Fourier-based interpolation bias prediction in digital image correlation

Yong Su; Qingchuan Zhang; Zeren Gao; Xiaohai Xu; Xiaoping Wu

doi:10.1364/OE.23.019242

1. Introduction

Digital image correlation has evolved into a reliable and flexible full-field optical metrology [1–3], and this technique has found numerous and ever-increasing applications in scientific and engineering fields [4–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–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.

2. 1-D interpolation bias analysis and numerical verification

2.1 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.

2.2 Principle of convolution-based interpolation

Interpolation is the process of constructing a function $f_{T} (x)$ to approximate the original function $f (x)$ though given sample values $f [n]$ . Nowadays, convolution-based interpolation is frequently employed because of its versatility and relatively low complexity. The essence of convolution-based interpolation is to represent $f (x)$ as a discrete sum of weighted and shifted basis functions ${r (x - k) | k \in ℤ}$ :

f_{T} (x) = \sum_{k = - \infty}^{\infty} c_{k} r (x - k),

where the interpolation coefficient

c_{k}

is chosen so that the original and reconstructed functions have identical values at the sample points. Equation (1) is equivalent to:

f_{T} (t) = \sum_{k = - \infty}^{\infty} f [k] φ (t - k) .

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].

2.3 1-D interpolation bias analysis

It is illuminating to analyze the 1-D situation first. Figure 1(a) illustrates a one-dimensional continuous function $f (x)$ whose domain is $ℝ$ (functions with limited domain can be extended by 0). $f (x)$ works similarly to the reference image in digital image correlation and is therefore called the reference function. $f (x)$ is sampled [see Fig. 1(b)] to produce a reference sequence $f [n] = f (n)$ . Figure 1(c) shows the process of translating $f (x)$ towards the positive direction of $x$ by $u_{0}$ units to obtain the translated function $f (x - u_{0})$ . The translated function is then sampled [see Fig. 1(d)] to obtain the target sequence $g [n] = f (n - n_{0})$ . Suppose that the interpolation basis is $φ (x)$ [Fig. 1(e) illustrates the basis related to linear interpolation]. According to Eq. (2) the reconstructed function is:

g (x) = \sum_{n = - \infty}^{\infty} g [n] φ (x - n) .

Figure 1(f) depicts the target function reconstructed by linear interpolation. If the sum of squared differences criterion is utilized, the calculated result

u

by digital image correlation satisfies:

u = argmin \sum_{n = - \infty}^{\infty} {[f (n) - g (n + u)]}^{2} .

The calculated result

u

is not equal to the real displacement

u_{0}

, the discrepancy between

u_{0}

and

u

is a result of imperfect interpolation [see Figs. 1(c) and 1(f)]. The objective of this section is to express the interpolation bias

u_{e} = u - u_{0}

in terms of

f (x)

φ (x)

.

Fig. 1 Translation simulation: (a) reference function; (b) sampled values of reference function; (c) translated function; (d) sampled values of translated function; (e) interpolation basis function; (f) reconstructed function.

Download Full Size | PDF

Parseval’s theorem applied to Eq. (4) yields:

\begin{array}{l} u = argmin \int_{- 1 / 2}^{1 / 2} {| \hat{F} (ν) - \hat{G} (ν) |}^{2} d ν \\ \hat{F} (ν) = \sum_{n = - \infty}^{\infty} f (n) e^{- j 2 π ν n}, \hat{G} (ν) = \sum_{n = - \infty}^{\infty} g (n + u) e^{- j 2 π ν n} \end{array}

where

\hat{F} (ν)

and

\hat{G} (ν)

are the discrete-time Fourier transforms (DTFTs) of

f (x)

and

g (x)

respectively. The continuous Fourier transform of

f (x)

is denoted as

\hat{f} (ν)

provided that

\hat{f} (ν) = \int_{- \infty}^{\infty} f (x) e^{- j 2 π ν x} d x

. The Poisson summation formulae are employed to present

\hat{F} (ν)

in terms of

\hat{f} (ν)

:

\hat{F} (ν) = \sum_{m = - \infty}^{\infty} \hat{f} (ν - m) .

To express

\hat{G} (ν)

in terms of

\hat{f} (ν)

, the relation between

f (x)

and

\hat{g} (ν)

should be obtained.

g (x)

is generated from

f (x)

through translation, sampling, and interpolation. Spatial translation [see Fig. 1(c)] corresponds to frequency-domain phase translation; spatial sampling [see Fig. 1(d)] corresponds to superposition of shifted replicas of the original spectrum in the frequency domain; and spatial interpolation [see Fig. 1(f)] corresponds to multiplying the original spectrum by the interpolation transfer function

\hat{φ} (ν)

(for the calculation algorithm for

\hat{φ} (ν)

, the reader is referred to [35].). Denote a Fourier transform pair by

\overset{F}{\leftrightarrow}

; then the process described above can be represented symbolically as:

\begin{array}{l} f (x) \overset{F}{\leftrightarrow} \hat{f} (ν) \\ f (x - u_{0}) \overset{F}{\leftrightarrow} e^{- j 2 π ν u_{0}} \hat{f} (ν) \\ f (x - u_{0}) comb (x) \overset{F}{\leftrightarrow} \sum_{k = - \infty}^{\infty} e^{- j 2 π (ν - k) u_{0}} \hat{f} (ν - k) \\ g (x) = φ (x) * [f (x - u_{0}) comb (x)] \overset{F}{\leftrightarrow} \hat{φ} (ν) \sum_{k = - \infty}^{\infty} e^{- j 2 π (ν - k) u_{0}} \hat{f} (ν - k) \end{array}

where

c o m b (x) = \sum_{k = - \infty}^{\infty} δ (x - k)

is the comb function; for a more detailed mathematical derivation of

\hat{g} (ν)

, the reader is referred to [38]. Using the Poisson summation formulae again yields:

\hat{G} (ν) = \sum_{m = - \infty}^{\infty} e^{j 2 π (ν - m) u} \hat{φ} (ν - m) \cdot \sum_{k = - \infty}^{\infty} e^{- j 2 π (ν - k) u_{0}} \hat{f} (ν - k) .

Substituting Eq. (6) and Eq. (8) into Eq. (5) gives:

\begin{array}{l} u = argmin Γ (u) \\ Γ (u) = \int_{- 1 / 2}^{1 / 2} [\begin{array}{l} - \sum_{k = - \infty}^{\infty} \hat{f} (ν - k) \cdot \sum_{n = - \infty}^{\infty} e^{j 2 π (ν - n) u_{0}} {\hat{f}}^{*} (ν - n) \cdot \sum_{m = - \infty}^{\infty} e^{- j 2 π (ν - m) u} {\hat{φ}}^{*} (ν - m) \\ - \sum_{k = - \infty}^{\infty} {\hat{f}}^{*} (ν - k) \cdot \sum_{n = - \infty}^{\infty} e^{- j 2 π (ν - n) u_{0}} \hat{f} (ν - n) \cdot \sum_{m = - \infty}^{\infty} e^{j 2 π (ν - m) u} \hat{φ} (ν - m) \\ + \sum_{m_{1} = - \infty}^{\infty} \sum_{m_{2} = - \infty}^{\infty} e^{j 2 π (m_{2} - m_{1}) u} \hat{φ} (ν - m_{1}) {\hat{φ}}^{*} (ν - m_{2}) \cdot \sum_{n_{1} = - \infty}^{\infty} \sum_{n_{2} = - \infty}^{\infty} e^{j 2 π (n_{1} - n_{2}) u_{0}} \hat{f} (ν - n_{1}) {\hat{f}}^{*} (ν - n_{2}) \end{array}] d ν \end{array}

Equation (9) is the frequency representation of Eq. (4); the first two terms in the integrand of Eq. (9) relate to

- 2 \sum_{n = - \infty}^{\infty} f (n) g (n + u)

, and the last term relates to

\sum_{n = - \infty}^{\infty} g^{2} (n + u)

.

Suppose that $Γ (u)$ is continuous enough; then the calculated result $u$ satisfies:

Γ^{'} (u) = 0.

The calculated result

u = u_{0} + u_{e}

. Recognizing that

u_{e}

is generally small, to a first-order approximation,

Γ^{'} (u_{0}) + {Γ^{'}}^{'} (u_{0}) u_{e} \approx 0.

The explicit dependence of the interpolation bias

u_{e}

upon the reference function

\hat{f} (ν)

and the interpolation transfer function

\hat{φ} (ν)

can be derived as follows:

\begin{array}{l} u_{e} \approx - \frac{Γ^{'} (u_{0})}{{Γ^{'}}^{'} (u_{0})} \\ Γ^{'} (u_{0}) = j 2 π \int_{- 1 / 2}^{1 / 2} [\begin{array}{l} \sum_{k = - \infty}^{\infty} \hat{f} (ν - k) \cdot \sum_{m = - \infty}^{\infty} \sum_{n = - \infty}^{\infty} e^{j 2 π (m - n) u_{0}} (ν - m) {\hat{f}}^{*} (ν - n) {\hat{φ}}^{*} (ν - m) \\ - \sum_{k = - \infty}^{\infty} {\hat{f}}^{*} (ν - k) \cdot \sum_{m = - \infty}^{\infty} \sum_{n = - \infty}^{\infty} e^{- j 2 π (m - n) u_{0}} (ν - m) \hat{f} (ν - n) \hat{φ} (ν - m) \\ + \sum_{n_{1} = - \infty}^{\infty} \sum_{n_{2} = - \infty}^{\infty} \sum_{m_{1} = - \infty}^{\infty} \sum_{m_{2} = - \infty}^{\infty} e^{j 2 π (n_{1} - n_{2} + m_{2} - m_{1}) u_{0}} (m_{2} - m_{1}) \hat{f} (ν - n_{1}) {\hat{f}}^{*} (ν - n_{2}) \hat{φ} (ν - m_{1}) {\hat{φ}}^{*} (ν - m_{2}) \end{array}] d ν \\ {Γ^{'}}^{'} (u_{0}) = 4 π^{2} \int_{- 1 / 2}^{1 / 2} [\begin{array}{l} \sum_{k = - \infty}^{\infty} \hat{f} (ν - k) \cdot \sum_{m = - \infty}^{\infty} \sum_{n = - \infty}^{\infty} e^{j 2 π (m - n) u_{0}} {(ν - m)}^{2} {\hat{f}}^{*} (ν - n) {\hat{φ}}^{*} (ν - m) \\ + \sum_{k = - \infty}^{\infty} {\hat{f}}^{*} (ν - k) \cdot \sum_{m = - \infty}^{\infty} \sum_{n = - \infty}^{\infty} e^{- j 2 π (m - n) u_{0}} {(ν - m)}^{2} \hat{f} (ν - n) \hat{φ} (ν - m) \\ - \sum_{n_{1} = - \infty}^{\infty} \sum_{n_{2} = - \infty}^{\infty} \sum_{m_{1} = - \infty}^{\infty} \sum_{m_{2} = - \infty}^{\infty} e^{j 2 π (n_{1} - n_{2} + m_{2} - m_{1}) u_{0}} {(m_{2} - m_{1})}^{2} \hat{f} (ν - n_{1}) {\hat{f}}^{*} (ν - n_{2}) \hat{φ} (ν - m_{1}) {\hat{φ}}^{*} (ν - m_{2}) \end{array}] d ν \end{array}

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 ${Γ^{'}}^{'} (u_{0})$ is a periodic function with period 1 and contains a direct component, and recognizing that $\hat{φ} (ν)$ is an approximation of an ideal low-pass filter [see Fig. 2(b)], ${Γ^{'}}^{'} (u_{0})$ can be approximated by (let $k = m = n = 0$ ):

{Γ^{'}}^{'} (u_{0}) \approx 4 π^{2} \int_{- 1 / 2}^{1 / 2} ν^{2} (\hat{φ} (ν) + {\hat{φ}}^{*} (ν)) {| \hat{f} (ν) |}^{2} d ν .

Γ^{'} (u_{0})

is a periodic function with period 1 analogous to

{Γ^{'}}^{'} (u_{0})

and therefore can be expanded into a Fourier series,

Γ^{'} (u_{0}) = \sum_{m = - \infty}^{\infty} a_{m} e^{j 2 π m u_{0}}

. 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 functions satisfy

\hat{φ} (ν) = {\hat{φ}}^{*} (ν)

. Utilizing the band-limited hypothesis and the fact that

\hat{φ} (ν)

is real in practice, Eq. (12) can be simplified to (let

k = n = 0

):

\begin{array}{l} u_{e} \approx - \frac{Γ^{'} (u_{0})}{{Γ^{'}}^{'} (u_{0})} \\ Γ^{'} (u_{0}) = - 4 π \int_{- 1 / 2}^{1 / 2} [\begin{array}{l} \sum_{m = - \infty}^{\infty} (ν - m) \hat{φ} (ν - m) \sin (2 π m u_{0}) \\ + \frac{1}{2} \sum_{m_{1} = - \infty}^{\infty} \sum_{m_{2} = - \infty}^{\infty} (m_{2} - m_{1}) \hat{φ} (ν - m_{1}) \hat{φ} (ν - m_{2}) \sin (2 π (m_{2} - m_{1}) u_{0}) \end{array}] {| \hat{f} (ν) |}^{2} d ν \\ {Γ^{'}}^{'} (u_{0}) \approx 8 π^{2} \int_{- 1 / 2}^{1 / 2} ν^{2} \hat{φ} (ν) {| \hat{f} (ν) |}^{2} d ν \end{array}

Equation (14) is referred to as the band-limited approximation of the interpolation bias.

{Γ^{'}}^{'} (u_{0})

associates with the squared sum of the gray gradient because by Parseval’s theorem

\int f^{'} {(x)}^{2} d x = 4 π^{2} \int ν^{2} {| \hat{f} (ν) |}^{2} d ν

. The integrand of

Γ^{'} (u_{0})

incorporates two parts: the first

\sum_{m = - \infty}^{\infty} (ν - m) \hat{φ} (ν - m) \sin (2 π m u_{0}) + \frac{1}{2} \sum_{m_{1} = - \infty}^{\infty} \sum_{m_{2} = - \infty}^{\infty} (m_{2} - m_{1}) \hat{φ} (ν - m_{1}) \hat{φ} (ν - m_{2}) \sin (2 π (m_{2} - m_{1}) u_{0})

, is determined by the interpolation method exclusively and is a superposition of sinusoidal waves. The second part

{| \hat{f} (ν) |}^{2}

is the power spectrum of

f (x)

. Recalling that the interpolation transfer function

\hat{φ} (ν)

is an approximation of an ideal low-pass filter and that the integration interval is [-0.5, 0.5], terms other than

\hat{φ} (ν)

and

\hat{φ} (ν \pm 1)

can be considered negligible. Moreover, if only the fundamental frequency is involved, further simplification can be accomplished:

\begin{array}{l} u_{e} \approx C \sin 2 π u_{0} \\ C = \frac{1}{2 π} \frac{\int_{- 1 / 2}^{1 / 2} [(ν - 1) \hat{φ} (ν - 1) - (ν + 1) \hat{φ} (ν + 1) + \hat{φ} (ν) \hat{φ} (ν + 1) + \hat{φ} (ν) \hat{φ} (ν - 1)] {| \hat{f} (ν) |}^{2} d ν}{\int_{- 1 / 2}^{1 / 2} ν^{2} \hat{φ} (ν) {| \hat{f} (ν) |}^{2} d ν} \end{array}

where

C

is a constant determined by the reference function

\hat{f} (ν)

and the interpolation algorithm

\hat{φ} (ν)

exclusively. Equation (15) will be referred to hereafter as the sinusoidal approximation of the interpolation bias.

Fig. 2 (a) Interpolation basis function; (b) interpolation transfer function; (c) interpolation bias kernels; (d) aliasing effect of Keys interpolation.

Download Full Size | PDF

Equation (15) explains the well-known sinusoidal-shaped curves of the interpolation bias $u_{e}$ with respect to the prescribed subpixel displacement $u_{0}$ and provides an intuitive physical interpretation of the interpolation bias. Within Eq. (15), ${| \hat{f} (ν) |}^{2}$ is the power spectrum of the reference function, and the factor

E_{i b} (ν) = (ν - 1) \hat{φ} (ν - 1) - (ν + 1) \hat{φ} (ν + 1) + \hat{φ} (ν) \hat{φ} (ν + 1) + \hat{φ} (ν) \hat{φ} (ν - 1)

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(b)], and the interpolation bias kernel [see Fig. 2(c)] of the Keys, cubic BSpline (4-tap), cubic OMOMS, quintic BSpline (6-tap), and septic BSpline (8-tap) interpolation algorithms. The domain of

E_{i b} (ν)

is limited to (−0.5, 0.5) due to the band-limited hypothesis (the frequency 0.5 corresponds to a period of two pixels). As shown in Fig. 2(c), in the low-frequency region

E_{i b} (ν) \approx 0

,

E_{i b} (ν)

increases with increasing

| ν |

until it reaches its peak value at a specific frequency, after which

E_{i b} (ν)

drops again. This tendency quantitatively interprets the result reported in other studies that high-frequency components are the major source of interpolation bias [20].

Because the energy in images is generally concentrated in the low-frequency region, an interpolation algorithm with $E_{i b} (ν)$ smaller in the low-frequency region tend to perform better. For instance, let us compare the Keys, cubic BSpline, and cubic OMOMS interpolation algorithms. Because

{\begin{cases} E_{ib}^{Keys} (ν) > E_{ib}^{BSpline} (ν), | ν | < 0.39 \\ E_{ib}^{B S p l i n e} (ν) > E_{ib}^{O M O M S} (ν), | ν | < 0.35 \end{cases}

so interpolation accuracy generally follows the order OMOMS > BSpline > Keys. Higher-order 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: $\hat{φ} (ν) \neq 1$ when $| ν | < 0.5$ relates to blurring effects; $\hat{φ} (ν) \neq 0$ when $| ν | > 0.5$ relates to aliasing effects. Because pre-filtering the reference function with a low-pass filter can suppress interpolation bias, it is reasonable to attribute interpolation bias to aliasing effects of interpolation, where the high-frequency components exceeding the Nyquist frequency are caused by sampling. From a mathematical perspective, this can be derived from Eq. (15). Figure 2(d) shows the Keys interpolation transfer function $\hat{φ} (ν)$ and its corresponding $\hat{φ} (ν \pm 1)$ . Because the aliasing effect, $\hat{φ} (ν \pm 1) \neq 0$ in the interval [-0.5, 0.5], thereby $E_{i b} (ν) \neq 0$ .

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.

2.4 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

f (x) = \sum_{k = 0}^{N} I_{k} \exp [- \frac{{(x - x_{k})}^{2}}{r^{2}}],

where

N

is the total count of Gaussian speckles,

x_{k}

is the center coordinate of speckle k,

r

is the Gaussian speckle radius, and

I_{k}

is the intensity level of speckle k. Its corresponding spectrum is:

\hat{f} (ν) = \sqrt{π} r e^{- π^{2} r^{2} ν^{2}} \sum_{k = 0}^{N} I_{k} e^{- j 2 π ν x_{k}} .

The spectrum consists of two parts:

\sqrt{π} r e^{- π^{2} r^{2} ν^{2}}

is the Fourier transform of the Gaussian speckle; the other part is determined by the speckle position

x_{k}

.

The intensity level was fixed as $I_{k} = 1$ ; 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⁻¹⁰. 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⁻⁶.

Fig. 3 1-D Synthetic speckle and interpolation bias: (a1)(b1)(c1) speckle patterns of radius 1.5, 2.0, 3.0 respectively; (a2)(b2)(c2) digital image correlation error, full estimate of interpolation bias [Eq. (12)] and sinusoidal approximation of interpolation bias [Eq. (15)] for radii 1.5, 2.0, 3.0.

Download Full Size | PDF

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 $u_{e}$ . For these representative speckle patterns, the interpolation bias follows the order Keys >> BSpline > OMOMS. These phenomena can be explained by the interpolation bias kernel $E_{i b} (ν)$ . 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 ${| \hat{f} (ν) |}^{2}$ and the interpolation bias kernel $E_{i b} (ν)$ . Figure 4 shows the power spectrum (on a logarithmic scale) for speckle radii 1.5, 2.0, and 3.0 and the interpolation bias kernels of Keys, cubic BSpline, and cubic OMOMS; only values in the interval [0,0.5] are displayed due to symmetry. As speckle radius decreases, the spatial period decreases, resulting in an expansion in the frequency domain. Because $E_{i b} (ν)$ is larger in the high-frequency region, $u_{e}$ increases. As illustrated in Fig. 4, the energy of a speckle pattern with radius 3.0 is concentrated in the frequency region $ν < 0.25$ where $E_{i b} (ν) \approx 0$ , and hence $u_{e}$ is about 1 × 10⁻⁴ pixel. The speckle pattern with radius 2.0 possesses more high-frequency energy, and hence $u_{e}$ increases to 0.001 for BSpline and OMOMS and 0.15 for Keys. The speckle pattern with radius 1.0 has a much stronger high-frequency component, leading to a further increase in $u_{e}$ . Although $E_{i b}^{O M O M S} (ν) > E_{i b}^{B S p l i n e} (ν)$ when $ν > 0.35$ , 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.

Fig. 4 Power spectrum (on a logarithmic scale) of speckle patterns with radii 1.5, 2.0, and 3.0 and the interpolation bias kernels of Keys, cubic BSpline and cubic OMOMS interpolation.

Download Full Size | PDF

This 1-D simulation has established the validity of the theoretical results presented earlier.

3. Interpolation bias prediction in digital image correlation

3.1 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 $N$ -dimensional space, suppose that the reference function is $f (x)$ , its corresponding Fourier transform is $\hat{f} (ν)$ , the real displacement is $u_{0}$ , the interpolation bias is $u_{e}$ , the interpolation basis is $φ (x)$ , its interpolation transfer function is $\hat{φ} (ν)$ . Taking an approach similar to the 1-D case and substituting variations for derivatives, the interpolation bias $u_{e}$ satisfies:

\begin{array}{l} \int_{Ω} [\begin{array}{l} \sum_{k} \sum_{m} \sum_{n} e^{j 2 π (m - n) \cdot u_{0}} {\hat{φ}}^{*} (ν - m) \hat{f} (ν - k) {\hat{f}}^{*} (ν - n) (ν - m) \\ - \sum_{k} \sum_{m} \sum_{n} e^{- j 2 π (m - n) \cdot u_{0}} \hat{φ} (ν - m) {\hat{f}}^{*} (ν - k) \hat{f} (ν - n) (ν - m) \\ + \sum_{n_{1}} \sum_{n_{2}} \sum_{m_{1}} \sum_{m_{2}} e^{j 2 π (m_{2} - m_{1} + n_{1} - n_{2}) \cdot u_{0}} \hat{φ} (ν - m_{1}) {\hat{φ}}^{*} (ν - m_{2}) \hat{f} (ν - n_{1}) {\hat{f}}^{*} (ν - n_{2}) (m_{2} - m_{1}) \end{array}] d ν \\ = j 2 π {\int_{Ω} [\begin{array}{l} \sum_{k} \sum_{m} \sum_{n} e^{j 2 π (m - n) \cdot u_{0}} {\hat{φ}}^{*} (ν - m) \hat{f} (ν - k) {\hat{f}}^{*} (ν - n) (ν - m) \otimes (ν - m) \\ + \sum_{k} \sum_{m} \sum_{n} e^{- j 2 π (m - n) \cdot u_{0}} \hat{φ} (ν - m) {\hat{f}}^{*} (ν - k) \hat{f} (ν - n) (ν - m) \otimes (ν - m) \\ - \sum_{n_{1}} \sum_{n_{2}} \sum_{m_{1}} \sum_{m_{2}} e^{j 2 π (m_{2} - m_{1} + n_{1} - n_{2}) \cdot u_{0}} \hat{φ} (ν - m_{1}) {\hat{φ}}^{*} (ν - m_{2}) \hat{f} (ν - n_{1}) {\hat{f}}^{*} (ν - n_{2}) (m_{2} - m_{1}) \otimes (m_{2} - m_{1}) \end{array}] d ν} \cdot u_{e} \end{array}

where

\otimes

denotes dyadic,

\sum_{k} f_{k}

denotes summation in the whole space

ℤ^{N}

, the integral domain is

Ω = {[- 1 / 2, 1 / 2]}^{N}

. High-dimensional interpolation basis can be constructed by tensor-product basis function [23]. Equation (20) is equivalent to a linear system and hence can be readily solved.

Digital image correlation is a particular case where N equals 2. If the translation is along the x-axis, then it is rational to impose the condition that the interpolation bias along the y-direction is 0. Introducing the band-limited hypothesis and taking into account the fact that the transfer function is real, Eq. (20) can be simplified to a band-limited form analogous to Eq. (14):

u_{e} \approx \frac{1}{2 π} \frac{\int_{- 1 / 2}^{1 / 2} \int_{- 1 / 2}^{1 / 2} [\begin{array}{l} \sum_{k = - \infty}^{\infty} \sum_{l = - \infty}^{\infty} (ν_{x} - k) \hat{φ} (ν_{x} - k, ν_{y} - l) \sin (2 π k u_{0}) + \\ \frac{1}{2} \sum_{k_{1} = - \infty}^{\infty} \sum_{k_{2} = - \infty}^{\infty} \sum_{l_{1} = - \infty}^{\infty} \sum_{l_{2} = - \infty}^{\infty} (k_{2} - k_{1}) \hat{φ} (ν_{x} - k_{1}, ν_{y} - l_{1}) \hat{φ} (ν_{x} - k_{2}, ν_{y} - l_{2}) \sin (2 π (k_{2} - k_{1}) u_{0}) \end{array}] {| \hat{f} (ν_{x}, ν_{y}) |}^{2} d ν_{x} d ν_{y}}{\int_{- 1 / 2}^{1 / 2} \int_{- 1 / 2}^{1 / 2} ν_{x}^{2} \hat{φ} (ν_{x}, ν_{y}) {| \hat{f} (ν_{x}, ν_{y}) |}^{2} d ν_{x} d ν_{y}}

If the fundamental frequency is considered exclusively:

\begin{array}{l} u_{e} \approx C \sin (2 π u_{0}) \\ C = \frac{1}{2 π} \frac{\int_{- 1 / 2}^{1 / 2} \int_{- 1 / 2}^{1 / 2} E_{i b} (ν_{x}, ν_{y}) {| \hat{f} (ν_{x}, ν_{y}) |}^{2} d ν_{x} d ν_{y}}{\int_{- 1 / 2}^{1 / 2} \int_{- 1 / 2}^{1 / 2} ν_{x}^{2} \hat{φ} (ν_{x}, ν_{y}) {| \hat{f} (ν_{x}, ν_{y}) |}^{2} d ν_{x} d ν_{y}} \\ E_{i b} (ν_{x}, ν_{y}) = (ν_{x} - 1) \hat{φ} (ν_{x} - 1, ν_{y}) - (ν_{x} + 1) \hat{φ} (ν_{x} + 1, ν_{y}) + \hat{φ} (ν_{x}, ν_{y}) \hat{φ} (ν_{x} - 1, ν_{y}) + \hat{φ} (ν_{x}, ν_{y}) \hat{φ} (ν_{x} + 1, ν_{y}) \end{array}

Equation (22) is similar to Eq. (15), where

E_{i b} (ν_{x}, ν_{y})

is the interpolation bias kernel in two dimensions and the aliasing effect along the y-direction is neglected.

3.2 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 1, Ref. [42].).

Fig. 5 Traditional spray-painted speckle patterns (a) autocorrelation functions of different subsets; (b) experiment speckle image and correlation subsets; (c) interpolation bias by digital image correlation and theoretical prediction by sinusoidal approximation [Eq. (22)] of different subsets.

Download Full Size | PDF

Fig. 6 Speckle patterns of different subsets and the corresponding discrete Fourier transforms; the speckle patterns are magnified to identical size for clarification: (a1)(b1)(c1)(d1) speckle patterns of subsets 51, 101, 251, 501; (a2)(b2)(c2)(d2) discrete Fourier transforms of subsets 51, 101, 251, 501.

Download Full Size | PDF

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 Numerically designed speckle patterns: (a1)(b1) speckle patterns and calculation subsets of coarse and fine patterns; (a2)(b2) discrete Fourier transforms of subsets of coarse and fine patterns; (a3) results of digital image correlation, band-limited approximation [Eq. (21)], and sinusoidal approximation [Eq. (22)] for coarse and fine patterns.

Download Full Size | PDF

Fig. 8 (a1)(b1) Close-ups of coarse and fine patterns; (a2)(b2) autocorrelation of coarse and fine patterns.

Download Full Size | PDF

Figures 7(a3)-7(b3) illustrates the results of digital image correlation and theoretical predictions. The theoretical predictions consist of band-limited and sinusoidal approximations. The digital image correlation results and the theoretical predictions show good agreement. As mentioned above, fine patterns have much more high-frequency energy than coarse patterns [see Figs. 7(a2)-7(b2)] and therefore have a larger interpolation bias. For example, the Keys interpolation bias is roughly 0.01 pixel for coarse patterns [see Fig. 7(a3)] and 0.02 pixel for fine patterns [see 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.

3.3 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.

Fig. 9 (a) Speckle pattern on the computer screen; (b) experimental setup.

Download Full Size | PDF

Fig. 10 (a1)(b1)(c1) Correlation subsets of fine, medium and coarse patterns (see Visualization 1, Visualization 2, and Visualization 3); (a2)(b2)(c2) translations calculated by digital image correlation; (a3)(b3)(c3) interpolation bias by digital image correlation and theoretical prediction by sinusoidal approximation [Eq. (22)]

Download Full Size | PDF

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.

4. 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.

Acknowledgments

This work was supported by the National Basic Research Program of China (2011CB302105) and the National Natural Science Foundation of China (Grant Nos. 11332010, 11472266, 11428206, and 11372300).

References and links

1. 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).

2. 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]

3. 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]

4. J. Goyens, J. Soons, P. Aerts, and J. Dirckx, “Finite-element modelling reveals force modulation of jaw adductors in stag beetles,” J. R. Soc. Interface 11(101), 20140908 (2014). [CrossRef] [PubMed]

5. R. Ghorbani, F. Matta, and M. A. Sutton, “Full-field deformation measurement and crack mapping on confined masonry walls using digital image correlation,” Exp. Mech. 55(1), 227–243 (2015). [CrossRef]

6. H. Zhang, D. Fu, H. Song, Y. Kang, G. Huang, G. Qi, and J. Li, “Damage and fracture investigation of three-point bending notched sandstone beams by DIC and AE techniques,” Rock Mech. Rock Eng. 48(3), 1297–1303 (2014). [CrossRef]

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. Z. Jiang, Q. Zhang, H. Jiang, Z. Chen, and X. Wu, “Spatial characteristics of the Portevin-Le Chatelier deformation bands in Al-4 at%Cu polycrystals,” Mater. Sci. Eng. A-Struct, Mater. Prop. Microstruct. Process. 403(1–2), 154–164 (2005). [CrossRef]

10. M. Pitter, C. See, J. Goh, and M. Somekh, “Focus errors and their correction in microscopic deformation analysis using correlation,” Opt. Express 10(23), 1361–1367 (2002). [CrossRef] [PubMed]

11. H. Yan and B. Pan, “Three-dimensional displacement measurement based on the combination of digital holography and digital image correlation,” Opt. Lett. 39(17), 5166–5169 (2014). [CrossRef] [PubMed]

12. M. Vo, Z. Y. Wang, B. Pan, and T. Y. Pan, “Hyper-accurate flexible calibration technique for fringe-projection-based three-dimensional imaging,” Opt. Express 20(15), 16926–16941 (2012). [CrossRef]

13. H. A. Bruck, S. R. Mcneill, M. A. Sutton, and W. H. Peters III, “Digital image correlation using Newton-Raphson method of partial differential correction,” Exp. Mech. 29(3), 261–267 (1989). [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. M. Bornert, F. Brémand, P. Doumalin, J. C. Dupré, M. Fazzini, M. Grédiac, 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]

16. B. Pan, H. Xie, Z. Wang, K. Qian, and Z. Wang, “Study on subset size selection in digital image correlation for speckle patterns,” Opt. Express 16(10), 7037–7048 (2008). [CrossRef] [PubMed]

17. 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]

18. W. Tong, “An evaluation of digital image correlation criteria for strain mapping applications,” Strain 41(4), 167–175 (2005). [CrossRef]

19. 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]

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. M. A. Sutton, S. R. Mcneill, J. S. Jang, and M. Babai, “Effects of subpixel image restoration on digital correlation error estimates,” Opt. Eng. 27(10), 870–877 (1988). [CrossRef]

22. S. Choi and S. P. Shah, “Measurement of deformations on concrete subjected to compression using image correlation,” Exp. Mech. 37(3), 307–313 (1997). [CrossRef]

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

24. G. K. Rohde, A. Aldroubi, and D. M. Healy Jr., “Interpolation artifacts in sub-pixel image registration,” IEEE Trans. Image Process. 18(2), 333–345 (2009). [CrossRef] [PubMed]

25. 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]

26. P. Mazzoleni, F. Matta, E. Zappa, M. A. Sutton, and A. Cigada, “Gaussian pre-filtering for uncertainty minimization in digital image correlation using numerically-designed speckle patterns,” Opt. Lasers Eng. 66, 19–33 (2015). [CrossRef]

27. B. Pan, “Bias error reduction of digital image correlation using Gaussian pre-filtering,” Opt. Lasers Eng. 51(10), 1161–1167 (2013). [CrossRef]

28. 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]

29. 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]

30. 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]

31. P. L. Reu, “Experimental and numerical methods for exact subpixel shifting,” Exp. Mech. 51(4), 443–452 (2011). [CrossRef]

32. M. Shimizu and M. Okutomi, “Sub-pixel estimation error cancellation on area-based matching,” Int. J. Comput. Vis. 63(3), 207–224 (2005). [CrossRef]

33. J. Inglada, V. Muron, D. Pichard, and T. Feuvrier, “Analysis of artifacts in subpixel remote sensing image registration,” IEEE Trans. Geosci. Rem. Sens. 45(1), 254–264 (2007). [CrossRef]

34. L. Svilainis, K. Lukoseviciute, V. Dumbrava, and A. Chaziachmetovas, “Subsample interpolation bias error in time of flight estimation by direct correlation in digital domain,” Measurement 46(10), 3950–3958 (2013). [CrossRef]

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

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

37. T. Blu, P. Thévenaz, and M. Unser, “MOMS: maximal-order interpolation of minimal support,” IEEE Trans. Image Process. 10(7), 1069–1080 (2001). [CrossRef] [PubMed]

38. S. K. Park and R. A. Schowengerdt, “Image sampling, reconstruction, and the effect of sample-scene phasing,” Appl. Opt. 21(17), 3142–3151 (1982). [CrossRef] [PubMed]

39. 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]

40. Y. Gao, T. Cheng, Y. Su, X. Xu, Y. Zhang, and Q. Zhang, “High-efficiency and high-accuracy digital image correlation for three-dimensional measurement,” Opt. Lasers Eng. 65, 73–80 (2015). [CrossRef]

41. L. Yu and B. Pan, “The errors in digital image correlation due to overmatched shape functions,” Meas. Sci. Technol. 26(4), 045202 (2015). [CrossRef]

42. Y. Su, “Interpolation bias prediction in digital image correlation,” figshare (2015) [retrieved 10 July 2015] http://dx.doi.org/ [CrossRef] .

Name	Description
Code 1	MATLAB implementation of the interpolation bias prediction algorithm
Visualization 1: MOV (7503 KB)	experiment of fine speckle pattern
Visualization 2: MOV (7503 KB)	experiment of medium speckle pattern
Visualization 3: MOV (7503 KB)	experiment of coarse speckle pattern

Fourier-based interpolation bias prediction in digital image correlation

Abstract

1. Introduction

2. 1-D interpolation bias analysis and numerical verification

2.1 Principle of digital image correlation

2.2 Principle of convolution-based interpolation

2.3 1-D interpolation bias analysis

2.4 1-D Simulation

3. Interpolation bias prediction in digital image correlation

3.1 Interpolation bias prediction in high-dimensional space

3.2 Numerical verification

3.3 Subpixel translation experiment

4. Conclusions

Acknowledgments

References and links

Supplementary Material (4)

Cited By

Figures (10)

Equations (22)

Optics Express