On the applicability of numerical image mapping for PIV image analysis near curved interfaces

Numerical meshes in computational fluid dynamics (CFD) are used to discretize the non-linear partial differential equations (PDEs) modelling the flow in order to allow the computation of an approximate solution. These meshes can be categorised as either structured or unstructured. While the former is defined by the intersection of curvilinear coordinate surfaces, the latter lacks any relation between coordinate directions. The current intention is for the mapping to describe the spatial rearrangement of image intensities when altering reference coordinate systems. As digital images are composed of intensities defined in pixel elements which are already arranged in an array format, structured meshes are consequently the most suitable grid type.

The generation of a structured mesh can be performed through several approaches of which the most common constitute analytical transformations, algebraic schemes and other based on partial derivative equations (PDE). Details regarding numerical grid generation can be found in Thompson et al (1999) but it presently suffices to note that while analytical transformations are the fastest and most accurate, their use is limited to more simplistic analytical geometries such as cylinders, sinusoidal walls, etc. In view of the intended general applicability of image mapping these will be left out of consideration. Grid generation based on partial differential equations is typically slower than its algebraic counterpart but is the preferred choice in the present work as they offer a higher degree of mesh smoothness and can concomitantly be expected to minimise image artefacts caused by image mapping. In fact, algebraic meshes are based on transfinite interpolation and are more likely to result in overlapping of the boundaries (Thompson et al 1999)

The numerical mesh is a discrete set of curvilinear coordinates (ξ, η) satisfying either an elliptic or hyperbolic system of PDE equations. Independent of their classifications, these equations describe the coordinate transformations when alternating between the physical (x, y) and logical (ξ, η) plane. Elliptic meshes are obtained by solving the system of elliptic equations presented in equation (1) through iterative schemes based on finite differences. The equation of this mesh comes from an analogy to a stretching membrane attached to the boundary of the object. These meshes are perfectly suitable for closed domains, where all the boundary geometries are specified and are considered best in terms of smoothness (Thompson 1987):

$$\begin{eqnarray}&&\frac{{{\partial}^{2}}\xi}{\partial {{x}^{2}}}+\frac{{{\partial}^{2}}\xi}{\partial {{y}^{2}}}=0\,\text{and}\,\frac{{{\partial}^{2}}\eta}{\partial {{x}^{2}}}+\frac{{{\partial}^{2}}\eta}{\partial {{y}^{2}}}=0.\end{eqnarray} \tag{ 1 }$$

Starting from the object geometry identified in the image recordings, boundary conditions are provided in terms of x and y. However, the solution to equation (1) expresses the unknown angular, ξ, and radial coordinates, η, as functions of the physical coordinates, x and y. System (1) is therefore usually rearranged such that the physical coordinates become the unknowns and the boundary condition can be defined with respect to the body geometry when η  =  0, as illustrated in figure 3. An additional user-defined non-homogenous term can be added to the equations to promote special grid properties such as line spacing, orthogonality, etc.

Alternatively a hyperbolic system of partial differential equations forms the basis for the intended mesh generation. Here the equations are analogous to the description of waves transporting disturbances on water surfaces:

$$\begin{eqnarray}&&\frac{\partial x}{\partial \xi}\frac{\partial x}{\partial \eta}+\frac{\partial y}{\partial \xi}\frac{\partial y}{\partial \eta}=0\,\text{and}\,\frac{\partial x}{\partial \xi}\frac{\partial y}{\partial \eta}-\frac{\partial y}{\partial \xi}\frac{\partial x}{\partial \eta}= \Delta A.\end{eqnarray} \tag{ 2 }$$

The first equality in (2) imposes orthogonality while the second controls the local cell area through ΔA. Parameter ΔA is user-defined and can subsequently take on any form; ΔA  =  constant, ΔA  =  ΔA(x, y) or ΔA  =  ΔA(ξ, η). The solution methodology starts from the geometry of the boundary conditions and propagates outwards until a user-defined extent. For this reason, this type of mesh is more suitable for open domains.

It is important to note that when attempting to generate a body-conforming mesh to map PIV images, the purpose of the grid, and therefore its inherent characteristics, will be slightly different from those used in CFD. Computational fluid dynamics utilises the mesh to obtain a discrete solution of the underlying flow equations at every grid node, therefore, an accrual of mesh points in vicinity of the surface is beneficial to the estimate of velocity near the object. In the current work, instead, the mesh is used to deform images that are thereafter analysed in sub-areas to evaluate the displacement of particle images. The choice of optimal mesh characteristics should thus introduce the least amount of artificial displacement gradients in the logical plane as these are well-known to hinder reliable cross-correlation (Westerweel 2008). For this reason, the optimal mesh should maintain a constant node density independent of the distance from the boundary or its curvature radius. When comparing elliptic and hyperbolic meshes, the former tends to automatically densify the number of points where surfaces have a lower radius of curvature and has the opposite behaviour where surfaces are smoother. This effect is to be avoided for image mapping purposes because of the aforementioned artificial gradients. Moreover, the iterative approach required to solve an elliptic system of partial differential equations renders the associated mesh generation slower compared to the hyperbolic type, where a faster marching technique can be adopted (Thompson et al 1999). These motivations lead the authors' grid choice to a hyperbolic mesh.

Having identified hyperbolic meshes as most suitable for PIV image analyses, this section highlights aspects in the mesh generation which require careful consideration. Utilising the hyperbolic partial differential equations expressed in (2) and a complementary boundary condition describing the pixel locations (x, y) of the interface geometry in terms of logical coordinates (ξ, η), the solution of the well-posed system of equations provides for each combination of curvilinear coordinates (ξ, η) the corresponding physical coordinates (x, y). This constitutes the image mapping. In this work, the authors opted for a simple first-order explicit solver because of its low computational effort. Details on the solver are provided in the appendix. Although the problem expressed in (2) is closed, three more parameters must be specified to generate the desired mesh; the local cell area ΔA, the wall-normal extent of the mesh and a smoothening factor. Maintaining a uniform cell area may be the most simplistic approach, yet it leads to severe stretching of the cell aspect ratio when increasing the radial extent of the mesh in regions of stronger boundary curvature as illustrated in figure 4(a). In section 5.2 the authors will show that such stretching introduces detrimental artificial displacement gradients in the mapping and should be categorically avoided. The area of the cells must therefore increase according to the spreading of the radial coordinate η. Here, an automatic adjustment of ΔA to the transformation's tangential derivatives is proposed:

$$\begin{eqnarray} &&\Delta A= \Delta A\left(\xi,\eta \right)= \Delta {{A}_{0}}\sqrt{{{\left(\frac{\partial x}{\partial \xi}\right)}^{2}}+{{\left(\frac{\partial y}{\partial \xi}\right)}^{2}}}\end{eqnarray} \tag{ 3 }$$

with 0 designating the initial condition. The effect of altering the cell-area according to the local Jacobian is exemplified in figure 4(b). The reader should bear in mind that the illustrated mesh in figure 4 is representative of the subsequent image deformation and is not related to the sizing or positioning of interrogation windows.

Zoom In Zoom Out Reset image size

The extent of the numerical mesh is defined by the user and should equal at least one initial correlation window size. The degree of grid smoothening must be selected considering the shape of the interface and is of pivotal importance both for numerical stability and mesh coherence. As an example, in absence of artificial smoothening, concave shaped interfaces are attributed meshes which intersect themselves as they propagate outwards. This creates a non-unique surjective transformation of the mesh as illustrated in figure 5(a). Incorporation of artificial smoothening by means of adding an explicit second-order dissipation term to the right hand side of equation (2) as proposed by Chan and Steger (1992) renders a more conducive mapping (figure 5(b)).

Zoom In Zoom Out Reset image size

When images will be transformed into the logical (ξ, η) plane and evaluated, corresponding particle image displacements dξ and dη must be converted into physical displacements dx and dy in order to be meaningful. The transformation of logical displacements is achieved using the Jacobian [J] of the transformation, expressed by equation (4):

$$\begin{eqnarray}&&\left[\begin{array}{c} \text{d}x \\ \text{d}y \end{array}\right]=\left[\begin{array}{c} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{array}\right]\left[\begin{array}{c} \text{d}\xi \\ \text{d}\eta \end{array}\right]=[J]\left[\begin{array}{c} \text{d}\xi \\ \text{d}\eta \end{array}\right].\end{eqnarray} \tag{ 4 }$$

Derivatives appearing in (4) are evaluated through an interpolation of the second-order finite differences of the physical coordinates (x, y) with respect to ξ and η.

The generation of the numerical mesh starts from the boundary condition provided in the form of a logical binary image identifying the object. As a direct consequence of pixelisation, edges cannot be straightforwardly used as mesh boundary conditions but need to be properly smoothened and rearranged (figure 6). The process of boundary smoothening can be performed in different ways and several techniques have been tested among which Gaussian convolution, discrete cosine transform (Garcia 2010) and Fourier descriptors (Feudjio et al 2014). In spite of the number of the techniques tested, all smoothening functions applied to the mask edges produce artificial oscillations due to the non-random nature of the quantization error of the discrete pixel positions making up the interface boundary. In order to overcome this problem, the authors propose a line simplification algorithm rather than smoothening for which the Douglas–Peucker algorithm (Douglas and Peucker 1973) has been efficiently adopted.

Zoom In Zoom Out Reset image size

The process of generating a boundary condition from the provided masking image for subsequent mesh generation is depicted in figure 6(b) and consists of four main steps. Edges are first extracted from the logical mask by means of classical edge detection routines. Since the edge of the mask is made of discrete pixels the suggested algorithm simplifies the step-shaped curve by iteratively reducing the number of points and selecting subsets with the aim of attaining a difference between the original and sub-sampled edge below a defined threshold. Following a spline interpolation of the simplified edges to repristinate the number of edge pixels, edge points are redistributed along the spline curve to be equidistant. The latter is of importance since a non-equally distributed mesh introduces artificial deformation into the mapped image, affecting the quality of the cross-correlation and therefore the displacement measurement accuracy. The number of points used to re-interpolate the simplified edge will affect the number of points of the mesh in the angular direction and consequently the accuracy of the transformation. For this reason the mesh density should be selected to ensure each image pixel is covered by at least one mesh node.

Another important aspect concerning the generation of the boundary conditions for a numeric mesh appertains the overall topology. The main obstacle is that the logical plane is always rectangular shaped. Although the transformation from a general shape to a rectangle can be quite straightforward (e.g. L- or C-shaped boundaries need only straightening) closed shapes need further attention. An annulus, for instance, cannot be transformed into a rectangular shape unless a branch cut is introduced on the physical plane. No general solution exists to overcome the inherent problem and the user is obliged to recursively adjust the mesh on a case-by-case basis to ensure continuity across the branch cut when interrogation windows overlap with it. Particular topological conditions could be needed also in case of shape singularities or strong curvature.

The hyperbolic numerical mesh is generated through the system of equations (2), which can be solved adopting e.g. the approach recorded in the appendix. Independent of the adopted solver, the obtained numerical solution consists of a discrete set of data points in the physical plane (x_ij, y_ij), which are located on the nodes of the numerical mesh and for which the corresponding equivalent in the logical plane (ξ_ij, η_ij) is known. Figures 7(a) and (b) show an example of those discrete points in the form of a 3D surface for a simple sinusoidal mesh. It should be noted that the locations (ξ_ij, η_ij) constitute the nodes of an equispaced rectangular grid and that the samples (x_ij, y_ij) do not need to coincide with original pixel locations. As will be shown in section 4.3, cross-correlation in the logical plane however requires image intensities defined at a higher resolution, i.e. at locations between the available (ξ_ij, η_ij) samples. The transformation between a general point on the logical (ξ, η) plane and the physical (x, y) plane thus requires an interpolation of the discrete point locations. In the current work the authors have adopted a quintic spline scheme (Astarita and Cardone 2004).

Zoom In Zoom Out Reset image size

The inverse transformation from the physical plane to the logical plane, which is depicted in figures 7(c) and (d), is also performed through interpolation, although additional attention is required since the coordinates (ξ, η) are distributed on the nodes of the mesh on the physical plane and, as such, do not constitute a structured grid as required for spline interpolation. For this reason, a non-structured interpolant is needed and in this work, a natural neighbours interpolation (Bobach et al 2009) has been selected accordingly.

Following the generation of the numerical mesh, interrogation window centres are distributed on the physical plane and, in line with the Chimera approach, analysed differently according to their distance from the interface boundary (see section 2). Interrogation areas of which the centroid is located outside the numerical mesh are analysed directly by means of cross-correlation on the physical plane as per standard PIV. Correlation windows whose centroid falls within the boundary fitted mesh will utilise the logical plane, where rectangular correlation windows are constructed and analysed by means of standard cross-correlation.

The relations established in section 4.2 describing image mapping between the different coordinate systems allow a complete reconstruction of the image recordings onto the logical plane. This is the common approach adopted in previous studies (Nguyen et al 2010, Jeon and Sung 2011, Park et al 2015), followed by cross-correlating intensities in the logical plane. However, digital images are not continuous but composed of discrete, square elements, i.e. physical pixels. Image mapping will consequently distort these pixel units. Due to the potentially strong non-linearity of the mesh shape and the consequent divergence (or convergence) of the radial lines in proximity of a curved surface, the use of a uniquely predefined correlation window size in the logical plane expressed in terms of physical pixels produces under-sampled (or over-sampled) particles at further distances from the boundary. The mapped particle intensity distributions consequently become strongly elliptical in the vertical (or horizontal) direction, biasing the outcome of the correlation operator as illustrated in figures 9(b), (c) and (e)².

This effect in itself is detrimental to the accuracy attainable with image mapping and any subsequent error analysis would as a result disfavour image mapping. Instead, to enable an objective evaluation of the overall applicability of image mapping, the authors present a fix to reduce this unfavourable effect. The authors propose a more robust approach involving image re-sampling of logical pixels prior to cross-correlating whereby the image mapping is performed on the basis of individual interrogation windows using an independent number of sampling points. As the process involves a transformation of the various units between the physical and logical plane, table 1 summarizes main related definitions to assist the reader, which are graphically clarified in figure 8.

Zoom In Zoom Out Reset image size

Strongly curved meshes distort the shape of the particle images as illustrated in figures 9(a) and (b). The stretched particle images (figure 9(c)) create a distorted correlation map (figure 9(e)) and negatively affect the mapping-related bias error as argued by Ronneberger et al (1998). To this extent, a re-sampling process in introduced. Starting from a logical rectangular interrogation window, whose sizing will be covered in section 4.4, the radial and angular dimensions of its mapped equivalent in the physical plane are evaluated in terms of physical pixel units. This yields a length γ_k (in physical pixels) for each curvilinear side k of the physical window. To maintain the circular shape of the particle images in the logical plane, the number of mapped pixels N_ξ and N_η along the tangential and radial directions respectively within each logical correlation window is then selected automatically considering these extents:

$$\begin{eqnarray}&&{{N}_{\xi}}=\frac{1}{2}\left({{\gamma}_{\text{top}}}+{{\gamma}_{\text{bottom}}}\right)\,\text{and}\,{{N}_{\eta}}=\frac{1}{2}\left({{\gamma}_{\text{left}}}+{{\gamma}_{\text{right}}}\right).\end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

The enhanced resampling thus re-evaluates the number of logical pixels dividing the logical, rectangular shaped correlation window ensuring a sufficient number of pixels to produce equally stretched particles in the horizontal and vertical direction (figure 9(d)). The shape of the particles in the mapped windows is consequently preserved yielding symmetric correlation peaks (figure 9(f)). The improvement can also be shown by evaluating a Gaussian fit of the correlation peaks, shown in figures 9(e) and (f), and evaluating the ratio between the horizontal and vertical widths. This ratio equals 5.2 for a constant sampling and 1.2 for the re-sampled correlation window, illustrating the improved correlation symmetry. Image intensities within the logical pixels are retrieved by quintic spline interpolation (Astarita and Cardone 2004) of the original intensity distribution in the corresponding physical locations obtained by forward mapping (Wolberg 1990). The reader should note that these logical pixels, which are arranged in a Cartesian manner in the logical plane (red circles in figure 8(a)) can consist of several logical units (figure 8(a)).

The outcome of the subsequent cross-correlation operation on the logical interrogation areas yields a displacement quantified in terms of logical pixels. These in turn depend on the number of imposed (re-) sampling points expressed in (5), necessitating a rescaling prior to the Jacobian-based conversion of displacements from the logical plane to the physical plane (see equation (4)). The conversion from logical pixels into logical units, for each displacement, is given by:

$$\begin{eqnarray}&&{{u}_{\xi,\eta}}=u_{\xi,\eta}^{\ast}\frac{\text{W}{{\text{S}}_{\xi,\eta}}}{{{N}_{\xi,\eta}}}.\end{eqnarray} \tag{ 6 }$$

Subscripts denote consideration of either the tangential ξ or radial direction η, u_ξ,η represents the displacement converted into logical units, $u_{\xi,\eta}^{\ast}$ symbolises the displacement in logical pixels resulting from the correlation in the logical plane and WS_ξ,η indicates the extent of the window in the logical plane in terms of logical units. The determination of WS_ξ,η will be discussed in section 4.4.

Contrary to standard PIV algorithms, interrogation window sizing during the process of image mapping is more complicated due to the non-linear correspondence between physical and logical interrogation window sizes:

$$\begin{eqnarray}&&\left[\begin{array}{*{35}{l}} \text{W}{{\text{S}}_{\xi}} \\ \text{W}{{\text{S}}_{\eta}} \end{array}\right]={{\left[J(x,y)\right]}^{-1}}\left[\begin{array}{*{35}{l}} \text{W}{{\text{S}}_{x}} \\ \text{W}{{\text{S}}_{y}} \end{array}\right]\end{eqnarray} \tag{ 7 }$$

where WS_x and WS_y are the horizontal and vertical window sizes on the physical plane. Variables WS_ξ and WS_η depict the angular and radial sizes respectively on the logical plane. This dependency was already highlighted in section 4.3 in the context of particle image distortions. As per equation (7), due to the spatial dependency of the Jacobian, equally sized correlation windows in the physical plane will yield non-uniform window sizes in the logical plane and vice versa. This problem has never been clearly stated in existing literature and the authors accordingly propose two different solutions to automatically adjust the size and the position of the correlation windows on the logical plane.

The first implementation, referred to as (straightforward) image mapping (IM), involves a recursive re-evaluation of interrogation area size and position on the logical plane in order to maintain the significance of the imposed window size. Once interrogation windows are distributed on the physical plane as per standard PIV (figure 10(a)) and selected spatial locations are transformed from physical to logical coordinates using the inverse mapping functions ξ(x, y) and η(x, y) (section 4.2, figure 10(b)), window sizes are converted from physical pixels into logical units using equation (7). This yields the logical window sizes WS_ξ and WS_η which appear in equation (6) (figure 10(c)). As per the original aim of applying image mapping to avoid any overlap between correlation windows and identified interfaces, logical interrogation windows are then shifted to avoid overlap with the horizontal line at η  =  0 referencing the interface boundary (figure 10(d)). The shifted logical positions correspond with new values of the Jacobian and therefore new window sizes on the physical plane. Window dimensions on the logical plane are then re-calculated and the process is repeated until none of the correlation windows extends below η  =  0. This is typically satisfied in only 2 iterations.

Zoom In Zoom Out Reset image size

The second implementation of the window sizing algorithm, known as image mapping with vector relocation (IMVR), is a modification of IM, taking advantage of the vector relocation technique (Usera et al 2004). Correlation windows in the logical plane are now allowed to exhibit a certain degree of overlap with η  =  0 and, contrary to shifting the window in its entirety, are cut-off followed by a repositioning of the resulting velocity vector. The first steps of the second approach thus remain identical to those of the IM method (figures 10(a)–(c)). In IMVR, the portion of those correlation windows extending below η  =  0 is now masked from the process of cross-correlation and the resulting vector is consequentially relocated in the centre of the remaining area (figure 10(d)).

It is clear that the combination of the additional operations underpinning image mapping must come at some elevated computational cost. A distinction can be made between one-off additional effort for the generation of the numerical mesh throughout the multi-grid interrogation process and a persistent one due to the mapping of particle images (figure 10). The mesh is generated once and remains invariant in case of stationary interfaces whereas the repeated effort related to the mapping, instead, cannot be prevented. The PIV image could be deformed in its entirety only once before the analysis. However the authors strongly discourage the reader to take this approach, especially when highly curved surfaces are considered due to the resulting strongly deformed particle images (figures 9(c)–(e)). The suggested solution of enhanced re-sampling (figures 9(d)–(f)), on the other hand, requires the interpolation of both numerical mesh data and image grey levels for every interrogation area. Although this process can be quickly accomplished by means of quintic splines, it nevertheless typically augments the overall processing time by 3 to 6 times the time spent on the cross-correlation operator, depending on the density of the numerical mesh and its wall-normal extent. This additional load, however, only applies to the Chimera region where image mapping is involved and does not affect all the correlation windows.

Independent of the adopted numerical mesh, the transformation from physical to logical plane must be carried out through an interpolation (see section 4.2). Without regard to the PIV algorithm in which image mapping is implemented, this inherently produces a numerical error for both vector components and locations when converted and is non-negligible although often overlooked in existing literature. Quantification of this error is not straightforward though because of the lack of explicit analytical expressions specifying the transformations. The authors have therefore performed a simplistic analysis by consecutive mapping of spatially distributed points from the physical plane to the logical plane (inverse transformation) and back (direct transformation) as shown in figure 7. The error heuristic, ε_transf, is then defined as the difference between the initial point locations and their doubly-converted position, as expressed by:

$$\begin{eqnarray}&&{{\varepsilon}_{\text{transf}}}=\left\Vert\begin{array}{*{35}{l}} {{x}^{\ast}}\left({{\xi}^{\ast}}\left({{x}_{k}},{{y}_{k}}\right),\,{{\eta}^{\ast}}\left({{x}_{k}},{{y}_{k}}\right)\right)-{{x}_{k}} \\ {{y}^{\ast}}\left({{\xi}^{\ast}}\left({{x}_{k}},{{y}_{k}}\right),\,{{\eta}^{\ast}}\left({{x}_{k}},{{y}_{k}}\right)\right)-{{y}_{k}} \end{array}\right\Vert\end{eqnarray} \tag{ 8 }$$

where x^*(·), y^*(·) and ξ^*(·), η^*(·) are respectively the interpolation of the numerical solutions for the physical and logical coordinates of the mesh, as described in section 4.2, and x_k and y_k are sampling points distributed over the physical plane. The contour of ε_transf is depicted for a hyperbolic mesh generated over a Joukowsky aerofoil in figure 12(a), in logarithmic scale.

Zoom In Zoom Out Reset image size

The observed error with regard to distribution and magnitude is conservative as it considers both errors originating from the direct and inverse transformations, with the latter being less accurate due to the scattered nature of the required interpolation (see section 4.2). In reality, the two transformations are never used consecutively, for which reason this result must be considered as an upper error limit. Figure 12(a) illustrates the maximum expected error in the vector position to be less than 10⁻³ pixels, i.e. one order of magnitude below the commonly accepted limit in displacement accuracy (~0.01 pixels) obtainable with standard PIV (Huang et al 1999). Nevertheless, figure 12(a) shows that highest errors in transformation are encountered in proximity of the leading edge. The reader is reminded that this error analysis is independent from the PIV analysis as it only originates from the mapping of points using the mesh transformation. However, this general tendency is independent of the mesh density (see section 4.1) and is related to the divergence of the characteristic mesh lines in presence of surface curvature. The digression of the lines is evidenced by the strong derivatives of the Jacobian shown in figure 12(b), which increase in magnitude as the distance from the aerofoil increases. The same effect can be observed on the trailing edge too, where mesh lines are instead converging. Unfortunately, this additional error will occur exactly where velocity gradients and image deformation are stronger, contributing overall to a deterioration of the measurement accuracy in these regions.

Mapping of velocity components, in theory, should be less prone to error than mapping of points, as it involves a simple rearrangement of equation (4). The Jacobian in this equation is evaluated using derivatives of the structured mesh data on the logical plane and this step is therefore more accurate than the inverse transformation of the points, with negligible error throughout the entire PIV algorithm. However, the conversion of velocity components is affected by a different source of error which is not related to the interpolation but rather to an amplification of the standard PIV error connected to the use of the Jacobian itself. This amplification can be easily illustrated by considering a simple polar to Cartesian transformation; x  =  ρ · cos(θ) and y  =  ρ · sin(θ). When the displacement field on the logical plane is subjected to an error ε_log, the transformed error on the physical ε_phy scales with radial distance

$$\begin{eqnarray}&&\left[\begin{array}{c} \text{d}x \\ \text{d}y \end{array}\right]+{{{\underline{\varepsilon}}}_{\text{phy}}}=\left[\begin{array}{c} \cos \left(\theta \right) & -\rho \sin \left(\theta \right) \\ \sin \left(\theta \right) & \rho \cos \left(\theta \right) \end{array}\right]\left(\left[\begin{array}{c} \text{d}\rho \\ \text{d}\theta \end{array}\right]+{{{\underline{\varepsilon}}}_{\log}}\right)\\ &&\quad =[J]\left[\begin{array}{c} \text{d}\rho \\ \text{d}\theta \end{array}\right]+[J]{{{\underline{\varepsilon}}}_{\log}}\to \|{{{\underline{\varepsilon}}}_{\text{phy}}}\|=\|{{{\underline{\varepsilon}}}_{\log}}\|\sqrt{1+{{\rho}^{2}}}.\end{eqnarray} \tag{ 9 }$$

In case of the Joukowsky aerofoil, the norm of the Jacobian is depicted in figure 12(b) and a similar amplification of the error is observed in figure 12(c). Here, an artificially imposed probability density function of the magnitude of normally distributed random noise, σ,with a standard deviation of 0.1 pixels was superimposed on to the ideal velocity field near the Joukowsky aerofoil's leading edge on the physical plane. When the same noise is superimposed to ideal displacements on the logical plane and converted back into the physical plane, noise is amplified by a factor 1.6 (figure 12(c)). The reader is reminded that this error amplification is independent of the PIV algorithm used, but is related to the geometry of the mesh adopted for image mapping. This error source may affect the quality of the measurements, especially in regions of stronger mesh curvature, thus limiting the general application of image mapping to PIV. Even if no other cause of mesh-related error was present, a random noise of a measurement will be amplified according to the mesh curvature simply because of the use of image mapping.

When implementing image mapping within a PIV image interrogation process, special attention must be paid to the distortion of particle images as these establish the signal for cross-correlation. In line with findings reported in the previous paragraph, this error source will be strongly dependent on the mesh arching. To isolate and quantify the source of measurement inaccuracy pertinent to the mesh curvature, a special test case was designed using synthetic images. A single correlation window was deformed using a Cartesian (x, y) to polar (ρ, θ) transformation around a cylinder of variable radius ρ, ranging from 10 (ρ_min) to 700 pixels (ρ_max). To nullify the influence of velocity gradients, a constant shift of 1.25 logical pixels was imposed along the θ-axis in the logical plane. Using the Jacobian in equation (4) with conventional polar coordinate transformations corresponding displacements in the physical (δx, δy) and logical plane (δρ, δθ) were expressed as

$$\begin{eqnarray}&&\left\{\begin{array}{*{35}{l}} \delta x=-{{U}_{\inf}}\rho \sin \left(\theta \right) \\ \delta y={{U}_{\inf}}\rho \cos \left(\theta \right) \end{array}\,\quad \text{and}\quad \,\left\{\begin{array}{*{35}{l}} \delta \rho =0 \\ \delta \theta ={{U}_{\inf}}=1.25 \end{array}.\right. \right. \end{eqnarray} \tag{ 10 }$$

To decouple the error due to particle deformation from other effects connected to image mapping, the size of the interrogation window on the logical plane was kept constant while testing multiple cylinder radii (see figure 13). The logical height of the window WS_ρ was arbitrarily set to 35 logical pixels, while the logical width WS_θ was based on the largest angular extent (i.e. the circumference) of the minimum radius tested; WS_θ  =  2πρ_min  ≈  63. Physical particle images were thus sampled using a constant number of logical pixels.

Zoom In Zoom Out Reset image size

However, as the physical extent of logical pixels depended on their radial location within the correlation window, particle images were resampled following a 1:1 ratio nearest to the cylinder radius (i.e. one logical pixel per physical pixel) and under-sampled by a factor 1  +  35/ρ furthest. Accordingly, this test allowed assessing solely the effect of elliptical warping due to mesh curvature. Bias and random errors in measured displacement, noted as β and σ respectively, for different particle diameters following Monte-Carlo simulations involving 10⁴ synthetic images are plotted against the cylinder radii in figure 14. The error for a correlation window of the same size, non-deformed (i.e. infinite cylinder radius), shifted by the same rigid displacement, is also depicted for comparison as a dashed line.

Zoom In Zoom Out Reset image size

A strong dependency between error and cylinder radius can be observed in both figures 14(a) and (b). Comparing the undeformed (dashed lines) with the deformed bias error (solid lines), figure 14(a) indicates the implementation of image mapping to augment the displacement bias error. Independent of the particle image diameters, this additional bias error can be up to 8 · 10⁻³ pixels in case of the strongest deformation tested (i.e. minimum cylinder radius) and decreases monotonically as the cylinder radius increases, reaching a minimum of 8 · 10⁻⁴ pixels for the maximum radius tested. This behaviour is to be expected as the size of the cylinder affects the deformation of the particle images, increasing their elliptical shape. For all the particle diameters analysed, results are in accordance with the existing literature, with smaller particle diameters presenting a higher bias error, as they tend to degrade the accuracy of correlation peak sub-pixel fitting (Westerweel 2000). Inspection of figure 14(b) reveals a similar behaviour for the random error, with image mapping introducing an additional error in the measurements inversely proportional to the cylinder radius. For cylinder radii exceeding 100px random errors increase with particle image diameter. In case of strong deformations a rapid increase in random error can be observed for the smaller particle image diameter, yielding error levels exceeding those for larger diameters. The discrepancy between the asymptotic behaviour of the curves shown in figure 14 and the corresponding reference cases is due to numerical errors originating from the different interpolation processes intrinsic to image mapping.

Zoom In Zoom Out Reset image size

The authors would like to stress the fact that the mapping was here defined analytically. Presented errors in this subsection therefore reflect the influence of mesh curvature only and not mesh type. When using image mapping, a decrease in the measurement accuracy can thus be expected where mesh curvature is stronger merely due to the particle image deformation. Notwithstanding, actual advantages of image mapping will be case dependent (geometry shape, mesh type, flow gradients, correlation window distribution, etc) and may outweigh the above error.

In this section, image mapping has been embedded within PIV image analysis processes and compared with alternative advanced PIV interrogation routines. Numerical simulations presented afore have indicated measurement errors in particle image displacement, inherent to image mapping, to be strongly affected by the mesh curvature. To assess the limits of applicability of image mapping and investigate the maximum extent of these errors, a Joukowsky aerofoil characterized by a strongly curved surface was therefore selected as a test case representative of aeronautical applications. Synthetic images were generated with seeding consisting of randomly distributed Gaussian shaped particle images with diameters selected from a normal probability (3 pixels mean, 1 pixel variance) and a concentration of 0.1 particles per pixel. Images had a bit depth of 8 with mean particle intensities of 120 and standard deviation of 40 grey levels. The analytical inviscid flow around a Joukowsky airfoil enabled the direct evaluation of resulting measurement errors in case of realistic aerodynamic shapes containing regions of high surface curvature.

Five combinations of image interrogation procedures (summarised in table 2) have been tested for different settings such as interrogation window size, number of predictor-corrector iterations and image noise. The algorithms combined different implementations of correlation window distributions, vector relocation and image mapping. Being the measurement in the near vicinity of the object the main interest of this paper, the authors focused their attention on the implementation of image mapping combined with an adaptive sampling distribution as proposed by Theunissen et al (2008b). For this implementation, formally referred to as adaptive sampling (AS), correlation windows are adaptively located to allow an automatic selection of the window properties. The density of correlation windows increases automatically towards the interface, with a consequent increase in windows overlap too, as depicted in figure 20(a). The combination of adaptive sampling and vector relocation (Usera et al 2004), whereby no use is made of image mapping, will be referred to as ASVR (adaptive sampling vector relocation). This approach will serve as a basis for the comparison of the performances with more advanced techniques. The image mapping algorithm, as described in section 4.4, can be implemented with or without vector relocation, yielding the ASIM (adaptive sampling image mapping, see figure 10(d)) and ASIMVR (adaptive sampling image mapping with vector relocation, see figure 10(e)) approaches. An adaptation of the algorithms on a standard Cartesian grid is also proposed for the sake of comparison, referred to as CSVR (cartesian sampling vector relocation) and CSIM (cartesian sampling image mapping). An overview of the conceptually different methodologies is provided in table 2.

In the following, plots depict the total measurement error ε averaged over 1000 synthetic images. For each velocity field the reconstruction error is quantified by interpolating the measured displacements, dx_m and dy_m, on a pixel-wise grid and evaluating the difference with the theoretical, analytical, displacements, dx_a and dy_a, in the corresponding location; ε  =  ((dx_m  −  dx_a)²  +  (dy_m  −  dy_a)²)^0.5. To ease comparison between the different interrogation approaches, the ratio between errors is presented on a logarithmic scale; log₁₀(ε_ASIM/ε_ASVR), log₁₀(ε_ASIMVR/ε_ASVR), and log₁₀(ε_CSIM/ε_CSVR). Values close to zero indicate the two methods to have the same performances (error near unity ratio). This can be expected outside the Chimera area where no image mapping is invoked and differences between methodologies are marginal.

In the absence of image noise, figure 15 shows the comparison of the error between ASIM and ASVR performed in 3 predictor-corrector iterations with window sizes reduced from 90 to 25 pixels. The spatial distribution in error presents patterns which can be attributed to the process of image mapping and to the characteristics of the numerical mesh. Starting from the leading edge, the highest error is concentrated in two agglomerations above and below the aerofoil. In these areas the Jacobian of the mesh reaches its maximum (figure 12(a)), amplifying random fluctuations of the velocity field and increasing the error as described in section 5.2. Moreover, due to the low radius of curvature of the aerofoil (10 pixels), this image area is affected by the strongest particle deformation and an additional error in terms of bias and random is expected, as already shown in figure 14.

On the edge of the airfoil ASIM persistently attains a higher error. This consistent behaviour highlights an obstacle to the application of image mapping which has never been addressed in literature before. In fact, ASIM maintains the interrogation window size constant in proximity of a boundary and, as a consequence, the minimum distance between a displacement vector and the aerofoil boundary is at least half the minimum window size. Conversely, with ASVR vectors are repositioned to the centre of gravity of the correlation window portion occupied by the seeded flow. This yields a minimum distance of one-quarter of the minimum window size. Moreover, due to the lack of noise in the synthetic images, smaller interrogation windows are always associated to a smaller error. ASVR, whereby parts of the correlation windows overlapping with the aerofoil are masked automatically, results in analysing correlation windows which are smaller than the minimum set by the user. Therefore, in this particular case, ASVR produces better vectors which are closer to the aerofoil. It will be shown that in the presence of image noise this reasoning no longer holds.

Moving to the aft of the aerofoil, a behaviour in error similar to the leading edge is noticeable on the trailing edge. In proximity of the edge, a strong reduction in error is achieved with ASIM, yet accuracy rapidly reduces. The observed behaviour can be explained considering three different phenomena; the process of window repositioning on the logical plane described in section 4.4; the adaptive point distribution and the interpolation effect of the predictor-corrector implemented in the algorithm. As already described in section 4.4, correlation windows on the logical plane are automatically repositioned in order to keep the seeding density constant and to ensure they all have exactly the same area. This process inevitably produces a clustering of sampling points on the edges of the Chimera region, as shown in figure 20. During the predictor-corrector process, the interpolation of velocity data in those areas of higher data density will produce a lower error of the velocity predictor and consequently a final error of ASIM which is lower than ASVR. The reader is reminded that this lower error is not related to the enhancement of the images produced by image mapping itself but rather to a mere effect of the interpolation due to the inherited clustering of data points.

The only difference between the implementation of image mapping in ASIM and ASIMVR is encountered near the airfoil edge (figure 16). ASIMVR does not shift the window position to keep the area constant, but in an analogous way to ASVR, masks the correlation window area overlapping with the aerofoil on the logical plane and relocates the vector in the centre of gravity of the captured flow region. The effect shows an error of ASIMVR on the edge of the aerofoil which is on par with ASVR. In this case the application of image mapping results in no improvement of the measurement accuracy. However, it is important to stress that this result is only valid for the ideal case of images without artificial noise. The ASIMVR algorithm does not produce any clustering of points since no repositioning of logical windows is performed. Nevertheless, comparison of figure 16 with figure 15 indicates a noticeable reduction of displacement error near the trailing edge. Despite the minimum interrogation window size imposed by the user, the actual area adopted by ASIMVR for the cross-correlation depends of many different conditions. On the trailing edge, the numerical mesh is logically cut along the aerofoil chord as shown in figure 11. Because of this discontinuity, interrogation windows close to the trailing edge analysed with ASIMVR will be smaller than windows positioned in the same area but analysed with ASVR. In absence of noise, this produces a smaller error for ASIMVR.

Zoom In Zoom Out Reset image size

With the purpose of understanding all the possible advantages and applications of image mapping, and proposing a more realistic case of images, the previous analysis has been repeated in the presence of artificial noise, in the form of a normal distribution with mean of 2.8% and standard deviation of 3.6% of the maximum intensity, added to the ideal images. The spatial distribution in displacement error produced by image mapping or standard adaptive PIV algorithms presents many similarities with the two previous cases. Now the relative importance of the error produced by the mesh transformation is dominated by the uncertainty brought by the artificial noise. Differences inherent to the image mapping technique are amplified accordingly.

Figure 17(a) presents the error ratios between ASIM and ASVR and shows features of high error similar to those in figure 15. These are again due to the strong deformations in the numerical mesh. However, a remarkable reduction in error with ASIM can be noticed in proximity of the aerofoil boundary. Because of the presence of artificial noise, a smaller window size, such is the case in ASVR, does not necessarily produce a better quality vector. Consequently, ASIM is able to yield measurements with higher accuracy by ensuring that all correlation windows respect the minimum window imposed by the user. This result is in accordance with the improvements already shown in literature (Jeon and Sung 2011) where a reduction of the error is expected in case of low curvature radius and images affected by artificial noise. Another area of smaller error in figure 17(a) is on the outer edge of Chimera which is ascribed to the higher density of sampling points.

Zoom In Zoom Out Reset image size

Figure 17(b) shows the error ratio for ASIMVR and presents similar features to ASIM. Because of artificial noise, the inherent characteristic of using smaller windows, which was beneficial in absence of noise, now becomes harmful. For this reason, the error of ASIMVR on the edge of the aerofoil reaches levels on par with ASVR. The same reduction in window size on the trailing edge due to the logical cut of the mask worsens the measurement precision. For the same reason, the error present on the outer area of the Chimera region is now negatively affected by the reduced window size and shows an increased level in figure 17(b) compared to figure 17(a).

With the objective of deepening the understanding in the conduciveness of image mapping in PIV algorithms, the analysis of the Joukowsky images has also been performed considering a single iteration thereby negating the predictor-corrector process. To allow the use of small correlation window sizes while maintaining a sufficient number of particle images per window, the maximum flow velocity magnitude was reduced to 1 px. Results are shown in figure 18 for both ASIM and ASIMVR and present once again the typical patterns of the numerical mesh error. However, given the very strong reduction of the velocity gradient magnitudes due to the very small displacement, there seem to be no actual advantages in using image mapping anymore. In fact, in spite of the artificial noise, quasi-static correlation windows are rotated, stretched and deformed by the image mapping without bringing any reduction in the gradients of the velocity field. This finding is in accordance with the results shown in figure 12(c); the mere application of the image mapping process to a static displacement field amplifies the random fluctuation increasing the displacement error.

Zoom In Zoom Out Reset image size

For the sake of completeness, a comparison with a structured Cartesian grid is also presented in figure 19. The figure shows similar patterns to the previous adaptive sampling, with the only differences being patchy artefacts due to the non-uniform distribution of the sampling points. The reader should note that the previous analysis, based on the adaptive point distribution, does not depend on the type of sampling distribution adopted, but rather on the peculiar aspects of image mapping and its algorithmic implementations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The numerical assessment has indicated the applicability of image mapping not to be as simple and straightforward as suggested in previous works (Nguyen et al 2010, Park et al 2015). Instead, the final outcome of the measurement strongly depends on a case-by-case basis, including the geometric shape of the interface and details of the algorithmic implementation. Essentially, two dominant effects determine the degree to wich image mapping can be beneficial; image deformation as a result of the imposed numerical mesh and the ability to maintain the image density in the vicinity of the object boundary.

Although the numerical mesh introduces a distortion of the particle images which is detrimental to the correlation accuracy as illustrated in figure 14, depending on the alignment of the numerical mesh with the flow streamlines, the resulting artificial displacement gradients can be either constructive or unfavourable. In case of the cylinder flow presented in figure 13 the mesh transformation straightens the curved streamlines in the logical plane, thereby minimising gradients. Contrarily, in the case of the Joukowsky airfoil near the leading and trailing edge, the misalignment between the governing flow and imposed numerical mesh combined with the strong gradients in the Jacobian (figure 12) led to a consistent deterioration of the PIV accuracy. Unfortunately, these effects cannot be quantified a priori as the underlying flow field is exactly the subject of the analysis.

The results presented in section 5.3 attest that image mapping does not generally yield an in increase in accuracy, even when incorporated in advanced algorithms. Instead, the results suggest more simplistic interrogation methodologies to outperform those including image mapping. Only for very specific imaging and flow conditions marginal improvements could be attained by implementing image mapping. In the near vicinity of the Joukowsky airfoil the ability of image mapping to avoid any overlap between the interrogation areas and the object meant that a more reliable cross-correlation could be achieved in the presence of image noise (figures 17 and 19). However, this improvement was irrelevant in the absence of image noise (figures 15 and 16). The potential gain in resolution consequently depends on two factors. First, the quality of the PIV images in proximity of the edges (i.e. severity of reflections, intensity of image noise, concentration of seeding, etc). This can be partly influenced by adequate pre-processing methodologies (e.g. Mendez et al (2017)). Second, the final correlation window size, which will typically involve trial and error. If the user selects a too large size to ensure a reliable correlation, no gain will be noticeable due to a too coarse resolution. Should the user opt for a too small correlation window size, independent of the analysis approach, cross-correlation will yield sporadic vectors.