Theoretical analysis on performance of digital speckle pattern: uniqueness, accuracy, precision, and spatial resolution

Yong Su; Zeren Gao; Zheng Fang; Yang Liu; Yaru Wang; Qingchuan Zhang; Shangquan Wu

doi:10.1364/OE.27.022439

1. Introduction

Digital image correlation (DIC) is a non-contacting and non-interferometric optical metrology [1, 2]. Because of the merits such as high precision, good versatility, and large measuring range, the DIC technique has been widely used in a diverse set of disciplines [3–8].

Metrological performance is of crucial importance for a measurement method. Grédiac et al. recommended characterizing the metrological performance by three parameters: the accuracy, the precision, and the spatial resolution [9]. Accuracy is a description of systematic error (or bias), precision is a description of random error (or displacement resolution), and spatial resolution is a reflection of the ability of a given technique to distinguish close features [10]. In addition to these three parameters, the subset uniqueness should be taken into account, because reduced uniqueness may give rise to mismatch, leading to erroneous measurements [11, 12].

Metrological performance of DIC depends highly on the quality of speckle pattern [13]. With the emergence of various pattern fabrication techniques, sought of optimized speckle pattern has come under the spotlight. Using genetic algorithm, Bomarito et al. optimized patterns that are composed of black and white pixels, with the consideration of both measurement precision and subset uniqueness [14]. Nevertheless, these pixel-composed patterns are less popular than digital speckle pattern (or numerically designed speckle pattern), which has advantages such as simpleness, controllable randomness, and convenient fabrication [15, 16]. From a mathematical standpoint, digital speckle pattern is a stochastic process controlled by generation parameters such as speckle profile, speckle coverage, and speckle variation. Pattern optimization can be achieved by optimizing these generation parameters. Recently, Chen et al. reported an optimized digital speckle pattern by calculating root-mean-square errors corresponding to various speckle radii and different speckle variations [16]. Lavatelli et al. developed a close-loop pattern optimization engine based on 3D experiment simulation [17]. These works, however, largely depend on numerical techniques, and thus have deficiencies such as long calculation time, lack of generality, and inability to reveal deep physical insight. Hence, there is a great need to develop a theoretical framework for pattern optimization.

Theoretical analysis on pattern optimization is seldom reported. Su et al. built a theoretical model to optimize speckle pattern being a filtered Poisson process (FPP) [18]. Nevertheless, that work has three drawbacks: (1) compared with digital speckle pattern, FPP-typed pattern is less practical, because of its long generation time; (2) the importance of subset uniqueness and spatial resolution was not realized, and therefore was overlooked; (3) optimal coverage could not be identified, because the quantization process is omitted.

The aim of this work is to offer a guide for selection of generation parameters of digital speckle pattern. To achieve this goal, the performance of digital speckle pattern is analyzed theoretically concerning four critical factors: uniqueness, accuracy, precision, and spatial resolution. This paper is organized as follows. In Section 2, digital speckle pattern and digital speckle image are introduced, and their differences are clarified. In Section 3, autocorrelation functions are researched, and formulas are derived to estimate the secondary autocorrelation peak height, which characterizes subset uniqueness. In Section 4, systematic error of 1-dimensional digital speckle pattern is analyzed; specifically, two approximate formulas are proposed. In Section 5, two-dimensional systematic error and random error are studied theoretically, and the correctness of proposed formulas are verified by numerical tests. In Section 6, empirical formulas are presented to feature the spatial resolution; moreover, a rudimentary model for spatial resolution is built. In Section 7, suggestions for selecting generation parameters of digital speckle pattern are put forward, and the limitations of this work are discussed. In Section 8, the conclusions of this work are drawn. To provide more implementation details, we make all raw data and source codes open source [19].

2. Digital speckle pattern and digital speckle image

2.1. 2D digital speckle pattern and 2D digital speckle image

Consider a grid of speckles: the speckle profile is $ψ (x, y)$ and the on-center spacing is a. Then, as shown in Fig. 1, the centroid of each speckle is displaced by imposing an independent random displacement $(ϵ, τ)$ . In this way, a random pattern $h (x, y)$ , termed digital speckle pattern (or numerically-designed speckle pattern), is produced [15, 16]. Essentially, digital speckle pattern is controlled by three factors: the grid spacing a, the speckle profile $ψ (x, y)$ , and the probability distribution of random displacement $(ϵ, τ)$ . These factors have received widespread attention, because they largely determine the measurement performance. By optimizing these pattern generation parameters, the metrological performance could be significantly improved.

Fig. 1 A two-dimensional digital speckle pattern is produced by imposing independent random offsets to a grid of speckles. Here, each disk represents a speckle.

Download Full Size | PDF

From a view of mathematics, the digital speckle pattern can be modeled as

h (x, y) = \sum_{i = - N}^{N} \sum_{j = - N}^{N} ψ (x - x_{i, j}, y - y_{i, j}),

where

x_{i, j} = a i + ϵ_{i, j}

and

y_{i, j} = a j + τ_{i, j}

represent the location of the speckle at the ith column and the jth row. In this paper, random variables, such as

x_{i, j}

and

y_{i, j}

, are written in boldface letters, following the convention in [20]. The total number of speckles is

{(2 N + 1)}^{2}

.

Digital speckle image $f [m, n]$ is produced by sampling and quantizing digital speckle pattern $h (x, y)$ . To prevent pixel saturation in real-case scenarios, the intensity pixel range is restricted from $q_{min}$ to $q_{max}$ . Denote the minimum and maximum values of pattern $h (x, y)$ as $h_{min}$ and $h_{max}$ , respectively. The image intensity $f [m, n]$ can be expressed as

f [m, n] = Round [q_{min} + (q_{max} - q_{min}) \frac{h (m, n) - h_{m i n}}{h_{m a x - h_{m i n}}}] .

It is worth noting that the quantization process have a subtle effect on the measurement precision for it controls the image contrast. This effect will be discussed in following sections.

2.2. 2D Gaussian digital speckle pattern and non-overlapping condition

Gaussian speckle is extensively used, because of its conceptual simplicity and computational efficiency [21]; the profile of a Gaussian speckle is

ψ (x, y) = exp (- \frac{x^{2} + y^{2}}{r^{2}}),

where r denotes the speckle radius [22]. It is generally assumed that the random offsets, both ϵ and τ, follow an uniform distribution over interval

[- b, b]

, where b is a parameter that controls the maximal offset iteChen:15,Mazzoleni:15; accordingly, the probability distribution functions of both ϵ and τ obey

p_{ϵ} (x) = p_{τ} (x) = \frac{1}{2 b} rect (\frac{x}{2 b}),

where

rect (x)

represents the rectangular function. Digital speckle pattern that satisfies Eqs. (3) and (4) is designated as Gaussian digital speckle pattern.

Gaussian digital speckle pattern can be characterized by three parameters: the grid spacing a, the radius r, and the maximum of random offsets b. The units of these three parameters are the same. The ratio of r to a is a characterization of the percent of area covered by speckles, and the ratio of b to a is a description of the randomness of speckle pattern. To facilitate the study, two dimensionless quantities, the speckle coverage ρ and the speckle variation ν, are defined as follows

ρ = \frac{2 r}{a}, ν = \frac{2 b}{a} .

Fig. 2 If the sum of radius r and maximum of imposed offset b is greater than one half of the grid spacing a, speckles might overlap.

Download Full Size | PDF

Speckle overlap might occur if coverage ρ or variation ν is relatively large. This overlap, however, is generally not desirable, because it tends to decrease the image contrast [12]. Approximate each speckle as a disk and let two speckles just touch, as shown Fig. 2; the speckle radius r, the grid spacing a, and the maximum of imposed offset b should satisfy

b + r = a / 2.

It follows from Eq. (5) that a sufficient condition for speckle non-overlap is

ρ + ν \leq 1.

This is designated as the non-overlapping condition.

2.3. 1D digital speckle pattern and 1D digital speckle image

It is advisable to investigate 1D digital speckle pattern first. Compared with its high-dimensional counterparts, the 1D digital speckle pattern is relatively simple. This simpleness will facilitate the drawing of brief physical picture and the reveal of deep insight. Additionally, it is generally straightforward to extend the 1D analysis to high-dimensional cases.

The 1D digital speckle pattern is mathematically modeled as

h (x) = \sum_{i = - N}^{N} ψ (x - x_{i}) .

Here, $2 N + 1$ is the total number of speckles; $ψ (x)$ is the profile of an individual speckle; $x_{i} = a i + ϵ_{i}$ is the position of the ith speckle, where parameter a characterizes the average distance between two successive speckles and ϵ_i denotes the random offset imposed on the ith speckle.

1D digital speckle image $f [m]$ is generated by directly sampling 1D digital speckle pattern $h (x)$ :

f [m] = h (m) .

It is worth mentioning that 1D digital speckle image is not quantized. The subtle features introduced by quantization process will researched in two-dimensional case.

2.4. 1D Gaussian digital speckle pattern

Analogous to 2D Gaussian speckle patterns, for a 1D Gaussian digital speckle pattern, the speckle profile $ψ (x)$ satisfies

ψ (x) = exp (- x^{2} / r^{2}),

where r denotes the speckle radius, and the probability distribution function of imposed offset ϵ satisfies

p_{ϵ} (ϵ) = {(2 b)}^{- 1} rect (ϵ / 2 b),

where b denotes the maximum of imposed offset.

Fig. 3 Schematic of 1-dimensional digital speckle pattern. Since the speckle centers are random variables, the digital speckle pattern is actually a stochastic process.

Download Full Size | PDF

Indeed, a digital speckle pattern $h (x)$ is a stochastic process, because the centroid of every speckle is random. It means, even the pattern generation parameters are fixed, each time a different pattern would be generated. As an example, four 1D Gaussian digital speckle patterns with identical generation parameters (N = 2, a = 8, $b = 1.2$ , and r = 2) are illustrated in Fig. 3, manifesting the inherent randomness of pattern generation. This spatial uncertainty should be considered and carefully handled [23].

3. Autocorrelation function

Autocorrelation function is of critical importance: it features pattern quality, influences convergence radius, and characterizes subset uniqueness [14]. However, analytical analysis on autocorrelation is scarcely reported. In this section, a theoretical model is built to formulate the autocorrelation function of digital speckle pattern, and approximation formulas are derived to estimate the secondary autocorrelation peak height, a key parameter that features subset uniqueness [11, 12]. The 1D case is investigated first, and then it is extended to 2D case.

3.1. Autocorrelation function of 1D digital speckle pattern

The autocorrelation function, which is a secondary property of a stochastic process, is defined as [20]

R_{h h} (ξ) = E {\int_{- \infty}^{\infty} h (x) h (x + ξ) d x} .

Let $p_{ϵ} (ϵ)$ be the probability distribution function of imposed offset ϵ. Define

ϕ (x) = ψ (x) \otimes p_{ϵ} (x),

where

\otimes

denotes convolution. We prove that the autocorrelation function of 1D digital speckle pattern [see Eq. (8)] is:

\begin{matrix} R_{h h} (ξ) = (2 N + 1) R_{ψ ψ} (ξ) + R_{ϕ ϕ} (ξ) \otimes χ (ξ), \end{matrix}

where

R_{ψ ψ} (ξ)

and

R_{ϕ ϕ} (ξ)

denote the autocorrelations of

ψ (x)

and

ϕ (x)

, respectively, and

χ (ξ) = \sum_{k = - N}^{N} \sum_{l = - N}^{N} δ (ξ + a k - a l) - (2 N + 1) δ (ξ) .

The derivation details can be found in Appendix A.

Fig. 4 An illustrative explanation of the formula of 1D autocorrelation function [Eq. (14)]. The pattern generation parameters are N = 3, a = 8, b = 1.2, and r = 2.

Download Full Size | PDF

The physical meaning of Eq. (14) is schematically illustrated in Fig. 4. Theoretical analysis in Appendix A indicates that (1) term $R_{ϕ ϕ} (ξ) \otimes χ (ξ)$ stems from the cross correlation between different speckles, and it constitutes the sub-peaks of autocorrelation function, (2) term $(2 N + 1) R_{ψ ψ} (ξ)$ comes from the autocorrelation of each speckle with itself, and it forms the primary peak. Function $χ (ξ)$ deserves special attention because it determines the positions of sub-peaks. As Eq. (15) can be rewritten as

\begin{matrix} χ (ξ) = 2 N δ (ξ \pm a) + (2 N - 1) δ (ξ \pm 2 a) + \dots + δ (ξ \pm 2 N a) \\ = \sum_{k = 1}^{2 N} (2 N + 1 - k) [δ (ξ - k a) + δ (ξ + k a)], \end{matrix}

the sub-peaks locate at

\pm a

,

\pm 2 a

,

\dots

,

\pm 2 N a

.

3.2. Secondary autocorrelation peak height of 1D digital speckle pattern

The ratio of secondary to primary autocorrelation peak heights is of essential interest, for it characterize the uniqueness of speckle pattern [14]. As autocorrelation is generally normalized, this ratio is usually known as secondary autocorrelation peak height. An observation of Fig. 4 yields that the primary and secondary peak heights could be approximated as $(2 N + 1) R_{ψ ψ} (0)$ and $2 N R_{ϕ ϕ} (0)$ respectively, if the speckles are not overcrowded. Consequently, the ratio of secondary to primary peak heights could be approximated as

A_{2} \approx \frac{2 N}{1 + 2 N} \frac{R_{ϕ ϕ} (0)}{R_{ψ ψ} (0)} .

Specifically, for 1D Gaussian digital speckle pattern, we prove the ratio of $R_{ϕ ϕ} (0)$ and $R_{ψ ψ} (0)$ is

Φ (k) = \frac{R_{ϕ ϕ} (0)}{R_{ψ ψ} (0)} = \sqrt{\frac{π}{2}} \frac{erf (\sqrt{2} k)}{k} - \frac{1 - e^{- 2 k^{2}}}{2 k^{2}},

where

k = b / r

is a dimensionless quantity that features the randomness of speckle pattern and

erf (x)

is the error function. The derivation details can be referred to Appendix B. Clearly, as

N \to \infty

,

A_{2}

approaches

Φ (k)

.

3.3. Numerical simulation for 1D digital speckle pattern

In order to verify the correctness of proposed theoretical model, the autocorrelation functions of 1D Gaussian digital speckle patterns were calculated numerically. Considering digital speckle pattern’s inherent randomness, $I =$ 1000 patterns were generated for each set of parameters. Denote the ith pattern as $h^{(i)} (x)$ , the autocorrelation function can be evaluated as

R_{h h} (ξ) = \frac{1}{I} \sum_{i = 1}^{I} \int_{- \infty}^{\infty} h^{(i)} (x) h^{(i)} (x + ξ) d x .

The integration in above equation was evaluated numerically using function trapz in MATLAB with an interval space of 0.01. The theoretical estimates of autocorrelation were obtained using Eq. (14).

Fig. 5 Numerical and theoretical results of normalized autocorrelation functions corresponding to (a1) b = 1.0, (a2) b = 2.0, and (a3) b = 3.0; other pattern generation parameters are r = 2, N = 2, and a = 8. (b) Function Φ(k) [Eq. (18)] describes the secondary autocorrelation peak height as N → ∞.

Download Full Size | PDF

The pattern generation parameters were selected as follows: the speckle radius r equaled to 2 pixels, following a suggestion in [16]; to reduce computation cost, parameter N was set as 2; the average in-center spacing a was 8 pixels to ensure a coverage of 50%; the parameter b, which features the maximum of random offsets, was set as 1, 2, and 3 pixels, respectively. The normalized autocorrelation functions ( $R_{h h} (ξ) / R_{h h} (0)$ ) by numerical simulation and by theoretical estimation are depicted in Figs. 5(a1)-5(a3). Numerical results coincide with the theoretical estimates, demonstrating the correctness of proposed formula.

The horizontal dash lines in Figs. 5(a1)-5(a3) represent the secondary autocorrelation peak heights estimated using Eqs. (17) and (18). For instance, consider b = 1 pixel; as the speckle radius r = 2 pixels, parameter $k = b / r$ equals to 0.5; in such a case, $Φ (k)$ [Eq. (18)] equals to 0.924 [see Fig. 5(b)]; as N = 2, factor $2 N / (2 N + 1) = 0.8$ ; hence, according to Eq. (17), the secondary peak height equals to 0.739 approximately [see Fig. 5(a1)]. The theoretical estimates of secondary peak heights agree well with the numerical results. Because $Φ (k)$ is a monotonically decreasing function [see Fig. 5(b)], the secondary peak height decreases with increment of b. From a physical standpoint, an increase in b causes an increase in pattern randomness; hence, each subset becomes more unique; as a result, the secondary peak height decreases.

3.4. Autocorrelation function of 2D digital speckle pattern

The autocorrelation function of a two-dimensional digital speckle pattern $h (ξ, ζ)$ is

R_{h h} (ξ, ζ) = E {\int_{- \infty}^{\infty} \int_{- \infty}^{\infty} h (x, y) h (x + ξ, y + ζ) d x d y} .

Let $p_{ϵ, τ} (ϵ, τ)$ be the probability distribution function of random offset $(ϵ, τ)$ . Define

ϕ (x, y) = ψ (x, y) \otimes p_{ϵ, τ} (x, y),

where

ψ (x, y)

is the profile of an individual speckle. Using a method similar to the 1D case, we prove the autocorrelation function of 2D digital speckle pattern [see Eq. (1)] is

R_{h h} (ξ, ζ) = {(2 N + 1)}^{2} R_{ψ ψ} (ξ, ζ) + R_{ϕ ϕ} (ξ, ζ) \otimes χ (ξ, ζ),

where

R_{ψ ψ} (ξ, ζ)

and

R_{ϕ ϕ} (ξ, ζ)

denote the autocorrelation functions of

ψ (x, y)

and

ϕ (x, y)

, respectively, and

\begin{matrix} χ (ξ, ζ) = \sum_{i = - N}^{N} \sum_{j = - N}^{N} \sum_{k = - N}^{N} \sum_{l = - N}^{N} δ (ξ + a i - a k, ζ + a j - a l) - {(2 N + 1)}^{2} δ (ξ, ζ) \\ = \sum_{s = - 2 N}^{2 N} \sum_{t = - 2 N}^{2 N} (2 N + 1 - | s |) (2 N + 1 - | t |) δ (ξ - s a, ζ - t a) - {(2 N + 1)}^{2} δ (ξ, ζ) . \end{matrix}

Research interests are generally concentrated on sub-peaks near the primary peak; these peaks correspond to small s and t, and therefore $(2 N + 1 - | s |) (2 N + 1 - | t |) \sim 4 N^{2}$ ; moreover, in practical situations, N is generally relatively large, for speckles usually distribute over a large area; hence, function $χ (ξ, ζ)$ can be approximated as

X (ξ, ζ) = \lim_{N \to \infty} χ (ξ, ζ) / {(2 N + 1)}^{2} = comb (ξ, ζ) - δ (ξ, ζ),

where

comb (ξ, ζ)

denotes the comb function. Thus, the autocorrelation function [Eq. (22)] can be approximated as

\lim_{N \to \infty} R_{h h} (ξ, ζ) / {(2 N + 1)}^{2} = R_{ψ ψ} (ξ, ζ) + R_{ϕ ϕ} (ξ, ζ) \otimes X (ξ, ζ) .

3.5. ZNCC-based autocorrelation function of 2D digital speckle image

In the field of image registration, the similarity between reference and target subsets is generally quantified by zero-mean normalized cross-correlation (ZNCC), because of its insensitiveness to offset and scale changes in the intensity of target subset. Given a reference subset $f [m, n]$ and a target subset $g [m, n]$ , the ZNCC is defined as

C_{ZNCC} = \frac{\sum_{m = - M}^{M} \sum_{n = - M}^{M} (f [m, n] - \bar{f}) (g [m, n] - \bar{g})}{\sqrt{\sum_{m = - M}^{M} \sum_{n = - M}^{M} {(f [m, n] - \bar{f})}^{2}} \sqrt{\sum_{m = - M}^{M} \sum_{n = - M}^{M} {(g [m, n] - \bar{g})}^{2}}},

where M denotes one half of subset size,

\bar{f}

and

\bar{g}

denote the intensity average in reference and target subsets, respectively.

An issue then arises: for given digital speckle image $f [m, n]$ , how to estimate its ZNCC-based autocorrelation function $A_{f f} [t, s]$ . Suppose $f [m, n]$ is originated from digital speckle pattern $h (x, y)$ [see Eq. (2)]; denote $h [m, n] = h (m, n)$ as the direct sample of pattern $h (x, y)$ ; since ZNCC is an invariant under linear transformation, ignoring quantization noise, the ZNCC-based autocorrelation functions of $f [m, n]$ and $h [m, n]$ are approximately equal:

A_{f f} [t, s] \approx A_{h h} [t, s] .

Thus, the issue is transformed into estimating $A_{h h} [t, s]$ . In practice, in order to overcome aperture problem, $M ≫ 1$ [24]; it follows from Eq. (26) that

\begin{matrix} A_{h h} [t, s] \approx \frac{\sum_{m = - M}^{M} \sum_{n = - M}^{M} (h [m, n] - \bar{h}) (h [m + t, n + s] - \bar{h})}{\sqrt{\sum_{m = - M}^{M} \sum_{n = - M}^{M} {(h [m, n] - \bar{h})}^{2}} \sqrt{\sum_{m = - M}^{M} \sum_{n = - M}^{M} {(h [m + t, n + s] - \bar{h})}^{2}}} \\ \approx \frac{\sum_{m = - M}^{M} \sum_{n = - M}^{M} h [m, n] h [m + t, n + s] - {(2 M + 1)}^{2} {\bar{h}}^{2}}{\sum_{m = - M}^{M} \sum_{n = - M}^{M} h^{2} [m, n] - {(2 M + 1)}^{2} {\bar{h}}^{2}}, \end{matrix}

where

\bar{h}

denotes the average intensity of image

h [m, n]

. Speckle pattern should be well sampled, and therefore

\begin{matrix} {(2 M + 1)}^{- 2} \sum_{m = - M}^{M} \sum_{n = - M}^{M} h [m, n] h [m + t, n + s] \\ \approx {(2 N a + 1)}^{- 2} \int_{- M}^{M} \int_{- M}^{M} h (x, y) h (x + t, y + s) d x d y \\ \approx a^{- 2} [R_{ψ ψ} (t, s) + (R_{ϕ ϕ} \otimes X) (t, s)] . \end{matrix}

A combination of Eqs. (27), (28), and (29) yields the ZNCC-based autocorrelation function of digital speckle image $f [m, n]$ can be estimated by

A_{f f} [t, s] \approx \frac{R_{ψ ψ} (t, s) + (R_{ϕ ϕ} \otimes X) (t, s) - a^{2} {\bar{h}}^{2}}{R_{ψ ψ} (0, 0) + (R_{ϕ ϕ} \otimes X) (0, 0) - a^{2} {\bar{h}}^{2}} .

3.6. Secondary autocorrelation peak height of 2D digital speckle image

If the speckles are not overcrowded, similar to the 1-dimensional case, the secondary autocorrelation peak height using ZNCC could be approximated as

A_{2} \approx \frac{R_{ϕ ϕ} (0, 0) - a^{2} {\bar{h}}^{2}}{R_{ψ ψ} (0, 0) - a^{2} {\bar{h}}^{2}} .

Specifically, for 2D Gaussian digital speckle image, we prove that above equation can be rewritten as

A_{2} \approx \frac{Φ {(k)}^{2} - π ρ^{2} / 2}{1 - π ρ^{2} / 2},

where

Φ (k)

is defined by Eq. (18) and ρ is speckle coverage defined by Eq. (5). The derivation details can be referred to Appendix C. The existence of expression

- π ρ^{2} / 2

in Eq. (32) is due to the mean value needs to be subtracted when calculating ZNCC.

3.7. Numerical simulation for 2D digital speckle image

Numerical simulations were conducted to verify the analysis on autocorrelation function. Following literature’s recommendation, the radius r was set as 2 pixels and the coverage ρ was set as 50% [16]; the speckle variation ranged from 0% to 100%, in increments of 5%; according to the process described in Section 2, Gaussian digital speckle images were generated, and some of them are shown in Fig. 6. The resolution of these images is 500 × 500; the depth is 8-bit; the minimum and maximum of intensities, $q_{min}$ and $q_{max}$ [see Eq. (2)], were set as 20 and 240 respectively. These parameters are summarized into Table 1; furthermore, every 2D digital speckle image in this article uses these parameters.

Fig. 6 Gaussian digital speckle images with speckle variations ranging from 0% to 100%, given a speckle radius of 2 pixels and a speckle coverage of 50%. The definitions of speckle variation and speckle coverage can be referred to Eq. (5).

Download Full Size | PDF

Table 1. Parameters of Digital Speckle Images.

View Table | View all tables in this article

As shown in Fig. 6(a), in case of variation ν = 0, the digital speckle image is periodic; as the variation increases, this periodicity fades out and the pattern becomes more random; excessive large variation would cause a reduction in image contrast, due to large number of speckles overlap each other.

The subsets were centered at the image center and had a side length of 451 pixels. The displacement along the x axis, t, ranged from -24 to 24 pixels, in increments of 1 pixel; meanwhile, the displacement along the y axis, s, was set as 0 pixel. Figure 7 illustrates the ZNCC-based autocorrelation functions of images shown in Fig. 6. If the variation ν = 0, the autocorrelation function is periodic, for in such a case the pattern is periodic [see Fig. 6(a)]; as variation ν increases, the secondary autocorrelation peak height decreases; for relatively large ν, sub-peaks vanish, probably due to the collapse of underlying periodicity. Sub-peaks, if exist, appear at integer multiples of grid spacing (here a = 8), in accord with the theoretical analysis.

Fig. 7 ZNCC-based autocorrelation functions of images shown in Fig. 6. Here, the speckle variation ranges from 0% to 100%, in increments of 10%.

Download Full Size | PDF

To investigate the relationship between secondary peak height and speckle variation, the correlation coefficients at shift $t = a$ were extracted and are depicted in Fig. 8, superimposed on the theoretical estimates using Eq. (32). Given coverage $ρ = 50 %$ , the critical value of variation $ν = 50 %$ , according to the non-overlapping condition [Eq. (7)]. As shown in Fig. 8, an excellent agreement can be found between the numerical simulation and the theoretical estimation, if the non-overlapping condition is satisfied; if sub-peaks exist, in most cases, the theory presents reasonable estimates. It is worth noting that excessive speckle overlap, which causes disappearance of sub-peaks, is not preferable in practical situations, because it decreases image contrast [see Fig. 6] and thus makes the image more sensitive to sensor noise [25, 26].

Fig. 8 Numerical results and theoretical estimates [Eq. (32)] of secondary autocorrelation peak heights, given a speckle radius of 2 pixels and a speckle coverage of 50%. According to the non-overlapping condition [Eq. (7)], the critical value of speckle variation is 50%.

Download Full Size | PDF

4. Systematic errors of 1D digital speckle pattern

Imperfect interpolation causes a systematic error in digital image correlation technique, termed interpolation bias [27, 28]. This interpolation bias receives extensive concerns, for it can induce virtual strain as large as the 40% of the actual strain level [27]. The aim of this section is to clarify the relationship between interpolation bias and pattern generation parameters; specially, two approximations are proposed.

4.1. Theoretical analysis of 1D interpolation bias

Consider a 1D pattern $h (x)$ . The reference sequence $f [m] = h (m)$ is obtained by directly sampling pattern $h (x)$ ; then, $h (x)$ is shifted along the x axis with u₀ unit, and then is sampled, producing a target sequence $g [m] = h (m - u_{0})$ . The shift patten can be reconstructed by interpolation; if the convolution-based interpolation is used [29], the target pattern can be approximated as

g (x) = \sum_{k = - \infty}^{\infty} g [k] φ (x - k),

where

φ (x)

denotes the interpolation basis.

The displacement can be estimated by

u = arg \min_{u} \sum_{m = - \infty}^{\infty} {f [m] - g (m + u)}^{2} .

If the reconstructed function $g (x)$ is identical to the shifted pattern $h (x - u_{0})$ , the solution of above equation is definitely the actual displacement u₀. However, in practice, interpolation error, namely the difference between $g (x)$ and $h (x - u_{0})$ , is unavoidable, because exact interpolation cannot be performed in practice. This interpolation error will introduce a systematic error $e_{b} = u - u_{0}$ , termed interpolation bias. For the widely used inverse compositional Gauss-Newton (IC-GN) method [30], the interpolation bias e_b could be approximated as [31]

\begin{array}{l} e_{b} \approx c_{b} \sin 2 π u_{0}, \\ c_{b} = \frac{\int_{- π}^{π} E_{b} (ω) | \hat{f} (ω) |^{2} d ω}{\int_{- π}^{π} ω \hat{d} (ω) \hat{φ} (ω) | \hat{f} (ω) |^{2} d ω}, \\ E_{b} (ω) = \hat{d} (ω) [\hat{φ} (ω - 2 π) - \hat{φ} (ω + 2 π)] . \end{array}

where c_b denotes the amplitude of interpolation bias,

E_{b} (ω)

denotes the interpolation bias kernel,

\hat{f} (ω)

denotes the spectrum of the speckle pattern,

\hat{d} (ω)

denotes the magnitude response of the gradient estimator, and

\hat{φ} (ω)

denotes the transfer function of the interpolation algorithm.

As the authors point out, digital speckle pattern is essentially a stochastic process; as a result, the interpolation bias is a random variable $e_{b}$ and so as its amplitude $c_{b}$ ; this randomness is termed spatial uncertainty [23]. To estimate the expectation of interpolation bias, the power spectrum $| \hat{f} (ω) |^{2}$ is substituted by its expectation ${| \hat{f} (ω) |^{2}}$ ; then, the spatial expectation of $c_{b}$ could be approximated as [18]

E {c_{b}} = \frac{\int_{- π}^{π} E_{b} (ω) E {| \hat{f} (ω) |^{2}} d ω}{\int_{- π}^{π} ω \hat{d} (ω) \hat{φ} (ω) E {| \hat{f} (ω) |^{2}} d ω} .

The spatial expectation of power spectrum of 1D digital speckle pattern [see Eq. (8)] is

E {| \hat{f} (ω) |^{2}} = | \hat{ψ} (ω) |^{2} E {| \hat{z} (ω) |^{2}},

where

\hat{ψ} (ω)

is the spectrum of an individual speckle and

\hat{z} (ω) = \sum_{i = - N}^{N} \exp (- j ω x_{i})

features the contribution of speckle positions. Since

\hat{ψ} (ω)

is deterministic and predefined, the key issue should be deriving the specific form of

E {| \hat{z} (ω) |^{2}}

. Through mathematical deduction, we prove for digital speckle pattern

E {| \hat{z} (ω) |^{2}} = (2 N + 1) + {\hat{p}}_{ϵ}^{2} (ω) [D_{N}^{2} (ω a) - (2 N + 1)],

where

{\hat{p}}_{ϵ} (ω)

is the characteristic function of random offset ϵ [20], and

D_{N} (x)

is the N

^{t h}

Dirichlet kernel whose closed form is

D_{N} (x) = \sin (N + \frac{1}{2}) x / \sin \frac{x}{2}

[32]. The derivation details can be referred to Appendix D.

4.2. Numerical simulation for interpolation bias

To demonstrate the validity of proposed theoretical analysis, interpolation biases of 1D digital speckle pattern were evaluated numerically and then compared with theory values.

As shown in Fig. 3, a given set of pattern generation parameters correspond to infinite patterns. To rate this inherent randomness, $I = 1000$ Gaussian speckle patterns were generated for given generation parameters (r, N, a, b). For each speckle pattern $h^{(i)} (x)$ , it was sampled to produce the reference sequence $f^{(i)} [m] = h^{(i)} (m)$ ; then, the pattern was shifted along the x axis in u₀ units, producing a shifted pattern $h^{(i)} (x - u_{0})$ ; the shifted pattern was sampled to produce the target sequence $g^{(i)} [m] = h^{(i)} (x - u_{0})$ ; here, the shift u₀ ranged from 0 to 1 pixel, with an increase of 0.05 pixels; finally, the displacement u was registered using IC-GN method, which is the most popular subpixel registration algorithm at present, because of its high computational efficiency and robustness to sensor noise [33]. The deformation within a subset was described by the first-order shape function, for it can balance the computation cost and the spatial resolution; the target sequence was interpolated using cubic B-spline algorithm, which provides good convergence and high interpolation accuracy; Barron gradient operator, whose truncation error is much smaller than commonly used Sobel and Prewitt operator, was employed to estimate the gradients of reference sequence. These parameters are summarized into Table 2; further more, every registration in this article uses these parameters. To remove boundary effects, the subset size was sufficient large ( $M = ⌈ N a + 2 b + 4 r ⌉$ ), so that all speckles were contained; an even larger subset size will not alter the result, for the newly added intensities are virtually zero, and therefore have no effect on the optimization of Eq. (34). Each pattern would give rise to a different measured displacement; denote $u^{(i)}$ as the measured displacement of the ith pattern; the expectation and variance of the interpolation bias can be evaluated as

E {e_{b}} = \sum_{i = 1}^{I} \frac{u^{(i)} - u_{0}}{I}, Var {e_{b}} = \frac{1}{I - 1} \sum_{i = 1}^{I} {[u^{(i)} - u_{0} - E {e_{b}}]}^{2} .

Table 2. Parameters for Image Registration.

View Table | View all tables in this article

Fig. 9 Interpolation biases corresponding to various pattern generation parameters. (a) r = 2, N = 1, a = 8, b = 2; (b) r = 2, N = 10, a = 8, b = 2; (c) r = 2, N = 1, a = 8, b = 10; (d) r = 3, N =1, a = 8, b = 2.

Download Full Size | PDF

To investigate the role that the generation parameters play, interpolation biases corresponding to various generation parameters were calculated. Since the optimal speckle size falls in the interval of 3 to 5 pixels [22, 24], speckle radius was set as 2 pixels at first, corresponding to a speckle size of 4 pixels; a speckle coverage ρ of 50% is generally preferred [34], and therefore the grid spacing a was set as 8 pixels ( $a = 2 r / ρ$ ); the speckle variation ν was set as 50% as well, and consequently $b =$ 2 pixels ( $b = a ν / 2$ ). Digital image correlation is computational intensive; moreover, in order to study spatial uncertainty, the displacements of I patterns should be registered, further increasing the amount of calculation; to reduce computation burden, N was set as 1 at first, corresponding to totally three speckles.

Figure 9(a) shows the interpolation bias corresponding to parameters mentioned above, namely r = 2, N = 1, a = 8, and b = 2, superimposed on the theoretical estimates given by Eq. (36). Clearly, an excellent agreement is found between the numerical simulation and the theoretical estimation. Then, the influences of generation parameters on interpolation bias are studied by changing parameter N, b, and r. First, to study the influence of speckle count, the speckle count was increased to 21 (N = 10), and the results are shown in Fig. 9(b); a comparison between Fig. 9(a) and Fig. 9(b) indicates that an increase in speckle count would cause a decrease in spatial variance $Var {e_{b}}$ , because patterns with more speckles have increased statistical stability. Second, the parameter b was increased to 10 to investigate the role that random offset plays; by comparing Fig. 9(a) and Fig. 9(c), it shows that the spatial variance $Var {e_{b}}$ increases with increment of b, because increasing b would increase the randomness of speckle pattern. Third, the radius r was increased to 3 pixels; as Fig. 9(a) and Fig. 9(d) indicate, the interpolation bias decreases with increment of speckle radius, due to high-frequency component is the major source of interpolation bias [27, 28]; pattern with relatively large speckle radius has less high frequency components, and therefore corresponds to small interpolation bias.

4.3. Influence of random offset on interpolation bias

Given a speckle radius of 2 pixels and a grid spacing of 8 pixels, interpolation biases at $u_{0} = 0.25$ pixels corresponding to various b were evaluated and are illustrated in Fig. 10, superimposed on the theoretical estimates given by Eqs. (36) and (38). Figures 10(a) and 10(b) correspond to N = 1 and $N = 10$ , respectively. According to Fig. 10, with increment of b, the interpolation bias increases firstly and then decreases; finally, when b is relatively large, the bias appears to approach a specific limit value. Indeed, this value can be estimated theoretically.

Fig. 10 Interpolation biases at u₀ = 0.25 pixels corresponding to various random offsets b when (a) N = 1 and (b) N = 10, given a speckle radius of 2 pixels and a grid spacing of 8 pixels.

Download Full Size | PDF

The imposed displacement ϵ follow an uniform distribution in interval $[- b, b]$ ; its characteristic function ${\hat{p}}_{ϵ} (ω) = sinc b ω$ ; thus, ${\hat{p}}_{ϵ}^{2} (ω) \sim b^{- 2}$ as $b \to \infty$ . It follows from Eq. (38) that

\lim_{b \to \infty} E {| \hat{z} (ω) |^{2}} / (2 N + 1) = 1, for any ω \neq 0.

The interpolation bias kernel $E_{b} (ω)$ equals to zero at ω = 0, because there is no interpolation error for a constant function [28]. Consequently, Eq. (36) can be rewritten as

\lim_{b \to \infty} E {c_{b}} = \frac{\int_{- π}^{π} E_{b} (ω) | \hat{ψ} (ω) |^{2} d ω}{\int_{- π}^{π} ω \hat{d} (ω) \hat{φ} (ω) | \hat{ψ} (ω) |^{2} d ω} .

The horizontal dash lines in Fig. 10 represent the values estimated by above equation. Clearly, Eq. (41) is a good description of the bias error in case of relatively large b.

Interestingly, Eq. (41) is identical to the formula of bias expectation of patterns that are composed of completely randomly positioned speckles [18]. This phenomenon can be explained as follows: in case of relatively large b, the inherent regularity controlled by the grid spacing a is overwhelmed; that is, digital speckle pattern tends to completely random pattern as $b \to \infty$ ; therefore, these two formulas are the same.

The theoretical estimates fit the numerical results better, when N is larger or b is small. This is due to the approximation used in derivation: consider a random variable $z = x / y$ , where both x and y

are random variables; because it is difficult to derive the specific form of $E {z}$ , it is approximated as $E {z} \approx E {x} / E {y}$ [from Eq. (35) to Eq. (36)]. This approximation is more accurate when the variations of x and y are less significant. Increasing b would increase pattern randomness; decreasing speckle count would decrease the statistical stability; therefore, our theoretical model is more accurate for large speckle count N and small random offset b.

4.4. Influence of speckle count on interpolation bias

Interpolation biases at $u_{0} = 0.25$ pixels corresponding to various speckle counts were calculated and are illustrated in Fig. 11, given a speckle radius of 2 pixels and a grid spacing of 8 pixels. Figures 11(a) and 11(b) correspond to b = 2 and $b = 10$ , respectively. The theoretical estimates, given by Eqs. (36) and (38), are consistent with the numerical results. With increment of speckle count, the systematic error quickly approaches a certain value. An issue is how to estimate this value.

Fig. 11 Interpolation biases at u₀ = 0.25 pixels corresponding to various speckle count N when (a) b = 2 and (b) b = 10, given a speckle radius of 2 pixels and a grid spacing of 8 pixels.

Download Full Size | PDF

As $N \to \infty$ , the Dirichlet kernel obeys

\lim_{N \to \infty} D_{N}^{2} (ω a) / (2 N + 1) = \frac{2 π}{a} comb (\frac{a ω}{2 π}) .

Substituting above equation into Eq. (38) yields

\lim_{N \to \infty} E {| \hat{z} (ω) |^{2}} / (2 N + 1) = 1 - {\hat{p}}_{ϵ}^{2} (ω) + \frac{2 π}{a} {\hat{p}}_{ϵ}^{2} (ω) comb (\frac{ω a}{2 π}) .

Hence, the limit of interpolation bias as N approaches infinity can be approximated as

\lim_{N \to \infty} E {c_{b}} = \frac{\int_{- π}^{π} E_{b} (ω) | \hat{ψ} (ω) |^{2} [1 - {\hat{p}}_{ϵ}^{2} (ω) + 2 π a^{- 1} {\hat{p}}_{ϵ}^{2} (ω) comb (ω a / 2 π)] d ω}{\int_{- π}^{π} ω \hat{d} (ω) \hat{φ} (ω) | \hat{ψ} (ω) |^{2} [1 - {\hat{p}}_{ϵ}^{2} (ω) + 2 π a^{- 1} {\hat{p}}_{ϵ}^{2} (ω) comb (ω a / 2 π)] d ω} .

The horizontal dash lines in Fig. 11 depict the estimates using above equation, demonstrating its validity.

5. Measurement errors of 2D digital speckle image

Sensor noise, which is inevitable in practice, would cause a measurement uncertainty, termed random error [25]. Hence, in practice, the measurement error is composed of two parts: the systematic error and the random error [26]. In this section, these two errors will be analyzed theoretically.

5.1. Theoretical analysis of 2D systematic error

It is straightforward to extend the 1D theoretical analysis on interpolation bias to the 2D case. The systematic error of 2D digital speckle pattern can be estimated as

E {c_{b}} = \frac{\int_{- π}^{π} \int_{- π}^{π} E_{b} (ω_{x}, ω_{y}) | \hat{ψ} (ω_{x}, ω_{y}) |^{2} E {| \hat{z} (ω_{x}, ω_{y}) |^{2}} d ω_{x} d ω_{y}}{\int_{- π}^{π} \int_{- π}^{π} ω_{x} \hat{d} (ω_{x}, ω_{y}) \hat{φ} (ω_{x}, ω_{y}) | \hat{ψ} (ω_{x}, ω_{y}) |^{2} E {| \hat{z} (ω_{x}, ω_{y}) |^{2}} d ω_{x} d ω_{y}},

where

E_{b} (ω_{x}, ω_{y})

denotes the interpolation bias kernel,

\hat{ψ} (ω_{x}, ω_{y})

is the spectrum of an individual speckle,

\hat{d} (ω_{x}, ω_{y})

is the magnitude response of gradient estimator,

\hat{φ} (ω_{x}, ω_{y})

is the interpolation transfer function, and

E {| \hat{z} (ω_{x}, ω_{y}) |^{2}} = {(2 N + 1)}^{2} + {\hat{p}}_{ϵ}^{2} (ω_{x}) {\hat{p}}_{τ}^{2} (ω_{y}) [D_{N}^{2} (ω_{x} a) D_{N}^{2} (ω_{y} a) - {(2 N + 1)}^{2}] .

Here, ${(2 N + 1)}^{2}$ is total number of speckles, ${\hat{p}}_{ϵ} (ω_{x})$ and ${\hat{p}}_{τ} (ω_{y})$ are the characteristic functions of random offsets ϵ and τ respectively, and $D_{N} (x)$ is the Nth Dirichlet kernel.

Equations (45) and (46), however, are hard to use, because in practice the subset is only a small part of the whole speckle area other than contains it. Inspired by the discussion in Section 4, the limits of $b \to \infty$ and $N \to \infty$ are suggested to estimate the interpolation bias, consequently

E {c_{b}} = \frac{\int_{- π}^{π} \int_{- π}^{π} E_{b} (ω_{x}, ω_{y}) | \hat{ψ} (ω_{x}, ω_{y}) |^{2} Θ (ω_{x}, ω_{y}) d ω_{x} d ω_{y}}{\int_{- π}^{π} \int_{- π}^{π} ω_{x} \hat{d} (ω_{x}, ω_{y}) \hat{φ} (ω_{x}, ω_{y}) | \hat{ψ} (ω_{x}, ω_{y}) |^{2} Θ (ω_{x}, ω_{y}) d ω_{x} d ω_{y}},

where

Θ (ω_{x}, ω_{y}) = 1, if b \to \infty

and

\begin{array}{l} Θ (ω_{x}, ω_{y}) = 1 & - \overset{2}{sinc} b ω_{x} \overset{2}{sinc} b ω_{y} \\ + \frac{4 π^{2}}{a^{2}} \overset{2}{sinc} b ω_{x} \overset{2}{sinc} b ω_{y} comb (\frac{ω_{x} a}{2 π}) comb (\frac{ω_{y} a}{2 π}), \end{array} if N \to \infty .

Here, it is assumed that random offsets, both ϵ and τ, follow an uniform distribution in interval $[- b, b]$ , so that the characteristic function ${\hat{p}}_{ϵ} (ω) = {\hat{p}}_{τ} (ω) = sinc (b ω)$ . Equations (48) and (49) are the two-dimensional counterparts of Eqs. (40) and (43), respectively. Since the interpolation bias does not change with the linear transformation of the gray scale, Eq. (47) also applies to 2D digital speckle image.

5.2. Theoretical analysis of 2D random error

The random error has been thoroughly studied [25, 26]. For the widely used IC-GN algorithm, the random error e_n could be estimated as [33]

e_{n}^{2} = \frac{σ_{f}^{2} + σ_{g}^{2}}{\sum_{m = - M}^{M} \sum_{n = - M}^{M} f_{x}^{2} [m, n]},

where σ_f and σ_g denote the standard deviation of sensor noise in reference and target images respectively, M is the one half of subset size, and

f_{x} [m, n]

denotes the image gradient along the x axis. Considering the spatial uncertainty, Eq. (50) could be rewritten as [18]

e_{n}^{2} = \frac{σ^{2}}{Q_{f}},

where

σ^{2} = σ_{f}^{2} + σ_{g}^{2}

characterizes the noise level, and

Q_{f} = \sum_{m = - M}^{M} \sum_{n = - M}^{M} f_{x}^{2} [m, n]

is the sum of square of subset intensity gradients (SSSIG) [25]. A large

Q_{f}

is generally preferred, because according to Eq. (51) increased

Q_{f}

is able to reduce the random error.

The key to estimate the random error is to derive the form of $Q_{f}$ . Suppose that digital speckle image $f [m, n]$ comes from digital speckle pattern $h (x, y)$ [see Eq. (2)]. Ignoring quantization noise, it follows from Eq. (2) that

Q_{f} \approx \frac{{(q_{max} - q_{min})}^{2}}{{(h_{max} - h_{min})}^{2} Q_{h}},

where

Q_{h} = \sum_{m = - M}^{M} \sum_{n = - M}^{M} h_{x}^{2} (m, n)

. Through mathematical derivation, we prove

Q_{h}

could be approximated as

Q_{h} \approx \frac{{(2 M + 1)}^{2}}{4 π^{2} a^{2}} \int_{- π}^{π} \int_{- π}^{π} | \hat{d} (ω_{x}, ω_{y}) |^{2} | \hat{ψ} (ω_{x}, ω_{y}) |^{2} Θ (ω_{x}, ω_{y}) d ω_{x} d ω_{y},

where

2 M + 1

is the subset size; a denotes the grid spacing;

\hat{d} (ω_{x}, ω_{y})

denotes the magnitude response of gradient estimator;

\hat{ψ} (ω_{x}, ω_{y})

denotes the spectrum of an individual speckle;

Θ (ω_{x}, ω_{y})

is given by Eq. (48) or Eq. (49), depending on which assumption is used. The derivation details can be referred to Appendix E.

A combination of Eqs. (51), (52) and (53) yields the random error can be estimated as

\begin{array}{l} E {e_{n}^{2}} \approx \frac{4 π^{2} a^{2} σ^{2} {(h_{m a x} - h_{m i n})}^{2}}{{(2 M + 1)}^{2} {(q_{m a x} - q_{m i n})}^{2}} \times {[\int_{- π}^{π} \int_{- π}^{π} | \hat{d} (ω_{x}, ω_{y}) |^{2} | \hat{ψ} (ω_{x}, ω_{y}) |^{2} Θ (ω_{x}, ω_{y}) d ω_{x} d ω_{y}]}^{- 1} . \end{array}

If speckle overlap occurs, it would be difficult to theoretically estimate the specific values of $h_{max}$ and $h_{min}$ , because of their inherent randomness. Conversely, if speckles do not overlap each other, the analysis is greatly simplified, because in such a case

h_{max} = ψ_{max}, h_{min} = ψ_{min} .

Here, $ψ_{min}$ and $ψ_{max}$ denote the the minimum and maximum of speckle $ψ (x, y)$ , respectively. In this way, Eq. (54) can be rewritten as

\begin{array}{l} E {e_{n}^{2}} \approx \frac{4 π^{2} a^{2} σ^{2} {(ψ_{max} - ψ_{min})}^{2}}{{(2 M + 1)}^{2} {(q_{max} - q_{min})}^{2}} \end{array} \times {[\int_{- π}^{π} \int_{- π}^{π} | \hat{d} (ω_{x}, ω_{y}) |^{2} | \hat{ψ} (ω_{x}, ω_{y}) |^{2} Θ (ω_{x}, ω_{y}) d ω_{x} d ω_{y}]}^{- 1} .

5.3. Numerical simulation for measurement errors

To verify the correctness of proposed theoretical formulas, systematic errors and random errors were evaluated numerically and then compared with the theoretical estimates. The calculation process of measurement errors is shown in Fig. 12: (1) generate a two-dimensional Gaussian speckle pattern $h (x, y)$ using given pattern generation parameters; (2) generate the translated pattern $h (x - u_{0}, y)$ by shifting the original pattern along the x axis with u₀ units; (3) produce the reference image $f [m, n]$ and target image $g [m, n]$ from the original pattern $h (x, y)$ and translated pattern $h (x - u_{0}, y)$ , respectively; (4) add random Gaussian noise to both the reference and the target images to simulate sensor noise; (5)register the full field displacement using IC-GN algorithm by comparing the noisy reference image and the noisy target image; (6) repeat step 4 and step 5 for J times to assess the temporal uncertainty.

Fig. 12 Flowchart of calculation of measurement errors for 2D digital speckle pattern.

Download Full Size | PDF

The full field displacement is registered (step 5) as follows: in the region of interest (ROI), with a given grid step, a grid of point of interest (POI) is identified, as shown in Fig. 13; then, for each POI, a subset centered at the POI is selected, and the displacement is registered using the IC-GN algorithm with parameters in Table 2. Parameters of images are shown in Table 1; other calculation parameters are shown in Table 3.

Fig. 13 The image resolution is 500× 500 pixels. The region of interest (ROI), indicated by the yellow frame, has a side size of 450 pixels. The grid step is 10 pixels. The point of interest (POI) is indicated by the red squares.

Download Full Size | PDF

Table 3. Parameters for Evaluation of 2D Measurement Errors.

View Table | View all tables in this article

Denote $u^{(i, j)}$ as the measured displacement of the ith POI corresponding to the jth noise addition. The systematic error $e_{b^{(i)}}$ and the random error $e_{n^{(i)}}$ of the ith POI could be estimated as

e_{b}^{(i)} = \frac{1}{J} \sum_{j = 1}^{J} [u^{(i, j)} - u_{0}], e_{n}^{(i)} = \sqrt{\frac{1}{J - 1} \sum_{j = 1}^{J} {[u^{(i, j)} - u_{0} - e_{b^{(i)}}]}^{2}},

where J denotes the total number of noise addition. Different POIs have different calculation errors due to spatial uncertainty. The spatial expectation and variance of systematic error and random error are

\begin{array}{l} E {e_{b}} = \frac{1}{I} \sum_{i = 1}^{I} e_{b}^{(i)}, Var {e_{b}} = \frac{1}{I - 1} \sum_{i = 1}^{I} {[e_{b}^{(i)} - E {e_{b}}]}^{2}; \\ E {e_{n}} = \frac{1}{I} \sum_{i = 1}^{I} e_{n}^{(i)}, Var {e_{n}} = \frac{1}{I - 1} \sum_{i = 1}^{I} {[e_{n}^{(i)} - E {e_{n}}]}^{2} . \end{array}

Fig. 14 Gaussian digital speckle images with radii (a) 1.0, (b) 1.5, (c) 2.0, (d) 2.5, (e) 3.0, and (f) 3.5 pixels, given a speckle coverage of 50% and a speckle variation of 50%.

Download Full Size | PDF

5.4. Influence of speckle radius on measurement errors

First, measurement errors of different radii were investigated. Excessive small speckle radius leads to significant interpolation bias; excessive large radius gives rise to large random error; the optimal radius falls within the range of 1.5 to 2.5 pixels. Figures 14(a)-14(f) shows digital speckle images with radii 1.0, 1.5, 2.0, 2.5, 3.0 and 3.5 pixels, given a speckle coverage of 50% and a variation of 50%. The measurement errors of these images were evaluated numerically. As $h_{min}$ and $h_{max}$ are difficult to estimate, Eq. (55) was employed as an approximation. Based on Eqs. (47) and (56), the measurement errors were estimated. Because function Θ(ω_x, ω_y) depends on the assumption adopted, there are two theoretical estimates, corresponding to $b \to \infty$ and $N \to \infty$ respectively.

Fig. 15 Systematic errors and random errors for speckle radii (a) 1.5, (b) 2.0, and (c) 3.0 pixels, corresponding to images shown in Fig. 14(b), 14(c), and 14(e).

Download Full Size | PDF

Figures 15(a)-15(c) illustrates systematic errors and random errors corresponding to radii 1.5, 2,0, and 3.0 pixels. A good agreement is found between the theory values and the numerical results. Both approximations, whether $b \to \infty$ or $N \to \infty$ , present valid estimates. For the systematic errors, compared with $b \to \infty$ , estimates of $N \to \infty$ seems closer to the numerical results, probably due to the bias converges quickly as N increases [see Fig. 11]. The difference between $b \to \infty$ and $N \to \infty$ becomes less significant as the radius increases. For the random errors, the estimates of $N \to \infty$ are slightly larger than that of $b \to \infty$ , but both of them provide reasonable estimates.

Fig. 16 (a) Systematic errors and (b) random errors corresponding to various speckle radii, given a speckle coverage of 50% and a speckle variation of 50%.

Download Full Size | PDF

To investigate the relationship between measurement errors and speckle radius, the measurement errors at $u_{0} = 0.25$ pixels were extracted and are illustrated in Fig. 16. As the speckle radius increases, the systematic error decreases and the random error increases, similar to the phenomenon observed in [18]. Theoretical estimates of systematic errors are in good agreement with the numerical results, except for r = 1.0 pixel, due to the sinusoidal form of systematic error is derived based on an assumption that sampling theorem is satisfied [31]. Theoretical estimates of random errors are consistent with the numerical results. Specifically, in the optimal radius interval $r \in [1.5, 2.5]$ , proposed theory provides rather good estimates for random errors.

With increment of speckle radius, the systematic error decreases and the random error increases. Hence, there is an optimal radius that minimizes the total error. This optimized value could be readily identified by making use of the theoretical analysis.

5.5. Influence of speckle variation on measurement errors

Second, the influence of speckle variation on measurement errors was investigated. The variation ranged from 10% to 80%, in increments of 10%. The radius was set as 2 pixels and the coverage was set as 50%. Gaussian digital speckle images generated are shown in Fig. 6. Measurement errors at $u_{0} =$ 0.25 pixels were calculated and are shown in Fig. 17, superimposed on the theoretical estimates using Eqs. (47) and (56).

Fig. 17 (a) Systematic errors and (b) random errors corresponding to various speckle variations, given a speckle radius of 2 pixels and a speckle coverage of 50%.

Download Full Size | PDF

The theoretical estimates of systematic errors agree well with the numerical results. If the non-overlapping condition is satisfied, theoretical estimates of random errors are in accordance with numerical results. Nevertheless, if a large number of speckles overlap each other, Eq. (55) fails to hold; consequently, Equation (56) becomes invalid. Unfortunately, it is difficult to estimate $h_{min}$ and $h_{max}$ in case of non-overlapping condition is not satisfied.

The estimates using assumption $b \to \infty$ are independent of the variation, because it describes the limit behavior as variation $ν \to \infty$ . By contrast, estimates of $N \to \infty$ characterize the role speckle variation plays. Compared with estimates using $b \to \infty$ , the estimates using $N \to \infty$ are in better accordance with the numerical results, probably due to its rapid convergence [see Fig. 11]. Nevertheless, the discrepancy between these two assumptions is not significant.

Figure 17 implies that if a regular pattern (ν = 0) is used, the measurement errors can be minimized. Although regular pattern can be used for image registration, there is a high probability that a mismatch may occur, leading to erroneous measurements [11, 12].

Fig. 18 Digital speckle images with coverages (a) 30%, (b) 40%, (b) 50%, (d) 60%, (e) 70%, and (f) 80%, given a speckle radius of 2 pixels and a speckle variation of 50%.

Download Full Size | PDF

5.6. Influence of speckle coverage on measurement errors

Third, the influence of speckle coverage on measurement errors was investigated. Figure 18 shows digital speckle images with coverages ranging from from 30% to 80%, given a speckle radius of 2 pixels and a variation of 50%. The measurement errors at $u_{0} = 0.25$ pixels were calculated and are illustrated in Fig. 19, superimposed by the estimates using Eqs. (47) and (56).

The systematic error increases with increment of speckle coverage. The estimates using $N \to \infty$ coincide well with the calculation. By contrast, the estimate using $b \to \infty$ is a constant because $Θ (ω_{x}, ω_{y}) = 1$ . As the speckle coverage increases, the random error decreases firstly and then increases. If the non-overlapping condition is satisfied, a good agreement could be found between the random errors by theoretical estimation and by numerical calculation. If overlap occurs, the estimates deviate from the numerical results, because the maximum and minimum of pattern, $h_{min}$ and $h_{max}$ , cannot be effectively estimated.

Fig. 19 (a) Systematic errors and (b) random errors corresponding to various speckle coverages, given a speckle radius of 2 pixels and a speckle variation of 50%.

Download Full Size | PDF

Here, compared with the random error, the systematic error is relatively small, because the speckle radius is 2 pixels, which is recommended by literature [16]. In such a case, the random error is the major source of the total error. Consequently, the optimal coverage is about 60% [see Fig. 19(b)], in which case overlaps do not occur in large numbers, and therefore the contrast increase caused by incremental speckle count overcomes the contrast lost due to overlaps. Nevertheless, the enhancement from $ν = 50 %$ to $ν = 60 %$ is not significant, and therefore 50% is also a good choice in practice.

6. Spatial resolution

Spatial resolution characterizes the ability of a given technique to distinguish close features [35]. A definition of spatial resolution is the period of the sine displacement beyond which the relative loss in amplitude of the displacement returned by the measuring technique is greater than a certain value [10]. Recently, the spatial resolution has evolved into a research hotspot [36, 37]. The aim of this section is to explore the relationship between spatial resolution and subset size.

6.1. Spatial resolution and Savitzky-Golay filter

Consider a sinusoidal deformation field

u = A_{0} sin ω x .

The measured amplitude $A (ω)$ is a monotonically decreasing function of ω. Given a threshold α, the limit frequency ω_α obeys [10]

A (ω_{α}) = α A_{0} .

Spatial resolution is the period corresponding to the critical frequency:

d_{u} = 2 π / ω_{α} .

Grédiac et al. proposed to identify the spatial resolution using the frequency response of Savitzky-Golay (S-G) filter [9], on the basis that the measured displacement field is approximately equal to the convolution of the actual displacement field and the S-G filter [38]. A S-G filter is characterized by two parameters: polynomial order d and half width M [39]. In terms of digital image correlation, these two parameters correspond to shape function order and half subset size, respectively [38].

The frequency responses of S-G filters, with various d and M, are illustrated in Fig. 20. The horizontal and vertical axes represent the normalized wave number $\tilde{k} = ω / π$ and the frequency response $\hat{H} (\tilde{k})$ of the S-G filter in dB, respectively. The critical value of normalized wave number ${\tilde{k}}_{α}$ obeys [9]

\hat{H} ({\tilde{k}}_{α}) = α,

where α is the threshold. Accordingly, the spatial resolution is [9]

d_{u} = 2 / {\tilde{k}}_{α} .

Grédiac et al. used a threshold of $α =$ 0.2, equal to -14 dB approximately [9]. As shown in Fig. 20, the critical value of normalized wave number, ${\tilde{k}}_{α}$ , increases with increment of shape function order d or with decrement of half subset size M.

Fig. 20 Frequency response of S-G filters for various subset sizes and shape function orders.

Download Full Size | PDF

6.2. Relationship between spatial resolution and subset size

In order to achieve desired spatial resolution, the subset size should not exceed a specific value. To identify this value, the relationship between spatial resolution and subset size needs to be clarified.

Fig. 21 Relation between spatial resolution d_u and half subset size M for both first- and second-order shape functions. With least squares fitting method, empirical formulas [Eq. (64)] are proposed.

Download Full Size | PDF

Figure 21 shows the spatial resolutions corresponding to various subset sizes for both first- and second-order shape functions, calculated by Eqs. (62) and (63) with a threshold of $α = 0.2$ . Interestingly, a good linear relationship is observed between the half subset size M and the spatial resolution d_u. With least squares fitting method, empirical formulas for spatial resolution are proposed:

d_{u} = {\begin{array}{l} 2.4211 M + 1.1853, & first - order shape function; \\ 1.2998 M + 0.5443, & second - order shape function . \end{array}

If a spatial resolution of d_u is desired, the required subset size $2 M + 1$ can be identified by above equation.

6.3. A rudimentary model for spatial resolution

It is rather difficult to deduce a closed form formula that can precisely characterize the spatial resolution, because the frequency-response of S-G filter has strong nonlinearity. Instead, we build a rudimentary model, based on an fact that the limited spatial resolution originates from the shape function under-match [9, 38].

For the first-order shape function, the error due to shape function under-match is approximately [40]

Δ u \approx \frac{M (M + 1)}{6} (\frac{\partial^{2} u}{\partial^{2} x} + \frac{\partial^{2} u}{\partial^{2} y}),

where u is the actual displacement field and M is half subset size. Hence, for the sinusoidal deformation field described by Eq. (59), the under-match error

Δ u \approx - \frac{M (M + 1)}{6} A_{0} ω^{2} sin ω x .

Then, the measured amplitude is

A (ω) \approx A_{0} - \frac{M (M + 1)}{6} A_{0} ω^{2} .

Given a threshold α, with Eqs. (60) and (61), the spatial resolution for the first-order shape function could be approximated as

d_{u}^{(1)} \approx 2 π \sqrt{\frac{M (M + 1)}{6 (1 - α)}} \sim \frac{2 π}{\sqrt{6 (1 - α)}} (M + \frac{1}{2} - \frac{1}{8 M} + \dots) .

Above equation might explain why in Eq. (64) the slope of first-order shape function is about two times of its intercept.

For the second-order shape function, the under-match error is [41]

Δ u = \frac{M^{2} {(M + 1)}^{2}}{2520} (9 \frac{\partial^{4} u}{\partial x^{4}} + 70 \frac{\partial^{4} u}{\partial x^{2} \partial y^{2}} + 9 \frac{\partial^{4} u}{\partial y^{4}}) .

Analogously, after a simple deduction, the spatial resolution for the second-order shape function could be approximated as

d_{u}^{(2)} \approx 2 π \frac{\sqrt{M (M + 1)}}{\sqrt[4]{280 (1 - α)}} \sim \frac{2 π}{\sqrt[4]{280 (1 - α)}} (M + \frac{1}{2} - \frac{1}{8 M} + \dots) .

It is worth mentioning that the under-match error formulas, both Eqs. (65) and (69), are deduced by ignoring high order terms of Taylor series [40, 41]. Therefore, these formulas are more accurate for smooth deformation and small subset size.

Fig. 22 Spatial resolutions by numerical calculation and by theory [Eqs. (68) and (70)] for threshold (a1) α = 0.99, (a2) α = 0.5, and (a3) α = 0.2. (b) Slopes by least squares fitting and by theory [Eq. (71)].

Download Full Size | PDF

To investigate the feasibility of the rudimentary model, the spatial resolutions corresponding to α = 0.99, 0.5, and 0.2 were calculated and are depicted in Figs. 22(a1)-22(a3) respectively, superimposed on the theoretical estimates given by Eqs. (68) and (70). In case of α = 0.99, as shown in Fig. 22(a1), an excellent agreement could be found between the theoretical estimates and the numerical results. As α decreases, the difference between theoretical and numerical results increases for both first- and second-order shape functions, probably due to for significant local deformations Eqs. (65) and (69) become invalid [40]. The value of spatial resolution is preferred to be small, because it means closer features could be figured out. Since the theory values are greater than the numerical results, the proposed model provides conservative estimates for spatial resolution.

The slope of the best fitting straight line is a key parameter characterizing the relationship between spatial resolution and subset size. According to Eqs. (68) and (70), the slopes could be approximated as

k^{(1)} = 2 π {[6 (1 - α)]}^{- 1 / 2} and k^{(2)} = 2 π {[280 (1 - α)]}^{- 1 / 4}

for first- and second-order shape functions, respectively. On the other hand, the slopes could be numerically evaluated using least squares fitting. Figure 22(b) shows the slopes for various thresholds by calculation and by Eq. (71). For relatively large α, the theoretical estimates of slopes agree with the numerical results excellently; for small α, the theoretical estimates are less than the numerical results, and thus the proposed model provides conservative estimates for spatial resolution, probably due to for sinusoidal deformations the under-match errors estimated using Eq. (65) are larger than the actual errors [40].

7. Discussion

7.1. Suggestions on selection of generation parameters of digital speckle pattern

Digital speckle pattern is characterized by three parameters: the speckle radius r, the speckle coverage ρ, and the speckle variation ν. In previous sections, four factors, including subset uniqueness, measurement accuracy, measurement precision, and spatial resolution, are researched. Considering these four factors, we proposed several suggestions on the selection of pattern generation parameters.

Based on the non-overlapping condition [Eq. (7)], it is suggested that the coverage ρ and variation ν should satisfy $ρ + ν = 1$ . Although numerical simulation suggests a slight larger coverage, the improvement is not significant, and that optimal coverage is hard to identify.
Secondary autocorrelation peak height $A_{2}$ , largely determined by k and ρ [Eq. (32)], features the subset uniqueness. Using the first suggestion, parameter $k = b / r = ν / ρ$ can be expressed as $k = (1 - ρ) / ρ = ρ^{- 1} - 1$ . Thus, $A_{2}$ is a function of ρ exclusively. Given preferred $A_{2}$ , the speckle coverage ρ can be numerically identified, and then the variation $ν = 1 - ρ$ .
The spatial resolution is mainly determined by the subset size. Given desired spatial resolution, the subset size can be identified using the empirical formulas [Eq. (64)] or by the rudimentary model [Eqs. (68) and (70)].
As the speckle radius r increases, the systematic error decreases and the random error increases [Fig. 16]. Hence, An optimal speckle radius exists. Given coverage ρ, variation ν, and subset size $2 M + 1$ , the overall error, which includes both systematic error and random error, is a function of radius r exclusively. With the theoretical analysis on systematic errors and random errors, the optimal radius can be determined.

7.2. Limits and outlook

Although we try to present a comprehensive study on the performance of digital speckle pattern, there are still two issues needing further research. (1) Some formulas, such as Eqs. (32) and (56), become invalid if large number of speckles overlap each other. Though speckle overlap is generally not preferable because it may cause potential decrease in image contrast, this is still a problem if we want to build a complete framework. (2) The influence of imaging system is not taken into account. Actually, the object of study of this work is the pattern on the sensor plane, not the original pattern on the specimen. Several preliminary experiments concerning imaging system have be carried out, and in our following work the influence of imaging system will be thoroughly discussed.

8. Conclusion

In this work, we build theoretical models for autocorrelation function, systematic error, random error, and spatial resolution. Considering these four factors, we propose suggestions on selection of generation parameters of digital speckle pattern. To provide more implementation details, all source codes and original data are uploaded and ready to be downloaded. The conclusions of this work are drawn as follows:

In terms of uniqueness, a theoretical model is built for autocorrelation function of digital speckle pattern, and formulas are derived to estimate the secondary autocorrelation peak height, which is a characteristic of pattern uniqueness. With these formulas, desired uniqueness can be obtained by choosing the appropriate pattern generation parameters.
In terms of measurement accuracy and precision, by deriving the spatial expectation of power spectrum of digital speckle pattern, approximate formulas are presented to characterize systematic errors and random errors. These formulas will help scholars evaluate accuracy and precision of digital image correlation technique.
In terms of spatial resolution, the linear relation between spatial resolution and subset size is featured by empirical formulas for both first- and second-order shape functions; beside, a rudiment model for spatial resolution, which provides conservative estimates, is presented. With these formulas, given desired spatial resolution, we can identify the required subset size.
Considering uniqueness, accuracy, precision, and spatial resolution, suggestions for selection of pattern generation parameters are proposed, which can eventually improve the measurement performance of digital image correlation technique.

These contributions will facilitate the standardization of speckle pattern, which is our primary motivation.

Appendix

A. Autocorrelation function of 1D digital speckle pattern

The autocorrelation function, which is a secondary property of a stochastic process, is defined as

R_{h h} (ξ) = E {\int_{- \infty}^{\infty} h (x) h (x + ξ) d x} .

Hence, the autocorrelation function of 1D digital speckle pattern [see Eq. (8)] is

\begin{array}{l} R_{h h} (ξ) = \sum_{k = - N}^{N} \sum_{l = - N}^{N} α_{k, l} (ξ), \\ α_{k, l} (ξ) = E {\int_{- \infty}^{\infty} ψ (x - x_{k}) ψ (x + ξ - x_{l}) d x} . \end{array}

Function $α_{k, l} (ξ)$ characterizes the correlation between the kth and the lth speckles. Particularly, if $k = l$ , it becomes the correlation of an individual speckle with itself:

α_{k, k} (ξ) = E {\int_{- \infty}^{\infty} ψ (x) ψ (x + ξ) d x} = R_{ψ ψ} (ξ) .

Let $p_{ϵ} (ϵ)$ be the probability distribution function of ϵ. Define $ϕ (x)$ as the convolution of $p_{ϵ} (ϵ)$ and speckle profile $ψ (x)$ :

ϕ (x) = ψ (x) \otimes p_{ϵ} (x),

then the expectation

E {ψ (x - ϵ)} = \int_{- \infty}^{\infty} ψ (x - ϵ) p_{ϵ} (ϵ) d ϵ = ϕ (x) .

If $k \neq l$ , since ϵ_k and ϵ_l are independent, $α_{k, l} (ξ)$ can be given by

\begin{matrix} α_{k, l} (ξ) = \int_{- \infty}^{\infty} E {ψ (x - x_{k})} E {ψ (x + ξ - x_{l})} d x \\ = \int_{- \infty}^{\infty} ϕ (x - a k) ϕ (x + ξ - a l) d x \\ = R_{ϕ ϕ} (ξ + a k - a l), \end{matrix}

where

R_{ϕ ϕ} (ξ)

is the autocorrelation function of

ϕ (x)

. Substituting Eqs. (74) and (77) into Eq. (73) yields the autocorrelation function

R_{h h} (ξ) = (2 N + 1) R_{ψ ψ} (ξ) + \sum_{k = - N}^{N} \sum_{\begin{matrix} l = - N \\ l \neq k \end{matrix}}^{N} R_{ϕ ϕ} (ξ + a k - a l) .

In above equation, the first term, associated with $R_{ψ ψ} (ξ)$ , characterizes the correlation of an individual speckle with itself; the second term, associated with $R_{ϕ ϕ} (ξ)$ , features the correlation of different speckles. To depart the influence of speckle positions, define

χ (ξ) = \sum_{k = - N}^{N} \sum_{l = - N}^{N} δ (ξ + a k - a l) - (2 N + 1) δ (ξ),

where

δ (ξ)

is the Dirac delta function. Thus, Eq. (78) can be rewritten as

R_{h h} (ξ) = (2 N + 1) R_{ψ ψ} (ξ) + R_{ϕ ϕ} (ξ) \otimes χ (ξ) .

B. Secondary peak height of 1D digital speckle pattern

The ratio of secondary peak height to the primary peak height can be approximated as

A_{2} \approx \frac{2 N}{1 + 2 N} \frac{R_{ϕ ϕ} (0)}{R_{ψ ψ} (0)} .

The main problem is to estimate the ratio $R_{ϕ ϕ} (0) / R_{ψ ψ} (0)$ . Since $ϕ (ξ) = ψ (ξ) \otimes p_{ϵ} (ξ)$ , the autocorrelation function of $ϕ (ξ)$ satisfies

R_{ϕ ϕ} (ξ) = R_{ψ ψ} (ξ) \otimes R_{p_{ϵ} p_{ϵ}} (ξ),

where

R_{p_{ϵ} p_{ϵ}} (ξ)

denotes the autocorrelation of

p_{ϵ} (ϵ)

. Hence,

R_{ϕ ϕ} (0) = \int_{- \infty}^{\infty} R_{ψ ψ} (ξ) R_{p_{ϵ} p_{ϵ}} (ξ) d ξ .

A very common situation in reality is that the speckle profile $ψ (x)$ can be characterized by a single parameter r and the probability distribution function $p_{ϵ} (ϵ)$ can be characterized by a single parameter c. In such a case, these two functions can be respectively normalized by theirs length scales:

ψ (x) = \bar{ψ} (x / r), p_{ϵ} (ϵ) = c^{- 1} {\bar{p}}_{ϵ} (ϵ / c) .

The factor $c^{- 1}$ ensures the integral of the probability distribution function to be unit. Accordingly, the autocorrelation function

R_{ψ ψ} (ξ) = r R_{\bar{ψ} \bar{ψ}} (ξ / r), R_{p_{ϵ} p_{ϵ}} (ξ) = c^{- 1} R_{{\bar{p}}_{ϵ} {\bar{p}}_{ϵ}} (ξ / c) .

A combination of Eq. (83) and Eq. (85) yields

\frac{R_{ϕ ϕ} (0)}{R_{ψ ψ} (0)} = \frac{\int_{- \infty}^{\infty} R_{\bar{ψ} \bar{ψ}} (ξ c / r) R_{{\bar{p}}_{ϵ} {\bar{p}}_{ϵ}} (ξ) d ξ}{R_{\bar{ψ} \bar{ψ}} (0)} .

The profile of a Gaussian speckle is $ψ (x) = \exp (- x^{2} / r^{2})$ . Hence,

\bar{ψ} (x) = exp (- x^{2}) .

The autocorrelation function of $\bar{ψ} (x)$ is

R_{\bar{ψ} \bar{ψ}} (ξ) = {(π / 2)}^{1 / 2} exp (- ξ^{2} / 2) .

Now let the random offset ϵ follow a uniform distribution over an interval $[- b, b]$ , then its probability distribution function is a rectangular function $p_{ϵ} (ξ) = {(2 b)}^{- 1} rect (ξ / 2 b)$ . Hence, the length scale parameter $c = 2 b$ and

{\bar{p}}_{ϵ} (ξ) = rect (ξ) .

The autocorrelation of a rectangular function is a triangle function, and thus

R_{{\bar{p}}_{ϵ} {\bar{p}}_{ϵ}} (ξ) = tri (ξ) .

In such a case,

\begin{matrix} \int_{- \infty}^{\infty} R_{\bar{ψ} \bar{ψ}} (ξ c / r) R_{{\bar{p}}_{c} {\bar{p}}_{c}_{ϵ}} (ξ) d ξ = \int_{- \infty}^{\infty} \sqrt{\frac{π}{2}} \exp (- \frac{2 b^{2} ξ^{2}}{r^{2}}) tri (ξ) d ξ \\ = \int_{0}^{1} \sqrt{2 π} (1 - ξ) \exp (- \frac{2 b^{2} ξ^{2}}{r^{2}}) d ξ \\ = \frac{π r}{2 b} erf (\frac{\sqrt{2} b}{r}) - \sqrt{\frac{π}{2}} \frac{r^{2}}{2 b^{2}} [1 - \exp (- \frac{2 b^{2}}{r^{2}})] . \end{matrix}

Hence, we obtain

Φ (k) = \frac{R_{ϕ ϕ} (0)}{R_{ψ ψ} (0)} = \sqrt{\frac{π}{2}} \frac{erf (\sqrt{2} k)}{k} - \frac{1 - e^{- 2 k^{2}}}{2 k^{2}},

where

k = b / r

is a dimensionless quantity that feature the randomness of the speckle pattern.

C. Secondary peak height of 2D digital speckle image

The secondary peak height of ZNCC-based autocorrealtion is

A_{2} \approx \frac{R_{ϕ ϕ} (0, 0) - a^{2} {\bar{h}}^{2}}{R_{ψ ψ} (0, 0) - a^{2} {\bar{h}}^{2}} .

The profile of an individual Gaussian speckle is

ψ (x, y) = exp (- \frac{x^{2} + y^{2}}{r^{2}}) .

The autocorrelation function of a Gaussian speckle is

R_{ψ ψ} (ξ, ζ) = \frac{π r^{2}}{2} exp (- \frac{ξ^{2} + ζ^{2}}{2 r^{2}}) .

The random displacements, both ϵ and τ, follow an uniform distribution in interval $[- b, b]$ . Thus, their joint probability distribution obeys

p_{ϵ, τ} (ϵ, τ) = \frac{1}{{(2 b)}^{2}} rect (\frac{ϵ}{2 b}) rect (\frac{τ}{2 b}) .

The autocorrelation function of $p_{ϵ, τ} (ϵ, τ)$ is

R_{p_{ϵ, τ} p_{ϵ, τ}} (ξ, ζ) = \frac{1}{{(2 b)}^{2}} tri (\frac{ξ}{2 b}) tri (\frac{ζ}{2 b}) .

Since $ϕ (ξ, ζ) = ψ (ξ, ζ) \otimes p_{ϵ, τ} (ξ, ζ)$ , the autocorrelation function of $ϕ (ξ, ζ)$ satisfies

R_{ϕ ϕ} (ξ, ζ) = R_{ψ ψ} (ξ, ζ) \otimes R_{p_{ϵ, τ} p_{ϵ, τ}} (ξ, ζ) .

Therefore,

\begin{matrix} R_{ϕ ϕ} (0, 0) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} R_{ψ ψ} (ξ, ζ) R_{p_{ϵ, τ} p_{ϵ, τ}} (ξ, ζ) d ξ d ζ \\ = \int_{0}^{2 b} \int_{0}^{2 b} \frac{π r^{2}}{2 b^{2}} \exp (- \frac{ξ^{2} + ζ^{2}}{2 r^{2}}) (1 - \frac{ξ}{2 b}) (1 - \frac{ζ}{2 b}) d ξ d ζ \\ = \frac{π r^{2}}{2} {\sqrt{\frac{π}{2}} \frac{r}{b} erf (\frac{\sqrt{2} b}{r}) - \frac{r^{2}}{2 b^{2}} [1 - \exp (- \frac{2 b^{2}}{r^{2}})]}^{2} . \end{matrix}

The average of $\bar{h}$ is

\bar{h} = a^{- 2} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} ψ (x, y) d x d y = \frac{π r^{2}}{a^{2}} .

Substituting Eqs. (95), (99), and (100) into Eq. (93) yields the ZNCC-based secondary peak height

A_{2} \approx \frac{Φ {(k)}^{2} - π ρ^{2} / 2}{1 - π ρ^{2} / 2},

where

ρ = 2 r / a

,

k = b / r

, and

Φ (k)

is given by Eq. (18). Indeed,

Φ (k)

provides an estimation of secondary autocorrelation peak height in the 1-dimensional case.

D. Power spectrum of 1D digital speckle pattern

The power spectrum of 1D digital speckle pattern $f (x)$ is

| \hat{f} (ω) |^{2} = | \hat{ψ} (ω) |^{2} | \hat{z} (ω) |^{2},

where

\hat{ψ} (ω)

denotes the spectrum of an individual speckle and

\hat{z} (ω) = \sum_{i = - N}^{N} exp (- j ω x_{i})

features the contribution of speckle position

x_{i}

. It follows from Eq. (103) that

\begin{matrix} | \hat{z} (ω) |^{2} = \sum_{m = - N}^{N} \sum_{n = - N}^{N} \exp [j ω (x_{n} - x_{m})] \\ = \sum_{m = - N}^{N} \sum_{n = - N}^{N} \exp [j ω a (n - m)] \exp [j ω (ϵ_{n} - ϵ_{m})] . \end{matrix}

Therefore, the expectation

E {| \hat{z} (ω) |^{2}} = \sum_{m = - N}^{N} \sum_{n = - N}^{N} β_{m, n} exp [j ω a (n - m)] .

where

β_{m, n} = E {exp [j ω (ϵ_{n} - ϵ_{m})]} .

The crucial issue is to identify the specific form of $β_{m, n}$ . If $m = n$ , it is clear that $β_{m, m} = 1$ . If $m \neq n$ , because of the mutually independence of imposed displacements,

β_{m, n} = E {exp (j ω ϵ_{n})} \cdot E {exp (- j ω ϵ_{m})} .

In probability theory,

E {exp (j ω x)} = \int_{- \infty}^{\infty} p_{x} (x) exp (j ω x) d x

is called characteristic function of random variable x [20]. Let

{\hat{p}}_{ϵ} (ω)

be the characteristic function of ϵ. As the probability distribution function of ϵ is real and even,

β_{m, n} = {\hat{p}}_{ϵ}^{2} (ω)

provided

m \neq n

. In summary,

β_{m, m} = {\begin{array}{l} 1, & m & = n; \\ {\hat{p}}_{ϵ}^{2} (ω), & m & \neq n . \end{array}

A combination of Eq. (105) and Eq. (109) yields

E {| \hat{z} (ω) |^{2}} = (2 N + 1) + {\hat{p}}_{ϵ}^{2} (ω) \sum_{m = - N}^{N} \sum_{\begin{matrix} n = - N \\ n \neq m \end{matrix}}^{N} exp [j ω a (n - m)] .

The second term in above equation

\begin{matrix} \sum_{m = - N}^{N} \sum_{\begin{matrix} n = - N \\ n \neq m \end{matrix}}^{N} \exp [j ω a (n - m)] = \sum_{m = - N}^{N} \exp (- j ω m a) \cdot \sum_{n = - N}^{N} \exp (j ω n a) - (2 N + 1) \\ = D_{N}^{2} (ω a) - (2 N + 1) . \end{matrix}

Here, $D_{N} (x) = \sum_{n = - N}^{N} \exp (j n x)$ denotes the N $^{t h}$ Dirichlet kernel whose closed form is [32]

D_{N} (x) = \frac{sin (N + \frac{1}{2}) x}{sin \frac{x}{2}} .

Hence, we can rewrite Eq. (110) as

E {| \hat{z} (ω) |^{2}} = (2 N + 1) + {\hat{p}}_{ϵ}^{2} (ω) [D_{N}^{2} (ω a) - (2 N + 1)] .

And thus the expectation of power spectrum can be derived.

E. Sum of square of subset gradients

Speckle pattern should be well sampled to avoid aliasing. Hence, for a pattern $h (x, y)$ , the sum of squared gradients of the whole pattern could be approximated as

{\tilde{Q}}_{h} = \sum_{m = - \infty}^{\infty} \sum_{n = - \infty}^{\infty} h_{x}^{2} (m, n) \approx \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} h_{x}^{2} (x, y) d ω_{x} d ω_{y} .

The gradient is generally estimated by convolving the pattern with a gradient filter $d [m, n]$ :

h_{x} (x, y) = \sum_{m = - \infty}^{\infty} \sum_{n = - \infty}^{\infty} d [m, n] h (x - m, y - n) .

With Parseval’s theorem, the integral in Eq. (114) could be transferred into the frequency domain:

{\tilde{Q}}_{h} \approx \frac{1}{4 π^{2}} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} | \hat{d} (ω_{x}, ω_{y}) |^{2} | \hat{h} (ω_{x}, ω_{y}) |^{2} d ω_{x} d ω_{y},

where

\hat{d} (ω_{x}, ω_{y})

is the magnitude response of gradient estimator

d [m, n]

and

\hat{h} (ω_{x}, ω_{y})

is the Fourier transform of

h (x, y)

. Then, for a digital speckle pattern,

E {{\tilde{Q}}_{h}} \approx \frac{1}{4 π^{2}} \int_{- π}^{π} \int_{- π}^{π} | \hat{d} (ω_{x}, ω_{y}) |^{2} | \hat{ψ} (ω_{x}, ω_{y}) |^{2} E {| \hat{z} (ω_{x}, ω_{y}) |^{2}} d ω_{x} d ω_{y},

where

E {| \hat{z} (ω_{x}, ω_{y}) |^{2}}

is given by Eq. (46).

In practical situations, the subset is a small part of the whole speckle area. The issue is to estimate the sum of square of intensity gradients within a subset:

Q_{h} = \sum_{m = - M}^{M} \sum_{n = - M}^{M} h_{x}^{2} (m, n),

where M is one-half of subset size and there are

{(2 M + 1)}^{2}

points in a subset. The digital speckle pattern are composed of

{(2 N + 1)}^{2}

speckles. The grid spacing is a, so that the area of each speckle could be seen as a². Therefore, it is believed that

\frac{E {Q_{h}}}{{(2 M + 1)}^{2}} \approx \frac{E {{\tilde{Q}}_{h}}}{{(2 N + 1)}^{2} a^{2}} .

Substituting Eq. (117) into above equation yields

\begin{matrix} E {Q_{h}} \approx \frac{{(2 M + 1)}^{2}}{4 π^{2} a^{2}} \int_{- π}^{π} \int_{- π}^{π} | \hat{d} (ω_{x}, ω_{y}) |^{2} | \hat{ψ} (ω_{x}, ω_{y}) |^{2} E {| \hat{z} (ω_{x}, ω_{y}) |^{2}} / {(2 N + 1)}^{2} d ω_{x} d ω_{y} \\ \approx \frac{{(2 M + 1)}^{2}}{4 π^{2} a^{2}} \int_{- π}^{π} \int_{- π}^{π} | \hat{d} (ω_{x}, ω_{y}) |^{2} | \hat{ψ} (ω_{x}, ω_{y}) |^{2} Θ (ω_{x}, ω_{y}) d ω_{x} d ω_{y} . \end{matrix}

Funding

National Natural Science Foundation of China (11702287, 11872354, 11627803); Strategic Priority Research Program of the Chinese Academy of Sciences (XDB22040502); Fundamental Research Funds for the Central Universities (WK2480000004).

References

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. C. Quan, C. J. Tay, W. Sun, and X. He, “Determination of three-dimensional displacement using two-dimensional digital image correlation,” Appl. Opt. 47(4), 583–593 (2008). [CrossRef] [PubMed]

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

5. X. Shao, M. M. Eisa, Z. Chen, S. Dong, and X. He, “Self-calibration single-lens 3D video extensometer for high-accuracy and real-time strain measurement,” Opt. Express 24(26), 30124–30138 (2016). [CrossRef]

6. X. Huang, Z. Liu, and H. Xie, “Recent progress in residual stress measurement techniques,” Acta Mech. Solida Sin. 26(6), 570–583 (2013). [CrossRef]

7. Z. Hu, T. Xu, H. Luo, R. Z. Gan, and H. Lu, “Measurement of thickness and profile of a transparent material using fluorescent stereo microscopy,” Opt. Express 24(26), 29822–29829 (2016). [CrossRef]

8. Z. Wang, H. Wu, K. Kang, S. Wang, Y. Li, W. Hou, H. Riaz, L. Li, and C. Li, “DIC/Moiré hybrid method based on regular patterns for deformation measurement,” Opt. Express 27(13), 18435–18444 (2019). [CrossRef] [PubMed]

9. M. Grédiac, B. Blaysat, and F. Sur, “A critical comparison of some metrological parameters characterizing local digital image correlation and grid method,” Exp. Mech. 57(6), 871–903 (2017). [CrossRef]

10. L. Wittevrongel, P. Lava, S. V. Lomov, and D. Debruyne, “A self adaptive global digital image correlation algorithm,” Exp. Mech. 55(2), 361–378 (2015). [CrossRef]

11. G. Stoilov, V. Kavardzhikov, and D. Pashkouleva, “A comparative study of random patterns for digital image correlation,” J. Theor. Appl. Mech. 42(2), 55–66 (2012). [CrossRef]

12. Z. Chen, C. Quan, F. Zhu, and X. He, “A method to transfer speckle patterns for digital image correlation,” Meas. Sci. Technol. 26(9), 095201 (2015). [CrossRef]

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

14. G. Bomarito, J. Hochhalter, T. Ruggles, and A. Cannon, “Increasing accuracy and precision of digital image correlation through pattern optimization,” Opt. Lasers Eng. 91, 73–85 (2017). [CrossRef]

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

16. Z. Chen, X. Shao, X. Xu, and X. He, “Optimized digital speckle patterns for digital image correlation by consideration of both accuracy and efficiency,” Appl. Opt. 57(4), 884–893 (2018). [CrossRef] [PubMed]

17. A. Lavatelli, R. Balcaen, E. Zappa, and D. Debruyne, “Closed-loop optimization of DIC speckle patterns based on simulated experiments,” IEEE Trans. Instrum. Meas. (to be published).

18. Y. Su, Q. Zhang, and Z. Gao, “Statistical model for speckle pattern optimization,” Opt. Express 25(24), 30259–30275 (2017). [CrossRef] [PubMed]

19. Y. Su, “Performance of digital speckle pattern,” figshare (2019) [retrieved 14 Jul 2019], http://dx.doi.org/10.6084/m9.figshare.8868221.v1.

20. A. Papoulis and S. U. Pillai, Probability, Random Variables, and Stochastic Processes (Tata McGraw-HillEducation, 2002).

21. Z. Hu, H. Xie, J. Lu, H. Wang, and J. Zhu, “Error evaluation technique for three-dimensional digital image correlation,” Appl. Opt. 50(33), 6239–6247 (2011). [CrossRef] [PubMed]

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

23. Y. Wang, P. Lava, P. Reu, and D. Debruyne, “Theoretical analysis on the measurement errors of local 2D DIC: part I temporal and spatial uncertainty quantification of displacement measurements,” Strain 52(2), 110–128 (2016). [CrossRef]

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

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

26. Y. Wang, M. Sutton, H. Bruck, and H. 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]

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

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

29. P. Thevenaz, T. Blu, and M. Unser, “Interpolation revisited [medical images application],” IEEE Trans. Med. Imaging 19(7), 739–758 (2000). [CrossRef] [PubMed]

30. S. Baker and I. Matthews, “Lucas-Kanade 20 years on: a unifying framework,” Int. J. Comput. Vis. 56(3), 221–255 (2004). [CrossRef]

31. Y. Su, Q. Zhang, X. Xu, Z. Gao, and S. Wu, “Interpolation bias for the inverse compositional Gauss-Newton algorithm in digital image correlation,” Opt. Lasers Eng. 100, 267–278 (2018). [CrossRef]

32. E. M. Stein and R. Shakarchi, Fourier Analysis: an Introduction (Princeton University Press, 2011).

33. X. Shao, X. Dai, and X. He, “Noise robustness and parallel computation of the inverse compositional Gauss-Newton algorithm in digital image correlation,” Opt. Lasers Eng. 71, 9–19 (2015). [CrossRef]

34. E. Jones and M. Iadicola, A Good Practices Guide for Digital Image Correlation (International Digital Image Correlation Society, 2018).

35. M. Grédiac and F. Sur, “Effect of sensor noise on the resolution and spatial resolution of displacement and strain maps estimated with the grid method,” Strain 50(1), 1–27 (2014). [CrossRef]

36. P. L. Reu, E. Toussaint, E. Jones, H. A. Bruck, M. Iadicola, R. Balcaen, D. Z. Turner, T. Siebert, P. Lava, and M. Simonsen, “DIC challenge: developing images and guidelines for evaluating accuracy and resolution of 2D analyses,” Exp. Mech. 58(7), 1067–1099 (2018). [CrossRef]

37. M. Grédiac, B. Blaysat, and F. Sur, “A robust-to-noise deconvolution algorithm to enhance displacement and strain maps obtained with local DIC and LSA,” Exp. Mech. 59(2), 219–243 (2019). [CrossRef]

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

39. A. Savitzky and M. J. E. Golay, “Smoothing and differentiation of data by simplified least squares procedures,” Anal. Chem. 36(8), 1627–1639 (1964). [CrossRef]

40. X. Xu, Y. Su, and Q. Zhang, “Theoretical estimation of systematic errors in local deformation measurements using digital image correlation,” Opt. Lasers Eng. 88, 265–279 (2017). [CrossRef]

41. B. Wang and B. Pan, “Self-adaptive digital volume correlation for unknown deformation fields,” Exp. Mech. 59(2), 149–162 (2019). [CrossRef]

Parameters
Image resolution	500 × 500 pixels
Bit depth	8 Bit
Minimum intensity $q_{min}$ [see Eq. (2)]	20
Maximum intensity $q_{max}$ [see Eq. (2)]	240

Parameters
Subpixel registration algorithm	IC-GN algorithm
Correlation criterion	sum of squared differences (SSD) criterion
Shape function	first-order shape function
Interpolation algorithm	cubic B-spline interpolation
Gradient estimator	Barron operator

Parameters
Region of interest	$[25, 475] \times [25, 475]$
Grid step	10 pixels
Subset size	31 pixels
Noise level of image	$σ_{f} = σ_{g} = 2.55$ (1% noise level)
Number of noise additions	100 times

Parameters
Image resolution	500 × 500 pixels
Bit depth	8 Bit
Minimum intensity $q_{min}$ [see Eq. (2)]	20
Maximum intensity $q_{max}$ [see Eq. (2)]	240

Parameters
Subpixel registration algorithm	IC-GN algorithm
Correlation criterion	sum of squared differences (SSD) criterion
Shape function	first-order shape function
Interpolation algorithm	cubic B-spline interpolation
Gradient estimator	Barron operator

Abstract

1. Introduction

2. Digital speckle pattern and digital speckle image

2.1. 2D digital speckle pattern and 2D digital speckle image

2.2. 2D Gaussian digital speckle pattern and non-overlapping condition

2.3. 1D digital speckle pattern and 1D digital speckle image

2.4. 1D Gaussian digital speckle pattern

3. Autocorrelation function

3.1. Autocorrelation function of 1D digital speckle pattern

3.2. Secondary autocorrelation peak height of 1D digital speckle pattern

3.3. Numerical simulation for 1D digital speckle pattern

3.4. Autocorrelation function of 2D digital speckle pattern

3.5. ZNCC-based autocorrelation function of 2D digital speckle image

3.6. Secondary autocorrelation peak height of 2D digital speckle image

3.7. Numerical simulation for 2D digital speckle image

4. Systematic errors of 1D digital speckle pattern

4.1. Theoretical analysis of 1D interpolation bias

4.2. Numerical simulation for interpolation bias

4.3. Influence of random offset on interpolation bias

4.4. Influence of speckle count on interpolation bias

5. Measurement errors of 2D digital speckle image

5.1. Theoretical analysis of 2D systematic error

5.2. Theoretical analysis of 2D random error

5.3. Numerical simulation for measurement errors

5.4. Influence of speckle radius on measurement errors

5.5. Influence of speckle variation on measurement errors

5.6. Influence of speckle coverage on measurement errors

6. Spatial resolution

6.1. Spatial resolution and Savitzky-Golay filter

6.2. Relationship between spatial resolution and subset size

6.3. A rudimentary model for spatial resolution

7. Discussion

7.1. Suggestions on selection of generation parameters of digital speckle pattern

7.2. Limits and outlook

8. Conclusion

Appendix

A. Autocorrelation function of 1D digital speckle pattern

B. Secondary peak height of 1D digital speckle pattern

C. Secondary peak height of 2D digital speckle image

D. Power spectrum of 1D digital speckle pattern

E. Sum of square of subset gradients

Funding

References

Cited By

Figures (22)

Tables (3)

Equations (120)

Optics Express