1. Introduction

DIC (Digital Image Correlation) ^[1,2,3] is increasingly being used in engineering fields and other fields such as medical image processing ^[4] for measuring displacements and deformations. The DIC measurements are basically provided by tracking the displacement of pixels by recording sequential images of an object before and after displacement using a pixel-wise search scheme. This process is computationally demanding and the probability of algorithmic mismatch or divergence is relatively high, especially for large deformations ^[5]. The accurate and efficient processing of the images therefore represents a challenge for the DIC algorithm.

More specifically, the DIC algorithm includes two computational steps: $(1)$ a coarse initial guess (IG) to find an approximation of the true displacements and $(2)$ a Newton-Raphson (NR) iteration process to successively find the accurate displacements and often their gradients ^[1,6,7].

Significant efforts were made in the past few decades to improve both computational processes in terms of the accuracy and computational efficiency. Earlier studies on subset size ^[8], interpolation scheme ^[9,10], shape functions ^[10], correlation criteria ^[11], speckle pattern ^[12] and convergence criteria ^[7] have considerably advanced the DIC. Experimental aspects such as the object surface ^[13], chemical etching of object surface ^[14], hand drawing ^[14], paint spraying ^[8,13,15], laser beam structuring ^[16] and even speckle projections ^[14] or specked uniqueness ^[17] were also systemically studied and led to important technological improvements. Other more recent studies successfully deployed advanced algorithms such seed-point initiation, feature detection or registration tracking techniques ^[5,17,18,19]. The latter studies used special function and transform (affine, Fourier-Mellin) for a more sophisticated tracking particularly for the objects undergoing large deformations and motions ^[20]. Zhang et al. (2015) ran a parallel DIC code on a graphics processing unit (GPU) for its parallel computing technology and successfully improved the algorithmic efficiency ^[21]. The DIC computational cost remains high and efficient solutions without compromising the accuracy are desired. In the present article, the focus is made on making the DIC algorithm more efficient. A new efficient algorithmic method is also effective in finding IG values that is closer to the true displacements.

The initial guess of the displacement is used to initiate the NR iteration and as such essential to the overall DIC convergence and efficiency. More precisely, the convergence radius $R = \max (\sqrt{(u_0-u)^2+(v_0-v)^2})$ ($u$, $v$) must be within 3 pixels from the true value for a high convergence success rate ^[7]. Here, $u_0$, $v_0$ represent the initial guess and $u$, $v$ represent the true values of displacement components. The accuracy of initial guess can be improved by increasing the subset size. However, the corresponding computational cost increases remarkably when increasing the subset size.

To obtain an enhanced convergence of NR iteration of DIC measurement with a more accurate IG approximation which is needed in our laboratory testing of porous acrylic cement specimens whose surface is porous and affect the mechanical testing ^[22], we used the Fuzzy Distance Transform (FDT) that was a method developed and routinely used in our laboratory in image treatment. FDT ^[23,24,25] was shown to be a particularly effective method to treat and study features in low-resolution micro-computed tomography images. As detailed below, FDT method is a different class of a transform that takes into account the grey value the neighboring pixels to more accurately approximate the distance between two pixels and as such to rebuild an object in the digital image. Therefore, it is hypothesized that FDT can help provide accurate IG values for the DIC calculation. In the present study, the efficiency and accuracy of a novel accelerated FDT method was accordingly introduced and verified using numerical experiments, with a wide range of motion conditions due to the noise-free advantage of numerically generated images.

---

2. Materials and method

In this study, the algorithm consists of three main steps: (1) accelerated integer-pixel level initial guess; (2) more accurate initial guess at sub-pixel level and (3) Newton-Raphson iteration for accurate deformation. Particularly, Step 1 is to calculate the integer-level displacement with an accelerated scheme and Step 2 is to increase the accuracy of initial guess to sub-pixel level in aim of reducing the computational cost in initial guess process and improving the accuracy of initial guess which potentially guarantees the convergence of NR iteration (Step 3) in complicated deformation conditions in experimental environment.

---

2.1. Initial guess algorithm

The efficiency of the new initial guess algorithm stems from reducing the computationally-demanding search to a smaller subset of the entire images. Specifically, the new efficient process proposed to accelerate the calculation of the initial guess for the displacements of all selected pixel of interest (POI) is broken into two steps: (1) an initial full-field search for the displacements of the first POI in the entire image after displacement and (2) a reduced search for the displacements of the remaining POIs using the displacements of the first POI, determined according to Step 1.

---

2.1.1. Full-field search for displacement of first POI

Figure 1 illustrated the scheme of the full-field search for the initial guess of $1^{st}$ POI. In the reference image, a grid of POIs was selected with a defined grid step (GS). The $1^{st}$ POI was located at $(x_{10}, y_{10})$, and a square of pixels $(x_{1i}, y_{1j})$ around $1^{st}$ POI was chosen to be the reference subset (size of $(2M+1)^2$, denoted as Ref. subset in 1). Thereafter, a same size of subset in the deformed image was also selected as the deformed subset (denoted as Def. subset in 1). The center and other pixels inside the deformed subset were written as $(x_{10}^{^{\prime }}, y_{10}^{\prime })$ and $(x_{1i}^{^{\prime }}, y_{1j}^{\prime })$, respectively. The intensity of two subsets were compared with the robust Zero-Normalized Cross Correlation (ZNCC) criteria^[12] as Eq 1 presented

$\begin{equation} C_{ZNCC} = \sum\limits_{i, j = -M}^{M}\frac{\left[ f( x_{i}, y_{j})-f_{m}\right] \times \left[ g(x_{i}^{\prime }, y_{j}^{\prime })-g_{m}\right] }{\Delta f\times \Delta g}. \label{c-zncc} \end{equation}$

(1)

where $f(x_{i}, y_{j} )$ and $g(x'_{i}, y'_{j})$ were the intensity (gray scale value) of each pixel inside the reference and deformed subsets, respectively; $f_{m}$ and $g_{m}$ were the averaged intensity of two subsets; and equations of $\Delta f$, $\Delta g$ were presented in Eq 2.

$\begin{equation}\label{deltaf_fm_gm} \begin{aligned} \Delta f & = &\sqrt{\sum\limits_{i, j = -M}^{M}\left[ {{f({{x}_{i}}, {{y}_{j})}-{{f}_{m}}} }\right] {^{2}}}, \\ \Delta g & = &\sqrt{\sum\limits_{i, j = -M}^{M}\left[ {{g(x_{i}^{^{\prime }}, y_{j}^{\prime })-{{g}_{m}}}}\right] {^{2}}} \end{aligned} \end{equation}$

(2)

By systematically shifting the location of $(x_{10}^{^{\prime}}, y_{10}^{\prime })$ over the displaced image, a matrix $M_{C}$ of correlation coefficients was obtained by correlating the reference subset with deformed subsets at all locations that $(x_{10}^{^{\prime }}, y_{10}^{\prime })$ has been shifted. As the table in Figure 1 shows, the correlation coefficient matrix $M_{C}$ had a size of $m\times n$ pixels which is approximately equal to the deformed image size. Each component of $M_{C}$ had a range of $[-1, 1]$ where the value 1 indicated the optimum, that is best match between images taken before and after deformation. By finding the maximum value in coefficient matrix $M_{C}$, the optimum location of $(x_{10}^{^{\prime }}, y_{10}^{\prime })$ where the deformed subset had the optimum correlation coefficient with the reference subset of $1^{st}$ POI was accordingly identified. The initial guess of the displacement components $(u_{1}, v_{1})$ could thus be obtained from Eq 3. The initial guess from this step was at integer-pixel level.

$\begin{equation} \begin{aligned}\label{dispU_V} u_1 & = & x'_{10}-x_{10}, \\ v_1 & = & y'_{10}-y_{10}. \end{aligned} \end{equation}$

(3)

---

2.1.2. Reduced IG search of the remaining POIs

From the previous step, the initial values of $1^{st}$ initial displacement guess was obtained by the full-field search algorithm. Yet for the $2^{nd}$ and rest of selected POIs, an accelerated search process was developed and used to estimate their respective locations using the coordinates of $2^{nd}$ POI and, the already determined $1^{st}$ POI initial guess. Therefore, the follow-up search is done in a considerably reduced area around the estimated locations instead of performing the computationally costly full-field search for $2^{nd}$ POI.

Take as per example the $2^{nd}$ POI whose coordinates were $(x_{20}, y_{20})$ and whose reference subset pixels were $(x_{2i}, y_{2j})$, the center $(x_{2c}^{^{\prime }}, y_{2c}^{\prime })$ of its estimated locations was determined using

$\begin{equation} \begin{aligned}\label{cx_cy_2poi} x'_{2c}& = & x_{20}+u_{1}, \\ y'_{2c}& = & y_{20}+v_{1}. \end{aligned} \end{equation}$

(4)

Then the reduced search zone was $[(x_{2c}^{\prime}-s):(x_{2c}^{\prime}+s), (y_{2c}^{\prime}-s):(y_{2c}^{\prime}+s)]$ with an offset of $s$ pixels in both $x$ and $y$ directions as Figure 2 illustrated. By shifting the location of deformed subset (in the reduced search zone), a reduced matrix $M_{C}$ of correlation coefficient was then obtained. Using the same scheme of calculating the deformed location of $1^{st}$ POI, the location at which the maximum correlation coefficient was found and was then considered as the new location $(x_{20}^{\prime}, y_{20}^{\prime})$ of POI $(x_{20}, y_{20})$. The initial guess $(u_{2}, v_{2})$ of this POI could be then obtained using Eq 3 by replacing the coordinate of $1^{st}$ POI with that of $2^{nd}$ POI.

Using this second step of the new process, the initial guess for displacements of all selected POIs except $1^{st}$ could be determined with a considerably reduced computational cost.

---

2.2. Fuzzy-logic based initial guess with sub-pixel resolution

It should be however noted that the displacement IG values resulting from the aforementioned accelerated/reduced search scheme are integer without a fractional part because of the discrete nature of the digital images. The accuracy of the tracking is therefore incremental by one pixel. The new fuzzy-logic based process particularly provided IG values of fractional nature (real numbers) at sub-pixel level by calculating a fractional increment over the integer initial guess values. Since this fuzzy-logic based process applies to all POIs, the $1^{st}$ POI $(x_{10}, y_{10})$ was still taken as an example to explain the structure of the novel fuzzy-based process.

To calculate the displacement increment, there were five steps to follow:

1. Identify the pixel of the optimum location $(x_{10}^{\prime}, y_{10}^{\prime})$ of the POI $(x_{10}, y_{10})$ and its correlation coefficient matrix $M_{C}$.

2. Identify the coordinate and correlation coefficient of 8 pixels closest to the optimum or index pixel $(x_{10}^{\prime}, y_{10}^{\prime})$. There were 9 pixels in total including the optimum $(x_{10}^{\prime}, y_{10}^{\prime})$. The coordinate and correlation coefficients of 8 closest pixels were labeled as $(x_{m}^{\prime}, y_{m}^{\prime})$ and $c_{m}$ $(m = 1, \cdots, 8)$ as Figure 3 illustrated.

3. Calculate, utilizing the fractional correlation values of the adjacent pixel of $M_{C}$, a weighting value $w_{m}$ $(m = 1, \cdots, 8)$ for each pixel according to

$\begin{equation} w_{m} = \frac{c_{m}-c_{\min }}{\sum_{m = 1}^{8}(c_{m}-c_{\min })}, \label{weight_fuzzy} \end{equation}$

(5)

where $c_{m}$ $(m = 1, \ldots, 8)$ was the correlation coefficient of $m^{th}$ pixel and $c_{min}$ was the minimum coefficient of the 8 pixels. Please note that the reason why Eq 5 was established in this form was because the correlation value of the index/optimum pixel in the $ZNCC$ correlation criteria was 1. The higher the coefficient, the higher the similarity between reference and displaced subset. Since $(x_{10}^{\prime}, y_{10}^{\prime})$ was the location with the biggest coefficient, its neighboring pixels whose coefficient was closest to the coefficient of $(x_{10}^{\prime}, y_{10}^{\prime})$ obtained biggest weight value.

4. Calculate the fractional displacement increment using the weighting value $w_{m}$ of the nine pixels and of the distance between an individual pixel and the identified optimum/index location $(x_{10}^{\prime}, y_{10}^{\prime})$. The increment $(\Delta u_{m}, \Delta v_{m})$ was calculated according to

$\begin{equation} \begin{aligned} \label{deltaU_V_fuzzy} \Delta u_m & = & w_m \times (x'_m-x'_{10}), \\ \Delta v_m & = & w_m \times (y'_m-y'_{10}). \end{aligned} \end{equation}$

(6)

where the pixel with the biggest weighting value contributed most to the total displacement increment. Moreover, $\Delta u_{m}$ and $\Delta v_{m}$ were negative when $x_{m}^{\prime}<x_{10}^{\prime}$ and $y_{m}^{\prime}<y_{10}^{\prime}$, respectively. Likewise, $\Delta u_{m}$ and $\Delta v_{m}$ were positive when $x_{m}^{\prime}> x_{1}^{\prime}$ and $y_{m}^{\prime}>y_{1}^{\prime}$, respectively. In brief, Eq 6 resulted in shifting the modified $(x'_{10}, y'_{10})$ towards the neighboring pixels with the higher correlation coefficients. The fractional increment was then added to the integer initial guess $(u_{1}, v_{1})$ for displacement.

5. Sum up the increments $(\Delta u_m, \Delta v_m)$ to obtain the total displacement increment from the weighting process,

$\begin{equation} \begin{aligned}\label{u_v_fuzzy} u_{fuzzy} & = &u_{1}+\sum\limits_{m = 1}^{8}\Delta u_{m}, \\ v_{fuzzy} & = &v_{1}+\sum\limits_{m = 1}^{8}\Delta v_{m}. \end{aligned} \end{equation}$

(7)

Through the above five steps, sub-pixel/fractional level increments of the displacements were obtained and can be added to the initial integer guess for the purpose of more accurate initial input in the subsequent NR iteration process. The more accurate initial guess $(u_{fuzzy}, v_{fuzzy})$ was expected to increase the computational speed of Newton-Raphson iteration process and enhance its convergence.

---

2.3. Newton-Raphson process for the accurate displacements and deformations

Once an accurate initial guess of displacement of all selected POI was obtained, the accurate displacement and its gradients could be computed using the NR iteration. If we defined the displacement and gradient components in the deformation vector ${p}$ as ${p} = \left[u;v;\partial u/\partial x;\partial u/\partial y;\partial v/\partial x;\partial v/\partial y\right]$, the initial guess obtained from fuzzy/weighting process of any POI $(x_0, y_0 )$ was identified as ${p}_n$ with $n$ was the iteration number. Then the next deformation vector ${p}_{n+1}$ could be obtained by

$\begin{equation} p_{n+1}-p_{n} = -\frac{\nabla C\left( {{p}_{n}}\right) }{\nabla ^{2}C\left( {{p }_{n}}\right) }, \label{nrIter} \end{equation}$

(8)

where $\nabla C(p_n )$ and $\nabla^2 C(p_n )$ were the first-order and second-order derivatives of the correlation function, which are the Jacobian or Hessian matrix. In accordance with previous studies ^{[1,2,7,10,11,26,27]}, a robust and fast Zero-Normalized Sum of Squared Differences (ZNSSD) correlation criteria was chosen to correlate the reference and deformed subset in NR iteration process,

$\begin{equation} C({p}) = \sum\limits_{i, j = -M}^{M}\, \left[ \frac{f({{x}_{i}}, {{y}_{j})}-f_{m}}{ \Delta f}-\frac{g(x'_{i}, y'_{j})-g_{m}}{\Delta g}\right]^{2} . \label{cpartPn} \end{equation}$

(9)

In Eq 9, the item $\frac{f({{x}_{i}}, {{y}_{j}})-f_{m}}{\Delta f}$ carried the information of the reference subset which was known, yet the item $\frac{g(x'_{i}, y'_{j})-g_{m}}{\Delta g}$ was unknown since the pixel $(x'_i, y'_j )$ inside the deformed subset may deform to a location between integer pixels in NR iteration process. The intensity $g(x'_i, y'_j)$ of all pixels inside deformed subset was not available in the digital image, thus the bicubic spline interpolation was used to interpolate the intensity of pixels at any location during whole NR iteration process,

$\begin{equation} g\left( x'_i, y'_j\right) = \sum\limits_{l, k = 0}^{3}{{a}_{lk}}\cdot (\delta x)^{l}\cdot (\delta y)^{k}, \label{basicbicub} \end{equation}$

(10)

where $\delta x$ and $\delta y$ were the distances between pixel $(x'_i, y'_j)$ and the starting border of interpolation interval where $(x'_i, y'_j)$ was located. $a_{lk}$ was the interpolation coefficient within the interpolation. There were $16$ of $a_{lk}$ for each interpolation interval, and they were computed in advance for the sake of low computational cost.

In fact, the coordinate of pixels $(x'_i, y'_j)$ in the deformed subset was related to the coordinate of pixels $(x_i, y_j)$ inside the reference subset by the shape function:

$\begin{equation} \begin{aligned}\label{1st_shapefunx_y} {x_{i}}^{\prime } = x_{i}+u+\frac{\partial u}{\partial x}(x_i-x_0)+\frac{\partial u}{ \partial y} (y_i-y_0), \\ {y_{j}}^{\prime } = y_{j}+v+\frac{\partial v}{\partial x}(x_i-x_0)+\frac{\partial v}{ \partial y} (y_i-y_0), \end{aligned} \end{equation}$

(11)

where only the six variables $u, v, \partial u/\partial x, \partial u/\partial y, \partial v/\partial x, \partial v/\partial y$ in the vector ${p}$ were unknown. Therefore, the correlation coefficient $C$ in Eq 9 becomes the function of vector ${p}$ only. To calculate the first-order and second-order derivatives of $C$ by the vector ${p}$, a simplification was made. We denoted each variable in the vector ${p}$ as $p_{t, r} = u, v, \cdots, \partial v/\partial y, (t, r = 1, \cdots, 6)$. The Jacobian matrices of the deformation vector ${p}$ thus was written as

$\begin{equation} \nabla C\left( {{p}}\right) = \left( \begin{array}{cccccc} \frac{\partial C}{\partial p_1}&\cdots&\frac{\partial C}{\partial p_t}&\cdots&\frac{\partial C}{\partial p_6} \end{array} \right) \label{jacob} \end{equation}$

(12)

$\begin{eqnarray} \nabla ^{2}C\left( p\right) = \left( \begin{array}{ccc} \frac{\partial ^{2}C}{\partial p_{1}\partial p_{1}}&\ldots &\frac{ \partial ^{2}C}{\partial p_{1}\partial p_{6}} \\ \vdots &\frac{\partial ^{2}C}{\partial p_{t}\partial p_{r}}&\vdots \\ \frac{\partial ^{2}C}{\partial p_{6}\partial p_{1}}&\ldots &\frac{ \partial ^{2}C}{\partial p_{6}\partial p_{6}} \end{array} \right) \label{hessn} \end{eqnarray}$

(13)

The optimization ^[27] was used to simplify the Jacobian matrix with a good accuracy ^[28,29]. Namely, the item $(f(x_i, y_j )-f_m)/\Delta f-(g(x'_i, y'_j )-g_m)/\Delta g \approx 0$ when ${p}$ was close enough to the actual value. Then

$\begin{equation} \frac{\partial^2 C}{\partial p_t \partial p_r } \approx 2\sum\limits_{i, j = -M}^{M} \prod\limits_{k = t, r} \frac{ \partial }{\partial p_{k}}\left[ \frac{g\left( x_{i}^{\prime }, y_{j}^{\prime }\right) -g_m }{\Delta g}\right] . \label{hessn_optimized} \end{equation}$

(14)

By substituting Eq 2, Eq 10 and Eq 11 into Eq 14, the Jacobian and Hessian matrices were then calculated. Eq 8 then iterated until the convergence criteria was met, and the optimum unknown vector ${p}$ were thus iteratively found. In this study, the convergence criteria were set stringently. Only if both of Eq 15 was met, then tell algorithm to output the vector ${p}$ to guarantee the success rate of convergence.

$\begin{equation} \begin{aligned}\label{converg1_2} \left\Vert \Delta p\right\Vert & = &\left\Vert p_{n+1}-p_{n}\right\Vert = \sqrt{\sum\limits_{t = 1}^{6}\left( \Delta p_{t}\right) ^{2}}<10^{-3}, \\ \frac{\left\vert \Delta C\right\vert }{\left\Vert \Delta p\right\Vert } &<&10^{-8}. \end{aligned} \end{equation}$

(15)

Comparing to the conventional algorithm which consists of two steps (integer initial guess and Newton-Raphson iteration), a novel algorithm was formed with two extra steps in initial guess calculation and thereafter referred as accelerated F-NR algorithm. Both of the acceleration and fuzzy processes as well as F-NR algorithm were verified as described in the following section.

---

3. Results and discussion

Since Digital Image Correlation is an image-based and light-intensive method to measure the deformation occurs in experimental environment, there is diversity of external error sources reported (e.g. lightening errors, lens distortion errors, camera self-heating errors, etc.) ^{[11,27,28,29,30]}. This study is mainly focused on verifying the new algorithm itself in terms of accuracy and computational efficiency regardless of external errors at this moment. Therefore, the numerical images with known deformation and free of various errors is chosen to verify the algorithm.

An algorithm for numerical image generation was developed according to Zhou ^[31] and then used to generate a series of image pairs with a prior known displacement and deformations. Based on the preliminary examination of image histogram-related pattern tests, the image size, speckle size and speckle number were chosen to be $256\times 256$ pixels, 4 pixels and 150000, respectively.

---

3.1. Analysis of the acceleration and fuzzy processes

Figure 4 showed an example of a pair of images generated with a rotation angle of $2^{\circ }$ counter-clockwise. To investigate the impact of the acceleration and fuzzy processes, three different versions of DIC algorithm were used to process the example image pair.

Specifically, Version 1 is the conventional algorithm which consists of integer initial guess (IG) search scheme which applies to all POIs as described in Section 2.1.1 and the NR process as described in Section 2.3. Version 2 consists of fast initial guess calculation (as described in Section 2.1) and the NR process described in Section 2.3. At last, Version 3 consists of fast initial guess calculation (as described in Section 2.1), improved initial at sub-pixel level (described in Section 2.2) and the NR process described in Section 2.3. The results of Version 1 algorithm were used as the control or reference for the comparison with Versions 2 and 3.

The image pair shown in Figure 4, with 441 POIs and a grid step of 5 pixels was used for comparing the computational cost of the three aforementioned versions of the DIC algorithm. The subset size used was of $31\times 31$ pixels. In the two-step acceleration IG search scheme, the offset $s$ in Figure 2 was set to 5 pixel which provided a sufficiently high accuracy while the computational cost was controlled. In our experience, the off-set $s$ is to selected in pilot tests, taking the level of displacement, the computational cost, and the accuracy required into consideration.

Figure 5 showed the displacement fields obtained from the DIC computation using the novel NR iteration procedure (Version 3). The transversal displacement increased to a range of $2$ to $5.5$ pixels as the $X$ axis increased. Whereas, the axial displacement decreased to a range of $-2.5$ to $-6.0$ pixels as the $Y$ axis increased. The results were expected and highly consistent with the a prior defined conditions.

Figure 6 showed the computation time related to running the Versions 1 to 3 on a Lenovo laptop (i5 processor, 4G RAM). Two group bars were shown the light blue bar denoted the running time used for IG calculation either from the standard full-field, or from reduced-field, or the reduced field and weighted/fuzzy IG procedure (the NR part are all same); the dark bar denoted the time consumed in the NR process (initial guess was obtained from three different ways as mentioned in previous sentence). The average time consumed for the initial guess calculation by Version 1 was $10.4$ seconds while the time consumed by Versions 2 (accelerated) and 3 (accelerated + weighted) were only $1.2$ seconds which indicated a time reduction of 88.5% over Version 1 (full-field search scheme).

The average time consumed in the NR iteration process in Versions 1 and 2 were $33.8$ seconds and $33.7$ seconds, respectively. This indicated the acceleration process mostly did not affect the computation time of NR iteration. Nevertheless, the average time consumed in the NR process in Version 3 was only 32.4 seconds, and the overall running of the novel procedure (Version 3) was reduced by $10.6$ seconds ($31.5$%) compared to the standard procedure (Version 1), which is substantial. Notably, the reduced and weighted IG search scheme (fuzzy-based initial guess calculation) reduced the NR running time from $33.7$ seconds to $32.4$ seconds, which showed the reduced running time in the new procedure stemmed from the use of the reduced IG search scheme and that fuzzy-based weighting scheme provided a more accurate initial guess. In a related note, the reason why the time consumed in this study was bigger than that reported in other studies ^[6,28,32] was the use of programming language-Matlab.

With respect to the accuracy issue of the novel procedure, Figure 7 showed the absolute average displacement error, standard deviation (SD), maximum and minimum errors in of the displacement measurements from: (1) the integer full-field IG search scheme (with no NR); (2) the novel fractional accelerated/fuzzy scheme in accordance with Section 2.2 (with no NR) and, (3) the final F-NR algorithm which combined (2) and the NR process. It is clear that novel scheme (2) largely enhanced the accuracy of the mean IG displacement value from 0.18 (full-field search scheme) to 0.11 pixel (novel scheme). Specifically, the average error in the displacement measurement from the standard scheme (1) was 0.18 pixel with an SD of 0.24 pixel. The maximum and minimum errors were 0.39 pixel and 0.07 pixel, respectively. It was noteworthy that the absolute error was as big as 18% of the smallest unit/increment, which was one pixel. In the novel accelerated fractional (fuzzy-weighted) scheme, (2) the average error was considerably reduced to 0.11 pixel which is equal to an error reduction by 39% in the IG values. This reduction has significant implications for the likelihood and the speed the convergence. The fuzzier the images or the speckle used, the more effective the novel procedure. Specifically, he maximum and minimum error using the fuzzy approximation process were 0.38 pixel and 0.04 pixel which indicated an error reduction by 2.6% and 43% compared to standard scheme (1) as well. In the NR process, the four types of error estimates were obviously reduced to 0.03, 0.002, 0.04 and 0.03 pixel, respectively. This showed the accurate and precise final measurement results of F-NR algorithm. The error can be even lower when a larger subset size is used and the speckle of the images are not clear. In our laboratory experience, the fuzzy-logics approximation provides a useful tool in dealing with low-resolution image features ^[23,24].

According to the structure of three versions of algorithm, the NR iteration part is the same for all three versions in aspect of algorithm structure. Version 1 used integer IG search scheme and NR iteration process. Version 2 used reduced field search schemed and NR iteration process. Version 3 used reduced field search scheme, fuzzy IG search scheme and NR iteration process. In another word, results in Version 1 is used as control group. Version 2 has an extra step of reducing the IG searching zone comparing to Version 1. Version 3 has an extra step of fuzzy IG searching scheme which gives out more accurate IG comparing to Version 3. And these two steps are exactly the new process that this study has proposed comparing to the conventional algorithm shown in Version 1 in terms of computational cost reduction. Therefore, Figure 6 has shown the results of targets that are of the interest in this study.

In addition, the full-field searching in DIC measurement is mainly calculating the deformation information pixel-wise. The reduce searching scheme proposed in this study reduces the searching zone of each pixel of interest. Although it is not shown in this study, the improvement in computational cost by this scheme is predictable in larger image size in case that a larger number of pixels of interest are to be calculated.

---

3.2. Investigating the measuring limits of F-NR algorithm

To further investigate the limits of the F-NR algorithm (Version 3), a series of image pairs simulating different classes of motions were generated in Table 1. Specifically, ten pairs of images simulating the rigid body translation, rotation and uniaxial tensile conditions were produced and used in this follow-up investigation. For the rigid body translation (RBT) images, translations of 0.1-1 pixel with an increment of 0.1 pixel were studied. Please be aware of the 1-pixel periodicity of errors ^[9,26]. For the images of simulating the rigid body rotation (RBR), images with clockwise rotation were studied due to a symmetric error curve found for images with counter-clockwise rotations in previous studies. For the image pairs of the uniaxial tensile (UAT) condition, both of the axial and transversal strains $\frac{\partial v}{\partial y}$ and $\frac{\partial u}{\partial x}$ were assigned to the deformed image yet the values of only $\frac{\partial v}{\partial y}$ was shown in Table 1. The Poisson ratio was set to 0.33 based on the testing of acrylic medical cements studied in our laboratory.

These image pairs were treated using the F-NR algorithm with the configuration, as described in Section 3.1. The axial strain applied ranged between $1\times10^5$ and $3\times10^5\mu \varepsilon$ with an increment of $0.5\times 10^5 \mu \varepsilon$. The corresponding transversal strains were also properly assigned to image pairs. This added a series of images with the tensile condition which simulated large strains of up to 30% for the sake of investigating the error of the F-NR algorithm for large deformations of interest for medical bone cements.

By applying the F-NR algorithm, the displacement components and gradients were calculated. For the RBR image pairs, the rotation angle was determined from the calculated displacement gradients according to Eq ^[1].

$\begin{equation} \theta = \frac{1}{2}\left( \frac{\partial v}{\partial x}-\frac{\partial u}{ \partial y}\right) , \label{thetaTrans} \end{equation}$

(16)

The displacement, rotation angle and strain results calculated using F-NR algorithm were further compared with the actual/real input data, which were a priori known. Figures 8a, b and c showed the displacement vector for one RBT image pairs ($u = v = 1$ pixel), for one RBR image pair $(\theta = 5^{\circ})$ and for UAT image pair ($\frac{\partial v}{\partial y} = 50000 \mu \varepsilon$ and $\frac{\partial u}{\partial x} = 16500 \mu \varepsilon$), respectively. The displacement vectors in Figure 8 were found to be highly consistent with the real vector orientation and absolute value, respectively.

Furthermore, Figures 9a, b and c showed the comparison of calculated rigid body displacement, rigid rotation angle and strain values to their respective exact values. The green diamond marks denoted the F-NR results at each level applied load and the blue line denoted the exact values. The red error bar was also plotted to show the standard deviation of the F-NR results. Figure 9 clearly demonstrated a high agreement between algorithmic results and exact values. The negative strains denoted the transversal strains while positive strains denoted axial strains in Figure 9c. It is clear that the strain calculated by the F-NR algorithm, even for large deformations of up to $50000 \mu \varepsilon$ were all accurate, fell all on the $45^{\circ}$ straight line with a rather narrow error bar in the studied strain range of up to $50000 \mu \varepsilon$. This agreement demonstrates that the F-NR algorithm was capable of measuring different types of motions and deformations in a rather wide range, as specified in Table 1.

To have a more precise view of algorithmic errors with increased motion, Figures 10a, b and c demonstrated the average absolute error for the three conditions listed in Table 1. Figure 10a shows that the absolute error of the transversal and axial displacement components approximately followed a sinusoidal curve as a systematic error due to the intensity interpolation ^[9,26]. And errors were close to zero at 0, 0.5, and 1 pixel positions which was reported by earlier study of Schreier, et al. (2000) ^[26]. The standard deviations at the 10 displacement increments were almost identical in value. The algorithmic overall absolute error fell between $-0.001$ and $0.001$ pixel which approximately equaled to 1% of relative error at its maximum.

Figure 10b shows a negatively increasing average error in algorithmically calculated rotation angle versus applied rotation angle. Specifically, the error ranged from 0 to approximately $0.006^{\circ}$ which corresponded to a maximum relative error of 1.2% at $5^{\circ}$ RBR. The standard deviations bar also remained in a small range which showed the stability of the F-NR algorithm.

Figure 10c show that the absolute average error of the calculated axial strain fluctuated between -60 $\mu \varepsilon$ and 60 $\mu \varepsilon$. The negative values in both axes denoted the transversal strain while the positive strain denoted axial strain (tensile). Please note that the increment of simulated strain in the narrow range of -16500~50000 $\mu \varepsilon$ was small. However, for a better appreciation of the results, an enlarged detail was shown in the subfigure of Figure 10c. It was seen that the error of the calculated strains still fell in the range of -60 to 60 $\mu \varepsilon$. This implied a relative algorithmic error of only 0.5% when processing images simulating strains of up to $3\times 10^5\mu \varepsilon$.

---

4. Conclusion

We were able to confirm the hypothesis of the study in that using the novel search scheme, the IG values became more accurate and the overall computational costs of running the DIC algorithm were reduced by 31.5%, which is substantial. Specifically, compared to the full-field search scheme, the accelerated scheme reduced the computational cost for the IG values by 88.5% (Figure 6). Furthermore, Figure 7 clearly showed that the average IG value was considerably closer to the average true displacement value if compared to the full-field search scheme. The F-NR algorithm was verified by a wide range of image pairs simulating rigid body translations and rotations as well as the uniaxial tensile strains of up to 30%. The overall error of three motion types remained smaller than 1.2% which indicated an algorithmic accuracy of 98.8%. The relative standard deviation was smaller than 1% which implied a precision rate of 99%. Our study showed that the novel F-NR algorithm is accurate, precise, and efficient. This algorithm was also robust for all types of displacement and deformation tested.

It is believed that the new scheme using reduced searching zone and fuzzy-logic based scheme for initial guess would perform better in images with larger image size and deformation due to the pixel-wise initial guess calculation. The results obtained in this study concluded that the novel algorithm is an efficient technique for accurately measuring full-field displacements and a wide range of deformations. The two extra steps and this complete algorithm is potentially useful in DIC measurement during the mechanical testing of bone cement specimens with sophisticated surface conditions during deformation.

---

Acknowledgments

The funding from Chinese Scholarship Council, Natural Sciences and Engineering Research Council of Canada (NSERC) are acknowledged. Dr Liang Wang and Dr. Ahmed Sweedy are acknowledged for the valuable discussions during the algorithm development and manuscript revision in this study, respectively.

---

Conflict of interest

The authors declare no conflict of interest.

---