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.

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 [

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 problems in relevant scientific and engineering fields. First, the principles of the optical flow method are briefly described, including the physics-based optical flow equation, the variational method, 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.

Flow visualization techniques are widely used based on the applications of observable tracers (such as particles and dyes) and change of certain physical properties of fluid (such as the density). Digital flow visualization images are usually obtained by using cameras detecting radiation with a certain wavelength bandwidth from a fluid medium in a flow. Liu and Shen [

where _{1}, _{2}_{β}^{2} = ∂^{2}/∂_{β}_{β}_{1}, _{2}, _{1}, _{2}

To determine the optical flow, a variational formulation with a smoothness constraint is typically used [

where

The standard finite difference method is used to solve Eq. (3) with the Neumann condition ∂

An error analysis of the optical flow method was given by Liu et al. [

Where the relevant parameters are the displacement magnitude ‖Δ_{1} and _{2} are coefficients to be determined. Since according to Eq. (4) the total error is proportional to ‖Δ_{p}

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:

where _{1}, _{2}, _{3} and _{4} are coefficients to be determined, and Δ^{m} is the _{p}_{p}_{p}_{p}_{p}_{p}_{p}_{p}

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 the main program, 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 “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 optical flow program.

The optical flow 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 optical flow computation are listed in Table _{1}, _{2}

Inputs.

Input files and parameters | Notation | Unit | Note |
---|---|---|---|

image pair | Im1, Im2 | – | tif, bmp, jpeg etc. |

Lagrange multipliers | lambda_1 for the Horn-Schunck estimator |
– | regularization parameter in variational solution |

Gaussian filter size | size_filter | pixel | removing random image noise |

Gaussian filter size | size_average | pixel | correction for local illumination intensity change |

scale factor for downsampling of original images | scale_im | – | reduction of initial image size in coarse-to-fine scheme |

number of iterations | no_iteration | – | iteration in coarse-to-fine scheme |

indicator for regional diagnostics (0 or 1) | index_region | – | “0” for whole image; |

Outputs.

Output files | Notation | Unit | Note |
---|---|---|---|

coarse-grained velocity | ux0, uy0 | pixels/unit-time | based on downsampled images |

refined velocity | ux_corr, uy_corr | pixels/unit-time | refined result with full spatial resolution in coarse-to-fine process |

The optical flow program uses 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). 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.

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 could be selected by a trial-and-error process. 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 [

Pre-processing of images is sometimes required to remove the random noise by using a Gaussian filter. 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.

The underlying assumption in optical flow computation is that the illumination light intensity keeps constant (or time-independent) 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 change 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.

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 the application of a Gaussian filter [

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 patterns/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.

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 from images of continuous patterns (typically less than 5 pixels depending on the size of image patterns). When the displacements in the image plane are large (for example more than 10 pixels), the error in optical flow computation would be significant. In this case, a coarse-to-fine iterative scheme can be used, which is implemented in a loop in the main program.

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). A scale factor is specified by the parameter called “scale_im” for downsampling. 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 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 in which a translation transformation is used for large displacements and the discretized optical flow equation is used 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. 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.

To demonstrate the application of this program, simulations are conducted on particle images (PIV images) in a synthetic flow 3:4 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 _{p}^{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 ^{2}/s and the vortex core radius is _{0} = 15 pixels. The uniform flow velocity is 10 pixels/s. The second image is generated by deforming the original synthetic particle image based on the given velocity field after a time step Δ

A pair of particle images of an Oseen vortex pair in a uniform flow,

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 Figure

Velocity vectors of an Oseen vortex pair in a uniform flow extracted from the particle images,

Streamlines of an Oseen vortex pair in uniform flow extracted from the particle images,

Vorticity field of an Oseen vortex pair in uniform flow extracted from the particle images.

The accuracy of the optical flow method applied to PIV images has been studied through simulations and experiments in comparison with the well-established cross-correlation method [

The White Ovals are distinct storms in the Jupiter’s atmosphere. Figure

A pair of images of the White Ovals on Jupiter,

Velocity vectors of the White Ovals on Jupiter,

Streamlines of the White Ovals on Jupiter,

Vorticity field of the White Ovals on Jupiter.

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 ^{2}/s. The sizes of the vortices are distributed in a Gaussian distribution with the mean of _{0} = 20 pixels and the standard deviation of 10 pixels. One of the pair of the cloud images is shown in Figure

Quasi 2D turbulence,

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 [

Wall-jet region of an impinging jet,

Based on Matlab (R2007a, or later versions): Windows.

Matlab (R2007a, or later versions).

None.

Several functions in Matlab image processing toolbox are required.

Tianshu Liu, Department of Mechanical and Aerospace Engineering, Western Michigan University.

Lixin Shen, Department of Mathematics, Syracuse University.

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 a spatial resolution of 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 [

The author has no competing interests to declare.