Initial Displacement Estimation Based on ESA in Improved Digital Volume Correlation and Its Application

In order to improve the algorithm efficiency in point matching and reduce the error of displacement field in deformation calculation, an Improved Digital Volume Correlation (IDVC) method is proposed. Aimed at the problem that the accuracy of initial displacement estimation of DVC has a great influence on the convergence calculation accuracy by the sub-voxel measurement algorithm, the coordinates of the estimated points in deformation are quickly matched by the Explosion Search Algorithm (ESA), with the displacement vectors obtained. Then, the non-physical fluctuation of displacement field and the wrong strain concentration are suppressed in regularization. Finally, the deformation information obtained by correlation search is substituted into the iterative least squares algorithm to attain the sub-voxel precision displacement. The results of the simulation experiment and concrete uni-axial compression experiment show that the computational efficiency by the proposed method is 20.44% higher than that by the traditional DVC method, and that the accuracy level of displacement field measurement reaches 0.001 voxels.


I. INTRODUCTION
High resolution computed tomography (CT) enables people to obtain three-dimensional volume images of micro-scale materials, thus providing an important way to analyze the mechanical properties of materials. As a novel and effective technique for measuring internal deformation fields of non-transparent solid materials, metal materials, and some internal tissues of organisms, the Digital Volume Correlation (DVC) detection technology has been applied to many fields [1]- [5] and has received growing attention. In 1999, Bay et al. [6] and Bay [7] proposed the DVC, which obtains the internal displacement and strain field by CT scanning and analyzing the image deformation of an object. As a result, it has become a major experimental means of characterizing the whole field deformation of an object. Similar The associate editor coordinating the review of this manuscript and approving it for publication was Danping He .
to the principle of DIC (Digital Image Correlation), DVC determines the displacement vector (unit: voxels) by tracking the correlation points of a matched object. The accurate estimation of initial displacement is the premise of obtaining reliable measurement results of a deformation field. The existing estimates of initial deformation include the coarsefine search method proposed by Sutton et al. [8]. To be specific, the search area is roughly searched with large step size to find the extremum position of a correlation coefficient, and then the search step size is reduced for fine search. The genetic algorithm proposed by Pilcha et al. [9], Ma et al. [10] is based on the mechanism of biological selection and the principle of randomized probability statistics to search the optimal solution. Chiang [11], Chen and Chiang [12], and Chen et al. [13] used FFT to perform correlation operations in the frequency domain. Most of the existing applications of initial displacement estimation are planar and require little computation. In the case of large deformation measurement inside an object, the displacement of a target point matched by traditional DVC method will rise sharply in the process of deformation, resulting in an increase in the number of iterations easy to fall into the local optimum.
Based on the Explosion Search Algorithm (ESA), an improved digital volume correlation (IDVC) method for estimating the initial deformation value will be proposed in this paper. ESA performs the steepest descent search in the process of point matching in deformation according to the intelligent search strategy. In the implementation of the algorithm, an adaptive factor is added to make the search area adaptively enlarged and achieve global search. The accurate initial value of deformation estimation is obtained, by ESA algorithm as the initial value of the iterative least squares algorithm. The quality of DVC results depends on the details (high frequency information), anomalies (including intrinsic noise, artifacts, etc.) and gray interpolation errors in the source image. We will introduce regularization to smooth the displacement field. In order to eliminate the non-physical strain concentration, the strain value will be calculated by Point Least Squares (PLS), and the high frequency noise in displacement data will be effectively resisted by local fitting.

II. BASICPRINCIPLE OF DVC
Speckle is used as a deformation information carrier of a measured object, and a natural material usually had a randomly distributed natural speckle texture. According to the three-dimensional volume image of the object, the DVC method calculated the three-dimensional displacement vector by matching the corresponding speckle feature points before and after the deformation of the object, and repeated the process to obtain the deformed three-dimensional displacement field. The strain field is extracted from the three-dimensional displacement field through an appropriate difference algorithm. The accuracy of the strain field calculation depended largely on the precision of the displacement field measurement. Therefore, the key to the DVC method is how to accurately acquire the three-dimensional displacement field of the object.
The subset containing n = (2M + 1) 3 voxels are selected in the three-dimensional volume image, and the target subset with the greatest similarity to the reference gray level is searched by coarse search combined with ESA to obtain the displacement of integer voxels u 0, v 0 ,w 0 . Zero-mean Normalized Cross Correlation (ZNCC) is introduced as an evaluation criterion [15] (1), as shown at the bottom of this page.
where Re(x i , y i , z i ) and D(x i + u 0 , y i + v 0 , z i + w 0 ) are the gray levels of point i in the reference subset before and after deformation. Re m and D m represent the average gray levels of the reference subset and the target subset respectively. ϑ zncc ∈ [−1, 1]. ϑ zncc = 0 means that the reference subset is completely unrelated to the target one. ϑ zncc = −1, which indicating a negative correlation between the two.

III. INITIAL DISPLACEMENT ESTIMATION BASED ON ESA
The ESA is inspired by bomb explosion and characterized by fast convergence, parallel search, high precision compared with genetics, particle swarm and the like. Firstly, explosion points with radii r t = (r t1 , r t2 , . . . , r tD ) are initialized. The initial position coordinate of the subset is set as the start point of search. The position change of an explosion point will cause its radius to be renewed correspondingly. The renewal of radius r can be divided into the following two cases: 1. When the position of the explosion point is updated, the formula for explosion radius r is as follows: where r tmax and r tmin are the maximum radius and minimum radius respectively, T represents the current iteration number, 96228 VOLUME 8, 2020 and MaxIt is the maximum number of iterations during the current ESA calculation. 2. If no new explosion point more suitable for the original explosion point can be found in this search scope, the location of explosion point will not be updated. The formula for explosion radius r is as follows: where ∂ is the scale factor, and ∂ > 1. In practical application, it is generally selected according to experience. In the process of ESA implementation in this paper, ∂ = 1.5.
According to the size of fitness value, all explosion points are reordered accordingly, and the explosion points with poor fitness values [µn] (with µ standing for the migration ratio, n for the total number of subset, and [ * ] for the rounding function) are shifted as: where X represents the location of the updated explosion point.
After the ESO and MO are executed, the operation is performed on MuO. Before the mutation operator is executed, a mutation probability P M is required.. The mutation probability in this paper is set to be 0.01. A random number is generated between [0,1]. If the value is less than P M , the MuO operation is performed on the explosion point. For X = (x l , x 2 , . . . x D ), the variation formula in each dimension contained in the explosion point is where σ = (σ 1 , σ 2 , . . . , σ D ) = 0.1 * (X max − X min ). X min and X max are the two corresponding pieces of boundary information on the search interval, whose dimension is same as that of the search space. In this paper, only three-dimensional space is involved. The calculation process of initial displacement estimation is as follows: n subsets are initialized as explosion points in the search space, the similarity values of gray distribution before and after deformation are determined, and the point with the smallest similarity value of gray distribution (min{−ϑ zncc }) is stored as Gbest.
According to the results of ESO, if the fitness (similarity) of the explosion point is less than the fitness of the original one, the position of the explosion point will be updated, and its radius would be updated according to formula (2). If the position of the explosion point stay unchanged, the search radius r will be updated according to formula (3).
The fitness values of all explosion points are calculated in parallel, with the corresponding points with poor fitness migrated according to formula (4).
The mutation probability P M and the size of the generated random number are determined, and the mutation operation is performed in line with formula (5) if necessary. t = t + 1, if t < MaxIt, then step 2 is followed. Otherwise, the iteration ends and the result will be output.

A. GLOBAL OPTIMIZATION
In order to optimize the susceptibility to local deformation differences [20], the global minimization function min {−ϑ zncc } is replaced by the minimum mean value of all subsets: where n is the sum of subsets, and is the corresponding displacement vector after deformation. It should be noted that the average deformation correlation of subset is different from the whole volume correlation.

B. REGULARIZATION
In the hoping of suppressing the inherent noise and alleviating the potential problems caused by the lack of details on the source image, regularization is added to the optimization process, and curvature penalty is added to the displacement field in the objective function: where e is the subset area, D is the reference volume image vector parameter, and the integral part of the equation used the 3 * 3 * 3 Gauss integral. The magnitude of the penalty depended on the displacement vector parameter. The objective function containing regularization is rewritten as follows: where θ denotes the penalty factor in regularization for adjusting the amount of regularization. One can attempt to relate the magnitude of this penalty factor to the wave length of spurious oscillations in the resulting displacement field, see the work of Leclerc et al. [23]., or minimize the condition number of the Hessian matrix of the optimization, similar to the work of Barber and Hose [24]. Regularization is a compromise measure for accuracy and precision of strain field. =[U,V,W] is the corresponding displacement vector obtained from the initial displacement estimation based on ESA.
C. SUB-VOXEL DISPLACEMENT CALCULATIO e Aimed at improving the accuracy of displacement measurement, the iterative least squares fitting [16] is used to solve the problem. Given that the gray contrast and shape of the target image might change except for the position update of the center point of the reference subset after deformation, a linear gray change model is adopted to reduce the displacement measurement error caused by the gray change of the target subset.
D(x i , y i , z i ) = Q * Re(x i , y i , z i ) + G, i = 1, 2, . . . , n (9) where Q and G represented the linear change of gray level and the offset value of the object before and after deformation, VOLUME 8, 2020 respectively. The gray level change model showed stronger anti-jamming ability. It is completely equivalent to the ϑ zncc function which evaluated the similarity degree.
Considering the possible rigid body translation, rotation, tensile deformation and shear deformation of the target subset, the linear displacement function [17] is used to express the parameter changes after deformation.
where x, y, z are the distances from point i (x i , y i , z i ) in the reference subset to the center point (x 0 , y 0 , z 0 ) respectively. u 0 , v 0 , w 0 , and u, v, w here are the integers and sub voxel level displacement in three-dimensional coordinates, respectively. u x,y,z , v x,y,z , w x,y,z are the displacement gradients of the subset image.
The objective function of point i is defined as L i (υ). υ is the parameter vector including sub-voxel displacements, displacement gradients, gray shift factors and scale change factors. The Taylor expansion of the objective function in D(x i + u 0 , y i + v 0 , z i + w 0 ) is as follows.
where υ N+1 , υ N are solutions after the N + 1 th and N th iterations. For a subset containing n = (2M + 1) 3 voxels points, the (N + 1) th parameter vector is obtained by iterative least squares calculation: In the above function: D xn x n D xn y n D xn z n D xn D yn x n D yn y n D yn z n D yn D zn x n D zn y n D zn z n D zn In each iteration process, the gray value and gray gradient of the target point position needed to be obtained in advance. The bicubic interpolation method with high accuracy is employed. The sub-voxel precision displacement is secured by calculating equation (12) repeatedly until the algorithm satisfied the preset convergence condition. The convergence condition of the proposed algorithm is that the displacement error fluctuation of two adjacent iterations should be less than or equal to 10 −3 voxels.

D. STRAIN CALCULATIO
To reduce the influence of error in displacement data on strain calculation, the point least squares fitting method (PLS) [14] is used to obtain strain values. High frequency noise in displacement data is effectively suppressed by local fitting.
where u, v, and w are the displacement data obtained by ESA. A i=0−3 , B i=0−3 , C i=0−3 are polynomial coefficients to be determined: Through these six Cauchy strain components, the principal strain, equivalent strain and volumetric strain can be obtained.
where ε 1 , ε 2 , and ε 3 are the first, second and third principal strains respectively. The positive strains of ε x , ε y , and ε z are along the directions of x, y, and z, respectively.

IV. VALIDATION OF IDVC ALGORITHM A. VIRTUAL DEFORMATION OF THREE-DIMENSIONAL VOLUME IMAGE BY COMPUTER SIMULATION
3D speckle volume images generated by computer simulation are used to verify the accuracy and computational performance of the algorithm. The computer configuration is an Intel i5-7500 CPU 3.40 GHz processor with an 8.00GB memory, and the software version of MATLAB is R2016b. The speckle pattern inspired by Zhou and Kenneth [18] is applied to simulating surface deformation process of objects, and the following formula are used for 3D speckle image in the process of deformation: where int denotes the rounding function, s represents the total number of speckle particles in the volume image generated by computer simulation, R is the radius of the speckle particles in the computer simulation volume image, and E 0 i is the random distribution intensity of the center of the first speckle particle. Besides, (x i , y i , z i ) represents the central position of the first spatial speckle particle before deformation, and (x i , y i , z i ) corresponds to the central position after deformation. The parameters used in this simulation experiment are as follows: s = 12000, R = 1, E 0 i = 2, and the size of speckle volume image reached 100 * 100 * 100 voxels.
The rigid body translation experiment is simulated by the computer simulation. In this experiment, the reference volume image is translated into 0.1 voxels in turn along Z axis, and the other deformation parameters (tension, compression & shear deformation) are set to be 0 to generate 10 consecutive sub-voxel translation volume images. In the calculation of sub-voxel displacement using IDVC, four different sizes of subset 11 3 , 21 3 , 31 3 and 41 3 are used respectively.  FIGURE 4.(a) shows that the measured displacement is basically the same as the pre-added sub-voxel displacement. FIGURE 4.(b) shows that the displacement is approximately sinusoidal. This error distribution is caused by sub-voxel interpolation algorithm [19] and the maximum mean error occurred at the displacements of 0.2 voxels and The displacement U in x direction is set to be 4.25 voxels, the displacement V in y direction is 3.26 voxels and W in z direction is 2.50 voxels. Other deformation parameters are set to be 0. The calculation results are compared with those based on the traditional DVC method, as shown in Compared with the gradient-based digital volume correlation method [22], the proposed method is more efficient but less precise. The average computing time is 2.725s when n = 31 3 voxels (the best size of subset in this paper), and the computing efficiency is improved by 20.44%.
We used coarse-fine search, fourier-transform-based cross correlation algorithm and the proposed algorithm based on ESA to match the center data points of random subsets in the reference image and the target image. The following table lists the calculation efficiency of point matching of different algorithms, the comparison results of measurement mean error and standard deviation. The deformation parameters are set as follows: U = 3.78, V= 3.82, W = 2.75.
According to the data analysis in FIGURE 5 & Table 3, it is concluded that the efficiency of Fourier-transform-based cross correlation is the highest, but the error of displacement measurement is the largest. This is because the algorithm can VOLUME 8, 2020 only solve the integer displacement in the matching process. It is very time-consuming for coarse-fine search, which puts forward higher requirements for DVC algorithm with more data points. The algorithm proposed in this paper takes a trade off between accuracy and calculation efficiency. Combining the advantages and disadvantages of the above three algorithms, the DVC algorithm based on ESA proposed in this paper has better application value in practice.

B. UN-AXIAL COMPRESSION EXPERIMENT OF CONCRETE
To verify the reliability and practical performance of the algorithm in real deformation, the uni-axial compression experiment of concrete is conducted [21]. The aggregate is made up of pebbles with a diameter of about 0.5-2cm. Standard  concrete cylinder specimens (height/diameter = 12cm/6cm) are chosen from these raw materials. The selected equipment used a Marconi M8000 spiral CT scanner from image center equipment of Xi'an Central Hospital. The concrete    before the peak compressive strength, and the concrete sample is scanned by CT, which is continuously circulated until the sample is destroyed. CT and loading equipment are described in FIGURE.6.
Due to the influence of the loading plate, the image ends are greatly affected by image distortion and artifacts. As a consequence, only the middle part is imported into the AVIZO software for 3D reconstruction.
The mortar threshold is 1103-1583 HU and the aggregate threshold is 1584-3095 HU based on the CT histogram of concrete. The three-dimensional volume image reconstructed by AVIZO is as follows: It could be inferred that MATLAB is slightly inadequate in three-dimensional display. In order to obtain detailed deformation information of materials, the concrete sections under different static compression loads are intercepted from the deformed volume images for analysis. As indicated by the transparent part of materials, no displacement is generated under loading conditions, though cracks are generated along the edge of aggregate and mortar. The specimen is loaded along Z-axis direction. In other words, axial loading is dominant, and the compression resulted in deformation and failure.  Judging from FIGURE 10 (a) and (b), the boundary between positive and negative displacements appeared. It could be seen from the equivalent strain nephogram that the local high strain region came into being before the failure of the specimens, and developed into macro cracks with the increase of load. n = 31 3 voxels is selected as the subset, and the crosssection strain value measured by the algorithm in this paper is |−0.06|, which basically coincided with 0.058 measured by the commercial software Avizo. It is shown that the algorithm remained reliable in practical application. The application of IDVC algorithm in concrete compression experiment could provide some technical means and theoretical basis for revealing the local generation and development of deformation before macro crack formation.

V. CONCLUSIONS
This paper has proposed a new IDVC method, which combines ESA and digital volume correlation method together to overcome the problem that the traditional DVC method is greatly affected by the initial value of iteration: • Through the parallel computing feature of ESA, the coordinates of the estimated points before and after deformation have been quickly matched, and then the displacement vectors have been obtained, improving the computational efficiency of DVC algorithm.
• Curvature penalty (regularization) has been used to suppress the non-physical fluctuation of displacement field and wrong strain concentration. The performance of the proposed IDVC method has been evaluated by computer simulation speckle experiment. The simulation results have shown that the proposed algorithm maintains the accuracy level of displacement field calculation and improves the calculation efficiency. The results of strain measurement in uni-axial compression test experiment have been basically consistent with that calculated by the commercial software, thereby verifying the feasibility of the algorithm in practical measurement.
• Regarding the IDVC method proposed in this paper, regular cubic subsets have been used to obtain the displacement field, with the internal deformation of subsets neglected. In actual calculation, the internal deformation of subset scan is considered to further improve the accuracy of the algorithm. At the same time, the complexity of the algorithm will be increased accordingly, thus increasing the calculation time.
Traditional initial displacement estimation like coarse-fine search has such characteristics that all search fields of the deformed image are searched point by point to find out the maximum position of the correlation coefficient C. Then, the fine search is carried out on the basis of the coarse search. It needs sub-pixel reconstruction near the integer pixel position of the deformed image. The coarse-fine search is easy to program, but it is quite time-consuming in search. The parameters obtained by Fourier-transform-based cross correlation algorithm are only zero order deformation vector parameters [u, v, w],which are difficult to deal with large deformation or complex non-uniform deformation. So this paper proposed this algorithm to solve the time-consuming problem of too many matching points in DVC, and simplifies the process of computer processing to find the extreme value of correlation coefficient. It is also suitable for large deformation. In the aspect of adaptive extended search, the algorithm proposed in this paper has some limitations. For example, when there are too many matching data points, bad point pairs need to be eliminated manually, which is not friendly for DVC processing flow. It requires continuous improvement of the algorithm