On the stability of iterative PIV image interrogation methods by F. Scarano Delft University of Technology Department of Aerospace Engineering PO Box 5058, 2600 GB Delft; The Netherlands E-Mail: f.scarano@lr.tudelft.nl ABSTRACT The stability of cross-correlation based PIV iterative interrogation methods is investigated. The iterative image interrogation equation is formulated and two different options are considered for the implementation of the predictorcorrector principle. The first is based on the simplest approach, which sums up the intermediate results of the crosscorrelation; the second applies the displacement predictor as the average value over the interrogation window. The numerical assessment is performed making use of computer-generated one-dimensional particle recordings where a known displacement distribution is applied. The difference between the result of each interrogation and the exact displacement is monitored as measurement residual. The convergence of the iterative process is monitored through the difference between subsequent results. The problem relevance relies in the broad adoption of iterative interrogation methods by the scientific community as well as by PIV systems manufacturers. Methods applying continuous window deformation are discussed. In general terms the stability analysis of predictor-corrector algorithm shows that an integral condition of the predictor must be fulfilled in order to obtain a stable and converging process. However in the latter case the spatial resolution is that of a moving average filter. Finally an analytical justification of the phenomenon of the most amplified unstable wavelengths is given, followed with a numerical assessment. 1. INTRODUCTION Since their appearance, iterative image interrogation methods have widely spread over the PIV community. Thanks to the availability of affordable computer power, researchers found out that iterative interrogation techniques can significantly improve the analysis of digital PIV recordings. Among the most important advantages is the increase of dynamic range when the multi-step interrogation is coupled with a multi-grid interrogation with window shift (Soria 1996, Westerweel et al. 1997; Scarano and Riethmuller, 1997). The introduction of the iterative sub-pixel window shift technique (Lecordier, 1997) showed that an important improvement could be obtained also in terms of measurement sub pixel precision. Finally, the introduction of the iterative window deformation method (Huang et al. 1993, Jambunathan et al. 1995) opened the way to a new generation of image processing methods able to maximise the number of particle image pairs even in flows with a large velocity gradient. The most important improvement brought to the technique was in terms of the velocity gradient dynamic range. However the properties of the iterative interrogation process were not known and the interrogation required a user determined criterion to end the processing. Fincham and Spedding (1997) discussed the accuracy of the iterative sub-pixel window shift technique and that of iterative window deformation as compared to the conventional (one-step correlation) interrogation method, concluding that the sub-pixel window shift compensates for the bias error term, while the window deformation acts mostly in a reduction of the random error. Scarano and Riethmuller (2000) extended the error analysis to the velocity gradient with results in agreement with Fincham and Spedding, also verifying that the spatial resolution of the interrogation method obeys to the rules of linear filters theory. Nogueira et al. (2001) reported about the possible unstable behaviour of iterative cross-correlation interrogation when associated with image deformation. In a recent discussion during the worldwide PIV Challenge II, the issue of the instability of iterative image deformation methods was raised repeatedly due to the large number of researchers and system manufacturers adopting the iterative window deformation interrogation method, which motivates the need for the present study. The first part of the paper is devoted to the analysis of the iterative process corresponding to cross-correlation predictor-corrector interrogation. The model is then applied to the most common schemes to displace and deform the interrogation windows returning the theoretically predicted behaviour of the residual error term. The assessment of the 1 proposed theory is performed by means of synthetic PIV recordings. The analysis is extended to the effect of vector data low-pass filtering. 2. ITERATIVE IMAGE PROCESSING The iterative analysis of PIV images to evaluate the tracers images displacement V = (u,v) can be schematically represented by the above equation:
 k
 k −1  k V =V p +r (1) Where Vk indicates the result of the evaluation at the kth iteration. The right hand side terms represent the displacement predictor Vp = (up,vp) from the previous iteration. The correction term r = (ru,rv) is the vector determined by means of cross-correlation analysis of the images after a proper transformation of coordinates has been applied. The correction (viz. residual) term is obtained from the cross correlation of the transformed images I ak and I bk obtained from the original images I a and I b according to the relation: I ak ( x, y ) = I a ( x − udk / 2, y − vdk / 2) (2) I bk ( x, y ) = I b ( x + udk / 2, y + vdk / 2) The term Vd = (ud,vd) is the image deformation displacement commonly obtained interpolating the predictor displacement onto a pixel grid. When the image deformation is preformed by a piecewise constant function, then we speak of window offset (by integer number of pixels or sub-pixel). If a polynomial base is used to obtain the image deformation field, then we enter the domain of window or image deformation methods. The most common choice for the displacement interpolation is a bi-linear interpolation (Huang et al., 1993, Jambunathan et al., 1995, Scarano and Riethmuller, 2000). Higher order functions have also been considered (Fincham and Delerce, 2000; Lecordier, 2003), which may in principle improve the accuracy of image deformation. However the topic will not be discussed since it goes beyond the scope of the present paper. In the cases of sub-pixel window shift and window deformation, the images have to be interpolated for the re-sampling at non-integer pixel position, which also opens another discussion about the optimal schemes for accurate image reconstruction. The analysis of the deformed images is assumed to return the (moving) average displacement VMA(x0,y0) between the two images over the selected interrogation area W. Therefore, for a given location P = (x0,y0) the expression of the residual displacement measured between the transformed images is:
1 r ( x0 , y0 ) = 2 l 

 V x y − V , d ( x, y )  dxdy ( ) ∫∫W   (3) Where the integral is extended over the interrogation window W of linear size l. The displacement Vd is then obtained from the integration of the chosen base function over the interrogation window. In case of simple window shift Vd coincides with the value in the middle Vmeas. If the transformation of the images is performed by means of a continuous function (linear or higher order interpolation) it can be demonstrated that Vd differs from Vmeas by an amount proportional to the window size and to the displacement second derivative (Scarano, 2003). The iterative equation can finally be rewritten as:
 k
 k −1

 k −1 1 V ( x0 , y0 ) = V p ( x0 , y0 ) + 2 ∫∫ V ( x, y ) − V d ( x, y )  dxdy      l W  current result predictor (4) corrector While the corrector term can be directly evaluated by cross-correlation of the deformed images, the predictor term has not been yet discussed with enough details. The most straightforward choice is to set the predictor equal to the value obtained from the previous interrogation (assuming that the iterative interrogation is performed on a constant grid), therefore the first choice reads as: 2 
 k −1
 k −1 V p ( x0 , y0 ) = V meas ( x0 , y0 ) (5) “local predictor” This approach appears very logical at a first sight, and the application is straightforward in a PIV interrogation algorithm, which probably justifies its broad adoption in the PIV community (Lecordier 2003, Nogueira et al. 2003, Wieneke 2003, Willert 2003, Wereley 2003). However it will be shown that the result of eq. 4 is intrinsically unstable and the iterative approach may diverge unless the image processing is interrupted at an early stage. 2.1 The integral equation for iterative correlation Previous studies showed that the result of cross-correlation could be modelled with a moving average linear filter (Scarano, 2003). The result to be expected at the end of the cross-correlation process is given by an expression similar to eq. 5.
 1 V meas ( x0 , y0 ) = 2 l
 ∫∫ V ( x, y ) dxdy (6) W The measured velocity at a given location P = (x0,y0) can be modelled as the window average value rather than the value at the centre of the window of the local velocity distribution. The consequence of this assumption is that the spatial resolution of the cross-correlation interrogation follows the same rules as linear filters (i.e. moving average). If the predictor term is rewritten to take into account the integral form of the iterative interrogation:
 k −1 1 V p ( x0 , y0 ) = 2 l ∫∫ W
 k −1 V meas ( x, y ) dxdy “integral predictor” (7) And inserted in the integral equation for iterative correlation, the result reads as:
 k
 k −1

 k −1 1 1 1 V ( x0 , y0 ) = 2 ∫∫ V meas ( x, y ) dxdy + 2 ∫∫ V ( x, y ) − V d ( x, y )  dxdy = 2   W W      l l l     current result predictor ∫∫ W
 V ( x, y ) dxdy (8) corrector The above equation is consistent with the definition of measurement by integral given in eq. 6, when the displacement field used to deform the images Vd is the same as that used to predict the displacement Vmeas from the previous iteration. 2.2 Stability of iterative interrogation An iterative process is defined in this context as convergent when a norm of the difference between the current result and the previous one (iteration residual |Vk – Vk-1| ) decreases monotonically. A convergent process does not necessarily lead to the exact value of the velocity. Therefore if the asymptotic value is that of the exact velocity we speak of absolute convergence Vmeas→V. Otherwise we denote it as relative or biased convergence. The iteration residual may follow in some cases an oscillatory behaviour. If the amplitude of the oscillations increases then the process is considered divergent or unstable. The result could be considered converged when the iteration residual falls below a given threshold (e.g. 1% of the measured velocity). 3. NUMERICAL ASSESSMENT The properties of the different iterative schemes are assessed and compared by means of computer simulated one-dimensional PIV recordings with a known displacement distribution. The particle image properties and the seeding density are chosen close to optimal level (particle image diameter dp = 2 pixels, particles density ppp = 0.5 particle per pixel) in order to avoid the interference of other effects in the simulation. The particle images are integrated over an array of length L = 4000 pixels, returning an intensity diagram as shown in the figure below. 3 400 300 Intensity 200 100 0 -100 -200 100 125 150 175 200 225 250 275 300 Position [pixels] Fig. 1. Sample of the one-dimensional particle image pattern (after mean subtraction). The one-component displacement field U(x) considered is a sinusoid with spatially varying wavelength and fixed amplitude of 2 pixels. 2 U ( x ) = 2sin ( 2π x / L ) ⋅ ( 2 x / L )  ;   x ∈ [1, L] (9) The reference result is taken as the output of a moving average top-hat filter UMA(x) with a filter length l equivalent to the interrogation window size. In the present case the interrogation window size is l = 201 pixels and an overlap factor of 75% between neighbouring interrogation windows is applied. Fig. 2 shows the result of the iterative interrogation technique based on the application of eq. 2-5. The first iteration returns the result (top-left) as commonly done by means of conventional cross-correlation. It can be seen that the measured displacement follows approximately the moving average of the applied displacement. After the first iteration, the measurement approaches more closely the applied displacement and it keeps close to it with a small overall error during several iterations. However proceeding in the iterative analysis the measured displacement shows an increasing noise level amplified iteration after iteration. For instance, at the twentieth iteration, noisy fluctuations have become clearly visible, and the velocity measurement is significantly degraded. Moreover the presence of this type of noise compromises even more seriously the measurement of the velocity spatial derivative. It may be noticed that the error grows more rapidly in the points where the curvature of the velocity profile is highest, which will be explained further in the present study. 4 2 2.5 2 1.5 1.5 1 0.5 U [pixels] U [pixels] 1 0 -0.5 0.5 0 -0.5 -1 -1 -1.5 -1.5 -2 0 1000 2000 Position [pixels] 3000 -2 0 4000 3 4 2 3 2000 Position [pixels] 3000 4000 1000 2000 Position [pixels] 3000 4000 2 U [pixels] 1 U [pixels] 1000 0 -1 1 0 -1 -2 -2 -3 0 1000 2000 Position [pixels] 3000 -3 0 4000 Fig. 2. Applied displacement field U(x) (solid black line), moving averaged displacement UMA (solid blue line) and measured displacement Umeas (solid red line and dots). Top-left: 1st iteration. Top-right: 2nd iteration. Bottom-left: 5th iteration. Bottom-right: 20th iteration. Fig. 3 describes the evolution of the measurement error and the iteration residual. The graph on the left side is related to the local predictor scheme according to eq. 5. It can be noticed that the measurement error decreases rapidly during the first four iterations and then it increases again. The iteration residual decreases for a larger number of iterations, reaching a minimum of 10-2 pixels. Fig. 4 displays the result obtained applying the integral form of the displacement predictor (eq. 7) in the iterative interrogation equation. One can clearly see that only a relatively small difference in terms of measurement error occurs between the second and the twentieth iteration. The measured displacement shows no oscillations and converges to the moving averaged displacement distribution. The convergence of the iteration residual within 1% is reached already at the 3rd iteration, however the measurement error does not reach values as low as in case of the local predictor scheme where a minimum below 10-1 pixels is obtained. This is essentially due to the lower spatial response of the integral method, which produces a larger absolute error in the high-frequency side of the diagram. 10 0 error 10 0 error Residual [pixels] Residual [pixels] 10 -1 residual 0 5 10 Iteration # 15 20 10 -1 10 -2 10 -3 10 -4 10 -5 0 residual 5 10 Iteration # 15 20 Fig. 3. Measurement RMS error (red solid line and +) and iteration residual (blue solid line and dots). Left: local velocity predictor (eq. 5). Right: integral velocity predictor (eq. 7). 5 2 1.5 1 1 0.5 0.5 U [pixels] U [pixels] 2 1.5 0 -0.5 0 -0.5 -1 -1 -1.5 -1.5 -2 0 1000 2000 Position [pixels] 3000 -2 0 4000 1000 2000 Position [pixels] 3000 4000 Fig. 4. Applied displacement field U(x) (solid black line), moving averaged displacement UMA (solid blue line) and measured at the first iteration Umeas (solid red line and dots). First iteration: see Fig. 2 top-left. Left: 2nd iteration. Right: 20th iteration. When the spatial derivatives of the displacement distribution are concerned, one observes that the effect of the unstable behaviour of the local predictor scheme is even more evident. The comparison of the two schemes after only four iterations is shown in Fig. 5. The diagram with solid black line is the exact spatial derivative of the signal. The data series with a solid blue line and dots represents the integral scheme and the red line diagram is the result from the local predictor scheme. One may conclude again that although the local predictor scheme follows better the exact values on average, however it exhibits a much larger data scatter, which makes the results more uncertain. 0.03 UDU/DX [pixels/pixel] 0.02 0.01 0 -0.01 -0.02 -0.03 0 500 1000 1500 2000 2500 Position [pixels] 3000 3500 4000 Fig. 5. Spatial derivative of the measured displacement. Solid black line: exact displacement derivative; red line and dots: measured displacement derivative at the fourth iteration with local predictor scheme; blue line and crosses: integral predictor scheme. A specific test is performed in order to explain the nature of the instability observed in the local predictor scheme. A null displacement distribution (initial noise level displacement 10-16 pixels) is considered. Starting from it, the iterative processing advances and the random correlation noise triggers the instability of the interrogation scheme, which acts amplifying selectively different wavelengths or a single wavelength. Fig. 6 shows the diagram of the measured displacement after 40 iterations, which corresponds directly to the measurement error since the applied displacement is null. In the specific case one can see that wavelengths of 2/3 l are amplified, corresponding to the situation where the local predictor and the average displacement over the window are in phase opposition. Denoting the error as ε(x) = U(x) – Umeas(x), we find that when ε(x) = asin[2π(x-1/3 l)/(2/3 l)], then ε(x0) = a. On the other side the measured displacement is given by U meas ( x0 ) = 1 l2 ∫ l 2 −l 2 ε ( x ) dx = −2a . Such conclusion is in agreement with the findings of Nogueira et al. (1999), who ascribed the instability to the sign reversal in the sinc shape of the response function of the top-hat window (Fig. 6-right). 6 When the above expressions are introduced in eqs. 4 and 5, it can be concluded that the right hand side term is increasing for any non-zero predictor displacement, which finally causes the growth of the unstable wavelength. 0.01 0.005 0 -0.005 -0.01 0 500 1000 1500 2000 2500 3000 3500 4000 Fig. 6. Measured displacement after 40 iterations. The initial random error is amplified at a selected wavelength of 2/3*l corresponding to the maximum of the first sign reversal in the frequency response shown on the right. 3.1 Stability vs. spatial resolution The spatial response of the different approaches to iterative processing is further assessed on a simulated PIV image (two-dimensional) of sinusoidal fluctuations with spatially varying wavelength. The normalised response is given as an indication of the modulation transfer function of the PIV technique. Moreover the level of the RMS fluctuations quantifies the measurement uncertainty associated with the different schemes. A computer-generated image of 1000×2000 pixels with a sinusoidal shear motion in the vertical direction is considered, which allows to draw direct quantitative conclusions on the properties of the measurement. The image is analysed with an interrogation window of l = 31 pixels and an overlap factor of 75%. The iterative scheme is applied over five iterations. Being the displacement constant along the y-coordinate, the measurement statistics can be evaluated (mean and RMS). The result is given in terms of mean displacement over the y-coordinate V(x) in Fig. 7-top. It can be clearly noticed that all methods are able to describe the fluctuations up to a minimum wavelength equal to the window size. Fig. 7-middle describes the amplitude ratio V/V0 versus the normalised window size l* = l/λ, where λ is the local wavelength of the sinusoid. The method based on the integral predictor once again returns the same result as that of a linear filter of equivalent size applied to the displacement field. The amplitude ratio drops below 50% at l* = 0.57. Conversely, the method based on the local predictor method, drawn in blue, shows the appealing feature of a significantly higher spatial response, which drops below 50% only at l* = 0.90. When the local predictor method is considered, however with the application of a smoothing operation between the iterations except for the last iteration, then the amplitude ratio drops below 50% at l* = 0.85. In terms of RMS uncertainty (Fig. 7-bottom), the local predictor method shows the highest noise level with a noise below 10-1 pixels up to l* = 0.4. When adopting the integral method or applying the local method with smoothing the noise level is kept almost one order of magnitude lower and stays below 10-1 pixels up to l* = 0.7. 7 V [pixels] 2 1 0 -1 -2 0 200 400 600 800 1000 X [pixels] 1200 1400 1600 1800 2000 0 0.2 0.4 0.6 0.8 1 l/lam 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 l/lam 1.2 1.4 1.6 1.8 2 V/V0 1 0.5 0 0 10 RMS(V) -1 10 -2 10 -3 10 Fig. 7. Results from the sinusoidal shear motion analysis. Top: applied displacement field (black); measurement with integral predictor scheme (red); measurement with local predictor scheme and smoothing except last iteration (green); measurement with local predictor without smooth (blue). Middle: amplitude ratio at sinusoid crests as a function of the normalised window size. Linear filter with window size (black) and half size (blue); colors same as above. Bottom: measurement RMS uncertainty as a function of the normalised window size. 4. CONCLUSIONS The study covered the aspect of stability of iterative cross-correlation analysis of PIV recordings from a theoretical and numerical viewpoint. It was shown that the predictor-corrector procedure based on the local value of the measured displacement leads to the growth of unstable waves. Opposed to that, the scheme that adopts the integral of the displacement over the interrogation window does not suffer from any type of instability and returns a measurement in good agreement with linear filters theory. Moreover convergence is reached almost immediately (three iterations are sufficient) with respect to the local predictor method, which reaches a minimum error after about 8 iterations and then starts to diverge. The major advantage of a local scheme is that it shows a higher spatial resolution, however at the cost of a higher noise level. An optimum is achieved if the unstable wavelengths are damped at each iteration except for the last iteration. In this case a high spatial response is also associated to the lowest measurement RMS uncertainty. REFERENCES Huang HT; Fielder HF; Wang JJ (1993) Limitation and improvement of PIV, part II. Particle image distortion, a novel technique. Exp Fluids 15: 263-273 Jambunathan K; Ju XY; Dobbins BN; Ashforth-Frost S; (1995) An improved cross correlation technique for particle image velocimetry. Meas Sci Technol 6: 507-514 8 Lecordier B; (1997) Etude de l’interaction de la propagation d’une flamme premelangee avec le champ aerodynamique, par association de la tomographie Laser et de la Velocimetrie par Images de particules. These de doctorat de l’Universite de Rouen Nogueira J; Lecuona A; Rodriguez PA (1999) Local field correction PIV: on the increase of accuracy of digital PIV systems. Exp Fluids 27: 107-116 Nogueira J; Lecuona A; Rodriguez PA; (2001) Identification of a new source of peak locking, analysis and its removal in conventional and super-resolution PIV techniques Exp. Fluids 30: 309–316 Scarano F; Riethmuller ML (2000) Advances in iterative multigrid PIV image processing. Exp Fluids 29: S051-S060 Scarano F; (2003) Theory of non-isotropic spatial resolution in PIV. Exp Fluids 35: 268-277 Soria J; (1996) An investigation of the near wake of a circular cylinder using a video-based digital cross-correlation particle image velocimetry technique. Experimental Thermal and Fluid Science 12: 221-233 Westeerweel J; Dabiri D; Gharib M (1997) The effect of a discrete window offset on the accuracy of cross-correlation analysis of digital PIV recordings. Exp Fluids 23: 20-28 Wieneke B; personal communication PIV Challenge II (2003) Willert C; personal communication PIV Challenge II (2003) 9 View publication stats