Surface reconstruction method for measurement data with outlier detection by using improved RANSAC and correction parameter

Abstract

The moving least squares (MLS) and moving total least squares (MTLS) are two of the most popular methods used for reconstructing measurement data, on account of their good local approximation accuracy. However, their reconstruction accuracy and robustness will be greatly reduced when there are outliers in measurement data. This article proposes an improved MTLS method (IMTLS), which introduces an improved random sample consensus (RANSAC) algorithm and a correction parameter in the support domain, to deal with the outliers and random errors. Based on the nodes within the support domain, firstly the improved RANSAC is used to generate a model for establishing the group of pre-interpolation and calculating the residual of each node. Subsequently, the abnormal degree of the node with the largest residual is evaluated by the correction parameter associated with the node residual and random errors. The node with certain abnormal degree will be eliminated and the remaining nodes are used to obtain the approximation coefficients. The correction parameter can be used for data reconstruction without insufficient or excessive elimination. The results of numerical simulation and measurement experiment show that the reconstruction accuracy and robustness of the IMTLS method is superior to the MLS and MTLS method.

Introduction

After decades of development, the reconstruction algorithms for discrete data have played an essential role in many engineering and scientific fields, especially in error analysis and data processing. The conventional numerical methods, such as finite element method (FEM), interpolate or approximate the nodes through defining a mesh based on known nodes. ¹ However, the accuracy of reconstruction will be greatly reduced and the fitting will even fail when the grid-based method is used to deal with large deformation and discontinuity problems. ² In addition, the human-labor and time cost of generating meshes in complex-shaped domains are not satisfactory. While meshless methods use node-based approximation without mesh discretization, the efficiency and accuracy of data processing are, therefore, greatly improved. ³ In order to meet the development requirement of different fields, a multitude of meshless methods have been presented and applied, such as the element-free galerkin (EFG), ⁴ diffuse elements method (DEM), ⁵ and moving least squares (MLS).

Mclain ⁶ proposed the weighted least squares method. On this basis, Lancaster and Salkauskas ⁷ introduced the moving concept and proposed the MLS method in 1981. To this day, the MLS has become an important method for constructing shape function. Unlike the traditional least square method using complete polynomials, the shape function of the MLS is composed of a coefficient vector and a basis function vector, which can obtain higher continuity under the condition with low order basis function. ⁸ The introduction of the weight function with compact support makes the reconstructed curve or surface accurate and smooth, which has contributed to the wide application of the MLS in various fields. For example, the MLS is applied to solve elasticity problems, ⁹ the compressible Navier-Stokes, ¹⁰ Kuramoto-Sivashinsky, ¹¹ and Burgers equation, ¹² and estimate mathematical model based on discrete points. ¹³ Due to its good performance, the MLS method is often used in combination with other methods to construct shape function. In finite element analysis, the introduction of the MLS method enhances the shape function of the active element. ¹⁴ In the smoothed particle hydrodynamics (SPH), the MLS method was used to construct kernel functions to obtain higher consistency. ¹⁵

The MLS method obtains the local fitting coefficients through the weighted least squares (LS) ¹⁶ method in the support domain ¹⁷ based on the Gauss Markov error model, in which only the dependent variable contains errors. ¹⁸ The total least squares (TLS) ¹⁹ method is an estimation method based on the errors-in-variables (EIV) model. Unlike the Gauss Markov model, the errors of both independent variable and dependent variable are considered in the EIV model. ²⁰ By replacing the LS estimation with the TLS estimation in the support domain, the MLS method is transformed into the moving total least squares (MTLS). ²¹

Nevertheless, due to the impact of factors, such as the disturbances in the measurement environment, the measurement data often contains outliers that seriously deviate from the actual value.^22,23 The LS and TLS respectively used in the support domains of the MLS and MTLS method are not robust estimation methods.^24,25 When there are outliers in the domain, larger deviations will exist in the fitted values around outliers. Many studies have been carried out to reduce the negative impact from the outliers, with the proposed robust algorithms dividing into two main forms. One is to select a subsample from the discrete points to obtain the regression coefficient.^26,27 If the outlier is in the subsample, it will be automatically eliminated. However, this method may eliminate some of non-outliers, in which case the accuracy of approximation will be significantly affected. ²⁸ The other type of algorithm is to identify outliers first and then weaken the influence of outliers by assigning weights to them. In this method, it is difficult to determine appropriate weights when multiple outliers with different levels exist in the discrete data.

In order to weaken the impact of outliers on reconstruction, we propose an improved MTLS (IMTLS) reconstruction method in this article. In the support domain, we deal with outliers by introducing an improved random sample consensus (RANSAC) and a correction parameter, and then the local approximation coefficients are obtained based on the TLS estimation. RANSAC is a robust model estimation algorithm, especially when measurement data contains high proportion of outliers. However, it has limitations in data fitting influencing its accuracy and stability. ²⁹ Therefore, RANSAC algorithm needs to be improved to estimate a relatively reliable initial model. The correction parameter associated with the random error and the node residual is introduced to detect and eliminate abnormal node.

The rest of this article is structured as follows: the second section is a brief introduction to the MLS, MTLS, and RANSAC algorithm, the third section explains the principle and procedure of the proposed algorithm in detail, and the fourth part verifies the performance of the IMTLS method through numerical simulation and experimental data.

Introduction to the basic algorithms

The MLS method

Consider that there are n discrete points x = [x₁, x₂, …, x_n] and corresponding y = [y₁, y₂, …, y_n] in a bounded area Ω. In the MLS method, the trial function f(x) for the fitting point x can be written as

$$f(x)=\sum _{j=1}^m{p}_j(x)a_j(x)={\mathbf{p}}^T(x)\mathbf{a}(x)$$

(1)

where p(x) =  [p₁(x), p₂(x), …, p_m(x)] ^T is a vector of the monomial basis function and the number of the basic is m, and a(x) = [a₁(x), a₂(x), …, a_m(x)] ^T is the coefficient vector to be solved. In this article, for 2D curve and 3D surface fitting, we choose p(x)=[1, x] ^T and p(x)=[1, x, y] ^T , respectively.

In order to solve the optimal a(x), the construct function E(x) is a quadratic function about a(x), that is,

$$E(x)=\sum _{J=1}^n\theta (\Vert x-x_J\Vert /r){(\sum _{j=1}^m{p}_j(x_J)a_j(x)-y_J)}^2$$

(2)

where r is used to control size of the support domain, and θ(||x–x_J||/r) is a non-negative and compactly supported weight function to attribute a weight to each node according to its position relative to x. The fitting property of MLS algorithm will be influenced by the weight function. For example, the fitting form will be interpolated when θ(0) = ∞. ³⁰ There are many kinds of weight functions with compact support such as exponential, Gaussian, and cubic spline weight functions. This article chooses exponential weight function, that is,

$$\theta (s)=\{\begin{array}{c}\frac{\exp (-\alpha s^2)-\exp (-\alpha )}{1-\exp (-\alpha )}\mathrm{if}s\le 1\hfill \\ \multicolumn{1}{c}{0\mathrm{else}}\end{array}$$

(3)

where α is a coefficient related to the convergence speed.

According to the principle of least squares, the coefficient vector a(x) can be solved by

$$\frac{\partial E(x)}{\partial a_j(x)}=2\sum _{i=1}^N\theta \left(s\right)[\sum _{J=1}^m{p}_J(x_i)a_J(x)-y_i]p_j(x_i)=0$$

(4)

where j = 1, 2, …, m, solving equation (4) to obtain

$$\mathbf{a}(x)={\mathbf{A}}^{-1}(x)\mathbf{C}(x)\mathbf{u}$$

(5)

where

$$\begin{array}{c}\mathbf{u}=[y_1,y_2,\cdots ,y_N{]}^T\hfill \\ \mathbf{A}\left(x\right)=\sum _{i=1}^n\theta (\Vert x-x_i\Vert /r)\mathbf{p}(x_i){\mathbf{p}}^T(x_i)\hfill \\ \mathbf{C}\left(x\right)=[\theta (\Vert x-x_1\Vert /r)\mathbf{p}(x_1),\theta (\Vert x-x_2\Vert /r)\mathbf{p}(x_2),\hfill \\ \cdots ,\theta (\Vert x-x_n\Vert /r)\mathbf{p}(x_n)]\hfill \end{array}$$

Taking equation (5) into equation (1), the function f(x) can be presented as

$$f(x)={\mathbf{p}}^T(x){\mathbf{A}}^{-1}(x)\mathbf{C}(x)\mathbf{u}$$

(6)

The MTLS method

Unlike the MLS using LS estimation to obtain the local fitting coefficients, the MTLS method applies the TLS estimation for those coefficients, in which the errors of both dependent variable and independent variables can be considered. In TLS, the solution of linear equation is considered, that is,

$$(\mathbf{B}+{\mathbf{E}}_B)\mathbf{X}=\mathbf{L}+\mathbf{e}$$

(7)

where Band L are coefficient matrix and observation vector respectively, E _B and e are random error matrix of B and L respectively, and X is the coefficient vector to be solved. According to the principle of total least squares, the solution of equation (7) can be transformed into a constrained optimization problem, that is,

$$\mathbf{X}=\arg min({\mathbf{e}}_B^T{\mathbf{e}}_B+{\mathbf{e}}^T\mathbf{e})$$

(8)

where e _B  = vec(E _B ).

Therefore, X can be solved by Lagrange multiplier method. The function is constructed as

$$\mathbf{J}(\mathbf{e},{\mathbf{e}}_B,\mathbf{\phi },\mathbf{X})={\mathbf{e}}^T\mathbf{e}+{\mathbf{e}}_B^T{\mathbf{e}}_B+2{\mathbf{\phi }}^T(\mathbf{L}+\mathbf{e}-\mathbf{BX}-{\mathbf{E}}_B\mathbf{X})$$

(9)

where φ is the vector of Lagrange multiplier. Let the first partial derivative of J with respect to e, e _B , φ, and X equal to 0, that is,

$$\{\begin{array}{c}\frac{1}{2}\frac{\partial \mathbf{J}}{\partial \mathbf{e}}=\mathbf{e}+\mathbf{\phi }=0\hfill \\ \multicolumn{1}{c}{\frac{1}{2}\frac{\partial \mathbf{J}}{\partial {\mathbf{e}}_B}={\mathbf{E}}_B-\mathbf{\phi }{\mathbf{X}}^T=0}\\ \multicolumn{1}{c}{\frac{1}{2}\frac{\partial \mathbf{J}}{\partial \mathbf{\phi }}=\mathbf{L}-\mathbf{BX}+\mathbf{e}+\mathbf{\phi }-{\mathbf{E}}_B\mathbf{X}=0}\\ \multicolumn{1}{c}{\frac{1}{2}\frac{\partial \mathbf{J}}{\partial \mathbf{X}}={\mathbf{B}}^T\mathbf{\phi }+{\mathbf{E}}_B^T\mathbf{\phi }=0}\end{array}$$

(10)

derived from equation (10), we can obtain

$$\{\begin{array}{c}{\mathbf{B}}^T\mathbf{AX}-{\mathbf{B}}^T\mathbf{L}=\mathbf{X}\mu \hfill \\ \multicolumn{1}{c}{{\mathbf{L}}^T\mathbf{L}-{\mathbf{L}}^T\mathbf{BX}=\mu }\end{array}$$

(11)

where μ =  ${\mathbf{e}}_B^T$ e _B  + e ^T e.

The equation (11) can be transformed into

$${\left[\mathbf{B}\mathbf{L}\right]}^T\left[\mathbf{B}\mathbf{L}\right]\left[\begin{array}{c}\mathbf{X}\hfill \\ \multicolumn{1}{c}{-1}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{B}}^T\mathbf{B}\hfill & {\mathbf{B}}^T\mathbf{L}\hfill \\ \multicolumn{1}{c}{{\mathbf{L}}^T\mathbf{B}}& {\mathbf{L}}^T\mathbf{L}\hfill \end{array}\right]\left[\begin{array}{c}\mathbf{X}\hfill \\ \multicolumn{1}{c}{-1}\end{array}\right]=\mu \left[\begin{array}{c}\mathbf{X}\hfill \\ \multicolumn{1}{c}{-1}\end{array}\right]$$

(12)

It can be seen from equation (12) that the problem of solving total least squares is transformed into obtaining the eigenvalues and eigenvectors of matrix [B L] ^T [B L]. When μ is taken as the smallest eigenvalue of the matrix [B L] ^T [B L], X in the corresponding eigenvector is the solution of TLS. For the MTLS method, the construction of the fitting function is consistent with the MLS method, and the same weight function is used for the MTLS. The acquisition of local regression coefficients can be expressed as

$$\mathbf{D}={\left[{\mathbf{B}}_W{\mathbf{L}}_W\right]}^T\left[{\mathbf{B}}_W{\mathbf{L}}_W\right]={\mathbf{M}}^{-1}\mathbf{\Lambda }\mathbf{M}$$

(13)

where

$$\begin{array}{c}{\mathbf{L}}_W=[\theta (\Vert x-x_1\Vert /r)y_1,\theta (\Vert x-x_2\Vert /r)y_2,\hfill \\ \cdots ,\theta (\Vert x-x_n\Vert /r)y_n{]}^T\hfill \\ {\mathbf{B}}_W=[\theta (\Vert x-x_1\Vert /r)\mathbf{p}(x_1),(\Vert x-x_2\Vert /r)\mathbf{p}(x_2),\hfill \\ \cdots ,\theta (\Vert x-x_n\Vert /r)\mathbf{p}(x_n){]}^T\hfill \\ \mathbf{\Lambda }=\mathrm{diag}({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_{m+1}),\hfill \\ \multicolumn{1}{c}{({\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_{m+1})}\\ \mathbf{M}=[{\mathbf{v}}_1,{\mathbf{v}}_2,\cdots ,{\mathbf{v}}_{m+1}],\mathbf{D}{\mathbf{v}}_i={\lambda }_i{\mathbf{v}}_i\hfill \\ {\mathbf{v}}_i=[v_{i,1},v_{i,2},\cdots ,v_{i,m+1}{]}^T(i=1,2,\cdots ,m+1)\hfill \end{array}$$

If λ_m+1≠λ_m, combined with equations (12) and (13), the regression coefficient vector a(x) is:

$$\mathbf{a}(x)=-\frac{1}{v_{m+1,m+1}}{[v_{m+1,1},v_{m+1,2},\cdots ,v_{m+1,m}]}^T$$

(14)

In the MLS method, the coefficient matrix B is error-free, while only the observation vector L is contaminated. Nevertheless, this assumption is usually unrealistic in some fields of engineering, because the errors, such as sampling errors and instrument errors, are always brought into matrix B. Therefore, the MTLS method can give a more reasonable approximation than the MLS method when both matrix B and the vector L contain errors. ³¹

The RANSAC algorithm

The RANSAC algorithm, presented by Fischler and Bolles, ³² can robustly estimate the model parameters. It has been applied to many fields, especially in computer vision,^33,34 due to good performance in handling the data with a tremendous level of outliers and remarkably simple structure. The main idea of the RANSAC is to calculate parameters of hypothesized model by randomly sampling the subsample from the entire dataset, ³⁵ and then the model is performed on the entire dataset. During this process, a distance threshold that is considered empirically ³⁶ is set to calculate the inlier rate of each model, and sampling is stopped when the number of iterations reaches the defined value. The parameters corresponding to the highest inlier rate (the highest consensus) are selected as the model parameters of the total sample. The brief calculation steps of the RANSAC algorithm are described as follows:

Algorithm 1. The RANSAC method.
Input: discrete data set C =  $\{x_i,y_i{\}}_{i=1}^N$ Output: model G(x) Set Q = 0 for (j = 1 to k) do  Randomly select a minimum subset from C to calculate model g(x) and get the residuals   Count the number of inliers q of the current model according to distance threshold  if (q > Q) then   G(x) = g(x) and Q = q end if end for

In order to get a subsample that is all inliers with probability p, the sampling times k can be obtained by the following formula, ³⁷ that is,

$$k=\frac{\log (1-p)}{\log (1-{(1-p)}^{\zeta })}$$

(15)

where w is the proportion of outliers in the entire dataset, and ζ is number of data nodes needed to calculate model parameters. For example, ζ = 2 is taken in the linear parameter estimation. Then a reliable model is obtained through multiple iterations, and the node with the residual outside the distance threshold is considered as outlier.

The improved MTLS method

When there are outliers in the measurement data, the robustness and accuracy of both the MLS and MTLS method are not satisfactory due to their construction principles. ³⁸ Therefore, this article proposes an IMTLS method to reduce the negative impact of outliers on reconstruction. In this method, an improved RANSAC and a correction parameter are introduced into the support domain of the MTLS method to detect and eliminate abnormal nodes, and then the local fitting coefficients based on the TLS estimation can be obtained. The improved RANSAC algorithm is different from the standard RANSAC algorithm in two aspects. Firstly, it considers every possible subsample to generate hypothesized models, which can enhance the stability of data reconstruction. Secondly, the threshold is automatically set to a value related to the random error of the data, which reduces the negative impact of the original method to set a threshold empirically.

Then, we mainly introduce the principle and procedure of the algorithm reconstruction in support domain, as shown in Figure 1. Assuming that there are N nodes in single support domain, the reconstruction process of the MTLS method can mainly be divided into three steps. Firstly, an estimation model is obtained by the improved RANSAC algorithm. Secondly, the abnormal node is eliminated by introducing a correction parameter δ. In this step, whether to eliminate the node with the largest residual is determined by the size of d_max/δ and l, where d_max is the maximum residual and l is automatically set according to the random error of the data. Lastly, local approximation coefficients are determined by TLS estimation. In the support domain, the calculation process of the IMTLS method is shown as follows:

Figure 1.

The fitting procedure in the support domain of IMTLS method.

Algorithm 2. The IMTLS method.

Input: discrete data C = {xi, yi ${\}}_{i=1}^N$ , parameter δ and l Output: local regression parameter a Set Q = 0 Select different combinations including ζ nodes from C,denoted as C1, C2, …, ${{\mathrm{C}}_C}_N^{\zeta }$ for (j = 1 to $C_N^{\zeta }$ ) do   Get model g(x) based on C j and residuals d = {di ${\}}_{i=1}^N$   Count the number of inliers p of the model according to l  if (q > Q) then   Q = q and D = d  end if end for Get maximum residual dmax = max(D) if (dmax/δ ≥ l) then   Eliminate the node with dmax  else   Reserve the node with dmax end if Remaining data is processed based on TLS to get a

After eliminating the abnormal nodes, assume that t nodes $\overline{\mathbf{x}}$  = [ $\overline{x}$ ₁, $\overline{x}$ ₂,…, $\overline{x}$ _t ] and $\overline{\mathbf{y}}$ =[ $\overline{y}$ ₁, $\overline{y}$ ₂, …, $\overline{y}$ _t ] are retained in the support domain. The solution of the TLS can also be obtained by the singular value decomposition (SVD). ³⁹ As shown in the linear equation (7), perform the SVD on the augmented matrix

$$\mathbf{F}=[\mathbf{B}\mathbf{L}]=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T=\underset{m+1t-m-1}{\left[{\mathbf{U}}_1{\mathbf{U}}_2\right]}\left[\begin{array}{c}\mathbf{\Gamma }\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]{\mathbf{V}}^T={\mathbf{U}}_1\mathbf{\Gamma }{\mathbf{V}}^T$$

(16)

in which U = [U₁U₂] _{t  ×  t} and V_{(m+1) × (m+1)} is an orthogonal matrix composed of eigenvectors of matrix FF ^T and matrix F ^T F respectively, and Γ = diag[ξ₁, ξ₂, …, ξ_m+1], ξ₁ ≥ ξ₂ ≥…≥ ξ_m+1 are the singular values of matrix F. Matrix Γ and matrix V can be presented as

$$\mathbf{\Gamma }=\underset{m1}{\left[\begin{array}{cc}{\mathbf{\Gamma }}_1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\mathbf{\Gamma }}_2\hfill \end{array}\right]}\begin{array}{c}m\hfill \\ \multicolumn{1}{c}{1}\end{array}\mathbf{V}=\underset{m1}{\left[\begin{array}{cc}{\mathbf{V}}_{11}\hfill & {\mathbf{V}}_{12}\hfill \\ \multicolumn{1}{c}{{\mathbf{V}}_{21}}& {\mathbf{V}}_{22}\hfill \end{array}\right]}\begin{array}{c}m\hfill \\ \multicolumn{1}{c}{1}\end{array}$$

When V₂₂ is non-singular and ξ_m ≠ ξ_m+1, the solution of TLS is unique, that is,

$$\mathbf{X}=-{\mathbf{V}}_{12}{\mathbf{V}}_{22}^{-1}$$

(17)

For the IMTLS method, the matrix F is weighted, that is,

$${\mathbf{F}}_W=[{\mathbf{B}}_W{\mathbf{L}}_W]={\mathbf{U}}_W{\mathbf{\Sigma }}_W{\mathbf{V}}_W^T$$

(18)

Where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{L}}_W=[\theta (\Vert x-{\overline{x}}_1\Vert /r){\overline{y}}_1,\theta (\Vert x-{\overline{x}}_2\Vert /r){\overline{y}}_2,\dots ,}\\ \multicolumn{1}{c}{\theta (\Vert x-{\overline{x}}_t\Vert /r){\overline{y}}_t{]}^T}\\ \multicolumn{1}{c}{{\mathbf{B}}_W=[\theta (\Vert x-{\overline{x}}_1\Vert /r)\mathbf{p}({\overline{x}}_1),\theta (\Vert x-{\overline{x}}_2\Vert /r)\mathbf{p}({\overline{x}}_2),\dots ,}\\ \multicolumn{1}{c}{\theta (\Vert x-{\overline{x}}_t\Vert /r)\mathbf{p}({\overline{x}}_t{)]}^T}\end{array}$$

The solution expression of equation (17) is rewritten as

$${\mathbf{X}}_W=-{\mathbf{V}}_{W12}{\mathbf{V}}_{W22}^{-1}$$

(19)

In the support domain of the IMTLS method, abnormal nodes are automatically detected and eliminated by the improved RANSAC method and the introduced correction parameter δ. The outliers are evidently included in the abnormal node and can be effectively processed. The local fitting coefficients are obtained by the TLS estimation. As the movement of fitting point x, outliers and random errors are processed in the entire domain.

Case verification

In this section, both numerical and experimental cases were considered by the MLS, MTLS, and IMTLS method, to verify the performance of the proposed method.

Case 1

Take the curve function

$$y=\sin (0.3\pi x)+0.8\sin (0.5\pi x)$$

(20)

as a model to generate uniformly distributed discrete points (x_0i, y_0i). Outliers (0, Δy_i) and random errors of normal distribution with zero mean are added to (x_0i, y_0i) to get simulation data (x_1i, y_1i). The reconstruction points (x_2i, y_2i) are obtained after reconstruction. The sum of absolute differences can be presented as

$$s=\sum _{i=1}^n|y_{0i}-y_{2i}|$$

(21)

which is considered as the index to evaluate three algorithms.

In this case, the parameter

$$l=\sqrt{{\sigma }_x^2+{\sigma }_y^2}$$

(22)

Let n = 201 and r = [max(x₀)−min(x₀)] × 4/100, where min(x₀) = 0 and max(x₀) = 5. Under the same random error condition (σ_x = 0.001 and σ_y = 0.001), we add four different outliers to (x_0i, y_0i), and the s values obtained by three algorithms are shown in Table 1. The fitting curves are shown in Figure 2.

Table 1.

The comparison of three methods in Case 1.

Error parameter	MLS	MTLS	IMTLS
s 1	1.974703	1.995617	0.177981
s 2	3.752185	3.852139	0.177980
s 3	3.854453	3.878887	0.179741
s 4	3.906223	4.180141	0.180511
$\overline{s}$	3.371891	3.476696	0.179053
σs	0.808569	0.864746	0.001107

Figure 2.

The curves fitting by three algorithms under different outlier conditions.

The figures and tables of the above cases show that outliers can be effectively detected and eliminated by the IMTLS method to reduce the incalculable negative impact on the reconstruction. Compared with the MLS and MTLS method, the error values of the reconstructed data processed by the IMTLS method is much smaller.

To further verify the performance of the IMTLS method, we take this case as an example, in which only random errors are added in the simulation data. In Table 2, we can see that the error values obtained by the IMTLS fitting is still minimum among the three reconstruction methods, which shows that the IMTLS method also has good approximation performance under this situation, and Figure 3 shows the bar graph of error values under different random errors.

Table 2.

The comparison of processing random errors in Case 1.

σx	σy	s
		MLS	MTLS	IMTLS
0.000001	0.001	0.274568	0.152626	0.151917
0.00001	0.001	0.274228	0.152063	0.151351
0.0001	0.001	0.275371	0.153263	0.152555
0.001	0.001	0.285897	0.176635	0.176257
0.001	0.0001	0.282678	0.162684	0.161941
0.001	0.00001	0.281492	0.162066	0.161338
0.001	0.000001	0.280397	0.162179	0.161444

Figure 3.

The s values under different random errors.

Furthermore, the CPU time of each algorithm is calculated by MATLAB based on the datasets with different size of n = 201 + i × 200 (i = 0, 1, …, 4). Two conditions are considered and presented in Tables 3 and 4 (r₀ = [max(x₀)−min(x₀)] × 4/100, where max(x₀) = 5 and min(x₀) = 0) respectively. One is to fix the number of nodes in support domain, and the other one is to fix the size of support domain. In the first condition, CPU time of IMTLS method changes steadily with the increase of n. While in the second condition, the CPU time of IMTLS method increases rapidly as n increases because more nodes are included in single support domain. These results illustrate that the number of nodes in the support domain has a significant impact on the CPU time of IMTLS method. Therefore, in order to ensure the efficiency of IMTLS method in data reconstructing, the size of support domain should be appropriately selected to control the number of nodes in support domain.

Table 3.

The CPU times of three methods with a fixed number of nodes in support domain (s).

n	201	401	601	801	1001
r	r0	(1/2)r0	(1/3)r0	(1/4)r0	(1/5)r0
MLS	0.0148	0.0375	0.0477	0.0602	0.0813
MTLS	0.0075	0.0164	0.0250	0.0422	0.0477
IMTLS	0.1473	0.3148	0.4531	0.5977	0.7516

Table 4.

The CPU times of three methods with a fixed radius of support domain (s).

n	201	401	601	801	1001
r	r0	r0	r0	r0	r0
MLS	0.0148	0.0563	0.1227	0.2407	0.4063
MTLS	0.0075	0.0273	0.0547	0.6000	1.1250
IMTLS	0.1473	0.9633	2.9602	45.3609	90.3750

Case 2

In this case, the surface data of standard ball was processed to verify the performance of the IMTLS method. The data is obtained by a commercial white light interferometer (WLI) – Taylor Hobson CCI 3000.

Figure 4 shows the schematic diagram of a WLI system. Firstly, a broadband illumination beam passes through an interferometric objective via a beam splitter. The beams that were reflected by reference mirror and measured surface were focused onto a camera. Interference fringes will generate when the optical path difference (OPD) between the reference and measurement arm is within the coherence length, and the visibility of the fringes increases as OPD decreases. A series of interferograms can be obtained through scanning the objective. Surface data can be obtained by tracking all coherence peaks or phase retrieval within the field of view of the objective.

Figure 4.

Schematic diagram of a WLI system.

The errors contained in experimental data will have a negative impact on obtaining the true profile of the ball measured, and the impact can be reduced by using appropriate algorithm for reconstruction. Three methods are used to reconstruct the experimental data respectively. Then, the reconstructed data are used for parameter regression based on simulated annealing algorithm ⁴⁰ and the regressed radius is used as the evaluation index of reconstruction method. The radius of standard ball is calibrated as 14.402 mm. Since the random error of measurement data cannot be calibrated technically, l is defined to the standard deviation of residuals fitted by the LS estimation in support domain. As shown in Figure 5, data selected in different locations are processed by three algorithms. As shown in Table 5, the radius corresponding to the IMTLS is more proximate to the calibrated value among the three algorithms.

Figure 5.

The processing of experimental data.

Table 5.

The regression radii of three methods in Case 2 (mm).

Radius parameter	MLS	MTLS	IMTLS
R 1	14.340	14.435	14.405
R 2	14.441	14.393	14.408
R 3	14.723	14.632	14.613
R 4	14.615	14.617	14.420
$\overline{R}$	14.530	14.519	14.462
σR	0.149	0.106	0.088

For the correction parameter δ introduced by the IMTLS method, we set different values to observe the change of reconstruction result. As shown in Figure 6(a), when there are outliers in the data, the reconstruction result of the IMTLS method is obviously better than the other two algorithms. And the result will be consistent with that of the MTLS method as δ increases, which shows that when δ increases to a certain value, the node will not be eliminated in the support domain. As shown in Figure 6(b), when there are only random errors, the result of the IMTLS method is still better among the three algorithms. It is clear that δ can be selected in this way for the curve and surface reconstruction.

Figure 6.

Trendlines of three algorithms in Case 1: (a) the data of the fourth group and (b) the data only containing random errors.

Furthermore, we set different number of nodes that can be eliminated in the support domain, to observe the number of points eliminated in each support domain. The data of the fourth group in Case 2 is taken as an example.

As shown in the Figure 7, the l calculated by the residuals fitted by the LS estimation is relatively large in the support domain with outliers, while in the support domain without outliers, the obtained l is relatively small. Therefore, compared with the support domain without outliers, fewer points are eliminated in the support domain with outliers.

Figure 7.

Number of points eliminated in each support domain under different conditions.

Figure 8 demonstrates the processing result with eliminated nodes, some outliers are not effectively eliminated, because the number of eliminated nodes is not sufficient to eliminate all outliers in some support domains. In addition, the processed data are regressed by simulated annealing to get the regression radii, as shown in Table 6. When the outliers are effectively processed, the corresponding regression radius is significantly closer to the calibrated value. This indicates that the outlier has a great negative influence on the parameter regression and can be processed effectively by setting an appropriate number of elimination nodes.

Figure 8.

The fitting surface under different conditions.

Table 6.

The regression radii under different conditions (mm).

R 44	R 43	R 42	R 41
14.420	14.394	14.464	14.484

Through numerical simulation and experimental data verification, the IMTLS method shows better performance. As seen from Case 1, compared with the MLS and MTLS method, the IMTLS method can effectively deal with errors, regardless of whether there are outliers in the data. In addition, the IMTLS method also inherits the good local approximation properties from the MTLS method. We further verify the algorithm by processing measurement data of a standard ball, and the performance of proposed method is evaluated by regression radius. As seen from Case 2, the processed result by the IMTLS algorithm is more proximate to calibrated value.

Conclusion

The MLS and MTLS method show good performance in the fitting of discrete data, such as realizing effective approximation for local geometry feature and obtaining high-order continuous approximation functions with low-order basis functions, etc. However, these two reconstruction methods are not robust, because the outliers in the measurement data will have extremely negative impact on the fitting results. In order to reduce this impact, we proposed an improved MTLS method, in which an improved RANSAC and a correction parameter are introduced into the support domain of the MTLS method to process the abnormal nodes, and then the local fitting coefficients are obtained by the TLS estimation. In this way, the improved MTLS method not only has the advantages of the MTLS method but also can effectively deal with outliers. Practically, we verified the proposed algorithm by dealing with experimental data obtained by CCI. The processing results of three cases show that the performance of the IMTLS method is significantly better than the other two algorithms.