1.

Introduction

Due to its incomparable advantages, such as low environment requirement and easy data processing, digital image correlation (DIC)^1,2 has been applied in many noncontact measurement areas.^3–6 Now with the development of high-speed cameras and computer performance, improving the computation efficiency of the DIC method has become a major problem. Continuous improvements have been made to the DIC algorithm.

Originally, the Newton–Raphson^7,8 algorithm was widely accepted and used to determine the subpixel displacement due to its high accuracy. This method requires several iterations to optimize coefficients and to get the correlational function maximum. Later, the inverse compositional Gauss–Newton^9,10 algorithm was proposed, omitting the repeated calculation of Hessian matrix and greatly improving the computation efficiency. Also, a reliability-guided displacement scanning strategy and a precomputed global interpolation coefficient look-up table are employed to accelerate the integer-pixel displacement searching and subpixel displacement calculation, respectively.¹¹ The calculation of a reliability-guided digital image correlation (RGDIC) method begins with a seed point and is then guided by ZNCC coefficients of the computed points.⁹ This ensures that the calculation path is always along the most reliable direction, and error propagation is avoided. Then, a seed point-based parallel method was proposed to improve the calculation speed, which can maximize the computing speed using an improved initial guess.¹²

In addition, the parallel computing based on graphics processing units (GPU)^13,14 has been found suitable for image processing, including in the DIC field. The compute unified device architecture (CUDA) as a parallel computing platform, developed by NVIDIA, has been widely used due to its flexibility. Through programming, GPU can be assigned different computational tasks. Points of interest can be divided into many groups, and points in a single group can be calculated at the same time. Considering that the path-dependent tracking strategy wastes the parallel computing power of modern computers with multicore processors, the RGDIC method is improved by using a two-section tracking scheme.¹⁵ These calculated points with correlation coefficients higher than a preset threshold are taken as reliable computed points and are given the same priority to extend the correlation analysis to their neighbors. Thus, DIC calculation can be initially executed in parallel at multiple points by separate independent threads. A parallel DIC (paDIC) method powered by GPU parallel computing is proposed.¹⁶ This method combines the IC-GN algorithm for subpixel registration with a path-independent fast Fourier transform-based cross correlation algorithm for integer-pixel initial guess estimation, achieving a superior computation efficiency over the DIC method, purely running on the CPU. Based on paDIC method, a pipelined system framework unifying five variants of combinations of CPU and GPU was proposed, which can be flexibly applied to various practical applications with different requirements of measurement scales and speeds.¹⁷ A fast method was introduced for estimation of dense two- and three-dimensional displacement fields from image correlation. It is based on a local, or window-based, optical flow algorithm, which is ideally suited for parallel processors. This method adopted a coarse-to-fine strategy on a multiresolution image pyramid helping to propagate a good initialization for the optimization throughout the scales; hence, convergence could be reached even for displacements significantly larger than the subset size.¹⁸ However, the repeated iteration is not avoided and is only simplified through these methods.

A DIC algorithm combining the spatial correlation with temporal continuity is proposed.¹⁹ The time-consuming iteration operation during subpixel displacement calculation is replaced by the fitting method based on the moving least-squares technique. Based on this method, a parallel temporal sequence DIC method is proposed. A fast integer-pixel searching algorithm based on the seed points is used to calculate the integer-pixel displacement. The integer-pixel displacement results of seed points are calculated along the time axis, and these results are taken as the initial displacement guess of other points. Thus, the searching field of other points can be greatly reduced and calculated together by separate independent threads. The fitting method based on the moving least-squares technique²⁰ is applied to determine the subpixel displacement, avoiding complex iterations and improving computation efficiency. Also, the fitting process of every point is fully independent, without any hinderance from the neighboring points and the overlapping subsets, such that the advantage of GPU parallel computing can be exploited.

In this paper, GPU parallel computing is used to speed up a temporal sequence DIC method, which chose seed points to determine the integer-pixel displacement and applied the moving least-squares fitting technique to acquire the subpixel displacement. The result shows a super-fast calculation speed under the equivalent precision.

2.

Principle of the Parallel Temporal Sequence DIC Method

The parallel temporal sequence DIC algorithm is divided into three steps, and it mainly considers a continuous deformation process. The first step is to calculate the integer-pixel displacement of seed points along the time axis. The second step is to use the results of seed points as the initial values to determine the integer-pixel displacement of other points, and the calculation of other points takes the advantage of parallel computing to accelerate. At last, the moving least-squares technique is applied to determine the subpixel displacement, also using parallel computing.

2.1.

Integer-Pixel Displacement Calculation

In order to decrease the computation time, the integer-pixel displacement of seed point at every deformation moment is calculated in advance. Along the time axis, the result of the last moment is also taken as the initial guess of the current moment considering the temporal continuity. As Fig. 1 shows, the point $P$ is chosen as the seed point. In each speckle image, the rectangles in different colors stand by the search field, and the search of every moment is centered on the result of the last moment. Since the deformation is continuous in space, once the integer-pixel displacement of one point $p(x,y)$ is acquired, the searching scale of its neighboring points $p_n(x_n,y_n)$ can be reduced greatly. The searching scales of the neighboring points are depended on the complexity of the deformation field. Considering the large deformation and deformation discontinuity in space, more seed points are needed, so the searching scale of neighboring points can be much smaller.

Fig. 1

Integer-pixel displacement calculation schematic diagram.

It is remarkable that the method is not suitable when there is a jump during the deformation process. Large deformation gradient will cause a large computation error, but this kind of computation error can be effectively decreased or avoided by increasing the search scale or by improving the exposure speed of the camera.

Other calculation points usually are chosen in a checkerboard pattern (if the surface of the specimen is irregular such as a hole, the calculation points can then be chosen according to the real situation). The whole points in different locations and in different deformation moments are calculated together. However, if the grid is too tight, the data conflict will still exist when reading gray data from the same speckle image. The strategy is to move the calculation grid. At first, a suitable vertical spacing and horizontal spacing are chosen to achieve a reasonable grid layout to avoid the data conflict, and then the grid is separately moved vertically and horizontally. Thus, the grid can now become tighter within several iterations and all the points of interest can be calculated.

Once the calculation task was determined, the CUDA programming is used to finish the job. As Fig. 1 shows, a multicore GPU usually consists of multiple steaming multiprocessors (SMs), and each SM includes a certain number of steaming processors (SPs). In the CUDA programming, the calculation task was assigned to grids of blocks, and each block contains a number of threads. In order to ensure maximum utilization of the GPU resource and hiding the latency effectively, the number of threads usually required is more than 64 (in this case around 512), and the number of blocks cannot be too small either,⁹ which is according to the quantity of the speckle image analyzed once (in this case around 100). During the calculation process, a kernel consisting of a series of data and instructions is assigned to blocks and then distributed to threads, which run on every SP independently within each block.

2.2.

Subpixel Displacement Calculation

Figure 2 shows the subpixel displacement calculation process at one point. A weighted moving least-squares method is used to fit the integer-pixel displacement along the time axis to obtain the subpixel displacement. The influence domain is defined by the central point and the influence radius, and the central point is the calculated point. A cubic polynomial is adopted as the basic function. However, the integer-pixel displacements happen to be stair-stepping due to the neighboring results being the same. In order to overcome the stair-stepping effect and improve the fitting accuracy, the integer-pixel displacement is weighted according to the correlation coefficient in advance. Previous studies have generally considered all integer-pixel displacements to have equal precision, hence the least squares algorithm is employed. However, accuracy differs among points in practice. The correlation coefficient is a simple representation of reliability, where larger correlation coefficient indicates that the integer-pixel displacement is a better approximation of the actual displacement. Therefore, a weighted function is generated based on the correlation coefficient, such that a point’s weight is larger for larger correlation coefficient, and vice versa. The weighted fitting principle is described as follows:

Eq. (1)

$$J(a_0,a_1,\dots ,a_n)=\sum _{i=1}^m{k}_i{(f-y)}^2=\sum _{i=1}^m{k}_i{(a_0+a_1{x}_i+\dots +a_n{x}_i^n-y_i)}^2,$$

where $J(a_0,a_1,\dots ,a_n)$ represents the distance function, $y$ represents the integer-pixel displacement results at different moments, and $k$ represents the weight coefficients. The function $f(a_0,a_1,\dots ,a_n)$ is the adopted polynomial, and $a_0,a_1,\dots ,a_n$ represents the unknown parameters. $m$ represents the number of integer-pixel displacement results involved in fitting. In order to get the minimum of $J(a_0,a_1,\dots ,a_n)$, $a_0,a_1,\dots ,a_n$ is substituted into Eq. (2):

Eq. (2)

$$\frac{\partial J}{\partial a_j}=2\sum _{i=1}^m{k}_i{(x_i)}^j(a_0+a_1{x}_i+\dots +a_n{x}_i^n-y_i)=0,$$

Eq. (3)

$$\left[\begin{array}{cccc}\sum _{i=1}^m{k}_i& \sum _{i=1}^m{k}_i{x}_i& \dots & \sum _{i=1}^m{k}_i{x}_i^n\\ \sum _{i=1}^m{k}_i{x}_i& \sum _{i=1}^m{k}_i{x}_i^2& \dots & \sum _{i=1}^m{k}_i{x}_i^{n+1}\\ ⋮& ⋮& \ddots & ⋮\\ \sum _{i=1}^m{k}_i{x}_i^n& \sum _{i=1}^m{k}_i{x}_i^{n+2}& \dots & \sum _{i=1}^m{k}_i{x}_i^{2n}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{c}\sum _{i=1}^m{k}_i{y}_i\\ \sum _{i=1}^m{k}_i{x}_i{y}_i\\ ⋮\\ \sum _{i=1}^m{k}_i{x}_i^n{y}_i\end{array}\right].$$

Fig. 2

Weighted moving least-squares method schematic diagram.

The parameters can be obtained according to Eq. (3), then the subpixel displacement of the central point can be calculated. It is remarkable that this method may not work well at the few first moments because the location in the fitting curve of these moments is at the beginning, which may cause large computation error. For example, when the fitting data number is 33, the central number is 16 ($t=16$), but the first to the fifth moment must be calculated through the same curve ($t=\mathrm{0,1},2,\dots ,15$) (Fig. 3).

Fig. 3

Subpixel displacement calculation schematic diagram.

Moreover, the fitting process of different points is fully independent, without any hinderance from the neighboring points and the overlapping subsets, such that the advantage of GPU parallel computing can be exploited. Similar to the integer-pixel displacement calculation, the computation task was divided into many groups and was distributed to every thread within every block. Also, the number of threads required is more than 64 (in this case around 512), and the number of blocks cannot be too small either (in this case around 100).

3.

Experimental Verification

The parallel fitting DIC method was programmed using the C++ language based on CUDA 8.0 and was run on a desktop computer equipped with Intel i7-7700 CPU and a GeForce GTX 1080Ti graphics card. A simulation duralumin fatigue tensile experiment was used to illustrate the accuracy and efficiency of the method. As Fig. 4 shows, the speckle image was obtained as the reference image from a real experiment. Other speckle images were generated with preset subpixel displacements. The subpixel displacements of the translation experiment including 800 deformation speckle images were set according to Eq. (4) as follows:

Fig. 4

Simulation speckle image.

Note that the errors of integer-pixel displacements are all within 1 pixel. In order to evaluate the calculation accuracy and speed of the proposed algorithm, the ICGN method was used to compare. As shown in Fig. 5, the green line was the real displacement curve of a random calculation point along the deformation time axis, and the blue and the red points stood by the calculation results of the proposed algorithm and the ICGN method at each deformation time, respectively. These points agreed well with the actual sine curve, thereby exhibiting that the calculation of the proposed algorithm was correct and highly accurate.

Fig. 5

Displacement result curve.

The average error⁷ and standard error⁷ at every deformation time were used to illustrate the computation accuracy. Figure 6 shows the average error of the proposed algorithm and the ICGN method at different deformation times, and the result showed that the average error of the proposed algorithm is almost at the same level with the ICGN method.

Fig. 6

Comparison between two methods: (a) average error curve and (b) standard error curve.

Tables 1 and 2 show that the calculation speed of the proposed algorithm increased when the calculation image number or the calculation point number at a single deformation time increased. With the calculation requirement increasing, the computing units of GPU device can be fully unitized and the calculation speed increases. At last, it will reach the limit, and the speed becomes stable. According to Table 2, the maximum calculation speed achieved 463,320 POI/s (points of interest per second) when the calculation speed of the ICGN method was 2700 POI/s under the same precision, which verified that the proposed algorithm could greatly improve the computation efficiency compared with the ICGN method. Also, the result showed that the subpixel fitting method was suitable for parallel computing, and the subpixel displacement calculation made up less than 1% of the entire calculation.

Table 1

Calculation speed under different calculated image numbers.

Table 2

Calculation speed under different calculated image numbers.

In order to verify the accuracy and efficiency of the proposed method in the actual experiment, a real low carbon steel tensile experiment was used. The tensile sample size is shown in Fig. 7(a). The loading rate was $0.02\mathrm{mm}/\mathrm{s}$, and 200 speckle images were obtained from the experiment, one of which is shown in Fig. 7(b). Also, the average error⁷ and standard error⁷ at every deformation time were used to illustrate the computation accuracy. (The real displacement of every point at every deformation time was obtained and calculated from the extensometer.)

Fig. 7

The experimental specimen: (a) specimen size and (b) speckle image.

Figure 8 shows the average error and standard error of the proposed algorithm and the ICGN method at different deformation times. The result showed that the average error and standard error of the proposed algorithm are little more than those of the ICGN method. But the average error was within 0.04 pixels, whereas the standard error was within 0.03 pixels, which shows that the computation accuracy of the proposed algorithm was at the same level with the accuracy of present algorithms (which were usually within 0.005 to 0.1 pixels).¹⁴

Fig. 8

Comparison between two methods: (a) average error curve and (b) standard error curve.

According to Table 3, in the actual tensile experiment, the calculation speed achieved 432,866 POI/s when the calculation speed of the ICGN method was 2052 POI/s under the same precision, which verified that the proposed algorithm could greatly improve the computation efficiency compared with the ICGN method. Also, in the simulation experiment, the subpixel displacement calculation made up less than 1% of the entire calculation. The computational efficiency could be further enhanced if a faster integer-pixel displacement calculation method is discovered. Moreover, it can be found that the calculation speed of the ICGN method clearly reduced in the actual experiment, whereas the speed of the proposed method almost remained unchanged. In the stage of subpixel displacement calculation using the ICGN method, more iterations are always needed in the actual experiment due to the lower quality speckle image compared with the simulation experiment, but its effect on the speed of the proposed method can be ignored because there were no iteration progress during the proposed calculation method.

Table 3

Calculation speed.

4.

Conclusion

In this study, a new parallel temporal sequence DIC method was described. This method chose seed points to determine the integer-pixel displacement and applied the weighted moving least-squares fitting technique to acquire the subpixel displacement. The complex iterations were completely avoided, and the computation efficiency was further improved by parallel computing. The computational accuracy and efficiency of this proposed method were analyzed, and the conclusions are as follows.