High Accurate Crack Reconstruction Based on an Improved Discontinuous Digital Image Correlation: Subset Restore and Adaptation Method

Digital image correlation (DIC) is an efficient nondestructive technique for measuring surface displacement in engineering. However, standard DIC is restricted to continuous deformation, and the existing discontinuous DIC (DDIC) techniques are only able to measure small-scale cracks. In this report, a novel subset restore model and a corresponding subset size adaptation algorithm are presented to overcome this limitation for crack-state and displacement field reconstruction for large-scale cracks. +e technique introduces a new subset restore method for splicing the segmented subset by tracing the motion trajectory caused by pure discontinuities. +e proposed model facilitates the calculation of the rotation angle and the pivot of the subset movement. +e subset size adaptation algorithm is designed based on an evaluation of the intensity gradient and correlation coefficient to allow themodel to achieve high accuracy. Validation of the approach was performed using two typical crackmodels, by deforming a numerically synthesized Gaussian speckle image according to the deformation data from finite element analysis (FEA) results and photographing a laboratory tensile test with a high-speed CCD camera, respectively. +e results validate the efficacy and high accuracy of the proposed approach compared to standard DIC in the reconstruction of the displacement fields in both continuous and discontinuous regions. +e accuracy of resultant displacement reconstruction achieves approximately 0.015 pixels and 0.05 pixels in continuous region and crack vicinity, respectively.


Introduction
Cracks are among the most common defects in engineered structures. Numerous nondestructive techniques, such as eddy current testing (ECT) [1], ultrasonic flaw detection [2], infrared thermal wave method [3], electromagnetic detection [4], and crack image recognition [5][6][7], have been widely used to detect and identify cracks in structures. e accurate measurement and reconstruction of surface displacement, deformation, and crack status are the key to the acquisition of qualitative and quantitative information to better understand the mechanism of load-resistance and the failure modes of cracked structure. A feasible noncontact surface displacement measurement method is the digital image correction (DIC) method, which was proposed by Peters and Ranson [8] and Yamaguchi [9] and tested by Chu et al. [10][11][12][13]. Nevertheless, the standard DIC method is based on the hypothesis that the pending surface is continuous everywhere, which will lead to serious errors in the vicinity of cracks.
To overcome this limitation, Réthoré et al. [14,15] utilized the extended finite element method (XFEM) [16] and proposed an extended DIC method by "enriching" the shape function of the elements to consider the presence of a crack.
e image is meshed with elements similar to the finite element method, and the shape function is enriched to yield nodal displacements. is technique is fundamentally different from the standard DIC. Moreover, it is computationally intensive. A novel pointwise DIC method was developed by Jin and Bruck [17,18]. ey considered the vertical and horizontal displacements of each pixel on the image as parameters and simultaneously optimized them within a subset. A large number of unknowns need to be optimized at one time using the genetic algorithms, which requires significant computational resources and may lead to poor universality. Helm [19] proposed a modified approach to identify a crack that exploits the consistency of quasiregular speckle patterns. However, due to the complexity of the generated surface, this method is still limited to the laboratory examples. Pan [20] suggested the use of a reliability guided (RG) DIC method to ensure that the calculation path was always along the points with the highest zero-mean normalized cross-correlation (ZNCC) coefficient and to avoid calculations for discontinuous areas. Based on this idea, Blaber et al. [21] developed an automatic program to combine an inverse compositional method with the RG DIC, which is theoretically capable of approaching discontinuous regions. Poissant and Barthelat [22] proposed a novel subset splitting technique and assumed that a crack is a straight line that divides the subregion into two parts. is is invalid in most engineering situations. Recently, Hassan et al. [23,24], Han and Pan [25], and Tang and Xiao [26] studied the relationship between the two parts of the divided subset when discontinuities pass through and established different discontinuous subset splitting models to improve the accuracy when the displacement and deformation of a discontinuous region are calculated. e results show that these approaches work well when discontinuous deformations are not strenuous. However, when large translations and rotations occur along the discontinuous path, they may lead to erroneous results.
Herein, we propose a novel improved DDIC method to handle crack discontinuities, with the aim of reducing errors and enhancing applicability. e technique is mainly based on a subset restore model to trace the motion trajectory of a subset caused by pure discontinuities, and a corresponding subset size adaptation algorithm to improve the correlation accuracy. e details are explained in Section 2. A special simulation experimental method for testing the proposed approach is introduced in Section 3 by deforming a reference image according to the deformation data from FEA. Two typical examples, a plate with a long inclined crack and a cracking process in a central hole plate, are investigated by the proposed simulation method and the laboratory tensile test, respectively, to evaluate the proposed technique in the latter part of Section 3. Finally, a brief conclusion and potential applications of the proposed technique are presented in Section 4.

Principles
e standard DIC method is a statistical process that correlates pixel values of two corresponding regions called subsets in a deformed and a reference image. Digitized grayscale images of the speckled specimen surface are acquired (using a CCD camera, for example) as pending images [27,28]. e kernel shape function is based on a continuous assumption that fails when discontinuities occur. An alternative model is needed to overcome this limitation.
Consider the displacement of pixels on the surface before and after deformation and cracking. A pixel region with size (2m + 1) × (2m + 1) is selected as a subset, and point P is the geometric center of the subset called the focus point. Due to the occurrence of a crack, the subset is divided into two regions by the crack face, defining a region that includes P as the master subset and the remaining region as the slave subset, as shown in Figure 1.
ere is a point A (x a , y a ) in the master subset far from the crack line that ensures a continuous subset containing A.
e position of this point in the deformed image can be obtained by standard DIC as A * (x * a , y * a ). Similarly, another point B (x b , y b ) can be found, and its position B * (x * b , y * b ) in the deformed image can be determined. It is assumed that the change in distance of the two points A and B is negligible after deformation within one subset compared to the discontinuous deformation, which can be expressed as |AB|≈| A * B * |. us, there must be a point O m (x m , y m ) and an angle θ m that allows AB and A * B * to be in mutual coincidence after the rotation of A * B * through an angle θ m with the pivot O m . e rotation equations of points A and B to points A * and B * through an angle θ m with pivot O m can be written as follows: e parameters x m , y m , and θ m can be solved using equation (1), which represent the rotation of the master subset by pure crack deformation. Similarly, the rotation of the slave subset by pure crack deformation can be obtained for the pivot point O s (x s , y s ) and a rotation angle θ s . us, the deformation of a discontinuous subset can be divided into two parts: a pure continuous deformation and a pure discontinuous deformation.
Assume the focus point P of a selected subset in the reference image is at (x p , y p ), and point Q (x q , y q ) is an arbitrary point in the subset. en, the movement to positions Q * (x * q , y * q ) and P * (x * p , y * p ) in the deformed image is as shown in Figure 2.
Consider that a crack exists in all the subsets, and for the continuous subset, the rotation angle θ m � θ s � 0.
us, a unified subset restore model can be used to deal with both the continuous and discontinuous subsets. e geometrical relationship can be established as follows: where u p v p and u q v q represent the pure continuous displacements of points P and Q, respectively. When the selected subset is continuous, let θ m � θ s � 0. Equation (2) can then be simplified as a pure continuous deformation: It can determine that the displacements u q v q of Q can be derived from the displacements u p v p of P when the subset is continuous as follows: Pending subset

A B
Master subset Slave subset

Mathematical Problems in Engineering 3
After substituting equation (4) into equation (2), the position of Q * (x * q , y * q ) can be expressed as follows: where pivot O s (x s , y s ) and the rotation angle θ s can be calculated from equation (1). us, the position Q * (x * q , y * q ) of any arbitrary point Q after deformation can be obtained from the pure continuous deformation vector p � u v zu/zx zu/zy zv/zx zv/zy . Let us define the gray intensity distribution of the reference image as f(x, y) and similarly f * (x, y) for the deformed image. en, the gray intensity of Q * can be written as follows: (6) proves that if the displacement p � u v zu/zx zu/zyzv/zxzv/zy] of the focus point P is known, then the position of any nearby point Q * can be determined. Conversely, if any point Q in the reference image can be found and the position can be computed in the deformed image as Q * , the displacement p of the focus point P can be estimated.

Reconstruction of the Displacement Field. Equation
To estimate the positions of P and any Q, the crosscorrelation coefficient (CC) and the sum of the squared difference coefficient (SSD) can be used to determine the degree of matching. eir zero normalized functions are preferred [29,30] and can be presented as follows: where S stands for the set of all pixel positions in the selected subset and the others are given as follows: It is easy to determine that C ZNCC ≤ 1 and C ZNSSD ≥ 0 from equation (7), and the closer the C ZNCC and C ZNSSD are to their limits, the higher the degree of matching of the pending subsets. e results from Pan [31] prove that C ZNCC and C ZNSSD are mathematically equivalent and can be quantitatively expressed as For simplicity, the C ZNSSD is selected as the correlation coefficient to determine the degree of matching in this report, and all correlation coefficients subsequently mentioned refer to C ZNSSD unless otherwise stated.
An appropriate interpolation technique must be utilized to achieve high accuracy. Several image interpolation functions could be used, e.g., B-spline functions, o-Moms functions, and Key's function [32], to estimate the gray intensity of any region between pixels. e bicubic B-spline function is selected because it is globally symmetric and second-time continuously differentiable. e equation for the B-spline kernel is given as follows: where by definition, we have By minimizing the correlation coefficient C ZNSSD , the overall displacements of the image can be estimated. For this purpose, let C ZNSSD assume its minimum value by computing the following: 4 Mathematical Problems in Engineering where i is the iteration. e classic NR minimization technique [23] is used to estimate and improve p in equation (12) as follows: where [H i ] is the Hessian matrix that contains the second derivatives of p i at each iteration. Finally, the reconstruction of the crack-state and full-field displacements can be acquired by combining the pure continuous and discontinuous deformation.

Subset Adaptation
Scheme. e selection of the subset and calculation order is important for the accuracy of the results. e results of Pan et al. [33] show that the standard deviation (SD) error of the displacement result decreases as the sum of the square of the subset intensity gradient (SSSIG) increases during the processing of the continuous DIC.
us, a subset adaptation scheme is performed to ensure that each subset can be calculated for a suitable position and size to minimize the error rate.
In this case, an intensity gradient function G(x, y) is introduced together with the correlation coefficient C to determine the position of the focus point P and the size of the selected subset as follows: where G x and G y are the sum of the square of the subset intensity gradients in the x and y directions, respectively, and can be calculated as follows: To implement the subset adaptation scheme, an initial focus point P is artificially selected as a "seed." Normally, a pixel point far away from the crack area with a tiny movement after deformation is suggested. us, the initial subset size M may be expressed as where l and w represent the length and width of the region of interest (ROI) in pixels. Before the process is initiated, the ROI needs to be determined and interpolated as previously indicated. e first part of the subset adaptation scheme is "seeding" of the reference image. According to the accuracy requirements of the calculation, the expectation of G(x, y) needs to be set as G th . During this process, as shown in Figure 3, the initial seed and subset size M are used to calculate its intensity gradient function G(x, y). If G(x, y) < G th , then we let M � M + 2 and use it to calculate the new G(x, y) for comparison with G th . e process is repeated, and M is accumulated until G(x, y) ≥ G th or G(x, y) does not increase further. e four vertexes of the subset are then selected as new pending seeds, the locations of these seeds are checked, and the pending seeds located in existing subsets are canceled. is process is repeated until the entire ROI is seeded, and then the subset and seed data are recorded. Subsequently, during the second part, the DDIC method is applied to each seed and the initial rotation angle is set as θ � 0 to reduce the required preliminary calculations. e correlation coefficient threshold for identification of the discontinuous subsets is set as C cr . e subsets that satisfy condition C ≥ C cr are identified as discontinuous subsets, and θ is no longer zero. e pivot O(x, y) plus θ can then be calculated using equation (1), and a new C can be obtained for the discontinuous subsets. A check is then performed to determine if all the Cs are in an acceptable range; otherwise, the "error" subsets are adjusted similar to the first part of the process.
Finally, the subpixel NR algorithm is applied to obtain the full-field displacements at the subpixel level; thus, the pure continuous deformation and pure discontinuous deformation of every subset can be obtained. en, for every subset, use the pure discontinuous deformation to calculate the crack line (ignore those subsets where pure discontinuous deformation is equal to zero), and summarize them to reconstruct the whole crack. A summary of the subset adaptation scheme is presented as a flow chart in Figure 4.

Validation and Verification
3.1. Acquisition of Test Images. Digital images for validation of the DIC method were mainly acquired from laboratory experiment [10][11][12][13]34] and digital simulation [35]. Although a laboratory DIC experiment is a complex procedure, a typical digital simulation to obtain reference and deformed images is much easier based on the assumption of deformation. Each approach has its advantages and disadvantages, as shown in the comparison in Table 1.
An improved solution is presented to simplify the acquisition of digital images and to overcome its limitations.
is entails the generation of sufficiently complex and random reference images via the accumulation of individual Gaussian speckle, to simulate real deformation via finite element analysis (FEA) and the deformation of reference images according to simulated deformations. e gray intensity distribution of the normal Gaussian speckle image can be computed as follows: where I 0 is the maximum light intensity of each speckle, S is the total number of speckles in the image, R is the size of each speckle, and [x k , y k ] represents the center location of the kth speckle.

Mathematical Problems in Engineering
More than 2.25 × 10 10 calculations of e ω (e's power calculation) are required to obtain a high-resolution speckled image with a resolution of 1000 × 1500 pixels and 15000 speckles. is significantly decreases the computation speed. It should be noted that, for the expression e −16 ≈1.125 × 10 −7 , the effect of the Gaussian speckle beyond a radius of 3R can be ignored. us, an additional condition can be added to equation (17) to reduce the computation as follows: ����������������� After the generation of the reference image, an FEA model that is consistent with a real situation in engineering is established and analyzed using an FEA software (Abaqus is suggested for discontinuous analysis). e deformation of the model can be obtained as U � u u x u y and V � v v x v y , and the gray intensity distribution of the deformed image can be derived as follows: It should be noted that the coordinates obtained using equation (19) are no longer integers. is indicates that the interpolation method should be applied to generate a valid deformed image after deformation. Comparison of laboratory experiments and normal digital simulation is listed in Table 1.

Simulation: Deformation of Tension Plate with a Long
Inclined Crack. Tension load is one of the most common engineering situations, and it is used in this case to test along an inclined precrack plate. A specimen model is established using carbon steel Q235 (E � 210 GPa and ] � 0.28) with a size of 100 mm × 200 mm × 10 mm, as shown in Figure 5(a). e inclined precrack starts at the middle of the vertical edge at an angle of 20°to the horizontal with a length of 30 mm. e displacement results are obtained using Abaqus and presented as shown in Figure 5(b). e reference image with a resolution of 1000 × 2000 pixels is obtained by applying equations (17) and (18)  Step 1 Step 2 Step n substituting the FEA displacements results into equation (19). ese two synthetic images are shown in Figures 5(c) and 5(d), respectively. It should be noted that an ROI with a size of 600 × 800 pixels is set as shown in Figures 5(c) and 5(d), and the unit conversion factor used to convert a length unit into a pixel unit is 10 pixel/mm. e synthetic reference image and deformed image are analyzed using the proposed discontinuous DIC method to reconstruct the displacements and crack state. In addition, the standard DIC method is used to obtain similar results for subsequent comparison. e reconstruction of the displacement field in both the horizontal and vertical directions is obtained using the proposed discontinuous DIC and standard DIC, respectively, and are presented as shown in Figure 6.
It is evident from Figure 6 that the proposed discontinuous DIC method works well in both the continuous area and the discontinuous area, whereas the standard DIC is invalid in the vicinity of the crack. Based on a comparison with Figure 5, it can be determined that the reconstructed displacement fields and crack state coincide with the FEA results. e correlation coefficients of the proposed discontinuous DIC and standard DIC were investigated for the ROI, and the results were collected and reorganized. e distribution of the correlation coefficient probability is shown in Figure 7(a), while the correlation coefficient values along the crack are shown in Figure 7(b).
It is evident that the distribution of the correlation coefficients obtained using the proposed method is more concentrated, and its values are smaller compared to results from standard DIC, which can result in a smoother displacement field, as shown in Figure 6. Furthermore, a comparison of the correlation coefficients along the crack indicates that the standard DIC is unsuccessful for the identification of displacements in the vicinity of the crack, while the proposed method facilitates reconstruction of the crack and the displacements.    Figures 8(a) and 8(b), respectively, along with the results obtained for the standard DIC and theoretical crack displacements acquired using FEA. e results indicate that the proposed method achieves the design goal.
A quantitative evaluation of the calculated displacements is adopted using the mean bias error and the standard deviation error as criteria. e mean bias error of the displacement result is defined as follows: where u and u real denote the measured and real displacements of each point in the ROI, respectively.  Mathematical Problems in Engineering e standard deviation error of the measured displacement is expressed as follows: e results of different mean bias errors and standard deviation errors calculated using the proposed method and standard DIC are listed in Table 2. e 20-pixel vicinity of the crack is defined as the area around the discontinuities to demonstrate the capability of appropriately addressing this region. e statistics in Table 2 show that the standard DIC fails to address the discontinuous problem, while the proposed method is capable of reconstructing the crack state with subpixel accuracy.

Laboratory Experiment: Cracking Process of the Central
Hole Tension Plate. Another experiment was performed in a central hole plate to evaluate the ability of the proposed technique to recognize and reconstruct the crack propagation process when a discontinuous area (the central hole) already exists. A thin plate model was made using high-density polyethylene (HDPE) with a through-hole in the geometric center of the plate. Two small notches are located at the edge of the through-hole, as shown in Figure 9  the major area which is "painted" are 80 mm × 140 mm × 2 mm, and the radius of the hole is 10 mm. A tension load of 30 μm/s is applied in the vertical direction of the specimen. A plot of the displacement and its corresponding load force on the edge of specimen is shown in Figure 9(b). e total test time is 100 s, and pictures are taken every 0.1 s by the high-speed CCD camera with a resolution of 960 × 1280 pixels. e first image is selected as the reference image. Ten moments of the tensile test process are selected as typical deformations of different cracking stages with crack lengths l c equal to 1 mm, 2 mm, and 3 mm-10 mm, as well as ten images are found according to their crack lengths as mentioned and are set as deformed images. Due to the symmetry of the model, an ROI is set with a size of 350 × 500 pixels, as shown in Figure 10; the ten deformed ROIs are presented in 10-time amplifications to illustrate the cracking process. It should be noted that the unit conversion factor for converting a length unit into a pixel unit is 10 pixel/mm. ese ten deformed images along with the reference image are substituted in the proposed method and the standard DIC, respectively. e preliminary run with the proposed method yields displacements fields, and the results from crack stage l c � 10 mm are selected and displayed, as shown in Figure 11. Similar results for the standard DIC are also presented for comparison. e reconstruction of the displacement fields demonstrates the capability of the proposed method to handle complex discontinuous deformations, while the standard DIC is limited in this respect. e reconstruction of the crack in each crack stage is calculated in both the horizontal and vertical directions. e results for three representative crack stages are selected for crack length l c of 2 mm, 6 mm, and 10 mm to illustrate the validity and accuracy of crack displacement reconstruction, as shown in Figures 12(a)-12(f ).
It is evident from Figure 12 that the proposed method achieves high accuracy for all the crack processes. e  fluctuation of the absolute error for each test point increases very little with crack propagation because the subset restore model is applied. e detailed correlation coefficients of the ten deformed images while processing the proposed method are collected to verify the quality of the correlation. In addition, the resultant standard deviation error of the displacement is adopted to quantitatively evaluate the reconstruction results and can be calculated as follows: Data are counted separately in the entire ROI, crack vicinity, and the continuous area to investigate the effect of crack propagation on the accuracy of the reconstruction, and the results are rearranged, as shown in Figures 13(a) and 13(b). It is evident that crack propagation affects the accuracy of the results, but this effect is limited. e correlation

Conclusion
is article presents an improved DDIC method based on the subset restore model and subset adaptation algorithm to address the discontinuous problems in standard DIC. e kernel of the proposed method is to trace the motion trajectory of the subset caused by pure discontinuities, to calculate the rotation angle plus the pivot, and to restore the separate master subset and slave subset. e deformation of the tension plate with a long inclined crack and the cracking process of the central hole tension plate were investigated by an improved simulation method and a laboratory tensile test, respectively. e results indicated that the proposed technique is superior for the reconstruction of the displacement fields for both continuous and discontinuous areas when compared to the standard DIC method. Furthermore, the proposed method is capable of handling not only small discontinuities but also cracks with large translations and rotations. e high reconstruction accuracy of the procedure facilitates the extraction of displacement in the vicinity of the crack and the crack tip, which is beneficial in fracture mechanics and can be applied to quantitative crack growth monitoring.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.