Synthetic aperture particle image velocimetry (SAPIV) provides a non-invasive means of revealing the physics of complex flows using a compact camera array to resolve the 3D flow field with high temporal and spatial resolution. Intensity-threshold-based methods of reconstructing the flow field are unsatisfactory in nonuniform illuminated fluid flows. This article investigates the characteristics of the focused particles in re-projected image stacks, and presents a convolutional neural network (CNN)-based particle field reconstruction method. The CNN architecture determines the likelihood of each area containing focused particles in the re-projected 3D image stacks. The structural similarity between the images projected by the reconstructed particle field and the images captured from the cameras is then computed, allowing in-focus particles to be extracted. The feasibility of our method is investigated through synthetic simulations and experiments. The results show that the proposed technique achieves remarkable performance, paving the way for non-uniformly illuminated particle field applications in 3D velocity measurements.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
Recent advances in sensor techniques and improved 3D reconstruction algorithms have led to renewed interest in 3D flow field measurements [1–6,2,3], Nonintrusive 3D particle imaging velocimetry (PIV) techniques have been widely applied to the flow fields in both water and air, whereby the fluid domain is captured from different perspectives and then reconstructed algorithmically [7,8,9]. Stereo PIV visualizes the 3D fluid in a 2D velocity field [10,11], and scanning PIV extends such single-slice two-dimensional three-component (2D3C) measurements to multiple planes by combining a scanning laser sheet with a high-speed camera. However, the measurable velocity is below 1 m/s because of the laser repetition rate and the speed of the scanning mirror [9,12]. Defocusing digital PIV recovers depth information from defocused images, which are normally produced by a three-aperture mask, but has a low measurable particle density . Holographic PIV was developed as a truly volumetric velocity measurement to characterize the flow field , but is limited by its cumbersome experimental setup and small measurement volume while holograms are recorded by CCD/CMOS sensors [14,15]. Elsinga et al. used four cameras to capture particle images from different viewing angles and reconstructed the 3D particle field via the multiplicative algebraic reconstruction technique (MART); this approach is called Tomographic PIV (Tomo-PIV) . This provided a significant step in the development of 3D velocity measurement techniques. Typically, the seeding density of Tomo-PIV ranges from 0.02 to 0.08 ppp [17,18]. The spatial resolution is high, and Tomo-PIV can be used with a relatively large measurement volume (measurable range along the optical axis is smaller than that in the lateral direction) [19,20]. Tomo-PIV is widely used in experimental fluid mechanics , but the cameras must be arranged for a total aperture of 40°–80° to maximize the accuracy of the particle field reconstruction . When the optical measurement window size is finite, such as diagnosis of the internal flow field structure inside an engine, it becomes difficult to implement Tomo-PIV. MART-based reconstruction methods [21,22,23,24,25] are also time-consuming and storage-intensive. Light field PIV (LFPIV), or plenoptic PIV, is a single-camera 3D velocity measurement technique that employs a microlens array and a plenoptic camera [26,27,28,29]. Its resolution is limited by the number of microlenses in the array, and cannot match the seeding density or spatial resolution of Tomo-PIV because of its smaller angular resolution [30,31,32,33,34,35]. Research on LFPIV is an active field of study. In 2013, Schanz et. al. developed a Shake-The-Box (STB) method to track the particles in 3D [36,37]. It features a considerable improvement compared to previous Tomo PIV methods in both accuracy and the applicable particle image densities. STB uses the available velocity information to create an estimated particle distribution. The errors introduced by the estimation are small enough to be easily corrected using image matching, and no further partner search is required for known particles. This method has been successfully applied in a transitional jet in a water tank. Lawson et. al. proposed a hybrid method combining scanning PIV with tomographic reconstruction to make spatially and temporally resolved measurements of the fine-scale motions in turbulent flows . Three high-speed cameras were used to record particle images as a laser sheet was rapidly traversed across a measurement volume. However, the accuracy is strongly affected by the scanning speed. In 2017, Xiong et. al. developed a Rainbow PIV method [39,40]. A liner color filter was used to generate a set of wavelength varied light planes. Further, a diffractive optical element was custom designed to be attached in front of the camera lens, and it provided wavelength-selective focus to ensure all the light planes were in focus simultaneously. Therefore, a single camera is capable to retrieve 3D particle locations by combining 2D spatial location and 1D color information. The results in a real stirred flow showed the method can reconstruct a significant part of the flow structures robustly and accurately. In 2018, Schneiders et. al. developed a coaxial volumetric velocimeter (CVV) for wind tunnel measurements . In CVV, laser illumination was provided along the same direction as that of the camera views, reducing the optical access requirements to a single viewing direction. Helium-filled soap bubbles were used to increase the intensity of scattered light. The CVV system was employed for two wind tunnel experiments. The results showed the capability of the CVV to access the near-wall velocity over a complex three-dimensional topology. Paciaroni et. al. proposed a backscatter particle image velocimetry via optical time-of-flight sectioning in 2018 . The flow field was illuminated with an ultrashort laser, and the time-gated backscattered light was collected. Only a single-optical-access port was used in this method. Experimental results indicate this technique has a potential to quantitative PIV measurements in challenging single-optical-access geometries. Very recently, a single-camera 3D velocimetry based on endoscopic tomography was proposed . Nine projections of the target 3D luminous field at consecutive time instants were registered from different orientations with one camera and customized fiber bundles. ART algorithm and optical flow method were used to reconstruct the velocity field. A proof-of-concept experiment was conducted to validate its performance. The results showed the feasibility of the proposed method.
Synthetic aperture PIV (SAPIV) was developed as a truly 3D, three component (3C) technique for measuring 3D flow fields . In SAPIV, images captured from a camera array are remapped onto virtual planes. Combining images from each camera using SA imaging gives a sharply focused image at an arbitrary depth by exploiting the parallax between the cameras. Unlike other 3D3C PIV methods, SA imaging can visualize self-occluding scenes, meaning that higher particle densities (up to 0.123 ppp) can be tolerated. SAPIV employs a low-cost camera array, with all camera centers of projection lying in the same plane, and omits the complicated Scheimpflug adapters. This compact structure is useful when the views of the scene are confined to small windows, such as portholes for nuclear reactors or engine internal flow field structure diagnosis . SAPIV has been successfully used in fish kinetics [46,47,48], free-flying painted lady butterfly kinetics , and vocal fold vibrations , and has also achieved promising performance with bubble flows .
The SAPIV imaging system contains a camera array (usually 8–16 cameras) [19,44], with the principal axes of all cameras intersecting at one point in space. SA imaging simulates images captured via a lens of an arbitrary sized aperture. The experimental setup of SAPIV is shown in Fig. 1.
In principle, as the aperture size of a lens increases, its depth of field decreases . SA imaging simulates a lens with an arbitrarily large aperture using an array of cameras looking at the same scene from slightly different viewpoints. The SAPIV algorithm consists of remapping and 3D particle field reconstruction [19,44,53].
In the remapping step, images acquired from different cameras are re-projected onto different virtual planes, which are called refocused planes. The re-projected images are then accumulated and averaged with the map-shift-average algorithm . This is expressed as19].
There are several methods for extracting 3D focused particles from the 3D remapped volume. Belden et al.  assumed a Gaussian intensity distribution of each refocused image, and used a threshold of three standard deviations above the mean image intensity on each refocused plane as the minimum brightness of a valid particle. The stack of all thresholded focal planes yielded the final 3D particle volume for PIV processing. Lynch et al. modified this method to dilate the thresholded images by one pixel in each direction, ensuring that each particle image contains neighboring low-intensity information . However, this process is relatively inefficient, because the threshold for every refocused image must be determined separately. In low-density fields, even if there are no particles in the remapped plane, “focused particles” can still be reconstructed. Bajpayee et al. developed an algorithm that searches for local intensity peaks in the remapped volume and converges to the local maximum of each peak . A binarization technique with an intensity-weighted average scheme determines the location of different particles, although the binarization threshold must be set manually. If the threshold value is too low, many discrete blurs will not be removed. Conversely, if the threshold is too high, many actual particles will be filtered out. Thus, the particle positions will not be accurate, and the reconstruction accuracy will decrease. Qu et al. developed an adaptive threshold method to reconstruct the 3D particle field. The threshold is calculated by the correlation between the re-projected images of the reconstructed particle field and the images captured by different cameras .
These threshold-based methods assume that the intensities of different tracer particles in the volume are similar. Actually, in experimental implementations, when the point light source expands to 3D volume illumination, the intensity distribution will be nonuniform (see Fig. 2) [9,56–59]. Further, in water flow, the laser intensity decays rapidly along the incident direction. Hence, the particle intensities vary significantly. Particles with low intensities may be filtered out, meaning the field cannot be reconstructed with high accuracy. The focused particles cannot be extracted using the difference in intensity values. More general and extensive characteristics of focused particles must be investigated to enable highly accurate reconstructions of the particle field.
In this article, the spatial characteristics near the particle regions in the remapped image stacks are investigated using multiple perspective imaging. We propose a convolutional neural network (CNN)-based method to reconstruct the 3D tracer particle field in SAPIV. The structural similarity (SSIM) is applied to filter out areas with a low probability of focusing particles and extract in-focus particles. The principle of CNN-based particle field reconstruction is detailed in section 2. A CNN is trained with synthesized data that include variations in the calibration parameters of the camera array and particle sizes. In Section 3, a nonuniform illuminated vortex ring field is simulated to evaluate the performance of our method. In section 4, the trained CNN is applied to a cylinder wake flow.
2. CNN-based SAPIV particle field reconstruction
Re-projection from a 16-camera array during the SAPIV remapping process is illustrated in Fig. 3. A1, A2, A3 …, A16 are the 16 images of particle A in the array. Particles that lie in remapped plane Zi are aligned and appear sharp, whereas particles off this surface are blurred due to the parallax between the cameras. The intensity distribution along the depth of the particle in the remapped planes has an hourglass shape, with the focused particle lying at the waist of the hourglass. This hourglass is the inherent spatial characteristic of the refocused particle image. The right part of Fig. 3 shows the refocused image of a particle in the experiment. The focused particle lies in plane Zi, which is at the waist of the hourglass. The focused particle has a regular 3D structure along the depth range in the refocused image stacks.
In the refocused particle image stacks, nine neighboring sub-images of size 7 × 7 are cropped and combined into one 21 × 21 image. Figure 4 shows this process in 3D.
In Fig. 4(a), Z represents the neighboring refocused planes along the depth direction, and W represents the window with 7 × 7 pixels, (i, j, k) represents the 3D coordinate value of the center pixel in the window Wi, j, k. Nine neighboring sub-images of size 7 × 7 are cropped from neighboring planes Zk-4, Zk-3, Zk-2, Zk-1, Zk, …, Zk + 4. The cropped sub-images are represented as Ii, j, k-4, Ii, j, k-3, Ii, j, k-2, Ii, j, k-1, Ii, j, k, …, Ii, j, k + 4 (see Fig. 4(b)). Each cropped image is 7 × 7 pixels. The coordinate of the center of Ii, j, k is (i, j, k). Then, these nine images are combined into one image, (see Fig. 4(c)), of size 21 × 21 pixels. The center of the combined image is (11,11).
If the combined image contains a focused particle, the intensity distribution of the 2D image is similar to that in Fig. 4(c), and the focused particle image center in is also (11,11). The corresponding 3D location of the focused particle in the refocused image stacks can also be obtained, and the focused particle can eventually be extracted. By doing so, the 3D spatial characteristics of the focused particles in refocused image stacks can be transformed into 2D images. This approach has been successfully used in the multi-focus microscopy techniques. Abrahamsson et. al. developed a phase pattern to project images from diﬀerent focal planes onto diﬀerent regions of a camera chip. And then a 2D Gaussian fit on the z-projection of each focal plane was performed to locate the beads in the different sub-images . Herein, A CNN architecture is used to extract the focused particles from the combined 2D images. We believe this CNN method also has a potential to locate the 3D positions of single molecules for nano-scale molecular tracking .
The proposed method is illustrated in Fig. 5. The refocused images are first cropped in the depth, and then combined into numerous 21 × 21-pixel images. A trained CNN is then used to recognize images containing focused particles. A 3D gray centroid method is applied to relocate the particle center, and an adaptable threshold method is used to refine the particles. Ultimately, the 3D particle field is reconstructed. Our CNN-based method performs excellently in flow fields under nonuniform illumination. The CNN training process and specific operations are detailed below.
2.1 CNN construction and training process
This section describes the internal architecture of our CNN (see Fig. 6). The CNN architecture was determined by the combination of the existed CNN architectures [62,63] and a series of experiments by trial and error.
The multilayer CNN consists of one input layer, six convolutional layers, one pooling layer, one dropout layer, one fully connected layer, and one output layer. The input image is 21 × 21. In each convolution layer, zero padding, batch normalization (BN), and a rectified linear unit (ReLU) activation function are used. BN allows the system to achieve much higher learning rates, reduces the sensitivity to the initialization conditions, and reduces the internal covariate shift . In each convolution layer, the stride is set to be 1. The pooling layer (max pooling) uses a regularizing dropout unit to avoid overfitting, and the probability of neurons being retained is 0.7. In the fully connected layer, nodes are fully connected with those in the dropout layer. Finally, the possibility of focused particles is output by three nodes through a softmax classifier.
To train the CNN, we simulated six particle fields with different calibration parameters and different particle sizes. The six simulated fields with different calibration parameters and different particle size are illustrated in Fig. 7, and the calibration parameters and particle sizes are listed in Table 1. In each simulated particle field, the particle size u ranges from 100 to 400 μm. d represents the distance between the measurement volume and camera projection center, f represents the focal length of the simulated lens, and D represents the interval between each camera. In each simulation, the intensity of each particle is in Gaussian shape. Six refocused image stacks can then be calculated, from which the training data sets are generated.
In the simulation, a particle is divided into 10000 small 3D points, and the positions of these 10000 points are in Gaussian distribution. From the center to the edge, the density of these points varies from high to low. Each small point represents 1 intensity unit. Each small point is projected into the simulated camera sensor. The number of the projected points that fall on each pixel is counted and is calculated as the weighting of each pixel intensity. The synthetic particle images captured by each simulated camera can be eventually obtained (shown in Fig. 8). In Fig. 8, d represents the distance between the measurement volume and camera projection center, f represents the focal length of the simulated lens.
The selection of the standard deviation relates to the particle size. In this paper, if the particle diameter is between 100 μm and 200 μm, the standard deviation is randomly selected from 15 to 30. Further, if the particle diameter is between 201 μm and 300 μm, the standard deviation is randomly selected from 31 to 45. If the particle diameter is between 301 μm and 400 μm, the standard deviation is randomly selected from 46 to 60. The intensity of simulated particle images ranges from 10 to 255. In this way, the simulated particle image is more similar to the particle images in the real experiment, and the trained CNN architecture has better generalization ability.
The combined 21 × 21 images were obtained using the method in Fig. 4. If the focused particle image is located at the cropped image center, i.e., the intensity center of the focused particle is (11,11), it is labeled (1,0,0). If the intensity center location has one pixel offset in the X or Y direction, it is labeled (0,1,0). If the intensity center location is offset by one pixel in both the X and Y directions, it is labeled (0,0,1). Otherwise, the image is labeled (0,0,0) (see Fig. 9). If one of the three output nodes is approximately equal to 1, the input image is assumed to contain a focused particle. The 3D location of the focused particle is given by the coordinates of the corresponding point in the center of the image in the 3D particle field. These images and labels are used to train the CNN.
The CNN architecture was implemented in Python using the TensorFlow/Keras framework. All weights were updated using the Adam optimizer  with a learning rate of 0.01 over 15 epochs; the Loss Cross Entropy function was used to measure the error. Each simulated particle field contained 1000 particles, and 14562 training samples and 5715 test samples were synthesized. Figure 10 shows the training images generated by the six synthetic fields with labels (1,0,0); the particle sizes in the figure are 200 μm. The ratio of training to cross-validation was set to 7:3. In the training process, the threshold of softmax was chosen to 0.7. The training and blind testing of the network were performed on a Intel E31225 CPU and took about 18 minutes to converge. The accuracy of the validation data sets was 0.8089 and the loss was 0.1747.
2.2 CNN-based particle field reconstruction
The softmax function outputs the probability of the combined image containing a focused particle. If one of the three softmax nodes has a value of 1, the image center is regarded as the focused particle center, and its 3D coordinates on the refocused image stacks can be obtained. The 3D location is considered as the focused particle center. A threshold T is used to approximate the softmax outputs. If the output is greater than T, the value is set to 1; otherwise, it is set to 0. This is similar to a one-hot code. Unlike intensity-based threshold methods, we filter out areas with low probabilities of containing focused particles and retain areas with high probabilities. The trained CNN determines whether a region contains particles. We just select an appropriate threshold value of T to extract the focused particles and filter out the unfocused speckles. An adaptable threshold method is used herein.
Firstly, an initial threshold value T0 was used to filter out parts with low softmax output probabilities. Secondly, the 3D gray centroid of the obtained particle center was calculated according to (2). The particle center Pc is in the xth refocused plane, and its image coordinates are (x, y) in the xth refocused plane. (x0, y0, z0) represents the 3D gray centroid of the particle center Pc, S represents the 7 × 7 × 7 pixels of point (x ± 4, y ± 4) in the z ± 4 refocused images. For (i, j, k)∈S, f(i, j, k) represents the intensity value of (i, j) in the kth refocused plane.
Each tracer particle was assumed to have the same mean diameter, and the focused particle image size on the refocused plane was calculated using the principle of pinhole imaging model. Next, the particle centers were expanded into 3D particles. Finally, the non-maximum suppression (NMS) method was employed to remove redundant particles in the reconstructed field. If neighboring particles overlap by more than 50%, the particle with the lower center intensity is removed. Figure 11 shows the process of expansion and NMS, where different colors represent the intensities of the relocated particle centers. The red dot has the highest intensity. After expansion and NMS, the refocused particle is reconstructed.
SSIM described in [66,67] was used to evaluate the reconstruction quality:54]. As the images are 8-bit grayscale images, L was chosen as 255. The reconstructed particle field is projected into different camera sensors and the projected image is denoted as U1. The image captured by the camera array is denoted as U2.
For a fixed threshold T, the corresponding SSIM can be calculated. To obtain the optimal threshold T, several thresholds are set incrementally. For each threshold Ti, the mean SSIM value between the projected images from the reconstructed particle field and the images captured by the cameras is calculated, and cubic curve fitting is applied to obtain the maxima of the curve. When the mean SSIM value is maximized, the particle field has been reconstructed with high accuracy. Ultimately, 3D cross-correlation is used to calculate the 3D velocity vectors in the flow field.
3. Performance in vortex ring field
A vortex ring field was synthetized to verify our proposed method. In the simulation experiment, the CNN is trained using the data sets mentioned above. The performance of our proposed method is detailed below.
3.1 Vortex ring field with 16 cameras
In this experiment, a 16-camera array was simulated. The centers of projection lay in the same plane, and the optical axes of all lenses intersected at the center of the measurement volume. The interval D between each camera was set to 50 mm, and the distance between the camera array and the center of the simulated particle volume d was 350 mm. The cameras were modeled with 12-mm focal length lenses. Camera imaging was simulated via the pinhole model, and neither image distortion nor optical media changes along the line of sight were considered. The extrinsic parameters of the camera array were set manually, and the mapping function between the camera array and the particle volume was easily obtained.
A volume measuring 30 mm × 30 mm × 30 mm was seeded with 6000 particles of size 100 μm. During the simulation, the particle intensity decreased along the laser illumination, and the particle intensities at the edges were higher than those in the center of the volume. Three views of the particle fields at adjacent times are shown in Fig. 12.
The fluid motion is prescribed by the vortex ring equation used in [16,19,44], where the velocity magnitude is given byFigure 13 shows the movement of the particles in the 3D particle field across two frames. Low-intensity noise was seeded randomly to the synthetic camera images obtained from the pinhole model.
In the remapping process, the refocused plane was set to be 1292 × 964 pixels, with the square pixel Δp set to 3.75 μm × 3.75 μm. The reconstruction area was set to 480 × 480 × 480 voxels, and the size of each voxel was 80 μm. The refocused image stacks were reshaped using various 21 × 21-pixel images. The method has been detailed in Section 2. Further, to generate the patches of 21 × 21-pixel images, the window W was offset by 1 pixel. Therefore, 6 × 7 × 9 pixels were overlapped between two adjacent 21 × 21-pixel images. The trained CNN was then used to obtain the focused particle center locations and reconstruct the 3D particle field.
The particle size was calculated to be 3 × 3 × 3 voxels . After applying the CNN, the initial threshold T0 was set to 0.80, and the increment ΔT was set to 0.02. Following the 3D gray centroid operation and NMS, the reconstructed particle fields were obtained.
The relationship between T and the mean SSIM is shown in Fig. 14. In Fig. 14(a), when T = 0.87, the mean SSIM reaches a maximum of 0.8235; in Fig. 14(b), when T = 0.88, the mean SSIM reaches a maximum of 0.8548. Hence, the threshold in the first frame is chosen as 0.87, and that in the second frame is chosen to be 0.88.
The threshold method (denoted as Method 1)  and the adaptive threshold method (denoted as Method 2)  were also used to reconstruct the simulated particle field for comparison. In Method 1, the threshold was chosen to be two standard deviations above the mean image intensity on each refocused plane by trial and error. In Method 2, the adaptive value was determined after image pre-processing and several iterations. The reconstructed particle fields are shown in Fig. 15.
From the adjacent frames reconstructed by Method 1, we can see that particles with low intensities are missing as they have been filtered out as unfocused speckles. In the case of Method 2, many of the low-intensity particles have again been filtered out. However, the reconstructed frames given by the proposed method still contain many low-intensity particles (right half of each field).
The reconstruction quality Q  was applied to quantify the particle field reconstruction:
The reconstruction quality Q of the three methods for two adjacent frames is presented in Table 2. If Q ≥ 0.75, the reconstruction is considered to be acceptable . In Table 2, Rec. No. denotes the number of reconstructed particles, Frame 1 and Frame 2 represent the particle fields in adjacent frames. We can see that Method 1 only reconstructs about 4700 particles, and the reconstruction quality Q is less than 0.75. In Method 2, about 5000 of the 6000 particles are reconstructed, however, the reconstruction quality is less than 0.75. With our proposed method, about 5700 particles are reconstructed in both fields, and the Q values are relatively high at 0.8471 (frame 1) and 0.8543 (frame 2), both are greater than 0.75. Thus, we believe that the particle fields have been accurately reconstructed.
The reconstructed particle fields were projected into each camera sensor in the camera array. Figure 16 shows partial screenshots of the projected images from camera 1 of the two frames. It is clear that many of the particles have not been reconstructed by Methods 1 and 2. Particles with relatively low intensities have been filtered out as unfocused speckles. With our method, most of the particles have been reconstructed, including those with low intensities. The projected images of the particle field reconstructed by our method are most similar to the synthetic images.
The SSIM values between each projected image and the synthetic image are shown in Fig. 17. Our method gives the highest SSIM index in both frames. The other two reconstruction methods, particularly the threshold method, score quite low. Therefore, we believe intensity threshold-based methods are not suitable for particle field reconstruction under nonuniform illumination, whereas our method can reconstruct the corresponding 3D particle field with high accuracy.
To fully compare the performance of the three methods, the displacement fields were obtained from the three reconstructed volumes using a 3D cross-correlation-based multipass algorithm . A final interrogation volume containing 64 × 64 × 32 voxels and 75% overlap generated 16,800 vectors (20 × 20 × 42 vectors). The obtained displacement fields are shown in Fig. 18, where the fields are plotted in voxel units. From Figs. 18(a) and 18(b), it is apparent that Method 1 does not correctly reconstruct the vortex ring vector field. Figures 18(c) and 18(d) indicate that Method 2 reconstructs a vortex ring, but the cross-sections show that the vectors contain significant errors. In contrast, our proposed method (shown in Figs. 18(e) and 18(f)) clearly generates a vortex that can also be seen in the cross-sections.
To give a quantitative evaluation of the vortex fields reconstructed by the different methods, the error in each reconstructed vector field, defined as the difference between the velocity fields of the synthesized and reconstructed particle fields at every vector location, was calculated. The results are shown in Fig. 19. In the case of Method 1, the errors are randomly distributed, and most of the vector errors are greater than 1 voxel. With Method 2, about 60% of the vectors have an absolute error of less than 0.5 voxels. Using our proposed method, more than 95% of the vectors have an absolute error of less than 0.5 voxels, and ~80% of the vectors have an absolute error of less than 0.3 voxels. From this synthetic 3D experiment, it can be concluded that our SAPIV particle field reconstruction method is capable of instantaneously measuring flow structures with high resolution.
4. Experiment in a cylinder wake flow
4.1 Experimental setup and camera calibration
An array of 16 charge-coupled device (CCD) cameras was developed based on two customized industrial computers with 16 interfaces. The system consists of lenses (Computer M1214) with a 12-mm focal length and monochrome CCD cameras (Allied Vision Technologies, Guppy Pro F-125B: 1292 × 964, pixel size 3.75 μm). The cameras record at ten frames per second (fps). The cameras were adjusted such that all centers of projection were approximately in the same plane and the optical axes of all lenses approximately intersected at one point in the measurement volume. The spacing between each camera was set to 35 mm. A 10 W laser emitting at 808 nm (LASEVER) was used to illuminate the particle volume. A 30° Powell prism and a cylindrical lens with 50-mm focal length created a laser illumination cross-sectional area of 6 mm × 25 mm. A surface mirror was placed at the end of the tank to reflect the beam back into the experimental volume. A 350 mm × 350 mm × 200 mm glass tank filled with water was placed 300 mm from the camera array plane. The experimental setup is shown in Fig. 20.
The cameras were calibrated by imaging a precision-machined calibration plate passing through the target volume, as shown in Fig. 21. A square calibration plate (7 × 7 points with 10-mm point spacing) was precisely maneuvered through 11 planes in the volume in 2-mm increments covering the entire measurement domain to generate the 3D calibration points. The accuracy of traversal was of order ± 0.0025 mm. The image of the calibration plate center point at the middle calibration plane lay at the center of the image captured by each camera. The world coordinate system was attached to the calibration point located at the left central calibration plane.
Third-order polynomial fits were used to map the image coordinates to the world coordinates in each Z calibration plane, and linear interpolation was used to find the polynomial fits in the Z planes between each calibration plane. This approach follows that of prior Tomo-PIV studies, in which polynomial fits were used to overcome distortion introduced when imaging through an air–glass–water interface [19,69].
Calibration errors were obtained by calculating the distance between calibration points in the captured images and the points re-projected by the calibration functions. The maximum calibration error (Errmax) and the root mean square error (ErrRMS) of each camera are listed in Table 3. Belden et al.  noted that all reconstruction errors are less than 0.45 pixels in SAPIV, and we believe our calibration results are adequate for remapping and particle field reconstruction.
4.2 Cylinder wake flow
A smooth cylinder glass rod was placed in the water tank, and a flow was generated by a tube with an outlet orifice diameter dt = 20 mm. The diameter of the cylinder glass rod dr = 8 mm, and the flow velocity at the nozzle was approximately 0.5 m/s. The distance between the cylinder glass rod and the optical observation area l was ~60 mm. The cylinder glass rod lay on top of the laser volume, and the angle between the direction of water flow and the direction of laser transmission was 150° (see Fig. 22). The water temperature was 25°C. Water in the tank was seeded with 100-μm neutrally buoyant polyamide seeding particles. The frame rate was determined to be 10 fps at full resolution through the dropped frame test. Two adjacent time events of the experimental flow field were reconstructed.
The discrete grid on the remapped plane was chosen to have a cell size of 100 μm. During the remapping procedure, the spacing of neighboring remapped planes was set to be 100 μm. The reconstructed particle size was calculated to be 3 × 3 × 3 voxels . The image captured by the camera array is shown in Fig. 23. The particles are brighter in the edge region of the laser volume, whereas many particles in the center exhibit relatively lower intensities.
Image pre-processing was performed to remove background noise . The reconstruction area was 65 mm × 24 mm × 10 mm. The volume was discretized with 650 × 240 × 100 voxels, resulting in a resolution of 10 voxels/mm (voxel pitch of 100 μm). The remapping was accomplished using third-order polynomial calibration. The generation of the patches of 21 × 21-pixel images was the same as in the simulated experiment. The particle field was reconstructed by using the CNN architecture trained by the simulated particles. Before CNN process, the normalization was used in the generated 21 × 21-pixel image.
During the reconstruction process, the trained CNN and adaptable threshold operation were used to extract focused particles. In each field, T0 = 0.65 and ΔT = 0.02. The relationship between SSIM and the threshold in the two adjacent frames is shown in Fig. 24. After cubic curve fitting, the optimal threshold T was set to 0.75 in the first frame, as the mean SSIM attains a maximum at this point. The optimal threshold T was also chosen to be 0.75 in the latter frame.
The proposed method, the threshold method described in  and  (Method 1), and the adaptive threshold proposed in  (Method 2) were used to reconstruct the particle fields of two adjacent time events. In the threshold method, the threshold used to filter out unfocused speckles was chosen to be the standard deviation above the mean image intensity on each focal plane. The reconstructed particle fields are shown in Fig. 25. With Methods 1 and 2, many particles in the edge region of the laser volume are extracted, but many of the lower-intensity particles are filtered out. With the proposed method, many low-intensity particles in the central region are also extracted.
The projected images of the reconstructed particle fields were also calculated (see Fig. 26). The threshold method only reconstructs particles with high intensity. Many of the relatively darker particles are filtered as unfocused speckles. Approximately 1050 particles are reconstructed in frame 1 and around 850 particles are constructed in frame 2 when using the Method 1. In Method 2, more low-intensity particles are reconstructed, with approximately 2800 particles reconstructed in frame 1 and 3000 particles reconstructed in frame 2. Using our proposed method, many of the focused particles with very low intensities are reconstructed, with around 6000 particles reconstructed in frame 1 and 6200 particles reconstructed in frame 2.
To evaluate the three reconstruction methods, the SSIM between each image captured from the camera array and the images projected by the reconstructed field were calculated. The results are shown in Fig. 27. In frame 1, the mean SSIM value of Method 1 is 0.5392, that of Method 2 is 0.6081, and that of our method is 0.7112. In frame 2, the mean SSIM values of Method 1, Method 2, and the proposed method are 0.5254, 0.6051, and 0.7206, respectively. The SSIM acquired from our method is much higher than that for the other methods. Hence, our method provides the best performance in this experiment.
The 3D velocity field was calculated using cross-correlation-based 3D PIV analysis. A multi-pass algorithm with one pass at an initial interrogation volume size of 64 × 32 × 64 voxels and two passes at a final interrogation volume size of 32 × 32 × 32 voxels and 75% overlap generates a 27 × 4 × 78 vector field. The resultant vector field resolution is 2.1 mm in X and Y and 3.2 mm in Z. Linear interpolation and smoothing with a 3 × 3 × 3 Gaussian filter were also applied based on the codes in . The instantaneous velocity distribution is shown in Fig. 28.
In Fig. 28, the Z-axis corresponds to the laser illumination direction, and the wake field has just moved into the field of view. Different colors represent the velocity magnitudes. In the left of the field, the velocity is nearly 0, and the turbulence increases to the right. As the direction of the flow field is not parallel to the laser illumination, the flow field flows into view from the upper-right. A boundary can be clearly seen between the flow field and the water in the tank.
In the practical experiment, our cameras record at 10 frames per second at full resolution through the dropped frame test. Losing frames is strictly prohibited in SAPIV, so all the cameras in the array need to record the measurement volume synchronously. While in the wake flow, the flow velocity is about 0.5m/s. During 0.1s, the particles in the flow have displaced a considerable distance. According, the reconstructed particle fields at two adjacent frames will have low correlation, and the results based on 3D cross-correlations will contain many outliners. Also, the 10 frames will lead to undersampling of the vortex field. As a result, the details of the velocity field cannot be clearly reconstructed. Thus, the precision of the velocity field is low, and the resulting velocity is unreliable. Consequently, we reconstructed the flow when the wake flow just enters the observation area in the practical experiment. Nevertheless, this comparison of three methods has proved the non-substitutable advantage of our CNN-based reconstruction method in SAPIV.
This paper described a CNN-based particle field reconstruction method for SAPIV to resolve nonuniform illuminated field flows. The 3D spatial characteristics of the focused particles in refocused image stacks were investigated and transformed into 2D images, and a CNN architecture was used to reconstruct the particle field. The CNN architecture was trained using synthetic images that incorporate the various postures of the camera array and particle sizes. An SSIM-based threshold method was used to filter out unfocused speckles with small probability values and retain focused particles. The mean SSIM values of the images projected from the reconstructed particle field and the images captured by the camera array were calculated to measure the particle field reconstruction quality. A simulated 3D vortex ring field under nonuniform illumination was exploited to test the performance of the proposed method, which exhibited clear advantages over previous intensity threshold based methods. The method was also applied to a cylinder wake flow in water.
The proposed CNN-based approach can identify focused particles and filter out unfocused particles with high accuracy. The reconstruction quality Q of our method was shown to be much greater than 0.75, which is considered to be adequate. The accuracy of the obtained velocity was computed, and the results show that the flow field reconstructed by our method produces a high-precision velocity field. Further, the SSIM values between the images projected from the reconstructed particle field and the images captured by the cameras were calculated for both the synthetic simulation and the practical experiment. The results show that our method is capable of reconstructing the particle field with high accuracy in SAPIV. In the future work, GPU-accelerated implementation with NVIDIA graphics card needs to be investigated. Considering the capabilities of the proposed technique for SAPIV, we believe it can be widely applied in 3D flow field measurements.
National Natural Science Foundation of China (NFSC) (61701239); Natural Science Foundation of Jiangsu Province (BK20170852).
1. T. Yu, H. Liu, and W. Cai, “On the quantification of spatial resolution for three-dimensional computed tomography of chemiluminescence,” Opt. Express 25(20), 24093–24108 (2017). [CrossRef] [PubMed]
2. W. Cai and C. F. Kaminski, “Tomographic absorption spectroscopy for the study of gas dynamics and reactive flows,” Pror. Energy Combust. Sci. 59, 1–31 (2017). [CrossRef]
3. A. Bauknecht, B. Ewers, C. Wolf, F. Leopold, J. Yin, and M. Raffel, “Three-dimensional reconstruction of helicopter blade–tip vortices using a multi-camera BOS system,” Exp. Fluids 56(1), 1866 (2015). [CrossRef]
4. K. Mohri, S. Görs, J. Schöler, A. Rittler, T. Dreier, C. Schulz, and A. Kempf, “Instantaneous 3D imaging of highly turbulent flames using computed tomography of chemiluminescence,” Appl. Opt. 56(26), 7385–7395 (2017). [CrossRef] [PubMed]
5. Y. Song, J. Wang, Y. Jin, Z. Guo, Y. Ji, A. He, and Z. Li, “Implementation of multidirectional moiré computerized tomography: multidirectional affine calibration,” J. Opt. Soc. Am. A 33(12), 2385–2395 (2016). [CrossRef] [PubMed]
6. Y. Jin, Y. Song, X. Qu, Z. Li, Y. Ji, and A. He, “Three-dimensional dynamic measurements of CH* and C2* concentrations in flame using simultaneous chemiluminescence tomography,” Opt. Express 25(5), 4640–4654 (2017). [CrossRef] [PubMed]
7. K. Lasinger, C. Vogel, T. Pock, and K. Schindler, “Variational 3D-PIV with sparse descriptors,” Meas. Sci. Technol. 29(6), 064010 (2018). [CrossRef]
8. S. Discetti and F. Coletti, “Volumetric velocimetry for fluid flows,” Meas. Sci. Technol. 29(4), 042001 (2018). [CrossRef]
9. F. Scarano, “Tomographic PIV: principles and practice,” Meas. Sci. Technol. 24(1), 012001 (2013). [CrossRef]
10. A. K. Prasad, “Stereoscopic particle image velocimetry,” Exp. Fluids 29(2), 103–116 (2000). [CrossRef]
11. M. Arroyo and C. Greated, “Stereoscopic particle image velocimetry,” Meas. Sci. Technol. 2(12), 1181–1186 (1991). [CrossRef]
12. T. Hori and J. Sakakibara, “High-speed scanning stereoscopic PIV for 3D vorticity measurement in liquids,” Meas. Sci. Technol. 15(6), 1067–1078 (2004). [CrossRef]
13. F. Pereira, M. Gharib, D. Dabiri, and D. Modarress, “Defocusing digital particle image velocimetry: a 3-component 3-dimensional DPIV measurement technique. Application to bubbly flows,” Exp. Fluids 29(7), S078–S084 (2000). [CrossRef]
14. K. D. Hinsch, “Holographic particle image velocimetry,” Meas. Sci. Technol. 13(7), 201 (2002). [CrossRef]
15. J. Katz and J. Sheng, “Applications of holography in fluid mechanics and particle dynamics,” Annu. Rev. Fluid Mech. 42(1), 531–555 (2010). [CrossRef]
16. G. E. Elsinga, F. Scarano, B. Wieneke, and B. W. van Oudheusden, “Tomographic particle image velocimetry,” Exp. Fluids 41(6), 933–947 (2006). [CrossRef]
17. G. E. Elsinga, B. Wieneke, F. Scarano, and A. Schröder, “Tomographic 3D-PIV and applications,” in Particle image velocimetry (Springer, 2007), pp. 103–125.
18. F. Scarano and C. Poelma, “Three-dimensional vorticity patterns of cylinder wakes,” Exp. Fluids 47(1), 69 (2009). [CrossRef]
19. X. Qu, Y. Song, Y. Jin, Z. Li, X. Wang, Z. Guo, Y. Ji, and A. He, “3D SAPIV particle field reconstruction method based on adaptive threshold,” Appl. Opt. 57(7), 1622–1633 (2018). [CrossRef] [PubMed]
21. C. Atkinson and J. Soria, “An efficient simultaneous reconstruction technique for tomographic particle image velocimetry,” Exp. Fluids 47(4-5), 553–568 (2009). [CrossRef]
22. S. Discetti and T. Astarita, “A fast multi-resolution approach to tomographic PIV,” Exp. Fluids 52(3), 765–777 (2012). [CrossRef]
23. M. Novara, K. J. Batenburg, and F. Scarano, “Motion tracking-enhanced MART for tomographic PIV,” Meas. Sci. Technol. 21(3), 035401 (2010). [CrossRef]
24. K. Lynch and F. Scarano, “An efficient and accurate approach to MTE-MART for time-resolved tomographic PIV,” Exp. Fluids 56(3), 66 (2015). [CrossRef]
25. F. J. Martins, J.-M. Foucaut, L. Thomas, L. F. Azevedo, and M. Stanislas, “Volume reconstruction optimization for tomo-PIV algorithms applied to experimental data,” Meas. Sci. Technol. 26(8), 085202 (2015). [CrossRef]
26. T. W. Fahringer, K. P. Lynch, and B. S. Thurow, “Volumetric particle image velocimetry with a single plenoptic camera,” Meas. Sci. Technol. 26(11), 115201 (2015). [CrossRef]
27. T. Nonn, J. Kitzhofer, D. Hess, and C. Brücker, “Measurements in an IC-engine flow using light-field volumetric velocimetry,” in 16th Int. Symp. Appl. Laser Tech. Fluid Mech., Lisbon, Portugal (July 9–12, 2012).
28. E. A. Deem, Y. Zhang, L. N. Cattafesta, T. W. Fahringer, and B. S. Thurow, “On the resolution of plenoptic PIV,” Meas. Sci. Technol. 27(8), 084003 (2016). [CrossRef]
29. T. W. Fahringer and B. S. Thurow, “Filtered refocusing: a volumetric reconstruction algorithm for plenoptic-PIV,” Meas. Sci. Technol. 27(9), 094005 (2016). [CrossRef]
30. E. M. Hall, B. S. Thurow, and D. R. Guildenbecher, “Comparison of three-dimensional particle tracking and sizing using plenoptic imaging and digital in-line holography,” Appl. Opt. 55(23), 6410–6420 (2016). [CrossRef] [PubMed]
31. S. Shi, J. Wang, J. Ding, Z. Zhao, and T. New, “Parametric study on light field volumetric particle image velocimetry,” Flow Meas. Instrum. 49, 70–88 (2016). [CrossRef]
32. B. S. Thurow and T. Fahringer, “Recent development of volumetric PIV with a plenoptic camera,” in 10th Int. Symp. PIV (PIV’13), Delft, The Netherlands (July 1–3, 2013).
33. E. M. Hall, D. R. Guildenbecher, and B. S. Thurow, “Uncertainty characterization of particle location from refocused plenoptic images,” Opt. Express 25(18), 21801–21814 (2017). [CrossRef] [PubMed]
34. T. Fahringer and B. Thurow, “On the development of filtered refocusing: a volumetric reconstruction algorithm for plenoptic-PIV,” in 11th Int. Symp. PIV (PIV15), Santa Barbara, California (2015).
35. S. Shi, J. Ding, C. Atkinson, J. Soria, and T. New, “A detailed comparison of single-camera light-field PIV and tomographic PIV,” Exp. Fluids 59(3), 46 (2018). [CrossRef]
36. D. Schanz, S. Gesemann, and A. Schröder, “Shake-The-Box: Lagrangian particle tracking at high particle image densities,” Exp. Fluids 57(5), 70 (2016). [CrossRef]
37. D. Schanz, A. Schröder, S. Gesemann, D. Michaelis, and B. Wieneke, “‘Shake The Box’: A highly efficient and accurate Tomographic Particle Tracking Velocimetry (TOMO-PTV) method using prediction of particle positions,” in International Symposium on Particle Image Velocimetry (2013).
38. J. M. Lawson and J. R. Dawson, “A scanning PIV method for fine-scale turbulence measurements,” Exp. Fluids 55(12), 1857 (2014). [CrossRef]
39. J. Xiong, R. Idoughi, A. A. Aguirre-Pablo, A. B. Aljedaani, X. Dun, Q. Fu, S. T. Thoroddsen, and W. Heidrich, “Rainbow particle imaging velocimetry for dense 3D fluid velocity imaging,” ACM Trans. Graph. 36(4), 1 (2017).
40. J. Xiong, Q. Fu, R. Idoughi, and W. Heidrich, “Reconfigurable rainbow PIV for 3D flow measurement,” in 2018 IEEE International Conference on Computational Photography (ICCP) (IEEE, 2018), pp. 1–9. [CrossRef]
41. J. F. Schneiders, F. Scarano, C. Jux, and A. Sciacchitano, “Coaxial volumetric velocimetry,” Meas. Sci. Technol. 29(6), 065201 (2018). [CrossRef]
42. M. E. Paciaroni, Y. Chen, K. P. Lynch, and D. R. Guildenbecher, “Backscatter particle image velocimetry via optical time-of-flight sectioning,” Opt. Lett. 43(2), 312–315 (2018). [CrossRef] [PubMed]
44. J. Belden, T. T. Truscott, M. C. Axiak, and A. H. Techet, “Three dimensional synthetic aperture particle image velocimetry,” Meas. Sci. Technol. 21(12), 125403 (2010). [CrossRef]
45. D. M. Kubaczyk, A. H. Techet, and D. P. Hart, “Assessment of radial image distortion and spherical aberration on 3D synthetic aperture PIV measurements,” Meas. Sci. Technol. 24(10), 105402 (2013). [CrossRef]
46. L. Mendelson and A. H. Techet, “3D synthetic aperture PIV of a swimming fish,” in 10th Int. Symp. PIV (PIV’13), Delft, The Netherlands (July 1–3, 2013), pp. 1–6.
47. L. Mendelson and A. H. Techet, “Quantitative wake analysis of a freely swimming fish using 3D synthetic aperture PIV,” Exp. Fluids 56(7), 135 (2015). [CrossRef]
48. L. Mendelson and A. H. Techet, “Multi-camera volumetric PIV for the study of jumping fish,” Exp. Fluids 59(1), 10 (2018). [CrossRef]
49. K. R. Langley, E. Hardester, S. L. Thomson, and T. T. Truscott, “Three-dimensional flow measurements on flapping wings using synthetic aperture PIV,” Exp. Fluids 55(10), 1831 (2014). [CrossRef]
50. J. R. Nielson, D. J. Daily, T. T. Truscott, G. Luegmair, M. Döllinger, and S. L. Thomson, “Simultaneous tracking of vocal fold superior surface motion and glottal jet dynamics,” in ASME 2013 Int. Mech. Eng. Congr. Expo., San Diego, California (2013), pp. V03BT03A039. [CrossRef]
51. J. Belden, S. Ravela, T. T. Truscott, and A. H. Techet, “Three dimensional bubble field resolution using synthetic aperture imaging: application to a plunging jet,” Exp. Fluids 53(3), 839–861 (2012). [CrossRef]
52. A. Bajpayee and A. H. Techet, “Fast volume reconstruction for 3D PIV,” Exp. Fluids 58(8), 95 (2017). [CrossRef]
53. V. Vaish, B. Wilburn, N. Joshi, and M. Levoy, “Using plane + parallax for calibrating dense camera arrays,” in Proc. IEEE Conf. Comput. Vision Pattern Recog. (IEEE, 2004), pp. I-2–I-9. [CrossRef]
54. K. Lynch and B. Thurow, “Preliminary development of a 3-D, 3-C PIV technique using light field imaging,” in AIAA Fluid Dyn. Conf. Exhib., Hawaii (2013).
55. A. Bajpayee, “3D particle tracking velocimetry using synthetic aperture imaging,” M.S. thesis (Massachusetts Institute of Technology, 2014).
56. F. Zhang, Y. Song, X. Qu, Y. Ji, Z. Li, and A. He, “Three-dimensional illumination system for tomographic particle image velocimetry,” in Optical Design and Testing VII (International Society for Optics and Photonics, 2016), p. 100211D.
57. A. Schröder, R. Geisler, G. E. Elsinga, F. Scarano, and U. Dierksheide, “Investigation of a turbulent spot and a tripped turbulent boundary layer flow using time-resolved tomographic PIV,” Exp. Fluids 44(2), 305–316 (2008). [CrossRef]
58. F. Scarano and C. Poelma, “Three-dimensional vorticity patterns of cylinder wakes,” Exp. Fluids 47(1), 69 (2009). [CrossRef]
59. X. Qu, Y. Song, L. Xu, Y. Jin, Z. Li, and A. He, “Analysis of laser intensity attenuation and compensation and the influence on imaging through particle field,” in AOPC 2017: Optical Sensing and Imaging Technology and Applications (International Society for Optics and Photonics, 2017), p. 104622N.
60. S. Abrahamsson, J. Chen, B. Hajj, S. Stallinga, A. Y. Katsov, J. Wisniewski, G. Mizuguchi, P. Soule, F. Mueller, C. Dugast Darzacq, X. Darzacq, C. Wu, C. I. Bargmann, D. A. Agard, M. Dahan, and M. G. Gustafsson, “Fast multicolor 3D imaging using aberration-corrected multifocus microscopy,” Nat. Methods 10(1), 60–63 (2013). [CrossRef] [PubMed]
61. S. Ram, J. Chao, P. Prabhat, E. S. Ward, and R. J. Ober, “A novel approach to determining the three-dimensional location of microscopic objects with applications to 3D particle tracking,” in Three-Dimensional and Multidimensional Microscopy: Image Acquisition and Processing XIV (International Society for Optics and Photonics, 2007), p. 64430D.
62. Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” Proc. IEEE 86(11), 2278–2324 (1998). [CrossRef]
63. K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” arXiv preprint arXiv:1409.1556 (2014)
64. S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” arXiv preprint arXiv:1502.03167 (2015).
65. D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980 (2014).
66. Y. Wu, Y. Rivenson, Y. Zhang, Z. Wei, H. Günaydin, X. Lin, and A. Ozcan, “Extended depth-of-field in holographic imaging using deep-learning-based autofocusing and phase recovery,” Optica 5(6), 704–710 (2018). [CrossRef]
67. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process. 13(4), 600–612 (2004). [CrossRef] [PubMed]
68. J. Belden, T. T. Truscott, M. C. Axiak, and A. H. Techet, “Synthetic aperture imaging core code package,” http://saimaging.byu.edu/2012/11/27/core-code-package-release-3.
69. S. M. Soloff, R. J. Adrian, and Z. C. Liu, “Distortion compensation for generalized stereoscopic particle image velocimetry,” Meas. Sci. Technol. 8(12), 1441–1454 (1997). [CrossRef]