流体中速度场,速度矢量场检测方案(粒子图像测速)

H←V"f(co)，gVf(co)d←-(H+XI)-gchoose aE 0，1)， so that f(co+ad) ≤f(co)return co + ad A Variational Approach to Adaptive Correlationfor Motion Estimation in Particle ImageVelocimetry* Florian Becker， Bernhard Wieneke， Jing Yuan， Christoph Schnorr Image and Pattern Analysis Group， Heidelberg Collaboratory for ImageProcessing， University of Heidelberg， Germany，becker，yuanjing，schnoerr}@math.uni-heidelberg.de LaVision GmbH， Gottingen，Germany， bwieneke@lavision.de Abstract. In particle image velocimetry (PIV) a temporally separatedimage pair of a gas or liquid seeded with small particles is recordedand analysed in order to measure fluid flows therein. We investigate avariational approach to cross-correlation， a robust and well-establishedmethod to determine displacement vectors from the image data. A“soft”Gaussian window function replaces the usual rectangular correlationframe. We propose a criterion to adapt the window size and shape thatdirectly formulates the goal to minimise the displacement estimation er-ror. In order to measure motion and adapt the window shapes at thesame time we combine both sub-problems into a bi-level optimisationproblem and solve it via continuous multiscale methods. Experimentswith synthetic and real PIV data demonstrate the ability of our ap-proach to solve the formulated problem. Moreover window adaptationyields significantly improved results. Introduction Overview. Particle image velocimetry is an important measurement techniquefor industrial fluid flow questions. Small particles are introduced into liquidsor gases and act as indicators for the movement of the investigated substancearound obstacles and in mixing zones. A 2D plane is illuminated by laser lightrendering the particles in there visible to a camera which records two images ofthe highlighted area within a short time interval. The analysis of the image data allows to determine the movement of particlesand with this to measure the speed， turbulence or other derived mechanicalproperties of the fluid. In contrast to particle tracking velocimetry where firstsingle objects are identified by their position and then matched between twoimage frames， algorithms for PIV determine patches from the first and secondframe that fit best to some similarity measure. Cross-correlation has developedthe state-of-the-art method for motion estimation in PIV and benefits from itsrobustness against noise and illumination disturbances. ( * This work was p artially financed by the EC project FLUID (FP6-513663). ) In this paper we describe a variational approach to cross-correlation by con-tinuously optimising over the displacement variables. In addition a criterion isdefined to locally adapt the correlation window in order to improve accuracy ofthe estimation. Related Work and Contribution. A vast number of literature exists onall aspects of the application of cross-correlation for PIV， here we only referto [1] for an excellent overview. Typically an exhaustive search over the integerdisplacements is performed to search for the highest correlation peak. The corre-lation function is interpolated to gain sub-pixel accuracy. In contrast we presenta variational approach to motion estimation based on continuously maximisingthe cross-correlation between two images. The correlation window is described bya “soft”Gaussian weighting function instead of a sharp rectangular mask. Thisidea is used both in a local [2] and global context [3] for smoothing the opticalflow constraint. However while most approaches use a fixed window size com-mon for all positions， we formulate a sound criterion for the location-dependentchoice of the window shape parameters (size， orientation， anisotropy) in wordsof a further optimisation problem. Both displacements and window parametersare determined as a solution of a combined， bi-level minimisation problem whichis being solved via a multiscale gradient-based algorithm. We test our approachwith synthetic and real particle images to demonstrate the ability to robustlydetermine displacements and that window shape adaptation can improve resultssignificantly. Organisation. In section 2 we explicitly define the cross-correlation in wordsof a continuous minimisation problem and introduce the utilised weighting func-tion. Our approach to window adaptation is motivated and described. Discreti-sation and optimisation of the defined approach is subject to section 3. Resultsof experiments with real and synthetic data are given in section 4. We concludeand describe further work in section 5. 2 Problem Statement 2.1 Variational Approach to Correlation The input data consists of two images defined on 2 C R2. However we definethem to vanish outside and thus obtain two infinite image functions g1，92 ：RR. For the continuous case， we define the negative cross-correlation func-tion at position a E by with a window function w(r，2) parametrised by 2. In order to estimate themovement between two image frames in the considered areas， the correlation function is minimised with respect to the displacement vector v e R2. The localestimation is extended to a global variational problem to determine a vectorfield u：2→ R. where2is defined on and describes the location-dependent window shape.This formulation allows to add regularisation terms such as physical priors，e.g.incompressibility constraints [4]， on the vector field. However here the integrationover the correlation window in (1) is the only spatial regularising mechanism. In this work we jointly determine the integer and fractional part of the dis-placement by continuously searching for an optimum of the correlation function.In addition we choose a “soft” window function w(z，) ：= exp(-c2-l c)which is basically a non-normalised Gaussian function， instead of a sharp， rect-angular window. The symmetric， positive definite two-by-two matrix 2e Sallows to continuously steer the size， anisotropy and orientation of the window，see Fig. 1 for some possible shapes. j\\-/， Fig.1. Possible shapes of the weighting function： varying size， anisotropy， orientation 2.2 Window Adaptation When cross-correlation is employed for estimating motion it is implicitly assumedthat the displacements within the considered window are homogeneous. Howeverthis only holds true in very simple cases and leads to estimation errors in areasof large motion gradients as the vector field is smoothed out. This effect could beavoided by reducing the window size， however with the harm of a smaller areaof support and number of particles and thus a higher influence of image noise. In order to improve accuracy we propose to adapt the window shape byminimising a function which models the expected error subject to the choice ofthe window parameter 2 at position aE . Given a fixed vector field u wedefine the energy function as The first term of (3) measures the deviation from the assumption u= const. Inaddition we assume a constant error eout if the correlation window incorporatesdata from outside the image domain. The second term describes the error caused by insufficient large supportfor the displacement estimation in the presence of a homogeneous movement.We assume that each vector u results from a weighted least-square estima-tion u= arg minufw(z，2)u-u(z) dr over independent measurements uof the true displacement u*， which are disturbed by Gaussian additive noise，i.e. u(z)~N (u*，oI). Then it is possible to show that the expected square er-ror is Ellu-u*l}=Jw(r，2)drdw 2mvdet ·The parameter o in E constitutean estimation of the influence of image noise on the measurement error. Note that our definition of the error measure can easily be extended to involvefurther expert knowledge about the local influence of experimental parameters，such as particle seeding density， on the error of the cross-correlation method. 2.3 Joint Optimisation In section 2.1 and 2.2 we proposed two concepts disregarding their dependenciesby assuming that the window shape is fixed during correlation respectively thata vector field is known to estimate the error caused by spatial displacementvariations. Now we combine both by defining a bi-level optimisation problem， The top-level optimisation estimates the displacements u at all positions c in theimage by maximising the correlation terms. The window shapes are adapted inthe underlying optimisation problems that constrains each 2(r) to a minimumof the error estimation function E and again depends on u. 3 Discretisation and Optimisation 3.1 Data and Variable Discretisation The discrete input data g1，g2 is assumed to be sampled at a regular grid Y withgrid spacing ay and is stored in a cubic spline representation which yields a twotimes continuously derivable representation. We use an efficient implementationbased on [5] to evaluate the function g， its gradient Vg and its second deriva-tives Vg also at non-integer positions. Grey values of gi and g2 are shiftedbeforehand so they have a mean of zero each. Values outside the image domainare defined to be zero. Displacement and window shape variables are discretised on a separate reg-ular grid X C with spacing a which is typically chosen to be coarser than Y.We denote the variables located at the coordinates siE X as ui ：=u(ci) re-spectively 2i：=(ai). For the discretisation of the integral expressions in (2)and (3) we use simple， piecewise constant element functions which incorporate a discrete basis for functions R2R if arranged on a regulargrid with spacing a. The discrete versions of the objective functions then read The derivative of the first function with respect to ui simplifies to For speed up evaluation is restricted to a rectangular area that encloses allpositions with w(c，2)≥10-. 3.2 Optimisation The joint formulation (4)-(5) is a bi-level optimisation problem with non-convexfunction in both layers. First we remove the explicit constraints 2； > 0 of thesecond level by including them into the objective function as a logarithmic barrierfunction and define with some constant p>0， here chosen as p=10-2. In the same manner itis possible to incorporate additional upper and lower bounds on the windowsize. We relax the optimisation objectives by considering only the first orderoptimality conditions， An initial solution is iteratively improved by updating the variables in parallel.initialise 2(0)， u0)，k+0 repeat k←k+1 until convergence The update step is a Levenberg-Marquardt method that just as the Newton stepmethod involves first and second order derivatives of the objective function f(a)，but additionally steers the step length by the parameters 入> 0. function update(f(c)，co，入) However due to the non-convexity of Cand E， a variable assignment that fulfilsequations (6)-(7) does not necessarily incorporate a local optima for (4)-(5). Forthis reason we line search along the step direction to reduce the value of theobjective functions and to avoid local maxima and saddle points. 3.3 Multiscale Framework The optimisation method described in the previous section finds a local optimumwhich however might be far from a globally optimal solution. Our approachavoids most of them by plugging the single-level optimisation into a multiscaleframework. Figure 2 illustrates how this allows to circumnavigate the many sub-optimal positions in the correlation function of two noisy images. 15 Fig. 2. Optimisation of the non-convex correlation function of noisy image data： Valueof -C(u，2，x) over u and evolution of the displacement over several iterations (brightline， starting in (0，0))； due to the multiscale framework the method does not get stuckin a local optima but finds the correct solution in (+8，-8). In order to gain a coarse-to-fine representation of the image frames， the datais repetitively low-pass filtered and sub-sampled. When moving in opposite direc-tion within the resolution pyramid， the variables are sampled to a finer grid and used as initialisation for the next iteration steps. We use cubic spline interpola-tion both for the re-sampling of data and vector variables. For the window shapevariables， bi-linear interpolation implicitly conserves the constraint 2e S2. 4 Experiments and Discussion 4.1 PIV Evaluation Data Set Our first experiment is based on an evaluation data set designed to verify theability of correlation-based motion estimation methods to cope with large gradi-ents in the vector field. Case A4 of the PIV Challenge 2005 6 contains an areanamed 1D Sinusoids I which consist of two synthetic particle images each 1000by 400 pixels in size. Their vertical displacements vary sinus-like with differentwavelength (400 pixels on the left down to 20 on the right) while the horizontalcomponent is zero everywhere， see Fig. 3(a)-(b). Six respectively two scale levels were used for the fixed-window experimentsand experiments with window adaptation. The scale factor between two suc-cessive levels is √2. The parameters in (3) were set to a=20 and eout =10.Windows were constrained so that their 50%-level curve lies within a radius ofabout 63 pixels，i.e.w(x，2)<0.5 for all ||c|l2≥63. The maximum displacementin data is about 2.7 pixels and 1.2 pixels in average. The results (see Fig. 3(c)) show that pure correlation with fixed windowshapes can capture the main structures of the images but fails to accuratelyestimate the vector fields especially in the presence of large displacements gra-dients. However when used as initialisation to adaptive correlation we obtaineda precise reconstruction (see Fig. 3(d)) of the displacement field， even for thesmallest wavelength. 4.2 Real Data In order to test the ability of our approach to cope with noise and disturbancesin real-life turbulent data we applied it to an image pair from a wind tunnelexperiment (wake behind a cylinder， see [7]). Figure 4(a) shows the resultingvector field. The multiscale framework used eight scale levels and a scaling factor of √2.Window adaptation parameters were chosen as o= 10 and eout = 10. Thewindow radius was constrained to a maximum of about 32 pixels. The averagemeasured displacement is about seven pixels. Although we used neither regularisation terms nor data validation algorithms(e.g. median filter) on the vector field we obtained a smooth solution withoutoutliers. Figure 4(b)-4(d) demonstrate how the window shapes align themselvesalong areas of equal displacements and avoid gradients as intended by the designof the window adaptation criterion. Fig.3. Experiments with synthetic data： (a) Groundtruth vector field (sub-sampled).(b) Vertical component (mapped to grey-values： bright =up， dark = down) of thegroundtruth vector field and highlighted detail， also shown on the right. (c) Correla-tion with fixed window shape； estimates inevitably deteriorate at small wavelengths.(d) Joint correlation and window adaptation can significantly improve accuracy despitespatially variant wavelengths. (d) detailed view (upper right) Fig. 4.Experiments with real data (wake behind a cylinder)： (a) Results (sub-sampled)of the correlation approach with window adaptation. The dark highlighting marks thearea of the (b)-(d) detailed views of the vector field with some adapted windows (con-tour line： w(，2)=0.5). Note that each window propagates into regions of homoge-neous movement and perpendicular to gradients in the vector field and not necessarilyalong the flow. 5 Conclusion and Further Work Conclusion. We proposed an approach to fluid flow estimation based on thecontinuous optimisation of the cross-correlation measure. The expected esti-mation error for the choice of the correlation window shape is modelled andminimised in order to adapt the windows to displacement gradients. Both sub-problems were combined in a bi-level optimisation problem. A multiscale gradi-ent-based approach was described that continuously searches for both optimaldisplacements and window shapes. Experiments with synthetic and real data showed that the approach can copewith large displacements and disturbances typical for real fluid flow experiments.The adaptation of window shapes by means of the error expectation model leadsto meaningful results in the presence of displacement gradients and improvederror significantly in the PIV-Challenge data set. Further Work. Our approach is the origin for further potential investigations：Due to its variational formulation it is possible to extend the correlation termto estimate also affine displacements within the window and to involve physicalpriors， such as incompressibility. Further expert knowledge and statistical information can be incorporatedinto window adaptation criterion to improve the estimation accuracy. Also morecomplex shapes for the weighting function should be considered. A comparisonto state-of-the-art correlation implementations is planned. References 1. Raffel， M.， Willert， C.， Kompenhans， J.：Particle Image Velocimetry. Volume 2.Springer-Verlag， Berlin (2001) 2. Lucas， B.， Kanade， T.： An Iterative Image Registration Technique with an Applica-tion to Stereo Vision. In： Proceedings of DARPA Imaging Understanding Workshop.(1981) 121-130 ( 3. Bruhn， A.， Weickert ， J . ， Schno r r， C .： Lucas/Kanade Meets Horn/Schunck ： Com-bining Local and Global Optic Flow Methods. I nternational Journal of ComputerVision 61(3)(2005) 211-231 ) ( Ruhnau， P.， Stahl， A .， S chnorr， C .： Variational Estimation of Experimental FluidFlows w ith P hysics-based S patio-temporal Regularization. Measurement Scienceand Technology 1 8 ( 2007) 755-763 ) ( 5. U L nser， M.， Aldroubi， A .， Eden， M.： B -Spline S ignal P rocessing： Part II-EfficientDesign and Applications . I IEEE Transactions on S i gnal P r ocessing 4 1 ( 2 ) (1993)834-848 ) 6. Stanislas， M.， Okamoto， K.， Kahler， C.， Westerweel， J.： Third International PIV-Challenge.(2006， to be published in Exp. in Fluids) ( 7. Carlier， J . ， Wieneke， B.：Deliverable 1.2： Report On P roduction and D i ffusionof Fluid Mechanic Images and Data. Activity Report， European Project FLUIDDeliverable 1 .2， Cemagref， LaVision (2005) )