Towards a new definition of areal surface texture parameters on freeform surface

Article history: Received 21 November 2016 Received in revised form 19 April 2017 Accepted 10 May 2017 Available online 17 May 2017


Introduction
Nowadays advanced manufacturing technologies empower complicated geometries to be designed and manufactured. For example additive manufacturing (AM) allows internal and external components with freeform geometries to be manufactured. This has led to a series of challenges in characterisation of complex surfaces. One challenge is the measurement of complex surfaces and the follow-on is the characterisation of the measured surface. While computed tomography (CT) has been considered as one of the promising measurement systems for certain complex surfaces (subject to its material), extracting surface information from CT measurement is an imperative operation.
To perform such extraction, the first necessary step is the reconstruction from volume to mesh that can produce an adaptive meshing based on the maximum allowable distance between the implicit function (surface with constant grey value) and the final triangular mesh [4]. A bonus of such operation is that the extracted mesh allows re-entrant features (undercuts). However one cannot directly compute areal surface texture parameters from such extracted mesh.
Areal parameters ('' Birmingham 14" project) where initially developed during a EU funded project leaded by Birmingham university [22]. There were some ambiguity in the definition of the parameters regarding the features characterisation and the material ratio parameters needed to be more indicative [14]. Further developments and standardisations, such as a detailed analysis of advanced and robust filtration techniques and subdivision of surface in hills and dales through the watershed decomposition, were proposed in the subsequent SurfStand project and published in Blunt and Jiang [3]. These definitions were later adopted in ISO 25178-2 [13]. The assumption on the calculation of the areal parameters is that the measured surface is represented by a planar mesh; it is assumed that the height coordinate of the measured surface can be described as a function of the others, z ¼ f ðx; yÞ. Current commercial measurement software packages compute height and hybrid parameters based on a discrete approximation on a rectangular grid. With irregular 3D meshes (such as the extracted surface from CT measurement) is not possible to perform the integrals based on the definitions from ISO 25178-2 otherwise the computation results of the parameters will have significant bias. It is then necessary to generalise those parameters to be suitable for both regular and irregular meshes. A study aiming to define the height parameters on a triangular mesh has been recently proposed in Abdul-Rahman et al. [1]. The authors proposed to firstly parameterise the surface form, represented through a triangular mesh, and then compute the areal parameters in the parametric space. The generalisation of the surface texture parameters is not restricted to triangular meshes, but can be applied to all surfaces. The extension to triangular meshes is performed through the finite element method and does not depend on the mesh parametrisation, whereby the proposed triangular mesh method is faster and no distortion is introduced.
In this paper a set of new definition of areal surface texture parameters (height parameters, i.e. Sa, Sq, Ssk and Sku, and hybrid parameters, i.e. Sdq and Sdr) for general parametric surface is proposed. One the advantages of this method is that it is possible to characterise a freeform surface without unwrapping it. It should be noted that if the form surface does not have a parametric representation it is not always possible to unwrap the mesh.
All computation of the parameters on the examples in this paper are carried out on the primary surface without any filtration, to eliminate any effect of filtration that bias the comparison of the parameters' result.
The paper is organised as follows: the proposed method is presented and discussed in Section 2 including two methods for the surface approximation. In Section 3 the proposed method is tested on six sets of different types of surfaces and the computed parameter results compared with the values of the parameters defined in the ISO 25178-2.

Parameter definitions
In this paper it is assumed that a manufactured surface can be described with a regular surface R & R 3 as rðu; vÞ ¼ with u ¼ ðu; vÞ T and U & R 2 ; U is called the parameters space.
Suppose that it is possible to decompose the surface into two parts rðu; vÞ ¼ r form ðu; vÞ þ r res ðu; vÞ ð 2Þ where R form : r form ðu; vÞ represents the form and R res : r res ðu; vÞ the residual surface [19].
The proposed method can be explained with a profile example shown in Fig. 1a. A portion of a profile is described by a parametric curve x; y ð Þ T ¼ rðtÞ. It is assumed that the form profile is y ¼ 0 and 20 points uniformly distributed on the curve are generated. If the curve is treated as a function y ¼ f ðxÞ and if the integral is computed with the adaptive rectangular rule, the computation of Ra can be carried out as illustrated in Fig. 1b. It is possible to observe that the y values oscillate. The proposed method is based on the computation of Ra of a parametric curve as be the area highlighted in Fig. 1c, where t is an abstract parameter. It can be, for example, arc length parametrisation [6]. The residual coincide with y because the form surface is y ¼ 0; the profiles in Figs. 1a and c represent the same profile, the only difference is the representation used for the abscissa. Although the discrete method is not conceptually correct, as the number of sampled points increases, it can approximate the value of the parameter. In fact the discrete method can be interpreted as a Monte Carlo integration of the function y ¼ yðtÞ, but the stretching of the form profile in the parameter space t is not considered. Here it is proposed to extend the roughness parameters based on the height of the surface, i.e. it is not related to the volume of the material (void if the value of the function is negative) between the scale limited surface and the form surface.
If a total least squares approach is implemented the last term in Eq. (2) could be rewritten as r res ðu; vÞ ¼ r res ðu; vÞn form ðu; vÞ ð 3Þ where r res ðu; vÞ is the distance between rðu; vÞ and its projection on the form surface r form ðu; vÞ and n form ðu; vÞ is the surface normal. r res ðu; vÞ can be interpreted as a scalar field on the surface r form ðu; vÞ. If it is possible to describe the value of the surface r res ðu; vÞ on the form surface without stretching and if no reentrant features appear on the residual surface, the parameters can be computed with the definition of the ISO 25178-2. Form surfaces developable to a plane are all the surfaces where the Gaussian curvature is null everywhere [6], such as cylinder, so it is possible to transform this surface to a plane without stretching. When the form surface cannot be described by a developable surface, such as a general freeform surface, the local stretching must be taken into account [10]. Let r form;i the partial derivative along the dimension i, the parameters on the primary surface can be computed by weighting the function values with the infinitesimal surface area dr form ¼ r form;u ðu; vÞ Â r form;v ðu; vÞ du dv: The scenario in Eq. (3) is analysed first, after that a method to estimate the areal parameters of a general surface is presented. According to the definition of the parameters in the ISO 25178-2, a generalisation of the arithmetic mean of the absolute value of the height can be computed as [19] Sa where A ¼ R R R form dr form is the area of the form surface, the root mean square height as the root mean square gradient as where J form ¼ J form ðu; vÞ is the Jacobian matrix of the form surface, form J form and $ u r res ðu; vÞ is the gradient of the residual surface.
The developed interfacial area ratio can be computed as where r u ðu; vÞ ¼ r form;u þ r res;u n form þ r res n form;u r v ðu; vÞ ¼ r form;v þ r res;v n form þ r res n form;v where n form ðu; vÞ is the normal vector of the form surface and E; F and G are the coefficients of the first fundamental form and e; f and g are the coefficients of the second fundamental form [6] If the vector function r res is available the parameters could be computed similarly to the flux integral [10]. The absolute value of the height can be computed as the root mean square gradient as If the unit normal of the form surface is not differentiable, an approximation of the Sdq parameter could be where J res ¼ J res ðu; vÞ is the jacobian matrix of the residual surface.
The developed interfacial area ratio can be computed as where r u ðu; vÞ ¼ r form;u ðu; vÞ þ r res;u ðu; vÞ r u ðu; vÞ ¼ r form;v ðu; vÞ þ r res;v ðu; vÞ: In order to compute the areal parameters the surface must be reconstructed first. In this work two reconstruction methods are analysed: a parametric surface represented by the locally refined B-spline and the triangular mesh. B-spline surface reconstruction has been used in Harris et al. [12] to approximate discrete data and to compare various software packages, while the triangular mesh is the common representation in finite element analysis. In the following sections the two methods are briefly introduced.

Surface reconstruction -LR B-spline
Since the reconstructed surface can be described by a triangular mesh, in order to transform it to a parametric function it has to be parameterised. Parametrisation of a mesh means attaching a coordinate system [5]. The aim of the parametrisation is to describe the points of the triangulated surface (x i ) with a function x ¼ f ðu; vÞ, where u and v are two abstract parameters. In the previous section, it has been assumed that the parameters space is known, but when with real data is used the parameter space must be constructed. Stretch minimising approach proposed in Yoshizawa et al. [25,24] is employed because it minimise the area distortion, it is therefore a good candidate for the surface reconstruction. The method has been implemented in OpenGI [20]. This parametrisation are used for all the surfaces illustrated in this paper.
The form surface is then estimated with a total least squares (TLS) approach. Plane, cylinder, sphere and a complex surfaces will be analysed to show that the proposed method does not depend on a specific form surface. The differences between the point cloud and the projections on the form surface are computed first. Both the two point clouds, form and residuals, are then approximated with the LR B-spline algorithm [9]. It is not possible to use the same algorithm used in Harris et al. [12] because the interpolation is possible only with samples available on a regular grid in the parameter space. The locally refined algorithm iteratively partitions the parameters domain adaptively, i.e. where the difference between the reconstructed and the measured surface is greater than a specified tolerance. The algorithm terminates when all the differences between the approximated and the measured points is below a threshold or when a maximum level of the refinements is reached. A maximum value of 10 iterations and a threshold value of 0.01 were set in the approximation algorithm. During the approximation stage all the values were coded between 0 and 1 to avoid the scale effect. The abstract parametrisation domain was 0; 1 ½ 2 . The areal parameters (Sa, Sq, Ssk, Sku, Std and Sdr) of the primary surface can be computed as the integral of the residuals on the form surface. Since both the involved surfaces, r form and r res , share the parameters domain, the integration can be performed with a numerical quadrature rule. The numerical integration is performed with the h-cubature method implemented in Johnson [15]. This method recursively partitions the integration domain into smaller sub-domains, the quadrature rule is applied to each, until a convergence criterion is reached.

Surface reconstruction -Triangular mesh
If the surface is approximated with a triangular mesh, the surface integral of a function f ðb 0 ; b 1 ; b 3 Þ (in barycentric coordinates system) can be computed as where n tri is the number of the triangles in the mesh, A i is the area of the i-th triangle of the form surface, n q is the number of the quadrature points, f ij is the value of the function in the i-th triangle at the jth quadrature point and w j is the quadrature weight. The area of each triangle can be computed using the definition of the cross product where v i is the i-th vertex of the triangle, while surface area is the sum of the triangles areas Linear elements are used to approximate the integrals. The parameters are computed as a per vertex scalar field (see Fig. 2). Sa, Sq, Ssk and Sku can be computed using Eq. (18). The quadrature points depends on the degree of the function to integrate. The efficient symmetrical Gaussian quadrature rule implemented in Shao [8] is used. It allows the computation of the minimum number of quadrature points to evaluate the integral on a triangular domain. The gradient of the residual function on the l-th triangle can be computed as [5] $r res ðbÞ ¼ ðr res ðe j Þ À r res ðe i ÞÞ T is a vector of barycentric coordinates, e i is the vector with one in the i-th entry and zeros otherwise and ? denotes a counterclockwise rotation by p 2 in the triangle plane. The mean square gradient can therefore be computed as where $r res ðiÞ is the gradient on the i-th triangle. There is only one quadrature point because the gradient is a constant function on each element. The value of the parameter Sdr can be approximated evaluating both the area of the measured or scale limited mesh and the area of the form surface.

Experimental results
In this section the proposed method is evaluated both on measured and simulated surfaces. The computation of parameters with the proposed method is first compared with the ISO method implemented in SurfStand v6.0 [23]. Four validation surfaces available in SurfStand are first analysed. The heights of the surfaces can be described as a function of two coordinates (z ¼ f ðx; yÞ).
The rest of this section focuses on the evaluation result of two measured surfaces: a casting sample (Rubert & Co. Ltd, UK) with nominal Ra of 25 lm and an AM component. The casting sample was measured both with a Focus Variation (FV) instrument (Alicona InfiniteFocus G4g [7]) and with a CT device (Nikon XT H 225 microfocus CT). The FV dataset presents an equally spaced sampling, the resulting mesh is isotropic (all the elements have similar area). It is used to evaluate the agreement of the proposed method with the ISO 25178-2 definition. Two simulated data sets are then constructed adding the residuals to a developable (cylinder) and a non developable (sphere) surfaces. The effect of the distortion during the unwrapping stage is discussed. CT sets of data for the casting sample and the AM component are analysed to investigate the effect of an unequally spacing of the points on the computation of the parameters. Finally, in order to prove that the proposed method does not depend on the form surface, a simulated complex surface is investigated. The following methods are analysed: SurfStand: it represents the parameters computed with the SurfStand software. LR B-spline SL: it represents the parameters computed with the scale limited (primary) surface, according to the ISO 25178-2 definition, approximating the surface with the LR B-spline algorithm. Mesh: it represents the parameters computed with the proposed method using the triangular mesh approximation. LR B-spline 1D: it represents the parameters computed with the proposed method based on the scalar field (Eqs. (5)-(10)), approximating the surfaces with the LR B-spline algorithm. LR B-spline 3D: it represents the parameters computed with the proposed method based on the vector field (Eqs. (11)-(17)), approximating the surfaces with the LR B-spline algorithm. Discrete: it represents the approximation of the surface texture parameters with a meshfree method.

Validation surfaces
In this section four surfaces are analysed to evaluate the performances of the proposed methods. As shown in Fig. 3, each surface has its own characteristics. The first surface presents two grooves along a direction, the second one has some regularly located peaks, the third one has some plateaus and the forth one is a typical random surface.
Each surface is levelled with a least squares plane in SurfStand and the primary surface analysed. Tables 1-4 show the computed values of the parameters. The discrepancies between the height parameters (Sa, Sq, Ssk and Sku) are negligible. The maximum differences is the estimation of Sq (2%) for the skin_1 surface, and is probably due to the limited amount of points of the surface compared to the other test cases. The method based on the triangular mesh tend to overestimate the value of the Sdq parameter. Sdr estimations agree between all the tested methods except for the stain_steel surface. The biggest difference is equal to 9.22% between SurfStand (28.15%) and LR B-spline 1D (18.93%). The difference is due to an approximation of the computation method implemented in SurfStand [2]. It should be noted that values based on the LR B-spline algorithm have slightly lower Sdr values, this is an effect of the smoothing during the surface approximation stage.

Rubert casting sample: FV reconstruction
A Casting sample with nominal Ra of 25 lm was measured with the Alicona InfiniteFocus G4g. The form surface was approximated with the total least square plane, that corresponds to the first two scores of the principal component analysis (PCA) of the point cloud covariance matrix. Fig. 4 shows a portion of the measured mesh of the sample (grey) and the estimated form surface (blue). Since no undercuts were present in the FV mesh, it is also possible to approximate the scores of the PCA with the LR B-spline method. In Table 5 the parameters computed according to the ISO 25178-2 and with the proposed methods are reported. Since the scores of the PCA are not located on a regular rectangular grid, it was not possible to use SurfStand. The discrete method was used to compute Sa, Sq, Ssk and Sku. There is no difference between the discrete, LR B-spline SL and mesh methods in the estimation of the parameters. The discrepancies of the proposed methods are bigger due to the surface approximation, but since the maximum percentage error is 0.6% the differences are negligible.
The above procedure was applied to the whole set of data of the measured casting sample (see Fig. 5). The evaluation of the parameters was applied to a bigger point cloud to check the robustness of the approximation. Table 6 shows the parameters computed on the Fig. 3. Demo surfaces (values in lm).
primary surface. As for the previous scenario the parameters estimated with all the analysed models obtain similar results.

Primary surface on cylinder
To check the stability of the proposed method a data set with a cylindrical form was analysed. The scores of the PCA of the previ-ous test case were approximated with the MBA algorithm [16] to predict the points on a regular grid. The LR B-spline algorithm was not applied in order to investigate the robustness of the reconstruction method. The approximate points on a regular grid of 985 Â 799 were then added, along the normal direction, to a portion of a cylinder. The angle of the cylinder ranges from Àp to 0, while the radius was computed as     q ¼

Dy Dh
where Dy is the resolution of the coordinate in radial direction and Dh is the angle resolution. With this radius the distances on the cylinder coincide with the distances on the plane. Fig. 6 shows the simulated mesh and the computed mesh parametrisation is shown in Fig. 7. The primary dataset is split between the y and z coordinates. It is not possible to see clusters of points, illustrative of the goodness of the parametrisation. The surfaces were reconstructed with the method described in the previous section. The computed parameters are reported in Table 7. From the table it is possible to observe that the proposed method can achieve good performance also when the form surface is a developable surface.
In order to check the robustness of the procedure when the nominal form is not a plane, the full FV point cloud was added to the half cylinder shape. Table 8 shows the computed height parameters. The differences are similar to the values of Table 7, illustra-tive of the robustness of the proposed method respect to the mesh sample size.

Primary surface on a sphere
In this section the data set of the previous section was added to a portion of a sphere to evaluate the performances on a non developable surface. Fig. 8 shows the simulate measured surface and the unwrapped one. The unwrapped coordinates are computed as where q, set to 1, is the radius, h is the longitude and / is the colatitude. Since it is not possible to unwrap a sphere without distortions the length on a curve parallel to the equator, especially near the pole, has a strong stretching. The estimated parameters are shown in Table 9. It is possible to observe that, since the residual are non structured on the surface, the differences between the parameters is negligible. To understand the effect of the stretching the values of the areas of the form and the measured surfaces are reported in Table 10. The difference of the ISO method, compared to the triangular mesh, is 50.99% and 48.09%, respectively, for the form and the measured surface. It is therefore possible to conclude that the small difference in the parameter Sdr is due to a compensation of errors.

Rubert casting sample: CT reconstruction
The casting sample with a nominal Ra of 25 lm was acquired also using a Nikon XT H 225 microfocus CT. Nikon CT-Pro software [18] was used to perform the volume reconstruction. CT voxel size for all coordinates was 12.9 lm (x; y; z). Mesh reconstruction was performed by an adaptive algorithm implemented in CGAL [4,21]. This algorithm allows the reconstruction of an implicit surface with the desired approximation; the maximum allowable   error was set to 5 lm, almost 1 3 of the voxel size. The output mesh is a manifold mesh, so no post processing is needed in order to compute the parametrisation. It should be noted that the mesh reconstructed with the marching cube algorithm [17] may have some non manifold vertices or edges. Fig. 9 shows the reconstructed mesh (grey) and the estimated form surface (blue). Two datasets are analysed, the small one is a subset of the other. These sets of data correspond to the meshes analysed in Section 3.2.
Considering the adaptive meshing applied the spacing of the points is not constant, this can lead to a biased estimation of the height parameters with the discrete method. Since the surface may present some undercuts, it is not possible to approximate the surface of the scores of the PCA and the Sdq parameters of Table 7 Areal parameters FV data subset on cylindrical shape.

Sa (lm)
Sq (   the 3D residual must be computed using Eq. (16). It should also be noted that the discrete method is not compliant with the ISO 25178-2 because it is not an approximation of integrals described in the standard. The values will be computed to evaluate the bias. After applying the parametrisation and the reconstruction the areal parameters were computed. Tables 11 and 12 show the parameters of the datasets. The difference between discrete and the other methods increase compared to the previous test case. This is the effect of the unevenly spaced points and the re-entrant features.

AM component: CT reconstruction
An additive manufactured component measured and the surface was reconstructed with the algorithm mentioned above. The CT voxel resolution was 17.5 lm in all x; y and z directions. The threshold was selected according to the ISO 50 method implemented in VGStudio Max software [11]. The analysed surface along with a magnification, where is possible to observe a reconstructed undercuts, these are shown in Fig. 10. The reconstructed surface was again parametrised and approximated with the LR B-spline approximation.
The computed parameters are reported in Table 13. Sa and Sq parameters computed with the discrete approximation are biased, while the differences between Ssk and Sku are negligible. The absolute value of the difference is comparable to the previous test case. But, since the estimation are smaller, the percentage differences correspond to 8.20% and 6.73% for Sa and Sq. These discrepancies should be taken into account because the values are computed on the same set of data.

Complex surface
In this section a complex form surface is tested. The aim of this section is to show that with the proposed method it is possible to analyse the surface texture on a freeform surface. The surface of a turbocharger was reconstructed with the CT device (the voxel size was 57.23 lm for x; y and z). A maximum allowable distance of 5 lm was set and the reconstructed surface is shown in Fig. 11a. A blade was manually segmented after the mesh reconstruction (blue mesh in Fig. 11a and mesh in Fig. 11b). After the parametrisation of the form surface the primary dataset of the casting sample data from CT reconstruction has been added to the form  surface. Fig. 11c and d shows the residual and the simulated meshes. During the parametrisation there are some distortions both on the left and right corners. Table 14 reports the computed parameters. In this test case it is possible to observe that the differences between the parameters are negligible. The goodness of the discrete method can be explained with the low value of the Sdr parameter and the isotropy of the reconstructed mesh. It should be noted that there is no other method to compute the Sdq and Sdr parameters.

Conclusion
A method to compute the height and hybrid parameters of areal surface texture for freeform surfaces has been proposed based on two types of surface reconstruction: LR B-spline and triangular mesh. The proposed approach has been compared first with the definition of the areal parameters described in the ISO 25178-2. The six test sets have been carried out on reconstructed irregular meshes from plane, cylinder, sphere and a complex surface. The test results show a good approximation to the ISO definitions. It also indicates that the Sa and Sq parameters computed with the discrete method are biased, while the difference between LR Bspline and triangular mesh is negligible.
Future work involves development of the generalisation of the definitions of other areal parameters if the analysed surface is complex or presents undercuts. The concept of scale-limited surface for complex surface has also to be defined and investigated. All the computed values in this paper are based on the primary surface, this is due to the fact that the S filter, L filter and F operator defined in the ISO standards are based on a planar surface. Therefore revising the set of filters based on irregular meshes is also part of the future work.