a. Main idea

The main procedures of the hybrid method are briefly described, as illustrated in Figure 1. From a pair of successive PIV images (image 1 and image 2), the cross-correlation method is used for initial estimation of a coarse-grained displacement field since it as an integral approach is less sensitive to the displacement magnitude, random noise and illumination change in PIV images. Then, a shifted image 1 is generated by using a shifting scheme based on the coarse-grained displacement field, where the spatial resolution of the original image 1 is restored by using interpolation. Further, from the shifted image 1 and the image 2, a residual displacement field is extracted by using the optical flow method since the residual displacement magnitude is very small. The corrected displacement field is reconstructed by superposing the coarse-grained displacement field and the residual displacement field. One iteration is required at least, and the procedures can be iterated for further improvement.

b. Cross-correlation method

The application of the correlation method is the first step of the hybrid method. The cross-correlation method is physically intuitive [1, 2]. An image domain is decomposed into many small subdomains called interrogation windows. A domain-averaged displacement vector in a window is determined by maximizing the cross-correlation function. The cross-correlation function can be directly computed. Since the cross-correlation function is essentially a convolution integral, it can be evaluated via the fast Fourier transform (FFT) algorithm. There are many research and commercial algorithms for PIV [3, 4, 5]. In this work, in order to make a stand-alone Matlab program of the hybrid method, we adopt the modified functions from the open-source Matlab software PIVlab [16]. The structure and performance of PIVlab were described by Thielicke and Stamhuis [17]. For operational simplicity, the main parameters retained in PIVlab are the window sizes in multiple paths and the window overlapping percentages.

c. Optical flow method

The optical flow equation in the image plane for PIV images is written as [6]

where I is the normalized image intensity, u = (u₁,u₂) is the optical flow, and ∇ = ∂/∂x_i (i = 1,2) is the gradient operator in the image plane (x₁, x₂). The boundary term $f=w_{\mathrm{ext}}{\langle C_{\mathrm{ext}}\rangle }_{a}n\cdot {\left(\mathrm{NU}\right)|}_{{\Gamma }_1}^{{\Gamma }_2}$M2 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ f = {w_{ext}}{\left\langle {{C_{ext}}} \right\rangle _a}{\boldsymbol n} \cdot \left. {\left({N{\boldsymbol U}} \right)} \right|_{{\Gamma _1}}^{{\Gamma _2}} \] \end{document} acts as a source/sink term representing the effect of particles accumulated within the illumination sheet confined by the two control planes Γ₁ and Γ₂, where U is the particle velocity in the 3D object space, N is the number of particles per unit total volume, n is the unit vector normal to Γ₁ and Γ₂, C_ext is the extinction cross section, and W_ext is the corresponding correlation coefficient. The optical flow in Eq. (1) is defined as u = (u₁,u₂) = γ ⟨U₁₂⟩_N, where ⟨U₁₂⟩_N is the light-path-averaged velocity of particles projected onto the image plane, and γ is a factor of projection. Essentially, Eq. (1) is a differential equation describing the projected motion of particles in the image plane. It is assumed that the net flux of particles through the control surfaces is zero, i.e., $n\cdot {\left(\mathrm{NU}\right)|}_{{\Gamma }_1}^{{\Gamma }_2}=0$M3 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\boldsymbol n} \cdot \left. {\left({N{\boldsymbol U}} \right)} \right|_{{\Gamma _1}}^{{\Gamma _2}} = 0 \] \end{document} .

To determine the optical flow as an inverse problem from Eq. (1), a variational formulation with a smoothness constraint is typically used [6, 18, 19]. Given I and f, we define a functional

where the first integral term is the equation functional, the second integral term is the first-order Tikhonov regularization functional, D is an image domain of interest, and α is the Lagrange multiplier. Minimization of the functional J(u) leads to the Euler-Lagrange equation

where ∇² = ∂²/∂x_i ∂x_i is the Laplace operator, and α is the Lagrange multiplier that controls the smoothness of the field. When a pair of sequential images is given and f is neglected in the first-order approximation, the standard finite difference method is used to solve this partial differential equation with the Neumann condition ∂u/∂n = 0 on the image domain boundary ∂D for the optical flow u. In theory, the optical flow method is able to extract a velocity field at the spatial resolution of one vector per pixel. However, the actual resolvable spatial scale is somewhat affected by the Lagrange multiplier that acts as an effective diffusion coefficient, i.e., the larger Lagrange multiplier tends to smooth out finer features. The optical flow algorithm used in this work is the same as the open source program “OpenOpticalFlow.v1” [14]. The program has the Horn-Schunck estimator for an initial solution [18] and the Liu-Shen estimator for a refined solution [6, 19]. Optical flow computation sensitivity to the Lagrange multiplier is not large. The hybrid method is described in [20, 21], which is recapitulated for the purpose of application.

d. Parameters related to error

A systematical error analysis for the optical flow method applied to PIV images is given by Liu et al. [13]. The accuracy of the optical flow method for PIV depends on the four parameters: particle displacement, particle velocity gradient, particle image density, and particle image diameter. The expression of error estimation is given below, which is applicable to not only the optical flow method but also the cross-correlation method. The total error of the optical flow method applied to PIV images is given by

where c₁, c₂, c₃ and c₄ are empirical coefficients. In Eq. (4), the main parameters related to the errors in PIV are the particle displacement magnitude ǁΔx_Pǁ, the particle image diameter d_p, the particle velocity gradient ǁ∇u_Pǁ, and the particle image density N_p.

The error is proportional to the particle displacement magnitude ǁΔx_Pǁ. The effect of the particle image diameter d_p in Eq. (4) is represented by the terms of $~d_p^2$M7 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\sim} {d_{p}^{2}} \] \end{document} and $~d_p^{-2}$M8 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\sim} {d_{p}^{-2}} \] \end{document} that have the opposite trends as d_p increases, which indicates that there could be the optimum.diameter. Generally, the error decrease with increasing the particle image density N_p. But when N_p becomes very large, the PIV images become more uniform, and thus extraction of the optical flow is not accurate due to the very small intensity gradient. Therefore, based on this observation, there would be an optimal value of N_p. For the convenience of application, the particle displacement ǁΔx_Pǁ is replaced by the maximum particle displacement max(ǀΔx_Pǀ), and the particle velocity gradient magnitude ǁ∇u_Pǁ is replaced by the averaged magnitude of the particle velocity gradient ǁ∇>u_Pǁ.

e. Shifting scheme

The shifting scheme in the image plane is a necessary element of the hybrid method, which is used to generate the shifted image 1 for optical flow computation. The shifting scheme has two procedures. The first procedure is pixel-to-pixel mapping of the intensity of a given image onto the shifted image denoted by I₀ based on a given velocity field in pixels per unit time.

To include the contribution of subpixel motion, the second procedure is described below. Since the coarse-grained velocity field obtained using the cross-correlation method has a much lower spatial resolution, data interpolation is used to generate a high-resolution velocity field with the same size as the original images. For the subpixel velocity components δu_x and δu_y in the image plane, they follow the equation ∂_tI + ∂_x (Iδu_x) + ∂_y (Iδu_y) = 0 [20]. Integration of this equation from 0 to 1 (one unit time) gives the intensity variation associated with the subpixel motion δI = –Δ_X(Iδu_X)–Δ_y(Iδu_Y) where Δ_x are Δ_y the differences over one pixel in the x and y coordinates in the image plane, and the spatial differences are Δx = Δy = 1 pixel. Therefore, the corrected mapping I(i,j) → I₀ + δI is used in the shifting scheme, where I₀ denotes the pixel-to-pixel shifted image obtained using the first procedure.

f. Pre-processing scheme

A pre-processing scheme is made to correct illumination intensity change. The laser illumination intensity may change between two successive PIV images acquired, which will affect optical flow computation. The illumination intensity change is contributed by an overall shift and a local change. The scheme has two steps. In the first step, the overall change in the whole images is corrected by shifting the averaged intensity difference between the two images.

The second step is to correct the local illumination intensity change. By applying a Gaussian filter with a suitable standard deviation to the two images, two sufficiently filtered (smoothed) images are obtained. The difference between the two filtered images represents the local illumination intensity change, and it is used to compensate the image intensity change by adding the difference to the original image to be corrected. The selection of the standard deviation σ of the Gaussian filter is a key, depending on the characteristic length scale of the local change pattern. When σ is too small, the optical flow could not be extracted correctly. On the other hand, when σ is too large, the local change could not be sufficiently corrected. Therefore, trial-and-error tests should be made to find a suitable value of σ for a particular application.

a. General description

The program of the hybrid method is written in Matlab, and the files including the main program, subroutines and sample images are contained in a folder named “OpenOpticalFlow_PIV.v1”. The main program is “OpenOpticalFlow_PIV_Run.m”. As shown in Figure 1, the main program has the following steps:

1. A pair of images (Image 1 and Image) is loaded;
2. The relevant parameters are set;
3. The images are pre-processed;
4. A coarse-grained displacement field is extracted from Image 1 and Image 2 by using a cross-correlation method;
5. A shifted image 1 is generated from Image 1 using the shifting scheme based on the coarse-grained displacement field that is interpreted at every pixel;
6. A residual displacement field is obtained from the shifted image 1 and Image 2 by using the optical flow method;
7. The corrected displacement field is obtained by superposing the interpreted coarse –grained displacement field and the residual displacement field;
8. The results including the velocity vectors, velocity magnitude, streamlines, vorticity and second invariant are plotted.

The coarse-grained displacement field for the initial approximation is obtained by using the subroutine “pivAnalyzeImagePair” that is adopted from the open-source Matlab software PIVlab [16, 17]. In optical flow computation, the subroutine for solving Eq. (3) is “liu_shen_estimator.m”, while the solution of the Horn-Schunck optical flow equation (“horn_schunck_estimator.m”) is used as an initial approximation for Eq. (3) for faster convergence. A combination of “liu_shen_estimator.m” and “horn_schunck_estimator.m” constitutes a key element of the optical flow program.

The processing can be made in a selected rectangular domain of interest for regional flow diagnostics when the index “index_region” is set at 1. Otherwise, the processing is conducted in the whole image domain when the index “index_region” is set at 0. The input files and relevant parameters for computation are listed in Table 1, which are set in the main program. The image files “Im1” and “Im2” could be in the conventional image format such as “tif”, “bmp”, or “jpeg”. The input parameters are discussed in the following subsections. The optical flow u = (u₁,u₂) is calculated by solving Eq. (3). Table 2 lists the outputs. The output data are “ux0” and uy0” that are the coarse-grained velocity components in the image plane (x, y), which are extracted by using the cross-correlation method. The refined velocity components “ux_corr” and “uy_corr” are obtained by using the optical flow method. The unit of the velocity given by the program is pixels/unit-time, where the unit time is the time interval between two input images. The velocity can be converted to the physical unit of m/s after the relation between the image plane and the 3D object space is established through camera calibration/orientation. The vorticity, strain rate and second invariant can be calculated based on a high-resolution velocity field. Furthermore, the main program can be adapted as a subroutine in a loop for processing of a sequence of image pairs such that the statistical quantities of the flow can be obtained.

b. The Lagrange multipliers

In optical flow computation, the Horn-Schunck estimator (“horn_schunck_estimator.m”) is used for an initial solution, and Liu-Shen estimator (“liu_shen_estimator.m”) is used for a refined solution of Eq. (1). In the main program, the Lagrange multipliers “lambda_1” and “lambda_2” are selected for the Horn-Schunck and Liu-Shen estimators, respectively. For example, “lambda_1” = 20 and “lambda_2” = 2000 for typical computations. Since no rigorous theory exists for determining the Lagrange multiplier, simulations based on a synthetic velocity field for a specific measurement are conducted to determine the Lagrange multiplier by using an optimization scheme [8]. The smallest Lagrange multiplier that still leads to a well-posed solution could be selected as a starting value in a trial-and-error process. Nevertheless, within a considerable range of the Lagrange multipliers, the solution is not significantly sensitive to its selection. Users can do some tests by adjusting the Lagrange multipliers for their applications.

c. Filtering

Pre-processing of images is sometimes required to remove the random noise by using a Gaussian filter, where the standard deviation (std) of a Gaussian filter is selected, depending on the noise level in a specific application (for example the std of a Gaussian filter is 4-6 pixels for images of 480×520 pixels). In the main program, the mask size of a Gaussian filter is given by the parameter called “size_filter”, and the std of the Gaussian filter is 0.6 of the mask size. In certain applications, the background noise in PIV images should be cleaned up by setting the cutting threshold value.

d. Correction for illumination change

In optical flow computation, it is assumed that the illumination light intensity is constant. However, the illumination intensity field could change in a time interval between two successively acquired PIV images. In this case, correction for this illumination intensity change is required before optical flow computation. This correction is done by calling the subroutine “correction_illumination”. As described in the previous section (f), the first step is the overall illumination change in the whole image is corrected by shifting an averaged intensity value in a selected domain. The parameter “window_shifting,size” defines the rectangular domain for calculating the averaged intensity value. The second step is the use of a simple scheme for correcting local illumination intensity changes based on the application of a Gaussian filter [11]. The selection of the std (or size) of a Gaussian filter is important to determine the local averaged intensity value for correction. This size of the Gaussian filter is given by the parameter called “size_average” in the main program. The selection of the filter size is a trial-and-error process depending on the pattern of illumination change in a specific measurement, and simulations on a synthetic velocity field are helpful to determine the suitable size.

e. Coarse-grained displacement field

The coarse-grained displacement field for the initial approximation is obtained by using the cross-correlation subroutine “pivAnalyzeImagePair”. The window sizes in the multiple-path process are given by the parameter “pivPar.iaSizeX”. For example, “pivPar.iaSizeX = [64 32 16]” defines the three-path process from 64×64 to 32×32 to 16×16 pixel² windows. The overlapping between the windows is defined by the parameter “pivPar.iaStepX”. For example, the 50% overlapping is defined by “pivPar.iaStepX = [32 16 8]” which corresponds to “pivPar.iaSizeX = [64 32 16]”. The parameter for the method of calculating the cross-correlation is “pivPar.ccMethod = ‘fft’”, indicating the FFT is used.

In principle, the displacement field can be successively improved by iterations to achieve a better accuracy. The iteration number for successive corrections is given by “no_iteration”. Since one iteration is required at least, “no_iteration” is set at 1 or a larger integer. Due to data interpretation, data near the left and top edges have large errors for certain images. In this case, these bad data are replaced by setting the data using the good data in the nearest interior region. The edge width is defined by the parameter “edge_width”. In most cases, this procedure is not required, and thus “edge_width” is usually set at 0.

a. Oseen vortex pair in uniform flow

An Oseen vortex pair in uniform flow has been used as a canonical case to evaluate the accuracy of the optical flow method [6, 13]. Here this flow is used as an example for application of the hybrid method. A pair of the particle image with the size of 500 × 500 pixels and the 8-bit dynamic range is generated for simulations, where 10000 particles with a Gaussian intensity distribution with the standard deviation of 2 pixels are uniformly distributed. Two Oseen vortices are placed at (m/3, n/2) and (2m/3, n/2) in an image, respectively, where m (500) and n (500) are the numbers of rows and columns of the image. The vortex strengths are Γ = ±7000 (pixel)²/s and the vortex core radius is r₀ = 15 pixels. The uniform flow velocity is 10 pixels/s. The second image is generated by using the shifting scheme based on the given velocity field after a time step. Figure 2 shows the generated image 1 and image 2, where the mean characteristic image diameter of particles is 4 pixels and the particle image density is 0.04 1/pixel².

The hybrid method is applied to Image 1 and Image 2. The cross-correlation method is applied with the three-path process from 32×32 to 16×16 to 8×8 and the 50% overlapping. In optical flow computation, the Lagrange multipliers in the Horn-Schunck estimator and the Liu-Shen estimator are set at 20 and 2000, respectively. Figure 3 shows the fields of the velocity vector, velocity magnitude and vorticity extracted from the particle images, where the streamlines are also shown and the velocity magnitude and vorticity field normalized by their maximum value.

The accuracy of the hybrid method has been studied through simulations of the Oseen vortex pair in in a uniform flow comparison with the optical flow method and cross-correlation methods [20]. As indicated before, the error depends on the four parameters: particle displacement, particle velocity gradient, particle image density, and particle image diameter. It is found that the hybrid method consistently performs better than the typical cross-correlation methods although the optical flow method is more accurate for PIV images with sufficiently small displacements (typically less than 3 pixels). To large extent, the hybrid method indeed combines the advantages of the cross-correlation method and the optical flow method for PIV measurements.

b. Vortex-induced separation

PIV images were obtained in a circular air jet normally impinging on a wall from a contoured circular nozzle in experiments [13]. Figure 4 shows a pair of PIV images focusing on the wall-jet region where strong interactions between vortices and boundary layer occur. Figure 5 shows the fields of the velocity vectors, velocity magnitude and vorticity obtained by using the hybrid method. It is revealed that large vortices generated in the shear layer of the wall jet due to the Kelvin-Helmholtz instability induce secondary vortices (secondary separations) in boundary layer near the wall. Compared to the cross-correlation method in PIV, the hybrid method yields the velocity field with much higher spatial resolution (one vector per a pixel), revealing more details of the near-wall flow structures in the thin wall jet (the thickness is about 3 mm in the physical scale, and 100 pixels in the image plane).

c. Twin jet

PIV images were obtained in twin jet excited acoustically where the ratio between the space between the two circular jets and exit diameters was 3 and the Reynolds number based on the diameter was 1800 [21]. Figure 6 shows a pair of particle images of acoustically excited twin jet. Figure 7 shows the fields of velocity vectors, velocity magnitude, and vorticity of acoustically excited twin jet extracted from the particle images.