RainbowPIV with improved depth resolution—design and comparative study with TomoPIV

In order to fully reconstruct dense 3D-3C velocity fields, prior work has utilized a pipeline approach that first estimates the particle distribution fields at successive time steps, and then reconstructs the corresponding flow fields using the estimated particle distributions. Separating particle distribution field and velocity field reconstructions neglects temporal coherence as a strong physical cue. Specifically, particles present at one time step should also be present at the next, as well as the previous time steps (excluding a small number of particles entering or leaving the observation volume), and their location should be consistent with the estimated 3D-3C flow fields. RainbowPIV was the first approach to utilize a joint optimization framework for the reconstruction of particle distribution fields and fluid velocity vector fields on sequentially captured video frames. As illustrated in [26], the benefit of this prior information is the refinement of any ambiguous particles computed from a single frame, resulting in more precise flow fields. This joint optimization framework has also been utilized for multi-camera 3D Fluid Flow Estimation [9] and x-ray computed tomographic applications [29], and has demonstrated improved reconstruction quality in the respective works. Specifically, the integrated optimization problem for RainbowPIV can be expressed as:

$$\begin{align} ( \mathbf{p}^{*}, \mathbf{u}^{*} ) & = \mathop{\mathrm{argmin}}_{\mathbf{p},\mathbf{u}}~ \frac{1}{2} \left\| \mathbf{A} \Bigg{[} \mathbf{p}_{1}|...|\mathbf{p}_{T} \Bigg] - \Bigg[ \mathbf{i}_{1}|...|\mathbf{i}_{T}\Bigg{]} \right\|_{2}^{2} \nonumber\\ & + \kappa_{1} \left\| diag(\boldsymbol{w})\left( \mathbf{p}_{1};...;\mathbf{p}_{T} \right) \right\|_{1} + \Pi_{[0,1]} \left(\mathbf{p}_{1};...;\mathbf{p}_{T} \right) \nonumber\\ & + \kappa_{2} \sum_{t = 1}^{T} \int_{\Omega} \left( \mathbf{p}_t - \mathbf{p}_{t+1}( \mathbf{u}_t, -\Delta t) \right) ^{\circ 2} d \Omega \nonumber \\ & + \kappa_{3} \sum_{t = 1}^{T} \left\|\nabla \mathbf{u}_{t}\right\|^{2}_{2} + \sum_{t = 1}^{T} \Pi_{\mathcal{C}_{DIV}}(\mathbf{u}_{t}) \nonumber \\ & + \kappa_{4} \sum_{t = 1}^{T} \Big{(} \left\| \mathbf{u}_{t}- \Pi_{\mathcal{C}_{DIV}}(\mathbf{u}_{t-1}(\mathbf{u}_{t-1}, \Delta t)) \right\|_{2}^{2} \nonumber \\ & + \left\| \mathbf{u}_t- \Pi_{\mathcal{C}_{DIV}}(\mathbf{u}_{t+1}(\mathbf{u}_{t}, -\Delta t)) \right\|_{2}^{2} \Big),\end{align} \tag{ 1 }$$

where $\mathbf{p}^{*}$ and $\mathbf{u}^{*}$ denote the target particle distribution field and the fluid velocity vector field, respectively. The above energy function is minimized using a coordinate descent method, alternating on these two variables. Specifically, we alternate between keeping $\mathbf{u}$ fixed, and solving for $\mathbf{p}$, then fixing $\mathbf{p}$ and solving for $\mathbf{u}$. The following sections provide intuitive motivations and explanations for the above optimization problem. More detailed derivations and solutions for this issue can be found in [26].

From a starting point of extracting particle distributions from a single RGB image, straightforward approaches have been proposed in [11, 23], which relate hue values of the captured image to particle depths in the volume. This type of approach, as discussed before, cannot correctly extract the locations of overlapped particles due to color mixing. Instead, we express the imaging system as a forward model (illustrated in figure 3), written by $\mathbf{A} \mathbf{p}_t$, with $\mathbf{A}$ being the coefficient matrix, and $\mathbf{p}_t$ being the particle distribution field at time step t (1 ≤ t ≤ T). Denoting the corresponding captured image as $\mathbf{i}_t$, ideally, we have $\mathbf{A} \mathbf{p}_t = \mathbf{i}_t$. However, this linear system is under-determined (or ill-posed), since we have fewer constraints than unknown variables.

Zoom In Zoom Out Reset image size

Introducing regularization terms and constructing a minimization problem are common strategies to tackle this kind of problem. A sparse distribution of particles in the volume is expected, and the value of the particle distribution fields should be within the range of 0 and 1. These two constraints account for the second and third terms in line 2 of equation (1), where $\| \cdot \|_1$ is the L₁ norm, and κ₁ is the penalty parameter enforcing the sparsity. Here, $diag(\boldsymbol{w})$ denotes the weights compensating the sensor sensitivity to different wavelengths [25]; Π_[0,1] is the operator projecting the value to the convex set [0, 1]. The fourth summation term in equation (1), referred to as the particle occupancy consistency, incorporates temporal coherence into the minimization problem, where $\mathbf{p}_{t+1}( \mathbf{u}_t, -\Delta t)$ indicates the advection of the particle distribution fields at time step t + 1 under flow field $\mathbf{u}_t$ by −Δt units of time. Note that this term is also used for reconstructing velocity fields.

Having captured the particle fields in sequential time steps, we can then perform velocity estimation, in order to generate time-resolved velocity fields. A common method for such velocity estimation is Digital Volume Correlation (DVC), based on the cross-correlation of volume neighborhoods, or windows [14]. However, the selection of an appropriate window size can be difficult, since large windows result in overly smooth flow fields, while small windows can result in incorrect matches, particularly in the case of noisy, low-light images. Furthermore, such algorithms require a sufficiently high density of particles, which cannot always be satisfied for every flow scenario, especially in 3D measurements. Moreover, the cross-correlation based framework does not easily support adding extra physical constraints, such as divergence free flows in the case of incompressible fluids.

An alternative method is therefore adopted from the computer vision community, which solves the variational 'optical flow' problem (corresponding to 2D digital image correlation), by finding the correspondence between two successive particle distribution fields. The Horn-Schunck algorithm is the most widely used global variational optical flow method; this assumes a consistency of brightness between successive images, together with the smoothness of the flow vectors in the spatial domain. Here, instead of assuming the photo-consistency of 2D image intensities (the brightness consistency constraint does not hold, since the color of particles will change when they traverse different depth levels in axial direction), our 3D version of Horn-Schunck assumes a consistency of particle distribution fields $\mathbf{p}$, as described in section 3.2. In this way, the particle occupancy consistency term ties particle and velocity reconstruction together into a joint optimization problem.

The framework described above is also amenable to the easy incorporation of the physical parameters of a given fluid. A number of works have already taken the physical constraints governed by Navier–Stokes equations into account when computing flows [6, 7, 16, 22, 26, 28]. In this work, we consider the fluids to be nonviscous and incompressible, which indicates zero divergence in the flow vectors. Moreover, the incompressible Navier–Stokes equation describes the time evolution of fluid flows. These two physical parameters account for the last two terms of equation (1). Here, $\Pi_{\mathcal{C}_{DIV}}$ is the operation to project arbitrary flow vector fields onto a space having a divergence-free velocity field, where $\mathbf{u}_{t+1}(\mathbf{u}_{t}, -\Delta t)$ and $\mathbf{u}_{t-1}(\mathbf{u}_{t-1}, \Delta t)$ approximate the time evolution of the flow vector fields backward and forward, respectively. Note that we are not using the Navier–Stokes Equations to simulate flows; instead, we are incorporating them into the optimization problem as regularization terms.

A coarse-to-fine strategy is applied in our flow computation framework to solve large displacements. Roughly computed flow fields in the coarsest level are transmitted into finer levels as an initial guess, and are then iteratively updated, based on the fixed-point theorem (here, a pyramid of four levels and 10 inner iterations are applied). The algorithm stops when it reaches the finest level,where one vector per voxel is achieved.

The depth resolution of the RainbowPIV system is dependent on the number of discretized depth levels, which is basically the number of PSFs. In our implementation, PSFs are expressed by a batch of 3 × 5 × 5 matrices, representing three color channels and a 5 × 5 spatial window. Owing to the relatively low sensitivity of the camera sensor to subtle wavelength shifts, only a limited number of PSFs is calibrated in the original RainbowPIV system, resulting in a relatively low depth resolution.

Specifically, the axial reconstruction accuracy relies on two factors: the camera's wavelength sensitivity, and the quantization error. The sparser the sampled signal, the larger the resulting quantization error. In this work, we show that it is feasible to reduce the quantization error by increasing the number of PSFs, using a simple linear interpolation scheme to generate in-between PSFs. As indicated in figure 4, in the prior implementation of RainbowPIV, the sampling rate in the axial direction is far lower than that in lateral directions, giving rise to a relatively large quantization error. With the proposed interpolation scheme, the same sampling rate can be achieved in all directions. While this change in PSF with wavelength is not linear over a large spectral range, piecewise linear approximations such as ours provide a good model of the true intensity change in each RGB color channel (see figure 5) over small spectral bands. With this simple digital adaptation (no hardware adjustment is required), we are able to generate super-depth particle distribution fields and flow vectors with reduced axial quantization errors. A similar interpolation scheme is also applied to reconstruct subpixel flow vectors between two images [8].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now assess the effectiveness of this interpolation scheme on simulated data. Here, we simulate a volume with a fixed lateral dimension of 20 mm × 20 mm, and various axial dimensions (2 mm, 4 mm, 8 mm, 16 mm), yielding different depth-to-width ratios. We use the calibrated PSFs to model the camera's response to different depth levels. The lateral dimension is discretized by a size of 200 × 200, yielding a spacing of 0.1 mm. A random particle distribution of ppp = 0.05 is generated, and is advected following simulated flow velocities, as described in [26]. We compare the reconstruction errors of the flow vectors in the axial direction (referring to $\mathbf{u}_z$) for both standard RainbowPIV and the proposed depth super-resolved RainbowPIV. The results are shown in figure 6. When the depth-to-width ratio is 0.1, RainbowPIV and the new proposed approach are identical, as they have the same axial spacing. When the depth-to-width ratio increases, RainbowPIV samples the volume more sparsely in the axial direction, yielding increased axial reconstruction errors (almost linear to the depth-to-width ratio). Using the interpolated scheme, the axial reconstruction accuracy is significantly improved, based on the simulated results; errorsare reduced by roughly half when the depth-to-width ratio reaches 0.8. The accuracy of the reconstructed flow vectors using this approach is further experimentally validated below.

Zoom In Zoom Out Reset image size

Here, we validate the proposed strategy using experimental data, including information relating to ground truth motion, based on the approach in [26]. We immerse the particles into a tank containing high viscosity liquid to "freeze" them. The tank is then put on a multi-dimensional rotation/translation stage. The true movement of particles can thus be managed by manipulating the stage. We move the tank in x and z directions, respectively, so as to compare the accuracy of both lateral and longitudinal reconstruction . A translation of 0.2 mm in the x direction and 0.5 mm in the z direction is applied. As a result, the standard RainbowPIV implementation [26] recovers a motion vector of 0.18 mm/time units, the mean value of the norm of the reconstructed velocity in the lateral translation, with a standard deviation of 0.015 (mm/time unit), and the proposed method delivers a similar reconstruction accuracy, 0.18 mm/time unit with standard deviation of 0.014 (mm/time unit). With regard to the longitudinal translation, 0.35 mm/time unit, with a standard deviation of 0.079 (mm/time unit) is achieved by the original RainbowPIV, while our super-resolved approach achieves 0.42 mm/time unit with a standard deviation of 0.045 (mm/time unit). We therefore observe a significant improvement with respect to axial reconstruction accuracy in this simple scenario, while the lateral reconstruction accuracy is unaffected.

In both implementations, the recovered lateral flows exhibit greater accuracy than the axial flows. As explained in section 3.4, quantization errors and camera wavelength sensitivity both affect the axial resolution. In the case of the original RainbowPIV, the sampling rate along the axial direction is lower than the lateral sampling rate, resulting in larger quantization errors. The super-resolved approach, in contrast, achieves an identical sampling rate in lateral and axial directions; however, the camera is more sensitive to lateral particle motions (changes in pixel positions) than axial motions (changes in wavelength). This fact brings in a great degree of uncertainty in relation to the estimation of axial flows.

We further validate our method by directly comparing it with a four-camera Tomo-PIV system in relation to practical fluid flow phenomena. The simultaneous experimental setup for four-camera Tomo-PIV and single-camera RainbowPIV is shown in figure 7. The effective resolution of the captured RainbowPIV images is 1600 × 1120 (pixels). We downsample the images by a factor of 8 in both dimensions (no further image preprocessing required), resulting in 200 × 140 (pixels) images, and each particle occupies roughly 3 × 3 pixels. The depth dimension is discretized into 100 levels in order to obtain regular voxels, where the PSFs of 20 levels (1.25 mm per level) are obtained by calibration, and the intermediate PSFs are digitally generated by linear interpolation, based on the calibrated neighbor PSFs. Therefore, the size of the reconstruction grid is 200 × 140 × 100, with a voxel spacing of 0.25 mm along all three axes. One vector per voxel will be generated by the proposed algorithm; therefore, the reconstructed velocity field has the same dimensional size and spatial resolution as the grids (200 × 140 × 100). The particle image density is roughly 0.05 ppp (particles per pixel), and one example of the captured RainbowPIV image at this particle density is shown in figure 8. A qualitative comparison of the reconstructed velocity fields for the introduced vortex ring at one time step is presented in figure 9. Two sliced views (one perpendicular and the other parallel to the image plane) of the reconstructed flow vectors and the color-coded vorticity magnitude are visualized. The overall flow structure obtained by RainbowPIV agrees well with that computed by Tomo-PIV. Nevertheless, with its multiple perspectives, Tomo-PIV demonstrates greater depth resolving capability than the proposed single view approach. We can observe that the reconstructed flow vectors from RainbowPIV with a large z component are noisier than flow vectors in the x − y plane. As discussed in section 3.4, this is due to the relatively lower sensitivity of the camera to the wavelength change, so that the camera is less sensitive to particle motions in the axial direction than to in-plane motions. Quantitative comparisons between Tomo-PIV and RainbowPIV show an average difference of about 0.05 m s⁻¹ for flow vector components in the x − y plane, and 0.1 m s⁻¹ for vector components in the z direction, with a maximum flow magnitude of 0.53 m s⁻¹. In comparison, the original RainbowPIV implementation has the same in-plane average difference, whereas the out-of plane difference is 0.21 m s⁻¹. Although the uncertainty of the axial flow vectors is still larger than for the lateral flow vectors, these results affirm the improved axial resolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An additional comparison is conducted by visualizing the isosurface of the vorticity magnitude, as shown at the top of figure 10. The figures reveal the similarity of the core structures reconstructed by these two measurement technologies. We further verify the mass conservation properties of the reconstructed flow fields, which should comply with their physical properties. The divergence of the computed flows is shown in the bottom of figure 10. Zero divergence is expected everywhere for incompressible fluids. This divergence-free property is explicitly enforced by our reconstruction method, whereas Tomo-PIV fails to generate flow fields obeying this physical property, and in addition does not readily support the addition of this constraint during computation. This result demonstrates that, despite the overall high quality and detail of the Tomo-PIV solution, it can in fact not be considered as a ground-truth solution.

Zoom In Zoom Out Reset image size

Next, we evaluate the RainbowPIV system in low particle density situations, which are far below the desired particle density required by Tomo-PIV. The correlation-based algorithms applied in Tomo-PIV require sufficiently dense seeding particles to extract accurate flow fields, and usually perform poorly at low particle densities. The variant of the optical flow method presented here, however, ameliorates this issue by exploiting a local constraint (particle occupancy consistency), which ensures that in those regions where particles are present, the reconstructed local flow vectors match with the practical particle motion, and two global constraints (global smoothness constraint and temporal coherence), transmitting the accurate local flow vectors to those regions where particles are not present. Moreover, the physical properties are always satisfied, regardless of the particle density.

The tested particle image densities are roughly calculated as 0.005 and 0.015 ppp, and the captured images for these are shown in figure 12. The reconstructed flow vectors and vorticity magnitude are also visualized in figure 11, which demonstrates that RainbowPIV successfully captures the expected flow structures for vortex rings at various low density levels. Owing to its use of both local and global constraints, the proposed method delivers decent results at rather low particle densities. Under such conditions, particle tracking systems would be preferable to Tomo-PIV. However, this generates Lagrangian flow vectors with sparse descriptors, rather than the desired Eulerian vector fields. As a result, we believe that RainbowPIV can also be very competitive in experiments where uniform particle seeding is difficult or impossible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size