An Optical Flow Algorithm with Automatic Parameter Adjustment for Fluid Velocimetry

1.1 Introduction

Accurate quantitative measurement of fluid motion is a fundamental problem in experimental fluid mechanics. One standard experimental method for fluid velocimetry is particle image velocimetry (PIV), which entails inferring fluid motion from the displacement of the particles in a seeded flow. A comprehensive treatment of the subject can be found in texts such as the one by Raffel et al. [28]. Briefly, small particles are seeded into a flow, and are illuminated using a laser sheet. The Mie scattering signal produced by the laser-particle interaction is captured using a digital camera, and the displacement of particles in two consecutive images is used to infer the velocity of the fluid. A widely popular way of determining particle displacements from PIV images is cross-correlation, which determines the particle displacement vectors by cross-correlating sub-regions of the images called interrogation windows (IW’s) of two consecutive particle image frames [36, 28]. However, cross-correlation-based PIV processing struggles to resolve sharp, small-scale velocity gradients accurately, and produces low-pass filtered velocity fields due to inherent spatial averaging over finite-sized interrogation windows [24, 30, 22, 12, 18, 13].

An alternative approach to cross-correlation-based PIV processing is optical flow velocimetry (OFV), which computes high spatial resolution velocity fields from the same particle image data. OFV has been shown to provide higher accuracy when compared to cross-correlation-based processing under certain conditions (i.e. high particle image density [32, 31]), and is summarized in the recent review paper by Jassal and Schmidt [19]. Mathematically, OFV inverts the displaced frame difference equation given in Eq. (1) for the displacement field $\mathrm{\Delta }\underset{¯}{x}=\underset{¯}{u}\mathrm{\Delta }t$, where I₀ and I₁ are two particle image pairs separated by some time interval Δt, and $\underset{¯}{u}$ is the unknown velocity field. The direct inversion of the DFD equation is an ill-posed inverse problem due to having only one equation per pixel, the intensity conservation from Eq. (1), but two unknowns, i.e. the two components of the velocity field. A typical way to remedy the ill-posedness of an under-determined problem is to apply regularization, which can be explicit, implicit, or both [34]. Such an approach is called a variational solution, and the interested reader is referred to the review article by Heitz et al. for a more complete discussion [10].

In a variational approach, the velocity field is determined by minimizing a functional based on Eq. (1). The functional used in this work to invert the DFD equation is given in Eq. (2).

To introduce flow physics into the displacement field estimated, a fluid physics-based regularization strategy is typically introduced to dampen the non-physical gradients that may emerge in the solution due to “over-fitting” to the DFD equation. The software presented in this work uses the viscosity-like regularization proposed by Schmidt and Sutton [34] given in Eq. (3). The role of the regularization term is to encourage smoothness of the velocity field while preserving vortical structures.

By combining the regularization term in Eq. (3) with the variational formulation of the DFD equation in Eq. (2), one can write the optical flow problem as a minimization problem, which aims to minimize the functional given in Eq. (4) for the velocity field $\underset{¯}{\overset{^}{u}}$. In this work, instead of determining the velocity field directly, its wavelet representation $\underset{¯}{\psi }$ is computed. This is done in order to remedy the ill-posedness of the problem by reducing the dimensionality of the unknowns by truncating the wavelet coefficients, hence achieving implicit regularization. The reader is referred to Refs. [8, 9, 34] for further details.

However, the optimization problem proposed in Eq. (4) has a distribution of user-tuned parameters $\theta \left(\underset{¯}{x}\right)$. This locally varying distribution specifies the relative weights of the regularization and the data terms on a per-pixel basis. Tuning this distribution manually is impractical as it would entail specifying over a million parameters for a one-megapixel image. The specification of parameters imposes a serious challenge when the ground truth is unknown, as in the case of experimental data, therefore a mathematical model is proposed to determine this distribution automatically. Briefly, it is assumed that $\theta \left(\underset{¯}{x}\right)$ are derived from a Γ-distribution, and more details on the subject can be found in Refs. [2, 3]. This results in a new optimization problem given in Eq. (5), which determines the regularization distribution θ along with the velocity field $\underset{¯}{u}$. Intuitively, the assumption about $\theta (\underset{¯}{x})$ belonging to a Γ-distribution is an additional constraint on the minimization problem in Eq. 4, and the goal becomes to determine the shape of this Γ-distribution instead of $\theta (\underset{¯}{x})$. The reader is referred to Jassal et al. [21, 12] for further details.

The new optimization problem given in Eq. (5) has only one global parameter which must be specified, η, which essentially decides the shape of the Γ-distribution that governs the spatial behavior of $\theta (\underset{¯}{x})$. In this work, it is proposed that η can be estimated automatically from the cross-correlation results using Eq. (6) as described in Sec. 1.3. Briefly, a coarse estimate of the velocity field is obtained using cross-correlation with the final vector spacing being ΔI pixels. The coarse velocity field is then linearly interpolated to the full image pixel grid, and is used to compute J_D and J_R, which in turn are used to obtain an estimate of η using Eq. (6). Details regarding the optimization problem posed can be found in Jassal et al. [21, 12].

The software described here solves the minimization problem given in Eq. (5) to obtain the velocity field from a pair of particle images. The solution to minimization problem is computed in the wavelet domain and is based on the work of Schmidt and Sutton [32, 33]. It should be noted that the presented approach is not a hybrid-optical flow approach. A hybrid-optical flow method typically uses the cross-correlation velocity field as the initial velocity estimate, and refines it further using an optical flow algorithm [19]. Such a hybrid approach is prone to biases from the initial guess. The presented method computes the velocity field independently and does not make use of the cross-correlation velocity field as an initial guess–it is only used to tune the parameter η.

1.2 Testing and comparison

The presented software is tested on synthetic data where the ground truth is known, and therefore the error can be accurately computed. For the first test case, three-dimensional homogeneous isotropic turbulence (3D-HIT) dataset described in Jassal and Schmidt [17] is used. Figure 1 shows the magnitude of the estimated velocity field using (b) wOFV^⋆ (the method presented) and (c) cross-correlation at an instant. For cross-correlation processing, a final interrogation window of size 32×32 with a 50% overlap was used. The cross-correlation method makes use of iterative window-deformation and outlier removal. It is observed that wOFV^⋆ is able to better resolve the small-scale turbulent structures compared to cross-correlation. For the velocity fields presented, the root mean square error (RMSE) normalized by the root mean square velocity is found to be 0.2079 and 0.2151 for wOFV^⋆ and cross-correlation respectively. Therefore, wOFV^⋆ not only leads to a gain in resolution but also increases the accuracy of the velocity estimates.

Figure 1

This gain in accuracy is observed for the whole dataset as evident from the results shown in Figure 2, which presents the normalized RMSE for (a) velocity fields and (b) vorticity fields estimated using wOFV^⋆ and cross-correlation for 999 image pairs for the entire 3D-HIT dataset. The mean velocity RMSE is found to be 0.194 for wOFV^⋆ and 0.2217 for wOFV^⋆ and cross-correlation respectively, leading to a 15% accuracy gain for the presented method. A similar gain in accuracy is observed for the vorticity fields as well, which requires an accurate estimate of the velocity gradients. The mean vorticity RMSE for wOFV^⋆ is found to be 0.0841, which is 16% lower than that of cross-correlation (0.0971). The presented software has been tested on multiple different datasets and further details on its accuracy and sensitivity can be found in Ref. [12] and Jassal et al. [21].

Figure 2

The software is also compared to the other open source optical flow methods, and it has been demonstrated to produce high-resolution, accurate velocity fields. Figure 3 shows the RMSE computed using the velocity fields obtained using wOFV^⋆ and other codes by Dérian [7], Corpetti et al. [5], Liu and Shen [26], and cross-correlation [36], for the widely used synthetic PIV benchmark test case of a direct numerical simulation (DNS) of 2D HIT flow produced by Carlier and Wieneke [4]. Liu and Shen’s optical flow method used can downloaded from https://github.com/Tianshu-Liu/OpenOpticalFlow˙PIV˙v1. For this method the following parameters were used: λ_HS = 0.7, λ_LS = 260, NumIter = 1, DownSampleScale = 0.5, AverageSize = 20, and GaussianFilterSize = 4. Dérian’s method can be accessed at https://allgo.inria.fr/app/typhoon. For processing data with Dérian’s method, 25% of the wavelet coefficients were retained. For cross-correlation processing, a final interrogation window of size 32×32 with 50% overlap was used. The cross-correlation method makes use of iterative window-deformation and outlier removal to enhance accuracy. It is evident that the presented software provides much more accurate velocity estimates compared to the other available contemporary methods, while not requiring any manual parameter tuning. However, the method (like all OFV methods [19]) needs high-quality particle images that contain lower noise and higher seeding density than what is acceptable cross-correlation images to produce high-resolution results. Obtaining such images could be a limitation when high-seeding density is not achievable in an experiment. For further details, the reader is referred to Refs. [12, 31].

Figure 3

1.3 Implementation and architecture

The method described in Sec. 1.1 is implemented as an additional processing module in an open-source, cross-correlation-based PIV processing code, PIVlab [36]. Hence it uses the same GUI interface of PIVlab, along with various other pre- and post-processing options. The core workflow of the software is described in Figure 4, and further details are provided in the subsequent discussion in Sec. 1.3.3. Briefly, the images are loaded using PIVlab and are passed to the RunMain function, which segments the images into smaller regions or “patches” if required, and passes each region to BayesPyd_sub. The optimization problem given in Eq. (5) is then solved iteratively for each image or patch to determine the velocity field, which is stitched together by RunMain to obtain the full flow field, that is passed back to PIVlab as an output. After the images are loaded into PIVlab and optical flow is selected as the processing module, the user is allowed change the four options described in Sec. 1.3.1.

Figure 4

1.3.1 End-user options

1.3.1.1 Patch Size (M)

As the algorithm is based on solving the optimization problem in the wavelet domain, it is necessary to work with square images, since a separable wavelet transform requires a square domain of size 2ⁿ, where n∈ℤ⁺. However, as PIV images can be of an arbitrary size and aspect ratio, the first step is to segment the image into smaller square sub-regions, referred to here as patches, as shown in Figure 5. The user can choose the size of these patches if desired from a drop down list offering ${128}^2,{256}^2,{512}^2$, and 1024² as the square patch size. By default, the largest patch size that can fit into the image domain is chosen. However, if parallel processing is enabled, the user can increase the computational efficiency by choosing a more computationally efficient domain segmentation. This can be done by using the suggest settings option, which suggests a computationally optimal patch size based on the largest displacement in the flow field using coarse cross-correlation results. Once the patch size is selected, the image domain is segmented such that there is a minimum of 15% spatial overlap between adjacent patches.

Figure 5

To mitigate any effects that may arise due to particle (data) loss at the edges of the patches and other edge artifacts, a two-dimensional Nuttall [27] weighting window is used when stitching together the full flow field from individual patches. Therefore, the velocity field near the center of each patch is weighted more strongly relative to the edges when combining velocities from various patches in an overlapping region. However, errors can still manifest at the edges of patches in the presence of large particle displacements in the flow field if the patch size is too small relative to the displacements. This is a consequence of a large number of particles crossing patch boundaries over the time interval between frames. Figure 6 illustrates this effect by showing the velocity magnitude estimated from the same synthetic particle image data processed with three different patch sizes, namely 64², 128², and 256² pixels. The original image is 256×256 pixels, and the flow has a maximum particle displacement of approximately 12 pixels. As evident, at a patch sizes of 64² and 128² pixels, the patch boundaries are very pronounced in the computed velocity field, and there are significant non-physical errors that are not observed for a patch size of 256² pixels. A general recommendation is to choose a patch size that is about 10–12 times larger than the maximum displacements in the flow field.

Figure 6

1.3.1.2 Pyramid Level (j)

Typically, optical flow algorithms can determine displacements accurately only when the displacements are at most 2–3 pixels for an image pair. However, in an experimental setting, PIV images commonly have displacements as large as 10–15 pixels. Therefore, to resolve large displacements, an image pyramid scheme originally proposed by Ruhnau et al. [29] is implemented at a patch level along with the wavelet multi-resolution framework. The pyramidal scheme is depicted in Figure 7. Each image patch is successively Gaussian-filtered and down-sampled j times, and OFV is performed on the coarsest level first. Note that each down-sampling reduces the displacement magnitudes by a factor of two. The result is then translated to the next finest level using bicubic interpolation and used to perform an initial distortion of the particle images at that level before performing OFV again, and so on until the resolution level of the original images is reached. Pyramidal decomposition not only makes OFV algorithms more robust, but also increases computational efficiency because the computation at the finest level, which is the most computationally intensive, is initiated close to the final solution from the previous pyramid level.

Figure 7

The efficacy of the image pyramiding scheme is evident from the results presented in Figure 8, which presents the magnitude of the velocity computed from the same PIV data processed at different final pyramid levels of j = 1 and j = 2, compared to the ground truth from the DNS. Clearly, the algorithm is unable to resolve the flow field for j = 1, that is when the image is not down-sampled. However, once the image is down-sampled, and the pyramid level j is set to 2, the algorithm successfully resolves the flow field. By default, as generally advised, j is set to 3, to increase robustness as well as computational efficiency.

Figure 8

1.3.1.3 Smoothness (η)

As evident from Eq. (5), the optimization problem requires the specification of a parameter η. This parameter can be determined automatically by using the suggest settings button, and this is the recommended choice in practice. However, the user can modify the parameter if desired. By default, the smoothness parameter is set to 40, and is typically between 0 to 100. The details regarding the automatic estimate of the smoothness parameter are given in Sec. 1.3.2.5. The effect of the value of η is demonstrated in Figure 9, which shows the vorticity field computed from the velocity field using the same image pair from synthetic DNS data with different values of η. Generally speaking, a smaller value of η allows for a large amount of high frequency variations in the flow, whereas a larger value of η typically leads to a smoother flow field.

Figure 9

1.3.1.4 Inter-pass median filtering (Δ_M)

In the case of noisy data with large displacements, the pyramidal scheme described in Sec. 1.3.1.2 may amplify the noise in the particle images [1, 29] when a large number of pyramids are used. Therefore, to further increase the robustness of the algorithm, an inter-pass median filter is implemented to mitigate the effect of spurious vectors resulting from image pyramiding on the final computed velocity field. The user may choose an inter-pass median filter of size 3×3, 5×5, 7×7, or 9×9 pixels, if desired, when processing noisy data that requires using a large number of down-sampling to resolve large displacements. By default, the inter-pass median filter is disabled, as it is typically not required.

1.3.2 Sub-routines

The algorithm uses multiple sub-routines to streamline the code and to reuse the same operation multiple times. The key sub-routines are described in this section.

1.3.2.1 Psitrunc_mat

Psitrunc_mat takes a square matrix A of size m×m along with a wavelet filter matrix F of the same size, to return the wavelet transform $\tilde{A}$ of matrix A. The function performs a triple matrix product written in Eq. 7.

7

$\tilde{A}=FA{F}^T$

The same function can be used to obtain the original matrix A form its wavelet transform $\tilde{A}$, by passing the inverse wavelet transform filter matrix $\tilde{F}$ instead of F, and $\tilde{A}$ instead of A. This routine is used multiple times to transform the velocity field components between the physical and the wavelet space.

1.3.2.2 warpS

As the problem being solved by the algorithm is an inverse problem, it necessitates an accurate model of the forward problem. The forward problem is essentially an image warping operation where a particle image is deformed by the two-dimensional velocity field, and it is performed by the sub-routine warpS in the code. The inputs to the routine are the gridded interpolant object formed using the image, and a vector displacement field. The algorithm uses a symmetric boundary condition at the image edges, along with 2D-spline interpolation to obtain an accurate representation of the deformed image. Figure 10 illustrates the warping operation. The transformation is non-rigid, i.e. the shapes within an image can be bent, stretched, or distorted. Therefore, the forward model is capable of handling complex motion fields, such as the ones prevalent in fluid flows.

Figure 10

1.3.2.3 Reg

The sub-routine Reg computes the gradient of J_R given in Eq. (3) with respect to the wavelet coefficients $\underset{¯}{\psi }$, i.e. the wavelet transform of the velocity field. For computational efficiency, the operation is performed directly in the wavelet domain via a series of matrix multiplications as described by Schmidt and Sutton [33]. By default, the algorithm uses the viscosity-based regularization proposed by Schmidt and Sutton [34], however other regularization schemes such as Horn and Schunck [11], div-curl [6], and Laplacian [35] are also implemented, and can be selected by changing the regscheme input flag for the sub-routine. Once the gradient is obtained, the magnitude J_R is obtained readily via the operation given in Eq. (8).

8

$J_R=\frac{1}{2}<\underset{¯}{\psi },\underset{¯}{{\nabla }_{\psi }{J}_R}>$

The operation $<\underset{¯}{\psi },\underset{¯}{{\nabla }_{\psi }{J}_R}>$ denotes the inner product between the the wavelet coefficients $\underset{¯}{\psi }$ and gradient of J_R with respect to the wavelet coefficients.

1.3.2.4 GenerateImPyramid

As discussed in Sec. 1.3.1.2, the user can choose to enable the pyramidal decomposition scheme, which can help in resolving large displacements using OFV. Pyramidal decomposition is achieved using the GenerateImPyramid sub-routine, which takes an image and the number of pyramids (j) as the input. To ensure spatial differentiability of the images, a 3×3-pixel Gaussian filter is used to smoothen the images at each pyramid level. The Gaussian filtering plays a major role in increasing the accuracy of the computations. Figure 11 shows the magnitude of the velocity field computed from the same PIV data processed at the same settings, both with (c) and without (b) a Gaussian filter applied to the particle images at each pyramid level, compared to the ground truth (a). It is observed that certain regions in the flow field can are negatively affected in the absence of Gaussian filtering of the images at each pyramid level, as the differentiability of the images can be sacrificed leading to erroneous gradients which mislead the optimization. For computational efficiency, the image pyramids are generated at each patch, which allows independent computations for each patch, hence enabling parallelization of the code.

Figure 11

1.3.2.5 PredictSmoothnessCoefficient

One of the core features of the software developed is that all of the parameters can be estimated automatically, without requiring any user input. The smoothness parameter η can be estimated using the suggest settings button. Once the button is pressed, PIVlab runs a coarse velocity estimation via cross-correlation on the images and provides the velocity estimate to the PredictSmoothnessCoefficient routine. A schematic of estimating η is shown in Figure 12. The cross-correlation results are first interpolated to the full image grid using bilinear interpolation to obtain ${\underset{¯}{u}}_{\mathit{\text{CC}}}$. This linearly interpolated velocity field is used to warp the I₀ image forward and I₁ image backward by a displacement ${\underset{¯}{u}}_{\mathit{\text{CC}}}/2$ to compute $I_0^F$ and $I_1^B$ respectively. An estimate of the data term, $J_D^{\mathit{\text{CC}}}$, is computed using Eq. (9).

9

$J_D^{\mathit{\text{CC}}}=\frac{1}{2}{\int }_{\mathrm{\Omega }}{\left[I_0^F–I_1^B\right]}^2\text{d}\mathrm{\Omega }.$

Figure 12

In the same spirit, ${\underset{¯}{u}}_{\mathit{\text{CC}}}$ is transformed to the wavelet space to estimate $J_R^{\mathit{\text{CC}}}$ using Eq. (3). Combining $J_D^{\mathit{\text{CC}}}$, $J_R^{\mathit{\text{CC}}}$, and the vector spacing in the cross-correlation relations ΔI, η can be readily estimated using Eq. (10), where 𝒪 represents the order of magnitude operation.

10

${\log }_{10}\eta =\mathcal{O}(J_D^{\mathit{\text{CC}}})–\mathcal{O}(J_R^{\mathit{\text{CC}}})–{\log }_2\mathrm{\Delta }I$

1.3.3 Algorithm

In this section, the full code workflow is described. The core steps are illustrated in Figure 4, and the details are discussed below. The computations are essentially handled at three levels. At the first level, i.e. the RunMain routine, the image is segmented and combined as necessary, while also choosing the processing options. The velocity field for each patch is computed by the routine BayesPyd_sub, which optimizes the caller function call_Bayes using a gradient-based optimizer making use of the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm [25]. A pseudo-code for RunMain is given in Code 1.

For each patch, the velocity computation occurs in the BayesPyd_Sub routine, and a pseudo-code is presented in Code 2. The optimizer options are set as follows.

• DERIVATIVECHECK=’off’: a derivative check with finite differencing is turned off for computational efficiency.

• METHOD=’lbfgs’: a low memory, gradient based optimizer is used. See Ref. [25] for details.

• MAXITER=100: the maximum number of iterations is capped to a reasonable value. Typically, the solution converges in less than 100 iterations.

• MAXFUNEVALS=300: the maximum number of function evaluations per iteration is capped to 500.

• Corr=8: a history of 8 steps is stored to decide the convergence of the solution.

• progTol=1e-5 and optTol=1e-5: converge tolerances are set as high as possible to obtain an accurate solution while optimizing for computational efficiency.

Briefly, each patch is processed at different pyramid levels, and a multi-resolution framework based on wavelets is implemented at each pyramid level. The functional and its gradient are provided to the call_Bayes routine, which calls Bayes_min to compute J_D, J_R, and their respective gradients with respect to the wavelet coefficients $\underset{¯}{\psi }$.

As the minimization problem posed in Eq. (5) involves optimizing both $\underset{¯}{u}$ and θ, an Iterative Alternating Sequential (IAS) algorithm for MAP estimation [2, 3] is used to split the optimization problem into two separate sequential minimization problems that can be solved iteratively. The energy functional for each of the problems are defined as follows.

11

${\mathcal{E}}_{\theta }\left(\underset{¯}{\psi }\right)={\displaystyle \frac{1}{2}}{\displaystyle {\int }_{\mathrm{\Omega }}}{\displaystyle \frac{{\left(I_0–I_1\right)}^2}{\theta }}+{\displaystyle \frac{1}{2}}{\displaystyle {\int }_{\mathrm{\Omega }}}{\left[{\nabla }^2\underset{¯}{u}+{\displaystyle \frac{1}{3}}\underset{¯}{\nabla }\left(\underset{¯}{\nabla }\cdot \underset{¯}{u}\right)\right]}^2,$

12

${\mathcal{E}}_{\underset{¯}{\psi }}(\theta )={\displaystyle \frac{1}{2}}{\displaystyle {\int }_{\mathrm{\Omega }}}{\displaystyle \frac{{\left(I_0–I_1\right)}^2}{\theta }}–\eta {\displaystyle \sum _{j=1}^n}\log {\theta }_j+{\displaystyle \sum _{j=1}^n}{\theta }_j$

The code first minimizes ${\mathcal{E}}_{\theta }\left(\underset{¯}{\psi }\right)$ at fixed regularization field θ, and then minimizes ${\mathcal{E}}_{\underset{¯}{\psi }}(\theta )$ at fixed $\underset{¯}{\psi }$ to update θ. The process is repeated until convergence. Note that ${\mathcal{E}}_{\underset{¯}{\psi }}(\theta )$ is quadratic in f(ψ) = (I₀–I₁(ψ)), so it can be solved analytically to obtain an update for θ at every iteration at fixed $\underset{¯}{\psi }$–as the values θ_j are assumed to be independent. The updating of θ_j can be performed via the following closed form expression given by

13

${\theta }_j=\left(\frac{\eta }{2}+\sqrt{\frac{{\eta }^2}{4}+\frac{{\left(I_0–I_1\right)}_j^2}{2}}\right).$

The computational details for the first minimization problem posed by ${\mathcal{E}}_{\theta }\left(\underset{¯}{\psi }\right)$ at fixed θ are similar to the one described by Schmidt and Sutton [33], however this work uses symmetric warping–that is both I₀ and I₁ are warped. This has been demonstrated in testing to enhance the temporal accuracy as well as the dynamic range of the proposed method. The functional and its gradient are evaluated in the Bayes_min routine, and the pseudo-code is provided in Code 3.

Pseudo-code 1

In summary, the wavelet coefficients are used to compute $\underset{¯}{{\nabla }_{\psi }{J}_R}$ and J_R using the sub-routine Reg as described in Sec. 1.3.2.3. The wavelet coefficients $\underset{¯}{\psi }$ are inverted to obtain $\underset{¯}{u}$, which is then used to warp both of the particle images. The warped images are then used to compute J_D and $\underset{¯}{{\nabla }_{\psi }{J}_D}$. The integrand of J_D is used to determine the local weights $\theta (\underset{¯}{x})$ using Eq. 13.

Quality control

The code has been extensively tested on various different synthetic as well as experimental particle image data [31, 14, 35, 23, 33, 21, 16, 15, 20]. For unit tests, multiple synthetic datasets have been made available by the authors [17] that contain particle images along with the ground truth to evaluate the accuracy of the algorithm.

Operating system

Based on MATLAB 7.10.0 (R2010a): Windows, UNIX/Linux, Macintosh

Pseudo-code 2

Pseudo-code 3

Programming language

MATLAB 7.10.0 (R2010a), upward compatible

Additional system requirements

MATLAB 7.10.0 (R2019b): 1 GB disk space, 1 GB RAM.

Dependencies

Matlab Image Processing Toolbox.

Software location

Archive An archive version of the software is made available as an OSF repository.

• Name: OSF

• Persistent identifier: https://www.doi.org/10.17605/OSF.IO/9SK6X

• Licence: MIT license

• Publisher: Gauresh Raj Jassal

• Version published: 3.08

• Date published: 16th of January 2025

Code repository The code is hosted in a Github repository.

• Name: Github

• Persistent identifier: https://shrediquette.github.io/PIVlab/

• Licence: MIT license

• Date published: 16th of December 2024

Standalone Windows version A standalone GUI version is also available for Windows users.

• Name: Standalone Windows

• Persistent identifier: https://www.pivlab.de

Language

MATLAB 7.10.0