## Abstract

In digital image correlation (DIC), the noise-induced bias is significant if the noise level is high or the contrast of the image is low. However, existing methods for the estimation of the noise-induced bias are merely applicable to traditional interpolation methods such as linear and cubic interpolation, but are not applicable to generalized interpolation methods such as BSpline and OMOMS. Both traditional interpolation and generalized interpolation belong to convolution-based interpolation. Considering the widely use of generalized interpolation, this paper presents a theoretical analysis of noise-induced bias for convolution-based interpolation. A sinusoidal approximate formula for noise-induced bias is derived; this formula motivates an estimating strategy which is with speed, ease, and accuracy; furthermore, based on this formula, the mechanism of sophisticated interpolation methods generally reducing noise-induced bias is revealed. The validity of the theoretical analysis is established by both numerical simulations and actual subpixel translation experiment. Compared to existing methods, formulae provided by this paper are simpler, briefer, and more general. In addition, a more intuitionistic explanation of the cause of noise-induced bias is provided by quantitatively characterized the position-dependence of noise variability in the spatial domain.

© 2016 Optical Society of America

## 1. Introduction

Digital image correlation (DIC) is a practical, flexible, and reliable optical metrology widely used for shape, motion, and deformation measurements [1–8]. This technique retrieves full-field displacements by matching subsets in different images. The matching process demands that the texture of the specimen contains abundant information content; therefore artificial speckle patterns are generally prepared. Since the qualities of the speckle patterns profoundly influence the metrological performance of DIC, both scholars and practitioners are concerned with the design and implementation of the optimal speckle patterns [9,10]. After long-term practice and research, several empirical rules have been summarized: the ideal speckle patterns should be isotropic, non-periodic, high information content, and good contrast [11]. However, for a precision metrology, the quantitative evaluation is more essential; quantitative quality assessment of speckle patterns is based on the evaluation of metrological performance of DIC, and requires exhaustive knowledge of DIC errors. Influenced by many factors [12,13], the displacement measurement errors of DIC are typically several hundredth pixels [14]; these errors can be classified in two categories: random errors (variability) and systematic errors (bias).

Systematic errors are mainly caused by under-matched shape function [15], interpolation error [16], and image noise [17]. Considering that the systematic errors due to under-matched shape function rely on the underlying displacement fields [15,18], this paper deals with the *ultimate error regime* exclusively [19], in which case the shape function is perfectly matched. In the absence of noise, imperfect interpolation will introduce sinusoidal-shaped systematic errors, which are referred to as *interpolation bias* [20]. In our previous work, a concept of *interpolation bias kernel* was introduced, the inherent nature of interpolation bias was revealed, and an interpolation bias prediction algorithm, which exploits both interpolation bias kernel and the information of power spectrum, was proposed [21]. Actual images are inevitably contaminated by image noise. Wang et al. indicated that, if image noise exists, the systematic errors of DIC comprise two terms: the first term is interpolation bias; the second term is *noise-induced bias*, which is due to the combination of non-ideal interpolation and image noise [17]. In the context of high noise level or low image contrast, the noise-induced bias is more significant than the interpolation bias. Unfortunately, the framework provided by Wang is merely available for *traditional interpolation* methods such as linear and cubic interpolation, but is not applicable to *generalized interpolation* methods such as commonly used BSpline [13,20] and OMOMS [22]. To overcome the limitations, Sutton et al. investigated the noise-induced bias through a concept called *interpolation phase error*. It is noted that the noise-induced bias is proportional to the noise energy and inverse proportional to the sum of squared intensity gradients [11]. Nevertheless, the dependence of noise-induced bias upon subpixel positions for generalized interpolation is still unknown.

The aim of the work is to present thorough analysis of noise-induced bias for convolution-based interpolation (both traditional interpolation and generalized interpolation belong to convolution-based interpolation). The paper is organized as follows. In part 2, theoretical analysis of noise-induced bias is presented. In part 3, the theoretical analysis is confirmed by numerical simulations. In part 4, actual subpixel translation experiment is performed to validate the correctness of the theoretical analysis. In part 5, the differences with previous finding are discussed, and a more intuitionistic explanation of the cause of noise-induced bias is proposed. In part 6, conclusions of this work are drawn.

## 2. Theoretical analysis of noise-induced bias

#### 2.1 Convolution-based interpolation, traditional interpolation, and generalized interpolation

Considering that the readers of digital image correlation community may not familiar with the interpolation theory developed by the signal processing community, this section provides a brief introduction to convolution-based interpolation, traditional interpolation, and generalized interpolation. The readers are referred to review papers such as [23–25] for details.

According to the Nyquist-Shannon sampling theorem, if the band-width *B* of a band-limited function *f*(*x*) satisfies *B*<0.5 (the Nyquist frequency), *f*(*x*) can be completely recovered by its samples *f*(*k*)

*sinc*function. Consequently, the

*sinc*function is generally approximated by a sinc-like function

*φ*

_{int}(

*x*), and accordingly the reconstructed function

*f*(

_{T}*x*) turns to

*φ*

_{int}(

*x*) is referred to as the

*interpolation basis*, and this interpolation model [Eq. (2)] is referred to as

*convolution-based interpolation*. An interpolation method is completely determined by its

*φ*

_{int}(

*x*). As

*φ*

_{int}(

*x*) is merely an approximation of the

*sinc*function, the reconstructed function

*f*(

_{T}*x*) is not identical to the original function

*f*(

*x*); the discrepancies between

*f*(

*x*) and

*f*(

_{T}*x*) depend on the choice of

*φ*

_{int}(

*x*).

In the early researches, it was common to choose *φ*_{int}(*x*) as a piecewise-polynomial with finite support for convenience. These interpolation methods are referred to as *traditional interpolation.* Two examples of traditional interpolation are linear and Keys interpolation [26]. Figure 1(a1) shows the interpolation bases corresponding to linear and Keys interpolation, superimposed on the desired *sinc* function. Due to the significant differences between the *sinc* function and these interpolation bases, the accuracies of these interpolation methods are frequently insufficient. In the frequency domain, interpolation transfer functions shown in Fig. 1(b) clearly demonstrate the poor accuracy of linear and Keys interpolation: it is evident that the differences between the ideal low-pass filter (the transfer function of *sinc* interpolation) and the transfers functions corresponding to linear and Keys interpolation are significant.

In order to enhance the interpolation accuracy, *generalized interpolation* was presented. The essential difference between the generalized interpolation and the traditional interpolation is the introduction of interpolation coefficients *c*(*k*). For generalized interpolation, the reconstruction function satisfies

*c*(

*k*) is the interpolation coefficients, and

*φ*(

*x*) is a finite-support piecewise-polynomial. To enforce exact interpolation, for integer position

*k*

_{0},

*c*(

*k*). Although these equations can be solved by traditional methods such as Gaussian elimination or LU factorization, a more efficient digital-filtering technique was advocated by Unser et al. [27]. Owing to their work, the time costs of calculating interpolation coefficients are generally negligible, and thus the computation speeds of traditional interpolation and generalized interpolation are roughly the same, if the orders and supports of the interpolation bases are the same.

For generalized interpolation [Eq. (3)], it can be represented in the form of Eq. (2), so that generalized interpolation belongs to convolution-based interpolation. At this time, *φ*_{int}(*x*) is referred to as *cardinal interpolation basis*. It should be emphasized that *φ*_{int}(*x*) and *φ*(*x*) are completely different for generalized interpolation; *φ*_{int}(*x*) does not have close-form and the support of *φ*_{int}(*x*) is infinite; thus the cardinal interpolation bases are mainly used for theoretical analysis. The relation between *φ*_{int}(*x*) and *φ*(*x*) is

*p*(

*k*) is an infinite sequence whose

*Z*transform

*P*(

*z*) satisfies

The introduction of interpolation coefficients offers new possibilities and thus enhances the interpolation accuracy significantly. Figure 1(a2) depicts the cardinal interpolation bases of several generalized interpolation methods, including cubic BSpline (BSpline3), cubic OMOMS (OMOMS3), quintic BSpline (BSpline5), and septic BSpline (BSpline7). These cardinal interpolation bases have infinite support. Due to they are more analogous to *sinc* function, the interpolation accuracies of methods shown in Fig. 1(b) are higher than methods shown in Fig. 1(a). From the view of frequency domain, the transfer functions shown in Fig. 1(b) clearly present the superiority of generalized interpolation.

Generalized interpolation is greatly superior to traditional interpolation and thus is widely used. In digital image correlation community, because traditional interpolation such as cubic interpolation introduces large bias errors, it is suggested to employ generalized interpolation, for instance, BSpline and OMOMS interpolation, to enhance the interpolation accuracy [20,22].

#### 2.2 Deficiencies of existing methods

In [17], the noise-induced bias was reported and the mathematical expressions of noise-induced bias for linear and cubic interpolation were derived. Nevertheless, mathematical framework provided by [17] needs to present a subpixel intensity using several nearby samples. This procedure can be performed for traditional interpolation whose supports are finite; however, it cannot be performed for generalized interpolation whose supports are infinite. Considering the widely use and high accuracy of generalized interpolation, it is necessary to extend the analysis to generalized interpolation, and this is the aim of this paper.

#### 2.3 Systematic errors of digital image correlation

At first, the 1-dimensional case is considered. Suppose the original function is *f*(*x*); the domain of *f*(*x*) is (-∞,∞). The samples of *f*(*x*) are *f*(*n*); the sequence *f*(*n*) can be regarded as the reference image. *f*(*x*) is translated along the *x* axis by *u*_{0} units, thereby the translated function *f*(*x*-*u*_{0}). The samples of *f*(*x*-*u*_{0}) are *f*(*n*-*u*_{0}); the sequence *f*(*n*-*u*_{0}) can be regarded as the deformed image. Interpolation is performed to reconstruct *f*(*x*-*u*_{0}) by its samples *f*(*n*-*u*_{0}), and the reconstructed function *g*(*x*) is given as [Eq. (2)]

*φ*

_{int}(

*x*) is the basis of convolution-based interpolation.

DIC measures the real displacement *u*_{0} by optimizing specific correlation criterion. If the sum of squared difference (SSD) criterion is employed, the measured displacement *u* is given as

Noise is inevitable in practice; the existence of noise will induce additional errors. Suppose the noise in the reference image is *ε _{f}*(

*n*); then the observed intensity

*f*

_{noise}(

*n*) is

*ε*(

_{g}*n*); since

*g*(

*n*+

*u*) is at a subpixel position, from Eq. (7), the influence of noise and interpolation will couple:

*g*(

*n*+

*u*) is the intensity in the absence of noise, and

*ε*(

_{g}*n + u*) is the intensity caused by noise. Substitute Eq. (9) and Eq. (10) into Eq. (8); then C

_{SSD}(

*u*) becomes

*u*) is a deterministic variable while Λ(

*u*) is a random variable.

The measured displacement *u* satisfies C_{SSD}’(*u*) = 0, and thus

*u*is not equal to the actual displacement

*u*

_{0}, but contains measurement error

*u*. Because

_{e}*u*=

*u*

_{0}+

*u*and

_{e}*u*is generally small, under first order approximation

_{e}*u*is small, the denominator of Eq. (13) is much larger than the numerator. In the denominator, the variation of Λ”(

_{e}*u*

_{0}) is generally much smaller than Γ” (

*u*

_{0}), which corresponds to the sum of squared intensity gradients. Therefore, the denominator of Eq. (13) is approximated by its expectation, and thus the expectation and variation of

*u*are

_{e}*N*points are contaminated by intensity noise with the standard deviation

*σ*. In this case,

*u*

_{e}is given as follows (the details are in the Appendix)

#### 2.4 Noise-induced bias for convolution-based interpolation

Equation (15) is complicated. Moreover, it contains terms about *g*(*x*), which are not convenient to calculate in practice. It is preferable that the estimation equation merely contains terms about *f*(*x*). Therefore, it is essential to simplify Eq. (15). If ${\sum}_{n=-\infty}^{\infty}{{g}^{\prime}}^{2}\left(n+{u}_{0}\right)}\gg N{\sigma}^{2$, the systematic errors can be approximated as

*u*is the interpolation bias; it has been thoroughly studied [17,20,21].

_{ib}*u*is the noise-induced bias; the aim of this paper is to derive mathematical formula of it and provide efficient technique to estimate it. Suppose ${\sum}_{n=-\infty}^{\infty}{{g}^{\prime}}^{2}\left(n+{u}_{0}\right)}\approx {\displaystyle {\sum}_{n=-\infty}^{\infty}{{f}^{\prime}}^{2}\left(n\right)$, namely the sum of squared intensity gradients of the reference and the deformed images are roughly the same (the validity of this assumption is discussed in the Appendix),

_{nb}*u*

_{0}) is referred to as

*interpolation-noise coupling function*in this work. The original contribution of this work is to derive the explicit expression of Φ(

*u*

_{0}) in terms of

*φ*

_{int}(

*x*) for both traditional interpolation and generalized interpolation.

#### 2.5 Sinusoidal approximation of noise-induced bias

Although Eq. (17) is simple, it still has some deficiencies in practice. As cardinal interpolation basis *φ*_{int}(*x*) does not have close-form for generalized interpolation, Φ(*x*) does not have close-form as well. Therefore, in practice, only numerical methods are available to calculate Φ(*x*), which is still quite inconvenient. In order to overcome this deficiency, in this paper, Φ(*x*) is approximated by its Fourier representation. Φ(*x*) is a periodic function with period 1, so that Φ(*x*) can be resolved into Fourier series. The Fourier coefficients of Φ(*x*) are expressed using the interpolation transfer function ${\widehat{\phi}}_{\mathrm{int}}\left(\nu \right)$.

Define $\Theta \left(x\right)=-{\phi}_{\mathrm{int}}\left(x\right){{\phi}^{\prime}}_{\mathrm{int}}\left(x\right)$, so that

*x*) is the

*comb*function [28]. Convolution in the spatial domain corresponds to multiplication in the frequency domain. Therefore, in the frequency domain, Eq. (18) corresponds to

*x*). Multiplication in the spatial domain corresponds to convolution in the frequency domain, so that the Fourier transform of $\Theta \left(\nu \right)$ is

*x*) is

*M*is an integer determined by the desired level of precision.

In practice, the subset size is finite. Although the theoretical analysis in this section assumes that the domain of the original function is (-∞,∞), if the subset is extended by adding zeros, this condition will be satisfied. Since the intensities outside the subset are zeros, the infinite sum ${\Sigma}_{n=-\infty}^{\infty}$ in Eq. (15), Eq. (17), and Eq. (22) can be substituted by ${\Sigma}_{n=-\left(N-1\right)/2}^{\left(N-1\right)/2}$, here *N* denotes the subset size.

## 3. Numerical results

#### 3.1 Systematic errors

To validate the correctness of the proposed theoretical analysis, numerical simulations were performed. In the numerical tests, three 1-dimensional speckle patterns with different speckle sizes were utilized. These speckle patterns are composed of individual Gaussian speckles [29]; the intensity *f*(*x*) at *x* is given as

*K*is the total count of Gaussian speckles,

*I*

_{0}is the peak intensity of the each speckle,

*x*

_{k}is the position of the

*k*th speckle with a rectangular distribution, and

*R*is the speckle radius. In this study,

*I*

_{0}is 1. The positions of the speckles

*x*

_{k}are restricted to the interval [-50, 50], so that the total interval length

*L*= 100. The speckle density is defined as 2

*KR/L*; in this simulation, the speckle density is 65%. The speckle radii are 1.5, 2.0, and 3.0 respectively. Figure 2 shows the speckle patterns utilized in the numerical simulations.

In order to evaluate the systematic errors, for given speckle pattern, 20 translated speckle patterns were generated by translating the speckle positions *x*_{k} 0.05 units each time numerically, corresponding to displacements range from 0.05 pixels to 1 pixel. Exact sample values were utilized to remove the quantization errors. Additive Gaussian white noise with standard deviations range from 0.01 to 0.05 were added. In order to evaluate the systematic errors and random errors, the noise addition repeated for *M* times. The corresponding *mean bias error* and *standard deviation error* are given as [30]

*e*

_{u}is the mean bias error,

*σ*

_{u}is the standard deviation error,

*u*

_{i}is the measured displacement corresponding to the

*i*th noise addition, and

*u*

_{0}is the real displacement. In this study, the repetition

*M*= 100000.

The DIC algorithm was implemented by a self-written program: Gauss-Newton method and SSD criterion were utilized; as the underlying deformation is rigid motion, zero-order shape function was employed; the convergence criterion was that the iterative increment $\left|\Delta u\right|$ was less than 1 × 10^{−7}; in order to explore the influence of interpolation, both cubic and quintic BSpline interpolation algorithms were employed. The subsets were chosen as the interval [-60, 60], so that there were 121 points in the subsets in total.

The theoretical estimations were obtained by Eq. (15). As the reasons pointed out in section 2, the infinite sum ${\Sigma}_{n=-\infty}^{\infty}$ in Eq. (15) was substituted by ${\Sigma}_{n=-\left(N-1\right)/2}^{\left(N-1\right)/2}$; in this simulation, as indicated above, the subset size *N* = 121. The cardinal interpolation basis *φ*_{int}(*x*) does not have close-form, so that *φ*_{int}(*x*) was obtained by numerical method; this method is based on the fact that *φ*_{int}(*x*) is the reconstruction function of *δ*(*x*) (Dirac delta function); the derivative *φ’*_{int}(*x*) and the second derivative *φ”*_{int}(*x*) were calculated in a similar way. *g*(*n* + *u*_{0}), *g*’(*n* + *u*_{0}), and *g*”(*n* + *u*_{0}) were calculated numerically from the reconstructed function *g*(*x*).

The systematic errors by DIC and by theory [Eq. (15)] are illustrated in Fig. 3. Figure 3(a1)-3(a2), 3(b1)-3(b2), and 3(c1)-3(c2) correspond to speckle patterns shown in Fig. 2(a), 2(b), and 2(c) respectively. Figure 3(a1)-3(c1) correspond to the cubic BSpline interpolation while Fig. 3(a2)-3(c2) correspond to the quintic BSpline.

Figure 3 demonstrates that the systematic errors by DIC and by theory show excellent agreement. Figure 3 indicates that (1) the noise-induced bias can be much larger than the interpolation bias (interpolation bias corresponds to the systematic errors in the absence of noise); (2) a speckle pattern with large speckle size introduces less interpolation bias, but is more sensitive to noise; (3) if the intensity noise is significant, the bias of a pattern with large speckle size is probably much larger than that of a pattern with small speckle size; (4) for given speckle pattern, high-order BSpline interpolation method reduces both interpolation bias and noise-induced bias.

#### 3.2 Noise-induced bias

This paper proceeds to validate the correctness of the theoretical results on noise-induced bias. The noise-induced bias by DIC was obtained by subtracting the interpolation bias (the systematics errors in the absence of noise) from the systematic errors in the presence of noise [see Eq. (16)]. The noise-induced bias by theoretical analysis was obtained by Eq. (17).

Figure 4 shows both noise-induced bias by DIC and by theoretical analysis. The DIC results are clearly consistent with our theoretical estimations. Figure 4 indicates that (1) the noise-induced bias increase as noise level increase; (2) the noise-induced bias of fine speckle patterns are less than that of coarse patterns; (3) for given speckle pattern and noise level, high order BSpline interpolation introduces less noise-induced bias.

Observation (1) and (2) have been reported and well explained [11]; observation (3) has been reported in [11], but its physical insight still remains unexplained. In the following section, this phenomenon will be explained.

#### 3.3 Sinusoidal approximation

Last section demonstrates the correctness of our theoretical results on noise-induced bias. Therefore, for the sinusoidal approximation, the issue is how many terms is enough to well represent the interpolation noise coupling function Φ(*x*). To tackle this problem, the coupling function Φ(*x*) corresponding to Keys, cubic BSpline, quintic BSpline, and septic BSpline were calculated numerically. The corresponding Fourier coefficients *a*_{m} were calculated as well. Because the Fourier coefficients merely depend on the interpolation transfer function, for given interpolation method, *a*_{m} are constants. The interpolation transfer functions of aforementioned interpolation methods are [23]

*integral*in MATLAB, the Fourier coefficients were calculated. These coefficients are shown in Table 1. Figure 5 shows the interpolation-noise coupling function Φ(

*x*) and corresponding fundamental frequency approximation

*a*

_{1}sin2π

*x*, indicating that the fundamental frequency approximation is adequate in practice. Therefore, for given noise level and speckle pattern, the noise-induced bias is primarily determined by the coefficient of the fundamental frequency,

*a*

_{1}. Table 1 indicates that,

*a*

_{1}of more sophisticated interpolation methods (high order BSpline are more sophisticated than the low order BSpline; cubic BSpline is more sophisticated than Keys method) are smaller, and therefore the corresponding noise-induced bias are smaller. The deep level reason lies in Eq. (21): since the transfer function of a sophisticated interpolation method is more analogous to the ideal low-pass filter, the overlap between ${\widehat{\phi}}_{\mathrm{int}}\left(\nu \right)$ and ${\widehat{\phi}}_{\mathrm{int}}\left(1-\nu \right)$ is smaller, so that

*a*

_{1}is smaller, and thus the noise-induced bias is smaller. Hence, why sophisticated interpolation generally reduce noise-induced bias is well explained.

## 4. Experimental verification

Actual subpixel translation experiment was performed to verify the correctness of the mathematical derivations of noise-induced bias.

The experiment method has been published in [21]. To clarify the situation, the principle is explained briefly. The experiment set-up is shown in Fig. 6(a). The idea of this method is to display a speckle pattern on the computer screen and then translate the speckle pattern by an image processing software [Fig. 6(b)]; the integer pixel translation on the computer screen will introduce a subpixel translation on the sensor plane [Fig. 6(c)]. The experiment was carried out as follows. (1) A camera (OK_IM1161, Beijing JoinHope Image Technology Ltd.) with 12mm lens was rigidly fixed to a variation-isolation table; then a tablet computer (Microsoft Surface Pro 3) was placed on the table. The distance was adjusted so that 10 pixel on the computer screen roughly corresponded to a pixel of the camera sensor. (2) A laser pen was placed on the camera and its direction was justified until the laser point was at the middle of the image. A CD disk was placed on the computer screen to reflect the laser beam; meanwhile, the position of the computer was adjusted until the reflected laser beam coincided with the incident laser beam, which ensured the perpendicular of the computer screen and the optical axis. Then the CD disk and the laser pen were carefully taken away. (3) A speckle pattern was displayed on the computer screen. The speckle pattern was translated using an image processing software. During the translation, a wireless keyboard was used. For each translation step, 100 frames were captured to analyze the noise distribution.

As before, the noise-induced bias was obtained by subtracting the interpolation bias (bias without noise) from the systematic errors (bias with noise). Therefore, it is necessary to obtain the bias both with and without noise. The data processing flows corresponding to bias with and without noise are different; they are shown in Fig. 7 schematically. The computation of bias without noise is shown in Fig. 7(a): the deformed images were averaged first to obtain a virtual noiseless deform image, and then the average deform image correlated with the reference image to obtain the measured displacement. The computation of bias with noise is shown in Fig. 7(b): the reference image correlated with 100 noisy images, thereby 100 measured displacements corrupted by image noise; these 100 measured displacements were averaged to obtain an average displacement. In both cases, the measured displacements were fitted by a linear model and the differences with the linear fit were the bias. The reference image was the average of 100 images and thus was virtually noiseless.

The DIC algorithm was implemented by a self-written program. Gauss-Newton method was utilized. In order to be consistent with the theoretical analysis, SSD criterion was employed. Zero order shape function was used. The subset size was chosen as 137 × 78. In order to investigate the influence of speckle size, three speckle patterns were utilized in this work as shown in Figs. 8(a1)-8(c1). The densities of these speckle patterns are 65%. The speckle sizes are roughly 4, 6, and 8 pixels respectively.

To estimate the noise-induced bias, the image noise should be known. The standard deviations of the image noise were evaluated by 100 frames and the results are depicted in Figs. 8(a2)-8(c2). In accord with [31], the image noise depends on corresponding image intensity. Hence, it is essential to extend previous theoretical analysis to cases of non-uniform noise. Besides, images are 2-dimensional signals. After some mathematical derivations, the formula of noise-induced bias for non-uniform noise is given as

*M*and

*N*are the height and the width of the subset respectively, $\overline{{\sigma}^{2}}$ is the mean of noise variance, and $\overline{{{f}^{\prime}}_{x}^{2}}$ is the mean of squared intensity gradients.

As indicated in [21], experimental results of BSpline are sensitive to environmental noise; thus Keys interpolation was employed. The interpolation basis of Keys interpolation is [26]

To estimate noise-induced bias by Eq. (26), three terms are required: the interpolation noise coupling function, the mean of noise variance, and the mean of squared intensity gradients. In this experiment, Φ(*u*_{0}) is shown in Eq. (28), the $\overline{{\sigma}^{2}}$ are 3.4053, 3.4175, 3.3956 respectively, and $\overline{{{f}^{\prime}}_{x}^{2}}$ are 1152.9, 983.37, and 809.38 respectively.

Figures 8(a3)-8(b3) illustrate the systematic errors with noise and without noise. The DIC results were obtained following the computational flow chart shown in Fig. 7. The theoretical estimation of the interpolation bias was based on the formula in [21]; the estimation of the noise-induced bias was based on Eq. (26); the estimation of the systematic error was the sum of estimations of interpolation bias and noise-induced bias. Figures 8(a4)-8(c4) show the noise-induced bias by experiment and theory, which show good agreement. In this experiment, the noise-induced bias is serval thousandth pixels, about one-tenths of the interpolation bias. As the speckle density is the same, the intensity gradient of a fine pattern is larger than that of a coarse pattern, resulting in a smaller noise-induced bias.

## 5. Discussion

#### 5.1 Comparison with findings in literatures

Findings in [17] are special cases of theoretical analysis presented in this paper. For example, if linear interpolation is used, the linear interpolation basis is

#### 5.2 Cause of noise-induced bias

In order to obtain insight into noise-induced bias, it is crucial to know its cause. In [11], Sutton et al. noted that the noise-induced bias is due to the position-dependence of noise caused by the low-pass effect of interpolation. This position-dependence means that, even though the variance of noise is${\sigma}_{0}^{2}$ at every integer position, the variance of interpolated intensity *σ*^{2}(*u*) is still a function of subpixel position *u*

*k*and the sub-pixel position

*u*[11]. However, Eq. (31) is seldom used, for $\widehat{h}\left(k,u\right)$ is generally difficult to obtain and the integral is generally difficult to calculate. In addition, this explanation is not intuitionistic. In this section, the position-dependence is characterized in the spatial domain directly, and a more intuitionistic explanation is given.

Suppose function *f*(*x*); its samples *f*(*k*) are corrupted by noise *ε*(*k*). If convolution-based interpolation is used, the interpolated intensity at a subpixel position *u* is

*f*

_{T}(

*u*) is the interpolated intensity in the absence of noise. It is evident that the expectation of

*f*

_{noise}(

*u*) is

*f*

_{T}(

*u*). If the variance of noise is ${\sigma}_{0}^{2}$ at each sample point, the variance of interpolated intensity

*σ*

^{2}(

*u*) depends on the subpixel position

*u*:

To validate the correctness of Eq. (33), numerical simulations were performed for the 1-dimentional speckle pattern shown in Fig. 2(a). Additive Gaussian white noise with standard deviation *σ*_{0} = 0.01 were added at integer pixel positions. In order to simulate the randomness of the noise, the addition repeated for 1000000 times. Keys, BSpline, and OMOMS interpolation were employed to interpolate intensities at subpixel positions in the interval [-1,1]. Figure 9 illustrates the standard deviations of the interpolated intensity *σ*(*u*) by numerical experiments and by Eq. (33); the simulation and theoretical results show excellent agreement. Remarkably, it can be found that (1) the fluctuation of high-order BSpline is more unobvious; (2) as 0 is the non-differentiable point of the interpolation basis of OMOMS [23], *σ*(*u*) of OMOMS is **C**^{0} continuity; (3) for Keys and BSpline interpolation methods, because of continuity and symmetry, *σ*(*u*) achieve their extreme values at integer and half-pixel positions.

On the basis of Eq. (33), the form of the interpolation-noise coupling function Φ(*u*) can be explained qualitatively (recalling that$\Phi \left(u\right)=-{\displaystyle \sum _{k=-\infty}^{\infty}{\phi}_{\mathrm{int}}\left(u-k\right){{\phi}^{\prime}}_{\mathrm{int}}\left(u-k\right)}$). If SSD criterion is employed [see Eq. (8)], the expectation of C_{SSD}(*u*) is

*u*, and it is this term that introduces noise-induce bias. Considering that the measured displacement satisfies C

_{SSD}’(

*u*) = 0, the noise-induced bias should depend on the derivative of ${\sum}_{k=-\infty}^{\infty}{\phi}_{\mathrm{int}}^{2}\left(u-k\right)$, so that the interpolation-noise coupling function should depend on ${\sum}_{k=-\infty}^{\infty}{\phi}_{\mathrm{int}}\left(u-k\right){{\phi}^{\prime}}_{\mathrm{int}}\left(u-k\right)$, which is accord with our previous analysis.

## 6. Conclusion

The purpose of this paper is to provide theoretical formulae to estimate noise-induce bias for both traditional interpolation and generalized interpolation. The conclusions are drawn as follows:

- 1) Considering that existing methods for estimating noise-induced bias are not applicable to the widely used generalized interpolation, this paper presents a theoretical analysis of noise-induced bias for convolution-based interpolation. The validity of the theoretical derivations is established by both numerical simulations and actual subpixel translation experiment. Analysis provided by this work can be applied to both traditional interpolation and generalized interpolation; findings in literatures are included as special cases. Compared to the existing methods, formulae provided by this paper are simpler, briefer, and more general.
- 2) In order to simplify the estimation of noise-induced bias, a sinusoidal approximate formula is derived which is capable of estimating the noise-induced bias with speed, ease, and accuracy. Based on the approximate formula, the mechanism of sophisticated interpolation methods generally reducing noise-induced bias is revealed.
- 3) In order to provide a more intuitionistic explanation of the cause of noise-induced bias, the position-dependence of noise is quantitatively characterized in the spatial domain, by which the cause of noise-induced bias is explained.

The primary motivation of this work is to provide a valid quality assessment of speckle patterns and thus facilitating the standardization of speckle patterns for DIC.

## Appendix

#### A. Detailed derivation of the systematic errors

Since

*N*points are contaminated by noise with standard deviation

*σ,*and hence

*u*) is

#### B. The dependence of sum of squared intensity gradients upon subpixel positions

The dependence of the sum of squared intensity gradients upon subpixel positions is characterized in this appendix. The Fourier transform of function *f*(*x*) is

*g*’(

*x*+

*u*

_{0}) is

*g*’(

*x*+

*u*

_{0}) is

*f*(

*x*) satisfies sampling theorem, and thus

*k*+

*m =*0. In this case, using the Parseval’s theorem, the sum of squared intensity gradients can be expressed in the frequency domain:

In order to verify the validity of Eq. (47), numerical experiments were performed. Cubic (4 tap), quintic (6 tap), and septic (8 tap) BSpline interpolation methods were utilized to obtain the sum of squared intensity gradients for speckle patterns shown in Fig. 2. Figures 10(a1)-10(a3) illustrate the sum of squared intensity gradients obtained by DIC and by Eq. (47) for speckle patterns shown in Figs. 2(a)-2(c). The theoretical estimations show good agreement with the results by DIC.

Equation (47) implies that SSSIG can be divided into two parts: a direct component, and a cosine-shaped component. The amplitude of the cosine-component is determined by the integral of the product of the power spectrum and the kernel *E*_{COS}(ν). *E*_{COS}(ν) related to cubic, quintic, and septic BSpline interpolation methods are illustrated in Fig. 10(b). A high order BSpline interpolation method has small frequency response *E*_{COS}(ν), and therefore the fluctuation is less obvious. Because small speckle size gives rise to significant high frequency component and the absolute value of *E*_{COS}(ν) is large at high frequency, the fluctuations of fine patterns are larger than that of coarse patterns. However, the cosine component is general negligible compared to the direct component, and therefore it’s rational to consider $\sum _{n=-\infty}^{\infty}{{g}^{\prime}}^{2}\left(n+{u}_{0}\right)}\approx {\displaystyle \sum _{n=-\infty}^{\infty}{{{g}^{\prime}}^{2}\left(n+{u}_{0}\right)|}_{{u}_{0}=0}}={\displaystyle \sum _{n=-\infty}^{\infty}{{f}^{\prime}}^{2}\left(n\right)$.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos. 11332010, 51271174, 11372300, 11127201, 11472266 and 11428206).

## References and links

**1. **W. H. Peters and W. F. Ranson, “Digital imaging techniques in experimental stress analysis,” Opt. Eng. **21**(3), 427–431 (1982). [CrossRef]

**2. **T. C. Chu, W. F. Ranson, and M. A. Sutton, “Applications of digital-image-correlation techniques to experimental mechanics,” Exp. Mech. **25**(3), 232–244 (1985). [CrossRef]

**3. **I. Yamaguchi, “A laser-speckle strain gauge,” J. Phys. E Sci. Instrum. **14**(11), 1270–1273 (1981). [CrossRef]

**4. **F. Hild and S. Roux, “Digital image correlation: from displacement measurement to identification of elastic properties – a review,” Strain **42**(2), 69–80 (2006). [CrossRef]

**5. **M. S. Kirugulige, H. V. Tippur, and T. S. Denney, “Measurement of transient deformations using digital image correlation method and high-speed photography: application to dynamic fracture,” Appl. Opt. **46**(22), 5083–5096 (2007). [CrossRef] [PubMed]

**6. **Z. Liu and J. Gao, “Deformation-pattern-based digital speckle correlation for coefficient of thermal expansion evaluation of film,” Opt. Express **19**(18), 17469–17479 (2011). [CrossRef] [PubMed]

**7. **G. F. Xiang, Q. C. Zhang, H. W. Liu, X. P. Wu, and X. Y. Ju, “Time-resolved deformation measurements of the Portevin–Le Chatelier bands,” Scr. Mater. **56**(8), 721–724 (2007). [CrossRef]

**8. **Q. Zhang, Z. Jiang, H. Jiang, Z. Chen, and X. Wu, “On the propagation and pulsation of Portevin-Le Chatelier deformation bands: an experimental study with digital speckle pattern metrology,” Int. J. Plast. **21**(11), 2150–2173 (2005). [CrossRef]

**9. **D. Lecompte, A. Smits, S. Bossuyt, H. Sol, J. Vantomme, D. Van Hemelrijck, and A. M. Habraken, “Quality assessment of speckle patterns for digital image correlation,” Opt. Lasers Eng. **44**(11), 1132–1145 (2006). [CrossRef]

**10. **H. Haddadi and S. Belhabib, “Use of rigid-body motion for the investigation and estimation of the measurement errors related to digital image correlation technique,” Opt. Lasers Eng. **46**(2), 185–196 (2008). [CrossRef]

**11. **M. A. Sutton, J. J. Orteu, and H. Schreier, *Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications* (Springer Science & Business Media, 2009).

**12. **B. Pan, K. Qian, H. Xie, and A. Asundi, “Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review,” Meas. Sci. Technol. **20**(6), 062001 (2009). [CrossRef]

**13. **M. Bornert, F. Bremand, P. Doumalin, J. C. Dupre, M. Fazzini, M. Grediac, F. Hild, S. Mistou, J. Molimard, J. J. Orteu, L. Robert, Y. Surrel, P. Vacher, and B. Wattrisse, “Assessment of digital image correlation measurement errors: methodology and results,” Exp. Mech. **49**(3), 353–370 (2009). [CrossRef]

**14. **M. A. Sutton, J. L. Turner, H. A. Bruck, and T. A. Chae, “Full-field representation of discretely sampled surface deformation for displacement and strain analysis,” Exp. Mech. **31**(2), 168–177 (1991). [CrossRef]

**15. **H. W. Schreier and M. A. Sutton, “Systematic errors in digital image correlation due to undermatched subset shape functions,” Exp. Mech. **42**(3), 303–310 (2002). [CrossRef]

**16. **M. A. Sutton, S. R. McNeill, J. Jang, and M. Babai, “Effects of subpixel image restoration on digital correlation error estimates,” Opt. Eng. **27**(10), 870–877 (1988). [CrossRef]

**17. **Y. Q. Wang, M. A. Sutton, H. A. Bruck, and H. W. Schreier, “Quantitative error assessment in pattern matching: effects of intensity pattern noise, interpolation, strain and image contrast on motion measurements,” Strain **45**(2), 160–178 (2009). [CrossRef]

**18. **X. Xu, Y. Su, Y. Cai, T. Cheng, and Q. Zhang, “Effects of various shape functions and subset size in local deformation measurements using DIC,” Exp. Mech. **55**(8), 1575–1590 (2015). [CrossRef]

**19. **F. Amiot, M. Bornert, P. Doumalin, J. C. Dupre, M. Fazzini, J. J. Orteu, C. Poilane, L. Robert, R. Rotinat, E. Toussaint, B. Wattrisse, and J. S. Wienin, “Assessment of digital image correlation measurement accuracy in the ultimate error regime: main results of a collaborative benchmark,” Strain **49**(6), 483–496 (2013). [CrossRef]

**20. **H. W. Schreier, J. R. Braasch, and M. A. Sutton, “Systematic errors in digital image correlation caused by intensity interpolation,” Opt. Eng. **39**(11), 2915–2921 (2000). [CrossRef]

**21. **Y. Su, Q. Zhang, Z. Gao, X. Xu, and X. Wu, “Fourier-based interpolation bias prediction in digital image correlation,” Opt. Express **23**(15), 19242–19260 (2015). [CrossRef] [PubMed]

**22. **L. Luu, Z. Wang, M. Vo, T. Hoang, and J. Ma, “Accuracy enhancement of digital image correlation with B-spline interpolation,” Opt. Lett. **36**(16), 3070–3072 (2011). [CrossRef] [PubMed]

**23. **P. Thévenaz, T. Blu, and M. Unser, “Interpolation revisited,” IEEE Trans. Med. Imaging **19**(7), 739–758 (2000). [CrossRef] [PubMed]

**24. **E. Meijering, “A chronology of interpolation: from ancient astronomy to modern signal and image processing,” Proc. IEEE **90**(3), 319–342 (2002). [CrossRef]

**25. **M. Unser, “Splines: a perfect fit for signal and image processing,” IEEE Signal Process. Mag. **16**(6), 22–38 (1999). [CrossRef]

**26. **R. G. Keys, “Cubic convolution interpolation for digital image-processing,” IEEE Trans. Acoust. Speech Signal Process. **29**(6), 1153–1160 (1981). [CrossRef]

**27. **M. Unser, A. Aldroubi, and M. Eden, “Fast B-spline transforms for continuous image representation and interpolation,” IEEE Trans. Pattern Anal. Mach. Intell. **13**(3), 277–285 (1991). [CrossRef]

**28. **J. W. Goodman, *Introduction to Fourier Optics* (Roberts and Company Publishers, 2005).

**29. **P. Zhou and K. E. Goodson, “Subpixel displacement and deformation gradient measurement using digital image/speckle correlation (DISC),” Opt. Eng. **40**(8), 1613–1620 (2001). [CrossRef]

**30. **B. Pan, H. Xie, B. Xu, and F. Dai, “Performance of sub-pixel registration algorithms in digital image correlation,” Meas. Sci. Technol. **17**(6), 1615–1621 (2006). [CrossRef]

**31. **X. D. Ke, H. W. Schreier, M. A. Sutton, and Y. Q. Wang, “Error assessment in stereo-based deformation measurements,” Exp. Mech. **51**(4), 423–441 (2011). [CrossRef]