Variational second order flow estimation for PIV sequences

RESEARCH ARTICLE Variational second order flow estimation for PIV sequences L. Alvarez ÆC. A. Castan ˜oÆM. Garcı ´aÆK. Krissian Æ L. Mazorra ÆA. Salgado ÆJ. Sa ´nchez Received: 14 July 2006 / Revised: 25 July 2007 / Accepted: 21 September 2007 / Published online: 10 October 2007 Springer-Verlag 2007 Abstract We present in this paper a variational approach to accurately estimate simultaneously the velocity field and its derivatives directly from PIV image sequences. Our method differs from other techniques that have been presented in the literature in the fact that the energy minimization used to estimate the particles motion depends on a second order Taylor development of the flow. In this way, we are not only able to compute the motion vector field, but we also obtain an accurate estimation of their derivatives. Hence, we avoid the use of numerical schemes to compute the derivatives from the estimated flow that usually yield to numerical amplification of the inherent uncertainty on the estimated flow. The performance of our approach is illustrated with the estimation of the motion vector field and the vorticity on both synthetic and real PIV datasets. 1 Introduction Particle image velocimetry (PIV) is a relatively new technique widely used to study the temporal evolution of all kinds of flows. From the initial work by Adrian (1988), where the theoretical framework for PIV was set, this technique has been a continuous research topic (Thomas et al. 2005). The basic idea behind this technique is to make visible fluid motion by adding small tracer particles to the fluid and infer the flow velocity field from the images of the position of these particles at two instances of time (Westerweel 1997). The interested reader is referred to the work of Keane and Adrian (1990), where the experimental design rules are detailed. In the context of PIV research, the standard procedure to infer the flow velocity field from the image pairs is to compute the cross-correlation function between local interrogation windows of consecutive frames (Keane and Adrian 1992). In this way, the local maximum of that function determines the displacement of the pixel values. Hierarchical and iterative schemes have been developed to overcome the limitations of correlation techniques. In this sense, a lot of work has been done since the pioneering work of Soria et al. (1996) and it is still an interesting research topic Westerweel et al. (1997), Scarano (2002), Ruhnau et al. (2005b), and Nogueira et al. (2005). Recently, variational motion estimation approaches widely used in computer vision or image processing domains have been proposed for PIV (Corpetti et al. 2002; Ruhnau et al. 2005b; Corpetti et al. 2006). Encouraging L. Alvarez (&)C. A. Castan ˜oM. Garcı ´aK. Krissian  L. Mazorra A. Salgado J. Sa ´nchez Departamento de Informa ´tica y Sistemas, Universidad de Las Palmas de Gran Canaria, Campus de Tafira s/n, 35017 Las Palmas de Gran Canaria, Spain e-mail: [email protected] C. A. Castan ˜o e-mail: [email protected] M. Garcı ´a e-mail: [email protected] K. Krissian e-mail: [email protected] L. Mazorra e-mail: [email protected] A. Salgado e-mail: [email protected] J. Sa ´nchez e-mail: [email protected] 123 Exp Fluids (2008) 44:291–304 DOI 10.1007/s00348-007-0402-3 results have also been obtained with these approaches, although some limitations have also to be overcome. For instance, variational methods are usually based on the assumption that light conditions do not change from frame to frame. However, real fluid images usually suffer from temporal intensity distortions, yielding to wrong results. In this sense, more robust methods with respect to light conditions have been proposed by Corpetti et al. (2002) based on an integrated version of the continuity equation of fluid mechanics. In any case, either with correlation-based methods or variational approaches, there is still a certain level of noise inherent to the estimation method (Foucaut and Stanislas 2002). This noise may introduce a significant error on the flow components, but it is even more important for the computation of the derivatives of the motion field due to numerical error amplification. Since flow derivatives are a great source of information in fluid dynamics (they appear in the Navier–Stokes equations and important parameters such as divergence, vorticity, the strain rate tensor and dissipation rate), it is particularly important to obtain an accurate estimation of flow derivatives. Usually the computation of the flow derivatives is a two steps procedure: First we compute a flow estimation using any method we want and second we compute flow derivatives using some kind of finite difference scheme. In this paper we present a variational method to compute simultaneously the flow and its derivatives directly from the image sequence. The proposed method is based in an energy minimization and a second order Taylor development of the flow. We illustrate the behavior of our approach by the estimation of the out-of-plane vorticity, which is an interesting measure to detect and to study the evolution of vortices in the flow (Hunt 1987). Several approaches have been proposed in the literature to achieve this task. On the one hand, discrete differential operators based on different numerical schemes are directly applied to the estimated velocity vector field. In general, their accuracy depends on the initial uncertainty inherent to the vector field, the sampling rate, neighborhood correlation, etc. (Etebari and Vlachos 2005; Luff et al. 1999; Raffel et al. 1998; Kim and Lee. 1996). A second technique was proposed by Fouras and Soria (1998) which consists on analytical derivation of the estimated vector field, once it has been fitted to a base of second order polynomials by means of av 2 fitting procedure. A third technique is presented by Raffel et al. (1998) and Abrahamson and Lonnes (1995) based on the computation on the circulation around an enclosed area, which is related to the vorticity by the Stokes theorem. The latter two methods are detailed in Sect. 3. The paper is organized as follows. In Sect. 2,we describe the details of our variational approach and the way we have adapted the energy function to be minimized in order to directly estimate significative fluid parameters such as vorticity and strain rate tensor components. Next, in Sect. 3we present an overview of standard techniques from the literature to compute those quantities. In Sect. 4 we present numerical experiments on synthetic and real data and we analyze the accuracy of the proposed method in comparison with the standard techniques described in Sect. 3, as well as we characterize the sensitivity of our method with respect to the main parameters. Finally in Sect. 5we present the main conclusions of the paper. 2 A variational approach for second order motion estimation To estimate the optic flow of a given sequence we propose a variational approach based on the minimization of an energy function Eð~ uÞdefined for each point x 0 with the special fact that it depends not only on the displacement vector components u=(u,v) T , but also on their partial derivatives ~ u¼ðu;v;ux;uy;vx;vyÞT;as shown in the following equation: Eð~ uÞ¼Eðu;v;ux;uy;vx;vyÞ ¼Z Xðx0Þ Krðxx0ÞI1ðxÞI2ðxþðu;vÞT  þuxuy vxvy  ðxx0ÞÞ2 dxð1Þ where I 1 (x) is the first image of the sequence and I 2 (x) is the following frame, where the time increment from frame to frame is assumed to be normalized (Dt= 1), K r (x–x 0 )isa Gaussian kernel with standard deviation rwhich weighs the pixels in the domain X(x 0 ) centered at point x 0 . Hence, the goal is to find the components of vector ~ u;that is, the displacement and its partial derivatives such that minimize the error between I 1 (x) and I 2 (x) displaced by the unknown flow. At a first glance, the dependence on the partial derivatives (u x ,u y ,v x ,v y ) might seem a redundant operation since they can be computed from the obtained motion vector field. However, it can be shown that numerical differentiation of the estimated motion vector field yields to inaccurate results mainly due to undesired numerical error amplification (Raffel et al. 1998). The partial derivatives of the motion vector field are the components of the velocity gradient tensor, also called the rate of deformation tensor, D=du/dxgiven by: D¼du dx¼uxuy vxvy  ð2Þ This tensor provides quantitative information of the unitary deformation rate of the flow. In a general context, 292 Exp Fluids (2008) 44:291–304 123 this tensor is not a symmetric matrix and it has four independent components. In the context of fluid flow analysis, this deformation rate tensor is usually decomposed into meaningful components related to physical phenomena. A common decomposition method is to divide the tensor into a symmetric and an antisymmetric matrix, as shown in Eq. 3. du dx¼ux1 2uyþvx  1 2vxþuy vy "# þ01 2uyvx  1 2vxuy 0 "# ð3Þ This decomposition is uniquely determined for any matrix A, where the symmetric part is given by AþAT 2and represents the strain rate tensor with the elongational strains on the diagonal and the shearing strains on the offdiagonal, whereas the antisymmetric part is given by AAT 2; which is the rate of rotation tensor, whose non-zero elements are the vorticity components. Since vorticity and the strain rate tensor have a more relevant physical meaning in the context of fluid analysis than just the flow gradient it is possible to reformulate Eq. 1as a function of the vorticity xand the strain rate tensor components e 1,1 ,e 1,2 ,e 2,2 performing the following change of variables: ux¼e1;1uy¼xþe1;2 vx¼e1;2xvy¼e2;2 ð4Þ and now the expression of the energy we want to minimize is: Eð~ uÞ¼Eðu;v;x;e1;1;e1;2;e2;2Þ¼ Z Xðx0Þ Krðxx0Þ I1ðxÞI2ðxþðu;vÞTþe1;1e1;2þx e1;2xe 2;2  ðxx0ÞÞ  2 dx ð5Þ 2.1 Energy minimization In order to be able to obtain the six components vector of unknowns ~ u¼ðu;v;x;e1;1;e1;2;e2;2ÞTthat minimize Eq. 5 we formulate the solution at step n+ 1 as a function of the solution at step nand the six components vector of residuals ~ h¼ðhu;hv;hx;he1;1;he1;2;he2;2ÞTcomputed at each step, as shows the following expression: ~ unþ1¼~ unþ~ hð6Þ Then, introducing Eq. 6in Eq 5we are able to obtain an iterative operator towards the local minimum of the energy function similar to a gradient descent algorithm Fletcher and Reeves (1964). But first, the following approximations are still necessary to simplify the term that depends on I 2 (x) in order to obtain our iterative operator: I2 xþðunþ1;vnþ1ÞT þenþ1 1;1enþ1 1;2þxnþ1 enþ1 1;2xnþ1enþ1 2;2 ! ðxx0Þ! ¼I2xþðun;vnÞTþen 1;1en 1;2þxn en 1;2xnen 2;2  ðxx0Þ  þðhu;hvÞTþhe1;1he1;2þhx he1;2hxhe2;2  ðxx0Þð7Þ To simplify notation, let us define I 2 n as: In 2ðxÞ¼I2 xþðun;vnÞT þen 1;1en 1;2þxn en 1;2xnen 2;2 ! xx0 yy0  !ð8Þ where the vector x–x 0 has been decomposed into their components (x–x 0 ,y–y 0 ) T . The partial derivatives are denoted as: In 2;xðxÞ¼oI2 ox xþðun;vnÞT þen 1;1en 1;2þxn en 1;2xnen 2;2 ! xx0 yy0  !ð9Þ In 2;yðxÞ¼oI2 oy xþðun;vnÞT þen 1;1en 1;2þxn en 1;2xnen 2;2 ! xx0 yy0  !ð10Þ Then, using the new notation, the solution at step n+1 can be approximated by its Taylor expansion over ~ unto provide: Inþ1 2ðxÞ’In 2ðxÞþ ðhu;hvÞT þhe1;1he1;2þhx he1;2hxhe2;2 ! xx0 yy0  !TIn 2;xðxÞ In 2;yðxÞ ! ð11Þ Finally, by using vectorial notation I 2 (x) n+1 can be expressed as: Exp Fluids (2008) 44:291–304 293 123 Inþ1 2ðxÞ’In 2ðxÞþ^ In 2ðxÞT~ hð12Þ where ^ In 2ðxÞ¼ðIn 2;xðxÞ;In 2;yðxÞ;ðyy0ÞIn 2;xðxÞðxx0ÞIn 2;yðxÞ; ðxx0ÞIn 2;xðxÞ;ðxx0ÞIn 2;yðxÞþðyy0ÞIn 2;xðxÞ; ðyy0ÞIn 2;yðxÞÞTð13Þ Hence, introducing Eq. 12 in Eq. 5we obtain the approximation for the energy function shown in Eq. 14, where, in addition, we introduce a second term in the energy function weighted by the parameter a. The role of this regularization term is to provide at every point a smooth vector field, by means of an additional constraint on the norm of the vector ~ h;which is forced to be small. In this sense, the parameter adetermines the importance of this additional constraint on the vector field. ~ Eð~ hÞ¼ Z Xðx0Þ Krðxx0ÞI1ðxÞIn 2ðxÞ^ In 2ðxÞT~ h  2dx þaZ Xðx0Þ k~ hk2dxð14Þ This formulation allow us to easily obtain an analytical expression to compute the local minimum of the energy as a function of ~ h: r~ Eð~ hÞ¼2Z Xðx0Þ Krðxx0ÞI1ðxÞIn 2ðxÞ^ In 2ðxÞT~ h  ^ In 2ðxÞdxþ2aZ Xðx0Þ ~ hdx¼0ð15Þ which is equivalent to: Z Xðx0Þ Krðxx0ÞI1ðxÞIn 2ðxÞ  ^ In 2ðxÞdx ¼Z Xðx0Þ Krðxx0Þð^ In 2ðxÞ^ In 2ðxÞTÞþaI  ~ hdxð16Þ This system of equations can be expressed with the standard matrix notation Ax =b, where the vector of unknowns in this case is the vector ~ h;while the system matrix and the independent term can be computed as: A¼Z Xðx0Þ Krðxx0Þð^ In 2ðxÞ^ In 2ðxÞTÞþaI  dxð17Þ b¼Z Xðx0Þ Krðxx0ÞI1ðxÞIn 2ðxÞ  ^ In 2ðxÞdxð18Þ The solution of the system of equations is given by ~ h¼A1b;so we only have to invert the 6 ·6 matrix Aor use any other algorithm to solve the system of equations. We observe that if a[0, matrix Ais positive definite and therefore the system of equations is well posed (it has a unique solution). It means that the parameter aavoids instabilities in the solution of the linear system of equations. 2.2 Numerical implementation Next, we describe the main steps we followed to derive an efficient algorithm of the method proposed in Sect. 2, including the numerical considerations involved on the computation of integrals and spatial derivatives in the proposed method. As it can be seen in Algorithm 1, our variational approach starts from an initial estimation of the motion vector field u 0 which can be obtained using any other optic flow estimation method. Then, we start the iterative procedure towards the local minimum of the energy function given in Eq. 5. At each iteration we check the convergence of our algorithm in order to discard the solutions that do not provide the optimal response. This algorithm is then applied for every point in the image (or a grid at a given scale) to obtain the desired first and second order flow parameters estimation, using information within a neighborhood of that point determined by the rparameter of the Gaussian kernel. The influence of the main parameters in this algorithm, rand a, is detailed in Sect. 4. Concerning other parameters like N iter , we remark that it should be as big as possible since it is automatically truncated when convergence is detected. Finally, the initialization of the energy at the first step requires the computation of the image derivatives. In our implementation, we have used central finite differences because it is a good compromise between an easy implementation, low time of computation and low error propagation Raffel et al. (1998). 3 Standard methods for vorticity estimation In order to obtain an estimation of the vorticity field, alternatives to finite differencing have also been proposed to compute such magnitude. In this section we describe two of these methods that we will use later to compare the 294 Exp Fluids (2008) 44:291–304 123 performance of standard vorticity estimation algorithms with the performance of our variational approach. 3.1 Circulation method By definition the vorticity is related to the circulation by Stokes theorem: C¼Iudl¼ZðruÞdS¼ZxdSð19Þ where ldescribes the path of integration around a surface S. The vorticity for a fluid element is found by reducing the surface Sand with it the path l, to zero: ^ nx¼^ nru¼lim S!0 1 SIudlð20Þ where the unit vector ^ nis normal to the surface S. Stokes theorem can also be applied to the two dimensional vector field: ðxZÞi;j¼1 ACi;j¼1 AIlðx;yÞðu;vÞdlð21Þ where ðxZÞi;jreflects the average vorticity within the enclosed area. In practice, Eq. 21 is implemented by choosing a small rectangular contour around which the circulation is calculated using a standard integration scheme such as the trapezoidal rule. The local circulation is then divided by the enclosed area to arrive at an average vorticity, as shows the following equation which uses a neighborhood of eight points: Exp Fluids (2008) 44:291–304 295 123 ðxzÞi;jffiCi;j 4DxDyð22Þ where Ci;j¼1 2Dxðui1;j1þ2ui;j1þuiþ1;j1Þ þ1 2Dyðviþ1;j1þ2viþ1;jþviþ1;jþ1Þ 1 2Dxðuiþ1;jþ1þ2ui;jþ1þui1;jþ1Þ 1 2Dyðvi1;jþ1þ2vi1;jþvi1;j1Þ ð23Þ 3.2 Second order polynomial fit and analytic differentiation The second standard method we use to compare our method with was initially proposed in Fouras and Soria (1998). For this method, it is necessary to provide an initial estimation of the motion vector components (u i,j ,v i,j ) within a certain neighborhood centered at the point (i 0 ,j 0 ) where we want to compute the vorticity. Then, the components (u i,j ,v i,j ) are separately fitted to a basis of polynomials {P k (x,y)} of K-th power using a v 2 procedure, as proposed in Press et al. (1992). The linear combination of the polynomials requires M=(K+1) 2 coefficients. Thus, the number of data points Nin the neighborhood required for the v 2 fitting process is N‡M. In other words, the v 2 fitting procedure provides the coefficients u k and v k we need to express our initial motion estimation vector as a linear combination of the elements of the basis of polynomials: ui;j¼X M1 k¼0 ukPkðx;yÞð24Þ vi;j¼X M1 k¼0 vkPkðx;yÞð25Þ Then, it is straightforward to analytically differentiate the previous expressions. The numerical value of the vorticity at the current position (i,j) is given by the following expression evaluated in (x=0,y= 0): ðxzÞi;j¼X M1 k¼0 vkoPkðx¼0;y¼0Þ oxukoPkðx¼0;y¼0Þ oy  ð26Þ In this method, the set of basis functions and the samples around the point of interest can be seen as parameters of the method since different possibilities are possible. In our study, as proposed in Fouras and Soria (1998), we have chosen the following set of second order polynomials in x and yas basis functions: Pðx;yÞ2f1;x;x2;y;y2;xy;x2y;xy2;x2y2gð27Þ This results in M= 9 coefficients to be fitted in the decomposition described by Eqs. 24 and 25. Thus, we need at least 9 samples from a 3 ·3 neighborhood centered at point (i,j), which is the neighborhood used in our experiments. 4 Numerical experiments In order to show the performance of the proposed approach, we present here the numerical experiments we have performed to evaluate the accuracy of the method. We perform experiments on synthetic and real data. In the case of synthetic experiments we use flow models with a known analytic expression, in this way, we are able to compare our results with the true solution we are looking for, also called ground truth in the field of computer vision. In this sense we focus our attention in the flow estimation improvement we obtain as well as in the vorticity estimation because its interest in flow analysis. In addition, we compare our vorticity estimation with the methods detailed in Sect. 3. In any case, the algorithm proposed in this paper starts from an initial estimation of the flow which can be obtained by any desired method. This initialization must be, at least, a rough approximation of the true motion vector field in order to the algorithm be able to converge. In our work, we use two different methods drawn from the state of the art to achieve such initialization in order to evaluate its influence on the obtained results. The first method (PDE) used in our experiments is based on a classical PDE scheme, presented in Alema ´n et al. (2005), composed by a data term that assumes that the image intensity is equal in two corresponding points and a regularization term introduced in Anandan (1989) and studied in Alvarez et al. (2000) which allows discontinuities preservation on the flow. In our implementation, we use a pyramidal approach in order to speed up the algorithm and to avoid the convergence towards spurious local minima. The second method (COR) used in our experiments is based on the multi-step iterative computation of crosscorrelation measures where the position of the interrogation window is displaced with subpixel precision at each iteration depending on the precomputed flow, as it is described in Scarano (2002). To find the peak location in the correlation map, we use a Gaussian interpolation scheme based on a 3 ·3 kernel Westerweel (1993), 296 Exp Fluids (2008) 44:291–304 123 whereas a bilinear interpolation scheme to perform the displacement of the interrogation window with subpixel precisition. 4.1 Experiment 1: Synthetic vortex flow In our first experiment, we have built a sequence of two PIV images of 1,024 ·1,024 pixels. The first image was synthesized by the research institute CEMAGREF (Rennes, France) using a uniform random distribution of the particles over the whole image, with a constant particle concentration of 256 particles each 32 ·32 window. The particle size follows a random normal distribution with an average of 1.25 pixels and a standard deviation of 0.25 pixels. The second image is computed by applying the synthetic vortex flow modeled in Eq. 28 and shown on the left side of Fig. 1. In this sense, we are able to obtain the analytic expression of the vorticity, given in Eq. 29 and shown on the right side of Fig. 1, which can be used as the ground truth to evaluate our approach. uðxÞ¼ u¼3y ffiffiffiffiffiffiffiffiffi x2þy2 p v¼3x ffiffiffiffiffiffiffiffiffi x2þy2 p 8 < :ð28Þ x¼vxuy¼3 ffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2þy2 pð29Þ Table 1shows some quantitative results that demonstrate the performance of our method in comparison with other methods drawn from the state-of-the-art, which were briefly described at the beginning of this section. In our experiments, we computed the average angular error and average Euclidean error of the flow components directly obtained from the PDE and COR methods. Then, we use those responses as an initialization of our method, yielding to a significant improvement of the error measures, as it can be seen in the table. Finally, to evaluate the noise introduced in our approach, we initialize our numerical scheme with the ground truth flow field, yielding to an error rate of the same order of magnitude as in the previous case. Figure 2shows a detail of the results obtained with the different algorithms under study using a PDE scheme to compute the initial estimation of the flow. On the left side, we show the vorticity obtained with our variational method. In the middle, the vorticity obtained with the circulation method detailed in Sect. 3.1 is shown, while the image on the right corresponds to the v 2 fitting method with analytic differentiation (v 2 method, from now on) detailed in Sect. 3.2. The results show that, in all cases, a data validation step is required in order to eliminate the outliers that appear on the images, which mainly appear by the fact that usually not all the particles are detected correctly in a PIV sequence, as it is justified in Ruhnau et al. (2005a). Hence, this effect leads to a bad convergence of the algorithm on spread points over the whole image. It is important to remember that the proposed energy minimization is resolved locally for every point of the image and, since no global regularization term is used to impose global regularity on the estimated vector field. Hence, the convergence from the input data may lead to a local minimum different from the global one. Fortunately, the noise pattern that appears can be easily removed using a median filter, as it is described in Senel et al. (2002), as a data Fig. 1 Left Synthetic vortex flow generated with Eq. 28. Right True Vorticity field obtained from the analytical expression (Eq. 29) Table 1 Angular and Euclidean error computed between the vortex flow in Eq. 28 and the estimations performed with the different approaches studied in this paper: PDE, COR and our method, where the initialization flow is specified in brackets Error Mean angular error Mean Euclidean error True, PDE 0.4308 3.1e-2 True, COR 0.2737 1.9e-2 True, Our method (PDE) 0.0328 9.3e-4 True, Our method (COR) 0.0347 9.7e-4 True, Our method (True) 0.0155 2.9e-4 Exp Fluids (2008) 44:291–304 297 123 validation post-processing step. In our experiments, to avoid the modification of correct signal values, we only substitute the current value by the median when we detect that it might be an outlier, since its value is quite different from the median. Figure 3shows the validated results after the data validation step. Now, it can be seen that the dispersed noise has successfully been removed improving the quality of the estimations. In order to check the influence of the initial flow estimation on the results, we use again the COR scheme as initialization to compute the vorticity. Figure 4shows a detail of the results obtained using the new initial estimation as input for the vorticity computation algorithms once the outliers have been removed with the median filter, as it was done in the last case. Fig. 2 Estimated vorticity before removing the outliers with a vortex flow model using PDE scheme for initial estimation. Left Results with our variational approach. Middle Results with the circulation method. Right Results with the v 2 method Fig. 3 Estimated vorticity after removing the outliers with a vortex flow model using PDE scheme for initial estimation. Left Results with our variational approach. Middle Results with the circulation method. Right Results with the v 2 method Fig. 4 Estimated vorticity after removing the outliers with a vortex flow model using a correlation scheme for initial estimation. Left Results with our variational approach. Middle Results with the circulation method. Right Results with the v 2 method Table 2 Mean error computed for the different vorticity estimation methods Init flow Our method Circulation v 2 method True 3.07e-5 3.45e-6 4.10e-6 PDE 1.35e-4 2.37e-3 2.39e-3 COR 1.38e-4 1.80e-3 1.80e-3 Our method (PDE) N/A 2.50e-4 1.91e-4 Our method (COR) N/A 2.75e-4 2.01e-4 298 Exp Fluids (2008) 44:291–304 123 Finally, in Table 2, we present quantitative error measurements, following Eq. 30, to evaluate the accuracy of our approach for vorticity estimation in comparison with the ground truth model. The first column indicates the estimated vector field from which the vorticity is computed with the corresponding schemes described in Sect. 3. In our experiments, we first compute the vorticity field from the ground truth vector field. In this case, the error rate obtained with the circulation and v 2 methods outperforms our results, since those schemes perform the computation directly from the true field and our method slightly modifies it. On the other situations, our method outperforms their results because of the numerical amplification of the initial error on the estimated flow. In this sense, even the vorticity given by the circulation or v 2 methods, using the flow provided by our method, yields to a lower error rate than using the PDE or COR flows. Error ¼PNx i¼1PNy j¼1jxTRUE i;jxest i;jj NxNyð30Þ 4.2 Experiment 2: Synthetic Lamb–Oseen flow To validate our algorithm, we have again built a sequence of two PIV images of 1,024 ·1,024 pixels using the same procedure detailed in the previous experiment, but following the mathematical flow model shown in Eq. 31, also known as Lamb–Oseen vortex: uðxÞ¼ u¼C0y 2pffiffiffiffiffiffiffiffiffi x2þy2 p1expðx2þy2 4mtÞ  v¼C0x 2pffiffiffiffiffiffiffiffiffi x2þy2 p1expðx2þy2 4mtÞ  8 > < > :ð31Þ where mis the cinematic viscosity and C 0 is the initial circulation. In our experiments we used C 0 = 0.05 m 2 /s and ffiffiffiffiffiffi 4mt p¼1 6;which can be seen as the radius of the Lamb–Oseen Vortex. Then, Eq. 32 presents the analytic expression of the vorticity. x¼vxuy¼C0 p4mtexp x2þy2 4mt  ð32Þ Figure 5shows in the left side the synthetic Lamb– Oseen vortex and in the right side the vorticity computed with the parameters used in the computation of the flow. In Table 3we show some quantitative results that demonstrate the behavior of our method. As it was done in the previous experiment, we computed the average angular error and average Euclidean error of the flow components directly obtained from the PDE and COR methods. Then, these flows are used as an initialization of our approach and again the flow obtained is improved. To emphasize the behavior of our approach with respect to the noise, we computed the error rate using the ground truth flow as initialization. In this case, we obtain a very low error rate, which shows that the noise introduced by our method is almost negligible. Figure 6shows the estimated out-of-plane vorticity obtained with our approach (image on the left) in comparison to the circulation method (in the middle) and the v 2 fitting method (on the right) once a median filter has been used to eliminate the outliers as a data validation postprocessing step. Fig. 5 Left Synthetic Lamb– Oseen vortex generated with Eq. 31.Right True vorticity field obtained from the analytical expression (Eq. 32) Table 3 Angular and Euclidean error computed between the Lamb– Oseen flow in Eq. 31 and the estimations performed with the different approaches studied in this paper: PDE, COR and our method, where the initialization flow is specified in brackets Error Mean angular error Mean Euclidean error True, PDE 0.1132 0.0574 True, COR 0.9642 0.4601 True, Our method (PDE) 0.0571 0.0083 True, Our method (COR) 0.8938 0.1968 True, Our method (True) 0.0112 0.0015 Exp Fluids (2008) 44:291–304 299 123