鲁棒光学中动态自洽的相关变分方法用于其流场分析检测方案(粒子图像测速)

Exp Fluids (2008) 45：595-608DOI 10.1007/s00348-008-0567-4 Exp Fluids (2008) 45：595-608596 Dynamic consistent correlation-variational approach for robustoptical flow estimation D. Heitz·P. Heas · E. Memin·J.Carlier Abstract We present in this paper a novel combinedscheme dedicated to the measurement of velocity in fluidexperimental flows through image sequences. The pro-posed technique satisfies the Navier-Stokes equations andcombines the robustness of correlation techniques with thehigh density of global variational methods. It can be con-sidered either as a reenforcement of fluid dedicated optical-flow methods towards robustness， or as an enhancement ofcorrelation approaches towards dense information. Thisresults in a physics-based technique that is robust undernoise and outliers， while providing a dense motion fieldThe method was applied on synthetic images and on realexperiments in turbulent flows carried out to allow athorough comparison with a state of the art variational andcorrelation methods 1 Introduction There has been significant progress during recent years inthe analysis of particle image sequence to estimate fluidvelocities. On one hand regions-based techniques produce D. Heitz· P. Heas·J. Carlier Cemagref， UR TERE， 17 Avenue de Cucille， CS64427，.35044 Rennes， France D. Heitz () P. Heas· J. CarlierUniversite Europeenne de Bretagne， Rennes， Francee-mail： dominique.heitz@cemagref.fr Present Address： P. Heas· E. Memin INRIA， Campus Universitaire de Beaulieu， 35042 Rennes， France velocities from maximum correlations between displacedinterrogation windows. These robust techniques providereliable but low-pass filtered vector fields， compared to thedensity of particles. On the other hand model based mea-surement techniques， known as optical flow approaches，rely on a model， the so-called data term or observationterm， which relates the observed image intensity to thedisplacements fields to estimate. Since the single use ofthe data term leads to an improperly posed problem—onlythe displacements along the direction of the intensity gra-dient can be measured-additional assumptions have beenproposed. The local approach of Lucas and Kanade (1981)tackles this problem while considering a local spatialconstancy assumptions on the optic flow field. Biginand Granlund (1988) extended this technique to localspatiotemporal constancy assumptions. Like correlationapproaches， these local schemes are region-based methods，providing robust but sparse estimations. Alternatively， a relatively simple important idea in theconceptual development of model-based measurementtechnique is to introduce a regularization of the velocityfield to devise a global approach (Tikhonov and Arsenin1977). Based on the minimization of an energy function，composed of a data term and a regularization term， thisglobal scheme provides dense vector fields with spatialcoherence. It should be noticed， that contrary to approa-ches using post-treatmentss to filter spurious； vectors.global approaches seek regular and coherent displace-ments fields. This promotion of regular and coherentsolutions is conducted from the beginning of the mini-mization. Hence， the key of global methods's accuracy isdirectly link with the regularization term. However， theseapproaches provides a framework to introduce physics-based constraints， thus giving opportunities to introducedprior knowledges of the dynamics of fluids to extract the motion through image sequences. IIn thispaperWeexamine'ss new physics-basedconstraintsleading todynamic consistency. The global scheme proposed by Horn and Schunck(1981)， using the optical flow constraint equation (OFC)，and a first-order regularization， was applied on particleimage by Ruhnau et al. (2005). However， this standardapproach was at the origin adapted to two-dimensional(2D) rigid motions. With particle image of turbulent fluidflows， the OFC equation-based on the variational for-mulation of the brightnessconstancy asasssumption-isviolated due to out of plane motions， and the regularizationpromotessmoothhvector fields.Thefluid dedicatedapproach of Corpetti et al. (2002) uses a data term based onthe continuity equation and a second-order regularizerconstraining the gradients of the flow components diver-gence and vorticity. This new method is able to extract theapparent 2D motion of 3D flow while preserving the spatialvariations of the apparent divergence and vorticity. Appliedto particle images it provides reliable results (Corpetti et al.2006). The optical flow method used by Corpetti et al.based on the continuity equation is physically consistentonly for transmittance images. Recently， Liu and Shen(2008) gave the mathematical definition and a physicalmeaning of this optical flow formulation for particle ima-ges. They showed that the continuity equation fortunatelyleads to the physics-based data term they derived， with theassumption that the boundary flux terms can be neglectedfor particle images. The accuracy (of second-order regularizationwasimproved with carefully designed Helmholtz decomposi-tion (Kohlberger et al. (2003) and discrete orthogonaldecomposition (Yuan et al. 2007). For 2D flows this lastmethod provides better results than advanced competitivecorrelation approaches (Heitz et al.2008). To obtain anaccurate estimation of the derivatives of the velocity，Alvarez et al. (2008) proposed an optical flow approachbased on a second-order Taylor development of the flow. Irthis technique the physical meaning was enhanced withdecomposition of the velocity gradient tensor into sym-metric and antisymmetric parts， thus exhibiting tlthevorticity and strain rate tensor. Recently， Rhunau andSchnorr (2007) devised a physically grounded regulariza-tion based on the Stokes equation， therefore restricted toviscous flow (i.e. for very small Reynolds numbers). Another way to introduce physics in the estimation ofmotion is to satisfy the Navier-Stokes equations. Thisapproach was first introduced by Okuno et al. (2000) topost-treat standard correlation measurements. Like themethod of Rhunau and Schnorr (2007)， this physic-basedmeasurement technique provides estimations of the veloc-ity and of the pressure. Using the incompressible velocity-vorticity formulation of the Navier-Stokes equations， Rhunau et al. (2007) developed a spatio-temporal regu-larizer for a global approach equipped with an additionalconstraint that mimics the small viscous term. In the con-text of data assimilation into Navier-Stokes equations，Cuzol et al. (2007) presented a Bayesian filtering techniqueto track vorticity fields， and Papadakis et al. (2007) pro-posed a promising variational data assimilation techniqueto estimate an optimal dynamically consistent vector fields.All the above-mentioned spatio-temporal techniques areonly valid for 2D flows， which represents a strong limita-tion for the analysis of most of real flows. Furthermore， variational approaches can exhibit otherlimitations with particle image sequences. The variationalformulation is only suited for displacements smaller thanthe shortest wavelength present in the image， i.e. smallerthen the size of the observed particles. To handle largedisplacements the data term is generally embedded within amultiresolution coarse-to-fine scheme (Bergen et al. 1992).Standard multiresolution approach consists in deriving theoriginal frame into a pyramid of images， by successivelow-pass filtering (Gaussian smoothing) and regular sub-sampling by a factor of two in each direction. The createdpyramid structure then allows multiscale motion estimationby incremental procedure. Low resolution components areestimated at coarsest level and then refined step by step： ata resolution level， an increment velocity field is estimatedaround the projection (by duplication or interpolation) ofthe final estimate at previous resolution level. In generalthe procedure of filtering and subsampling leads to a loss ofinformation. In the case of particle images of turbulentfluid flows， small particles with large velocities， can besmoothed out by the multiresolution procedure. Therefore，due to this loss of information， the variational approachwill not extract reliable velocities. In presence of noise inthe image the multiresolution procedure enhances thisphenomenon. At the coarsest level， the subsampling pro-cess keeps statistically more noise than particle informationwhen there is poor particle density and noise in the originalimage. Bruhn et al. (2005) proposed a combined local-global method which has the robustness of local methodswith the density of global approaches. Unfortunately， thistechnique， not specifically suited for fluid flow， was nottested on particle images. In the present paper we describe a novel dynamic con-sistent global approach combining the robustness ofcorrelation techniques with the high density of globalvariational rmethods.. constrainedwithNavier-Stokesequations extending the estimations towards 3D flowsdriven by shallow water equations. The proposed approachis an adaptation of Heas et al. (2007) scheme-developedfor the measurement of atmospheric dynamics from satel-lite image sequences-for laser sheet illuminated particleimage. This results in a physics-based technique that is robust under noise and outliers， while providing a densespatio-temporal coherent motion field. This article is organized as follows. In Sect. 1， wepresent multiscale motion estimation focusing on the fluiddedicated variational methods. Then， in Sect. 2， we presentthe collaborative correlation-variationalaapproach， fol-lowed by the description of spatio-temporal constraint. InSect. 6， results from synthetic images based on a DNS of2D turbulence， and real images in the wake of a circularcylinder， are presented and analyzed. We provide someelements of comparison of our method with standardoptical flow techniques. 2 Multiscale motion estimation 2.1 Regions-based or local approach A classical approach used with PIV is correlation-basedmatching technique. The velocity u at point s is obtainedfor continuous image fields by minimizing where W(s) denotes windows centered on different pointsof the images and w is the displacement vector to estimate.For local optical flow approaches， which belong to regions-based methods， the negative product of the images isreplaced by the absolute value or the square of the DFDequation (see Sect. 2.2). Theseapproaches suffer from several deficiencies：traceable features must be sufficiently contrasted and mustpersist over time on consecutive images. Furthermore， theestimation prone to erroneous spatial variability， which canbe reduced with the use of filters. On the other hand， thesetechniques constitute very fast methods and are generallylocally very robust to noise. These techniques are based ondisjoint local estimation， and thus produce sparse vectorfields estimated locally and independently. 2.2 Global or dense approach 2.2.1 Classical optical flow method The apparent motion， perceived through variations ofimage intensity I is called optical flow in the computervision community. All optical flow estimation methods relyon the temporal conservation of some invariants. The mostcommon are photometric invariant which can easily beextracted and may lead to dense measurements. In the caseof PIV of bi-dimensional flows， the gray level conservationcan be assumed， leading to an integrated non-linear formulation called the displaced frame difference (DFD)equation I(s + u(S)， t+At)=I(s， t)， or at time t， to alinear differential formulation called the OFC equationu(s)· VI(s) +I(s)=0， where s denotes the spatial coor-dinates (x， y) and u(s)=(u，v) the apparent motion field atthis point and at a given time t. These two formulations cannot be used alone， as theyprovide only one equation for two unknowns at each spa-tio-temporal location (S，t). This is the well-known apertureproblem where the normal flow to the brightness gradient isestimated while the tangential velocity component remainsundetermined. In order to remove this ambiguity one mustrely on other assumptions. The mainly used assumption islocal spatial or spatio-temporal constancy of the velocityfields (Lucas and Kanade 1981； Bigiin and Granlund 1988). Alternate approaches， derived in the framework pro-posed by Horn and Schunck (1981)， have the greatadvantage of producing dense vector fields with spatialglobal coherence. They are based on the minimization ofan energy function J= Ja+J. composed of two terms.The first term Ja is called the data term and implements theOFC equation The penalty function d is usually the L2 norm but it maybe changed to a robust function attenuating the effectof data that deviate significantly from the OFC (Black andAnandanni1996)..Theieseconddterm jJr isiscalled theregularization term. It is usually a standard first-orderspatial smoothness term where a >0 is a parameter controlling the balance betweenthe smoothness and the global adequacy to the brightnessconstancy assumption..stands for the Euclidean norm.Function d may be the quadratic penalty if the searchedsolution is smooth everywhere or a robust norm function ifone wants to handle implicitly the spatial discontinuities ofthe field (Black and Anandan 1996) 2.2.2 Fluid dedicated optical flow method With most of real turbulent flows， the OFC equation dosenot hold. To estimate the apparent motion (i.e. the velocityin the image plane) of 3D flows， Corpetti et al. (2002，2006)proposed a data term based on the continuity equation. Thisnew constraint although derived for transmittance images isfortunately suited for laser-sheet illuminated particle ima-ges (Liu and Shen 2008). Liu and Shen demonstrated thatthe optical flow is proportional to the path averaged velocity of particles across the laser sheet and proposed the fol-lowing physics-based optical flow equation， where D is a diffusion coefficient， c iss a coefficientof particle scattering/absorption， B=-n·VN|r；-v.(Nr， VT2 +N|rVIi) is a boundary term that is relatedto N， the particle number per unit total volume， and itsderivatives coupled with the derivatives of the controlsurfaces T1， I2 of the laser sheet illuminated volume. Sincethe control surfaces are planar， there is no particle diffusionby molecular process， and the rate of accumulation of theparticle in laser sheet illuminated volume is neglected， theterm f(s， I)= 0 and the physics-based OFC can be written Furthermore in the case of fluid flows， it canbedemonstrated1thatafirst-order regularizationis notadapted as it favors the estimation of velocity fields with1OWvorticity (Corpetti et al. 2006). Asecond-orderregularizationnccanadvantageouslybee considered asproposed by Suter (1994)： where =curlu=u-8 and =divu=+ standfor the vorticity and divergence of the apparent velocityu=(u， v). Note that in the case of bi-dimensionalincompressible flows， divergence of the apparent motion isequal to zero. High-order regularization leads to non-trivialminimization problems. As a consequence more sophisti-cated discretization schemes have to be used. Using finitemimetic differences， Yuan et al. (2007) proposed a high-order motion estimator providing accurate measurements(Heitz et al. 2008). In the present study， to circumvent the difficulty ofimplementing second-order smoothness constraint， thesecond-order regularization term was simplified in a com-putational point of view in two interleaved first-order div-curl regularizations based on auxiliary variable i and approximating the vorticity and the divergence of the flow(Corpetti et al. 2002，2006). Thus we have： where A is a positive regularization parameter. AlthoughYuan et al. (2007) numerical scheme provides betterresults， here the Corpetti et al. regularization is used forsimplicity. 2.3 Multiresolution approach A problem with variational approaches is large displace-ment estimation. The intensity function must be locallysufficiently close to a linear function. Since the larger thedisplacement the more narrow the linearity domain is， largedisplacements cannot be recovered directly. A commonway to overcome this limitation is to create a image pyr-amid， constructed by successive low-pass filtering anddown sampling of the original images. In this framework，principal components of displacements are first estimatedat coarse resolution where motion amplitude should besufficiently reduced in order to satisfy the linearityrequirement. Incremental displacement are then estimatedwhile going down the pyramid， and the solution is refined(Bergen et al. 1992) (Fig.1). Themultiresolution schemes estimates principalcomponent displacements only at coarse resolutions wheresmall particles are rubbed out. Hence， this approachenables the characterization of large displacements of smallparticles only in the case when their motion are closeenough to the principal component’s one. Since at thecoarsest level， the subsampling process keeps statisticallymore noise than particle information， especially when thereis poor particle density or small particles (around 1-2 pixeldiameter) and noise in the image， the multiresolutionapproach leads to erroneous estimations. 3 Collaborative correlation-optical flow scheme To overcome the multiresolution limitations， we proposean alternative approach relying on a unique representationof the full resolution image. Thus， the proposed method Fig. 1 Coarse to fine resolution with multiresolution representationof the images tackles the non-linear estimation problem without makingsuccessive approximations in the calculation of interpo-lated images nor restricting itself to the characterization oflarge displacements of sufficiently large structures. Theproposed method takes advantage of the non-linear for-mulation of the motion estimation problem within adifferential framework appropriated for globalized localsmoothing. More explicitly， large displacements of both，fine or large structures， can be recovered by correlation-based methods for sufficiently contrasted and persistingimage regions. Thus， the idea of the method is to replacethe coarse estimates of the multiresolution scheme by adense large scale displacement estimate derived from acollection of correlation-based vectors ue obtained byEq. 1. A new functional was defined to estimate the vectorfields where J() is the energy function constraining estimateddisplacements u to be close to a sparse correlation-basedvector field uc. Thus， functional J() is defined as aquadratic distance between the solution and a collection ofcorrelation-based vectors u=(u’，v) located at the points=(x’，y) and influencing their neighborhood accordingto a bi-dimensional Gaussian N'(s) of mean s'：and ofvariance c where g and y denote respectively the confidence and thefunctional weighting factors. 4 Spatio-temporal smoothing 4.1 Energy function We propose to incorporate a priori physical knowledge onfluid dynamical evolution in the combined correlation-optical flow estimation scheme. In order to preservespatio-temporal consistency of displacement estimates， asimplified Navier-Stokes dynamical model is adapted toimages depicting 3D fluid flows. A dense displacementfield is predicted by time integration of a physicaldynamical model. The propagated field is then introducedin the estimation process as a spatio-temporal smoother. Apriori information introduced with this dynamic model mayenhance significantly the quality of estimates especially inthe case of noisy observations. Keeping notations of the previous section， we define forthe estimation of variable u the new functional where is an energy functionnCconstrainingdisplacements u to be consistent with a physically soundprediction up relying on Navier-Stokes equations. Asproposed in Heas et al. (2007)， we define this functional asa quadratic distance between the estimated field u and thedense propagated field up =(up，Vp)： where p denotes a weighting factor. This approach con-stitutes an alternative to the spatio-temporal smootherdefined in Weickert and Schnorr (2001) and is to someextent similar to the temporal constraint introduced inRhunau et al. (2007). It is important to distinguish thisa priori dynamical model from the model used for the data-term. Indeed， from the probabilistic point of view， mini-mizing the functional is analogouss to searching themaximum a posteriori estimate obtained by maximizing theproduct of the likelihood probability distribution repre-sented by the data-term and a priori probability distributionrepresented by spatio-temporal smoothing. 4.2 Dynamical model Dynamical models describing motion field evolution areneeded here for the prediction at time t + At of a motionfield u， using the previous estimation performed betweentime t- At and t. In Weickert and Schnorr (2001)， thetemporal derivative of the velocity vectors are constrained tobe weak， which is not consistent with fluid flow dynamics.The smoothing proposed in Rhunau et al. (2007) constrainsthe estimated vorticity to be close to a prediction. Relying onNavier-Stokes equations adapted to shallow flows， thepropagation model in Heas et al. (2007) extends this 2D caseto flows driven by a shallow water evolution law. In the present paper we chose to use the latter simplifieddynamical model for short time propagation of shallowfluid motion and adapt it to laser-sheet illuminated particleimage. Let us present this evolution model. We assumeincompressibility of the flow of constant densityy Po.Denoting respectively by F=(Fu，F)，T=(Tu，T)，the viscous forces and turbulent dissipation at subgridscales the horizontal momentum equations read： where the unknown quantities are the pressure function pand horizontal motion field u=(u， v). An integratedform of the momentum equations is produced by verticalintegration. It yields to the shallow water momentum equations of Saint-Venant (1871). However， a simplifiedmodel of these equations can be obtained when relying onthe approximation that the vertical profile of horizontalmotion is constant within the flow. In other words， weneglect horizontal motion vertical derivatives. By addingthe mass conservation model， we form a simplified shallowwater system which reads： where h denotes the depth function， here the thickness of thelaser sheet. Consequently， denoting by v and vs the kinematicviscosity and a subgrid turbulent viscosity， F and T arereplaced by diffusion terms in the above momentumequations. Let us denote the vorticity by = curl(u) andthe divergence by=div(u). The previous system may beexpressed in its vorticity-divergence form： where II is the determinant of the Jacobian matrix ofvariables (u， v)， vs=(CA，)1tl is the enstrophy-basedsubgrid scale model proposed by Mansour et al. (1978)，and C the Lilly’s universal constant equal to 0.17. For laser-sheet illuminated images of turbulent flows，the thickness of the laser sheet is constant in space andtime. Hence， the vorticity-divergence equations provide adynamical model independent of variable h. Based on therealistic assumption that in 3D flows the divergence of theapparent motion is not weak due to out of plane motions，the simplified vorticity-divergence model based on shallowflows reads： This model for shallow flows reducessininthe 2Dincompressible case to the vorticity equation： 4.3 Algorithm of resolution The curl and divergence completely determine the under-lying horizontal velocity field and the current velocityestimate can be recovered from these quantities up to alaminar flow. Indeed， denoting the orthogonal gradient byV-=(-0/0y，0/Ox)， the Helmholtz decomposition ofthe field into a sum of gradients of two potential functionsis expressed as where Uhar is a harmonic transportation part(div Uhar = curl uhar =0) of the field u and where thestream function Y and the velocity potential @ correspond tothe solenoidal and the irrotational part of the field. The latterare linked to divergence and vorticity through two Poissonequations. Expressing the solution of both equations asa convolution product with the 2D Green kernel Gassociatedd withhthe Laplacian operator：Y= G * k，①=G*， the whole velocity field can be recovered withthe equation which can be efficiently solved in the Fourier domain. Theharmonic transportation component uhar is recovered bysubtracting to the field u its solenoidal and irrotationalparts. Let us sum up this prediction process. The vorticity andthe divergence fields are developed in time in betweenconsecutive image frame using a discretized form of Eq. 17and time increments 8t. After each time increment， assuminguhar constant within each frame interval， Eq. 20 is used toupdate the velocity u needed by Eq. 17， with the currentvorticity and divergence estimates. To avoid instability， asemi-implicit first-order upwind time discretization schemeis used to integrate forward Eq. 17. Classical second-ordercenteredfinite difference schemesare used for thecurl and divergence discretization. To solve the linearsystem associated with the semi-implicit discretizationscheme， the matrix has been constrained to be diago-nally dominant， which is a sufficient condition for awell-conditioned inversion problem. This condition reads1/ot max，(lul +lvl-Il). Finally， the dynamical modeltime integration procedure results in the predicted horizontalmotion field u. The global algorithm to estimate dynamic consistentmultiscale motions is presented in Fig. 2 5 Numerical and experimental evaluation This section compares results on synthetic and real imagesequences obtained with the following approaches： Correlation To estimate correlation-based velocityfields， the commercial software DaVis 7.2 from LaVi-sion was used. A multipass algorithm with a finalinterrogation window size of 16 x 16 pixels and 50%overlapping was applied. Image deformation and roundGaussianweighting functionnwere used. Spurious.velocities3were identified withmedian filter andreplaced by the median. The vector fields obtained Fig. 2 Synthetic algorithm for dynamic consistent multiscale motionestimation with thiscorrelation technique： vwereuseddin thecollaborative approach. ●Dynamic consistent correlation (this paper) Thisapproach is a filtering process of the correlation via asecond-order regularization and a dynamic consistencysatisfying the Navier-Stokes equations. This physics-based post-treatment， obtained without the data term inEq. 12， provides a dense vector field (1 vector per pixel). Multiresolution (Corpetti et al. 2006) This fluid dedi-catedl optical ffllowapproachh uuses second-orderregularization and multiresolution to handle largedisplacements (see Sect. 2.2 for details). LDynamic consistent multiresolution (this paper) Thisapproach is an extension of the fluid dedicatedapproach0of Corpetti et al. (2006)) withh a spatio-temporal c constraint satisfyingg ttrhe Navier-Stokesequations. Correlation variational (this paper) This collaborativeapproach is an adaptation of Corpetti et al. (2006)approach to handle large displacements with correlationinstead of multiresolution. LDynamic consistent correlation-variational (this paper)As described above， this approach uses second-orderregularization， correlation to handle large displace-ments and dynamic consistency satisfying Navier-Stokes equations. For all approaches the same set of parameters was used，providing an optimal performance： α=1，000， Blo = 1with o=8 pixels and y/o=0.05. 5.1 Synthetic image sequence To evaluate the performance of the collaborative approachcompare to state of the art optical-flow and correlation approaches， a synthetic particle image sequence was gen-erated based on the direct numerical simulation (DNS) of2D turbulent flow (Carlier and Wieneke 2005).The presentflow contains typical difficulties for image-based mea-surements techniques， like high velocity gradients andlarge dynamic range. The Reynolds number based on thelength scale of domain of computation 2n and on the rms ofthe forcing velocity was equal to 30，000. Figure 3 presentsthe whole flow domain and the region considered for athorough comparison. The image sequence was generated with a home-madeparticle image generator involving comparable methodsthan those developed for the EUROPIV Synthetic ImageGenerator (Lecordier and Westerweel 2003). A velocity-vorticity formulation of the Navier-Stokes equations wasadopted for the DNS. The vorticity equation was solved inFourier space using desaliased Fourier expansions in twodirections with periodic boundary conditions. The timeintegration was third-order/three steps with a Runge-Kuttascheme. The code is called pseudo-spectral. The coordinatesof each particles were calculated in the physical space inside x [px] Fig. 3 Two-dimensional turbulence test case-colormap of thenormalized vorticity with superposed velocity field the vorticity equation computation using the Runge-Kuttascheme of the DNS in order to use the intermediate velocitycomponents. However， a bicubic interpolation of thevelocity field was applied since the particles were not sup-posed to be located at the nodes of the grid In the present evaluation the image pair considered ledto a maximum particle image displacement up to 6 pixelsThis particle image displacement was globally optimal tominimize the RMS error provided by the variationaltechnique. This relatively small maximum displacementnicecessitates a multiresolution pyramid with only two lev-els. However， the addition of noise in the image pairsenables to put variational techniques， involving multires-olution scheme， on the wrong track. In order to provide more realistic conditions the signal-to-noise ratio (SNR) was increased， simulating a reductionof the power of the virtual laser (Hain and Kahler 2007)The noise was quantified with peak SNR (PSNR) which isthe ratio between the maximum possible power of a signaland the power of corrupting noise that affects the fidelity ofits representation. Because many signals have a very widedynamic range， PSNR is usually expressed in terms of thelogarithmic decibel scale. The PSNR is most commonlyused as a measure of quality of reconstruction in imagecompression. The PSNR is defined as where d is the dynamic of the grey level of the signal， heretheeparticles. Foran image where theparticleaarerepresented over 8 bits d=2-1=255. MSE is themean squared error， which for two m× n images I， and Irwhere oneoftheeimages isconsidered：danoisy.approximation of the other， is defined as Typical PSNR values for images of good quality liebetween 30 and 40 dB. Figure 4 shows particle imagesinfluenced by the different levels of noise used for theevaluation. In comparison with these synthetic particleimages， typical real particle images recorded with a highspeed camera (Photron APX) exhibit relatively low PSNRvalues lying between 12 and 18 dB (see Fig. 5). 5.2 Real image sequence The collaborative approach was also applied to a particleimage sequence recorded in one of the wind tunnels of theRennes regional Center of the Cemagref. The sequenceshows the near wake flow of a circular cylinder at Reynoldsnumber Re= 3，900. Fig. 4 Synthetic particle images (50 ×50 pixels) influenced by thereduction of the power of a virtual laser. From left to right and top tobottom， PSNR values equal 36， 24，18 and 12 dB， respectively Fig.5 Particle image pair (50 x50 pixels) recorded with a highspeed camera (Photron APX) The circular cylinder had a length L= 280 mm and adiameter D= 12 mm. It was equipped with two thinrectangular end plates with the specification recommendedby Stansby (1974). The distance between the end-plateswas 240 mm providing an aspect ratio L/D=20. Theblockage ratio was4.3%. The circular cylinder wasmounted horizontally at 3.5D from the entrance of thetesting zone. The free-stream velocity was adjusted at4.8 m s . 2D2C PIV experiments were carried out with a LaVi-sion commercial system including a New Wave laser Solo 3Nd-YAG (Energy by pulse of 50 mJ) and 2 PCO camerasSensiCam (CCD size of 1，280× 1，024 pixels， pixel size of6.7 ×6.7 um’and dynamics of 12 bits). A lens with focallength of 50 mm and aperture of 5.6 was mounted on thecamera. The field of view was 3.6Dx 2.9D. The diameter of the particle seeding (diluted polyglycol in water) wereestimated to be less than 10 um. Given the image magni-fication (0.2)， lens aperture and light wavelength， theimages of the seeding particles are diffraction limited witha diameter of 1-2 pixels in diameter. The correlation-based velocity fields were calculatedwith the commercial software DaVis 6.2 from LaVision. Amultipass algorithm with a final interrogation window sizeof 16×16 pixels and 50% overlapping was applied.Image deformation and round Gaussian weighting functionwere used. Spurious velocities were identified with medianfilter， each vector component is checked independently andreplaced by the median following Westerweel (1994). This investigation of the near wake flow is described indetail in Parnaudeau et al. (2008). 5.3 Multiresolution limitations In Fig. 6 the total RMS error of the approaches is given asa function of the time for five image sequences character-ized by different SNR. The same velocity field was used togenerate the five image sequences. For all methods theRMS error was increasing with decreasing values of PSNR(increasing noise). For images of poor quality (i.ePSNR=12 dB)， the large error observed for the multiresolutionapproach Wassmainly explainedbyyttheinaccurate estimation of the large displacements. Figure 7clearly illustrates this behaviour. It shows the vector fieldof the exact solution for region of large displacements(larger than 2 pixels) and the colormap of the deviation ofoptical-flow from the exact velocity modulus. For regionsof large velocities the deviation from the exact solution waslarger than equal to 1 pixel. The corresponding RMS error，out of the range of the Fig. 6， was equal to 1.2 pixel. Smallparticles with large displacements in the images were PSNR Fig. 6 Effect of the noise on the total RMS error for the differentapproaches and during the image sequences x x [px] x Fig. 7 Colormap of the deviation of the multiresolution approachfrom the exact velocity modulus-vector field of the exact solutionfor displacements larger than 2 pixels. Top PSNR=12 dB， bottomPSNR = 18 dB smoothed out by the down sampling of the multiresolutionscheme. In these regions the real motion information waslost. As illustrated in Fig. 7， through the results of PSNRequal to 18 and 12 dB， this effect was further promoted bythe addition of noise in the images. The inaccurate estimation of the large displacements ledto a poor extraction of the dynamics especially at largescales. Energy spectra reported in Fig. 8 give a quick viewof the distribution of energy at different scales for differentlevel of noise in the images. For PSNR= 12 dB， themultiresolution approach completely failed to estimate thedynamics in the flow. For PSNR=18 dB and for largervalues of the SNR， the multiresolution showed consistentenergy distributions across the scales. Spectra of the error(difference between the multiresolution's estimation andthe DNS solution) are also shown in Fig. 8 indicating an Fig. 8 Spectral response of the multiresolution approach for differentlevels of noise in the images. Spectra ofthe error for the same data areshown in inset Fig. 9 Spectral response of different approaches for PSNR = 12 dB.Spectra of the error for the same data are shown in inset improvement of the estimations with decreasing noise inthe images. Inspection shows that the variational approachdoes not yields any windowing effect as it is observed forthe correlation in the form of a sinc function， visible atsmall scales (for k=1/16=0.0625). The monotonicallyvanishing error spectra provided by the variationalapproach clearly exhibit that the dense informationestimated was consistent with the DNS reference down tothe smallest scales. These behaviours can also be seen withthe results of the third international PIV Challenge for thevariational approach of Corpetti et al. (2006) and thecorrelation techniques (see Figs. 1-17 in Stanislas et al.2008) (Fig. 9). The accuracy of the small scales dynamics is directlylinked with the regularization involved to tackle the aper-ture problem. The regularization is conducted from the beginning of the minimization process and complement theinformation of the data term with a spatial coherence.Hence， this global variational technique extracts densevector fields (one vector per pixel)， beyond the densityassociated to the particles. We will see in the following thatthieeintroductionof physics-basedregularizationnccanimprove the estimation of the vortex dynamics at smallscales. Further informations about this subject can be foundin Heitz et al.(2008). 5.4 Correlation-variational approach As shown in Fig. 6 the combined correlation-variationalapproach resulted in lower RMS error. Compared to thecorrelation method， the proposed technique allowed theestimation of the large displacements with a lower devia-tion from the exact velocity modulus (see Fig. 10). Theimprovement that could be gained with the proposedapproach is exhibited in Fig. 9， giving the distribution ofthe energy estimated across the spatial scales. The errorspectra plotted in inset shows that the combined schemesmoothed out the window effect due to the correlation. Thebenefit of the combined technique was especially observedfor image sequences with large noise. However， for imageswith less noise， the colormap of the deviation from theexact velocity modulus plotted in Fig. 11， indicates that theproposed approach yielded a finer description of the vortexfilaments. In the present study， as given by the parameterswe manually set for best performance， the contribution ofthe data term of the variational approach was poor and theproposed correlation-variational global scheme behavedmainly as a regularizer of the correlation vector fields，denoising the vector fields and giving an enlargement ofthe dynamics. As already discussed above， the second-order variational technique used in this study gives resultscomparable with the most performing correlation tech-niques. Better performances， i.e. enlarged dynamic rangewith high accuracy， could be provided by the global vari-ational scheme proposed by Yuan et al. (2007)， usingmimetic finite difference to compute second-order regu-larization， or other physics-based regularization techniques(Heitz et al. 2008). The combined variational-correlation scheme was alsoevaluated with real images of turbulent flows. Figure 12shows the estimated velocity fields using variationaltechnique with multiresolution and the proposed correla-tion-variational approach. The multiresolution methodunder-estimated the external high velocity which wasapproximately two times smaller than the ground truth.This large error was reduced significantly by using theproposed collaborative approach. It should be noticed thatin Corpetti et al. (2006) the same multiresolution techniquesucceeded in estimating the fluid motions in comparable xa x [px] x [px] x [px] Fig.10 Deviation of different approaches from the exact velocitymodulus and corresponding， vector field-from top to bottommilitarisation， correlation， correlation-variational， dynamic consistentcorrelation-variational. PSNR = 12 dB xa x [px] x xa x [px] xa x [px] xa x [px] Fig. 11 Deviation of different approaches from the exact velocitymodulus and corresponding vector field-From topp to bottommultiresolution， correlation， correlation-variational， dynamic consis-tent correlation-variational. PSNR =36 dB x2 -2 488x [px] 6 x2a -2 x[px] Fig. 12 Instantaneous vector field with horizontal velocity colormap， in near the wake of a circular cylinder at Re=3，900. Topoptical-flow approachh((Corpetti et al.2006)， bottom proposedcollaborative approach wake flow. However， in that experiment the observedparticle diameters were around 3-4pixels， whereas in thepresent case they were around 1-2 pixels. In Fig. 13， thevorticity colormaps of correlation-based and combinedapproaches are shown. Although the vorticity was com-puted on the finest grid provided by each methods (8 pixelsfor the correlation and 1 pixel for the collaborative tech-nique)， less noisy vorticity maps are obtained by thecombined scheme， associated with a finer description of theflow. This evaluation with experimental data indicates thatthe proposed collaborative approach yielded reliable resultswhere classical optical flow techniques failed， with anenlargement towards small spatial scales compared tocorrelation-based technique. 5.5 Spatio-temporal regularization As shown in Fig. 9 for images with large noise， theimprovement of the estimation with the addition of thedynamic consistency is quite efficient. The enlargementof the dynamic range resolved was almost of one decade.In addition， the computational cost associated to thisspatio-temporal regularization was rather low (four timesless) compared to the computation cost of the correlation-variational step. However， this computational cost couldbe larger if the time separation between the images was 488x [px] Fig. 13 Instantaneous vector field with vorticity colormap， in thenear wake of a circular cylinder at Re =3，900. Top correlationapproach， bottom proposed collaborative approach Fig. 14 Spectral response of the dynamic consistent correlation-variational approach during the image sequence for PSNR= 12 dB.Spectra of the error for the same data are shown in inset increased， like in some real-time PIV experiments. Thedecay of the RMS observed in Fig. 6， indicated that theproposed approach was able to enhance the estimationsover several frames， the main improvement being pro-vided from the first to the second frame. During thefollowing time steps the spatio-temporal regularizationenhanced the descriptions of the spatial scales rangingbetween 30 and 10 pixels， corresponding to the size of the vortices shown in Fig. 3. This phenomenon was alsoobserved in the distribution of spectra during the timesequence， plotted in Fig. 14 and in the low deviationfrom the exact velocity given in the Fig. 10 where lessspuriousvectors can be seen compared to the otherapproaches. 6 Conclusion In this paper， we have proposed and evaluated a newmethod for robust and dense estimations of instantaneousvelocities of fluid flows from image sequences.. Thismethod is a collaborative approach satisfying Navier-Stokes equations and combining variational methods withcorrelation technique. The novel scheme has been specifi-cally designed to provide physics-based， robust over noiseand dense estimations. For laser sheet illuminated particleimages of 3D turbulent flows a dynamic model has beenderived from Navier-Stokes equations following shallowflow assumptions. The developed approach has been tested on syntheticparticle images based on 2D turbulence and on realimage sequence recorded in the wake of a circular cyl-inder. In each cases， we compared our results with theones issued from PIV and optic-flow.It was pointed outthat for image sequence with large noise the collabora-tiveschemere-enforced fluid dedicated variationalmethods towardsSrobustness， enlarging the range ofscales that can beresolved withcorrelation-basedtechniques. AcknowledgmentsThe financial support by the Region of Bretagnein France under Grant No.20048347 and by the European Unionunder Grant No. FP6-513663 are gratefully acknowledged. References ( Alvarez L ， Castano C， Garcia M ， K rissian K， Mazorra L， Salgado A， S anchez J ( 2 008) Variational s e cond or d er fl o w estimation for P IV s equences. E xp F l uids 44(2)：291-304 ) ( Bergen J， Anadan P ， Anna K， Hingorani R ( 1 992) H i erarchicalmodel-based m otion estimation. In： Sandin i G(ed) Proceedingsof E uropean conference on c o mputer vi s ion. S p ringer， Heidel- berg， pp 237-252 ) ( Bigiin J， Granlund G (1988) Optical f low based on the inertia matrixin the frequency domain. In： Proceedings of SSAB symposiumon p icture processing. Lund， S weden ) ( Black M ， Anandan P ( 1996) The rob u st est i mation of multiple motions： p arametric a nd piecewise-smooth flow fields 63(1)：75-104 ) Bruhn A， Weickert J， Schnorr C (2005) Lucas/Kanade meets Horn/Schunck： combining local and global optic flow methods. Int JComput Vis 63(3)：211-231 ( Carlier J， Wieneke B ( 2 005) Report 1 o n pr o duction and dif f usion offluid mechanics i mages and d a ta， Fluid project deliverable 1 .2. ) European Project ‘Fluiddimage analisys anddescription’(FLUID). http：//www.fluid.irisa.fr/ Corpetti T， Memin E， Perez P (2002) Dense estimation of fluid flows.IEEE Trans Pattern Anal Mach Intell 24(3)：365-380 Corpetti T， Heitz D， Arroyo G， Memin E， Santa-Cruz A (2006) Fluidexperimental flow estimation based on an optical-flow scheme.Exp Fluids 40：80-97 Cuzol A， Hellier P， Memin E (2007)A low dimensional fluid motionestimator. Int J Comput Vis 75(3)：329-349 ( Hain R ， K ahler C (2007) F u ndamentals of multiframe particle image velocimetry (PIV). Exp Fluids 42：575-587 ) ( Heas P ， Memin E，Papadakis N ， Szantai A (2007) Layered estimationof atmospheric mesoscale dynamics f r om s a tellite imagery. IEEE Trans Geosci Remote S e ns ( i n press) ) Heitz D，Memin E， Schnorr C (2008) Fluid flow motion measure-ments from image sequence： a review and some perspectives.ExpFluids (submitted) Horn B， Schunck B (1981) Determining optical flow. Artif Intell17：185-203 Kohlberger T， Memin E， Schnorr C (2003) Variational dense motionestimation using the Helmholtz decomposition. In： Griffin L，Lillholm M (eds) Scale space methods in computer visionLNCS， vol 2695. Springer， Berlin， pp 432-448 ( Lecordier B， Westerweel J (2003) The EUROPIV synthetic imagegenerator (SIG). In： Stanislas M， Westerweel J ， Kompenhans J(eds) Particl e imag e velocimetr y ： recent improvements.EURO- P IV 2 workshop， Springer， p p 145-162 ) Liu T， Shen L (2008) Fluid flow and optical flow. J Fluid Mech (inpress) ( Lucas B， K a nade T (1981) An i terative ima g e registration tec h nique with a pplication t o s tereo v ision. I n： P roceedings o f s e venthinternational joint conference on artificial i ntelligence， V ancou- ver， Canada， pp 674-679 ) ( Mansour NN， F e rziger JH， Reynolds W C (19 7 8) Lar g e-eddysimulation of a t urbulent mixing layer. Technical report. R eportTF-11， Thermosciences D ivision， Department of MechanicalEngineering， Standford University ) Okuno T， Sugii Y， Nishio S (2000) Image measurement of flow fieldusing physics-based dynamic model. Meas Sci Technol 11：667-676 Papadakis N， Corpetti T， Memin E (2007) Dynamically consistentoptical flow estimation. In： IEEE international conference oncomputer vision， ICCV’07. Rio de Janeiro， Brazil Rhunau P， Schnorr C (2007) Optical Stokes flow estimation： animaging based control approach. Exp Fluids 42：61-78 Rhunau P， Stahl A， Schnorr C (2007) Variational estimation ofexperimental fluid flows with physycs-based spatio-temporalregularization. Meas Sci Technol 18：755-763 Ruhnau P， Kohlberger T， Schnorr C， Nobach H (2005) Variationaloptical flow estimation for particle image velocimetry. ExpFluids 38：21-32 Saint-Venant A (1871) Theorie du mouvement non-permanent deseaux， avec application aux crues des rivieres et l'introduction desmarees dans leur lit. C R Acad Sci Paris 73：147-154 Stanislas M， Okamoto K， Kahler CJ， Westerweel J， Scarano F (2008)Main results of the third international PIV challenge. Exp Fluids45(1)：27-71 Stansby PK (1974) The effect of end plates on the pressure coefficientof a circular cylinder. R Aeronaut J 78：36-37 Suter D (1994) Motion estimation and vector splines. Seattle， USA，pp 939-942 ( Tikhonov AN， Arsenin VY (19 7 7) Solutions of i ll-posed problems. In： Scripta series in mathematics， chap II . ISBN 0470991240 ) ( Weickert J， Schnorr C (2001) Variational optic f low computation w i th a spatio-temporal smoothness constraint. JMath Imaging Vis 14(3)：245-255 ) ( Westerweel J ( 1994) Efficient detection of spurious vectors in particleimage velocimetry data sets. E xp Fluids 1 6：236-247 u1d ) ( Yuan J， Schnorr C， Memin E ( 2 007) D i screte o rthogonal decompo-sition and variational fluid flow estimatio n .J Math Imaging Vis28：67-80 ) Springer