Performance analysis and application of grid interpolation techniques for fluid flows

Although it is common for automated image processing techniques to claim subpixel accuracy in the identification of particles, or centroids of displacements of groups of particles, additional errors are inevitably introduced when and if these data are reinterpolated back onto a grid mesh whose nodes lie at different locations from the original data. Moreover, these errors can be large compared to the errors introduced in the original image processing step.Two different techniques, convolution with an adaptive Gaussian window (AGW), and a two-dimensional thin-shell spline (STS), have been compared and contrasted for interpolating irregularly spaced data onto a regular grid. Both techniques are global interpolators; the Gaussian kernel applies an ad hoc choice of smooth function, while the thin-shell spline minimises a global functional proportional to the Laplacian of the velocity field. In this way, the smoothness constraint on the spline coefficients may be thought of as akin to a viscous smoothing of the fluid flow.Performance curves are given, enabling the investigator to make an informed choice of interpolating routine and grid interpolation parameters to minimise the interpolation errors, given various external constraints. Some illustrative example applications on real experimental data are described. In general, the importance of matching the interpolation technique to the characteristics of the original data is stressed. It is also pointed out that a correct interpretation of grid interpolated data must be based on a basic knowledge of the performance characteristics of that interpolator. Finally, recommendations are made concerning the development of surface spline techniques for problems involving large numbers of data points.


Introduction
The design and application of image processing and particle image velocimetry (PIV) techniques to fluid mechanics problems has been reviewed by Hesselink (1988), and Adrian (1991).While many new reported techniques actually involve comparatively straightforward application of simple and proven image processing algorithms, the *Current address: Jet Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, CA 91109 Correspondence to: G. R. Spedding MS 300-235, accuracy, precision and robustness in the presence of noise of these techniques is quite often unknown and overlooked.In the absence of such a complete error analysis, the reliability and utility of PIV methods must be viewed with some scepticism.Pioneering efforts towards detailing relevant error analyses were described by Imaichi and Ohmi (1983) and by Agfii and Jimenez (1987)  .The former concentrated primarily on image digitising errors, while the latter were the first to give a thoughtful consideration to interpolation errors that arise when irregularly spaced data reduced to some regular grid.
The progress of quantitative imaging techniques for velocity field measurements can be monitored in the pages of Experiments in Fluids, and three broad classes of techniques seem to have gained a foothold in the community.Methods involving the estimation of velocity from coded image streaks (e.g.Imaichi and Ohmi, op. cit.;Gharib et al. 1985Gharib et al. , 1986)), and those requiring tracking of individual particles from frame-to-frame of digitised images that are closely-spaced in time (e.g.Perkins and Hunt, 1989;Malik et al. 1993), both require at least one interpolation step before results can be represented on a regular grid, where further analysis usually takes place.By contrast, the popular cross-correlation methods (e.g.Adrian, 1988;Willert and Gharib, 1991;Utami et al. 1990Utami et al. , 1991) ) result in vectors that apparently are already given on the regular grid generated by the discrete interrogation steps (in both analog and digital implementations).This is not strictly true, however, as the best estimate of the location of the velocity vector from any particular correlation, is not at the centre of the correlation window, but at a point halfway along the computed displacement.Consequently these DPIV methods should also use an interpolation scheme, and are in fact obliged to do so to reduce certain systematic errors that will give spurious divergent flows where the vorticity is large.
It appears, then, that the issue of avoiding large interpolation errors remains pertinent in the majority of PIV applications.It is also an issue of some importance in numerical simulations involving transformation from Eulerian to Lagrangian reference frames in fluid mechanics (see Yeung and Pope, 1988).This paper is based on more extensive results in a detailed report on automated image processing and grid interpolation techniques that has been available for some time (Rignot and Spedding, 1988 [RS88]).The focus here is on the interpolation part of that report, where the presentation and discussion of the spline interpolation methods, in particular, has been considerably refined and updated, following further research into the literature and four years of practical experience using these techniques.The original report is available on request and includes a PC-readable floppy disk with source code for all the algorithms used here.

Interpolation
may be forced to reconstruct exactly the original data, and the condition of continuous second derivatives will guarantee some minimum degree of smoothness of the surface fit.
The performance of these different interpolators will be examined on real and simulated velocity fields, and an attempt will be made to describe some of the strengths and weaknesses of each technique, with a view towards offering practical guidelines concerning the choice of method for experiments in fluid mechanics.

Inverse-distance weighting methods -(AGW)
One of the simplest interpolation schemes computes the value of the interpolation function F (x, y) at a point from the weighted sum of each of the k known data points, fk, Lancaster and Salkauskas (1986) provide an excellent introduction to techniques in curve and surface fitting, and describe a broad array of polynomial, spline and finite element techniques.Comparisons and brief descriptions of a large number of interpolating functions on standard test data sets have been reported by Franke (1979Franke ( , 1982)), although the results were described in qualitative terms and the dependence on various of the parameter ranges was not clearly identified, nor would this have been feasible on such a range of techniques.In more specific reference to fluid mechanics applications, Imaichi and Ohmi (1983) have clearly described the implementation of a simple linear least squares technique and its application to two dimensional vortex flows, and a variation on this scheme was successfully applied by Spedding and Maxworthy (1986) to 2D unsteady separated flows.In two most interesting papers concerning polynomial interpolation and particle tracking techniques, Jimenez and Agfii (1986, [AJ87]) reported small advantages in certain lower order polynomial interpolators and in kriging techniques for samples in the presence of noise, but opted for a simple Gaussian convolution operator in their particle tracking experiments, on grounds of reduced computational effort and complexity.Nguyen Duc and Sommeria (1987) described the use of a 2D spline interpolation on particle streak data from vortex couple experiments, and this technique, in turn, was based on the work of Paihua Montes (1978).

Introduction to interpolation techniques
These last two interpolation techniques serve as examples of two different classes of interpolator.Typically, those involving convolution with operators, or windows, of different kinds, will be simple to implement and fast to compute.There will be some kind of inherent smoothing, whose exact nature and the particular side effects will vary according to type.Spline interpolators, on the other hand,

Wk = dE 1
as originally discussed by Shepard (1968).Another is where a is the width of the Gaussian function.
After using ( 1) and (2) to compute the two components of velocity (u, v) in the plane (x, y), the vertical component Ov Ou of vorticity, ~Oz = 0x 0y can be computed using finite differences: h is a relative displacement, which may or may not equal the original grid spacing s.The main advantages of this technique are its simplicity and its speed.

Global Basis Functions -(STS)
The basic idea is very simple in essence.One simply seeks a function Gk(X, y) and determines a set of coefficients, 2k such that F (x, y) = ~k 2kGk(X, y), (4) interpolates the data.The trick, of course, lies in the appropriate choice of Gk(X, y).A number of possibilities have been considered, as described quite clearly by Franke (1982).In the multiquadric (MQ) method of Hardy (1971) (also Hardy and Gopfert, 1975), for example, where r is a parameter chosen by the user.Alternatively, a class of spline interpolating functions known as thin-plate splines, or thin-shell splines (with rotational symmetry in 91z), discussed in rigorous detail by Duchon (1977) and Meinguet (1979a, b), requires no such free parameter, and can be shown to minimise an energy function, roughly analogous to the bending energy of a thin plate, or in this context, the Laplacian of the fluid velocity vector field.Specifically, one searches for functions that minimise the functional, (5) For O=9t 2, Duchon showed that there exists a single continuous function, a, that minimises (5), written in the form: o'(t)= ~ ~,kGk(tk, t)-I-~lXq-~2yq-fl, One then obtains the coefficients 21,--., 2,,el,0~2, fl by solving the following system: where K is the reproducing kernel matrix, K=[kJ .... kij=Gk(tk, t) The bulk of the part of the thesis of Paihua Montes (1978) dealing with 2D splines concerns the design, evaluation and performance of numerical techniques for solving (8), based on careful investigation of the properties of K. Three alternative techniques were proposed.The first (method A) involved a change of basis function based on the data themselves, assuming the first three points in the data to be non-aligned.This is the method that will be used in the following comparisons, and is the one that is discussed in further detail, together with FORTRAN source code, in [RS88].Montes' methods, B and C involve direct solution of the system using a Cholesky decomposition of a re-ordered form of (8), and an iterative, gradient descent technique.The original thesis should be consulted for details.In all cases (including the software distributed in [RS88]), we actually use a smoothed spline interpolation, where the functional to be minimised has the form: k=l c.f. ( 5).This expression represents a balance between the global smoothness and the sum of the squares of the residuals, weighted by p.The system (8) now becomes, E A F where For p > 0, ~ behaves like a smoothing parameter.
The matrices K (or /~) and E are the same for both components of velocity (substitute u or v for f in the equations above), and are quite independent of the particular choice of coordinate system (x, y).A further advantage of the spline solution is that spatial derivatives of fin (x, y) are given analytically from differentiation of (6), so, for example, e)=, can be computed directly, since for the ith grid node, one immediately has k=l where the ,Ix, 2 coefficients are computed from the UR and vk velocity components respectively.

Smoothing
The rational choice of smoothing parameter is an entire topic unto itself and it will not be addressed explicitly here.A number of internal tests, described more fully in [RS88], showed that a reasonable choice of smoothing for the AGW technique was to smooth all grid data (velocities and their spatial derivatives) with an 8-neighbour mean smoothing filter.For both the test and experimental data, a fixed value of ~ = 1.0 was used in the STS interpolations.It is not claimed that this is the best choice (see, for example Wahba (1979Wahba ( , 1982Wahba ( , 1990 pp 45-65) pp 45-65) for possible approaches in computing optimum smoothing in thinplate splines), only that it allowed reasonable solutions to the stystem (10) to be computed that were not greatly in error, and were not oversmoothed either so that the effect of changing the other interpolation parameters could be clearly observed.Subsequently, in practical use, with real, noisy data, an appropriate and careful choice of ~-=~([v[ .... Inoisel, L/6), has been found to be essential.

Error functions
Consider, to begin with, the formulation of the Adaptive Gaussian Window (AGW) interpolation in Eqs.(1)(2)(3).The performance of this scheme depends, a priori, on four parameters: L: characteristic length of the flow 6: mean nearest neighbour distance between samples.h: displacement in the numerical differentiation.a: Gaussian window width.
If a global measure of the interpolation error is taken, such as the root mean square of the difference between the true data field and its interpolated values at the grid nodes, this quantity will be a function of these four parameters.Hence, RMSE=F (L,6,h,a).
Although, in principle, h may be chosen independent of the real grid spacing, s, for convenience, we set h = s for all the results reported here, so the grid spacing and the differentiation length are equal.A further simplification can be introduced by using the result of [A J87], who concluded that an optimum ratio of HI6 was around 1.24, independent of h/6 and L/6.We have retained this magic number of H/6, so the RMSE is now a function of three variables ~.This may be expressed in two intuitive forms: and RMSE = g(L/6, h/6).

(13b)
These formulations have quite practical interpretations, which have some influence on how one might x6 in [A J87] is defined as the mean particle separation, not the mean nearest neighbour distance, as used here, which diverges from the former as the particle distribution becomes less uniform.view the following performance analysis measurements.Assume that L is known or can be reasonably approximated.There are therefore just two alternative scenarios to be considered.In the first case, suppose the data must be interpolated on to a grid with fixed, predetermined dimensions, h and L/h are fixed, and the relation (13a) will give minimum requirements for 6 (i.e. the minimum number of data points to be obtained experimentally) to satisfy a given degree of accuracy of reconstruction of the field.Alternatively, in the second case, an image (or any other source of measurement) has already given a certain number of particles or data points, 6 is thus fixed, and an appropriate choice for the grid spacing, h will be given by (13b).
Since the window width in AGW, a, has been eliminated from the functions f and g, they may equally be applied to the Thin Shell Spline (STS) interpolation, which is a function of the same three parameters L, 6 and h.

Error analysis procedure -Burgers' vortex
The interpolation errors committed by AGW and STS were measured in tests on artificial velocity fields, and also estimated in real data.For economy of presentation, the errors are described by single RMS values, computed at the nodes of the grid and averaged over the whole field.
The artificial velocity was prescribed by Burgers' vortex, which is a closed-form solution of the Navier-Stokes equations, providing a simple canonical vorticity distribution.In a 2D flow, the tangential velocity, vt is: where x is the circulation, r the radial distance from the core centre and v is the core size.The vertical vorticity, COz, is then Varying numbers of velocity vectors, n, were distributed about this velocity field, at locations determined by two independent random number generators.The distribution of vorticity og~(x, y) and an example of randomly-located velocity samples in this field are shown in Figs.la and lb, respectively.
In real data, the true error is always unknown, and is usually estimated from the variance of resampled data.In unsteady and/or turbulent flows, this is rarely practicable; instead, one may estimate the error by resampling the same data, synthesising new samples with random number generators.This technique is known as bootstrapping (Efron, 1982), and is discussed more thoroughly by [A J87] and in I-RS88].Here the bootstrap techniques were implemented so as to calculate errors at the grid nodes (rather than at the data points) for both AGW and STS interpolators.Some general conclusions and operating principles for the AGW interpolation may be drawn: 9 A single vortex field may be reconstructed with an average error of around 5% for the velocity, and 10% for the vorticity, given appropriate choices of h/6 = 1 and L/h = 20 for this case.It is difficult to see how this error can be reduced below these values in most foreseeable practical applications.9 In real applications, given L, or an estimate for its likely value, and constraints on either h or 3, one may exercise one's remaining degree of freedom with an appropriate choice according to the curves in Fig. 2.

AGW Interpolation Errors
Figure 2a shows the mean RMSE function (9 in Eq. ( 13b)) of L/6, for different values of hi6, for the velocity field of the pure vortex of Fig. la, computed from a number of iterations of AGW on randomly distributed velocity samples.For this flow there are originally no small scales to be detected, and the relative grid spacing hi6 has little effect on the results.On the other hand, plotting the same function, but for the vorticity field (Fig. 2b), shows a strong influence of hi6, which should be greater than 1 and probably less than 2 (for a typical range of L/6), to keep errors down to acceptable levels.Expressing the RMSE velocity as a function of hi6, for different L/h (fin Eq. ( 13a)) in Fig. 2c, the rapid increase in the RMSE for hi6 < 1 is clear, as is the expected trend for the RMSE to decrease, for fixed hi6, with increasing L/h.The same effect is evident in Fig. 2d for the vorticity.

STS Interpolation Errors
A similar performance analysis for STS is displayed in Fig. 3. Figures 3a and b show the RMSE function (g in Eq. ( 13b)) of L/6, for different h/6, for the velocity and vorticity fields, respectively.In both cases, g is quite insensitive to changes in hi6, while the importance of achieving a satisfactory L/6 is obvious.STS does not have the same dependence on h/6 for the vorticity as was observed for AGW, owing to the different techniques for computing derivatives of u.Expressing the RMSE as given by f in Eq. (13a), Figs.3c and d show this function for the velocity and vorticity fields, varying L/h and h/6, as with AGW in Fig. 2. The results are similar too, with the relative effects of L/h and h/6 being quite clear.In general, it is these curves which will be the most useful for making appropriate choices of parameters and interpolators in practical applications.To summarize: 9 Interpolation errors of approximately 2.5% for the velocity, and 5% for the vorticity might be expected, given suitable choices of h/d and L/h. 9 STS is not especially sensitive to the particular choice of h, and so for a fixed or known L, one is simply left with a choice of 6 in order to interpolate within a certain error.

AGW vs. STS
In terms of the interpolation error performance on artificial velocity fields composed of pure vortices, inspection of Figs. 2 and 3 reveals some differences in the two techniques.STS is much less sensitive to h than AGW and, in fact, for this type of data, is more or less independent of the grid spacing, within reasonable limits.On the other hand, h has a significant influence in AGW as it governs the differentiation displacement length.It is true that we have chosen to couple h and s (h = s) in this analysis, and similarly the fixed ratio of a= 1.246 determines the Gaussian window width.There is some degree of flexibility in AGW, then, which has not been fully utilised, but the general point remains that this scheme is more sensitive to aspects of the grid spacing and geometry.Operating under the same conditions, the RMS errors in both velocities and their first spatial derivatives were roughly halved by STS in comparison with AGW, and on these grounds alone, STS is clearly superior.Sometimes small differences in RMSE values belie rather more obvious and dramatic differences in the distribution of errors across the whole data field.To illustrate this point, isometric surfaces of e~(x, y) interpolated from data of Fig. lb by STS and AGW appear in Figs.4a and  b respectively.
Both interpolations in Fig. 4 give acceptably accurate values for COm,, (the true value of m~,ax = 1) and for F, the circulation (measured by integrating over the surfaces of Fig. 4), but from a qualitative standpoint one is much more likely to be satisfied with the STS result.There are many spurious oscillations in the AGW co field, including false peaks.On the other hand, AGW does have the property that o) drops to zero where there is no original data.Compare Fig. 4b with Fig. lb at the edges of the grid to see that it will do this regardless of the remainder of the data field; in short, where the data are sparse, there is no implicit extrapolation or interpolation from surrounding areas aside from the Gaussian itself.This point is really rather obvious, but it is all that is required to explain the perhaps initially surprising unsmoothness of Fig. 4b.By contrast, STS displays the classical spline curve characteristics where the extrapolated end points are dictated by smoothness conditions of the interior points.If one knows these end conditions a priori, they may be specified, but in this case they were not, and the errors are the corners of the STS interpolation are the result.

Analysis of Bootstrap Error
The ability to accurately estimate unknown errors in interpolated velocity fields would be a valuable asset, and the bootstrapping algorithms were thoroughly tested on artificial velocity fields to check the agreement with the true errors.The solid curves of Figs 9 5a and b  Fig. 3a-d.STS performance analysis.RMS errors in u and r as functions of h/6, L/6 and L/h, as for AGW, above values of L/h.The dashed lines are the bootstrap error estimates under the same conditions.The trend for the bootstrap error in the velocity field is always the same as the true error, with both L/h and h/6, but the absolute magnitude of the error may itself be in error by up to 100%.Fig. 5b shows the same general trends for the vorticity field.
Once again, in Fig. 5b the bootstrap error always underestimates the true value.This is because both interpolators have a nonzero smoothing, so a bootstrap error estimate computed at the grid nodes will be an underestimate due to repeated smoothed samples having a lower variance than if the smoothing were absent.This is quite important in correctly interpreting the bootstrap result.The absolute value of the bootstrap error may underestimate the true error by a factor of two 2, but the relative change in bootstrap error with L/h and hi6 agrees well with the real error curves, and Fig. 5 suggests it could be a useful measure of comparative errors across the spatial domain.
2 This factor will depend on the degree of smoothing, and on the interpolating method.
Bootstrap algorithms were successfully implemented for both interpolators, but a single example using STS will suffice to demonstrate.Figure 6a is an isometric surface of the true error distribution in a,(x, y), O:u= U --l~l i where u and us are the real and interpolated velocities, and Fig. 6b shows the bootstrap estimate of the same quantity, computed only from the measured velocities.Fin the general notation of section 2, af =fk--Fk over the k known data points.]Note how % has its greatest magnitude at the corners, as previously observed (Fig. 4a), and how the bootstrap estimate recovers this general trend quite satisfactorily.For a~,,

~o,= I~o-co~f
In Fig. 7, both the true error (7a), and its bootstrap estimate (7b) are confined mostly to the corners and edges of the grid.The difference in scales in the two figures exaggerates the apparent extra noise in the bootstrap error field, which would be further reduced with more iterations.It appears that a realistic distribution of errors a(x, y) can be generated from the data, in which case these values might be used as limiting constraints in smoothing filter operations ([AJ87], see also Wahba (1990, pp 67-71)).This was not attempted here though, owing to the considerable uncertainty in the magnitude of the bootstrap error estimates.

Example Applications
It is important to verify the conclusions from numerical simulations in real experimental data, where sources and forms of noise are usually less uniform and predictable.Two examples of unsteady vortex shedding experiments follow, one in which the flow is close to completely two-dimensional, and data come from one instantaneous realisation.The second flow is strongly threedimensional, with flow out of the plane of measurement, and measurements have been phase-averaged from multiple realisations.

Oscillating 2D Flap
A two dimensional oscillating flap experiment was conducted by hinging a flat plate at the leading edge of a horizontal flat plate spanning a water channel.The streamwise length of the flap, c, was 5 cm, the freestream velocity, U, was 11 cm.s -1, and the flap was oscillated continuously at a frequency of .33Hz.The reduced frequency, k, defined by, k=nfc/U, was 0.47.The crossstream vorticity thus dominates the local flow field, which is characterised by the shedding of strong vorticity at the trailing edge of the flap.At the maximum flap opening angle of around 70 ~ this vorticity has rolled up into a single coherent vortex, shown in Fig. 8a.The flow has been seeded with polystyrene beads, and a modulated laser light source is directed along a streamwise, vertical sheet, with a thickness of 2-3 mm.Below this photograph, Fig. 8b shows the velocity vectors for this flow, computed from the digitized image of Fig. 8a.The theoretical accuracy of these vectors is approximately 0.3 pixels (see [RS88] for details).Figure 10 compares the AGW-and STS-interpolated results.The choices of L/h= 14.5 and h/6=0.8 were determined by the performance curves of Figs. 2 and 3, given the constraint that the particle density is a little low.The 15 x 15 grids in the top of Fig. 10 are on restricted areas over the vortex (the origin of the x, y coordinate system is at the hinge point of the flap on the plate) and, qualitatively, the velocity fields look quite similar.The vorticity distribution (lower figures) is rather different in the two cases, however, with the peak e)m,x being sharper and higher in the STS interpolation.
The difference in the two interpolations is not so great when averaged over the whole field so that, while ratio of co .... AOW/e) .... STS = 0.55, the circulation, F, computed by integrating co over all (x, y), is -11.3 cmZ.s -1 according to AGW, and -13.6crn2.s-a from STS, a ratio Fa~w/FsTs=0.83.This reflects the tendency of AGW to smear out a peak over a broader area.According to bootstrap error estimates, the RMSE, is 11.1% for AGW and 6.4% for STS.These values are quite consistent with the predictions from the performance curves of Fig. 2 and 3 for this range of L/h and h/6.In this example where n (=44) is small, and the velocity vectors are known with good accuracy, STS is clearly the best choice of interpolator.The difference in CPU time requirements is negligible but the superior accuracy of STS in determining o)(x, y) is quite noticeable.

3D Separated Flow over a Delta Wing with Flap
The delta wing experiment with a leading-edge flap is described more fully in Spedding et al (1987), and in [RS88], and is shown schematically in Fig. 12.The delta wing has a flap of local span 0.45b, where b is the semispan, and the flap oscillates sinusoidally about the leading edge, through an angle of 70 ~ A light slice was placed normal to the freestream, at the wing trailing edge, and particle trajectories in the (y, z) plane were digitised by hand and phase-averaged over several frames.Due to the manual digitising technique and the phase-averaging of the data, one might expect there to be more errors in the original velocity vectors, which are distributed very unevenly over (y, z).The strong axial flow along the vortex core meant that few particles remained there in the light slice during the time of exposure, and the sparse data in the far field is of somewhat dubious quality.In practice, some of these problems can be avoided by filtering and zero, unless some constant, uniform offset is added), and so the STS interpolation cannot be expected, in general, to correctly reproduce the edges of flow fields where this occurs.Either some zero velocity data points are inserted at the appropriate locations in the flow field, or an incorrect uniform flow must be tolerated.AGW grid values decay to zero by default, in the absence of data, but at the expense of generating artificial and probably erroneous vorticity data at those grid points.x(cm) 12

Bootstrap error estimates on real data
Detailed results are available in [RS88], and only one example will be shown here.Figures 13 and 14, show, on the upper surfaces, the original interpolated data of u 2 and o9 for the STS interpolation, and on the lower surface are the bootstrap distributions of (tru) 2 and tro,.The units of each pair of surfaces are therefore the same in each case, except the bootstrap error scales have been expanded by a factor of 100 and 10 for (au) 2 and ~r,0 respectively.The curl-up of the error surface (au) 2 in Fig. 13 at the edges close to the vortex is reminiscent of the artificiallyreconstructed Burgers vortex result.Overall, the small magnitude of these errors is encouraging, and the vorticity surfaces (Fig. 14) present a similar picture.The shape of the vorticity distribution in the top of Fig. 14 is quite likely to be correct, within error bounds given by the lower surface.The RMS value of the relative error, ~,, over all the m x n grid nodes, When the particle distribution is so non-uniform, the optimum grid spacing for AGW will vary across the field, and globally-averaged values of 6 may not yield the best distribution of mesh points.This matter is discussed further in [RS88].A reasonable solution was arrived at with a value of h/g= 1.27, slightly higher than the optimum value of 1 derived from tests with uniform random particle distributions.The result is shown in Fig. 13.The circulation for this separation vortex, F =--34.1 cm 2. s-1.
The choice of h is not clearly so crucial for STS, and h/6 = 1 can be safely used, resulting in a denser distribution of grid points (and possibly superior spatial resolution).The velocity and vorticity fields are plotted in Fig. 14.In contrast to AGW, the velocities do not decay to zero at the edges of the grid, but instead tend to remain close to interior neighbouring values.Strictly speaking, most particle streak detection schemes are incapable of resolving a zero velocity vector (when the displacement is is 2.5%.Noting that this could be an underestimate by a factor of two, the upper bound on the true interpolation error is thus unlikely to exceed 5%; this number is consistent with the quantitative predictions of Fig. 3 for this interpolation scheme.

Conclusions
Given the constraints imposed, STS can be expected to interpolate scattered data with errors about one half of the magnitude incurred by AGW.Moreover, its comparative insensitivity to small changes in the grid spacing, h, and the mean nearest neighbour distance, 6, make STS simpler to use.In both cases, it may be possible to significantly improve the performance by careful adjustment of local window width, tr (in AGW), and smoothing parameter, (in STS), in response to noise and/or spatial nonuniformity in the data distribution.
The principal justifications for using a scheme of the AGW type for interpolation are its simplicity and the convenience of implicit smoothing.However, in almost all circumstances tested, given an appropriate choice of free parameters, STS can significantly outperform AGW in terms of mean and local error in velocity and vorticity.The report referred to as [RS88] here provides source code for the spline interpolation algorithm and so complexity of implementation should no longer be a hindrance.Furthermore, now that CPU and RAM limitations are less severe, there seems to be little reason to use AGW-type routines for all but rough solutions to problems involving very large n.In these cases of large n (n > 2048 can be considered large by current workstation limitations) the problem is essentially the requirement in Eq. ( 9) of inverting an O(n2), double precision matrix.This in turn stems from insisting on a global solution to the spline surface coefficients, and can be removed if some local, patched spline interpolation is used instead.
In practical terms, the procedure is to subdivide the data into smaller domains, where local subrectangles are spline interpolated, and glued together at common edges by polynomials that can be computed using information in the neighbouring spline coefficients themselves.Now, the choice of numerical solution should reflect the fact that it is more difficult to track and check on the alignment of the first three data points in each subrectangle, and so the solution of (10) proceeds by a direct solution of a re-ordered projection of (10) on a basis computed from all the data.This is Montes' method B, which is a little more time consuming than A, and slightly less accurate, although the differences in accuracy, following a simple , y), tk=(Xk, Yk), and Gk(tk, t)= I&l 2 log[dk[ 2 .(7)

Fig. 1 .
Fig. 1. a The analytic function ~(x, y) for Burgers' vortex, b A randomly-chosen sample (N = 80) of velocity data points in the flow of a.
Fig.2a-d, AGW performance analysis.RMS error functionsfand g, in u and m as function of h/6, L& and L/h.Each point is an average of 4 or 5 tests, a g.; b go; el,; df~, are reproduced from Figs 9 2c and d and show the RMSE function of velocity and vorticity for AGW with h/5 for four different

Fig
Fig. 4a, b.Interpolations of the velocity vectors in Fig. lb, onto a 20 x 20 grid by STS a, and AGW b Fig.5a, b.The RMSE functions of u, a and co, b, for AGW (solid lines) and from the AGW bootstrap estimate (dashed lines).The RMSE values are shown as functions of L/h and h/6 and the solid curves are those already shown for AGW in Figs.2c and d Fig. 6a, b. a.(x,y), a and its bootsrap estimate; b for the STS interpolation of Fig. lb . 7a, b. a,o(x, y), a and its bootstrap estimate, b for STS

Fig. 8
Fig. 8. a Original particle streak image of the 2D flap, at an open-, tro~ ing angle, 0 = 76 ~ close to the maximum.The mean flow is from left o9 to right b Velocity vectors assigned to the flow above

Fig. 9 .
Fig. 9. Grid velocity vectors and vorticity distributions from the AGW (left), and STS (right) interpolations of the 2D flap data from Fig. 8b

Fig. 10 .
Fig. 10.Cross sections of the 3D flow over a delta wing/flap configuration are imaged in an inclined mirror, downstream of the trailing edge.The light slice is normal to the flow and deliberately thickened by defocussing the laser