Mathematical approach to the validation of field surface texture parameter software

To ensure mathematical traceability to analytical surface representations, the surface texture parameter calculations must also be performed analytically. To perform analytical, or symbolic, calculations computationally, a computer algebra system (CAS) is required. A CAS is a mathematical software package or module that can manipulate symbolic mathematical expressions directly, typically using a series of rules or look-up tables dictating mathematical procedures, without relying on numerical approximation. Such an approach delivers seemingly infinite precision results, providing the system is successfully able to obtain a mathematical expression, and ensures mathematical traceability throughout the operation. There are a variety of CASs available, and while some are open-source such as SageMath [20], the systems that are most popular are closed-source and commercial products. For the majority of the work in this paper, Wolfram Mathematica 11.1 was used as the primary CAS, with MAPLE 2017 also used to verify the results of Mathematica for first time calculations (subsequent calculations using a similar mathematical form did not require MAPLE verification) [21, 22]. Both Mathematica and MAPLE are among the most capable CASs available at the time.

The calculation of each parameter requires adhering to the definitions given in ISO 25 178-2 [7, 23]. In an ideal situation, this means substituting the surface equation for z(x, y) in the parameter definition and evaluating the result. The following sections give the mathematical operations used to calculate each parameter, and, where necessary, any alternative steps taken to circumvent limitations of the CAS used.

Before any parameter calculation can be performed, the surface function must be adjusted so that it has a mean-plane of zero. This is achieved by simply finding the mean height across the entire evaluation area of the surface and subtracting it from the surface function, given by

$$\begin{eqnarray}&&z(x,y)={z}_{0}(x,y)-\displaystyle \frac{1}{{A}_{\mathrm{eval}}}{\iint }_{{A}_{{\rm{R}}}}{z}_{0}(x,y){\rm{d}}x{\rm{d}}y,\end{eqnarray} \tag{ 1 }$$

where z₀(x, y) is used to represent the unadjusted surface function, the integration region, A_R, is the defined evaluation area of the surface, A_eval is the total evaluation area, and z(x, y) will represent the adjusted surface function in all future instances.

2.1.1. Sq, Ssk and Sku

The first three parameters, Sq, Ssk and Sku, represent the root mean square, skewness and kurtosis of the surface, respectively [7]. The definitions for each of these have similar forms, relying on a integral of some power function of z(x, y). Such an operation is well supported by most modern CASs, and so these parameters can be evaluated directly using the original definitions,

$$\begin{eqnarray}&&{S}_{{\rm{q}}}=\sqrt{\displaystyle \frac{1}{{A}_{\mathrm{eval}}}{\iint }_{{A}_{{\rm{R}}}}{z}^{2}(x,y){\rm{d}}x{\rm{d}}y},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{S}_{\mathrm{sk}}=\displaystyle \frac{1}{{S}_{\mathrm{sq}}^{3}}\left(\displaystyle \frac{1}{{A}_{\mathrm{eval}}}{\iint }_{{A}_{{\rm{R}}}}{z}^{3}(x,y){\rm{d}}x{\rm{d}}y\right)\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&{S}_{\mathrm{ku}}=\displaystyle \frac{1}{{S}_{\mathrm{sq}}^{4}}\left(\displaystyle \frac{1}{{A}_{\mathrm{eval}}}{\iint }_{{A}_{{\rm{R}}}}{z}^{4}(x,y){\rm{d}}x{\rm{d}}y\right).\end{eqnarray} \tag{ 4 }$$

In all surface texture parameter equations throughout this document, parameter sysmbols are given as a single letter with subscripts as needed, which differs from the parameter form given in text. The reason for this differentiation, in accordance with the ISO 25 178-2 definitions, is to avoid misinterpretation of compound letters as an indication of multiplication between quantities in equations [7]. As both Ssk and Sku rely on Sq in their definitions, Sq is calculated first, and the obtained expression substituted into the expression for Ssk and Sku. In the definitions given above, the double integral is performed over the defined area, A_R, which is the chosen domain of the continuous analytical expression z(x, y) over which the surface area is defined. A_R is assumed to be a rectangle defined by the interval [x_l, x_h] in the x direction and [y_l, y_h] in the y direction.

2.1.2. Sp, Sv and Sz

The next three parameters, Sp, Sv and Sz, are conceptually simple in that they relate to the largest peak/valley height values on the surface [7]. Sp corresponds to the largest peak height within the defined area, Sv corresponds to the magnitude of largest valley depth, and Sz is the total distance between the largest peak height and largest valley depth. For a discrete dataset, maximum peak/valley heights are obtained trivially by sorting the discrete heights and obtaining the highest and lowest values.

For a continuous analytical surface, obtaining maximum peak/valley heights symbolically becomes more complex. The parameters require calculating the global maxima and minima of the function within the specified region. Several optimisation techniques can potentially obtain values corresponding to only local minima, not global minima, and so care must be taken to avoid this. Modern CASs have many global optimisation functions available, each with strengths and weaknesses, which can be adopted to find these global maxima and minima [24].

A more manual approach to obtaining the largest peak and valley can be performed by finding all local minima and maxima points within the defined area using first and second derivatives, and sorting the results. Local minima/maxima are found at values of z(x, y) with zero gradient values in both x and y,

$$\begin{eqnarray}&&\displaystyle \frac{\partial z(x,y)}{\partial x}=0,\quad \displaystyle \frac{\partial z(x,y)}{\partial y}=0.\end{eqnarray} \tag{ 5 }$$

Furthermore, local maxima are found at locations where the second derivative of the function z(x, y) is less than zero, and local minima are found where the second derivative is greater than zero, given by

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}z(x,y)}{\partial {x}^{2}}\lt 0,\quad \displaystyle \frac{{\partial }^{2}z(x,y)}{\partial {y}^{2}}\lt 0\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}z(x,y)}{\partial {x}^{2}}\gt 0,\quad \displaystyle \frac{{\partial }^{2}z(x,y)}{\partial {y}^{2}}\gt 0,\end{eqnarray} \tag{ 7 }$$

respectively. By solving these equations simultaneously, values for x and y can be found that correspond to each local maxima and minima. By evaluating each of these points in z(x, y), expressions for surface height can be found, which can then be sorted to obtain the largest peak and valley.

Using the method of minimising the derivatives of a surface, any values that are not located at a zero-gradient in both x and y will not be found. This can miss the true Sp or Sv values if they happen to be located on the boundary of the defined area, where the zero-gradient for both x and y is not necessarily satisfied. Visual inspection of the surface can help to identify if this is the case for any individual scenario. If so, simple one-dimensional gradient analysis along the boundary line, at z(x, y_l) for example, can be performed for each of the four boundary edges to find any zero-gradient locations, which can then be evaluated to find the expressions for the heights. Finally, direct evaluation of z(x, y) at each for the four boundary corners can be performed to ensure any potential maxima/minima are identified there, as these points may not satisfy the boundary gradient analysis conditions, for similar reasons as before.

2.1.3. Sal and Str

Sal and Str are part of a subset of field parameters called spatial parameters [7]. Spatial parameters address the similarities of a surface with itself if spatially translated by a certain amount in the x and y directions. Such an analysis is valuable in identifying the uniformity of the surface. Sal assesses how abruptly a surface's height changes, and Str gives the ratio of fast-changing directions to slow-changing directions, indicating the presence of any directionality of the surface texture, such as grinding marks.

These parameters are calculated on the autocorrelation function of the surface, f_ACF(t_x, t_y). The autocorrelation function describes the degree of agreement between the surface function, z(x, y), and the same surface translated by t_x and t_y in the x and y directions, respectively. The definition for the autocorrelation function is given in ISO 25 178-2 as

$$\begin{eqnarray}&&{f}_{{ACF}}({t}_{x},{t}_{y})=\displaystyle \frac{{\iint }_{{A}_{{\rm{R}}}}z(x,y)\,z(x-{t}_{x},y-{t}_{y})\,{\rm{d}}x{\rm{d}}y}{{\iint }_{{A}_{{\rm{R}}}}z(x,y)\,z(x,y)\,{\rm{d}}x{\rm{d}}y}.\end{eqnarray} \tag{ 8 }$$

The denominator of equation (8) is used to normalise the autocorrelation function to give a value of 1 at t_x = 0 and t_y = 0. Calculation of the autocorrelation function requires evaluation of integrals similar to Sq, Ssk and Sku, and thus can be performed using a CAS.

Sal and Str are obtained by finding the distance it takes for f_ACF to decay from 1 to a specified value, s, defaulting at $s=1/5$ in ISO 25 178-3 [25]. This is achieved by finding the expression for the contour line of

$$\begin{eqnarray}&&{f}_{{ACF}}({t}_{x},{t}_{y})=1/5.\end{eqnarray} \tag{ 9 }$$

Sal is defined by the shortest distance from the origin, (0,0), to the contour line. The distance from the origin to a point on the contour is given by

$$\begin{eqnarray}&&\mathrm{Distance}=\sqrt{{t}_{x}^{2}+{t}_{y}^{2}}.\end{eqnarray} \tag{ 10 }$$

By rearranging equation (9) in terms of t_y, and substituting into equation (10), it is possible to obtain an expression in one dimension that gives the distance from the contour to the origin. This distance can then be minimised within the defined surface area to find the shortest distance to the origin, corresponding to Sal. Maximum points can be found in a similar fashion, and the ratio of shortest distance to longest distance gives the value for Str.

2.1.4. Sdq and Sdr

Sdq and Sdr are hybrid parameters, and related to the slope of the surface [7]. These parameters assess the steepness of features on the surface and can give valuable information for surfaces that could otherwise present similar Sa or Sq values. Hybrid parameters can be useful for the assessment of sealing performance and wettability, or any application that requires understanding of average surface gradient values.

These parameters require calculation of the surface gradients in the x and y directions, which can be obtained using partial derivatives as shown in section 2.1.2 for calculating Sp and Sv. Once the partial derivatives are calculated, Sdq can be found by evaluating the definition

$$\begin{eqnarray}&&\begin{array}{l}{S}_{{dq}}\\ =\,\sqrt{\tfrac{1}{{A}_{\mathrm{eval}}}{\iint }_{{A}_{{\rm{R}}}}\left[{\left(\tfrac{\partial z(x,y)}{\partial x}\right)}^{2}+{\left(\tfrac{\partial z(x,y)}{\partial y}\right)}^{2}\right]{\rm{d}}x{\rm{d}}y},\end{array}\end{eqnarray} \tag{ 11 }$$

which can be performed similarly to that discussed in section 2.1.1 for Sq, Ssk and Sku. Sdr has the definition

$$\begin{eqnarray}&&\begin{array}{l}{S}_{{dr}}\,=\,\displaystyle \frac{1}{{A}_{\mathrm{eval}}}\left[{\iint }_{{A}_{{\rm{R}}}}\left(\sqrt{\left[1+{\left(\tfrac{\partial z(x,y)}{\partial x}\right)}^{2}+{\left(\tfrac{\partial z(x,y)}{\partial y}\right)}^{2}\right]}-1\right){\rm{d}}x{\rm{d}}y\right],\end{array}\end{eqnarray} \tag{ 12 }$$

which, while appearing similar in composition to Sdq at first glance, poses a problem to CASs due to the square root of disparate terms within the integral being inseparable. Calculating Sdr using this definition, for the general case of an analytical surface expression of any form, requires the object within the integral to be evaluated as a single object, as it is inseparable. This particular inseparable evaluation leads to incalculable symbolic expressions and hanging run-times, due to seemingly endless cycles of techniques such as integration by parts, wherein an integral can be simplified by performing an integral on the derivative of some part of the original function. In this instance, the derived function would still be as complex as the original function, requiring further integration by parts and causing the recursive cycles. Equation (12) takes the form of a nonelementary integral, which is an integral of an elementary function whose solution is not elementary (an elementary function is a function comprised of a finite number of arithmetic operations on traditional functions such as powers, roots, trigonometric functions, exponential functions and logarithmic functions). The concept of a nonelementary integral is proven in Louiville's theorem [26, 27], and evaluations of nonelementary integrals require the use of either Taylor series expansion (to a required order), numerical integration or nonelementary 'special' functions (such as incomplete elliptical integrals or error functions). Both Taylor expansion and numerical integration are approximations, giving the result to a finite precision. Nonelementary special functions are advanced mathematical concepts that cannot be expressed using regular elementary functions, making it challenging to promote widespread use and adoption, and so will not be used here.

Practically, Sdr is the ratio between the planar area of evaluation of the surface, A_eval, and the total topographical area of the surface itself. The topographical area can be approximated by splitting the surface into a series of triangles, BCD, and summing the area of each of the triangles using the formula

$$\begin{eqnarray}&&\mathrm{Area}=\displaystyle \frac{1}{2}| \hat{u}\times \hat{v}| ,\end{eqnarray} \tag{ 13 }$$

where $\hat{u}$ is the vector that joins point B to point C, and $\hat{v}$ is the vector that joins point B to point D. By sampling a grid of height values from the surface expression, splitting each square defined by four height values into two triangles where the height values describe the vertices, and evaluating the cross product defined in equation (13), an expression can be derived for the approximation of the total topographical area of the surface,

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{Area} & = & \sum _{m=1}^{M}\sum _{n=1}^{N}\displaystyle \frac{\delta x\delta y}{2}\left(\sqrt{1+\tfrac{{\left({Z}_{n+1,m}-{Z}_{n,m}\right)}^{2}}{\delta {x}^{2}}+\tfrac{{\left({Z}_{n,m+1}-{Z}_{n,m}\right)}^{2}}{\delta {y}^{2}}}\right.\,+\left.\sqrt{1+\tfrac{{\left({Z}_{n+1,m+1}-{Z}_{n,m+1}\right)}^{2}}{\delta {x}^{2}}+\tfrac{{\left({Z}_{n+1,m+1}-{Z}_{n+1,m}\right)}^{2}}{\delta {y}^{2}}}\right),\end{array}\end{eqnarray} \tag{ 14 }$$

where N and M are the total number of vertices in the x and y directions, respectively, and δx and δy are the widths between adjacent vertices in the x and y directions, respectively. ${Z}_{n,m}$ describes the surface height values at each vertex,

$$\begin{eqnarray}&&{Z}_{n,m}=Z({x}_{n},{y}_{m}),\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{x}_{n}={x}_{0}+(n-1)\delta x\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{y}_{m}={y}_{0}+(m-1)\delta y.\end{eqnarray} \tag{ 17 }$$

Here, x₀ and y₀ are the coordinates of the starting corner of the surface area of interest. By using this method with an analytical surface expression, it is possible to define δx and δy to be any desired value. δx and δy can be reduced to the infinitesimal case, allowing a value for the total topographical area to be found that is equivalent to a symbolic evaluation. In practice, this would require seemingly infinite computational resources, and so compromises must be made, instead calculating to a precision that is appropriate for the required application. This approach is an alternative workaround to the limitations of calculating Sdr in the general case for any surface, however, it sacrifices the mathematical traceability by resorting to discrete, numerical methods. An alternative method, that calculates Sdr analytically for a specific case, in presented in section 4.4.

2.1.5. Sa

Sa is one of the most popular areal surface texture parameters in industry [28]. Sa is simply the mean deviation of height values on the surface from the mean plane and is broadly used in industry to define the 'roughness' of a surface, as a high Sa value would correspond to a larger spread between high and low points on a surface. However, due to the simplicity of the definition, the practical application of this metric, for example to understand lubricant retention, is limited, as surfaces with deterministic elements can skew the Sa value.

Mathematically, Sa is defined as the arithmetic mean of the absolute of the surface heights within the evaluation area,

$$\begin{eqnarray}&&{S}_{a}=\displaystyle \frac{1}{{A}_{\mathrm{eval}}}{\int }_{{A}_{{\rm{R}}}}| z(x,y)| {\rm{d}}x{\rm{d}}y.\end{eqnarray} \tag{ 18 }$$

Unfortunately, due to the absolute function, direct integration of the analytical surface expression is not possible.

All areas of the surface that lie below the mean plane are effectively mirrored to the positive z-axis, causing discontinuous edges at the points where the surface meets the mean plane. These discontinuities cannot be integrated as a single function. For simple cases, it is possible to separate the full discontinuous expression into a series of region-bound expressions and evaluate them all individually, however, a new expression is required for each region, along with information about the boundaries for each region (which are not necessarily linear, thus requiring line integrals). For more complex surfaces with larger numbers of mean plane crossing points, the formulation of each discrete term becomes prohibitively complex, particularly as the process would be bespoke for each individual surface expression, and not something easily implementable within a CAS. An alternative solution for simple parametric surfaces has been implemented in section 4.

A selection of ten simple analytical surfaces were defined for use in showcasing the calculation of field surface texture parameters. In order to highlight the variety of surface functions that can be used, three different types of function were selected:

• Three cosine functions, using the technique introduced in [19].
• Five polynomial functions, comprising of a cubic polynomial and four Bernoulli polynomials (two fourth order and two fifth order).
• Two Gaussian functions, using the natural exponential function, e^x.

To showcase the validity of this new approach, comparisons were performed with a selection of third-party parameter calculation software packages. A total of four different software packages were used, labelled A to D, which included a combination of commercial software and software developed by NMIs. Each software package was given a 700 × 700 high resolution SF surface. SDF dataset sampled from the surface expressions defined in above, where resolution is defined as the number of individual height values included in the surface dataset. This resolution was chosen because software D suffered from an upper limit on the file size of datasets of 10 MB, and a resolution of 700 × 700 was the highest resolution that reliably came under this limit, ensuring comparisons to be possible across all four software packages. No additional form removal or filtration operations were performed to enable sole assessment of the parameter value calculation methods.

The results of each of the third-party software tests are normalised relative to the mathematically-obtained parameter values to allow for more direct comparisons. Figure 1 shows the normalised values obtained for the Sq parameter, relative to the mathematical value, for a selection of eight analytical surfaces. Upon initial investigation, all software packages perform well in calculating Sq, however, closer inspection shows deviation from the mathematical value of the order of tenths of a percent, particularly for software B and D. This level of deviation could prove important when dealing with high precision applications, and helps metrology software developers to identify the aspects of their software that perform worse than competitors, highlighting key areas where their software requires improvement.

Zoom In Zoom Out Reset image size

Figures 2 and 3 show the results for the Sp and Sdq parameters, respectively. All software packages under test show good agreement with the mathematically-obtained values for Sp, which is expected due its simple definition of maximum peak height within the evaluation area. Discrete-based software can obtain the largest height value in a surface dataset easily with simple sorting algorithms for the height values (although this may not correspond to the true peak, depending on where the continuous surface was sampled). In contrast, the Sdq parameter values show more significant differences from the mathematically-defined values, although these variations are still within a tenth of a percent of the mathemtically-obtained value. Despite these deviations, three of the four software packages perform similarly. The definition of the Sdq parameter is based on the average gradient of the surface, which can be highly dependent on the scale at which the surface is measured [29]. This effect can also be affected by the resolution of the surface, and so it should be considered that the 700 × 700 dataset sampling resolution of the analytical surfaces is a primary reason for the deviations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results for autocorrelation parameter Sal are given in figure 4, for a selection of six analytical surfaces for which the autocorrelation parameters were calculable. In comparison with the previous results, the scale of deviations from the mathematically-obtained value is immediately apparent, with differences of approximately 10% common. This suggests the increased complexity and higher number of algorithmic steps required to calculate the parameter, relative to the simpler parameters, is contributing to the deviations in the values obtained by the software under test. It should be noted here that only software packages A and B provided the option to calculate autocorrelation parameters, and so results for only two of the four software packages under test are presented here.

Zoom In Zoom Out Reset image size

This novel method of software assessment by using mathematically-obtained reference values has demonstrated an ability to identify areas of the software under test that would benefit the most from improvement. By identifying which parameters deviate the most from a mathematical reference, in comparison to other software packages, resources can be allocated more efficiently, ensuring the weakest areas of the software are addressed first. The amount of information given to the software via a discrete dataset input can put a limit on the software's ability to obtain a parameter value close to the mathematical reference value. However, by performing an assessment at several resolutions of input surface dataset, software developers can identify to what extent the deviation from the mathematical parameter value is due to discretisation, and what is due to shortcomings in the software. This idea is addressed further in section 5.

To perform the extended precision calculations, the MATLAB plug-in Chebfun was used [30, 31]. Chebfun is an open source software package that utilises the idea that smooth functions can be effectively represented by expansions of Chebyshev polynomials (which are a series of recursively defined polynomials whose roots are used in polynomial interpolation techniques). Analytical functions are stored using the roots of Chebyshev polynomials, and are able to be represented accurately to machine precision, that is, fifteen to sixteen significant decimal digits. The use of Chebyshev polynomials to represent smooth non-periodic functions has been seen as an equivalent to Fourier series and its ability to represent smooth periodic functions and is a major component of the branch of mathematics known as approximation theory [32, 33]. The nature of the process is stable for even very large numbers of Chebyshev points, corresponding to complex analytical surfaces.

The aim of the Chebfun software is to achieve for functions what floating-point arithmetic achieves for numbers: rapid computation in which each successive operation is carried out exactly apart from a rounding error that is very small in relative terms [34]. This precise, numerical approach to computation with functions is an effective solution for overcoming the current computational limitations of purely symbolic calculations using CASs.

Before applying the Chebfun methods to new applications, such as complex surfaces and analytically challenging parameters, it is first worthwhile to perform a test to gain confidence in Chebfun's performance. This can be achieved by using the Chebfun method to calculate the same parameters as can be obtained symbolically using a CAS.

A selection of seven parameters were calculated for eight analytical surfaces mathematically, using the methods described in section 2.1. Each of the analytical results were then evaluated to fifteen significant decimal figures and stored numerically. Each analytical surface expression was then input into MATLAB and converted into a two-dimensional Chebfun object over the required evaluation area. Algorithms for the calculation of each surface texture parameter were performed, utilising the relevant Chebfun-specific functions. The results for each are given in tables 1 and 2 and show good agreement between the two methods. The Chebfun method is consistently able to match the mathematical result to within fourteen significant decimal figures, with the final fifteenth significant figure often only deviating by at most one digit.

With the Chebfun approach successfully able to obtain parameter values that agree with the mathematically obtained values to machine precision (minus rounding error in some cases), the next step is to extend this approach to more complex surfaces that better simulate realistic surfaces.

A selection of five surfaces were created with increased complexity, three pseudo-random surfaces with increasing high spatial frequency components, and two modified versions of simple surfaces from section 2.2 that have had pseudo-random elements added. Each of these surfaces were created utilising the Chebfun framework. For each, an analytical expression is converted into a Chebfun object and represented as a matrix of Chebyshev coefficients. This matrix is then manipulated, adding pseudo-random components to each element to increase the complexity of the surface and add realistic height variations. The Chebyshev coefficients in each matrix are then combined to form polynomial expressions using the recursive Chebyshev polynomial definitions

$$\begin{eqnarray}\begin{array}{rcl}{T}_{0}(x) & = & 1\\ {T}_{1}(x) & = & x\\ {T}_{n+1}(x) & = & 2{{xT}}_{n}(x)-{T}_{n-1}(x),\end{array}\end{eqnarray} \tag{ 19 }$$

and substituting them into the following equation

$$\begin{eqnarray}&&z(x,y)=\sum _{m=0}^{M}\sum _{n=0}^{N}{a}_{n,m}{T}_{n}(x){T}_{m}(y)\end{eqnarray} \tag{ 20 }$$

to create the final polynomial expression, where ${a}_{n,m}$ denotes the Chebyshev coefficient matrix elements, and N and M are the size of the Chebyshev coefficient matrix in the x and y directions, respectively. Surface topography images of the five surfaces created this method are shown in figure 5.

Zoom In Zoom Out Reset image size

The five complex surfaces were converted into Chebfun objects and then used to obtain surface texture parameters.

As the limitations of some symbolic calculations have been alleviated due to the numerical-based approach, parameters Sa and Sdr were also calculated in addition to the original set shown in section 2.2. Whilst the Chebfun method's calculation of the original set of parameters was verified against symbolic mathematical calculation in section 3.2, no such verification could be performed for Sa and Sdr. Therefore, although Chebfun may be able to calculate these parameters as accurately as the original set, caution must be taken when using the Chebfun parameter values as a reference, as errors may be present. For this reason, the calculation methods used should be explicitly stated alongside the obtained value to inform users of the potential sources of error. This is good practice that promotes transparency and should be carried out for any method of calculating reference values. In addition, alternative numerical calculation methods can be used alongside the Chebfun methods, and the number of significant figures that agree between both methods can be used as the reference value, increasing the confidence in the presented value.

With these reference parameter values obtained, it is now possible to use these values to assess the performance of third-party metrology software packages. The process was performed in a similar method to section 2.2; the five complex surfaces were sampled to create high resolution 700 × 700. SDF datasets which were then input into each third-party software package to calculate parameters. The results were then normalised to the mathematically-obtained values to enable easier comparisons.

Figure 6 shows the results for the Ssk parameter obtained by each of the software packages under test for all five complex surfaces. An interesting result to note here is the apparent similarities and differences between the underlying calculation methods employed by each of the software packages. Software A, C and D obtain similar results for each surface; even those for which they deviate from the reference value by as much as 3%. This suggests software A,C and D each calculate Ssk using a similar approach. Software B, however, clearly uses a different calculation method, allowing it to remain more consistently close to the reference value, but does not perform as well as the best cases of the other software packages.

Zoom In Zoom Out Reset image size

Figure 7 shows the results for the Sa parameter, newly calculable using the Chebfun method. All software packages under test show good agreement with the reference value, deviating by no more than one tenth of a percent. Here, software B shows the greatest amount of variation, suggesting the methods employed to calculate the Sa parameter are not as accurate as those employed by the other three software packages. Again, there is very close agreement between software A, C and D.

Zoom In Zoom Out Reset image size

Figure 8 shows the results for the Sdr parameter. Here, there are significant differences in the values obtained by software A and C, suggesting Sdr is one parameter that the two software packages calculate differently. Similar to the results for the Sdq parameter shown in figure 3, the Sdr parameter is affected by the resolution of the dataset, and the results here show that surfaces with more high spatial frequency components are the more significantly affected, with each of the software packages performing worse for those cases. It should be noted here that software D only gave Sdr results to two decimal places, making meaningful comparisons with the other software package results difficult.

Zoom In Zoom Out Reset image size

Through the use of the numerical methods made available using Chebfun, high precision reference parameter values can be obtained not only for parameters whose calculations are prohibitive for symbolic calculations using CASs, but also for surfaces with an increased degree of complexity that would otherwise be too computationally intensive. Using these new methods, a greater understanding of the strengths and weaknesses of the software under test can be performed using a wider range of surface texture parameters. In addition, the use of high complexity surfaces with high spatial frequency components delivers more information about the software under test, and how the calculated parameter values can vary.

As mentioned previously, Sal and Str are based on the autocorrelation function, and are calculated using the longest and shortest distances from the origin to a specified autocorrelation value. For the general case, difficulty can occur when applying the autocorrelation definition to the surface function, leading to an autocorrelation function of unknown shape which must be analysed to find the required decay lengths. This can be avoided by using a two-dimensional Gaussian function, whose autocorrelation function is also a Gaussian function. The Gaussian function about the origin is given by

$$\begin{eqnarray}&&{f}_{{ACF}}({t}_{x},{t}_{y})={e}^{-\left({{at}}_{x}^{2}+{{bt}}_{x}{t}_{y}+{{ct}}_{y}^{2}\right)},\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&a=\displaystyle \frac{{\cos }^{2}(\theta )}{{\sigma }_{x}^{2}}+\displaystyle \frac{{\sin }^{2}(\theta )}{{\sigma }_{y}^{2}},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&b=\displaystyle \frac{-\sin (2\theta )}{{\sigma }_{x}^{2}}+\displaystyle \frac{{\sin }^{2}(2\theta )}{{\sigma }_{y}^{2}},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&c=\displaystyle \frac{{\sin }^{2}(\theta )}{{\sigma }_{x}^{2}}+\displaystyle \frac{{\cos }^{2}(\theta )}{{\sigma }_{y}^{2}},\end{eqnarray} \tag{ 24 }$$

where σ_x and σ_y are the standard deviations, or widths, of the two-dimensional bell, and θ is the rotation of the Gaussian function about its centre. Using the default specified height of 0.2, the equation for the contour of interest becomes

$$\begin{eqnarray}&&{e}^{-\left({{at}}_{x}^{2}+{{bt}}_{x}{t}_{y}+{{ct}}_{y}^{2}\right)}=\displaystyle \frac{1}{5},\end{eqnarray} \tag{ 25 }$$

which can be rearranged to obtain

$$\begin{eqnarray}&&\left(\displaystyle \frac{{{at}}_{x}^{2}}{\mathrm{ln}5}+\displaystyle \frac{{{bt}}_{x}{t}_{y}}{\mathrm{ln}5}+\displaystyle \frac{{{ct}}_{y}^{2}}{\mathrm{ln}5}\right)=1,\end{eqnarray} \tag{ 26 }$$

which is the form of an equation for an ellipse, where the semi-major and semi-minor axes are given by ${\sigma }_{x}\sqrt{\mathrm{ln}5}$ and ${\sigma }_{y}\sqrt{\mathrm{ln}5}$ respectively. Therefore, the values of Sal and Str can be calculated directly using the semi-major and semi-minor axes from a surface defined in terms of σ_x, σ_y and θ.

Section 2.1.5 explains how the definition of Sa leads to a discontinuous surface function, causing problems when trying to evaluate a surface in the general case. Discontinuities occur after the application of the absolute function at the points where the surface crosses the zero-plane with a non-zero gradient. The absolute function causes the sign of the gradient to switch, leading to an abrupt change in surface slope. This can be avoided by defining a surface in such a way that the gradient at the zero-line is always zero.

One example of this is a modified version of a simple cosine wave, such that the sign of the wave is inverted for every second period. This example is given by an equation of the form

$$\begin{eqnarray}&&\begin{array}{l}z(x,y)\\ =\,\left\{\begin{array}{ll}A\left(1-\cos \left({f}_{1}\pi x\right)\right), & 0\leqslant x\mathrm{mod}\tfrac{4}{{f}_{1}}\lt \tfrac{2}{{f}_{1}},\\ -A\left(1-\cos \left({f}_{1}\pi x\right)\right), & \tfrac{2}{{f}_{1}}\leqslant x\mathrm{mod}\tfrac{4}{{f}_{1}}\lt \tfrac{4}{{f}_{1}},\end{array}\right.\end{array}\end{eqnarray} \tag{ 27 }$$

where A is the amplitude of the cosine wave and f₁ is the frequency. A profile example is given in figure 9. When the absolute function is applied to this surface, the parts of the surface below the zero-plane, which are the inverted sections of the surface, return above the zero-plane and recreate the unmodified cosine wave, for which Sa can easily be found.

Zoom In Zoom Out Reset image size

The Std parameter is defined as the texture direction of the scale-limited surface and gives the angle (with respect to a specified direction) of the most prominent texture on a surface. Mathematically, this is given as the absolute maximum value of the angular power spectrum of the surface for a given direction, s, given by

$$\begin{eqnarray}&&{f}_{{APS}}(s)={\int }_{{R}_{2}}^{{R}_{1}}r| F\left[r\sin \left(s-\theta \right),r\cos \left(s-\theta \right)\right]{| }^{2}{\rm{d}}r,\end{eqnarray} \tag{ 28 }$$

where R₁ and R₂ denote the range of integration in the radial direction, and also requires the Fourier transformation

$$\begin{eqnarray}&&F(u,v)={\iint }_{{A}_{{\rm{R}}}}z(x,y){e}^{-\left({iux}+{ivy}\right)}{\rm{d}}x{\rm{d}}y,\end{eqnarray} \tag{ 29 }$$

where u and v are spatial frequencies in the x and y directions, respectively. This calculation, involving two transformations and an absolute function (which has already been explained to be symbolically challenging for the case of Sa, see section 2.1.5), poses a complex symbolic challenge for the general case. Instead, it is possible to create reference analytical expressions with a known Std value that operate based on the parameter's more practical definition: the angle of the most prominent directional texture on a surface.

Creating a surface with a known, specific direction of texture can be achieved with a sinusoid with a known direction of propagation. By defining a surface of the form

$$\begin{eqnarray}&&z(x,y)=\cos \left(\left[A\pi x\right]+\left[B\pi y\right]\right),\end{eqnarray} \tag{ 30 }$$

the angle perpendicular to the direction of wave propagation, which gives the texture direction, is given by

$$\begin{eqnarray}&&{S}_{{td}}={\tan }^{-1}\left(B/A\right)+90^\circ \end{eqnarray} \tag{ 31 }$$

in degrees.

As mentioned previously, Sdr is calculated by finding the ratio of the evaluation area of the surface to the topographical area of the surface and subtracting 1. For the general case, the calculation of topographical area takes the form of a nonelementary integral. However, Sdr is readily calculable for more specific surfaces, a simple example of which is a linear slope, defined by

$$\begin{eqnarray}&&z(x,y)={Ax}+{By}.\end{eqnarray} \tag{ 32 }$$

A schematic of this is shown in figure 10. The topographical area of a linear gradient surface can be found using the Pythagorean theorem to find the lengths of the two perpendicular sides of the rectangle of the form

$$\begin{eqnarray}&&\sqrt{{x}^{2}+{A}^{2}{x}^{2}}=x\sqrt{1+{A}^{2}},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&\sqrt{{y}^{2}+{B}^{2}{y}^{2}}=y\sqrt{1+{B}^{2}},\end{eqnarray} \tag{ 34 }$$

and multiplying them together. The ratio of the topographical area to the evaluation area is then easily obtained, giving an expression for Sdr of

$$\begin{eqnarray}&&{S}_{{dr}}=\sqrt{1+{A}^{2}}\sqrt{1+{B}^{2}}-1.\end{eqnarray} \tag{ 35 }$$

Zoom In Zoom Out Reset image size

Similar to the assessments performed in sections 2.2 and 3.4, four surface texture parameter calculation software packages were used to showcase the validity of the parametric reference surfaces as a means to assess the performance of software. A total of five surfaces were created for parameters Sdr and Std, substituting different variable values into the parametric surface definitions. This approach demonstrates the flexibility of parametrically-defined surface/parameter reference pairs, as it enables a theoretically infinite variety of surfaces to be created that have mathematically traceable associated parameter values by adjusting the variables in the surface definition.

Figure 11 shows the results obtained by the software packages of the Sdr parameter for five surfaces obtained using a parametric surface function. Software C and D both show good agreement with the mathematical reference value for all surfaces, with software C obtaining normalised parameter values of the order of 1.000 000 05 for each surface. Software D shows more deviation than the other software packages for the fourth and fifth surface, which correspond to the surfaces with the lowest gradients. Conversely, software A and B performed poorly for the three surfaces that have the highest gradients and improve for the low gradient surfaces. In addition, software A and B obtain similar results for all five surfaces, suggesting a similar implementation of the Sdr calculation is performed for both software packages.

Zoom In Zoom Out Reset image size

Figure 12 shows the results obtained for the Std parameter, for five surfaces obtained using a parametric surface function. It should be noted here that only software packages A and C offered the ability to calculate the Std parameter. Initial results, given in the upper graph, show a significant difference between the performance of software A and software C in their ability to obtain parameter values comparable to the mathematical reference value. Further investigation revealed the values obtained by software C were all less than those obtained by software A by around 90°, suggesting a difference in interpretation of the definition of the Std parameter given in ISO 25 178-2 [7], specifically, relative to what axis should the most prominent angle of texture should be measured. Manually rotating the values obtained by software C gives the results shown in the lower graph, which show significantly better agreement with the mathematical reference values. With the corrections applied, software A is able to obtain a value for Std closer to the mathematical reference than software C for four of the five surfaces tested.

Zoom In Zoom Out Reset image size

In this section, a method for the assessment of surface texture parameter software performance is introduced. This approach quantifies the number of significant figures given by the software-obtained value that agree with the mathematical reference value. This number is then given alongside the resolution of the dataset used to obtain the software value, giving a performance metric describing of the quality the software's performance relative to the amount of discretisation error present. In addition, this performance metric allows for direct comparisons between third-party software.

Table 3 showcases the assessment of significant figure agreement for two analytically defined surfaces introduced in section 2. Here, the level of numerical precision that the software-obtained values agree with the reference values is presented for each parameter, allowing direct comparisons between software in a meaningful way that is directly applicable to the use of the software under test: to understand the extent to which the software gives an accurate result in comparison to a traceable reference. In addition, the total number of agreed significant figures is also calculated for each parameter set, providing a simple performance metric for the software under test for the reference surface as a whole. The results of the totals show software C gives the overall best agreement with the reference values for a dataset resolution of 700 × 700, narrowly improving upon software A for the '4th Bernoulli' analytical surface.

It should be noted that software D, as discussed in section 2, is limited to 10 MB dataset file sizes, restricting the maximum resolution of the tested datasets to 700 × 700. In addition, software D often gave parameter values to a reduced precision in comparison to the other test software, restricting its ability to perform well in this assessment.

A stage of the discretisation process that can affect the dataset and, therefore, the resulting parameter values is the sampling method used to obtain discrete height values from the analytical surface. The continuous surface can have different height values at very close x,y locations on the surface, and so small variations in the sampling locations can obtain different height values for use in the dataset. Additionally, averaging methods can be used that obtain single height values that describe areas of the analytical surface instead of specific sampling locations, which could also produce a different dataset for the same evaluation area.

There is no 'correct' sampling method, as different choices can be used to represent the height capture methods of different measurement processes. For example, point autofocus and contact stylus measurement instruments sample at discrete lateral intervals, whereas areal techniques, such as focus variation, obtain a full area image wherein each pixel represents a small area of the surface. In any case, sampling method variations can be addressed by detailing the method used when performing the software assessment. By explicitly defining the sampling method used, any further tests can be made using the same sampling methods and inter-software comparisons can still be made.

To highlight this variation, a second sampling implementation was used to compare to that used in all previous sections. For the previous sections, each height value data point was treated as a pixel with a width equal to the spacing between x, y sampling points. Therefore, to span the appropriate area, one sampling point was removed to account for the extra half spaces either side of the first and last data points in each direction. For example, for a 1000 × 1000 resolution surface spanning 10 mm by 10 mm, height values were sampled between 0 mm and 9.99 mm in each direction, accounting for the pixel widths which extend the surface size to between −0.005 mm and 9.995 mm. The second sampling implementation uses the same spacing between points, but includes an extra point on the end, giving a 1001 × 1001 resolution dataset with height values sampled between 0 mm and 10 mm. This method treats the dataset as a collection of discretely sampled height values with zero pixel width, and so still spans the same evaluation area as the previous implementation.

Table 4 gives the same significant figure analysis as performed in table 3, but for the second sampling method. Software B shows a significant improvement in its agreement with the reference values, while software A, C and D all show worse results. These performance metrics suggest this style of sampling is well interpreted by software B, whereas the other software packages are more suitable for the previous implementation. Again, neither implementation is 'correct', and clear specification of the sampling choices used when creating a dataset are necessary for proper assessment of metrology software.