Accurate and simple digital volume correlation using pre-interpolation

Existing incremental digital volume correlation methods can reduce the number of errors introduced by interpolation calculations in the inverse-compositional Gauss–Newton (IC-GN) algorithm iteration. However, the accuracy of these existing methods is insufficient for some conditions as the curve-fitting method has high computational efficiency but lacks accuracy. A simple pre-interpolation method is proposed to improve the accuracy and computational efficiency of digital volume correlation. First, the pretreatment of a deformed volume image is calculated by the cubic spline interpolation method with the most often chosen interpolation step of 1/2 sub-voxel. Next, the pre-interpolation is calculated only once and the block calculation techniques solve the memory problem. Then, the reference sub-volume in the updated reference volume image is translated into the nearest half-integer voxel position instead of the integer voxel position or other sub-voxel positions. The pre-interpolation method is applied to both the IC-GN algorithm and the curve-fitting method. Experimental results show that the maximum mean bias error and the maximum standard deviation of the improved IC-GN algorithm are reduced by 34% and 75%, respectively. The improved curve-fitting has better accuracy and computational efficiency than the IC-GN algorithm under small strain and the curve-fitting method can achieve about 3.2 times speedup than the IC-GN algorithm.


Introduction
The digital volume correlation method was developed from mature two-dimensional digital image correlation (DIC) technologies by Bay et al [1]. Since 1999, the DVC methodhas been widely used to quantify the deformation field of materials under external loading [2]. DVC can contribute to not only a better understanding of the mechanical behavior of materials from the perspective of a microstructural study [3], but also to the verification of the results of 3D finite element simulations [4]. In the past 20 years, many researchers have been committed to improving the accuracy and computational efficiency of DVC algorithms [5][6][7][8][9]. The inverse-compositional Gauss-Newton (IC-GN) algorithm is an advanced matching algorithm [10] that can accurately estimate displacement fields when large deformation occurs.
Although the Hessian matrix only needs to be calculated once in this method, its iterative calculation is still time-consuming. Therefore, Pan et al [6] used the interpolation coefficient lookup table approach to reduce the amount of calculation when using the IC-GN algorithm in order to improve the computational efficiency; they obtained an improvement in computational speed, but the approach was memory-consuming. The fast Fourier transform algorithm was also introduced into the DVC algorithm for cases of small deformation to improve computational efficiency [7]. In recent years, with the rapid development of graphics processing unit (GPU) hardware and software, GPUs have been gradually employed to accelerate the DVC algorithm and computational speed was significantly improved compared to that of a central processing unit (CPU) [11]. A simple and effective incremental DVC method was proposed by Wang et al [8] to reduce redundant interpolation calculations and errors introduced by interpolation calculations in IC-GN iterations. The deformation of different regions in a real material was calculated by Wang et al [9] using the self-adaptive DVC approach; the large deformation is calculated using a small sub-volume size while the small deformation is calculated using a large sub-volume size.
However, there is still space for improvement in the accuracy and computational efficiency of these algorithms. For example, incremental algorithms that perform well in terms of computational speed may have problems with accuracy and system stability. The self-adaptive DVC approach has the disadvantage of programming complexity. The accuracy and stability of the curve-fitting algorithm are relatively poor compared to the IC-GN algorithm; however, its computational efficiency is outstanding [12]. The curve-fitting algorithm would still be a good choice if its accuracy could be improved.
In this work, a simple and effective pre-interpolation method is proposed for the DVC algorithm. First, a 1/2 subvoxel interpolation calculation is performed in advance, using the proposed pre-interpolation method, on a deformed volume image matrix. Then, the sub-voxel displacement is estimated by the DVC incremental algorithm and the curve-fitting algorithm. The pre-interpolation calculation only needs to be carried out once, which reduces the error caused by intensity interpolation of renewed reference correlation points. Finally, the efficiency and capability of the proposed improved DVC approach is demonstrated in practical applications.

IC-GN: large strain
Interpolation calculation is a time-consuming task in DVC calculations; it significantly influences the computational speed. As mentioned earlier, Pan et al [6] effectively reduced redundant calculations by using the interpolation coefficient lookup table approach to perform sub-voxel interpolation; however, the deformed sub-volume must be recalculated in each iteration. Subsequently, incremental deformation was employed by Wang et al [8] to solve the interpolation problem wherein each incremental position of the deformation is mapped to the nearest integer point. Their simulation results show that the calculation efficiency is increased by approximately 2.5 times the initial value and the error caused by intensity interpolation of renewed reference correlation points is reduced. However, according to the results of this algorithm, in cases such as a sub-voxel displacement calculation, the random errors cannot be ignored when the sub-voxel displacement is about 0.5 voxels (assuming a sub-volume size of 41 × 41 × 41 the mean bias error is 0.003). Moreover, the systematic error of the incremental algorithm cannot be ignored when the sub-voxel displacement is about 0.25 voxels or 0.75 voxels (assume a sub-volume size is 41 × 41 × 41, stand deviation error is 0.007). The accuracy of the incremental algorithm is far lower than that of the normal IC-GN algorithm in some cases. In this study, the combination of the pre-interpolation method and the incremental algorithm is proposed in order to solve these problems.
A basic 3D IC-GN algorithm flowchart is illustrated in figure 1. For each measurement point, the robust zero-mean normalized sum of squared difference criterion can be optimized through the sub-voxel displacement determined by the IC-GN algorithm. It can be seen from equation (1) that the deformation submatrix needs to be constantly updated in each iteration; the corresponding mathematical operation is the interpolation calculation (more details regarding the DVC algorithm can be found in [6]). (1) where f (x) and g(x) are the gray intensity values at point x= (x,y,z) T in the reference and target sub-volumes; ξ = (∆x,∆y,∆z) T is the local coordinates of integer voxel points in the reference sub-volume; f m and g m represent the mean  intensity values of reference and deformed sub-volumes; p is the linear deformation vector; W(ξ;p) is the linear displacement mapping function used to describe the deformation of the target sub-volume; ∆p is the incremental deformation vector; W(ξ; ∆p) is the incremental mapping function exerted on the reference sub-volume.
Since the general CT volume image matrix exceeds 1000 × 1000 × 1000 voxels and reduces the memory requirement, the pre-interpolation is usually set to 1/2 sub-voxel interpolation. In this work, the interpolation coefficient α is employed to represent the number of interpolation segments between two points; α is a positive integer. In figures 2(a) and (b), the pre-interpolation calculation when α = 2 is briefly explained. For example, as shown in figures 2(a) and (b), when taking a matrix with 2 × 2 × 2 voxel size in the deformed volume image and performing cubic spline interpolation, 19 interpolation points are added at the 1/2 grid nodes. Then, the 2 × 2 × 2 voxels matrix is converted into a 3 × 3 × 3 voxels matrix and the number of nodes changes from 8 to 27.
The interpolated deformation matrix is 2 3 = 8 times larger than the original matrix when α = 2. Although it needs to occupy a massive amount of memory, this problem can be overcome using block calculation techniques (as shown in figure 3). According to the memory of a computer, the original matrix can be partitioned in different blocks. s is the maximum displacement, M is the sub-volume size, and the deformed block volume image is bigger than the reference block volume image.
A schematic illustrating the brief flow of the DVC increment algorithm when α = 2 is shown in figure 2(c). According to the firstst-order shape function and the position relationship between reference and target subvolume is: where (x, y, z) are the coordinates of the updated reference calculation point; (x 0 , y 0 , z 0 ) are the coordinates of the reference calculation point; ∆x, ∆y, ∆z are the local coordinates; u,v,w are three displacement components; u x , u y , u z , v x , v y , v z , w x , w y and w z are nine components of the displacement gradient.
For the incremental displacement vector of the updated measurement point of the incremental IC-GN method, the displacement (x, y, z) are translated to the nearest integer point positions, thereby avoiding interpolation calculation and improving the calculation speed. The incremental displacement vector of the proposed method needs to be modified on the basis of the incremental IC-GN method, and the displacement (x, y, z) are translated to the nearest half-integer point positions: where x 1/2int = round(x), y 1/2int = round(y), and z 1/2int = round(z) are coordinates providing the location of the nearest half-integer voxel point of (x, y, z), and function round refers to rounding operation; 1/2 int denotes half-integer.
In the conventional IC-GN algorithm, the gradients of the intensity within the reference sub-volume ▽f (f x , f y , f z ) are gradients under integer voxel conditions. Since the deformed volume image matrix is interpolated, the differential formula needs to be corrected as: where α is the interpolation coefficient.

Curve-fitting: small strain
In real DVC strain analysis, the strain error is difficult to control if the deformation is too large. Generally, the average strain level of the material is small (usually less than 5%); however, the local strain may be large (for example, the shear strain of some shear bands may exceed 20%). It is time-consuming to fully utilize the IC-GN algorithm in terms of calculating small strains. In the iterative process of the IC-GN algorithm, the calculation of the multiplication and inversion of the remaining matrices is still very large, even though the Hessian matrix only needs to be calculated once. In the DIC algorithm, the computational speed of the curve-fitting algorithm is seven times that of the IC-GN algorithm [12,13]; the computational speed difference between these two algorithms will be more disparate in the DVC algorithm. Therefore, the curve-fitting method is more suitable for calculating small deformation regions. Three algorithms, including the curve-fitting method, gradient-based method and the Newton method, were compared in terms of calculation accuracy and stability [12]. The Newton method has the best calculation accuracy and stability; the gradientbased method is the second; the curve-fitting method is the worst (the accuracies and stabilities of the Newton method and IC-GN are almost the same [13]). The accuracy and stability of the curve-fitting method can barely meet the calculation requirements. In addition, the pre-interpolation calculation method described in section 2.1 is also applied to the curve-fitting algorithm in order to improve the accuracy and stability. First, the location of the extreme point of the zeronormalized cross-correlation (ZNCC) is calculated by the integer search algorithm; the ZNCC function is established at the extreme point. The ZNCC function can be set as C(x,y,z) and the general curve-fitting equation is a ternary quadratic function [14]: C(x, y, z) = a 0 + a 1 x + a 2 y + a 3 z + a 4 xy + a 5 xz + a 6 yz In some research, the least squares method is used to solve the equation parameters [12,15]; however, the least squares method involves complicated and time-consuming matrix operations. The explicit method is used for acceleration while the accuracy and stability of the existing explicit calculation method are insufficient [14]. Therefore, this method needs to be improved and solved using equation (13), which can be found in the appendix where detailed derivation processes are included.

Computational complexities of the proposed DVC method
In order to evaluate the computational efficiency of the proposed DVC method, we referred to the algorithm complexity calculation method in [6]. As the DVC algorithm is relatively complicated, it is difficult to accurately calculate each calculation step and there may be some differences in programming implementations by different people. Therefore, we focused on the main time-consuming steps in the calculation process, which are usually directly related to the subvolume size and the remaining calculation steps were not considered. Taking the traditional DVC algorithm as a reference, and then calculating the computational complexity of the improved incremental DVC algorithm and the curve-fitting algorithm, the statistical results are shown in table 1.
The main difference between the incremental DVC algorithm (or incremental IC-GN algorithm) and the traditional DVC algorithm is the interpolation calculation. There is no interpolation calculation in the process of incremental IC-GN algorithm and there are only a few calculations related to the extraction of deformed subvolumes. This difference becomes more significant as the number of iterations of the IC-GN algorithm increases.
Comparing with the improved incremental DVC algorithm, there are no operations about the Hessian matrix and subvoxel interpolation calculation in the process of curve-fitting algorithm. The sub-voxel displacement calculation of the curve-fitting algorithm only needs to be done once. The curvefitting algorithm reduces a lot of calculations compared to the incremental IC-GN algorithm in the case of small strain.

Accuracy of the proposed DVC method
In this work, a 100 × 100 × 100 voxel size referencing a 3D speckle pattern (R = 4 voxels, s = 12 000) was first generated, as shown in the inset of figure 4(a). Then, according to the Fourier shift theorem [16], ten pure rigid body translation volume images were generated as the deformed volume images. Along the z direction, the displacements range from 0 to 1 voxel (0.0, 0.1, 0.2, 0.3, …, 0.9, 1.0 voxel); the displacements of the x direction and y direction are zero. The random Gaussian noise with a mean value of zero and a variance of four was added to the previous 11 volume images. All DVC analyses were calculated by DVC software written using MATLAB 2018b language on a desktop computer (i7-6700 CPU 3.40 GHz and 8 GB RAM).
The result of the IC-GN algorithm represents the existing incremental DVC algorithm result when α = 1; this means that the incremental DVC algorithm proposed in [8] is a special case of the proposed method in this work. As shown in figures 4(a) and (b), for the existing incremental algorithm, the accuracy of the algorithm will be significantly reduced to be far from the accuracy of the conventional IC-GN algorithm when the sub-voxel displacement is in a certain range. Moreover, 1/2 sub-voxel interpolation is performed on the deformed volume image when α = 2; the result of IC-GN calculation shows that the insufficiency of the incremental algorithm can be effectively improved. Particularly, the maximum mean bias error is reduced by 34% and the maximum random errors are reduced by 75% when the sub-voxel displacement is between 0.4 and 0.6 voxels. The improvement of the mean bias error and random error is very significant. According to the trend of the curves, the period and amplitude of the mean bias error curves of α = 2 are about half of those of the mean bias error curves of α = 1. Since the interpolation calculation is completed before the formal DVC calculation, the interpolation calculation time can be almost ignored; however, the improvement of the accuracy and stability of the algorithm is particularly remarkable. Furthermore, the accuracy of the DVC algorithm is very close to that of the traditional IC-GN algorithm when α = 2. The specific comparison of algorithms can be found in [6].
As shown in figures 4(c) and (d), the mean bias error of the curve-fitting method is close to that of the IC-GN method when α = 1; the maximum error is about 0.003 voxels. However, the difference between the random errors of these two algorithms approaches one order of magnitude. Therefore, it is unreasonable to use the curve-fitting method to estimate the sub-voxel displacement directly. The curve-fitting calculation is much more accurate compared to the existing IC-GN algorithm when the deformed volume is pre-interpolated by a 1/2 sub-voxel (α = 2). Not only can the mean bias error of sub-voxel displacement be reduced nearly ten times, but the random errors can also be well controlled, which is very close to the random errors of the IC-GN algorithm and fully meets the requirements of a real calculation. The pre-interpolation method can effectively improve the calculation accuracy and stability of the curve-fitting method while retaining the outstanding advantages of high computational efficiency and simple programming.
The interpolation coefficients are set as α = 1,2,3,4 and the window sub-matrix is set to the sub-volume size of 41 × 41 × 41 in order to verify the accuracy and stability of the algorithm under different interpolation coefficients α. Calculation results are shown in figure 5 and it can be seen that the IC-GN algorithm has better accuracy and stability. The mean bias error and standard deviation (SD) error of the IC-GN algorithm gradually decrease when the interpolation coefficient α increases from 1 to 3.
However, when α = 4, although the mean bias error continues to decrease, the SD does not continue to decrease significantly, and the standard deviation is close to that when α = 3. The reason for this is that the IC-GN algorithm has a remarkable ability to predict the extreme values. When α increases to a certain value, that is, when the accuracy of initial value reaches a certain value, the initial value is not the main factor that affects the accuracy of the IC-GN algorithm. In addition to memory consumption, the accuracy of the algorithm can be improved by appropriately increasing the interpolation coefficient α.
The accuracy and stability of the algorithm are improved when the interpolation coefficient α increases from one to two. However, the curve-fitting algorithm is less successful than the IC-GN algorithm when the interpolation coefficient α > 2; meanwhile, the mean bias error of the algorithm becomes very unstable. Moreover, the SD error increases by nearly seven times its initial value.
In these cases, the curve-fitting algorithm cannot meet the calculation accuracy requirements. It can be illustrated from the analysis above that: (1) some interpolation errors could be easily caused if the interpolation coefficient is too large; (2) if α is too large, the size of the cubic lattice used by curve-fitting is too small (the dimension of the lattice is 1/α), which means the nodes of a lattice are very close to the extreme position, which may cause numerical errors in the ZNCC calculation. As a result, the curve-fitting method becomes more sensitive to interpolation errors. These two aspects may affect the accuracy and system stability of the curve-fitting algorithm.
Since the accuracy and stability of the algorithm could not be improved by the increase in the interpolation coefficient, it is necessary to find the optimal α. The memory requirement is increased by 27 times its initial value when the interpolation coefficient α = 3. The memory requirement is sharply increased to 64 times its initial value when the interpolation coefficient α = 4, seriously reducing the applicability of the algorithm. Considering the calculation accuracy and memory consumption of the algorithm, the interpolation coefficient α =2 is the best choice. This not only ensures calculation accuracy but also does not occupy a massive amount of memory.
Moreover, as shown in figure 6, the curve-fitting algorithm has better accuracy than the IC-GN algorithm when α = 2. The SD of those two algorithms are similar; however, the mean bias error differs. The mean bias error of curve-fitting is 15 times less than that of IC-GN. Although curve-fitting can only deal with small strain deformation of materials, the improved curve-fitting not only has good accuracy but also has very high calculation efficiency. Therefore, it has a very high utilization value in small strain situations.

Computational efficiency of the proposed DVC method
The proposed algorithm has better accuracy and stability when the interpolation coefficient α = 2. In this section, the computational efficiencies of the two algorithms under different conditions are compared. Selecting α = 1 and 2, the computational speeds of the curve-fitting method and the IC-GN method are compared. The IC-GN algorithm has the same computational speed as existing incremental DVC algorithms when α = 1. Due to the use of different programming  languages and computer configurations, the computational speed in this paper is slower than that in [8] (MATLAB programs are much slower than C + + in general). As shown in figure 7, the computational speed of the curve-fitting method can reach 757.38 points s -1 and the improved IC-GN algorithm can reach 237.35 points s -1 when α = 2 and sub-volume size is 41 × 41 × 41. The computational speed of the curve-fitting algorithm is 3.2 times that of the IC-GN algorithm. It can be seen that the curve-fitting algorithm is much faster than the IC-GN algorithm. The computational speed of the small deformation region in the material can be greatly improved and the displacement of the large deformation region can also be accurately estimated when calculating the strain field of a real material.

Discussion and conclusions
In this work, a pre-interpolation method is proposed that typically uses a 1/2 sub-voxel interpolation. It was only necessary to calculate the pre-interpolation once and the reference sub-volumes in the updated reference volume image automatically translated into the nearest half-integer voxel position. Therefore, the redundant sub-voxel interpolation calculation was avoided completely. Compared with the integer increment algorithm, the proposed method can improve the accuracy and retain high computational efficiency. Additionally, the curve-fitting algorithm was improved by the preinterpolation method, whose accuracy and computational efficiency improved significantly. It was also found that if the interpolation coefficient larger than two, the error of the IC-GN algorithm and the curve-fitting algorithm increases; simultaneously, the memory requirement also increases significantly.
The simulation results show that the proposed improved DVC algorithm has better performance in accuracy, computational efficiency and ease of implementation compared to the existing DVC algorithm. In the MATLAB programming environment, the computational speed is increased by 1.0-3.2 times the original value, which can effectively improve the accuracy and efficiency of the DVC algorithm. Future studies will investigate how to use IC-GN and curvefitting algorithms more effectively. This is the curve-fitting method.
In this appendix, the steps of solving the curve-fitting equation (6) are illustrated in detail. First, for the convenience of writing, the 19 nodes in the cube lattice are marked as follows: Solving equation (A6) and the sub-voxel displacement is given as: