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

Authors: Wenzhi Tang; Hanbin Xiao; Chentong Chen
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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. The technique introduces a new subset restore method for splicing the segmented subset by tracing the motion trajectory caused by pure discontinuities. The proposed model facilitates the calculation of the rotation angle and the pivot of the subset movement. The subset size adaptation algorithm is designed based on an evaluation of the intensity gradient and correlation coefficient to allow the model to achieve high accuracy. Validation of the approach was performed using two typical crack models, 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. The 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. The accuracy of resultant displacement reconstruction achieves approximately 0.015 pixels and 0.05 pixels in continuous region and crack vicinity, respectively.


---

## Body

## 1. 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–7], have been widely used to detect and identify cracks in structures. The 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–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. The image is meshed with elements similar to the finite element method, and the shape function is enriched to yield nodal displacements. This 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]. They 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. This 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. The 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. The 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. The details are explained in Section2. 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.


## 2. Principles

### 2.1. Subset Restore Model

The 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]. The 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.Figure 1
Deformation of a cracked subset.There is a pointA (xa, ya) in the master subset far from the crack line that ensures a continuous subset containing A. The position of this point in the deformed image can be obtained by standard DIC as A∗xa∗,ya∗. Similarly, another point B (xb, yb) can be found, and its position B∗xb∗,yb∗ 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∗|. Thus, there must be a point Om (xm, ym) 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 Om. The rotation equations of points A and B to points A∗ and B∗ through an angle θm with pivot Om can be written as follows:(1)xa∗−xmya∗−ym=xa−xmya−ymcosθmsinθm−sinθmcosθm,xb∗−xmyb∗−ym=xb−xmyb−ymcosθmsinθm−sinθmcosθm.The parametersxm, ym, 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 Os (xs, ys) and a rotation angle θs. Thus, the deformation of a discontinuous subset can be divided into two parts: a pure continuous deformation and a pure discontinuous deformation.Assume the focus pointP of a selected subset in the reference image is at (xp, yp), and point Q (xq, yq) is an arbitrary point in the subset. Then, the movement to positions Q∗xq∗,yq∗ and P∗xp∗,yp∗ in the deformed image is as shown in Figure 2.Figure 2
The relationship of pixels in a cracked subset.Consider that a crack exists in all the subsets, and for the continuous subset, the rotation angleθm = θs = 0. Thus, a unified subset restore model can be used to deal with both the continuous and discontinuous subsets. The geometrical relationship can be established as follows:(2)xp∗−xmyp∗−ym=xp−xmyp−ymcosθmsinθm−sinθmcosθm+upvp,xq∗−xsyq∗−ys=xq−xsyq−yscosθssinθs−sinθscosθs+uqvq,where upvp and uqvq 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:(3)xp∗yp∗=xpyp+upvp,xq∗yq∗=xqyq+uqvq.It can determine that the displacementsuqvq of Q can be derived from the displacements upvp of P when the subset is continuous as follows:(4)uq=up+∂up∂xdx+∂up∂ydy,vq=vp+∂vp∂xdx+∂vp∂ydy.After substituting equation (4) into equation (2), the position of Q∗xq∗,yq∗ can be expressed as follows:(5)xq∗=xq−xscosθs−yq−yssinθs+up+∂up∂xdx+∂up∂ydy,yq∗=xq−xssinθs+yq−yscosθs+vp+∂vp∂xdx+∂vp∂ydy,where pivot Os (xs, ys) and the rotation angle θs can be calculated from equation (1). Thus, the position Q∗xq∗,yq∗ of any arbitrary point Q after deformation can be obtained from the pure continuous deformation vector p=uv∂u/∂x∂u/∂y∂v/∂x∂v/∂y. Let us define the gray intensity distribution of the reference image as f(x, y) and similarly f∗(x, y) for the deformed image. Then, the gray intensity of Q∗ can be written as follows:(6)f∗Q∗=f∗xq−xscosθs−yq−yssinθs+u+∂u∂xdx+∂u∂ydyxq−xssinθs+yq−yscosθs+v+∂v∂xdx+∂v∂ydy.


### 2.2. Reconstruction of the Displacement Field

Equation (6) proves that if the displacement p=uv∂u/∂x∂u/∂y∂v/∂x∂v/∂y 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.To estimate the positions ofP and any Q, the cross-correlation coefficient (CC) and the sum of the squared difference coefficient (SSD) can be used to determine the degree of matching. Their zero normalized functions are preferred [29, 30] and can be presented as follows:(7)CZNCC=∑Sfx,y−f¯f∗x∗,y∗−f∗¯ΔfΔf∗,CZNSSD=∑Sfx,y−f¯Δf−f∗x∗,y∗−f∗¯Δf∗2,where S stands for the set of all pixel positions in the selected subset and the others are given as follows:(8)f¯=12m+12∑Sfx,y;f∗¯=12m+12∑Sf∗x∗,y∗,Δf=∑Sfx,y−f¯2;Δf∗=∑Sf∗xi∗,yi∗−f∗¯2.It is easy to determine thatCZNCC ≤ 1 and CZNSSD ≥ 0 from equation (7), and the closer the CZNCC and CZNSSD are to their limits, the higher the degree of matching of the pending subsets. The results from Pan [31] prove that CZNCC and CZNSSD are mathematically equivalent and can be quantitatively expressed as(9)CZNSSD=2−2CZNCC.For simplicity, theCZNSSD is selected as the correlation coefficient to determine the degree of matching in this report, and all correlation coefficients subsequently mentioned refer to CZNSSD 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. The bicubic B-spline function is selected because it is globally symmetric and second-time continuously differentiable. The equation for the B-spline kernel is given as follows:(10)βnx=∑k=0n+1−1kn+1n+1−k!k!n+12+x−k+n,∀x∈R,∀n∈N∗where by definition, we have(11)x+n=max0,xn,n>0.By minimizing the correlation coefficientCZNSSD, the overall displacements of the image can be estimated. For this purpose, let CZNSSD assume its minimum value by computing the following:(12)∇Cp=∂CZNSSD∂pii=1,n=0,where i is the iteration. The classic NR minimization technique [23] is used to estimate and improve p in equation (12) as follows:(13)pi+1=pi−Hi−1∇Cpi,where Hi is the Hessian matrix that contains the second derivatives of pi 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.


### 2.3. Subset Adaptation Scheme

The selection of the subset and calculation order is important for the accuracy of the results. The 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. Thus, 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 functionGx,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:(14)Gx,y=GxGy,where Gx and Gy are the sum of the square of the subset intensity gradients in the x and y directions, respectively, and can be calculated as follows:(15)Gx=∑Sfx,y∂x2,Gy=∑Sfx,y∂y2.To implement the subset adaptation scheme, an initial focus pointP is artificially selected as a “seed.” Normally, a pixel point far away from the crack area with a tiny movement after deformation is suggested. Thus, the initial subset size M may be expressed as(16)M=2m+1,m=intminl,w50,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. The first part of the subset adaptation scheme is “seeding” of the reference image. According to the accuracy requirements of the calculation, the expectation of Gx,y needs to be set as Gth. During this process, as shown in Figure 3, the initial seed and subset size M are used to calculate its intensity gradient function Gx,y. If Gx,y<Gth, then we let M = M + 2 and use it to calculate the new Gx,y for comparison with Gth. The process is repeated, and M is accumulated until Gx,y≥Gth or Gx,y does not increase further. The 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. This process is repeated until the entire ROI is seeded, and then the subset and seed data are recorded.Figure 3
The “seeding” procedure for adapting the subset size.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. The correlation coefficient threshold for identification of the discontinuous subsets is set as Ccr. The subsets that satisfy condition C ≥ Ccr are identified as discontinuous subsets, and θ is no longer zero. The 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. Then, 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 Figure4.Figure 4
A flow chart summary of the subset adaptation scheme.


## 2.1. Subset Restore Model

The 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]. The 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.Figure 1
Deformation of a cracked subset.There is a pointA (xa, ya) in the master subset far from the crack line that ensures a continuous subset containing A. The position of this point in the deformed image can be obtained by standard DIC as A∗xa∗,ya∗. Similarly, another point B (xb, yb) can be found, and its position B∗xb∗,yb∗ 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∗|. Thus, there must be a point Om (xm, ym) 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 Om. The rotation equations of points A and B to points A∗ and B∗ through an angle θm with pivot Om can be written as follows:(1)xa∗−xmya∗−ym=xa−xmya−ymcosθmsinθm−sinθmcosθm,xb∗−xmyb∗−ym=xb−xmyb−ymcosθmsinθm−sinθmcosθm.The parametersxm, ym, 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 Os (xs, ys) and a rotation angle θs. Thus, the deformation of a discontinuous subset can be divided into two parts: a pure continuous deformation and a pure discontinuous deformation.Assume the focus pointP of a selected subset in the reference image is at (xp, yp), and point Q (xq, yq) is an arbitrary point in the subset. Then, the movement to positions Q∗xq∗,yq∗ and P∗xp∗,yp∗ in the deformed image is as shown in Figure 2.Figure 2
The relationship of pixels in a cracked subset.Consider that a crack exists in all the subsets, and for the continuous subset, the rotation angleθm = θs = 0. Thus, a unified subset restore model can be used to deal with both the continuous and discontinuous subsets. The geometrical relationship can be established as follows:(2)xp∗−xmyp∗−ym=xp−xmyp−ymcosθmsinθm−sinθmcosθm+upvp,xq∗−xsyq∗−ys=xq−xsyq−yscosθssinθs−sinθscosθs+uqvq,where upvp and uqvq 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:(3)xp∗yp∗=xpyp+upvp,xq∗yq∗=xqyq+uqvq.It can determine that the displacementsuqvq of Q can be derived from the displacements upvp of P when the subset is continuous as follows:(4)uq=up+∂up∂xdx+∂up∂ydy,vq=vp+∂vp∂xdx+∂vp∂ydy.After substituting equation (4) into equation (2), the position of Q∗xq∗,yq∗ can be expressed as follows:(5)xq∗=xq−xscosθs−yq−yssinθs+up+∂up∂xdx+∂up∂ydy,yq∗=xq−xssinθs+yq−yscosθs+vp+∂vp∂xdx+∂vp∂ydy,where pivot Os (xs, ys) and the rotation angle θs can be calculated from equation (1). Thus, the position Q∗xq∗,yq∗ of any arbitrary point Q after deformation can be obtained from the pure continuous deformation vector p=uv∂u/∂x∂u/∂y∂v/∂x∂v/∂y. Let us define the gray intensity distribution of the reference image as f(x, y) and similarly f∗(x, y) for the deformed image. Then, the gray intensity of Q∗ can be written as follows:(6)f∗Q∗=f∗xq−xscosθs−yq−yssinθs+u+∂u∂xdx+∂u∂ydyxq−xssinθs+yq−yscosθs+v+∂v∂xdx+∂v∂ydy.


## 2.2. Reconstruction of the Displacement Field

Equation (6) proves that if the displacement p=uv∂u/∂x∂u/∂y∂v/∂x∂v/∂y 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.To estimate the positions ofP and any Q, the cross-correlation coefficient (CC) and the sum of the squared difference coefficient (SSD) can be used to determine the degree of matching. Their zero normalized functions are preferred [29, 30] and can be presented as follows:(7)CZNCC=∑Sfx,y−f¯f∗x∗,y∗−f∗¯ΔfΔf∗,CZNSSD=∑Sfx,y−f¯Δf−f∗x∗,y∗−f∗¯Δf∗2,where S stands for the set of all pixel positions in the selected subset and the others are given as follows:(8)f¯=12m+12∑Sfx,y;f∗¯=12m+12∑Sf∗x∗,y∗,Δf=∑Sfx,y−f¯2;Δf∗=∑Sf∗xi∗,yi∗−f∗¯2.It is easy to determine thatCZNCC ≤ 1 and CZNSSD ≥ 0 from equation (7), and the closer the CZNCC and CZNSSD are to their limits, the higher the degree of matching of the pending subsets. The results from Pan [31] prove that CZNCC and CZNSSD are mathematically equivalent and can be quantitatively expressed as(9)CZNSSD=2−2CZNCC.For simplicity, theCZNSSD is selected as the correlation coefficient to determine the degree of matching in this report, and all correlation coefficients subsequently mentioned refer to CZNSSD 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. The bicubic B-spline function is selected because it is globally symmetric and second-time continuously differentiable. The equation for the B-spline kernel is given as follows:(10)βnx=∑k=0n+1−1kn+1n+1−k!k!n+12+x−k+n,∀x∈R,∀n∈N∗where by definition, we have(11)x+n=max0,xn,n>0.By minimizing the correlation coefficientCZNSSD, the overall displacements of the image can be estimated. For this purpose, let CZNSSD assume its minimum value by computing the following:(12)∇Cp=∂CZNSSD∂pii=1,n=0,where i is the iteration. The classic NR minimization technique [23] is used to estimate and improve p in equation (12) as follows:(13)pi+1=pi−Hi−1∇Cpi,where Hi is the Hessian matrix that contains the second derivatives of pi 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.


## 2.3. Subset Adaptation Scheme

The selection of the subset and calculation order is important for the accuracy of the results. The 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. Thus, 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 functionGx,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:(14)Gx,y=GxGy,where Gx and Gy are the sum of the square of the subset intensity gradients in the x and y directions, respectively, and can be calculated as follows:(15)Gx=∑Sfx,y∂x2,Gy=∑Sfx,y∂y2.To implement the subset adaptation scheme, an initial focus pointP is artificially selected as a “seed.” Normally, a pixel point far away from the crack area with a tiny movement after deformation is suggested. Thus, the initial subset size M may be expressed as(16)M=2m+1,m=intminl,w50,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. The first part of the subset adaptation scheme is “seeding” of the reference image. According to the accuracy requirements of the calculation, the expectation of Gx,y needs to be set as Gth. During this process, as shown in Figure 3, the initial seed and subset size M are used to calculate its intensity gradient function Gx,y. If Gx,y<Gth, then we let M = M + 2 and use it to calculate the new Gx,y for comparison with Gth. The process is repeated, and M is accumulated until Gx,y≥Gth or Gx,y does not increase further. The 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. This process is repeated until the entire ROI is seeded, and then the subset and seed data are recorded.Figure 3
The “seeding” procedure for adapting the subset size.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. The correlation coefficient threshold for identification of the discontinuous subsets is set as Ccr. The subsets that satisfy condition C ≥ Ccr are identified as discontinuous subsets, and θ is no longer zero. The 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. Then, 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 Figure4.Figure 4
A flow chart summary of the subset adaptation scheme.


## 3. Validation and Verification

### 3.1. Acquisition of Test Images

Digital images for validation of the DIC method were mainly acquired from laboratory experiment [10–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.Table 1
Comparison of different image acquisition methods.
Description
Acquisition methods for digital images
Laboratory experiment
Normal digital simulation
Improved digital simulation
Total time spend
Weeks
Hours
Hours
Cost
Yes
No
No
Image quality
High
Medium
High
Noise
Yes
No
No
Obtaining disposed surface
Hard
Easy
Easy
Proximity to reality
Agree
Approximate
Disagree
Acquisition of displacement results
Complex for limited position
Easy for arbitrary position
Easy for arbitrary positionAn improved solution is presented to simplify the acquisition of digital images and to overcome its limitations. This 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.The gray intensity distribution of the normal Gaussian speckle image can be computed as follows:(17)fx,y=int∑k=1sI0exp−x−xk2+y−yk2R2,where I0 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 [xk, yk] represents the center location of the kth speckle.More than 2.25 × 1010 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. This 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. Thus, an additional condition can be added to equation (17) to reduce the computation as follows:(18)x−xk2+y−yk2≤3R.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). The deformation of the model can be obtained asU=uuxuy and V=vvxvy, and the gray intensity distribution of the deformed image can be derived as follows:(19)f∗X+U,Y+V=fX,Y.It should be noted that the coordinates obtained using equation (19) are no longer integers. This 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.


### 3.2. 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). The inclined precrack starts at the middle of the vertical edge at an angle of 20° to the horizontal with a length of 30 mm. The displacement results are obtained using Abaqus and presented as shown in Figure 5(b). The reference image with a resolution of 1000 × 2000 pixels is obtained by applying equations (17) and (18) using Matlab (The MathWorks, Natick, MA, USA), and the deformed image is generated by substituting the FEA displacements results into equation (19). These 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.Figure 5
The FEA model and images used for testing: (a) design of the test model; (b) total displacement result; (c) reference image; (d) deformed image.The 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. The 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 Figure6.Figure 6
The reconstructed displacement fields using standard DIC and the proposed technique: (a) horizontal displacement using the proposed technique; (b) horizontal displacement using standard DIC; (c) vertical displacement using the proposed technique; (d) vertical displacement using standard DIC.
(a)
(b)
(c)
(d)It is evident from Figure6 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.The correlation coefficients of the proposed discontinuous DIC and standard DIC were investigated for the ROI, and the results were collected and reorganized. The distribution of the correlation coefficient probability is shown in Figure7(a), while the correlation coefficient values along the crack are shown in Figure 7(b).Figure 7
The correlation coefficient (a) distribution and (b) values along the crack.
(a)
(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 Figure6. 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.The reconstruction of the horizontal and vertical displacements of the crack is displayed in Figures8(a) and 8(b), respectively, along with the results obtained for the standard DIC and theoretical crack displacements acquired using FEA. The results indicate that the proposed method achieves the design goal.Figure 8
Reconstructions of crack displacements: (a) horizontal displacement; (b) vertical displacement.
(a)
(b)A quantitative evaluation of the calculated displacements is adopted using the mean bias error and the standard deviation error as criteria. The mean bias error of the displacement result is defined as follows:(20)δ=1N∑Nu−ureal,where u and ureal denote the measured and real displacements of each point in the ROI, respectively.The standard deviation error of the measured displacement is expressed as follows:(21)σ=1N−1∑Nu−ureal2.The results of different mean bias errors and standard deviation errors calculated using the proposed method and standard DIC are listed in Table2. The 20-pixel vicinity of the crack is defined as the area around the discontinuities to demonstrate the capability of appropriately addressing this region. The 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.Table 2
Mean bias error and standard deviation error obtained using different methodsa.
Result range
The proposed discontinuous DIC
Standard DIC
δ
u
σ
u
δ
v
σ
v
δ
u
σ
u
δ
v
σ
v
Full ROI
0.0105
0.0213
0.0143
0.0175
0.3472
0.5142
0.4543
0.5398
Vicinity of crack
−0.1302
0.1916
−0.2133
0.1583
−3.7218
4.2046
−5.4437
4.9795
Continuous area
0.0198
0.0177
0.0227
0.0169
0.0271
0.0196
0.0218
0.0215
aIn pixels.


### 3.3. 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 Figure9(a). The dimensions of 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).Figure 9
(a) Specimen configuration; (b) loading process.
(a)
(b)The 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. The 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 lengthslc 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.Figure 10
The tensile test and deformed ROIs.These ten deformed images along with the reference image are substituted in the proposed method and the standard DIC, respectively. The preliminary run with the proposed method yields displacements fields, and the results from crack stagelc = 10 mm are selected and displayed, as shown in Figure 11. Similar results for the standard DIC are also presented for comparison. The 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.Figure 11
The displacement reconstruction results for the proposed technique along with the standard DIC: (a) horizontal displacement obtained using the proposed technique; (b) horizontal displacement obtained using the standard DIC; (c) vertical displacement obtained for the proposed technique; (d) vertical displacement obtained for the standard DIC.
(a)
(b)
(c)
(d)The reconstruction of the crack in each crack stage is calculated in both the horizontal and vertical directions. The results for three representative crack stages are selected for crack lengthlc 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).Figure 12
The reconstructions of crack displacement in the horizontal direction: (a)lc = 2 mm, (c) lc = 6 mm, and (e) lc = 10 mm, and vertical direction: (b) lc = 2 mm, (d) lc = 6 mm, and (f) lc = 10 mm.
(a)
(b)
(c)
(d)
(e)
(f)It is evident from Figure12 that the proposed method achieves high accuracy for all the crack processes. The fluctuation of the absolute error for each test point increases very little with crack propagation because the subset restore model is applied.The 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:(22)σs=σu2+σv2.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 Figures13(a) and 13(b). It is evident that crack propagation affects the accuracy of the results, but this effect is limited. The correlation coefficients increase by approximately 0.1, while the resultant standard deviation errors increase by approximately 0.08 pixels.Figure 13
The effect of crack propagation on (a) correlation coefficients and (b) resultant standard deviation error.
(a)
(b)


## 3.1. Acquisition of Test Images

Digital images for validation of the DIC method were mainly acquired from laboratory experiment [10–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.Table 1
Comparison of different image acquisition methods.
Description
Acquisition methods for digital images
Laboratory experiment
Normal digital simulation
Improved digital simulation
Total time spend
Weeks
Hours
Hours
Cost
Yes
No
No
Image quality
High
Medium
High
Noise
Yes
No
No
Obtaining disposed surface
Hard
Easy
Easy
Proximity to reality
Agree
Approximate
Disagree
Acquisition of displacement results
Complex for limited position
Easy for arbitrary position
Easy for arbitrary positionAn improved solution is presented to simplify the acquisition of digital images and to overcome its limitations. This 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.The gray intensity distribution of the normal Gaussian speckle image can be computed as follows:(17)fx,y=int∑k=1sI0exp−x−xk2+y−yk2R2,where I0 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 [xk, yk] represents the center location of the kth speckle.More than 2.25 × 1010 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. This 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. Thus, an additional condition can be added to equation (17) to reduce the computation as follows:(18)x−xk2+y−yk2≤3R.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). The deformation of the model can be obtained asU=uuxuy and V=vvxvy, and the gray intensity distribution of the deformed image can be derived as follows:(19)f∗X+U,Y+V=fX,Y.It should be noted that the coordinates obtained using equation (19) are no longer integers. This 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.


## 3.2. 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). The inclined precrack starts at the middle of the vertical edge at an angle of 20° to the horizontal with a length of 30 mm. The displacement results are obtained using Abaqus and presented as shown in Figure 5(b). The reference image with a resolution of 1000 × 2000 pixels is obtained by applying equations (17) and (18) using Matlab (The MathWorks, Natick, MA, USA), and the deformed image is generated by substituting the FEA displacements results into equation (19). These 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.Figure 5
The FEA model and images used for testing: (a) design of the test model; (b) total displacement result; (c) reference image; (d) deformed image.The 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. The 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 Figure6.Figure 6
The reconstructed displacement fields using standard DIC and the proposed technique: (a) horizontal displacement using the proposed technique; (b) horizontal displacement using standard DIC; (c) vertical displacement using the proposed technique; (d) vertical displacement using standard DIC.
(a)
(b)
(c)
(d)It is evident from Figure6 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.The correlation coefficients of the proposed discontinuous DIC and standard DIC were investigated for the ROI, and the results were collected and reorganized. The distribution of the correlation coefficient probability is shown in Figure7(a), while the correlation coefficient values along the crack are shown in Figure 7(b).Figure 7
The correlation coefficient (a) distribution and (b) values along the crack.
(a)
(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 Figure6. 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.The reconstruction of the horizontal and vertical displacements of the crack is displayed in Figures8(a) and 8(b), respectively, along with the results obtained for the standard DIC and theoretical crack displacements acquired using FEA. The results indicate that the proposed method achieves the design goal.Figure 8
Reconstructions of crack displacements: (a) horizontal displacement; (b) vertical displacement.
(a)
(b)A quantitative evaluation of the calculated displacements is adopted using the mean bias error and the standard deviation error as criteria. The mean bias error of the displacement result is defined as follows:(20)δ=1N∑Nu−ureal,where u and ureal denote the measured and real displacements of each point in the ROI, respectively.The standard deviation error of the measured displacement is expressed as follows:(21)σ=1N−1∑Nu−ureal2.The results of different mean bias errors and standard deviation errors calculated using the proposed method and standard DIC are listed in Table2. The 20-pixel vicinity of the crack is defined as the area around the discontinuities to demonstrate the capability of appropriately addressing this region. The 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.Table 2
Mean bias error and standard deviation error obtained using different methodsa.
Result range
The proposed discontinuous DIC
Standard DIC
δ
u
σ
u
δ
v
σ
v
δ
u
σ
u
δ
v
σ
v
Full ROI
0.0105
0.0213
0.0143
0.0175
0.3472
0.5142
0.4543
0.5398
Vicinity of crack
−0.1302
0.1916
−0.2133
0.1583
−3.7218
4.2046
−5.4437
4.9795
Continuous area
0.0198
0.0177
0.0227
0.0169
0.0271
0.0196
0.0218
0.0215
aIn pixels.


## 3.3. 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 Figure9(a). The dimensions of 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).Figure 9
(a) Specimen configuration; (b) loading process.
(a)
(b)The 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. The 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 lengthslc 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.Figure 10
The tensile test and deformed ROIs.These ten deformed images along with the reference image are substituted in the proposed method and the standard DIC, respectively. The preliminary run with the proposed method yields displacements fields, and the results from crack stagelc = 10 mm are selected and displayed, as shown in Figure 11. Similar results for the standard DIC are also presented for comparison. The 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.Figure 11
The displacement reconstruction results for the proposed technique along with the standard DIC: (a) horizontal displacement obtained using the proposed technique; (b) horizontal displacement obtained using the standard DIC; (c) vertical displacement obtained for the proposed technique; (d) vertical displacement obtained for the standard DIC.
(a)
(b)
(c)
(d)The reconstruction of the crack in each crack stage is calculated in both the horizontal and vertical directions. The results for three representative crack stages are selected for crack lengthlc 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).Figure 12
The reconstructions of crack displacement in the horizontal direction: (a)lc = 2 mm, (c) lc = 6 mm, and (e) lc = 10 mm, and vertical direction: (b) lc = 2 mm, (d) lc = 6 mm, and (f) lc = 10 mm.
(a)
(b)
(c)
(d)
(e)
(f)It is evident from Figure12 that the proposed method achieves high accuracy for all the crack processes. The fluctuation of the absolute error for each test point increases very little with crack propagation because the subset restore model is applied.The 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:(22)σs=σu2+σv2.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 Figures13(a) and 13(b). It is evident that crack propagation affects the accuracy of the results, but this effect is limited. The correlation coefficients increase by approximately 0.1, while the resultant standard deviation errors increase by approximately 0.08 pixels.Figure 13
The effect of crack propagation on (a) correlation coefficients and (b) resultant standard deviation error.
(a)
(b)


## 4. Conclusion

This article presents an improved DDIC method based on the subset restore model and subset adaptation algorithm to address the discontinuous problems in standard DIC. The 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. The 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. The 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. The 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.



---

*Source: 5271962-2020-05-05.xml*