Noise robustness and parallel computation of the inverse compositional Gauss–Newton algorithm in digital image correlation

Optics and Lasers in Engineering 71 (2015) 9–19

Contents lists available at ScienceDirect

Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng

Noise robustness and parallel computation of the inverse compositional Gauss–Newton algorithm in digital image correlation Xinxing Shao, Xiangjun Dai, Xiaoyuan He n Department of Engineering Mechanics, Southeast University, Nanjing 210096, PR China

art ic l e i nf o

a b s t r a c t

Article history: Received 15 November 2014 Received in revised form 18 January 2015 Accepted 5 March 2015 Available online 31 March 2015

The inverse compositional Gauss–Newton (IC-GN) algorithm is one of the most popular sub-pixel registration algorithms in digital image correlation (DIC). The IC-GN algorithm, compared with the traditional forward additive Newton–Raphson (FA-NR) algorithm, can achieve the same accuracy in less time. However, there are no clear results regarding the noise robustness of IC-GN algorithm and the computational efficiency is still in need of further improvements. In this paper, a theoretical model of the IC-GN algorithm was derived based on the sum of squared differences correlation criterion and linear interpolation. The model indicates that the IC-GN algorithm has better noise robustness than the FA-NR algorithm, and shows no noise-induced bias if the gray gradient operator is chosen properly. Both numerical simulations and experiments show good agreements with the theoretical predictions. Furthermore, a seed point-based parallel method is proposed to improve the calculation speed. Compared with the recently proposed path-independent method, our model is feasible and practical, and it can maximize the computing speed using an improved initial guess. Moreover, we compared the computational efficiency of our method with that of the reliability-guided method using a four-point bending experiment, and the results show that the computational efficiency is greatly improved. This proposed parallel IC-GN algorithm has good noise robustness and is expected to be a practical option for real-time DIC. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Noise robustness Parallel computation Inverse compositional Gauss–Newton algorithm Digital image correlation

1. Introduction The technique of digital image correlation (DIC) [1–4] has been widely used for non-contact deformation measurement in the field of experimental mechanics. Improving the performance of DIC has been extensively investigated in different ways, including sub-pixel registration algorithms [5–7], sub-pixel interpolation schemes [8,9], subset size [10], subset shape functions [11,12], camera lens distortion [13,14], and image noise [15,16]. Among these aspects, the selection of the sub-pixel registration algorithm is vitally important. Pan et al. showed that the Newton–Raphson (NR) algorithm would be an optimal approach given its measurement accuracy and stability, if it were not for its slow calculation speed [7]. To speed the calculation of the NR algorithm, an inverse compositional Gauss–Newton (IC-GN) algorithm combined with a zero-mean normalized sum of squared difference (ZNSSD) correlation criterion was recently introduced for DIC [17]. Compared with the traditional forward additive Newton– Raphson (FA-NR) algorithm, the IC-GN algorithm can achieve the same accuracy in less time [17,18]. In terms of accuracy and stability, the IC-GN algorithm has been proven to be mathematically equivalent to the FA-NR algorithm [18]. In

n

Corresponding author. E-mail address: [email protected] (X. He).

http://dx.doi.org/10.1016/j.optlaseng.2015.03.005 0143-8166/& 2015 Elsevier Ltd. All rights reserved.

addition to accuracy and stability, robustness is an important standard for evaluating algorithms, especially the robustness in the presence of noise [19,20]. The noise robustness of the FA-NR algorithm has been thoroughly investigated by many researchers [15,16]. For example, Wang et al. assessed the quantitative error of the FA-NR algorithm caused by image noise [15]. However, few studies have systematically and quantitatively evaluated the noise robustness of the IC-GN algorithm. Furthermore, which algorithm has better noise robustness is an interesting but confusing question for the users of DIC. The computation efficiency of IC-GN is still in need of further improvements, although its calculation speed is faster than FA-NR [17,18]. As we know, DIC is a block-by-block processing technique, and parallel computing can help to calculate these blocks concurrently. Recently, a path-independent method combining a fast Fourier transform-based cross correlation (FFT-CC) algorithm and an IC-GN algorithm was proposed [21]. Parallel computing can be achieved based on this method, but the initial guess obtained by the FFT-CC algorithm is limited to translation, and an imperfect initial guess could lead to redundant iterations. In this paper, a novel one-dimensional (1D) theoretical model of the IC-GN algorithm is proposed. This model, based on the sum of squared differences (SSD) correlation criterion and linear interpolation, clearly indicates that the IC-GN algorithm has better noise robustness than FA-NR algorithm, and shows no noise-induced bias if the gray gradient operator is chosen properly. Both numerical

10

X. Shao et al. / Optics and Lasers in Engineering 71 (2015) 9–19

simulations and experiments were conducted to validate this, and the results show good agreement with the theoretical predictions. Moreover, a seed point-based parallel method is proposed to improve the calculation efficiency in Section 4. The implementation and advantages of the proposed method are discussed in detail. The accuracy and computational efficiency of the proposed parallel IC-GN algorithm are demonstrated by real experiments.

In IC-GN, the Hessian matrix need not be recalculated after precalculation because the reference subset remains the same [18]. Furthermore, there is another outstanding characteristic of IC-GN algorithm, that the first-order Taylor expansion of the grayscale intensities caused by the deformation increment is performed on the integer-pixel pointsWðx; y; 0Þ. Thus, the noise-induced bias error can be suppressed by choosing a proper gray gradient operator, as discussed in Section 3.

2. Brief introduction to the IC-GN algorithm 3. Noise robustness of IC-GN The IC-GN algorithm was first proposed by Baker and Matthews as an improvement of the Lucas–Kanade algorithm with respect to reducing computing time [18]. The algorithm was initially used as an image alignment method combined with the SSD correlation criterion. Tong [22] suggested using a noise-free reference image as well as additional reference images (generated by offsetting the reference image) in inverse update algorithm to improve its precision. Pan et al. [17] skillfully combined the IC-GN algorithm with the ZNSSD correlation criterion and used it as a sub-pixel DIC registration method. Gao et al. [23] analyzed the performance of IC-GN with first- order and second-order shape functions. Fig. 1 schematically shows the principles of the IC-GN and FANR algorithms, where Wðx; y; pÞ is the shape function that depicts the shape of the target subset relative to the reference subset, x, y are local coordinates, and p is the parameter vector. The practical and robust ZNSSD correlation criterion is used to evaluate the similarity between the reference and target subsets. In each iteration of the FA-NR algorithm, the incremental warp WðWðx; y; pÞ; ΔpÞ is exerted on the target subset, hence the parameter increment Δp can be added to the parameter vector directly (namely, pn þ 1 ¼ pn þ Δp), n is the number of iterations. However, the incremental warp Wðx; y; ΔpÞ is first exerted on the reference subset in each iteration of the IC-GN algorithm. Hence, it is subsequently inverted and composed with the deformation of the target subset. If a first-order shape function is used, the operation can be expressed as follows: 2 6 Wðx; y; pÞ’Wðx; y; pÞ 3 W  1 ðx; y; ΔpÞ ¼ 4

1 þ ux

uy

vx 0

1 þvy 0

u

32

76 v 54 1

1 þ Δux Δvx

Δuy 1 þ Δvy

0

0

3

Δu  1 Δv 7 5 1

ð1Þ where u; v are the displacement components and ux ; uy ; vx ; vy are the displacement gradients. This iteration is repeated until the pre-set convergence condition is satisfied. In this work, the convergence condition, ‖Δp‖ r 0:001, is predefined.

For the convenience of analysis, a 1D SSD correlation criterion is used to derive a theoretical model of the IC-GN and FA-NR algorithms. Consider a subset of 2Mþ1 pixels, the 1D SSD correlation function is defined as follows: ðuÞopt ¼ arg min

M X

½f ðxi Þ  gðxi þ uÞ2

ð2Þ

i ¼ M

where u is the deformation between the reference and target images, f ðxi Þ is the gray intensity value at point xi of the reference image, and gðxi þ uÞ is the gray intensity value at point xi þ u of the target image. To obtain the optimal value for u, non-linear iterative algorithms such as FA-NR, IC-GN, Steepest Descent, or Levenberg– Marquardt can be used. FA-NR and IC-GN are the mostly widely used in DIC. 3.1. IC-GN theoretical error model In IC-GN algorithm, we define u0 to be the initial value and u0 to be the measured translation after iteration. Thus, we have u0 ¼ u0  Δu

ð3Þ

where Δu is the deformation increment. Here we let u0 be the same with the exact translation, so the error of the measured translation due to iteration can be obtained. Minimization of the SSD correlation criterion function with respect to the variable Δu gives ! M P ∂ ½f ðxi þ ΔuÞ  gðxi þ u0 Þ2 i ¼ M

∂ðΔuÞ

¼0

ð4Þ

Here, xi is the coordinate of an integer-pixel point. Neglecting the high-order terms, the first-order Taylor expansion of f ðxi þ ΔuÞ

Fig. 1. Schematic of (a) IC-GN and (b) FA-NR.

X. Shao et al. / Optics and Lasers in Engineering 71 (2015) 9–19

yields f ðxi þ ΔuÞ ¼ f ðxi Þ þ Δuf x

ð5Þ

Then, the solution of Eq. (4) can be expressed in the following closed form M P

Δu ¼ i ¼  M

½gðxi þ u0 Þ f ðxi Þf x M P i ¼ M

ðf x Þ2

gðxi þ u0 Þ ¼ gðxi Þ þ u0 ½gðxi þ 1 Þ  gðxi Þ

ð7Þ

g x ¼ gðxi þ 1 Þ  gðxi Þ

ð8Þ

f x ¼ f ðxi þ 1 Þ  f ðxi Þ

ð9Þ

where g x and f x are the first-order derivatives of the grayscale intensities of the deformed and reference images, respectively. Various types of noise (e.g., thermal noise, readout noise, and shot noise) are unavoidably introduced into the images during recording. Obviously, the sub-pixel displacements calculated from Eqs. (7)–(9) will deviate from their actual values because of these types of noise. Thus, the image noise should be considered during the theoretical derivation. The grayscale intensity distribution of the noisy images can be written as follows: 0

f ðxi Þ ¼ f ðxi Þ þ ξ1 ðxi Þ

ð10Þ

0

f x ¼ f ðxi þ 1 Þ þ ξ1 ðxi þ 1 Þ  f ðxi Þ  ξ1 ðxi Þ

ð11Þ

g 0 ðxi Þ ¼ gðxi Þ þ ξ2 ðxi Þ

ð12Þ

g ðxi þ u0 Þ ¼ gðxi Þ þ ξ2 ðxi Þ þ u0 ½gðxi þ 1 Þ þ ξ2 ðxi þ 1 Þ  gðxi Þ  ξ2 ðxi Þ

ð13Þ

0

0

g x ¼ gðxi þ 1 Þ þ ξ2 ðxi þ 1 Þ gðxi Þ  ξ2 ðxi Þ

0

Δu ¼ i ¼  M

M P

0

M P

Δu ¼

ð19Þ

As shown in Eq. (17), the expectation value EðΔuÞ of the bias Δu has two terms. The first term only depends on the difference between the interpolated gray values and the true gray values. The second term is the bias error due to image noise. It is clear from Eq. (17) that the second term of bias is closely related to the noise of the reference image. Comparing Eq. (16) and Eq. (17) carefully, it can be further found that the impact of image noise on expectation is because Eðξ1 ðxi Þ2 Þ ¼ Varðξ1 ðxi ÞÞ þ Eðξ1 ðxi ÞÞ2 ¼ Varðξ1 ðxi ÞÞ ¼ σ 1 2

ð20Þ

The noise-induced bias error is introduced because the gray gradient operator contains f ðxi Þ, so it leads to the existence of ξ1 ðxi Þ in intensity gradients when considering image noise, as shown in Eq. (11). If a gray gradient operator is used that does not contain f ðxi Þ, such as the Barron gray operator [24]. The Barron gray operator is a four-point central differences operator and can be deduced from a cubic spline interpolation. the  1 In8 this8 research,  1 Barron gradient operator with a mask of 12 ;  12; 0; 12;  12 is used. 0 Then, f x can be rewritten as follows: 0

1 1 8 8 8 f ðx Þ þ ξ ðx Þ  f ðx Þ  ξ ðx Þ þ f ðx Þ 12 i  2 12 1 i  2 12 i  1 12 1 i  1 12 i þ 1 8 1 1 ð21Þ þ ξ1 ðxi þ 1 Þ  f ðxi þ 2 Þ  ξ1 ðxi þ 2 Þ 12 12 12

fx ¼

Substituting Eq. (10), Eq. (13) and Eq. (21) into Eq. (15), the expectation and variance of Δu can be obtained as follows: M P

EðΔ

uÞ ffi i ¼  M M P

½hðxi Þf x  ð22Þ

0

i ¼ M

ðf x Þ2

σ 21 þ σ 22

M P

ð23Þ

0

i ¼ M

ðf x Þ2

ðf x Þ2

Therefore, from Eq. (22) and Eq. (23), the expectation and variance of the measured translation u0 are obtained: M P

Eðu0 Þ ffi u0  i ¼

M M P

½hðxi Þf x  ð24Þ

0

i ¼ M

ðf x Þ2

½gðxi Þ þ ξ2 ðxi Þ þ u0 ðgðxi þ 1 Þ þ ξ2 ðxi þ 1 Þ  gðxi Þ  ξ2 ðxi ÞÞ  f ðxi Þ  ξ1 ðxi Þ½f ðxi þ 1 Þ þ ξ1 ðxi þ 1 Þ  f ðxi Þ  ξ1 ðxi Þ M P i ¼ M

M P M

½hðxi Þf x 

M P

i ¼ M

0

ðf x Þ2

ð16Þ

½f ðxi þ 1 Þ þ ξ1 ðxi þ 1 Þ f ðxi Þ  ξ1 ðxi Þ2

Applying the method described in Ref. [15], the expectation and variance of Δu can be written as follows:

EðΔuÞ ffi i ¼

hðxi Þ ¼ g 0 ðxi þu0 Þ  f ðxi Þ

ð15Þ

Note that this theoretical analysis is based on the following two assumptions. First, the image noise is independent at each pixel. Second, ξ1 and ξ2 are both white Gaussian noise with zero mean and standard deviations σ 1 and σ 2 , respectively. Substituting Eq. (10), Eq. (11) and Eq. (13) into Eq. (15), the expanded form can be obtained:

i ¼ M

ð18Þ

0

ðf x Þ2

where hðxi Þ is the gray interpolation error and it can be written as follows:

VarðΔuÞ ffi

0

i ¼ M

M P

i ¼ M

ð14Þ

where ξ1 and ξ2 denote the noise of the reference and target 0 0 images, respectively. Replacing f , f x and g in Eq. (6) withf , f x and 0 g , respectively, the following equation can be obtained: ½g 0 ðxi þ u0 Þ  f ðxi Þf x

σ 21 þ σ 22

VarðΔuÞ ffi

0

ð6Þ

After linear interpolation, we have

M P

11

Varðu0 Þ ffi

σ 21 þ σ 22

M P

i ¼ M

ð2M þ 1Þσ 21 þ M P 0 ðf x Þ2 i ¼ M

ð17Þ

0

; σ u0 ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Varðu0 Þ

ð25Þ

ðf x Þ2

Eq. (24) shows that the noise-induced bias of IC-GN is eliminated, although it is affected by the gray interpolation error and sum of square of subset intensity gradients (SSSIG). It can be seen

12

X. Shao et al. / Optics and Lasers in Engineering 71 (2015) 9–19

from Eq. (25) that there is no improvement in the displacement measurement accuracy of IC-GN; the variance remains unchanged. Further work to improve the accuracy of DIC from the perspective of this algorithm can be expected.

VarðΔuÞ ffi

u 0 ¼ u 0 þ Δu

¼0

ð27Þ

Here, xi þ u0 is the coordinate of a sub-pixel point. Neglecting the high-order terms, the first-order Taylor expansion of gðxi þ u0 þ ΔuÞ yields gðxi þ u0 þ ΔuÞ ¼ gðxi þ u0 Þ þ Δug x

ð28Þ

The solution of Eq. (27) can be then expressed as follows: M P

Δu ¼ i ¼  M

½f ðxi Þ  gðxi þ u0 Þg x ð29Þ

M P i ¼ M

ðg x Þ2 0

Δu ¼

0

i ¼ M

½f ðxi Þ  g 0 ðxi þ u0 Þg x

0

ð30Þ

M P i ¼ M

ðg x 0 Þ2

Substituting Eq. (10), Eq. (13) and Eq. (14) into Eq. (30), the expectation and variance of Δu can be written as follows [15]: M P

EðΔ

½  hðxi Þdg x 

uÞ ffi i ¼  M M P

i ¼ M

þf ðu0 Þ ðg x 0 Þ2

ð2M þ 1Þσ 22 M P ðg x 0 Þ2

ð32Þ

σ 21 þ σ 22

M P

i ¼ M

ð31Þ

i ¼ M

f ðu0 Þ ¼ 1  2u0 VarðΔuÞ ffi

8 1 1 ξ ðx Þ  gðx Þ  ξ ðx Þ 12 2 i þ 1 12 i þ 2 12 2 i þ 2

ð34Þ

Substituting Eq. (10), Eq. (13) and Eq. (34) into Eq. (30), the expectation and variance of Δu can be obtained as follows: M P

EðΔuÞ ffi i ¼

½  hðxi Þdg x 

M M P

i ¼ M

 ðg x 0 Þ2

Varðu0 Þ ffi

ðg x 0 Þ2

σ 21 þ σ 22

M P

i ¼ M

; σ u0 ¼

8 ð2M þ1Þσ 22 u0 M 12 P ðg x 0 Þ2

ð37Þ

i ¼ M

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V arðu0 Þ

ð38Þ

ðg x 0 Þ2

It can be seen from Eq. (37) that the noise-induced term still exists although the Barron gray operator was used. The noise-induced bias error of FA-NR algorithm cannot be eliminated in the same way 0 because of the existence of g 0 ðxi þ u0 Þg x as given in Eq. (30). The 0 0 intensity g ðxi þ u0 Þ and the gradient g x are estimated following Eq. (13) and Eq. (34), respectively. Inevitably, neighboring points contaminated by image noise are used to calculate the two parameters. Then, the effect of noise item, which has a quadratic form, cannot be avoided in the deformation increment. Hence, no matter which type of gray gradient operator is chosen to calculate derivatives of the grayscale intensity g x , the impact of noise on the expectation cannot be eliminated unless u0 is an inter-pixel movement.

In DIC, the bias error of measured displacement is a sum of two bias terms. The first term is the interpolation bias and it depends on the difference between the interpolated gray values and the true gray values. The second term is the noise-induced bias and it affects the precision a lot. Furthermore, it is worth mentioning that the first term (namely interpolation bias) is also affected by image noise. As shown in Section 3.1, the noise-induced bias of IC-GN can be eliminated. This is because the first-order Taylor expansion of grayscale intensities caused by the deformation increment is performed on integer-pixel points. Hence, the noise-induced bias can be eliminated by choosing a gray gradient operator that does not contain the factor f ðxi Þ. However, the bias error cannot be eliminated in the same way from the FA-NR algorithm because the first-order Taylor expansion of grayscale intensities is performed on sub-pixel points, as described in Section 3.2. Hence, we conclude that compared with FA-NR, IC-GN can achieve a faster calculation speed and better noise robustness with the same accuracy.

4. Seed point-based parallel method

1 1 8 8 8 gðx Þ þ ξ ðx Þ  gðx Þ  ξ ðx Þ þ gðx Þ 12 i  2 12 2 i  2 12 i  1 12 2 i  1 12 i þ 1

þ

i ¼ M



ð33Þ

ðg x 0 Þ2

Here, Barron gray operator is used to calculate intensity 0 gradients, sog x can be rewritten as follows: gx ' ¼

½hðxi Þdg x 

M M P

3.3. Noise robustness comparison between IC-GN and FA-NR 0

Replacing f , g, and g x in Eq. (28) with f , g 0 , and g x , respectively, the following equation can be obtained: M P

Eðu0 Þ ffi u0  i ¼

ð26Þ

Minimization of the SSD correlation criterion function with respect to the variable Δu gives ! M P ∂ ½f ðxi Þ  gðxi þ u0 þ ΔuÞ2

ð36Þ

ðg x 0 Þ2

Therefore, from Eq. (35) and Eq. (36), the expectation and variance of the measured translation u0 are obtained: M P

The iterative function of the FA-NR algorithm is defined as follows:

∂ðΔuÞ

M P

i ¼ M

3.2. FA-NR theoretical error model

i ¼ M

σ 21 þ σ 22

8 ð2M þ 1Þσ 22 u0 M 12 P ðg x 0 Þ2 i ¼ M

ð35Þ

Parallel computing can be used to improve the calculation speed of DIC. However, a parallel version of DIC has received very little attention, and few studies about it have been published. Many points need to be simultaneously calculated during a parallel implementation of DIC, hence the reliability-guided method [17] is not practical approach in this case. Jiang et al. [21] used the FFT-CC algorithm to estimate the initial guess for ICGN and proposed a path-independent method for parallel computation. Nevertheless, there are two weaknesses in the pathindependent method. First, the initial guess obtained by the FFTCC algorithm is limited to translation; hence it is not suitable for large rotations and deformations. Second, compared with the initial guess that is inherited from the neighboring points [25],

X. Shao et al. / Optics and Lasers in Engineering 71 (2015) 9–19

Fig. 2. Schematic of (a) improved initial guess transfer scheme and (b) seed point-based method.

Fig. 3. Flowchart of seed point-based parallel method.

13

14

X. Shao et al. / Optics and Lasers in Engineering 71 (2015) 9–19

the initial guess obtained by FFT-CC costs much more in redundant iterations. To solve these problems, a seed point-based parallel method is proposed in this paper. The deformation of the seed point is used to provide an improved initial guess for its neighboring points and every successfully calculated point is used as a new seed point for its neighboring points. The threshold of zero-mean normalized cross-correlation (ZNCC) coefficients can be used to distinguish the successfully calculated points. If the ZNCC coefficient of a computed point is more than the threshold value, the point is regarded as a successfully calculated point. The threshold of the ZNCC coefficient is set to 0.8 in this study. These unsuccessfully calculated points would not be considered as new seed points. Furthermore, the computed results of these points should not be used because the convergence and correctness cannot be verified by such low ZNCC coefficient. Fig. 2 shows a schematic of the improved initial guess transfer scheme and the seed point-based method. The initial guess of the neighboring points can be computed from the computed deformation parameter of the seed point [17].Once a seed point is determined, the neighboring four points can be calculated at the same time. The four calculated points are then used as new seed points, and another eight points can be computed. In this way, more and more points can be calculated simultaneously. A flowchart of the seed point-based parallel method is illustrated in Fig. 3. Before calculation, a queue Q and binary mask M c with the same size as the uncalculated points should be built. In addition, threads, which will be triggered during the calculation process, need

Fig. 4. Speckle pattern used in the numerical experiment.

to be created. During calculation, for a computed point (also known as a seed point), if its neighboring points have not been calculated, the neighboring points and the index of the seed point are inserted into Q. The index of the seed point can be used to retrieve its corresponding deformation from the deformation matrix and offer an improved initial guess. The binary mask M c is used to label the computed points and avoid recalculation of the same point. The initial value of each point in M c is set to “0,” and the value is changed to “1” if the point has been computed. The implementation of the proposed method is summarized as follows. Step 1: Specify a seed point in the region of interest (ROI) and calculate it using an automatic search scheme [26,27] and the IC-GN algorithm. The seed point is then marked as computed in the binary mask M c . Next, check the four neighboring points around the seed point. Uncalculated points that are labeled as not computed in M c and the index of the seed point are inserted into Q. Finally, trigger all threads. Step 2: The threads, which are used to retrieve an uncalculated point from Q and calculate it using the improved initial guess, can be truly concurrent. The improved initial guess [17,25] provides more newly computed points as seed points after a successful calculation. Next, analyze each of the neighboring four points around the seed point. If any of them have not been calculated, they are then inserted into Q along with the index of the seed point. Step 3: Repeat step 2 until Q is empty and all threads have finish calculating, indicating that the image correlation computation is complete.

From the above description, it is clear that points can be calculated in different threads simultaneously after the calculation of the first specified seed point. Furthermore, each computed point is used as a new seed point to calculate its neighboring points. In other words, the improved initial guess is for all points. Hence, the computing speed can be maximized with fewer iterations. Note that the initial guess of the first seed point can be calculated using an automatic search scheme so that the initial guess is not limited to translation. Compared to the path-independent method, the proposed method is more feasible and practical. In some critical cases, the ZNCC coefficients of most points are below the pre-set threshold. On this occasion, the existing reliability-guided method is recommended to be used if the results are really necessary.

Fig. 5. Mean bias error of measured u displacements for FA-NR (a) and IC-GN (b) using bicubic spline interpolation.

X. Shao et al. / Optics and Lasers in Engineering 71 (2015) 9–19

15

Fig. 6. Mean bias errors (a) and standard deviation errors (b) of measured u displacements for FA-NR and IC-GN at 0 noise level using bicubic spline interpolation.

Fig. 7. Mean bias errors (a) and standard deviation errors (b) of measured u displacements for FA-NR and IC-GN at 3% noise level using bicubic spline interpolation.

target subset relative to the reference subset in the following calculation.

5.1. Numerical simulations

Fig. 8. Experimental setup for cantilever beam bending.

5. Experimental verification To quantitatively compare the noise robustness of IC-GN to FANR, both numerical simulations and experiments were carried out in this study. The computational efficiency of the parallel IC-GN algorithm was demonstrated and compared to the reliability-guided IC-GN algorithm presented in [17]. All the algorithms were implemented using the Cþ þ language and tested on a laptop computer (Inter(R) Quad-Core (TM) i7-3720QM CPU with a main frequency of 2.6 GHz, 8.0 GB RAM). Because of Intel Hyper-Threading Technology, the Quad-Core CPU has eight logical threads in total. The First-order displacement shape function is used to depict the shape of the

Numerical analyses of six sets of synthetic images were carried out in this study. Fig. 4 shows the speckle pattern used in the following numerical experiments. The 8-bit speckle pattern with a resolution of 1280  1024 pixels was taken from our previous experiment. In the first image set, translated speckle images for each speckle pattern were generated by applying the appropriate shift in the Fourier domain according to the shift theorem [8]. The sub-pixel displacements applied in the x direction varied from 0 to 0.95 pixels, corresponding to a shift of 0.05 pixels between two successive images. The remaining five image sets were generated by adding various levels of white Gaussian noise to all the images in the first image set. The standard deviations of the added noise levels were 1%, 3%, 5%, 8% and 10% of the full 8-bit grayscale for image sets 2–6, respectively. The displacement of each translated speckle pattern was computed at 22,704 (¼ 172  132) regularly distributed points (the distance between neighboring points was five pixels) using 31  31 pixel subsets. To quantitatively evaluate the calculated displacements of the two different iterative algorithms, the error of the computed displacements were decomposed into two components: mean bias error and standard deviation error. The mean bias error of

16

X. Shao et al. / Optics and Lasers in Engineering 71 (2015) 9–19

Fig. 9. Calculation of the cantilever beam bending using FA-NR and IC-GN with a subset size of 51  51 pixels: (a) FA-NR v-displacement field, (b) IC-GN v-displacement field, (c) ZNCC coefficient distributions using FA-NR, and (d) ZNCC coefficient distributions using IC-GN.

Fig. 10. Real experimental images of a plastic specimen subjected to four-point bending: (a) reference image and (b) deformed image.

the measured displacement is defined as follows: Eu ¼ umean  uimp where umean ¼ N1

ð39Þ N P i¼1

ui represents the mean of the N estimated

displacements and uimp denotes the actual imposed sub-pixel displacement. The standard deviation error of the measured displacement can be defined as follows: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u 1 X σu ¼ t ðu  umean Þ ð40Þ N 1 i ¼ 1 i

Fig. 5 shows the mean bias error of the measured u displacement for FA-NR and IC-GN using the same bicubic spline interpolation, respectively. It can be concluded from Fig. 5 that: (1) bias error of the FA-NR algorithm increases with the increase of noise level. When the added noise level increases from 1% to 10%, the maximum bias error rises from 0.00522 to 0.148865 pixels. (2) Bias error of the IC-GN algorithm is only affected by the gray interpolation error and remains nearly unchanged while the image noise level increases. Furthermore, the maximum bias error is no more than 0.00679 pixels. Fig. 5 (b) shows that the gray interpolation error is also affected by image noise, but it is much smaller than the noise-induced bias error. When

X. Shao et al. / Optics and Lasers in Engineering 71 (2015) 9–19

using the same bicubic spline interpolation in both algorithms, the ICGN algorithm shows significantly less noise-induced bias. Figs. 6 and 7 show the mean bias and standard deviation errors of the measured u displacement for FA-NR and IC-GN at 0 and 3% noise level. For noiseless images, the standard deviation and bias error of both algorithms are close to each other. The negligible differences can be attributed to the different algorithm implementation of the two algorithms [17]. The result verifies that the two algorithms can achieve the same performance for noiseless images. When the standard deviation of the noise is 3%, the standard deviation of both IC-GN and FA-NR are similar while the bias error of IC-GN is much smaller than that of FA-NR. In summary, the noise robustness of IC-GN is better than FA-NR and IC-GN shows no noise-induced bias if the correct gray gradient operator is chosen.

17

aluminum specimen. The geometry of the specimen is shown in Fig. 8. The width is 78 mm, while the height and thickness are 15 and 2 mm, respectively. To increase the noise level of the images, a low-quality camera was used for image acquisition. Reference and target images with a size of 2048  2048 pixels and 8-bit gray

5.2. Experiments on real data To further compare the noise robustness of IC-GN with FA-NR, a real cantilever beam bending experiment was performed on an

Fig. 13. Comparison of the computational speed of parallel and reliability-guided IC-GN algorithms for various subset sizes.

Table 1 Comparison of the average number of iterations of parallel and reliability-guided IC-GN algorithms for various subset sizes. Subset size (pixels)

Fig. 11. Computational speed of parallel IC-GN for various numbers of created threads.

21 31 41 51

Average number of iterations Proposed parallel IC-GN

Reliability-guided IC-GN

2.68 2.11 1.99 1.94

2.68 2.11 1.99 1.93

Fig. 12. Calculation of the plastic specimen using parallel IC-GN with a subset of 31  31 pixels: (a) u-displacement field, (b) v-displacement field, and (c) ZNCC coefficient distributions.

18

X. Shao et al. / Optics and Lasers in Engineering 71 (2015) 9–19

levels were taken before and after deformation. During the DIC analysis of these images, a rectangular area was chosen to be the ROI, as shown in Fig. 8. Furthermore, the displacements were calculated at a mesh of 16,020 (267  60) points with a grid step of five pixels and subset size of 51  51 pixels. Fig. 9(a) and (b) presents the v-displacement field computed by FA-NR and IC-GN using the same bicubic spline interpolation, respectively. In addition, the ZNCC coefficient distributions are shown in Fig. 9(c) and (d). It is clear that the ZNCC coefficient of most points is larger than 0.96 for both algorithms, verifying their convergence. It can be seen from Fig. 9(a) and (b) that the vdisplacement field calculated by IC-GN is smoother than FA-NR. The distribution of Fig. 9(a) is similar to steps with a step width of nearly 1 pixel displacement. As the actual displacement field should be smooth, this result clearly shows the difference in noise robustness of the FA-NR and IC-GN algorithms. 5.3. Computational efficiency To verify the computational efficiency of the proposed parallel IC-GN implementation, a real four-point bending experiment was performed on a plastic specimen. As shown in Fig. 10, images with a size of 2048  2048 pixels at 8-bit gray levels were taken before and after deformation. A rectangular ROI containing 18,841 (83  227) points was calculated using the proposed parallel and reliability-guided IC-GN algorithms. It is quite important to note that to achieve a good calculation efficiency using the parallel IC-GN algorithm, the number of threads created before calculation should be considered. Fig. 11 shows the computational speed of the parallel IC-GN algorithm using various numbers of threads, with a range of 1–8. It can be seen from Fig. 11 that the computational speed of the parallel ICGN algorithm increases with the number of threads. Therefore, the computational speed is a linear function of the number of created threads. Because there are a total of eight threads in the computer used in the experiment, all eight threads were used in the following calculation. Displacement fields in the x and y directions computed by the parallel IC-GN algorithm are presented in Fig. 12(a) and (b) using a subset of 31  31 pixels and a grid step of 5 pixels. Fig. 12(c) shows the ZNCC coefficient distributions. It is clear that the ZNCC coefficients of most points are larger than 0.97, verifying the correctness of the calculated results. Although not shown here, we note that the reliability-guided IC-GN algorithm produces identical results. In terms of accuracy, there are no differences between the parallel implementation and reliability-guided method. Fig. 13 compares the computational efficiency of the parallel IC-GN algorithm with the recently proposed reliability-guided IC-GN algorithm. Compared with reliability-guided IC-GN, parallel IC-GN increases the computational speed by approximately six to seven times, depending on the subset size used. The gain in computational efficiency is mainly attributable to the parallel computation. If there are more than eight threads on the computer, the computational speed of parallel ICGN would be even faster. The results shown in Fig. 13 prove that a seed point-based parallel method can significantly increase the computational speed of DIC. In addition, Fig. 13 shows that the computational speed of both algorithms decreases steadily as the subset size increases. This can be explained by the fact that a larger subset requires more pixels to be analyzed. Table 1 further compares the average number of iterations of the parallel and reliability-guided IC-GN algorithms using various subset sizes, ranging from 21 pixels to 51pixels, for a fixed grid step 5pixels. The average number of iterations for both algorithms is almost the same and decreases with subset size. The results show that the improved initial guess was applied to all points and redundant iterations were avoided.

6. Conclusion DIC methods using an inverse compositional matching strategy, typically the reliability-guided IC-GN algorithm, have been widely used and commonly accepted for deformation analysis. In this study, the noise robustness of the IC-GN algorithm was investigated and compared with the FA-NR algorithm. Theoretical analysis and experimental results indicate that IC-GN has better noise robustness than FA-NR, and shows no noise-induced bias if the gray gradient operator is chosen properly. It is worth mentioning that gray interpolation error also suffers the impact of image noise in both IC-GN and FA-NR, but it is much smaller than the noiseinduced bias error. Additionally, a seed point-based parallel implementation of ICGN was proposed to improve the calculation speed of the method. An automatic search scheme was used to obtain the initial guess of the first seed point, and the improved initial guess was applied to all points via their neighboring points. Parallel DIC was implemented based on this proposed method. Experimental results show that the proposed parallel IC-GN algorithm is feasible, practical, and efficient. Using a computer with eight threads, the computation speed is approximately six to seven times faster than the existing reliability-guided IC-GN algorithm [17], depending on the subset size used. With the development of multi-thread technology, the calculation speed of the parallel DIC method could be further improved. Hence, the proposed parallel IC-GN algorithm is expected to be a practical choice for real-time DIC.

Acknowledgments This work is supported by the National Natural Science Foundation of China (under Grants 11272089 and 11327201).

References [1] Sutton MA, Orteu JJ, Schreier H. Image correlation for shape, motion and deformation measurements. New York: Springer; 2009. [2] Pan B, Qian K, Xie H, Asundi A. Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review. Meas Sci Technol 2009;20:062001. [3] Sutton MA, Wolters WJ, Peters WH, Ranson WF, McNeill SR. Determination of displacements using an improved digital correlation method. Image Vis Comput 1983;1:133–9. [4] Yamaguchi I. A laser-speckle strain gauge. J Phys E: Sci Instrum 1981;14:1270–3. [5] Davis CQ, Freeman DM. Statistics of subpixel registration algorithms based on spatiotemporal gradients or block matching. Opt Eng 1998;37:1290–8. [6] Bruck HA, McNeill SR, Sutton MA, Peters Iii WH. Digital image correlation using Newton–Raphson method of partial differential correction. Exp Mech 1989;29:261–7. [7] Pan B, Hui-Min X, Bo-Qin X, Fu-Long D. Performance of sub-pixel registration algorithms in digital image correlation. Meas Sci Technol 2006;17:1615–21. [8] Schreier HW, Braasch JR, Sutton MA. Systematic errors in digital image correlation caused by intensity interpolation. Opt Eng 2000;39:2915–21. [9] Luu L, Wang Z, Vo M, Hoang T, Ma J. Accuracy enhancement of digital image correlation with B-spline interpolation. Opt Lett 2011;36:3070–2. [10] Pan B, Xie H, Wang Z, Qian K, Wang Z. Study on subset size selection in digital image correlation for speckle patterns. Opt Express 2008;16:7037–48. [11] Schreier HW, Sutton MA. Systematic errors in digital image correlation due to undermatched subset shape functions. Exp Mech 2002;42:303–10. [12] Lu H, Cary PD. Deformation measurements by digital image correlation: implementation of a second-order displacement gradient. Exp Mech 2000;40:393–400. [13] Lava P, Van Paepegem W, Coppieters S, De Baere I, Wang Y, Debruyne D. Impact of lens distortions on strain measurements obtained with 2D digital image correlation. Opt Lasers Eng 2013;51:576–84. [14] Yoneyama S, Kitagawa A, Kitamura K, Kikuta H. Lens distortion correction for digital image correlation by measuring rigid body displacement. Opt Eng 2006;45:023602. [15] Wang YQ, Sutton MA, Bruck HA, Schreier HW. Quantitative error assessment in pattern matching: effects of intensity pattern noise, interpolation, strain and image contrast on motion measurements. Strain 2009;45:160–78.

X. Shao et al. / Optics and Lasers in Engineering 71 (2015) 9–19

[16] Wang ZY, Li HQ, Tong JW, Ruan JT. Statistical analysis of the effect of intensity pattern noise on the displacement measurement precision of digital image correlation using self-correlated images. Exp Mech 2007;47:701–7. [17] Pan B, Li K, Tong W. Fast, robust and accurate digital image correlation calculation without redundant computations. Exp Mech 2013;53:1277–89. [18] Baker S, Matthews I. Lucas-kanade 20 years on: a unifying framework. International journal of computer vision 2004;56:221–55. [19] Réfrégier Ph. Optimal trade-off filters for noise robustness, sharpness of the correlation peak, and Horner efficiency. Opt Lett 1991;16:829–31. [20] Kim Jong Min, Lee Ok Kyun, Ye Jong Chul. Improving noise robustness in subspacebased joint sparse recovery. IEEE Trans Signal Process 2012;60:5799–809. [21] Jiang Z, Kemao Q, Miao H, Yang J, Tang L. Path-independent digital image correlation with high accuracy, speed and robustness. Opt Lasers Eng 2015;65:93–102.

19

[22] Tong W. Subpixel image registration with reduced bias. Opt Lett 2011;36:763–5. [23] Gao Y, Cheng T, Su Y, Xu X, Zhang Y, Zhang Q. High-efficiency and highaccuracy digital image correlation for three-dimensional measurement. Opt Lasers Eng 2015;65:73–80. [24] Barron JL, Fleet DJ, Beauchemin SS. Systems and experiment performance of optical flow techniques. Int J Comput Vis 1994;12:43–77. [25] Zhou Y, Chen YQ. Propagation function for accurate initialization and efficiency enhancement of digital image correlation. Opt Lasers Eng 2012;50:1789–97. [26] Zhao JQ, Zeng P, Lei LP, Ma Y. Initial guess by improved population-based intelligent algorithms for large inter-frame deformation measurement using digital image correlation. Opt Lasers Eng 2012;50:473–90. [27] Zhou Y, Pan B, Chen YQ. Large deformation measurement using digital image correlation: a fully automated approach. Appl Opt 2012;51:7674–83.