1 OpenOpticalFlow : An Open Source Program for Extraction of Velocity Fields from Flow Visualization Images

This paper describes OpenOpticalFlow, an open source optical flow program in Matlab for extraction of high-resolution velocity fields from various flow visualization images. This program is a useful tool for researchers to use the optical flow method in various flow measurements. The principles of the optical flow method are concisely described, including the physics-based optical flow equation, the variational solution, and errors. The central part of this paper is the descriptions of the main program, relevant subroutines and selection of the relevant parameters in optical flow computation. Examples are given to demonstrate the applications of the optical flow method.


Introduction
In addition to particle image velocimetry (PIV) widely used for global velocity measurements, the optical flow method has been developed for extraction of high-resolution velocity fields from various images of continuous patterns from flow visualization images obtained in laboratories to cloud and ocean images taken by satellites/spacecraft [1][2][3][4][5][6].The rational foundation for application of the optical flow method to fluid flow measurements is the quantitative connection between the optical flow and the fluid flow velocity for various flow visualizations.Liu and Shen [7] have derived the projected motion equations for various flow visualizations including laser-sheet-induced fluorescence images, transmittance images of passive scalar transport, schlieren, shadowgraph and transmittance images of density-varying flows, transmittance and scattering images of particulate flows, and laser-sheet-illuminated particle images.Further, these equations are recast into the physics-based optical flow equation in the image plane, where the optical flow is proportional to the light-ray-path-averaged velocity of fluid (or particles) weighted in a relevant field quantity like dye concentration, fluid density or particle concentration.The optical flow method has been used to study the flow structures of Jupiter's Great Red Spot (GRS), impinging jets, and laser-induced underwater shock wave [8][9][10][11].A mathematical analysis of the variational solution of the optical flow and an iterative numerical algorithm are given by Wang et al. [12].The systematical error analysis of the optical flow method in velocity measurements is given by Liu et al. [11] in comparison with the crosscorrelation method in PIV. where is the velocity in the image plane referred to as the optical flow, g is the normalized image intensity that is proportional to the radiance received by a camera, is the Laplace operator, and is proportional to the pathaveraged velocity weighted with the field quantity  related to a visualizing medium.In the special case where 0 g    u and 0 f  , Eq. ( 1) is reduced to the Horn-Schunck brightness constraint equation [13].In general, the optical flow is not divergence-free, i.e., 0    u .

b. Variational solution
To determine the optical flow, a variational formulation with a smoothness constraint is typically used [7,13], which in fact is the first-order form of the Tikhonov's formulation for illposed problems.Given g and f , we define a functional where  is the Lagrange multiplier and  is an image domain.By minimizing the functional, i.e., min ) ( J  u , we obtain the Euler-Lagrange equation The standard finite difference method is used to solve Eq. ( 3) with the Neumann condition on the image domain boundary   for the optical flow [7, 12].c.Errors Liu et al. [11] gave an error analysis of the optical flow method.In general, the estimate for the total error of optical flow computation is given by where the relevant parameters are the displacement x  , the image intensity gradient g  , the velocity gradient u  and the velocity magnitude u , and 1 c and 2 c are coefficients to be determined.According to Eq. (4), it is indicated that the error is proportional to x  , and small displacements are generally required for a good accuracy in optical flow computation.The error would be larger in the regions where u  is large and/or g  is smaller.
For particle images in PIV, Eq. ( 4) can be further evaluated by including the effects of the particle image diameter and the particle image density.Therefore, the estimate for the total error of the optical flow method applied to PIV images is given by 6  c are coefficients to be determined, and m  is the mth-order difference operator.In Eq. ( 5), the main parameters are the particle displacement

Program and selection of parameters a. General description
The optical flow program is written in Matlab, and the files are contained in a folder named "OpenOpticalFlow.v1".The main program is "Flow_Diagnostics_Run.m", in which a pair of images is loaded, the relevant parameters are set, images are pre-processed, optical flow computation with a coarse-to-fine iteration is conducted, and the results are plotted.The subroutine for solving Eq. ( 3) is called "liu_shen_estimator.m".In computations, 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.The combination of "liu_shen_estimator.m" and "horn_schunck_estimator.m" constitutes the key elements of the subroutine "OpticalFlowPhysics_fun.m" used in the main program.
The optical flow processing can be done in a selected rectangular region 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 optical flow computation are listed in Table 1.The image files "Im1" and "Im2" could be in the conventional image format such as "tif", "bmp", "jpeg" etc.The input parameters are discussed in the following subsections.The optical flow  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 from the dowmsampled images in the coarse-to-fine process.The refined velocity components "ux_corr" and "uy_corr" are obtained by using the coarse-to-fine scheme.The unit of the velocity given by the program is pixels/unit-time, where the unit time is the time interval between two sequent 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, this program can be adapted for processing of a sequence of image pairs such that the statistical quantities of the flow can be obtained.

b. The Lagrange multipliers
The optical flow program has the Horn-Schunck estimator ("horn_schunck_estimator.m") for an initial solution and Liu-Shen estimator ("liu_shen_estimator.m") for a refined solution of Eq. ( 1).
There is no rigorous theory for determining the Lagrange multiplier in the variational formulation of the optical flow equation.The Lagrange multiplier acts like a diffusion coefficient in the corresponding Euler-Lagrange equations.Therefore, a larger Lagrange multiplier tends to smooth out finer flow structures.In general, the smallest Lagrange multiplier that still leads to a well-posed solution is selected.However, within a considerable range of the Lagrange multipliers, the solution is not significantly sensitive to its selection.Simulations based on a synthetic velocity field for a specific measurement could be carried out to determine the Lagrange multiplier by using an optimization scheme [8].

c. Filtering
Pre-processing of images is sometimes required to remove the random noise by using a Gaussian filter (the subroutine "pre_processing_a.m").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 "Flow_Diagnostics_Run.m", 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.

d. Correction for illumination change
The underlying assumption in optical flow computation is that the illumination light intensity keeps constant in flow visualization, which is valid in the well-controlled laboratory conditions.In some situations, however, the illumination intensity field could change in a time interval between two successively acquired images.For example, when the Jupiter's atmosphere structures were imaged by spacecraft, the illumination intensity field from the Sun could be changed in a relatively long time interval (hours) during image acquisition depending on the relative movement between the Sun, Jupiter and spacecraft.In this case, correction for this illumination intensity change is required before applying the optical flow method to these images.This correction is made by using the subroutine "correction_illumination.m".
In the first step in this subroutine, the overall illumination change in the whole image is corrected by normalizing both the images.In the second step, a simple scheme for correcting local illumination intensity changes is also based on 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 "Flow_Diagnostics_Run.m".
When the size of a filter is too large, the local variation associated with the illumination intensity change in images cannot be corrected.On the other hand, when the filter size is too small, the apparent motion of features in the two images would be artificially reduced after the procedure is applied.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 used to determine the suitable size.

e. Coarse-to-fine scheme
As indicated in Section 2.3, the optical flow method as a differential approach is more suitable for extraction of high-resolution velocity fields with small displacement vectors (typically less than 5 pixels depending on image patterns) from images of continuous patterns.
When the displacements in the image plane are large (for example more than 10 pixels), the error in optical flow computation would be large as.In this case, a coarse-to-fine iterative scheme can be used.This coarse-to-fine iterative scheme is implemented as a loop in the main program "Flow_Diagnostics_Run.m".
First, the original images are suitably downsampled by a suitable scale factor (such as 0.5) so that the displacements in pixels are small enough (1-5 pixels).The parameter called "scale_im" gives a scale factor for downsampling, and for example "scale_im" = 0.5 means that the original images are reduced to 50% of the original size.A coarse-grained velocity field is obtained by applying the optical flow algorithm to the downsampled images.
Then, the resulting coarse-grained velocity field is then used to generate a synthetic shifted image with the same spatial resolution as the original image #1 (i.e., the first one in the two successive images) by using an image-shifting (or image-warping) algorithm with an embedded spatial interpolation scheme.This subroutine "shift_image_fun_refine_1.m" uses a translation transformation for large displacements and the discretized optical flow equation for sub-pixel correction.
Next, the velocity difference field between the synthetically shifted image and the original image #2 is determined by using the optical flow algorithm, and it is added on the initial velocity field for correction or improvement.Thus, a refined velocity field is successively recovered by iterations to achieve a better accuracy.This process can be repeated iteratively.The iteration number is given by "no_iteration".Usually, one or two iteration is sufficient.When "no_iteration" is set at 0, the coarse-to-fine scheme is disabled particularly when the displacements in the original images are sufficiently small.

Examples a. Oseen vortex pair in uniform flow
To demonstrate the application of this program, simulations are conducted on particle images (PIV images) in a synthetic flow  an Oseen vortex pair in uniform flow.A sample particle image with the size of 500 500  pixels and the 8-bit dynamic range is generated, where 10000 particles with a Gaussian intensity distribution with the standard deviation of  = 2 pixels are uniformly distributed.Figure 1 shows a pair of the generated particle image that is used for simulations, where the mean characteristic image diameter of particles is The second image is generated by deforming the original synthetic particle image based on the given velocity field after a time step t  .This processing is made by applying an image-shifting algorithm, which faithfully describes the motion of image patterns for a given velocity field since the physics-based optical flow equation is derived from the governing transport equations for flow visualizations.The time step t  is used to control the maximum displacement in synthetic particle images.This flow with particle images has been used as a typical case in the simulations to evaluate the accuracy of the optical flow method [11].
In optical flow computation, the Lagrange multiplier in the Horn-Schunck estimator is set at 20, which provides a good initial solution for refinement by the Liu-Shen estimator based on Eq. ( 1).The Lagrange multiplier in the Liu-Shen estimator is fixed at 2000 for a refined velocity field, and it does not significantly affect the velocity profile in a range of 1000-20000 except the peak velocity near the vortex cores in this flow.A pair of particle images with 10000 particles in Fig. 1 is processed as a typical case, where the maximum displacement is 2.6 pixels.Figures 2(a) and 2(b) show the coarse-grained and refined velocity vector fields extracted from the particle images, respectively.The optical flow method gives 500 500  vectors from the original images.Figure 3 shows the corresponding streamlines.Figure 4 shows the refined vorticity field normalized by its maximum value.

b. White Ovals on Jupiter
The White Ovals are distinct storms in the Jupiter's atmosphere.PIV images were obtained in an air jet normally impinging on a wall from a contoured circular nozzle in experiments, focusing on the wall-jet region where strong interactions between vortices and boundary layer occur [11].Figure 10(a) shows one sample PIV image of the walljet region with relatively high particle density.Figure 10(b) shows the velocity vectors overlaid on the normalized vorticity field.The optical flow method is able to provide the velocity field near the wall at the resolution of one vector per one pixel.The optical flow is able to reveals large vortices generated in the shear layer of the free jet and wall jet due to the Kelvin-Helmholtz instability, and induced secondary vortices (secondary separations) in boundary layer near the wall.Compared to the traditional cross-correlation method in PIV, the optical flow 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). (

2) Availability
The folder called "OpenOpticalFlow_v1" is available for downloading at a website of Western Michigan University, which contains a set of Matlab programs, sample images and a document describing the optical flow method and the relevant programs.The program can be run in Matlab (R2007a, and newer versions) in Windows on a PC.Several functions in Matlab image processing toolbox are required.Questions on the use of this open source program could be sent to the email address of the author for discussion.

(3) Reuse Potential
Global velocity diagnostics is of fundamental importance in the study of fluid mechanics in order to understand the physics of complex flows.In addition to particle image velocimetry (PIV), the optical flow method allows extraction of high-resolution velocity fields from flow visualization images obtained in various measurements in different scientific and engineering fields.These images could be obtained by using inexpensive CCD and CMOS cameras under various illumination conditions.The physical and mathematical foundations of the optical flow method are well established.Optical flow computation can be conducted in a PC with Matlab to extract velocity fields at one vector per a pixel.The optical flow method has been used in some challenging flow measurements in which the traditional cross-correlation method in PIV could not provide sufficient high-resolution velocity fields in large domains [7][8][9][10][11][12] The objective of this paper is to describe an open source optical flow program in Matlab for extraction of high-resolution velocity fields from various flow visualization images.It would help researchers to apply the optical flow method to their specific applications in relevant scientific and engineering fields.First, the principles of the optical flow method are briefly described, including the physics-based optical flow equation, its variational solution, and a short error analysis.The physical meaning of the optical flow in fluid-mechanic measurements is interpreted.Then, the main program and the key subroutines are described, particularly on the selection of the relevant parameters in optical flow computation.Next, to demonstrate the applications of this program, several examples based on simulated and real flow visualization images are presented.Principles of optical flow method a. Physics-based optical flow equation Flow visualization techniques use tracers (such as particles and dyes) or on the change of certain physical properties of fluid (such as the density).Digital images of flow visualization are obtained by using cameras that detect radiation with a certain wavelength bandwidth from a fluid medium in a flow.Liu and Shen [7] have derived the physics-based optical flow equation for different flow visualizations, which is written in in the image coordinates as


should be small in optical flow computation using particle images, particularly in regions of large velocity gradients.The particle image diameter p d could have the optimal value for optical flow computation.The particle image density p N should be suitably large, and there would be the optimal value of p N .

p d = 4
pixels and the particle image density is 0.04 1/pixel 2 .A synthetic velocity field, generated by superposing an Oseen vortex pair on a uniform flow, is used as a canonical flow.Two Oseen vortices are placed at in an image, respectively, where m (500) and n (500) are the numbers of rows and columns of the image.The circumferential velocity of an Oseen vortex 2 /s and the vortex core radius is 15 r 0  pixels.The uniform flow velocity is 10 pixels/s.

Figure 5
shows two successive near-infrared continuum filtered (756 nm) images of the White Ovals.Three sets of images were taken on 02/19/1997 at a range of 1.1 million kilometers by the Solid State Imaging (CCD) system aboard NASA's Galileo spacecraft.Each is taken one hour apart.Unlike discrete particle images, these images have continuous patterns that are particularly suitable for the application of the optical flow method.It is noticed that the illumination from the Sun was considerably changed in a local and non-uniform fashion.Thus, the correction for this change is made by using the subroutine "correction_illumination.m" before applying the optical flow method.Figures6(a) and 6(b) show the coarse-grained and refined velocity vector fields extracted from the White Ovals images, respectively.Figure7shows the corresponding streamlines.Figure8shows the refined vorticity field normalized by its maximum value.It is indicated the interactions between the three cyclonic vortices.The balloon-shaped vortex is seen between the two well-formed ovals.In optical flow computation in this case and the following cases, the Lagrange multipliers in the Horn-Schunck and Liu-Shen estimators are set at 20 and 2000, respectively.c.Other examples Quasi 2D turbulence is simulated by randomly distributing 50 Oseen vortices in the uniform distribution in a cloud image.The distribution of the strengths of the vortices is the Gaussian distribution with zero mean and the standard deviation of 3000   (pixel) 2 /s.The sizes of the vortices are distributed in a Gaussian distribution with the mean of 20 r 0  pixels and the standard deviation of 10 pixels.One of the pair of the cloud images is shown in Fig. 9(a).
Figure 9(b) shows the velocity vectors overlaid on the normalized vorticity field extracted by the optical flow algorithm from the cloud image pair.The optical flow method can give better estimation in the distribution of velocity particularly in the regions with rapid changes in comparison with the traditional cross-correlation method applied to such images.The mean error in the whole image is about 0.03 pixels/unit-time in this case.

Figure 4 .Figure 5 .Figure 6 .Figure 7 .
Figure 4. Vorticity field of an Oseen vortex pair in uniform flow extracted from the particle images.