Study on the Performances of an Approximating Spline Filter Based on the ADRF Function in Surface Roughness Evaluation

Abstract

Isotropy is an important feature of an area filter in the three-dimensional surface roughness evaluation. First, the transmission characteristic deviation between the approximating spline filter and the Gaussian filter is reduced by cascading approximating. Second, the approximating spline filter is superimposed on the orthogonal direction to obtain an isotropic areal filter. Then, four direct methods for the solving approximating spline matrix are applied. Based on the profile filtering and areal filtering, the computational efficiency and accuracy are compared. The experimental results show that the improved square root method (LDLT decomposition) combines both computational efficiency and filtering precision, and is a good choice for solving the approximating spline matrix. Finally, six kinds of robust approximating spline filters are constructed. Taking the output value of robust Gaussian regression filter (RGRF) as reference, and the honing profile and step surface with deep valley characteristics were used as test surfaces to compare their robustness and iteration time. The experimental results show that the approximating spline filter based on the ADRF function has the shortest iteration times, while its roughness is close to the robust Gaussian regression filter.

1. Introduction

Surface roughness is an important predictor of crack, corrosion, and fatigue damage of mechanical parts, and it represents a tradeoff between the manufacturing cost and performance of mechanical parts [1]. For example, Laura et al. studied the relationship between surface roughness and friction resistance of fluids enhanced with halloysite nanotubes [2]. Therefore, it is very important to accurately evaluate the surface roughness. In precision engineering and production measurement, filters have often been used to suppress unwanted components in measurement data [3]. The surface texture includes three main components: form deviation, waviness, and roughness. Thus, before evaluating the surface roughness, it is necessary to separate the information on roughness from the surface topography signal [4]. The existing filtering technologies can effectively separate different wavelength components. The filter decomposes the measured surface topography signal in the frequency domain and separates different wavelength components into well-defined bandwidths. Besides, there is no need to express the surface morphology as a specific function [5]. At present, the most commonly used filtering technology in the surface roughness evaluation is the Gaussian filter because it has the transmission characteristics of isotropy and phase correction [6]. However, the Gaussian filter produces severe end effect, which results in a large distortion of the filter reference at the boundary [7]. In order to overcome this shortcoming, Krystek constructed a cubic spline filter by combining cubic spline function with the variational method [8]. This spline filter can suppress the end effect by selecting natural boundary conditions. In 2006, ISO/TC213 formally brought the cubic spline filter into the international standard [9]. Although the cubic spline filter has the advantage of suppressing the end effect, there is a great difference in the transmission characteristics between the spline filter and the Gaussian filter. In other words, when these filters are applied to the same surface, their surface metrological evaluation results cannot be compared. Moreover, the areal spline filter formed by superposition has anisotropic amplitude frequency characteristics in the orthogonal direction [10], which limits its wide application in the three-dimensional surface topography. Therefore, many studies on the isotropic areal spline filters have been conducted. For instance, in 2008, Hanada et al. proposed a new frequency-domain areal filtering method, which can obtain transmission characteristics of an arbitrary amplitude [11]. In 2015, Zhang et al. designed a high-order spline profile filter and constructed a spline areal filter with isotropic amplitude characteristics [12]. In 1989, D’Haeyer proposed the approximating spline theory model in order to approximate the transmission characteristics of the Gaussian filters with high accuracy [13]. The difference of transmission characteristics of spline filter and Gaussian filter can be further reduced by cascading approximating. Approximating spline is a further improvement of cubic spline by adding the first derivative to the bending energy of cubic spline. The cascade approximating takes the single step of approximating spline filter as the basic process prototype, and then uses approximating spline filter in series, which is essentially to achieve further approximation of Gaussian transmission characteristics through iteration. In 2009, Zhang et al. further studied the approximating spline theory and realized cascade approximating spline filter [14].

Krystek first applied the numerical solution method to the spline function matrix, which greatly improved the computational efficiency of spline filters [15]. Muralikrishnan et al. applied Gaussian elimination to the cubic spline matrix [16]. Similar to the spline filters, the approximating splines also need to be converted to the matrix form. Therefore, matrix decomposition efficiency and numerical stability are important factors to be considered in the surface roughness evaluation. At present, there are many methods to solve matrix equation, including the direct method [17] and iterative method in numerical analysis. The iterative method generally requires enough number of iterations to achieve the desired accuracy, and it can also sometimes diverge. The direct rule is an effective method to solve large matrix equations (e.g., large banded equations) [18]. The exact solution of spline matrix can be obtained by a finite step operation. Common direct methods include LU decomposition (Doolittle) method [19], LDLT method [20], and Thomas method [21]. LU decomposition method is an improvement of Gauss elimination method, and its essence is to decompose the coefficient matrix into the product of a unit lower triangular matrix and a unit upper triangular matrix. The LDLT method is to decompose the matrix into the unit lower triangular matrix (i.e., all the elements of the main diagonal are 1), the transpose matrix of the lower triangular matrix, and the diagonal matrix (Only the main diagonal of the matrix has elements, and the rest are zeroes). Thomas method is a formula derivation of band preserving decomposition for banded matrix (especially tridiagonal matrix), which is essentially a kind of LU decomposition. In this work, the direct method is applied to approximating spline filter to compare their performance in the profile and areal filtering.

Under the comprehensive action of many factors, including different processing technology and material properties, the surface of parts can produce various shapes and sizes of microscopic geometry, such as scratches, grooves, and other outliers. These outliers do not follow the strict normal distribution strictly, and the interference of outliers can cause the distortion of the filtered mean line. Thus, the authenticity of filtered data is affected. In addition, data distortion caused by surface defect and instrument limitation often appears in the process of surface measurement. These distorted data form outliers on the surface topography, for example, missing data will appear as a blank area on the surface topography. The existence of outliers also affects filtering results. Therefore, robust spline filters have been widely studied. For instance, in 2005, Goto et al. proposed a robust algorithm based on the L2-norm, which has good robustness to outliers [22]. In 2013, Zhang et al. proposed a new robust spline profile filter based on Tukey function [23]. Compared with the robust spline filter specified by international standards, this filter has not only better robustness, but also higher computational efficiency. The robust Gaussian regression filter and robust spline filter have also been included in international standards [24,25]. However, the computational cost of the robust Gaussian regression filter is high. Moreover, the precision of robust spline filter is low and there is no corresponding international standard of three-dimensional. More importantly, the filtering benchmark given by the international standard for robust splines is not satisfactory in all the cases. Therefore, it is necessary to design a robust filter that can not only maintain the precision of the filter but also improve the computational efficiency to replace the standard robust spline filter. At the same time, the difference between the robust Gaussian filter and the robust Gaussian filter should be reduced to obtain a better comparison of measurement results.

In this work, the transmission characteristic deviation between the approximating spline and the Gaussian filter is reduced by cascade approximating, and the transmission characteristic deviation value corresponding to a different approximation order is determined, thus providing reference data for approximation order selection. Also, the cascade approximating method also realizes the isotropic areal spline filter. Besides, direct methods for solving an approximate spline matrix are derived and compared in terms of computational efficiency and accuracy, and relevant suggestions are provided. In addition, six robust approximating spline filters are constructed and used in the experiments, which take honing profile and step surface with deep valley characteristics as experimental objects. The robust Gaussian regression filter is used as a reference filter to compare the robustness and computational efficiency of the constructed filters. The robust approximating spline filter based on ADRF function not only has higher computational efficiency, but also has robustness close to that of the robust Gaussian regression filter. The weight function of ADRF is a descending and cut-tailed weight function. The median absolute deviation which is a more robust scale parameter is used instead of median or standard deviation. The smoothness of Tukey and Andrews weight functions is considered in ADRF weight function. With an increase of deviation, the descending amplitude of weight values is gradually accelerated by stages until it is truncated and abandoned. To sum up, the LDLT decomposition approximating spline filter was realized in this work, which gives consideration to both computational efficiency and accuracy. Based on ADRF robust function, a robust approximating spline filter is designed. This work studies and improves the performance of approximating spline filter, which is helpful to promote the further application of approximating spline filter.

2. Transmission Characteristics

In order to reduce the difference of transmission characteristics between spline filter and Gaussian filter, the international standard ISO 16610-22 also recommends the approximating spline filter with a tension parameter β = 0.62524. As shown in Figure 1, the approximating spline filter in the standard is closer to the transmission characteristics of Gaussian filter compared with the cubic spline filter. The maximum deviation is 4.35%, however, the deviation of 4.35% cannot meet all the engineering requirements. According to [26], the spline function consists of two parts. One of them is the residual norm, which ensures that the result is close to the evaluation of profile. The other part is bending energy, which controls the smoothness of the spline curve. Lagrange parameters control the tradeoff between the curve fitting and bending. According to the principle of the variational method [27], adding the first derivative term to the bending energy can improve the transmission characteristics of a spatial frequency function. This method can be used to approximate the transmission characteristics to those of the Gaussian filter. The approximating spline function is based on the cubic spline function and adds the first derivative term to the bending energy. The approximating coefficient is introduced to provide better control of the approximating degree of the transmission characteristics. The filtering characteristic of the approximating spline can be expressed as follows.

$$\begin{array}{l}G\left(\omega \right)=\frac{1}{1+2\tau \mu \left(1-\cos \omega \right)+4\mu {\left(1-\cos \omega \right)}^2}\\ \\ \mu ={(\sqrt{5}-1)}^2/64{\sin }^4(\pi \Delta x/{\lambda }_c)\end{array}$$

(1)

where τ denotes the approximating coefficient, and μ denotes the Lagrangian parameter, and Δx denotes the sampling interval, and λ_c denotes the cut-off wavelength. As shown in Figure 2a, the transmission characteristics of the approximating spline and Gaussian functions is enumerated when $\tau =1/\sqrt{\mu }$, $\tau =2/\sqrt{\mu }$, or $\tau =3/\sqrt{\mu }$. From Figure 2b, it can be seen that the transmission characteristics of the approximating spline filter are the closest to those of the Gaussian filter when $\tau =1/\sqrt{\mu }$ [13]. The maximum transmission characteristic deviation is 4.263%. This value is in line with the deviation range (−5.0%, 5.0%) of the simplified Gaussian filtering algorithm specified in the ISO 11562 [28]. Compared with the approximating spline recommended by international standard, the improved approximating spline is closer to the transmission characteristics of Gaussian filter as shown in Figure 3.

In order to obtain higher approximation accuracy, many approximating spline filter are connected in series. The transmission characteristic of the cascade approximating spline filter can be expressed as

$$\begin{array}{l}G\left(\omega \right)=\frac{1}{{\left(1+2\tau \mu \left(1-\cos \omega \right)+4\mu {\left(1-\cos \omega \right)}^2\right)}^n}\\ \\ \mu ={(\sqrt{4\sqrt[n]{2}-3}-1)}^2/64{\sin }^4(\pi \Delta x/{\lambda }_c)\end{array}$$

(2)

where n denotes the cascade order. The transmission characteristic deviation for different cascade orders is presented in Figure 4, where it can be seen that with the increase in the cascade order, the approximation degree also increases. However, the improvement in the approximation accuracy results in higher calculation cost. Therefore, considering the effectiveness, the cascade order should be controlled in a reasonable range. The deviation is less than 1.0% when n = 2, so it is enough to use the second-order joint approximating spline filter in practice.

In the three-dimensional surface measurement, a method similar to the separability of the Gaussian filter is used to realize a two-dimensional approximating spline filter. Two groups of one-dimensional approximating spline filters process the original data along the x-axis and y-axis. First, the filter is applied m times in the x-direction (the same as the original data row number), and then n times in the y-direction (the same as the original data column number) [29]. Two-dimensional approximating spline filtering characteristics are defined as follows:

$$G\left(\omega \right)=\frac{1}{\left\{1+2\tau \mu \left(1-\cos {\omega }_x\right)+4\mu {\left(1-\cos {\omega }_x\right)}^2\right\}\times \left\{1+2\tau \mu \left(1-\cos {\omega }_y\right)+4\mu {\left(1-\cos {\omega }_y\right)}^2\right\}}$$

(3)

The anisotropy of the normal spline filter in orthogonal direction limits its application in 3D surface roughness evaluation. A two-dimensional isotropic spline filter can be built by superimposing a one-dimensional cascade approximating spline filter in the x and y directions, and its amplitude-frequency characteristics are close to those of the Gaussian filter. In the one-dimensional cascade approximating filter, when n = 2, the maximum deviation is less than 1.0%, which represents an ideal balance between the approximation accuracy and the computational efficiency. Similarly, as shown in Figure 5, when n = 2, the maximum deviation between the amplitude-frequency characteristics of the two-dimensional cascaded approximating spline filter and the two-dimensional Gaussian filter is less than 1.5%, which can meet most of the application requirements.

3. Filtering Algorithm

3.1. Gauss Elimination

An accurate and efficient solution to the approximating spline equation plays an important role in surface topography filtering. However, the solution to the spline function in engineering technology is often obtained by solving linear algebraic equations. As shown in Equation (4), the approximate spline partial differential equation is transformed into a system of linear algebraic equations by means of the difference method and derivation. As shown in Equation (5), A represents the coefficient matrix of the approximating spline. It can be found that the coefficient matrix of approximating spline filter is a kind of sparse matrix (that is, the number of non-zero elements in the matrix is far less than the total number of matrix elements). A is consistent with the order of the original data. The principle of the Gauss elimination method is to eliminate the unknown number step-by-step and transform the linear equation system Aw = z into an equivalent triangular linear equation system. In order to save the computer storage capacity, three one-dimensional arrays that are denoted as a, b, and c were used in this work to store five diagonal elements of matrix A in MATLAB as shown in Equation (6). A is transformed into an upper triangular matrix by row transformation; namely, the second row represents a product of the difference between the second and first rows and b₁/a₁ and the first-column element in the second row becomes zero. The third row represents a product of the difference between the third and first rows and c₁/a₁, and the first-column element of the third row becomes zero. The same operation is performed on all matrix rows, and the same operation is performed on all rows:

$$\begin{array}{c}{A}_{1,1}{w}_1+A_{1,2}{w}_2+\cdots \cdots +A_{1,n}{w}_n=z\\ A_{2,1}{w}_1+A_{2,2}{w}_2+\cdots \cdots +A_{2,n}{w}_n=z\\ ⋮\\ A_{n,1}{w}_1+A_{n,2}{w}_2+\cdots \cdots +A_{n,n}{w}_n=z\end{array},w=\left[\begin{array}{c}{w}_1\\ w_2\\ w_3\\ ⋮\\ ⋮\\ w_{n-1}\\ w_n\end{array}\right],z=\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\\ ⋮\\ ⋮\\ z_{n-1}\\ z_n\end{array}\right],$$

(4)

$$A=\left[\begin{array}{ccccccc}1+\mu +\sqrt{\mu }& -2\mu -\sqrt{\mu }& \mu & & & & \\ -2\mu -\sqrt{\mu }& 1+5\mu +2\sqrt{\mu }& -4\mu -\sqrt{\mu }& \mu & & & \\ \mu & -4\mu -\sqrt{\mu }& 1+6\mu +2\sqrt{\mu }& -4\mu -\sqrt{\mu }& \mu & & \\ & \ddots & \ddots & \ddots & \ddots & \ddots & \\ & & \mu & -4\mu -\sqrt{\mu }& 1+6\mu +2\sqrt{\mu }& -4\mu -\sqrt{\mu }& \mu \\ & & & \mu & -4\mu -\sqrt{\mu }& 1+5\mu +2\sqrt{\mu }& -2\mu -\sqrt{\mu }\\ & & & & \mu & -2\mu -\sqrt{\mu }& 1+\mu +\sqrt{\mu }\end{array}\right],$$

(5)

$$A=\left[\begin{array}{ccccccc}{a}_1& b_1& c_1& & & & \\ b_1& a_2& b_2& c_2& & & \\ c_1& b_2& a_3& b_3& c_3& & \\ & \ddots & \ddots & \ddots & \ddots & \ddots & \\ & & c_{n-4}& b_{n-3}& a_{n-2}& b_{n-2}& c_{n-2}\\ & & & c_{n-3}& b_{n-2}& a_{n-1}& b_{n-1}\\ & & & & c_{n-2}& b_{n-1}& a_n\end{array}\right],$$

(6)

In short, the Gauss order elimination method is used to transform the pentadiagonal matrix A into an upper triangular matrix as shown in Equation (7). The last waviness w_n is solved, and then the remaining waviness is obtained by back substitution. As shown in Figure 6, Gaussian elimination method includes two steps which are elimination and back substitution. This work estimates the calculation times of the two steps. The elimination of lower triangle elements b_i requires one division, two multiplications, and two subtractions for each element. For b_n−1, there are only two elements in the (n − 1)th row, so computer repeats it three times to eliminate b_n−1. However, when c_i is eliminated, it needs one division, one multiplication, and one subtraction. At the same time, the original data z participates in the elimination process, but the first element of these data is not included in the elimination process, while the second element is eliminated only once, and each following element starting from the third element is eliminated twice. Each time the original data are eliminated, the process includes one division, one multiplication, and one subtraction. In fact, in a computer implementation, the elimination process of the original data and the elimination process of matrix A share the division operation at the same time. Therefore, the elimination of original data considers only the latter two operations. In general, the calculation time of the Gaussian elimination is mainly the time of the elimination process. When a structure (Aw = z (1), Aw = z (2), … , Aw = z (n)) exists in a system, the elimination process needs to be repeatedly calculated for the change in z. This means that the Gaussian elimination method for 3D surface filtering includes complicated computational steps. When a small diagonal element is divisor, it will lead to the order of magnitude growth of other elements and the diffusion of rounding error, the order of magnitude of other elements increases, and the diffusion of rounding error is caused, making the filter unstable.

$$A=\left[\begin{array}{ccccccc}{A}_{1,1}& A_{1,2}& A_{1,3}& & & & \\ & A_{2,2}& A_{2,3}& A_{2,4}& & & \\ & & A_{3,3}& A_{3,4}& A_{3,5}& & \\ & & & \ddots & \ddots & \ddots & \\ & & & & A_{n-2,n-2}& A_{n-2,n-1}& A_{n-2,n}\\ & & & & & A_{n-1,n-1}& A_{n-1,n}\\ & & & & & & A_{n,n}\end{array}\right],$$

(7)

3.2. LU Decomposition

The coefficient matrix of the approximating spline equation has symmetric positive definite property. As shown in Equation (8), a positive definite matrix A can be decomposed into unit lower triangular matrix L and unit upper triangular matrix U by rewriting the Gauss elimination method into a compact form. Thus, the problem of solving Aw = z is transformed into the problem of solving Ly = z, Uw = y. The advantage of transforming the system of equations into the matrix form is that the matrix algorithm can be used to decompose the spline equations into simple parts, thus improving the computational efficiency. The flow chart of LU decomposition method is realized on the computer as shown in Figure 7. The elimination process of LU decomposition method uses the Gaussian elimination method, and generates the upper triangular matrix U. At the same time, when the multipliers are eliminated, they are placed in a lower triangular matrix, and the main diagonal elements are set to one to generate the corresponding lower triangular matrix L:

$$A=\left[\begin{array}{ccccc}1& & & & \\ L_{2,1}& 1& & & \\ L_{3,1}& L_{3,2}& 1& & \\ & \ddots & \ddots & \ddots & \\ & & L_{n,n-2}& L_{n,n-1}& 1\end{array}\right]\left[\begin{array}{ccccc}{U}_{1,1}& U_{1,2}& U_{1,3}& & \\ & U_{2,2}& U_{2,3}& \ddots & \\ & & U_{3,3}& \ddots & U_{n-2,n}\\ & & & \ddots & U_{n-1,n}\\ & & & & U_{n,n}\end{array}\right],$$

(8)

When the computer realizes matrix U, the corresponding matrix L is also formed. Moreover, the original data z do not participate in the elimination process, and the algorithm performs elimination only once. Thus, the whole elimination process becomes mainly the formation of matrix U. The elimination includes approximately 11 operations. Then, computer need to do two back substitution processes for z. Compared to the Gaussian elimination, although there is an additional back-substitution, this additional calculation replaces the part of the calculation that is saved by the absence of z in the elimination process. For a set of raw data z, the LU decomposition method and the Gaussian elimination method have the same calculation cost. If A = LU decomposition is adopted, it is quite convenient to use LU decomposition method to solve some equations such as Aw = z (1), z (1), …, z(1). There is no need for a double elimination process, and each solution Aw = z requires only two back substitution. However, the LU decomposition still adopts the Gaussian elimination process, which can cause numerical instability.

3.3. LDLT Decomposition

At present, the square root method L^TL that is the LU decomposition method is commonly used in computers. In this method, symmetric positive definite equations are solved by trigonometric decomposition of symmetric positive definite matrices. This method splits the coefficient matrix A into two symmetric trigonometric matrices A = L^TL. When elements of L are determined, the elements of L^T can be obtained. In the process of computer processing, only non-zero elements of L can be stored. Additionally, the order of magnitude of L does not increase during the decomposition, and the diagonal elements are always positive. Therefore, the square root method is a numerically stable decomposition method. However, through the analysis of the calculation process, it is found that the back-generation process of the algorithm takes almost half of the whole solving-process time, and the symmetry and sparsity of the decomposed matrix are not fully utilized. In addition, the decomposition process requires the square root calculation, which can cause precision loss. In order to avoid calculating the square root, the improved square root method (LDLT) can be used to decompose the sparse matrix A. The decomposition formula is shown in Equation (9).

$$A=\left[\begin{array}{ccccc}1& & & & \\ L{b}_1& 1& & & \\ L{c}_1& L{b}_2& 1& & \\ & \ddots & \ddots & \ddots & \\ & & L{c}_{n-2}& L{b}_{n-1}& 1\end{array}\right]\left[\begin{array}{ccccc}{d}_1& & & & \\ & d_2& & & \\ & & d_3& & \\ & & & \ddots & \\ & & & & d_n\end{array}\right]\left[\begin{array}{ccccc}1& L{b}_1& L{c}_1& & \\ & 1& L{b}_2& \ddots & \\ & & 1& \ddots & L{c}_{n-2}\\ & & & \ddots & L{b}_{n-1}\\ & & & & 1\end{array}\right],$$

(9)

The LDLT decomposition method decomposes matrix A into two symmetric triangular matrices and one diagonal matrix. Therefore, the interference caused by rooting operation is avoided and improves the stability of the filtering result. In addition, the sparse matrix A is a symmetric matrix. Compared with the square root method, the LDLT decomposition method simplifies the program design. As shown in Figure 8, it is the program design flow chart of LDLT decomposition method.

3.4. Thomas Method

In engineering practice, Thomas method has been often used to solve ordinary differential equation, heat conduction equation, and many other equations [21]. These problems tend to be tridiagonal matrices, and the spline matrix is a pentadiagonal matrix. In this work, the spline matrix equations are derived by the direct triangulation of matrices. The decomposition is as follows:

$$A=\left[\begin{array}{ccccc}{l}_1& & & & \\ m_1& l_2& & & \\ s_1& m_2& l_3& & \\ & \ddots & \ddots & \ddots & \\ & & s_{n-2}& m_{n-1}& l_n\end{array}\right]\left[\begin{array}{ccccc}1& p_1& q_1& & \\ & 1& p_2& \ddots & \\ & & 1& \ddots & q_{n-2}\\ & & & \ddots & p_{n-1}\\ & & & & 1\end{array}\right],$$

(10)

Similar to LU decomposition method, pursuit method decomposes sparse matrix A into upper triangular matrix and lower triangular matrix. The program design flow chart of Thomas method is shown in Figure 9.

The computational efficiency of the above matrix algorithm was evaluated. The test surfaces with different data lengths (Gaussian random distribution surfaces) were generated by MATLAB software, and their consumption times of profile and areal filtering were compared. The two-level approximating spline filter was used as an experimental object. The cut-off wavelength of the long wave was 0.8 mm, and the sampling length was 0.01 mm. In order to reduce the error, under the condition of keeping the same data length, the surface profiles of five groups of different data were generated, and the average of the five filtering times was taken for comparison value. As previously mentioned, different matrix algorithms are written into functions and are called by MATLAB for filtering, and the time consumption is calculated by the running time of each function by the profiler of MATLAB. The uncertainty of a computer itself is not taken into account. The filtering times of 1000, 10,000, 100,000, 1,000,000, and 5,000,000 data point by the matrix algorithm are given in

Table 1. Time consumption results of the profile filtering test.

Matrix Algorithm	Profile Size (Points)
	1000	10,000	100,000	1,000,000	5,000,000
Gauss elimination	0.0202 s	0.0286 s	0.105 s	0.887 s	4.411 s
LU decomposition	0.0202 s	0.0442 s	0.150 s	1.576 s	7.197 s
LDLT	0.0216 s	0.0392 s	0.160 s	1.182 s	5.283 s
Thomas	0.0246 s	0.0396 s	0.179 s	1.327 s	5.939 s

. The results in Table 1 show that the Gaussian elimination method has the shortest time and the highest computational efficiency at the same data length among all the methods. However, when the amount of data is small, the time difference between the methods is small. With the increase in the data amount, the advantages of the Gaussian elimination and improved square root method become more apparent. Especially when facing millions and more data points, the Gaussian elimination method shows higher computational efficiency than LU method. In areal filtering, 0.8 mm wavelength was used in the orthogonal direction, and the sampling interval was 0.01 mm. The obtained filtering times are given in

Table 2. Time consumption results of the areal filtering test.

Matrix Algorithm	Surface Size (mm)
	500 × 500	2000 × 2000	6000 × 6000	10,000 × 10,000
Gauss elimination	0.468 s	5.221 s	46.239 s	200.224 s
LU decomposition	0.164 s	1.772 s	17.803 s	55.935 s
LDLT	0.183 s	1.754 s	16.691 s	50.506 s
Thomas	0.178 s	1.715 s	16.589 s	49.953 s

, where it can be seen that the Gaussian elimination method, which performed well in profile filtering, takes too long time in areal filtering. This is because the Gaussian elimination method needs to recalculate the elimination process in the face of a new column of data, which leads to complex and time-consuming calculations. The LDLT method still achieves high computational efficiency. In terms of computational efficiency, the matrix algorithm achieves a satisfactory result. Thus, the Gauss elimination is time-consuming and increases the calculation cost in areal filtering, while the LDLT and Thomas methods achieve high computational efficiency.

Another important factor that should be considered is the computational accuracy. In the matrix decomposition process, rounding errors accumulate. Therefore, the decomposition method of the approximate spline matrix directly affects the error of roughness parameters. As shown in Figure 10a, the turning surface was selected for the experiment. The Form Talysurf PGI^® 1500S from Taylor Hobson was used to obtain the turning surface profile. The type of Form Talysurf stylus is P57340, and tip radius of the stylus is 2 μm. The sampling interval of Form Talysurf is 0.125 μm and the measuring speed is 0.25 mm/s. LS-line fitting was conducted to remove the form deviation, as shown in Figure 10b. The roughness parameter Ra and Rq were chosen for characterization and calculated according to the ISO 4288 [30]. The sampling length was set as 2.5 mm in this experiment, the evaluation length was 12.5 mm, and the cut-off wavelength was 2.5 mm which was numerically equivalent to the sampling length. In order to better compare the differences brought by different algorithms, four algorithms were applied to the standard cubic spline. Gaussian filter and spline filter in commercial software Mountains^® were taken as the reference in this experiment. The roughness parameter Ra, Rq were obtained as shown in

Table 3. Ra, Rq obtained by the Mountains^® software and the four algorithms.

Methods	Ra (μm)	Rq (μm)
Gaussian filter (Mountains®)	1.490326514	1.759749942
Spline filter (Mountains®)	1.489708718	1.759108632
Gauss elimination	1.490295906	1.759401134
LU decomposition	1.490295906	1.759401134
LDLT	1.490299845	1.759401523
Thomas	1.490296286	1.759401445

. It can be calculated that the errors of parameter values Ra and parameter values Rq obtained by Gaussian filter and spline filter in the evaluation software were 0.06178% and 0.06413% respectively. This was due to the difference of filtering algorithm and transmission characteristics between the two filters. As shown in

Table 4. The errors between the roughness values obtained by the four algorithms and the roughness values obtained by the spline filter in the Mountains^® software.

Methods	The Errors of Ra	The Errors of Rq
Gauss elimination	0.05872%	0.02925%
LU decomposition	0.05872%	0.02925%
LDLT	0.05911%	0.02929%
Thomas	0.05876%	0.02928%

, the roughness parameter value obtained by the cubic spline filter in the software was different from the results of the cubic spline filter realized by the four algorithms. In fact, the difference was caused by the implementation algorithm of spline filter. In 2011, Mathia et al. pointed out that the calculation errors of roughness parameters after using a spline filter was bigger than using a Gaussian regression filter [31], which further proved the evaluation difference between the Gaussian filter and spline filter in the software Mountains^®. Therefore, this study takes the output value of the Gaussian filter as the final reference. Compared with the spline filter in the software, the implementation algorithm reduces the evaluation difference between the spline filter and the Gaussian filter as shown in

Table 5. The errors between the roughness values obtained by the four algorithms and the roughness values obtained by the Gaussian filter in the Mountains^® software.

Methods	The Errors of Ra	The Errors of Rq
Gauss elimination	0.003061%	0.03488%
LU decomposition	0.003061%	0.03488%
LDLT	0.002667%	0.03484%
Thomas	0.003023%	0.03485%

. The roughness value obtained by the LDLT algorithm is the closest to the Gaussian filter, and the error were only 0.002667% and 0.03484%.

The roughness values obtained by Gauss elimination method and LU decomposition method are the same. In fact, the LU decomposition is based on the idea of the Gaussian elimination, which is highly affected by the rounding error. Therefore, with the increase in the matrix order, the error accumulates significantly. The LDLT avoids the square root operation and further reduces the error accumulation. The result shows that the roughness value obtained by the LDLT is the closest to the reference value. In essence, Thomas method is the Gaussian elimination method without selecting the principal elements. When the main diagonal elements are dominant, the calculation process of the Thomas method will not introduce a large increase in the order of magnitude of the intermediate results and will not cause a severe accumulation of rounding errors. The approximate spline matrix is not a strictly diagonally dominant matrix. Therefore, the catch-up method will cause serious error accumulation in the face of massive data.

With the development of machining technology, the precision of surface roughness measurement has been becoming more and more demanding. Therefore, the detection accuracy and filtering accuracy of an instrument are highly important. Considering the computational efficiency and accuracy, the Gaussian elimination is a good choice for profile filtering when processing a small amount of data. However, LU decomposition and Gaussian elimination cannot meet the precision requirement in the face of massive data. In that case, the LDLT method should be preferred. In the areal filtering, the LDLT method represents a good choice for solving an approximating spline matrix because it has high computational efficiency and filtering precision.

4. Robust Treatment

The reliability of reference data is an important index to characterize surface roughness. When the surface contains outliers, such as scratches, grooves, and missing surface data, filter reference is distorted. The robustness estimation weight function is introduced into the filter to enhance the robustness against the abnormal signals. Also, different weight functions may produce completely different results when processing the same measured signal, thus increasing the difference in surface evaluation between different filters. In 2010, the robust Gaussian regression filter was included in the ISO 16610-31. Also, it has been shown that the combination of the Tukey estimation function and the Gaussian regression filter can improve the robustness of the Gaussian regression filter [18]. However, it increases the time and calculation costs. In the same year, the ISO 16610-32 stipulated a robust spline filter with p = 1(the L1-norm based Robust Spline Filter). Based on the idea of a variational robust spline filter, a compromise function between the sum of the absolute values of the resists and the approximate bending energy of the data has been constructed. However, this function requires many cyclic iterative operations to be solved. It has been proved that the standard robust spline filter cannot achieve satisfactory results in the surface evaluation [23]. Therefore, it is necessary to design a robust filter, which not only has high filtering accuracy and computational efficiency, but also can reduce the evaluation difference between robust spline filter and robust Gaussian regression filter, so as to realize the unification of the two evaluation standards.

In this study, the approximating spline is used as a base filter, and several typical weight functions are selected for comparison. The matrix definition of the robust approximating spline is as follows:

$$\left(\rho +\tau \mu P+\mu Q\right)w=\rho z,$$

(11)

where ρ denotes a robust weight function, and τ denotes the approximating coefficient, and μ denotes the Lagrangian parameter. P and Q are matrices related to boundary conditions. First, the filter results are obtained by using the approximating spline filter, and then the results are substituted into the weight function to obtain the weights of each point. Then, the calculated weights are assigned to the points with small residual errors, and small weights are assigned to the points with large residual errors. Finally, the weights of each point are substituted into the robust approximating spline, and the robust filter values are iteratively calculated.

According to M-estimation, different robust weight functions may lead to diametrically different results when one and the same signal observed is processed. An impertinent choice of weight functions can sometimes make the estimating method invalid. In practice, it is rather difficult to grasp the exact distribution of the measured signal. To make approximating spline filter adaptable to widely practical application, several typical weight functions are selected and compared with each other [32].

(1)

WLAV

$$\rho \left(\upsilon \right)=\frac{1}{\upsilon },$$

(12)

where v = r/cMAD, MAD = |med(r)|.

(2)

Tukey

$$\rho \left(\upsilon \right)=\{\begin{array}{ll}\left({\left(1-{\upsilon }^2\right)}^2\right),& \left|\upsilon \right|\le 1\\ 0,& \left|\upsilon \right|>1\end{array},$$

(13)

where v = r/cMAD, and c denotes the modulation constant used to balance robustness and efficiency; c = 4.44 is adopted as an international standard.

(3)

ADRF

$$\rho \left(\upsilon \right)=\{\begin{array}{ll}1,& \left|\upsilon \right|\le a\\ a/\left|\upsilon \right|,& \left|\upsilon \right|\in \left(a,b\right]\\ b/{\upsilon }^2,& \left|\upsilon \right|\in \left(b,c\right]\\ 0,& \left|\upsilon \right|>c\end{array},$$

(14)

where v = r/s, s denotes a scale parameter, and a, b, and c denote the fine-tuning parameters.

(4)

Andrews

$$\rho \left(\upsilon \right)=\{\begin{array}{ll}\frac{1}{\pi \upsilon }\sin \left(\pi \upsilon \right),& \left|\upsilon \right|\le 1\\ 0,& \left|\upsilon \right|>1\end{array},$$

(15)

(5)

Hampel

$$\rho \left(\upsilon \right)=\{\begin{array}{ll}1,& \left|\upsilon \right|\le a\\ a/\left|\upsilon \right|,& \left|\upsilon \right|\in \left(a,b\right]\\ a\frac{c-\left|\upsilon \right|}{\left(c-b\right)\left|\upsilon \right|},& \left|\upsilon \right|\in \left(b,c\right]\\ 0,& \left|\upsilon \right|>c\end{array},$$

(16)

(6)

Huber

$$\rho \left(\upsilon \right)=\{\begin{array}{ll}1,& \left|\upsilon \right|\le c\\ \frac{c}{\upsilon },& \left|\upsilon \right|>c\end{array},$$

(17)

where c is the threshold value, and its value is 1.345.

5. Experimental Results

A honing surface with an obvious deformity was used in the experiments as shown in Figure 11a. The data were sampled at 4-μm interval, and 3854 data points were obtained. In order to evaluate the surface roughness accurately, the polynomial fitting method was used to remove the large shape components of the surface. The same cut-off wavelength (0.8 mm) was used in all filtering processes. The results of the robust Gaussian regression filter were used as a reference. As shown in Figure 11b, the waviness curve after using approximating spline was distorted at a large deep valley, resulting in a large distortion of the roughness. After different robust processing, the waviness profile of the filter output was smooth and natural, and the influence of abnormal features was suppressed. It can be seen that robust processing can obtain more realistic and reliable filtering values even in the case of surface profiles with large outliers, such as deep valleys or high peaks. The profile processed by the robust Gaussian regression filter was used as the reference data. A robust second-order Gaussian regression filter with a cut-off wavelength of 0.8 mm was used to perform the filtering operation. The evaluation length was set to be the same as the total length of the profile. Figure 12 shows the similarity between the waviness profile processed by different robust approximating spline filters and robust Gaussian regression filter. It can be seen that for most of the profiles, the profiles processed by the Tukey, ADRF, and Andrews are in good agreement with the reference curve, while those of the WLAV, Hampel, and Huber show significant differences.

The typical roughness parameters values of using different filtering methods are presented in

Table 6. Comparison of roughness parameter values of using different filtering methods.

	Not Robust	Tukey	ADRF	WLAV	Hampel	Andrews	Huber	Robust Gaussian
Ra (μm)	1.0945	1.1652	1.1618	1.1422	1.0763	1.1590	1.1325	1.1518
Rq (μm)	1.5434	1.6665	1.6626	1.6375	1.6650	1.6588	1.6210	1.6532
Rp (μm)	2.7425	3.2163	3.4880	3.2760	2.5425	2.9162	3.1268	3.6277
Rv (μm)	9.0952	10.2168	10.1310	9.9840	10.4705	10.2113	9.9364	10.1526
Rt (μm)	11.8377	13.4331	13.6191	13.2601	13.0130	13.1275	13.0632	13.7802
Rku	7.7288	8.8127	8.6766	8.7273	10.2548	8.9193	8.7533	8.8800
Rsk	−1.7914	−2.0305	−1.9900	−1.9887	−2.5561	−2.0567	−2.0004	−2.0167
Rsm (mm)	0.3310	0.3509	0.3594	0.3508	0.5251	0.3384	0.3344	0.3604
k		3	2	17	4	4	4	3
time	0.008 s	0.027 s	0.017 s	0.054 s	0.028 s	0.024 s	0.022 s	0.276 s

. First, the amplitude parameters were investigated. Ra (the arithmetical mean deviation of the profile) is the most commonly used 2D roughness parameter. After robust treatment, Rp (the maximum peak height of the profile) and Rv (the maximum valley depth of the profile) increased. Rt is the total height of the profile and equal to the sum of Rp and Rv. Meanwhile, Rq (the root mean square deviation of the profile) and Rku (the kurtosis of the profile) also increased accordingly. Rsk (the skewness of the profile) showed a negative decreased trend. The results show that compared to the unrobust approximating spline filter, the robust processing can effectively suppress the influence and propagation of the malformed characteristic signal and reduce the distortion of the median line of the filter. The roughness profile retains the characteristics of the deep valley and high peak well. Therefore, robust filtering can more truly reflect the low-frequency waviness and high-frequency roughness components and is more conducive to the evaluation of surface morphology. The roughness parameters processed by Huber and Hampel methods differed greatly from the robust Gaussian processing values. Furthermore, the spacing parameter Rsm (mean width of the profile elements) increased after robust treatment as well. The results show that the robust processing can more truly reflect the density of the peak and valley of the surface profile compared to the unrobust approximating spline filter. The results were further verified, and the roughness parameters of Tukey, ADRF, and Andrews were further analyzed as shown in Figure 12. The relative errors between the roughness parameters obtained by Tukey, ADRF, and Andrews and the reference data are given in

Table 7. Relative errors of roughness parameters obtained by different filtering methods.

Filtering Methods	Relative Errors of the Roughness Parameters
	Ra	Rq	Rp	Rv	Rt	Rku	Rsm
Tukey	0.012%	0.80%	11.3%	0.63%	2.58%	12.9%	2.64%
ADRF	0.008%	0.57%	3.85%	0.21%	1.17%	0.76%	0.28%
Andrews	0.006%	0.34%	19.6%	0.58%	1.40%	1.40%	6.10%

. The relative errors of Rp parameters between Tukey, ADRF, and Andrews methods and the robust Gaussian regression filter were 11.3%,3.85%, and 19.6%, respectively; the relative errors of Rv parameters were 0.63%, 0.21%, and 0.58%, respectively; lastly, the relative errors of Rt parameters were 2.58%, 1.17%, and 4.74%, respectively. Generally speaking, the filtering effect of the ADRF method on the deep valley and high peak was more similar to that of the robust Gaussian regression filter. The relative errors of Rq parameters were 0.80%, 0.57%, and 0.34%, while the relative errors of profile kurtosis Rku parameters were 12.9%, 0.76%, and 1.40%. Thus, the results of the Tukey, ADRF, and Andrews methods in reflecting peak-valley symmetry are close to that of the robust Gaussian regression filter. The relative errors of Rsm parameters were 2.64%, 0.28%, and 6.10%. ADRF was closer to the robust Gaussian regression filter in reflecting peak-valley symmetry. Another important factor for filter is the iteration time. The numbers of iterations of the Tukey and Andrews methods were k = 3 (0.027 s) and k = 4 (0.024 s) respectively, and the calculation efficiencies of these methods were relatively low. The number of iterations of the ADRF method was k = 2 (0.017 s), but it achieved high computational efficiency. The experimental results show that the ADRF method has good filtering performance for honing profile. The evaluation result of the surface roughness is also close to that of using the robust Gaussian regression filter. Through a large number of experiments, it is found that the ADRF method is neutral when it is applied to the profile with regular form and Gaussian distribution. It shows good robustness when ADRF method is applied to the surface profile with outliers such as scratches and grooves, which further verifies the applicability of ADRF method in surface roughness extraction.

The constructed ADRF robust approximating spline filter was then applied to areal filtering. The step surface obtained by the Talysurf CCI measurement system is presented in Figure 13. The existing surface shape was deleted, and the surface details were retained. The waviness surface and rough surface obtained by the approximating spline filtering without robust treatment is presented in Figure 14. According to the output results, there was an obvious distortion in the step area, and a false peak phenomenon occurred in the vicinity of the deep valley. Therefore, the presence of steps will seriously affect the accurate evaluation of roughness.

The waviness profiles of the RGRF and ADRF outputs are almost the same, and both suppress the influence of the step as shown in Figure 15a and Figure 16a. It can be seen that the robust filtering not only produced small distortion in the step region and around it but also had small deformation at the boundary. The robust processing effectively extracted the low-frequency waviness component and eliminated the influence of deep valleys features, and the obtained roughness surface was more real and reliable. The corresponding height parameter values are given in

Table 8. Roughness parameter values after using the ADRF and RGRF filter.

Filtering Method	Sa (μm)	Sq (μm)	Sv (μm)	Sku	Ssk	time
RGRF	8.49	16.10	46.76	5.27	−2.26	49.492 s
ADRF	8.43	16.19	46.67	5.30	−2.28	0.316 s

. The parameters of Sa and Sq have little difference, therefore the overall height deviation of the surface was almost the same. The error of Sv was 9%. Both of them had the same treatment effect at the maximum deep valley and could eliminate the influence of the step. The Ssk parameters were all negative, indicating that there were corresponding deep valley area, and the relative error between them was small. From Figure 17, it can be seen that the material ratio curve was almost the same, but there was a small difference between the profile peak and the profile valley. In addition, the computational time of ADRF is only 0.316 s, and the computational efficiency is much higher than that of the robust Gaussian regression filter. The experimental results show that the ADRF robust approximating spline has good performance in areal filtering with deep valleys.

6. Conclusions

In this work, the performance of approximating spline filtering were compared and improved. The transmission characteristic deviation between the approximating spline and Gaussian filters is further reduced by the cascade approximating method. It is found that the deviation of the second-order approximating spline can meet the general measurement requirement of 5%. Further research shows that the two-dimensional approximating spline can produce approximately isotropic amplitude characteristics, which ensures that the attenuation degree of the signal is consistent in all directions when the signal is passing through the filter. This is beneficial to the application of approximating spline in areal filtering. In fact, the filtering equation of approximating spline can be transformed into the form of a matrix equation. Matrix factorization involves computational efficiency and numerical accuracy. The error caused by matrix factorization can directly affect the roughness parameter values. In order to find a suitable matrix algorithm, four direct methods of matrix decomposition were realized by the numerical analysis theory. Through profile and areal filtering, it is found that LDLT algorithm has the highest computational efficiency. The Gaussian filter in the commercial Mountains^® software was used as a reference in the accuracy comparison experiment. Compared with Gauss elimination method, LU decomposition method, and Thomas method, it is found that the roughness obtained by LDLT method is the closest to that of using Gaussian filter. At the same time, the algorithm can further reduce the evaluation difference between the spline filter and the Gaussian filter in the Mountains^® software. Finally, in order to improve the robustness of the approximating spline filter and reduce the evaluation difference with the robust Gaussian regression filter, six types of robust approximating spline filters were constructed and analyzed. The results of the experiment of honing profile and step surface show that the robustness of the ADRF filter is closer to that of the robust Gaussian filter when dealing with the outliers. At the same time, the number of iterations of the ADRF filter is also less than those of the other function filters. The ADRF robust approximating spline filter and the robust Gaussian regression filter have similar robust performance, and the evaluation results for the same surface have a little difference, which is helpful for realizing the unification of the two sets of international standards. In particular, the robust filter is neutral when facing a uniform surface, such as that of obeying the Gaussian distribution. Therefore, in practical application, the appropriate filter should be selected according to the surface characteristics.