Displacement field determination using an iterative optical flow strategy

This study presents a new strategy of retrieving displacement fields from two different images or two video frames obtained at different times using an iterative optical flow method. First, illumination modeling and image preprocessing are required before the displacement calculations can be performed. Due to the displacement continuity assumption, the displacement differences between the neighboring points are small. Then, the initial displacements of the unprocessed points are determined by the displacements of the processed neighboring points. Next, an iterative optical flow method is performed with a small window size to obtain the displacement offsets between the matching point pairs. The calculation is performed at the subpixel level in each iteration, and thus image interpolation is necessary. In this paper, image upsampling of the deformed image is applied using bi-cubic spline interpolation to achieve the subpixel intensity before all the calculations. This improvement avoids the requirement for image interpolation to be in all iterations. Several numerical tests and real image datasets were used to validate the performance of the proposed method. The results were compared with those obtained by a traditional digital image correlation method using the traditional forward additive Newton–Raphson algorithm and the iterative least square algorithm. Both tests showed that the proposed algorithm exhibits excellent performance regarding computational accuracy and speed.


Introduction
Among the full-field optical techniques available to measure the displacement and deformation of a surface [1,2], 2D digital image correlation (DIC) is considered a well-established method that has been widely accepted and is commonly used in many engineering applications such as mechanical analysis [3][4][5], materialogy [6][7][8] and biomedicine [9,10]. Because it is superior to the 2D DIC method applied for in-plane displacement measurement, the 3D DIC method can address the curve surface deformation and out-of-plane displacement measurement [11,12] and has recently gained increased attention in industrial analysis [13,14]. The DIC method was first proposed by Peters and Ranson in the 1980s [15], and the primary principle underlying this method is the determination of the point in a deformed image that matches the given point in a reference image by cross-correlation of the two image regions surrounding the two points. Due to its simple theory and promising performance for displacement field measurements, the DIC method has attracted significant attention during recent decades. Therefore, it has been modified and improved to enhance its measurement accuracy and computational efficiency. The classical DIC algorithm was proposed by Bruck et al in 1989 [16] and used a forward additive Newton-Raphson (FA-N-R) algorithm to optimize the correlation function. Vendroux et al used an approximate Hessian matrix for the Newton-Raphson (N-R) scheme, greatly simplifying the iteration process and leading to substantial improvement in computational speed without loss of accuracy [17]. A reliability-guided DIC method proposed by Bing et al provided more precise initial predictions for the N-R iteration scheme, which increased the speed of convergence and reduced the iteration times [18]. Instead of the classical FA-N-R iteration process, an efficient inverse compositional (IC) matching strategy combined with the Gauss-Newton algorithm (IC-GN) was implemented for DIC [19] to greatly reduce the computational time without calculating and inversing the Hessian matrix at each iteration. In contrast to the local subset-based DIC method introduced above, in which the computational subsets are unrelated to each other, a finite-element-based (FE-based) global DIC method [20][21][22][23] globally describes the displacement field by discretizing the region of interest (ROI) into finite elements connected with nodes and tracking all elements and nodes simultaneously. The global FE-based strategy solves the problem of displacement discontinuity among the subsets, which cannot be guaranteed by the subset-based DIC method. Bing performed a detailed comparison between the two DIC methods and analyzed their merits and deficiencies regarding efficiency, accuracy, and continuity [24].
Optical flow describes the relative motion between the observer and the moving object in the image plane [25]. The optical flow method has received significant attention from researchers regarding enhancing its performance due to its importance in motion analysis. According to an additional constraint, the optical flow method has two branches: local and global. The most common local constraint is the Lucas-Kanade constraint [26], i.e. the Lucas-Kanade method obtains the optical flow using a least squares method in local regions. The global methods are usually based on the Horn-Schunck constraint [27]; the Horn-Schunck method computes the velocity field by minimizing the brightness constraint and field smoothness terms. The performances of the two different constraints were evaluated in [28].
In this paper, we propose an iterative optical flow method to obtain the displacement field accurately and rapidly. First, local linear intensity modeling between the reference and deformed image is performed. Then, a coarse estimation of the matching points is applied, and the image upsampling of the deformed image is pre-calculated for the subpixel intensity acquisition. After the determination of the matching points in the image pairs, the Lucas-Kanade optical flow method is employed to achieve the displacement offsets and the local intensity change parameters. If the results satisfy the terminal conditions, the calculation of the current point ends. Otherwise, the position of the matching point in the deformed image is updated by the result of the offset calculation, and the offset calculation is looped until the terminal conditions are satisfied. The results were tested on simulated images and standard test images provided by the Society for Experimental Mechanics (SEM) DIC Challenge which was clearly introduced in reference [29]; the results obtained by the proposed method were compared with those calculated by the conventional N-R method in [16] and the iterative least square (ILS) algorithm introduced in Zhou's paper [30]. The comparisons show that the proposed method exhibits excellent performance regarding computing speed and precision.

Description of the proposed algorithm
The proposed method is mainly divided into three steps: image preprocessing and digital upsampling, initial matching points estimation, and displacement offsets calculation. The image preprocessing procedure is necessary to guarantee constant constraints of the pixel intensity and to enhance the signalto-noise ratio (SNR) of the input images. Local linear illumination variation modeling is applied to decrease the negative impact of intensity changes between the images before and after deformation. Digital upsampling of the deformed image is performed to obtain the subpixel intensity, which is needed in the following procedures. In this paper, the upsampling rate is 10, and we have discussed the relationships between the performance of the proposed method and the upsampling rate σ u . Then, a random point in the ROI is chosen as a starting point. The displacement components of the starting point are precisely determined using the time-consuming normalized cross correlation (NCC) method [31]. The initial displacement components of the unprocessed points are determined by the average displacement parameters of the calculated neighboring points. After the initial matching points estimation, an optical flow method is applied to obtain the displacement offsets between the matching point pairs and the illumination modeling parameters using weighted least squares. The optical flow procedure is looped until the differences of the displacement offsets between two iterative results are smaller than 0.01 pixels or the displacement offsets are both less than 1/(2σ u ). A flowchart of the proposed algorithm is shown in figure 1.

Preprocessing of the input images
The pixel intensity constant is one of the hypotheses of the optical flow method due to the brightness constancy constraint equation. Because pictures obtained at different time usually exhibit illumination variations, a linear grayscale intensity model between the reference image and deformed image is needed. The linear model is described as where I 1 and I 2 are the intensity values of the reference image and deformed image, and (x, y) and (x * , y * ) are matching points in the image pairs. Due to the existence of local illumination variation, the parameters k and m cannot be determined before calculation, and they are treated as unknown parameters and calculated in the displacement offset calculation procedure. Smoothing using a low-or bandpass filter is necessary to enhance the SNR of the images. In the proposed method, a Gaussian filter with a standard deviation of 0.5 is applied to smooth the images.

Image upsampling of the deformed image and initial matching points estimation
Initially, in the proposed method, the displacement parameters of the seed point are computed using the NCC method. Due to the displacement continuity assumption, the initial displacement components (u 0 , v 0 ) of the unprocessed points (x, y) are determined by the neighboring points for which displacements have been computed and are shown as follows: where n is total number of the processed neighboring points, and (u i , v i ) is the displacement of the ith computed neighboring points. Since (u 0 , v 0 ) usually are not integer values, the subpixel intensity is needed. In this paper, image upsampling (rate σ u ) is employed to obtain the subpixel intensity values, and the upsampled matrix is denoted as I ud . Bicubic spline interpolation is implemented to obtain the upsampled matrix I ud . Note that digital upsampling is performed before all the calculations, and the subpixel intensity values used in the following calculations are all acquired from I ud .

Subpixel level estimation using the optical flow method
The pixel intensity constant assumption for the linear modeling image pairs is expressed as equation (1). The equation can be expressed as a Taylor series as follows: (3) H.O.T denotes the high-order terms in the Taylor series. Due to the displacement continuity assumption, the differences of the displacements between the initial matching points (x, y) and (x * , y * ) are very small, and the high-order terms of ∆x, ∆y can be neglected. Thus, equation (3) can be written as follows: (4) where I 1x and I 1y are the spatiotemporal image brightness derivatives of the reference image, and ∆x, ∆y are the u and v displacement offsets of the point (x, y). To solve the equation for the unknown parameters, another constraint is needed. The Lucas-Kanade optical flow method assumes a constant displacement in the small region Ω centered by the point (x, y). Then, we can obtain an overdetermined equation set with four parameters: displacement offsets ∆x and ∆y and linear illumination modeling parameters k and m. The equation set is expressed as follows: where I 1 1x , I 1 1y , …, I n 1x and I n 1y are the intensity derivatives of the n points in Ω in the reference image. I 1 1 , …, I n 1 are the intensity values of the n points in Ω in the reference image, and I 1 2 , …, I n 2 are the intensity values of the points in the counter part of Ω in the deformed image. To emphasize the current matching point, a weighted least squares method is applied to solve the equation set [32]. The weighted solution of the equation set is shown as follows: where W is the window function that gives each equation a weighting value, and the point at the center (the current matching point) has the largest value. The 11 × 11-pixel image region Ω centered by the current matching point is applied, and the window function W is separable and isotropic; the 1D weighting values in the horizontal and vertical directions are all described as a Hann window function. The horizontal weighting vector is calculated in figure 2.
The displacement offsets are then obtained after the weighted least squares are determined. If the offsets ∆x and ∆y are both smaller than 1/(2σ u ) or the differences of displacement offsets between two iterations are both smaller than 0.01 pixels, the calculation of the current point ends, and preparation for computation of the next unprocessed point begins. Otherwise, the position of the matching point in the deformed image is updated by the following equations: Then, the position of the matching point in the deformed image is quantized, and the intensity values of the points in the counter part of Ω are acquired from upsampled matrix I ud .
Next, we can obtain ∆x k+1 and ∆y k+1 using equations (5) and (6). The calculation is looped until the terminal condition is satisfied. A brief description of the iterations is shown in figure 3.

Displacement field acquisition
When the iteration of the current computing point is terminated, the displacement (u, v) of the point (x, y) can be expressed as follows:

Verification of the accuracy of the proposed method using numerical images
Computer-simulated speckle images can provide appropriate image features and accurate deformation information. Thus, the displacement parameters of the simulation can be used as benchmark data to examine the performance of the proposed method.

Performance estimation using numerical speckle images
To evaluate the performance of the proposed method, numerical speckle patterns generated by given parameters were utilized in the simulation, and the test results were compared with those obtained by the traditional DIC method and the ILS algorithm. A series of 8-bit grayscale images with only rigid inplane translation were tested. The image pairs were obtained by applying the appropriate shift in the Fourier domain according to the shift theorem [33]. An image pair is shown in figure 4. The displacement components in the vertical direction were set to 0 pixels as a constant. The displacement in the horizontal direction ranged from −1 to 0 pixels with an interval of 0.1 pixel, totaling 11 speckle image pairs. The size of the ROI was 150 × 150 pixels, and the subset for the correlation procedure was a 23 × 23-pixel rectangle for the DIC method. The initial values for the N-R and ILS methods were both obtained by a coarse-to-fine search. In this paper, the upsampling rate was 10, and the interpolation step was 0.1 pixels. This upsampling rate guarantees sufficient precision as well as fast computational speed, which is discussed in the following section.
To quantitatively demonstrate the computational accuracy of the two methods, two types of calculated displacement errors were applied: the mean bias error and standard deviation error. The mean bias error reflects the accuracy of the method, which is defined as follows: where u represents the N measured displacement values, and u ref is the imposed displacement. The standard deviation error reflects the smoothness of the displacement fields obtained, and its mathematical formula is The mean bias errors and standard deviation errors of the u displacement field obtained by the three different algorithms are plotted in figure 5. The proposed method clearly achieved better accuracy than did the conventional N-R method and the ILS method in the numerical image simulations. Since the conventional N-R method is highly precise, the proposed method has a little advantage over the N-R method in calculation accuracy, while the calculation accuracy is highly superior to that of the ILS method. The proposed method avoids redundant image interpolation and Hessian matrix regeneration and is very fast when a suitable terminal condition is imposed. In this experiment, the upsampling rate was set to 10, and the number of iterations was usually 2 or 3 for the offset calculations. Table 1 compares the average elapsed time of the three methods when processing the 11 image pairs. Implementation of the three different methods was accomplished in the MATLAB software environment (MathWorks, USA). The calculations were performed on a desktop computer (four Intel ® Core ™ i5-4590 3.30 GHz CPUs with 8 GB of RAM).
The proposed method achieves a high computing speed that is almost 4 times faster than the conventional method and 2 times faster than the ILS method. The results obtained in the numerical simulations all demonstrate the excellent performance of the proposed method regarding precision and computational efficiency.

Discussion of the impact of upsampling rate σ u on the accuracy and speed of the proposed method using numer ical images
The key factor in determining the accuracy and calculation speed of the proposed method is the upsampling rate σ u , which is the terminal condition of the iterations. When σ u is too large, the calculation accuracy is not guaranteed. However, a too small value of σ u may lead to redundant iterations or even endless loops. Here, we impose different σ u values in the numerical simulation to determine the optimum value to terminate the iterations. The results are shown in figure 6. The average accuracy (AA) is determined by the following equation: where the function abs(x) indicates the absolute value of variable x. According to figure 6, the suitable σ u of the proposed method ranged from 5 to 10, which is consistent with high precision and fast computing. In this experiment, we choose 10 as the upsampling rate to guarantee the efficiency and especially the precision of the proposed method.

Robustness analysis and computational efficiency estimation of the proposed method using SEM DIC Challenge image samples
The DIC datasets provide diverse image samples captured in experiments in which illumination variation and large deformations always occurred. To further estimate the performance of the proposed method, these experimental image samples were used to analyze the robustness of the proposed algorithm.

Verification of the proposed method using image pairs with severe illumination variations
Intensity constancy is the basic premise of the optical flow method; while in most experimental applications,  illumination variations frequently occur between the images captured before and after deformation. The proposed method applies linear modeling and determines multiple intensity change parameters in a small image sub-region, and this method is also suitable for the local illumination changes. In this test, the image size was 487 × 325 pixels, and a rectangular ROI measuring 125 × 200 points with a spacing of 1 pixel was processed using the proposed method. The ground truths of the u and v displacements were −0.3 pixels and −0.6 pixels, respectively. The reference image and the deformed image along with their gray level histograms are shown in figure 7, which demonstrates the serious illumination variations between the images before and after deformations.
The u and v displacement fields obtained by the proposed method are shown in figure 8.
The mean bias errors of the u and v displacements are 0.00069 pixels and −0.0001 pixels, respectively. The standard deviation errors of the u and v displacements are 0.00047 pixels and 0.00037 pixels, respectively. The proposed method achieves high accuracy for applications with variations in illumination.

Verification of the proposed method in rotated images
An application involving a small rotation angle was used to estimate the effectiveness of the proposed method. The size of the image pairs was 512 × 512 pixels, and the ground truth of the rotation angle was 1°. A total of 151 × 151 points in the ROI with an image size of 300 × 300 pixels were processed. The small rotation angle can be expressed by the partial derivatives of displacement components as follows: The images before and after deformation are shown in figure 9; the ROI is also marked in the reference image. The u and v displacement fields were calculated by the proposed method. Since the differential operation amplifies the noise and the calculation errors, it is not reasonable to obtain the rotating angle directly by equation (12). Here, we use pointwise least squares to suppress the noise in the displacement field and obtain precise displacement gradients. The   displacement field in the local region centered on the current point is approximately treated as a linear distribution: Then all points in the local rectangle region satisfy equation (13), and we can obtain an overdetermined equation set: . . .
The rectangle size is (2M + 1) × (2M + 1), and u(x, y) is the u displacement calculated at position (x, y). Thus, the parameters (a 0 , a 1 , a 2 ) can be calculated using least squares for equation solving, and the parameters (b 0 , b 1 , b 2 ) can be obtained in the same way. Then, the rotation angle of each point was obtained using the following equation: The local rotating angle of each point using the displacement filed computed by the proposed method and the pointwise least squares is shown in figure 10.
The maximum value of the rotating angle was 1.005°, and the minimum value was 0.995°, and the calculation error was ±0.5%. The results obtained in this test indicate that the proposed method exhibits high precision for applications with small angles of rotation.

Verification of the proposed method by a tensile beam experiment
To further verify the robustness of the proposed method, a tensile beam test was employed. We calculated 226 × 113 points in the ROI, and the spacing between two neighboring points was 3 pixels. The results obtained by the proposed method were compared with those obtained by the high-precision N-R method. The image pairs and a comparison of the results are shown in figures 10(a)-(c).
The u and v displacement fields of the tensile beam are easily recognized in figure 11(b), as expected. This finding indicates the effectiveness of the proposed method for applications with large deformations. The average differences of the u and v displacement fields between the proposed method and the conventional N-R method in this test were 0.0017 and 0.0021 pixels, respectively. These negligible differences demonstrate the high accuracy of the proposed method. In this application, a large displacement gradient occurred in the y direction. Although the differences of the u and v displacement fields obtained by the two methods in the tensile beam experiment were larger than those of the in the translation and small angle rotating applications, the proposed method achieved high computational accuracy, mainly due to the weighted solution of the linear equation set. This verification demonstrated that the proposed method has excellent performance in terms of the computational accuracy in these applications with large deformations.
The proposed method achieved high accuracy in the three experimental applications with different deformation styles. These results confirmed that the proposed method exhibits good robustness for speckle images.

Computing speed estimation of the proposed method in the real experimental image verification
The accuracy and robustness of the proposed method in three real experimental applications was verified in the previous sections. In this section, we estimated the computational speed of the proposed method and compare it with that of the traditional N-R method and the ILS method in the three verifications. Here, we use computing points per second to describe the calculation speed. The size of the calculation  region for the ILS method was set to 11 × 11 pixels, which is the same as that of the proposed method. The subset sizes for the conventional N-R method in the three applications were 23 × 23, 23 × 23 and 35 × 35 pixels, respectively. The initial values for the ILS method and the conventional N-R method were all obtained by coarse-fine search. The computing speed of the three different methods is shown in figure 12.
The computing speed of the proposed method decreased as the complexity of the displacement field increased. The computing speed advantage of the proposed method over the conventional N-R DIC method and the ILS method was substantial.

Conclusions
In this paper, a novel displacement field determination algorithm using an iterative optical flow method is proposed. The proposed method calculates the displacement offsets by updating the position of the matching points in deformed images and solving a linear equation set. The proposed method calculates the interpolation coefficients before the iterations are performed. Thus, the proposed method avoids redundant image interpolations and Hessian matrix updating, and it is faster than both the conventional N-R method and the ILS method. We also discuss the influence of the key factor on the performance of the proposed method and choose the optimal value to limit the number of iterations. Numerical simulations and experimental estimations were conducted to verify the accuracy, computational efficiency and robustness of the proposed method. All the results showed that the proposed method has excellent performance regarding accuracy, computational speed and robustness.