Meshfree Digital Image Correlation Using Element Free Galerkin Method: Theory, Algorithm and Validation

The association of advanced digital image correlation (DIC) and numerical simulation has been widely used for inverse parameter identification. It is attractive to develop an accurate DIC method sharing the common features with numerical simulation, which can lead to better synergy between experiments and simulations. A new meshfree digital image correlation (MF-DIC) using element free Galerkin method (EFGM) is proposed for deformation measurement. The EFGM is a classical meshfree method in numerical studies, and it is directly used to construct the shape function in MF-DIC from a set of scattered nodes for image matching. The MF-DIC is principally different from the classical local DIC and global DIC since it does not rely on the concept of a subset or an element. In MF-DIC, the C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}^{1}$$\end{document}-continuous displacement for every point is constructed based on a group of scattered nodes in a small support domain surrounding it. The continuous strain map can then be directly derived from the displacement, instead of using an additional smoothing technique as required in classical local DIC or post-processing used in global DIC. A performance assessment based on the Metrological Efficiency Indicator (MEI), as defined in DIC Challenge 2.0, shows that the proposed MF-DIC yields an excellent balance between spatial resolution and measurement resolution for both displacement and strain measurements. Given that the proposed MF-DIC shares common features with the classical meshfree method in computational mechanics, it paves the way for an enhanced synergy between experiments and simulations required for robust inverse parameter identification methods.


Introduction
Digital image correlation (DIC) [1,2] is a powerful and widely accepted optical method for full-field surface deformation measurements. In DIC, displacement/deformation is measured using recorded image series of a sample surface decorated with random features [3]. By comparing the images before and after deformation, using given algorithms, the full-field displacement and strain maps on the sample surface can be retrieved. The DIC technique inherits the attractive advantages of non-contact and full-field measurement, easy implementation and good adaptivity. It has been successfully used to study the mechanical behavior of biological materials [4], composites [5], metals [6], polymers [7], etc. at various challenging conditions, such as high-temperature [8,9], underwater [10], transient deformation [11] environments.
Generally, DIC methods have been developed under two major branches of formulations: local DIC (or subset-based DIC) and global DIC [12,13]. In local DIC, a grid of points of interest (POIs) is specified on the image before deformation (termed as reference image). A subset (commonly rectangular) centered at each POI is matched to the deformed image to retrieve the displacement at each POI. The matching of each POI is independent and fast, and the computational efficiency can be further enhanced using parallelization [14,15]. Local DIC is currently the most popular DIC method because it holds a good balance between spatial resolution, precision and efficiency. However, it still has some limitations. First, displacement continuity is only guaranteed in each subset. The overall strain field is not compatible due to the independent subset matching. Additionally, a large uncertainty arises when a small subset is selected in the local DIC. Moreover, as one of the most important applications of DIC, the displacement field retrieved by DIC can be used to inversely identify the parameters of constitutive material models [16][17][18]. In this regard, it is deemed to be more attractive to use the same kinematic relations in the involved experimental and computational mechanics. There is a natural gap between local DIC and simulation [19] regarding kinematic relations, which can be alleviated by the so-called DIC-leveling approach [20,21]. The key point of the latter approach is the creation of synthetically deformed images based on the nodal displacements provided by the finite element simulation. Although this method enables a fair comparison between DIC and finite element method (FEM), it comes with a large computational effort which hampers inverse identification strategies.
Finite-element-based global DIC [19,22] that integrates the idea of DIC and the FEM is an attempt to alleviate the difference between the displacement description in the DIC measurement and numerical simulation. This method was originally proposed by Sun et al. [23] in 2005. It meshes the region of interest (ROI) into a group of connected elements. All nodes in the mesh are tracked simultaneously instead of independently. The displacement field measured by global-DIC is C 0 -continuous. Additionally, it holds advantages over local DIC for measuring small strain with high strain gradients [22]. However, similar to FEM, the mesh in global DIC arises some typical limitations [24]: • Costly to generate a high-quality mesh for arbitrary geometry; • Hard to guarantee an arbitrary order of continuity for the displacement; In this work, we propose a novel meshfree DIC (MF-DIC) method that directly integrates DIC and meshfree (or meshless) method [24,25]. This method can effectively deal with the aforementioned challenges related to global DIC, while preserving a good balance between spatial resolution and measurement resolution. The meshfree method is an alternative to the popular FEM widely used when simulating the mechanical behavior of solids. It eliminates the constraint of the mesh and approximates the investigated physical fields from a group of arbitrarily distributed nodes. The first meshfree method, the smoothed particle hydrodynamics [26], can be traced back to 1977, close to the year in which DIC emerged [27,28]. Meshfree methods [29] gained a lot of attention since this pioneering work because it effectively eliminates the aforementioned issues. Beyond that, the meshfree method is also a good tool to smooth the displacement fields. Several researchers in the DIC community have employed the meshfree method as a post-processing tool to calculate the strain fields from noisy displacement data or unstructured calculation points [30][31][32]. However, the meshfree method was simply adopted as a displacement smoothing technique instead of DIC matching method. Hence, such a post-processing procedure based on the meshfree method cannot recover the deformation information lost in image matching. In 2016, a meshfree method [33] called "partition of unity" was first integrated into DIC technology to construct the displacement field. This work, however, shows only very brief principal details and primary results.
The MF-DIC proposed in this work is an approach that directly integrates the advanced meshfree idea, i.e. the Element Free Galerkin Method (EFGM) [34], into image matching. We need first to specify a group of scattered nodes, and then the displacement of every pixel on the reference image will be formulated using EFGM from these scattered nodes. For any given point, the displacement is constructed based on the scattered nodes in a small support domain surrounding it. After that, MF-DIC will iteratively retrieve the displacement vector at all nodes simultaneously based on the optical flow conservation [12]. To reduce the redundant computation involved in the repeated evaluation of the Hessian matrix during the optimization, the high-efficiency inverse-compositional Gaussian-Newton (IC-GN) algorithm [35,36] is also employed. Finally, the smoothed displacements and strains at every point can be directly calculated from the retrieved displacement vector of the nodes. The proposed MF-DIC holds some potential advantages compared to conventional local DIC and global DIC. It can measure the displacement map on the whole field with C 1 -continuity [37], enabling direct calculation of continuous strain maps for all points.
The major contributions of this work are: • The EFGM-based MF-DIC that directly integrate EFGM into image matching is proposed. This method is accurate and holds a good balance between spatial resolution and measurement resolution. Additionally, it holds potential in dealing with large deformation and discontinuous physical fields due to the elimination of node connectivity. • The principle and algorithm of the proposed MF-DIC are presented in detail. The performance of MF-DIC is well validated following the method in DIC challenge 2.0 [37]. • IC-GN is used to alleviate the repeated evaluation of the Hessian matrix involved in the forward-additive Newton-Raphson method.
The paper is organized as follows. "General Digital Image Correlation Methods" presents the general principle of local DIC and global DIC. "Principle of Meshfree DIC" first briefly overviews the MF-DIC approach, which is followed by an introduction of the EFGM. The DIC parameter optimization using the IC-GN algorithm and the calculation of smooth displacement and strain fields are described. In the end, a summary of the implementation procedures of the MF-DIC method is outlined. "Experiments and Results" compares the performance of the proposed MF-DIC algorithm with 26 state-of-the-art DIC methods by using the benchmark images provided in the DIC Challenge 2.0 [38]. In "Discussion and Conclusions", the discussion and conclusions about MF-DIC are presented.

General Digital Image Correlation Methods
DIC retrieves the full-field deformation by processing the image series of the object surface. The fundamental assumption in DIC is that the image intensity will be conserved after deformation following the optical flow conservation [12].
The assumption can be formulated by where f and g are the image intensity on the reference image and deformed image, respectively. The displacement at a given point = (x, y) T is ( ) = (u( ), v( )) T . Due to the limited dynamic range and inevitable image noise, the optical flow conservation cannot be perfectly satisfied for all pixels. It is impossible to directly determine the displacement ( ) from a single pixel. Instead, in DIC, the displacement ( ) will be estimated by utilizing the intensity information in a domain containing several pixels. The displacement map will be estimated by minimizing a correlation criterion, such as the sum of squared differences (SSD) over domain Ω (subset in local DIC or ROI in global DIC), where the displacement in the domain can be interpolated with where i ( ) is the shape function related to local coordinate components, u i the unknown related to the displacements, i the index of the displacement component. Both local DIC and global DIC follow this principle. The major differences are the domain size, the form of the shape function and the unknowns.
Local DIC calculates the displacement of each POI independently, see the upper part of Fig. 1. To this end, a subset is selected centered at the POI. The displacement within the subset is assumed to be continuous. It can be formulated from the deformation vector (displacements and their first-order derivatives or even higher-order derivatives) of the POI and the local coordinate in the subset. The displacement at this POI can be determined by minimizing a correlation criterion over this subset. After retrieving the displacement at every POI by independently matching each subset, the strain field can be calculated from the displacement map with some post-processing techniques, such as linear least squares [39,40]. Global DIC is another DIC method that discretizes the reference image into a mesh, see the lower part of Fig. 1. The displacement of a given point can be interpolated from the nodal displacements of the corresponding element covering it. The matching of each element is not solely based on a single element. Instead, it also needs the information from its neighbored elements. Consequently, the global DIC method requires the optimization of the displacements at all nodes over the whole ROI simultaneously.

Basic Principle
In MF-DIC, a group of scattered nodes (see the black dots, see Fig. 2) is first specified on the reference image, just similar to the node selection in conventional global DIC. The MF-DIC does not use the concept of subset or element. Instead, a support domain containing several scattered nodes is used to construct the shape function. As is shown in Fig. 2, for a given point (see the red cross) at location , a surrounding circular or rectangular support domain is selected. In this work, we focus only on the circular support domain. Inside this support domain (the circular region in Fig. 2), there are several nodes, see the purple dots covered in the support domain. The displacement of point can be calculated according to the displacements of the nodes enclosed in its support domain. Every point in the ROI has its own support domain and a group of nodes around it. Then, the cost function in Eq. (2) can be reformulated as

Moving least square approximation of the displacement
The key issue in equation (5) is to construct the shape function ( ) in the support domain. In this work, the EFGM [34] based on the moving least square (MLS) method [41] is used to determine the shape functions. Suppose that the displacement in equation (5) can be expressed as where ( ) = (1, x, y, x 2 , y 2 , xy, ..., ) T is the polynomial basis, and ( ) = a 1 ( ), a 2 ( ), ..., a K ( ) T is the coefficient vector with K the terms of the basis in the polynomial basis ( ) . Commonly, a first-order polynomial basis ( ) = (1, x, y) T or a second-order polynomial basis ( ) = (1, x, y, x 2 , y 2 , xy) T is used, offering satisfactory results for most cases. To determine the coefficient vector ( ) , a weighted residual of the nodes inside the domain is constructed as where W is the so-called weight function. The weight function controls the importance of the nodes inside the support domain of the POI. Here, we adopt the quadric spline weight function and the radius of the support domain, see Fig. 2. To minimize the weighted residual J with respect to ( ) , J is differentiated with respect to ( ): which leads to where Then, the coefficient vector can be determined as Substituting ( ) into equation (6), the MLS approximation of the interpolant displacement vector is estimated by Comparing equations (14) and (5), the shape function matrix is estimated to be with a given ( ) , ( ) can be readily determined from this equation. However, note that ( ) needs to be evaluated for every pixel only during the first iteration. In subsequent iterations the initial evaluation can be reused. Then, we can calculate the displacement at any point using equation (14) if S = (u 1 , u 2 , ..., u M ) T in the corresponding support domain is known.

Strain calculation
The MLS approximation of the displacement field is C 1 -continuous. We can directly calculate the continuous strain where the spatial derivative of the shape function is written as The derivatives T ,x ( ) and W ,x ( − i ) depend on the selection of the polynomial basis and weight function, which can be easily calculated. Then ,x ( ) can be computed from equation (17).
From the former derivation, we can easily calculate the displacement and strain at a given point . We can further simplify the calculation by translating the origin of the coordinate system to the examined point . This can be done by directly replacing with − = 0 and i with i − in the former derivation.

Parameter Optimization Using IC-GN Algorithm
For computational efficiency, the IC-GN algorithm is used here to iteratively minimize equation (4) to determine S . In this method, the displacement vector S in the support domain is updated on the reference image instead of the deformed image. This operation leads to the benefit of only evaluating the Hessian matrix once. Suppose N scattered nodes are selected inside the ROI. Here, we consider the displacement on both x and y dimensions simultaneously. At k-th iteration, the nodal displacement vector (including displacement on both dimensions) at all nodes in ROI is estimated to be k ∈ 2N × 1 , while the displacement vector increment on the reference image is Δ k ∈ 2N × 1 . Substituting the displacement into the cost function in equation (4) leads to (16) where S, k ∈ 2M × 1 and Δ S, k ∈ 2M × 1 are the displacement vector and displacement vector increment of the nodes in the support domain of , respectively. Then, ( ) is updated to Only the parameters in the support domain are included in S, k , and S, k varies with . However, we need to optimize k for all nodes in the ROI simultaneously from equation (21). Similar to the principle of stiffness matrix assembly in finite element method, we reformulate equation (21) to be where matrix ∈ 2M × 2N is the nodal assembly matrix. It is used to transform the displacement vector k ∈ 2N × 1 to the displacement vector S, k ∈ 2M × 1 for the nodes in a support domain. Considering the s-th displacement element in vector S, k , if it is corresponding to the t-th displacement element in k , the s-th row and t-th column element in matrix should be 1. All the elements in matrix are 0 except for the 2M elements with value of 1.
Calculating the derivative of C SSD with respect to Δ k leads to: where ∇ T f ( ) = f x ( ), f y ( ) is the reference image gradient with respect to T . The nodal displacement vector increment can be estimated to with The point + ( ) k is commonly located at the subpixel location, thus the image intensity g( + ( ) k ) on the deformed image needs to be retrieved using image interpolation. Here, Bicubic B-spline interpolation [42] is applied.
The calculated displacement increment is imposed on the reference image, which should be applied to the deformed image instead. The updated displacement vector k+1 on deformed image could be determined either by equation (A6) or improves the accuracy by including an extra correction term but requires more computational resource. If the displacement increments and displacement gradient are small enough, equation (A7) can also lead to fast and accurate convergence, which can be rewritten to equation (28) for global displacement vector updating After a few iterations, if the stop criterion is satisfied, such as the maximum displacement increment at all nodes is smaller than a threshold, i.e., or the maximum iteration number (e.g., 20) is reached, then the optimized k+1 is the final node displacement vector .
Note that on the nodes is not the final displacement vector we need since it cannot be used directly. It is necessary to further calculate u h ( ) and strain for all points in the ROI, even the nodes, using equation (14) and equation (16), respectively. The shape function ( ) does not have the Kronecker delta functions property. The u h ( ) is not necessarily equal to the corresponding one in the iteratively derived .

Summary of Meshfree DIC Method
To facilitate the implementation, the flow chart of the proposed MF-DIC technique is shown in Fig. 3. The procedures can be roughly divided into pre-processing, iterative optimization and post-processing. The pre-processing stage specifies some parameters at first and then calculates ( ) and its derivatives and H. All pre-processing procedures operate only once and the resulting parameters can be reused. The iterative optimization procedures are used to obtain k via the IC-GN algorithm. The displacement of every pixel needs to be computed at each iteration. The image intensity at the corresponding sub-pixel position on the deformed image is retrieved using Bicubic B-spline interpolation [42]. Once the stop criterion is reached, the parameter vector k will be transferred to the post-processing procedure. In this stage, the displacement and strain at all points in the ROI will be calculated.

Experiments and Results
To fairly assess the performance of the proposed MF-DIC, numerical experiments were conducted to evaluate the proposed method by resorting to the methodology used in the DIC Challenge 2.0 [38]. The proposed method is assessed by comparing its measurement results with the results of 26 state-of-the-art algorithms (25 algorithms in DIC Challenge 2.0 plus the inverse-compositional Levenberg-Marquardt (IC-LM) algorithm [43] with first-order shape function).

Test Configuration
In DIC Challenge 2.0 [38], a set of so-called "star" images was numerically generated to evaluate the metrology performance, including spatial resolution, measurement resolution and the Metrological Efficiency Indicator (MEI) [44]. Based on these images, over 20 algorithms were benchmarked with their results available to the public. In this work, we will also make use of the "star" images used in the DIC Challenge 2.0, enabling a fair comparison between the proposed method and the state-of-the-art algorithms. The "star 5" dataset in the DIC Challenge 2.0 is used to evaluate the displacement measurement performance. "Star 5" includes three 501 × 4000 pixels images consisting of a noisy reference image, a noise floor image, and a sinusoidal displacement "star" pattern with linearly varied spatial frequency and a displacement amplitude of 0.5 pixels. The "star 6" dataset is used to assess the performance of the DIC algorithms for strain measurement. Three synthetic 501× 4000 pixels images consisting of one reference image, one noise floor image, and one deformed image with a star strain pattern are also involved in "star 6" dataset. The star strain pattern has the inherent strain field with a linearly varied spatial frequency and a constant amplitude of 0.05. The noise floor images do not have any deformation. They aim to assess the displacement/strain measurement resolution of the methods, while the star images are used to estimate the displacement/ strain spatial resolution of the methods.
The noise floor image and the deformed image were matched to the noisy reference image using a local DIC algorithm (an open-source IC-LM DIC algorithm [43]) and the proposed MF-DIC approach. The local DIC uses an affine shape function. A grid of points with a constant space of 1 pixel was selected on the reference image. These points were matched to the deformed images with subset sizes varying from 9 × 9 pixels to 59 × 59 pixels in increments of 10 pixels. In the MF-DIC approach, a grid of scattered nodes with constant node space increasing from 5 to 33 pixels with an increment of 2 pixels was specified on the reference image. Both first-order and second-order polynomial basis were employed together with quadric spline weight function to approximate the deformation field. The support domain radius follows the node space, which is in this study taken 2.5 times of the node space. To ensure a fair comparison, the uniformly distributed nodes and constant support domain size were specified in this work. However, it must be noted that the nodes can also be randomly distributed, and the support domain size does not necessarily need to be fixed in real applications of the MF-DIC algorithm. Figure 4 shows the vertical displacement maps retrieved by local DIC (upper) and MF-DIC (lower). In local DIC, affine shape function and a subset size of 19 × 19 pixels were used. In MF-DIC, the first-order polynomial basis and a constant node space of 19 pixels were specified. From the left to the right side, the spatial frequency of the vertical displacement decreases linearly. Both methods have better performance on the right side, while failing on the left high-displacementfrequency region. The result from the local DIC method is noisier than MF-DIC.

Comparison with Local DIC
To better visualize the displacement distribution, the vertical displacement curves along the central horizontal line in the measured displacement maps were selected and plotted in Fig. 5. The results indicate that the proposed MF-DIC can measure the displacement accurately and precisely if the displacement period is larger than a given threshold ( ≈ 50 pixels in this case). When the period decreases, the error of MF-DIC will suddenly increase. Local DIC with affine shape function has good accuracy only when the period is large. With the decrease of the period, the measurement of local DIC deteriorates slowly. It must be noted that, the lengthscale parameter is not exactly the same in local DIC and MF-DIC. Consequently, the performance of the two DIC methods cannot be directly compared from a single strain map.

Metrological Performance Evaluation
In this section, the metrological performance of MF-DIC is compared with 26 state-of-the-art algorithms, including 9 local methods with affine shape function, 4 local methods with quadratic shape function, 7 global methods, and 6 miscellaneous approaches [45]. For simplicity, we denote the MF-DIC with first-order polynomial basis and secondorder polynomial basis as "MF-DIC.Order1" and "MF-DIC. Order2", respectively. Different subset sizes or element sizes were used to measure the displacements. For simplicity, we use a length-scale parameter to designate subset size in local DIC, element size in global DIC and the constant node space in MF-DIC (note that for the MF-DIC it is not necessary to have a constant node space). The curves at the central line of the measured displacement and strain maps are publicly available for 25 methods [38], while the remaining one method is the IC-LM method [43] with an affine shape function. We will evaluate the measurement resolution versus length-scale parameters, spatial resolution versus length-scale parameters, and measurement resolution versus spatial resolution.

Displacement measurement evaluation
The measurement resolution of displacement is defined by the standard deviation of the measured displacement of the undeformed noise floor image. Figure 6 shows the displacement measurement resolution versus length-scale parameters for the 26 reference algorithms and MF-DIC. The measurement resolution decreases almost linearly with the lengthscale parameter. The order of the polynomial basis has negligible influence on the measurement resolution of MF-DIC.
The spatial resolution [44,46] is the ability of the method to distinguish high-frequency deformation. To measure it, we use a 12th-order polynomial [38] to fit the central horizontal displacement curves in Fig. 5. The period corresponding to a 10% error of the signal on the fitted curve is defined to be the spatial resolution, see the red square and triangle in Fig. 5. Note that for better fitting the displacement curves in MF-DIC with 12th order polynomial, we removed some data points which are obviously wrongly matched on the left side. The discarding of these points leads to improved detection accuracy for the spatial resolution. Figure 7 presents the central displacement curves and the corresponding fitted curves for MF-DIC with first-order polynomial basis when different node spaces are used. The spatial resolution is measured to be the period at the red square marker. It can be observed that the spatial resolution is well detected in all cases. The spatial resolution versus length-scale parameter curves is presented in Fig. 8. The MF-DIC with first-order polynomial basis has better spatial resolution than MF-DIC with second-order polynomial basis in some cases. In most cases, they have very similar performance.
The former comparisons of spatial resolution and measurement resolution may not be fair because the length-scale parameter in local DIC, global DIC, miscellaneous methods, and MF-DIC cannot be the same. Both spatial resolution and measurement resolution are varied following the lengthscale parameters. However, it is questionable to accurately assess the method performance directly from the spatial resolution or measurement resolution. Therefore, we plot the measurement resolution versus the spatial resolution in Fig. 9. Both spatial resolution and measurement resolution      Fig. 9 Measurement resolution versus spatial resolution for 26 benchmark methods and MF-DIC method are expected to be small. However, there is always a tradeoff between them. In Fig. 9, we found that MF-DIC with first-order polynomial basis and the local DIC with quadradic shape functions are among the best. The MF-DIC with second-order polynomial basis, however, exhibits a poorer balance between spatial resolution and measurement resolution compared to the results obtained with a first-order polynomial basis.
To quantitively compare these methods, we focus on a single indicator independent of the length-scale parameters. To do so, we resort to the MEI, i.e., the multiplication of measurement resolution and spatial resolution. MEI is proven to be almost constant to parameter selection in DIC algorithms [44]. We use the average of the three minimum values of the MEI as the performance indicator of each method following Ref. [38]. The displacement MEIs for all methods are shown in Fig. 10. The proposed MF-DIC shows the best performance among all the evaluated methods. It can be inferred that MF-DIC is comparable with the local DIC methods involving a quadratic shape function. The proposed MF-DIC is obviously better than the global DIC methods, miscellaneous approaches, and the local DIC approaches with an affine shape function.

Strain measurement evaluation
The metrological performance of the proposed MF-DIC for strain measurement is evaluated using "star 6" dataset in the DIC Challenge 2.0. The strain maps measured by MF-DIC.Order1 with node space of 11 pixels and 21 pixels are shown in Fig. 11. Similar to the displacement measurements, the high-frequency strain information is lost on the left, while the low-frequency strain is well measured. With the increase of the node space, the measurement resolution increases while the spatial resolution decreases.
To keep the paper concise, the curves of for spatial resolution versus length-scale parameter, measurement resolution versus length-scale parameter, and measurement resolution versus spatial resolution are not given here for the strain measurement. As shown in Fig. 12, the strain MEIs of these methods were evaluated on the "star 6" dataset in the DIC Challenge 2.0. The strain MEI is defined as the multiplication of the square of the spatial resolution and the the measurement resolution of the algorithms [38]. Similar to the MEIs for displacement measurement obtained with MF-DIC, the MEIs for strain measurement with a first-order polynomial basis are superior to the ones with a second-order polynomial basis. MF-DIC also holds almost the lowest MEI score compared to the state-of-the-art algorithms (it must be noted the MEIs for three algorithms are missing in the publicly available data).

Discussion and Conclusions
A novel meshfree digital image correlation (MF-DIC) method integrating the EFGM-based meshfree idea into image matching is proposed. The IC-GN algorithm is used to speed up the optimization by reducing the redundant calculation of the Hessian matrix in the regular forward additive Newton-Raphson algorithm. The MF-DIC and the meshfree method in computational mechanics use the same node generation strategy and shape function. The latter offers great potential for better bridging the gap between experimental mechanics and computational mechanics, which is important for improving mixed experimental-numerical inverse identification strategies in the framework of Material Testing 2.0 [47].
The MF-DIC is substantially different from the classical local DIC and global DIC, see the implementation procedures summarized in Fig. 13 Fig. 12 Comparison of strain MEI between MF-DIC and 26 state-ofthe-art methods concept of a subset or an element. Instead, it first selects a group of scattered nodes (without mesh connection) on the reference image. For every point on the reference image, a support domain containing several scattered nodes is constructed to determine the displacement for this point using EFGM [34] based on the MLS approximation. The displacement for all nodes is simultaneously determined by minimizing the SSD of pixel intensity for all pixels in the ROI. Simulation tests on the "star" image used in DIC Challenge 2.0 show that the MEIs for both displacement and strain measurements of MF-DIC are among the smallest compared with the state-of-the-art algorithms. It is demonstrated that the proposed MF-DIC holds a good balance between spatial resolution and measurement resolution. The results also show that the method performance is insensitive to the order of the polynomial basis. This might be attributed to the fact that the number (2N) of the unknown displacement components and the available information are the same for both first and second order polynomial basis. Additionally, with the first-order polynomial basis, the weight function offers the method great ability to describe complex deformation.
However, the good accuracy and attractive features of this method are realized at the cost of the more computational burden. The shape function matrix needs to be evaluated for every point in the ROI. This computational burden is partially relieved since it is constant during the iterative optimization and only requires evaluation in the first evaluation. Additionally, the method suffers a lack of accuracy close to the boundary of the ROI, similarly as the conventional local DIC and global DIC. The reason is that the support domain size is reduced at the boundary, leading to an increase of the measurement uncertainty. Figure 13 schematically compares the flowcharts of local DIC, global DIC, and MF-DIC. MF-DIC has some unique features compared with classical local DIC and global DIC. First, it is still a global DIC method but without using a mesh or an element, leading to a series of potential features: • For complex structures, it is much easier to generate a group of scattered nodes than a high-quality mesh; • It holds the potential to deal with large deformation and discontinuity due to the elimination of a mesh compared with global DIC.
Additionally, in the MF-DIC method, continuous strain maps can be directly calculated from the C 1 -continuous displacement map. That means additional smoothing is not required for strain calculation. It is also beneficial in constitutive parameters identification relying on the minimization of a cost function composed of strain components.
As an attractive research area in numerical studies, meshfree methods have been flourishing in the past decades. Several powerful meshfree algorithms, such as the point interpolation method [48], radial point interpolation method [49] and reproducing kernel particle method [50] have been proposed. These methods will continue inspiring DIC research in the future. Furthermore, MF-DIC as a sister method to meshfree methods in computational mechanics, is also promising to be integrated with meshfree methods for constitutive parameters identification.

Appendix A
Suppose that we have iteratively derived the displacement increment Δ k ( ) at point on the reference image based on the displacement vector k ( ) in k-th iteration, where Δ k ( ) and k ( ) is imposed on the reference image and deformed image, respectively. We need to update the displacement vector k+1 ( ) on the deformed image for the next iteration. According to equation (21), we have It is equivalent to The Taylor expansion of Δ k ( ) and k ( ) at position + Δ k ( ) are where ∇ T k ( ) ∈ 2 × 2 denotes the gradient matrix of k ( ) with respect to T . Neglecting the second and higher-order terms, we have (A.1) f ( + Δ k ( )) = g( + k ( )).