AN ADVANCED COARSE-FINE SEARCH APPROACH FOR DIGITAL IMAGE CORRELATION APPLICATIONS FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 14, N o 1, 2016, pp. 63 - 73 Original scientific paper AN ADVANCED COARSE-FINE SEARCH APPROACH FOR DIGITAL IMAGE CORRELATION APPLICATIONS  UDC 004.93 Samo Simončič1, Melita Kompolšek2, Primož Podržaj1 1Faculty of Mechanical Engineering, University of Ljubljana, Slovenia 2Upper-Secondary School of Electrical and Computer Engineering and Technical Gymnasium, Ljubljana, Slovenia Abstract. The paper presents a newly developed fine search algorithm used in the application of digital correlation. In order to evaluate its performance a special purpose application was developed using C# programming language. The algorithm was then tested on a pre-prepared set of the computer generated speckled images. It turned out to be much faster than the conventional fine search algorithm. Consequently, it is a major step forward in a never ending quest for a fast digital correlation execution with sub pixel accuracy. Key Words: Digital Image Correlation, Motion Detection, Deformation, Sub-pixel Registration, Coarse-fine search Scheme 1. INTRODUCTION In the last two decades, the research field of digital image correlation (DIC) [1, 2] has been growing at an ever increasing pace. This can be attributed to its distinct advantages such as a simple experimental set-up, low requirements regarding experimental setup and a wide range of applicability. In addition, the digital image correlation has been widely used for deformation and shape measurement, mechanical parameter characterization and numerical-experimental cross validations [3-5]. The basic idea behind the standard and most widely used subset-based DIC method [6] is to track the subsets (or sub-images) specified in a reference image through the sequence of deformed images. The result of this procedure gives us the information about the full-field motion and deformation. It is evident that the proposed approach is rather Received August 20, 2015 / Accepted November 8, 2015 Corresponding author: Samo Simončič Faculty of Mechanical Engineering, University of Ljubljana, Aškerčeva cesta 6, 1000 Ljubljana, Slovenia E-mail: samo.simoncic@fs.uni-lj.si mailto:samo.simoncic@fs.uni-lj.si 64 S. SIMONČIČ, M. KOMPOLŠEK, P. PODRŢAJ simple and intuitive, but sub-pixel registration accuracy and computational efficiency are two important aspects which need to be considered if the approach is to be applied. In the first case (accuracy of DIC measurements) the obtained accuracy depends on various factors, such as speckle pattern [7], subset size [8], correlation criterion [9], shape function [10], sub- pixel interpolation scheme [11] and sub-pixel registration algorithm [12] where the latter one has the major influence on the registration accuracy of DIC. For this purpose different types of sub-pixel registration methods have been developed. Currently, an iterative spatial domain cross-correlation approach (for example a Newton-Raphson [13] (NR) approach) is one of the most widely used sub-pixel registration algorithms. In the last decade, original NR approach [14] has been improved significantly [13, 15] by reducing its computation complexity [16], improving its robustness [17] and extending its applicability. In this field, it is considered as a golden standard for accurate and robust sub-pixel motion detection in digital image. Its ability to take into the account the relative deformation and rotation of the target subset is one of the main reasons for its wide applicability in various fields of science. It is also capable of providing the best sub-pixel registration accuracy compared to the other methods. The NR approach, however, needs the correlation function, which by nature is nonlinear with respect to the desired mapping parameters vector. This of course implies nonlinear numerical optimization. In such a case the initial value of mapping parameters vector is required, which presents initial guess in the correlation procedure. It should be noted that the initial guess must be defined as accurately as possible because only in this way the convergence of the NR approach is guaranteed [15]. As already mentioned, the NR approach is very appropriate in the cases where the deformation and/or rotation of the target subset are present. On the other hand, if only rigid body translation of the target subset is considered, then it is possible to use the sub-pixel registration approaches which are not so demanding from the mathematical point of view. In many situations, this fact provides the possibility to use the methods which are very simple from the theoretical point of view and very straightforward for implementation which also means that they are computationally efficient. In our opinion, the coarse-fine search algorithm [18] is well suited for such a task. Because of that we have decided to implement it together with some improvements which will be presented in detail. The basic idea behind the coarse-fine search strategy can be described by the following steps. In the first step, it calculates the predefined correlation coefficient for all points of interest in the searching area with a 1 pixel step. In order to improve its accuracy it is logical to reduce the searching step, for example to 0.1 pixel or 0.01 pixel. From practical point of view, the value of the searching step depends on accuracy, which is needed in actual application. In many cases, the coarse-fine searching strategy is able to handle only rigid body motion and consequently the shape changes of the deformed subset cannot be evaluated as in the case of the NR approach. If the searching step is less than 1 pixel, the gray level at sub-pixel locations must be reconstructed and for this purpose a certain interpolation scheme is needed. From the execution time point of view this is the most demanding part of the coarse-fine searching approach. Consequently, it should be treated with some care (see reference [19] for example). The conventional fine search methods usually search for the best matching point by a given searching step in the x and y direction, respectively. Thus the computation complexity of such searching method is proportional to the given searching step (x_step × y_step) where x_step and y_step present number of steps in the x and y direction, respectively. In this An Advanced Coarse-Fine Search Approach for Digital Image Correlation Applications 65 manner, the fine search calculation is executed (x_step -1 + 1) x (y_step -1 + 1) times for each sample point. To be more precise, if the searching step is 0.01 pixel in both directions, the algorithm has to be executed 101 x 101 (10,201) times for each sample pixel. It is obvious that the presented scheme is very time consuming since for each sample pixel a great number of the correlation coefficients must be calculated. In similar manner, the interpolation coefficients must be calculated at each execution step as well. Based on this fact, it is evident that such searching method has a high computational complexity. Some work has been done to reduce the computation cost without decreasing the accuracy. The scheme presented in [20] needs to be executed n x 11 x 11 (121n) times for a searching step of 10 -n pixel in both directions. For instance, if the searching step is 0.01 pixel in both directions the fine search approach only needs to be calculated 242 times for the actual sample pixel. This improvement was our inspiration to develop a novel course-fine searching strategy in which the computation cost is significantly reduced compared to the other existing methods in the literature. The proposed approach needs to be executed n x 2 x 2 (4n) times for each sample point with a searching step of 2 -n in both directions. If a searching step is assumed to be 2 -7 (= 0.007812), then the proposed scheme needs to be executed only 28 times for each pixel. From this point of view, it is evident that the proposed approach outperforms the conventional one by a great margin (265 times) in the case when the searching steps are 0.01 and 0.0078 pixels for the conventional and the proposed scheme, respectively. The performance of a novel coarse-fine searching approach is tested and evaluated by a sequence of computer-generated speckled images, where each of them was translated compared to the previous one for some known value. To evaluate the measurement accuracy and effectiveness of the proposed method, we developed a Windows application which was developed specially for this purpose in the Visual Studio development environment using C# programming language. The results obtained by the conventional and the novel approach are also compared and evaluated from the computational complexity point of view. 2. PRINCIPLES OF DIGITAL IMAGE CORRELATION After digital images of the object surface before and after deformation are obtained, the DIC method can be used to calculate the movement of each image pixel. If full-field deformation is required, a ROI in the reference image must be specified within which image pixels are evenly spaced by virtual grids. The basic idea of the DIC method is to match a reference mask (subset) in the original image with a deformed mask (subset) in the image after deformation (deformed image) as illustrated in Fig. 1. We have assumed that the reference subset is a square with a reference pixel in its center. Once the location of the target subset in the deformed image is found, the displacement components of the reference and target subset centers can be determined. 66 S. SIMONČIČ, M. KOMPOLŠEK, P. PODRŢAJ Fig. 1 Schematic representation of the subset before and after deformation In order to estimate the degree of similarity between the reference and the deformed subsets, a certain correlation criterion should be defined first. As already mentioned, the evaluation of the proposed approach is tested on the set of computer-generated speckled images where the illumination is controlled. Due to this fact, we have only used a normalized sum of squared differences (NSSD) correlation criterion and not a more robust zero- normalized sum of squared differences (ZNSSD) one. The main drawback of the ZNSSD, however, lies in the fact that it is computationally expensive comparing to the applied NSSD correlation criterion. It should be noted that the used NSSD correlation criterion is insensitive to the linearly varying illumination intensity but sensitive to the offset in the illumination intensity. The latter one is not true for ZNSSD correlation criterion. However, the NSSD correlation criterion is defined by the following equation [1], 2 ),(),(            M Mi M Mj iiii NSSD g yxg f yxf C (1) where f and g are mean intensities of reference and deformed subsets, respectively. It is well known that the digital image is represented by a finite number of pixels. Consequently, one might think that accuracy is limited to one pixel, but, as already mentioned, it is possible to find registration methods with accuracy better than one pixel. Even in this case, however, the first step is the application of an algorithm with one pixel accuracy. This (one pixel accuracy) can be made by a simple searching scheme within the deformed image or with a more advanced approach [21]. In Ref. [21], the authors developed a fast recursive scheme to mathematically reduce the computational complexity of the traditional DIC technique with one pixel accuracy. In general, the sub-pixel methods are only used afterwards to get a more accurate displacement evaluation and this step was done by a novel searching scheme which is presented in more detail hereafter. The proposed approach can be used directly (without simple searching scheme or coarse scheme) if the actual displacement at measured image point between two consecutive images is not larger than one pixel. An Advanced Coarse-Fine Search Approach for Digital Image Correlation Applications 67 3. A NOVEL FINE SEARCHING STRATEGY The starting point for a fine searching method is a square subset of 1 × 1 pixel, where its location in the deformed image is determined by a coarse search. This step is the same as in the other fine searching methods. The obtained square subset is used for searching the best matching point by a given searching step in x and y directions, which is presented by x_step and y_step, respectively. Due to the fact that the correlation coefficient for each sub-pixel location in the square subset must be calculated, the computation of these fine search methods is quite time-consuming. More precisely, if the searching step is assumed to be 0.01 pixel in both directions, then the fine searching scheme will be executed 101 x 101 times for each sample point. This fact was the main reason to develop a novel fine search method in which computational complexity would be significantly reduced. For simplicity, we suppose that the searching steps are of the form x_step = y_step = 2 -n (n = 1, 2,…), where parameter n presents its accuracy. If parameter n is for example equal to 5, the obtained motion error would not be greater than 2 -5 (= 0.03125) pixel. Based on the known searching step, a novel searching scheme can be defined as follows. In the first step, searching for the matching point at searching steps )5.0(2 1   is performed. The square subset of 1 × 1 pixels centered at the location given by the coarse search is divided in four subsets centered at four sub-pixel locations. For each of them the NSSD correlation coefficient is calculated and the location of the point with the best match is denoted as ),( 11 yx . In the second step, the searching step is reduced from 1 2  to 2 2  pixel and the searching procedure is performed in the same manner as in the first step, but in this case a square subset of 0.5 × 0.5 pixels centered at ),( 11 yx is used for searching area. Thus, the location of the best match when searching with 2 2  pixel step is denoted as ),( 22 yx . This procedure is repeated until the value of the searching step is sufficiently small compared to some predefined threshold. If the searching step is assumed to be n 2 pixels for both directions, then one of the four sub-pixel locations in a square subset of 11 22   nn pixels centered at ),( 11  nn yx presents the best match which is denoted as ),( nn yx . To be clear, the points from ),( 11 yx to ),( nn yx present sub-pixel locations where their gray levels are calculated by predefined sub-pixel interpolation algorithm. The idea behind a novel fine searching scheme is clearly presented in Fig. 2. It is evident that the proposed scheme needs to be executed 22n times, if the searching steps are assumed to be n 2 pixels in both directions. Based on the presented theoretical facts, the numbers of iterations which need to be calculated for each sample point at different searching steps are presented in Fig. 3. They are denoted by red and blue bars for the proposed and the conventional approach, respectively. It can clearly be seen that the obtained number of iterations is increasing much slowly than in the case of the conventional approach. For example, if the searching step is 0.0156 pixel in both directions, then the conventional approach needs to be executed 4225 times for each sample point. This means 4201 iterations more than in the case of the novel fine searching scheme. This fact confirms that the proposed approach drastically reduces computational cost which results in a much wider range of applicability. 68 S. SIMONČIČ, M. KOMPOLŠEK, P. PODRŢAJ Fig. 2 Schematic diagram of a novel fine search scheme Fig. 3 Comparison of the number of iterations needed for calculation of the best match at different searching steps between the conventional and the proposed fine searching scheme for each sample point As shown, the number of iterations is obviously reduced for each sample pixel, but the execution time for each iteration has not been modified in the sense of an improved performance at that point. As mentioned, at each iteration the interpolation coefficients must be calculated. This part of the algorithm has a major impact on execution time. Because of that, the calculation of sub-pixel interpolation coefficients must be done carefully. The gray level at sub-pixel locations must be reconstructed and for this propose some kind of sub-pixel interpolation algorithm must be used. The sub-pixel interpolation calculation of a pixel point of a certain reference subset is not only performed in each iteration, but it also needs to be carried out for the same pixel point that appeared in adjacent reference subsets. The repeated interpolation calculation performed at sub-pixel locations consumes a lot of execution time. To eliminate such unnecessary computation, the authors in [19] present an An Advanced Coarse-Fine Search Approach for Digital Image Correlation Applications 69 elegant solution, which provides an efficient calculation of the interpolation coefficients. A more detailed explanation of this idea and its experimental verification are given in [19]. In our experiments, the bicubic interpolation is used to estimate gray level and first- order gray gradient at sub-pixel locations. They are defined by the following equations:     3 0 3 0 ),( m n nm mn yxyxg       3 1 3 0 1 ),( m n nm mnx ymxyxg  (2)      3 0 3 1 1 ),( m n nm mny ynxyxg  where ),( yxg  is gray level at sub-pixel location ),( yx  , ),( yxg x  and ),( yxg y  are first-order gray gradient at the same sub-pixel location ),( yx  in x and y direction, respectively. In bicubic interpolation, the unknown 16 coefficients, i.e., 330100 ,...,, aaa can be calculated by gray intensity of the neighboring 44 pixels centered at the sub-pixel location. It should be noted, that the 16 interpolation coefficients have the same values for each interpolation block, irrespective of the sub-pixel locations within it. In our application, the idea described in [19] was, however, used for efficient calculation of the interpolation coefficients. The basic idea behind it is presented in the following section. 4. PRE-COMPUTED INTERPOLATION COEFFICIENT FOR EFFICIENT BICUBIC INTERPOLATION As presented in the previous section, the proposed fine searching approach needs to be implemented by sub-pixel interpolation method which calculates intensity ),( yxg  at each sub-pixel location when that is necessary. This procedure is the same for conventional approach as well and it is extremely time-consuming, because it takes most of the computational time at each iteration step of the fine searching approach. More precisely, if the searching step is assumed to be n 2 pixels for both directions, then the proposed approach will be executed n4 times at each sample point. In this case, each pixel within the considered subset will be interpolated repeatedly n4 times. Therefore, it is evident that this procedure has a big overlap with the adjacent subsets because the pixel point within one reference subset may also appear in its neighboring reference subset as shown in Fig.4. This means that repetitive interpolations should be performed on the same pixel point. Considering a reference subset of )12()12(  MM pixels dimension and a grid step of L pixels, each pixel in this subset is then also used in the adjacent 2 [floor((2 1) / ) 1]M L   reference subsets. If the considered pixel point is displaced to a sub-pixel location, this pixel point will be interpolated approximately 2 4 [floor((2 1) / ) 1]n M L     times using the conventional method (parameter n presents the desired accuracy of the proposed approach). The study of other 70 S. SIMONČIČ, M. KOMPOLŠEK, P. PODRŢAJ references has shown that the 16 interpolation coefficients which need to be determined, if Equation 2 is to be applied, are repeatedly calculated. As already noted in the case of bicubic interpolation, the 16 coefficients are not modified in each interpolation block regardless of the sub-pixel location. Due to this fact it is possible to form a global look-up table for each interpolation block in the deformed image before a novel fine search approach is executed. Consequently it is possible to determine all the 16 coefficients in advance. As soon as a certain pixel falls in a certain block these coefficients are used to directly reconstruct the intensity at sub-pixel locations. As a consequence, the construction of interpolation coefficients must be performed only once and not 2 4 [floor((2 1) / ) 1]n M L     times for each interpolation block. For example, if we take a subset of 21×21 pixels and the grid step is set to 5 pixels, the number of interpolations is reduced from 75 to just 1. Obviously, a larger subset size and or a smaller grid step will further reduce the computation time of the newly presented approach. Fig. 4 Schematic representation of the redundant interpolation in the newly presented fine search approach. The displaced pixel point (denoted with a blue dot) is repeatedly appearing in the adjacent subsets (denoted with dashed lines) as the reference subsets defined in reference image are overlapping each other. This means that there are redundant interpolations during the optimization of subsets An Advanced Coarse-Fine Search Approach for Digital Image Correlation Applications 71 5. EXPERIMENTAL VERIFICATION To evaluate the measurement accuracy and effectiveness of the proposed method computer-generated speckled images were generated. Each image was translated compared to the previous one for some known value of the displacement vector. As already mentioned the calculation is performed by the conventional and the proposed fine search approaches together with a normalized sum of squared differences (NSSD) correlation criterion and zero-order displacement mapping function. The intensity of reference I1(x,y) and deformed ),( 2 yxI images were calculated by the following expressions 2 2 0 1 2 1 ( ) ( ) ( , ) exp s k k k k x x y y I x y I d           (3) and 2 2 0 0 0 2 2 1 ( ) ( ) ( , ) exp s k k k k x x u y y v I x y I d             , (4) where s is the total number of speckle granules, R is the size of the speckle granule, (xk, yk) are the positions of each speckle granule with random distribution and I 0 k is the peak intensity of each speckle granule. All computations are executed within an application which has been developed especially for this purpose and is written in C# programming language on a personal computer with an Intel Pentium i7, 2.30GHz processor and 8 GB memory. The graphical user interface of the developed program is presented in Fig. 5. The validation of the proposed and the conventional fine searching schemes is performed using the developed application. Fig. 5 Graphical user interface of the developed program for measurements displacements by the proposed fine searching scheme 72 S. SIMONČIČ, M. KOMPOLŠEK, P. PODRŢAJ Fig. 6 compares the computation time of the conventional and the proposed fine searching scheme at different searching steps, ranging from 0.25 to 0.0156 pixels for a fixed subset of 41 x 41 pixel size. As shown, the computation time of the conventional approach begins to increase rapidly as the searching step decreases. Fig. 6 Comparison of the computation time needed for calculation the best match at different searching steps between the conventional and the proposed fine searching scheme for each sample point For example, if the value of the searching step is assumed to be 0.0156 pixels, the conventional approach needs more than 16 seconds for each sample point to find the best match. On the other hand, the proposed approach only needs few milliseconds at the same searching step. From this fact, it is evident that the computation time is extremely reduced. It should be noted that the results presented in Fig. 6 have similar distribution with respect to the searching step as to the number of iterations for each sample point presented in Fig. 3. This is mainly due to the fact that the computation time is directly related to the number of iterations. 6. CONCLUDING REMARKS In order to test our idea of a newly developed fine searching method an application was made in Visual Studio using C# programming language. Within this application the algorithm was tested on a pre-prepared set of computer generated images, which were translated for some predetermined displacement. As expected from the theoretical analysis of both the approaches the results have clearly shown that the proposed fine searching scheme clearly outperforms the conventional method. As the newly developed method turns out to be roughly 1000 times faster for a typical searching step, it vastly widened the set of applicability of digital correlation for real life problems. An Advanced Coarse-Fine Search Approach for Digital Image Correlation Applications 73 REFERENCES 1. Pan, B., Qian, K., Xie, H., Asundi, A., 2009, Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review, Measurement Science and Technology, 20(6), pp. 1-17. 2. Schreier, H., Orteu, J.-J., Sutton, M. A., 2009, Image Corelation for Shape, Motion and Deformation Measurements, Springer. 3. Razi, H., Birkhold, A. I., Zehn, M., Duda, G. N., Willie, B. M., Checa, S., 2014, A finite element model of in vivo mouse tibial compression loading: influence of boundary condition s, Facta Univesitatis Series Mechanical Engineering, 12(3), pp. 195-207. 4. Simončič, S., Klobčar, D., Podrţaj, P., 2015, Kalman filter based initial guess estimation for digital image correlation, Optics and Lasers in Engineering, 73, pp. 80-88. 5. Badaloni, M., Rossi, M., Chiappini, G., Lava, P., Debruyne, D., 2015, Impact of Experimental Uncertainties on the Identification of Mechanical Material Properties using DIC , Experimental Mechanics, 55(8), pp. 1411-1426. 6. Tong, W., 2013, Formulation of Lucas-Kanade Digital Image Correlation Algorithms for Non-contact Deformation Measurements: A Review, Strain, 49, pp. 313-334. 7. Pan, B., Lu, Z., Huimin, X., 2010, Mean intensity gradient: An effective global parameter for quality assessment of the speckle patterns used in digital image correlation, Optics and Lasers in Engineering, 48(4), pp. 469-477. 8. Pan, B., Huimin, X., Wang, Z., Qian, K., Wang, Z., 2008, Study on subset size selection in digital image correlation for speckle patterns, Optics Express, 16(10), pp. 7037-7048. 9. Pan, B., Xie, H., Wang, Z., 2010, Equivalence of digital image correlation criteria for pattern matching, Applied Optics, 49(28), pp. 5501-5509. 10. Schreier, H. W., Sutton, M. A., 2002, Systematic errors in digital image correlation due to undermatched subset shape functions, Experimental Mechanics, 42(3), pp. 303-310. 11. Schreier, H. W., Braasch, J. R., Sutton, M. A., 2000, Systematic errors in digital image correlation caused by intensity interpolation, Optical Engineering, 39(11), pp. 2915-2921. 12. Bing, P., Hui-min, X., Bo-qin, X., Fu-lonh, D., 2006, Performance of sub-pixel registration algorithms in digital image correlation, Measurement Science and Technology, 17(6). 13. Wang, Z., Wang, S., Wang, Z., 2014, An analysis on computational load of DIC based on Newton- Raphson scheme, Optics and Lasers in Engineering, 52, pp. 61-65. 14. Bruck, H., McNeill, S., Sutton, M., Peters III, W., 1989, Digital image correlation using Newton- Raphson method of partial differential correlation, Experimental Mechanics, 29(3), pp. 261-267. 15. Vendroux, G., Knauss, W. G., 1998, Submicron deformation field measurements: Part 2. Improved digital image correlation, Experimental Mechanics, 38(2), pp. 86-92. 16. Pan, B., Li, K., Tong, W., 2013, Fast, Robust and Accurate Digital Image Correlation Calculation Without Redundant Computations, Experimental Mechanics, 53, pp. 1277–1289. 17. Huang, J., et al., 2013, Digital image correlation with self-adaptive Gaussian windows, Experimental Mechanics, 53(3), pp. 505-512. 18. Peters, W., Ranson, W., 1982, Digital Imaging Techniques in Experimental Stress Analysis, Optical Engineering, 21(3), pp. 427-431. 19. Pan, B., Li, K., 2011, A fast digital image correlation method for determination measurement, Optics and Lasers in Engineering, 49(7), pp. 841-847. 20. Zhang, Z.-F., Kang, Y.-L., Wang, H.-W., Qin, Q.-H., Qiu, Y., Li, X.-Q., 2006, A novel coarse-fine search scheme for digital image correlation method, Measurement, 9(8), pp. 710-718. 21. Jianyong, H., Tao, Z., Xiaochang, P., Lei, Q., Xiaoling, P., Chunyang, X., Jing, F., 2010, A high- efficiency digital image correlation method based on a fast recursive scheme, Meas. Sci. Technol., 21, pp. 1-12.