Fourier-based interpolation bias prediction in digital image correlation

Abstract: Based on the Fourier method, this paper deduces analytic formulae for interpolation bias in digital image correlation, explains the well-known sinusoidal-shaped curves of interpolation bias, and introduces the concept of interpolation bias kernel, which characterizes the frequency response of the interpolation bias and thus provides a measure of the subset matching quality of the interpolation algorithm. The interpolation bias kernel attributes the interpolation bias to aliasing effect of interpolation and indicates that high-frequency components are the major source of interpolation bias. Based on our theoretical results, a simple and effective interpolation bias prediction approach, which exploits the speckle spectrum and the interpolation transfer function, is proposed. Significant acceleration is attained, the effect of subset size is analyzed, and both numerical simulations and experimental results are found to agree with theoretical predictions. During the experiment, a novel experimental translation technique was developed that implements subpixel translation of a captured image through integer pixel translation on a computer screen. Owing to this remarkable technique, the influences of mechanical error and out-of-plane motion are eliminated, and complete interpolation bias curves as accurate as 0.01 pixel are attained by subpixel translation experiments.


Introduction
Digital image correlation has evolved into a reliable and flexible full-field optical metrology [1][2][3], and this technique has found numerous and ever-increasing applications in scientific and engineering fields [4][5][6][7][8][9]. The combination of digital image correlation with other techniques such as microscopy, holography, and fringe projection has further extended the scope of its application and highlighted its potential [10][11][12].
Digital image correlation retrieves full-field deformation by matching subsets in reference and target images. The displacement and strain errors of digital image correlation were first reported by Bruck et al. [13]. Detailed evaluation presented by Sutton et al. indicated that the displacement measurement error of digital image correlation is approximately 0.01 pixels [14]. This metrological performance is influenced by many factors [15], such as image noise [16], shape function [17], correlation criterion [18], optimization algorithm [19], and interpolation accuracy [20]. Among these issues, interpolation plays an essential role. In the 1980s and 1990s, sinusoidal-shaped systematic error in digital image correlation was observed, which was periodic with a period of 1 pixel and an amplitude possibly exceeding 0.06 pixel [21,22]. Schreier et al. attributed this periodic error to imperfect interpolation [20]. To achieve subpixel accuracies, grayscale must be evaluated at non-integer locations in digital image correlation. Non-ideal interpolation will lead to this systematic error, which is called interpolation bias. As a result of sinc function's infinite support and slow decay, ideal interpolation cannot be performed in practice [23], and thus interpolation bias is inevitable. Schreier et al. also pointed out that interpolation bias can be evaluated by correlating subpixel shifted images [20]. Since then, diverse approaches have been presented to reduce interpolation bias. These approaches can be classified into three categories: more sophisticated interpolation algorithms, low-pass image filtering, and stochastic integration [24]. High-order interpolation algorithms generally produce better results but require long execution time, and therefore it is preferable to optimize the interpolation basis [1]. Luu et al. introduced an inverse gradient weighting version of the BSpline interpolation algorithm, thereby enhanced accuracy [25]. Pre-filtering the speckle image with a low-pass filter has also been suggested [20,26,27]. The low-pass filter can suppress interpolation bias dramatically, but consequent blur will reduce the spatial resolution of digital image correlation. In the digital image processing community, Rohde et al. proposed the usage of an integral form of the correlation function instead of the summation form; a stochastic integral was employed to approximate the continuous integral of the correlation function, thus enhancing accuracy [24].
Despite the large number of approaches that have been suggested to decrease interpolation bias, the explicit dependence of this bias upon the speckle pattern and the interpolation algorithm is still unknown, and the sinusoidal-shaped curves remain unexplained. In addition, reducing interpolation bias by speckle pattern optimization demands knowledge of interpolation bias [28,29]. However, quantitative discussions are rare in the literature because interpolation is difficult to tackle mathematically and correlation is inherently nonlinear. In Wang's pioneering work [30], the interpolation bias for linear and cubic interpolation was predicted using sample values. However, his theory is inapplicable to high-accuracy interpolation algorithms such as BSpline [20], OMOMS [25], which have an infinite convolution kernel. Moreover, because his formulae include the difference between interpolated and real grayscales, they are hard to employ in theoretical analyses such as interpolation quality assessment and speckle pattern optimization.
Another challenge is the purely experimental measurement of interpolation bias. Subpixel translation experiments have seldom been reported because the interpolation bias can easily be overwhelmed by other experimental uncertainties such as stage error, out-of-plane motion, and image noise. Recently, Mazzoleni et al. reported the best result obtained to date, but nevertheless they failed to attain the complete curve due to accidental and variable vibrations of the camera and its non-isolated support [26]. Measuring interpolation bias accurately by purely experimental techniques is still a challenge.
At present, numerical methods are prevalently utilized to investigate interpolation bias. The conventional approach is to correlate a series of subpixel-shifted speckle images generated by synthetic speckle [19], FFT [20] or binning methods [31], which are inconvenient and time-consuming. If the exact dependence of interpolation bias upon the speckle pattern and the interpolation algorithm can be determined analytically, the interpolation bias can be predicted readily, considerable physical insight into the problem will be gained and studies of speckle pattern assessment and interpolation algorithm optimization will be benefited.
Similar problems arise in image processing [24], computer vision [32], remote sensing [33], and flight time estimation [34]. Despite distinct backgrounds and distinct correlation criterions, imperfect interpolation inevitably introduces systematic error, which occurs in various research fields which are not confined to digital image correlation.
Previous studies on interpolation bias generally utilized spatial methods. Considering the sinusoidal nature of the interpolation bias, a frequency method is probably preferable. Facilitated by the Fourier method, this paper provides a thorough discussion of interpolation bias.
This paper is organized as follows. Section 2 introduces the principles of digital image correlation and convolution-based interpolation, after which the analytical formulae for 1-D interpolation bias are deduced through Fourier analysis and confirmed by numerical experiment. Section 3 extends the theoretical analysis to high-dimensional space. Exploiting the speckle spectrum and interpolation transfer functions, a novel interpolation bias prediction approach is proposed. Subset size effects are analyzed, and theoretical predictions are confirmed by numerical simulations and actual experiments. In the course of the experiments, a novel subpixel translation technique was developed. Section 4 draws conclusions from this work.

Principle of digital image correlation
The principle of digital image correlation is to match subsets in reference and target images using a correlation criterion. The choice of correlation criterion can be loosely divided into two classes: cross-correlation and sum of squared differences [2]. This work uses the sum of squared differences criterion.

Principle of convolution-based interpolation
Interpolation is the process of constructing a function where the interpolation coefficient k c is chosen so that the original and reconstructed functions have identical values at the sample points. Equation (1) is equivalent to: The relationship between ( ) r x and ( ) x ϕ is discussed in [35]. The Keys, BSpline and OMOMS interpolation algorithms are all convolution-based interpolation methods [23,36,37].

1-D interpolation bias analysis
It is illuminating to analyze the 1-D situation first. Figure 1 [35].). Denote a Fourier transform pair by ←⎯→ F ; then the process described above can be represented symbolically as:  is the comb function; for a more detailed mathematical derivation of ( ) g ν , the reader is referred to [38]. Using the Poisson summation formulae again yields: Substituting Eq. (6) and Eq. (8) into Eq. (5) gives: Equation (9) The calculated result 0 e u u u = + . Recognizing that e u is generally small, to a first-order approximation, The explicit dependence of the interpolation bias e u upon the reference function ( ) f ν and the interpolation transfer function ( ) ϕ ν can be derived as follows: Equation (12) is referred to as the full estimate of interpolation bias in this work. Because the full estimate of the interpolation bias is complex and cumbersome, it is difficult to derive a physical meaning from the above expression. Moreover, it is computationally expensive due to its complexity, and therefore simplification is desirable. Recognizing that is a periodic function with period 1 and contains a direct component, and recognizing that ( ) ϕ ν is an approximation of an ideal low-pass filter [see can be approximated by (let 0 k m n = = = ): is a periodic function with period 1 analogous to and therefore can be expanded into a Fourier series, ( ) However, this function is still too complex, and hence a better strategy might be to simplify it under specific conditions. Digital systems demand a sufficient sampling frequency to meet the requirements of the sampling theorem, and therefore the band-limited hypothesis is a reasonable assumption. To eliminate phase distortion, the interpolation basis is invariably an even function so that the transfer Utilizing the band-limited hypothesis and the fact that ( ) ϕ ν is real in practice, Eq. (12) can be simplified to (let 0 k n = = ): Equation (14) is referred to as the band-limited approximation of the interpolation bias.
associates with the squared sum of the gray gradient because by Parseval's theorem incorporates two parts: the first where C is a constant determined by the reference function ( ) f ν and the interpolation algorithm ( ) ϕ ν exclusively. Equation (15) will be referred to hereafter as the sinusoidal approximation of the interpolation bias. Equation (15) is exclusively determined by the interpolation algorithm; it has been called the interpolation bias kernel by the author. This interpolation bias kernel plays a central role in this work and is of significant importance because it characterizes the bias response at specific frequencies, thus providing a measure of the subset matching quality of various interpolation algorithms. Figure 2 illustrates the interpolation basis [see Fig. 2(a)], the transfer function [see Fig. 2 so interpolation accuracy generally follows the order OMOMS > BSpline > Keys. Higherorder BSpline algorithms show less bias response in the low-frequency region and therefore yield better results. These phenomena have been noted before [1,20,25], but this work provides a quantitative explanation and reveals their inherent nature. The preceding discussion indicates that the interpolation bias kernel quantifies whether an interpolation algorithm is appropriate for digital image correlation.
Imperfect interpolation involves blurring and aliasing effects: ( ) Theoretical formulae can predict interpolation bias using the reference function and the interpolation transfer function, thus avoiding the generation and correlation of a series of subpixel-shifted speckle images. This theory can be employed to optimize the interpolation algorithm further. OMOMS is optimum in the asymptotic sense [37], but this does not imply that it is the best choice for digital image correlation. Moreover, different interpolation methods can be optimized for different types of speckle patterns. This theory can explain why a low-pass filter can decrease interpolation bias and can be used to direct the selection of filter size. Furthermore, it can be used for speckle pattern assessment and design.

1-D Simulation
To demonstrate the validity of the analytical results, numerical experiments were carried out. 1-D speckle patterns were generated by Zhou's algorithm [39]. The reference function where N is the total count of Gaussian speckles, k x is the center coordinate of speckle k, r is the Gaussian speckle radius, and k I is the intensity level of speckle k. Its corresponding spectrum is: The spectrum consists of two parts:  The intensity level was fixed as 1 k I = ; speckles were generated in the interval (−50,50), with a speckle density of 65%, and speckle patterns with radii 1.5, 2.0, and 3.0 were produced. Figures 3(a1)-3(c1) illustrates the reference function and sample values. Exact sampled values were used to remove quantization error. The speckle centers were subjected to a translation of 0.05 unit each time, and 20 shifted speckle patterns were obtained. Correlations were implemented using a zero-order shape function by the forward Gauss-Newton method [40], zero order shape function is not fundamental because high order shape function will not induce additional systematic error [41]. The convergence criterion was that the successive iterative increase be less than 1 × 10 −10 . The subset was chosen as [-60, 60] and consisted of a total of 121 points. Keys, cubic BSpline, and cubic OMOMS interpolations were used; the theoretical predictions incorporated full estimation [Eq. (12)] and sinusoidal approximation [Eq. (15)]. The infinite sum was approximated by a discrete sum from −10 to 10, and the numerical integral was implemented using the trapezoidal method with step 1 × 10 −6 . Figures 3(a2)-3(c2) shows the digital image correlation error, the full estimate, and the sinusoidal approximation of the interpolation bias for Keys, cubic BSpline, and cubic OMOMS. Because the interpolation bias of Keys was much larger than the others, the subfigures show results only for BSpline and OMOMS. It is evident that for the Keys method, the full estimate shows excellent agreement with the digital image correlation error, whereas the sinusoidal approximation loses details of the interpolation bias curves, but shares the same order of magnitude. For the cubic BSpline and OMOMS methods, the interpolation bias curves show excellent agreement with all theoretical results. This demonstrates that the method proposed here can be pronounced effective. Obviously, decreasing r induces an increase in e u . For these representative speckle patterns, the interpolation bias follows the order Keys >> BSpline > OMOMS. These phenomena can be explained by the interpolation bias kernel ( ) ib E ν . Because speckles with sharper edges tend to larger interpolation bias [20], and the denominator of the constant C in Eq. (15) is associated with the sum of squared grayscale gradients, the numerator of C dominates its behavior. Hence, the interpolation bias was determined mainly by the integral of the product of the power spectrum ( )   , the vast majority of the energy is concentrated in the low-frequency region; because OMOMS has a better bias response at lower frequencies, the performance of OMOMS is superior to that of BSpline. Nevertheless, this is not the truth when high frequencies dominate the energy. This 1-D simulation has established the validity of the theoretical results presented earlier.

Interpolation bias prediction in high-dimensional space
Digital image correlation and digital volume correlation, which are both of practical importance, correspond respectively to two and three dimensions, and therefore it is essential to expand the preceding 1-D discussion to higher-dimensional space. In higher dimensions, coordinates, displacements, and frequencies are vectors rather than scalars. For Ndimensional space, suppose that the reference function is ( ) f x , its corresponding Fourier transform is ( ) f ν , the real displacement is 0 u , the interpolation bias is e u , the interpolation basis is ( ) ϕ x , its interpolation transfer function is ( ) ϕ ν . Taking an approach similar to the 1-D case and substituting variations for derivatives, the interpolation bias e u satisfies: If the fundamental frequency is considered exclusively: Equation (22) is similar to Eq. (15), where

( )
, ib x y E ν ν is the interpolation bias kernel in two dimensions and the aliasing effect along the y-direction is neglected.

Numerical verification
The preceding theoretical results can be employed to predict the interpolation bias. Nevertheless, experimental speckle images are captured by digital cameras, and therefore their original spectrum is unknown. To approximate the original spectrum, it is proposed to utilize the discrete Fourier transform of the subsets. Clearly, the accuracy of the spectrum approximation depends on the choice of subset size. This section analyzes the effect of subset size and compares this novel approach with the conventional one. Experimental speckle image of a concrete cylinder was captured with a resolution of 2448 × 2048. A series of 20 subpixel-shifted images were produced numerically from the original image using the FFT method. Successive images correspond to a shift of 0.05 pixel. Square subsets with side lengths of 51, 101, 251, 501 pixels were chosen as illustrated in Fig. 5(b). The speckle pattern of each subset is shown in Figs. 6(a1)-6(d1), magnified to identical size for clarification. The corresponding autocorrelation functions are shown in Fig. 5(a). They indicate that the fluctuation of larger subsets is less obvious and that the autocorrelation width is largely independent of subset size. Digital image correlation was implemented by the Gauss-Newton method with a zero order shape function. The Keys and cubic BSpline interpolation methods were used. Interpolation bias curves were obtained through correlation. A sinusoidal approximation [Eq. (22)] was used to predict the interpolation bias, and the original spectrum was approximated by the discrete Fourier transform. Figures 6(a2)-6(d2) shows the discrete Fourier transform of each subset. The continuous integral in Eq. (22) is approximated by the numerical integration function trapz in MATLAB as we show in Code File 1 (Ref. [42].). Figure 5(c) implies that the interpolation bias remains roughly invariant as the subset size varies. The predictions using Eq. (22) are larger than the correlation results, but the discrepancy decreases as the subset size increases. The phenomenon that theoretical predictions for larger subsets are more accurate can be explained as follows: first, it can be inferred from Figs. 6(a2)-6(d2) that discrete Fourier transforms of larger subsets approximate the original spectrum better; second, a larger subset size leads to a smaller numerical integration step and more integral points, so that the numerical integral is more accurate. Because the interpolation bias remains roughly invariant as subset size varies, larger subsets are recommended for theoretical prediction. However, when small subsets are chosen, safe predictions are obtained.  Recently, the use of numerically designed speckle patterns was proposed; these patterns consist of randomly positioned circles [5]. Numerically designed speckle is more suitable for large fields of view, and a Gaussian pre-filter can decrease both interpolation bias and random error [26]. Numerically designed speckle patterns were subjected to further verification. Speckle patterns were generated, printed, and affixed onto the surface of a rigid plate. The plate was placed near to, then far from the camera to capture speckles with different radii. An IDS camera with resolution 2048 × 2048 was used. Figures 7(a1)-7(b1) shows the captured image, in which the frames indicate the correlation subsets; close-ups can be found in Figs. 8(a1)-8(b1). The speckle density is roughly 65%. The speckle radius is about ten pixels for coarse patterns and four pixels for fine patterns. The autocorrelation functions of subsets are shown in Figs. 8(a2)-8(b2). They suggest that the distance between the crest and the first trough is approximately equal to the speckle diameter, meaning that finer speckles tend to produce a sharper autocorrelation function. Comparing the autocorrelation functions for numerically designed speckles [see Figs. 8(a2)-8(b2)] with those of traditional spray-painted random patterns [see Fig. 5(a)], the fluctuations of numerically designed speckles are more significant. The discrete Fourier transforms of the subsets are illustrated in Figs. 7(a2)-7(b2) and indicate that the spectra of numerically designed speckle patterns show a ring-like shape which can be explained by the Fourier transform of a circle function, quite unlike traditional speckles. Traditional speckles are highly random patterns without well-defined particles. It is also apparent that that the energy in coarse patterns is concentrated at low frequencies, whereas fine patterns have more energy in the high-frequency domain.    Fig. 7(b3)]. The interpolation bias of BSpline is much smaller than that of Keys, and the interpolation bias increases as the high-frequency components of the speckle patterns increase.
Rigid motions were mainly concerned in this work. If the deformations are not rigid motion, a number of digital image correlation techniques use the displacement gradient to describe the kinematics of the subset [15]. The existence of displacement gradients will lead to the coupling of the interpolation errors and the displacement fields. This coupling effects need further research.

Subpixel translation experiment
This subsection describes the subpixel translation experiment which was conducted to measure interpolation bias experimentally. To eliminate mechanical error and out-of-plane motion, a novel subpixel translation technique, implementing subpixel translation of the captured image through integer pixel translation on a computer screen, has been developed. Microsoft Surface Pro3 was used to visualize speckle patterns. The PPI of Surface Pro3′s screen is 216, so that a screen pixel corresponds to 117.6 μm. A JoinHope camera with 640 × 480 resolution was used. Its sensor pixel size is 7.1 μm. The focal length of the lens was 12 mm. The experimental setup is shown in Fig. 9(b). The experimental procedure was the following. First, the camera and the computer were placed on a vibration-isolated table, and the distance between the camera and the computer screen was adjusted so that 1 pixel on the screen corresponded to roughly 0.1 pixel in the image. Second, a laser pen was positioned just above the camera. To confirm that the laser direction was parallel to the camera optical axis, the laser pen was justified until the corresponding laser point was in the middle of the image. After this, a CD was affixed on the screen, and the screen was adjusted so that the incident light coincided with the reflected light, which guaranteed the perpendicular of the screen and the optical axis. Then the laser pen and CD could be removed. Third, speckle patterns were displayed on the screen; these speckle patterns were translated by 1 pixel each time [see Fig.  9(a)]. To decrease the influence of image noise, 100 frames were captured and averaged each time. Fine, medium, and coarse patterns were captured as illustrated in Figs. 10(a1)-10(c1); the corresponding speckle diameters were approximately 4, 6, and 8 pixels.
The Keys and cubic BSpline interpolation algorithms were used to retrieve the speckle displacement. Translation curves for different speckle radii were measured as shown in Figs. 10(a2)-10(c2), and the interpolation bias was evaluated by linear fit. Figures 10(a3)-10(c3) illustrates the digital image correlation results and the sinusoidal approximations of different speckle patterns. Because the environmental noise is relatively small compared to the bias of Keys (it can exceed 0.01 pixel), the interpolation bias of the Keys algorithm showed good agreement with theoretical predictions. For BSpline, the interpolation bias is small, and therefore it is susceptible to noise. However, the theoretical predictions of BSpline were of the same order of magnitude as the experimental results. These results also indicated that fine patterns result in large interpolation bias, which is consistent with the preceding discussion. For instance, the bias of Keys interpolation was about 0.02 pixel for fine patterns and decreased as the speckle radius increased. For medium patterns, the bias was 0.015 pixel, and the bias for coarse patterns was even smaller. In summary, the experimental results in this research have shown good agreement with theoretical predictions, indicating that the proposed method represents a reliable prediction approach.

Conclusions
The emphasis throughout this paper has been on analysis of interpolation bias in digital image correlation. An interpolation-bias prediction algorithm has been presented and verified by numerical simulation and actual experiments. The following conclusions can be drawn: 1) The explicit dependence of interpolation bias upon speckle pattern and interpolation algorithm has been characterized, analytic formulae for interpolation bias have been deduced, and the well-known sinusoidal-shaped curves of interpolation bias have been explained.
2) The concept of interpolation bias kernel has been introduced. An interpolation bias kernel characterizes the bias response at specific speckle frequencies, thus providing a measure of the subset matching quality of the interpolation algorithm. The properties of the interpolation bias kernel indicate that high-frequency components are the major source of interpolation bias. Further investigation attributed the interpolation bias to aliasing effect of interpolation.
3) A simple and effective interpolation bias prediction approach has been proposed. This approach exploits speckle spectrum and interpolation transfer function to predict interpolation bias. Significant acceleration has been accomplished compared to traditional methods, the effect of subset size has been analyzed, and both numerical simulations and experimental results show good agreement with theoretical predictions.

4)
A novel experimental translation technique has been developed which implements subpixel translation of a captured image by integer pixel translation on a computer screen. Owing to this noteworthy technique, the influences of mechanical error and out-of-plane motion are eliminated and a complete interpolation bias curve with an accuracy of 0.01 pixel can be attained by experimental subpixel translation. The primary motivation for this work is to standardize speckle pattern. This work can be applied not only to interpolation bias prediction, interpolation algorithm optimization, and speckle pattern optimization, but also to image processing and computer vision.