Gaussian process regression‐based detection and correction of disturbances in surface topography measurements

Modern smart and intelligent manufacturing is characterised by an increasing use of highly engineered surfaces and quasi‐free form geometries, for example, by additive manufacturing, and the requirement for fast and informative measurement tools for quality controls in production. These have lately pulled towards the adoption of optical surface topography measuring instruments to qualify technological surfaces, which are core to being assessed to characterise the product and optimise the manufacturing process. Surface topography measurements performed by optical instruments are known to suffer from disturbances, such as non‐measured points and spikes. These are due to complex interactions between the measurand and the instrument and may result in poor measurement quality and biased characterisation results. Currently, the correction of measurement disturbances is carried out by empirical approaches which mostly involve thresholding and interpolations, whose sensitivity is often devolved to the operator. In this work, a Gaussian process regression‐based approach is outlined to identify and correct measurement disturbances exploiting spatial correlation properties of the measurements to provide a formal, robust, and univocally defined approach to manage measurement disturbances. The proposed approach, differently from currently available alternatives, is capable of managing at once the different type of disturbances by means of a supervised machine learning technique. The formal and practical advantages of the proposed method are discussed exploiting case studies of industrial interest applied on quasi‐flat surfaces to stress the disturbances effect.


F I G U R E 1 (A) surface topography with evident and isolated spikes. (B) topography with some voids or non-measured points INTRODUCTION
To satisfy the continuous ambitious demands for enhanced performances, customisation and more sustainable use of resources, manufacturing is facing the challenging development of novel processes, for example, additive manufacturing, 1 and advanced materials, for example, innovative composites 2 and coatings. 3 The evolutionary context of modern manufacturing, enlivened by big data, intelligent and interconnected cyber-physical systems, requires flexible and fast quality inspections that relies on thorough, accurate and precise characterisation methods. [4][5][6][7][8] The measurement and characterisation of technological surfaces, that is, the extraction of some synthetic indexes, are core to describe a topography quantitatively. These find application in quality and statistical process control to set up, optimise and control manufacturing processes, 9,10 which are necessary to provide and control proper functionalisation of the surface of a component. [11][12][13] As the relationship among technological surface topography, manufacturing process, via the concept of manufacturing signature, [14][15][16] and surface properties 17,18 has been established well, academia and industry have been increasingly focussing on design, 19,20 modification, employing surface technology, 18,21 and characterisation of technological surfaces. [22][23][24] Advanced surfaces, which accommodate aggregated functions, necessitate of complex structures often layered over different length scales. For instance, the enhanced possibilities in surface generation have made very attractive emulating the behaviour of living creatures in nature. 25,26 Such bio-inspired surfaces are intrinsically free form and can exhibit diffraction and/or antireflection properties, 27 hydrophobicity and self-cleansing properties, 28 plasmonic colouring, 29 friction control for grip improvement of robots and cobots. 30 Freeform structuring of surfaces requires reliable and versatile surface metrology capabilities to characterise the dimensional and geometrical features that achieve the surfaces' functions. Surface topography measuring instruments are intrinsically 'areal', as they map the height z(x,y) of a surface topography with respect to a reference measurement surface. In particular, surface topography measuring instruments acquire a data set of pixels mapping a surface that can be represented mathematically as a function F in the Euclidean space defined on an open subset of the Euclidean plane F: ℝ 2 → ℝ 3 . They can offer several advantages, ranging from better resolution, flexibility and allowing for non-contact operations, which is of great importance when dealing with soft surfaces. 31 Nonetheless, noise, and disturbances in general, affect these instruments, typically with negligible uncertainty contribution on the horizontal axes A recent international comparison of optical instruments, evaluating areal roughness parameters, highlighted how the presence of noise developed into an additional contribution to the measurement uncertainty. 32 The measurement noise can be increased by a non-optimal choice of the instrument settings (e.g., the light intensity), but also by spurious interactions between the instrument's light and the material of the component under measurement. Thus, defects or irregularities on the surface may trigger high-light scattering, and anomalous changes of the back-reflected light may be enhanced by the structuring itself of a component's surface. These complex interactions lead to two major issues: the absence of measured pixels, that is, pixels characterised by the absence of the corresponding recorded height values in the micrographs-called voids or non-measured points-and anomalous recorded values, which are significantly discordant with respect to the neighbouring pixels-usually referred to as spikes. 31,33,34 Figure 1 shows two examples of the mentioned disturbances in surface topographies that will be considered further on in the paper.
Both spikes and voids counteract micrographs evaluation, preventing the correct detection of specific surface features and significantly changing the assessment of roughness parameters. [35][36][37] Consequently, such discordant observations cannot be included in the evaluation of the measured set of pixels mapping a portion of a surface. Therefore, proper identification and correction, that is, the management, of the measurement disturbances are necessary to guarantee the correctness and unbiasedness of the characterisation for reliable application in quality and statistical process control. Surprisingly, although acknowledging the issue, the literature is scarce of solutions to address the management of these measurement disturbances. 32,35,36,38 As far as non-measured points (NMs) are concerned, these are inherently identified by their definition and typically corrected by filling, that is, imputation, strategies. These can either exploit a smooth interpolation or the selection of a particular value, for example, the average, the maximum, the minimum, of neighbouring points height. 36,39 Because the latter option is clearly arbitrary, it is little exploited. 36 Conversely, more alternatives are available to manage spikes. These have been recently reviewed by Podulka et al. 40 and can be split into two groups: (i) thresholding methods and (ii) filtering methods. The former identifies as spikes all the measured values beyond some certain upper and lower limits UL and LL respectively, selected by the operator. Thus, the set of spikes results as = { ( , ) ∀ ( , ) ∈ ∶ ( , )⟨ ∨ ( , )⟩ }. Giusca et al. 38 suggested as thresholding criterion = − = 3 (Sq is the root mean square error of the measured z(x,y), see Section 2.2 for more details). This choice, in line with other empirical quality rule of thumb, 41 has the advantage of being univocally defined, and hence is independent from the operator experience. However, when addressing complex surfaces, such as freeform and, thus, inherently 3D, this criterion is inadequate as it may identify as spikes topographical features actually belonging to the topography. Replicated measurements can be exploited, to remove the form and reduce the problem back to simple flat surfaces. In this case, the thresholding is applied on the micrograph resulting from the subtraction between the i-th measurement and the average micrograph of the replicated measurements. 38 Relying on replicated measurements, alternative thresholding limits can be evaluated pixel-wise by computing 3-sigma confidence intervals. 40 Although effective, replicated measurements may not be practical for current industrial application involving in-line quality control because complex surfaces require long measurement time, which is not compliant with the request of fast and informative inspection methods.
Differently, filtering methods identify spikes as anomalous values, that is, outliers, of the residuals and operate the correction by substituting the spike with the filtered value. Several filters are available, Gaussian, splines, opening and closing filters, with their robust counterpart [42][43][44][45] and median filter. 39 The literature points out that filtering may be less severe than thresholding. 35,40 However, filtering is a non-trivial operation which may result in a wrongful identification of some topographical scales as spikes, or denoising, 40,46 for example, by erroneously dimensioning the nesting index.
Therefore, although several methods are available in the literature for managing measurement disturbances, they lack generality, and, in some cases, they are not suitable for current industrial challenges and feature, as discussed, some criticalities. Moreover, even though the literature demonstrated that measurement noise is strongly spatially correlated, 20,47-50 none of the methods mentioned above takes into account the correlation neither in the detection nor in the correction phase. Furthermore, the literature alternatives consist of separated approaches to manage the different types of disturbances. This might expose the practitioners to correcting first the non-measured points, which though would result in biased corrections if any spike is present in the neighbourhood considered for the interpolation.
Therefore, in this work, we propose a methodology to unify the management of spikes and non-measured points, which is also capable of considering the spatial correlation of the measured topography. This paper proposes a supervised learning model, based on Gaussian Process Regression (GPR), for locating and correcting voids and spikes. The model is successively compared with a commercial software, 39 and, for what concerns spikes, also with a hands-on 'thresholding' method proposed in the literature. 51 Specifically, the paper is organised as follows: Section 2 discusses the GPR-based approach formulated for detecting and correcting the disturbances. Section 3 presents some case studies to test the presented methodology with respect to approaches available in the literature. Section 4 discusses the obtained results relevant to the case studies and Section 5 eventually draws conclusions.

GAUSSIAN PROCESS REGRESSION MODEL
GPR, also known in geostatistics as kriging and in computer experiments as surrogates or emulators, 52 is a stochastic interpolation technique apt to predict values from a sparse and noisy dataset, provided that the model responses are correlated. It was first conceived for geostatistics application by D. G. Krige in the early 1950s as a heuristic approach. Later on Matheron provided a theoretical structure 53 that was modified further for the application to the field of computer experiments by Sacks et al. 54,55 In the years GPR has been applied in manufacturing field to predict manufacturing macro-geometrical 56,57 and micro-geometrical errors, 58 assembling errors for compliant 59 and composites 60 materials, to support in the definition of adaptive and sequential sparse and economic of inspection plans. 12,61,62 More recently, GPR modelling has been classified as a machine learning technique, successfully applied for reliability and life prediction of cutting tools 63,64 and to support multi sensor data fusion. 65,66 Although other regression models based on spatial correlations are available to predict data arranged according to a structured grid, for example, CAR, SAR, Markov Random Field, GPR is preferred here. In fact, other solutions hypothesise a discrete domain. Conversely, the measurand in the case at hand, that is, a surface topography, underlies a latent continuous spatial domain, although discretised by the measurement instrument. This feature supports the adoption of GPR to model the measured surface based on a set of finite measured data. 53,[67][68][69] Let D ⊂ ℝ 2 a domain where the response z has been experimentally observed in a set of points = { 1 , … , , … , }, with = ( 1 , ) ∈ , ∀ ∈ {1, 2, … , }, that is, D is the definition area A, that is, the measured field of view, so that the set of tried points is = ( 1 , … , ) . The Gaussian process model assumes the empirical response z to be a realization of the process Z(x): where = ( 1 , … , ) contains the set of coefficients to linearly combinate some specified trend functions in ( ) = ( 1 ( ), … , ( )) . Ψ( ) ∼ (0, 2 (ℎ; )) models the regression error as a Gaussian random process with zero mean and stationary covariance, [Ψ( ), Ψ( )] = 2 (ℎ; ), with , ∈ {1, 2, … , }. 2 is the process variance and R is the correlation matrix, whose entries = (ℎ; ) are the spatial correlation function and depend only on the displacement vector h between any pair of points in and on a vector parameter . If the value (ℎ; ) depends only on the length ‖ℎ‖ 2 of the vector ℎ, and not on the direction along which it is computed, then the stochastic process is isotropic. From the former discussion follows that: GPR aims to predict the model response at an untried location 0 , so that 0 = ( 0 ). It is demonstrated that by min- , the predictor of 0 , that is,0, is: where 0 = ( ( 0 − 1 ), … , ( 0 − )) and F is the matrix with entries = ( ), ∈ {1, 2, … , }, ∈ {1, 2, … , }. 70 The formulation in Equation (3) is the most general one and is often referred to as universal kriging. 53 It is worth pointing out that the prediction is conditioned to the tried points, that is,0 = [ 0 | ]; moreover, it is the best linear unbiased predictor (BLUP), and it is unique. 53,62 The prediction can also be interpreted as made up of two contributions: a regressive term, ( 0 ) , which depends on (or equivalently its estimate), and a correction, 0 −1 ( − ). The meaning of the latter is particularly interesting. In fact, it is a linear combination of the regression residuals weighted for the spatial correlation function value. Specifically, the correction term controls the efficiency of the prediction and allows obtaining acceptable predictions even if the regression model is not the best one, that is, even if the problem is ill-stated or ill-conditioned, 71 to the point that, often, the regression function is chosen as a constant, that is, F = 1. 54,55 This choice is referred to as ordinary kriging. 53 Therefore, the prediction is highly dependent on the spatial correlation of the process, which motivates why it is necessary to adequately model the spatial correlation function, or kernel, R. Provided the relationship with the process covariance, R is a correlation matrix, which is symmetric, positively definite and with ones on the main diagonal. 72 Therefore, (ℎ; ) is an even function with respect to h and in particular (0; ) = 1 and lim ℎ→∞ (ℎ; ) = 0. From a practical perspective, these properties implicate a dependence of the predicted response,0, on the tried locations, , which is decreasing at increasing distance between the predicted location, 0 , and the tried one, . Its mathematical properties have made GPR extremely interesting as supervised machine learning technique, in which ( , ) represents a training set. Moreover, the properties of the predictor in Equation (3) and of the kernel introduce a component of artificial intelligence in the algorithm. In fact, differently from conventional regression methods, GPR disregards in estimating the prediction the effect of little impacting (i.e. uncorrelated) input data. 73,74 The functional form of the kernel has to be chosen based on the physics of the process. In this regards, literature, mostly for machine learning application, often resorts to a Bayesian approach, choosing the functional form a priori. Alternatively, as suggested by geostatiscians, 53 an educated guess of the function form can be obtained by the analysis of the variogram. The variogram, γ, was introduced by Matheron and is defined as half the variance of the difference between any pair of realisation at two locations 53,75 : Given the set of tried locations , the scatter plot of the variogram, is usually referred to as variogram cloud. Its relationship with the spatial correlation function suggests that a fitting of the variogram can aid in selecting the best functional form of the kernel R. However, the fitting is usually not carried out on the whole variogram cloud, but on some estimator. One of the most widely adopted is the Matheron estimator, which is unbiased 53 : where (ℎ) = {( , ) ∶ − = ℎ; , ∈ {1, 2, … , }} and the operator # is the cardinality.

Algorithm operational flow
Provided the GPR capability of modelling the spatial correlation-the known spatial correlation properties of measured topographies, and the highlighted shortcomings of the available methods to manage measurement disturbances-GPR seems suitable to be the statistical basis to devise a supervised machine learning technique to overcome these issues. In the case under study, the regressive part of the model in Equation (1) models the deterministic component of the measured electromagnetic 76,77 surface, whilst the stochastic part models the regression error that includes the spatially correlated measurement noise. 54,78 The algorithm developed is the following: 1. Sampling of the measured dataset to build the training data set 2. Selection of the most suitable kernel model 3. Validation of the kernel through the variogram 4. Application of GPR to predict data on the whole surface 5. Evaluation of interpolation residual 6. Identification of spikes, exploiting the residuals (non-measured points are inherently identified) 7. Correction of identified spikes and imputation of non-measured points based on their GPR prediction.
Subsampling of the measured dataset to build the training set has twofold advantages. First, because GPR is a computationally heavy procedure, it reduces computational time and cost. 68 Despite big data analytics and high-performance computing are increasingly resorted to in academia and industry, they require expensive facilities, while subsampling allows running the algorithm on personal computers. Second, it makes the procedure robust. In fact, exploiting the whole original dataset, which presents measurement disturbances, can distort the resulting model to fit the GPR. Conversely, provided that spikes are typically a small fraction of the measured data, subsampling prevents them from negatively affecting the fitting procedure without information loss. 79 The kernel selection can be performed comparing different typical functions, for example, Matèrn, squared exponential, Gaussian, 80,81 with the optimisation criterion aimed at minimising the root mean square error of the residuals. 74 The choice can be validated by exploiting the empirical variogram. Provided a preliminary plane correction is performed, ordinary kriging can be adopted to predict the whole surface, so that no particular assumptions in the regressive part of the prediction are bounded to describe the shape of the topography. 79,82 As far as the identification of measurement disturbances is concerned, non-measured points are inherently identified. On the other hand, spikes require an identification criterion. Provided the discussion in Section 2, residuals of the GPR model has to be normally distributed. 53 Because spikes are anomalous values, typically greatly different, with respect to the measured topography, these can be identified as the outliers of the residuals, identified by the modified interquartile method. 83 As it can be seen, in addition to providing better agreement with the physics of the measurements, the proposed GPR algorithm for managing measurement disturbances is univocally defined, and it does not rely upon possible operator interpretations. Moreover, the procedure manages at once both spikes and non-measured points to avoid biased corrections (cfr. Section 1).

Assessment of the prediction model
The assessment of the GPR is performed by considering some of the most widely adopted surface parameters. 46 They are the root mean square height Sq, the skewness Ssk, the kurtosis Sku and the maximum height Sz. Specifically, Sq and Sz are, respectively, the square root of the second statistical moment and the maximum peak-to-valley distance of the measured surface height z(x,y). They are directly relatable to the measure of the quality and irregularity of optical surfaces-namely root mean square error and peak-to-valley error. 46 Ssk is the normalised third statistical moment of the surface height and indicates a measure of the symmetry of the surface topography with respect to the reference surface (e.g., a negative value indicates few spikes but relatively deep valleys). Lastly, the kurtosis Sku is the normalised fourth statistical moment of the surface height, and such that a spiky surface is characterised by a high kurtosis value, while a bumpy surface by a low kurtosis value. Normalisation is performed to obtain dimensionless parameters. 46 Most of the areal roughness parameters (except for Sz) are defined as surface integral over a definition area, A, considering the surface topography height z(x,y) as a continuous function of the coordinates (x,y) ∈ A. Nevertheless, the measurements from optical instruments contain discrete information, being the measurand a set of pixels mapping the surface topography in the definition area. Therefore, the operational formulas are reported in Equation (6), while the formal definitions can be found elsewhere 76 : where and are the number of pixels in the measured field of view, respectively along the x and y directions. Sp is the maximum peak height and Sv the maximum pit height. It is clear, from their definition, how both the presence of spikes and voids may affect the evaluation of these parameters, and consequently are suitable to assess the performance of the proposed methodology with respect to state of the art.
Additionally, the prediction model is assessed according to the micrographs Topography Fidelity (T FI ). T FI is a metrological characteristic, that is, a characteristic of the measuring instrument with a direct connection to measurement uncertainty, 84,85 expressing the agreement amongst surface topographies when the related uncertainties are negligible by comparison. Considering the abstract definition of Topography Fidelity in the related ISO standard, 84,86 the comparison was performed on the uncertainties STR , 87 evaluated as surface topography repeatability (STR) of the w replicated topography measurements, ( , ), ∈ {1, … , }, that is:

MATERIALS AND METHODS-CASE STUDIES
The comparison will exploit some case studies of industrial interest with different types of technological surfaces. Specifically, quasi-flat surfaces were considered to enhance the comparison of the different disturbances management methods. Moreover, the case studies were chosen so that both the disturbances were not present at the same time, to avoid possible effect of superimposition between spikes and non-measured points corrections, and therefore highlight the suitability of the methodology in both cases. Spikes management were tested on two samples that have been polished and lapped, namely a copper (Cu) and a stainless steel (SS) sample. This type of surface finish is normally required for applications in optics, precision manufacturing and metrology (often to calibrate surface topography measuring instruments). The measurement of these surfaces is inherently more sensitive to spikes due to their nature. Examples of the acquired micrographs are provided in Figure 2. The spikes management was addressed employing: (1) The thresholding method, at 3 . (2) By application of a median filter, as it is based on a robust statistics and with straightforward implementation. 39 (3) Through the proposed GPR-based approach.
The non-measured points are due to complex interactions between the measurand and the instrument, depending on the working principle of the surface topography measuring instrument. For example, in the case of focus variation microscope they are due to poor contrast, for interferometers most typically they are locations that cannot be reached by the light or whose reflection does not return into the objective, which is characteristic of steep, complex and rough surfaces. Non-measured points management by the proposed GPR-based method was tested in comparison with the conventional smooth interpolation 39 on measurements of a Ti6Al4V surface manufactured by electron beam melting (EBM-already shown in Figure 1B). EBM is an additive manufacturing technique, 45,88 which enables great design flexibility, and, thus, has been largely adopted for highly customised parts in several sectors, for example, automotive, aerospace, tooling. Consequently, the characterisation of the related topography is of considerable interest for process and product control. 45 Measurements were performed by a Coherence Scanning Interferometre, the Zygo NewView 9000 shown in Figure 3, hosted at the Technological Surface Metrology Laboratory of Politecnico di Torino, Italy. Measurements were performed on a single field of view, with a Mirau 20 × objective (field of view of 429 μm × 429 μm and pixel size of 429 nm × 429 nm).
Each sample was measured fifteen times, to provide sufficient degrees of freedom, in repeatability conditions. The repeated measurements were then considered as independent for the characterisation. This involved a prior levelling, the management of measurement disturbances and the evaluation of the surface topography parameter presented in Equation (6). No filtering or standard operator was applied, to avoid the removal of structure and enhance disturbances management method differences. 89 Average and confidence interval, at 95% confidence level, of the parameters were computed, to test F I G U R E 3 The Coherence Scanning Interferometre, Zygo NewView 9000, exploited to perform the measurements both the robustness of the disturbance management method to replicated measurements and to assess the presence of significant differences amongst the methods.
As far as the Gaussian process regression-approach is concerned, in order to manage measurement disturbances, the training set was defined by a random subsampling of the surfaces. The ratio between sampled and measured points was set to 0.4% from prior experiences. 79,82 Conventional disturbances management methods and surface topography parameter evaluations have been performed by means of state-of-the-art commercial software, MountainsMap v8.0. 39 The software allows computing standard topographical parameters defined in Section 2.2 and the application of disturbances management technique that are reported in literature, as discussed in Section 1. All the calculations performed by the software are strictly according to the indications provided in relevant ISO standards of the Geometrical Product Specification system, and other international documents in the field. The GPR management method was implemented in MATLAB 2019b and run on a Windows 10 operating system running on a machine with 32 GB of RAM.

Spikes management
The management of spikes by the GPR-based approach requires choosing a kernel. The kernel was chosen after resampling the surfaces, as explained by the method outline in Section 2.1. In particular, the choice depends on the specific surface to be modelled. 74 In the presented case study, for both samples, the spatial covariance structure was suitably modelled by an exponential function, that is, a Matérn function with parameter = 0.5. Figure 4 shows the comparison of the fitting with the chosen kernel and the empirical variogram, validating the adequacy of the choice. As a result, a characteristic correlation length was defined as the distance at which the spatial correlation became negligible, and correspondingly the variogram assumed the value of the process variance. The characteristic correlation length was found of 40 pixels and 257 pixels, respectively for Cu and SS.  Therefore, the whole topography was predicted according to the GPR model, and the residuals of the prediction were computed accordingly. Figure 5A and B shows the normal probability plots (NPP) of the residuals, which highlight the spikes in the long and significant tails. By applying the modified interquartile range method on the residuals, spikes are correctly identified as shown by Figure 5C and D. Figure 6 and Figure 7 show a comparison of the location of the identified spikes by the GPR-based algorithm and the conventional 3 -thresholding and median filter approaches, respectively for Cu and SS samples. Table 1 reports the number of the identified spikes. As it can be noticed from Figure 6B and Figure 7B, the median filter is too severe and detects as spikes also topographical features, as the scratches that are generated during the polishing and lapping. Conversely, the GPR-based algorithm is, in general, less severe, but it is sensitive to local variations of minor extent, as it identifies smaller and isolated spikes on the topography, this is particularly evident in Figure 7A.
Once identified, the spikes can be corrected. The effect of the correction is tested on the areal texture parameters and on the . 76 The texture parameters are shown in Figure 8 and Figure 9 (respectively for Cu and SS sample) as average and 95% confidence interval of the replicated measurements, while the values are reported in Table 2. The considered texture parameters' dispersion represents the sole measurement repeatability contribution, and is computed with a coverage factor equal to a quantile of Student's t distribution with 14 degrees of freedom.
Sz is, by definition, the parameter most sensitive to spikes and allows to show that the thresholding is the most severe approach to correct this type of disturbance, even introducing a systematic difference in the estimation in the SS case study. Because Sq describes the dispersion of the topographical values, it highlights that both the GPR-interpolation

F I G U R E 8
Average and 95% confidence interval of topographical parameters evaluated on the Cu sample after the management of spikes according to the three considered management methods and the median filter introduce a smoothing effect on the topography. However, the median filter is too severe and introduces a systematic difference in the parameter evaluation in the case of the Cu sample. These considerations are reflected in the values of the asymmetry, kurtosis and , even though with less evident effects. The dispersion of the evaluated parameters is informative of the robustness and repeatability of the method. Hence, it is evident that the GPR-based approach outperforms the others. The comparison is also addressed qualitatively, by visual inspection of the corrected topographies, shown in Figure 10 and Figure 11. The greater severity of the median filter, as formerly noted, removes the manufacturing signature, and in its worst spike-managing performance some spikes are still present at the edges of the measured field of view (see Figure 10), increasing the dispersion of Sz, as shown in Figure 8.

F I G U R E 9
Average and 95% confidence interval of topographical parameters evaluated on the SS sample after the management of spikes according to the three considered management methods F I G U R E 1 0 Qualitative comparison of the topography corrected from spikes according to different methods; measurement on Cu sample. Notice that the median filter removes the scratch and is less robust to edge effects.

Non-measured points management
Similarly to the case study of spikes, the correction of non-measured points by the GPR-based method was assessed comparing with a smoothing spline correction, as discussed in Section 3.1. Nonetheless, in this case, the greater complexity of the topography required the choice of a squared exponential kernel, with a characteristic correlation length of 43 pixels, and validated by the variogram shown in Figure 12. Given the greater scales of the topography, the correction was addressed directly by a quantitative comparison of surface topography parameters, as shown in Figure 13, and by the summarised in Table 3. Systematic differences between the two approaches can be appreciated examining Sz and Sq (95% confidence level). In particular, the spline interpolation led to a systematic overestimation of both parameters, and consequently also of . As expected, skewness and kurtosis were less suceptible to highlight differences between the two methods. Nonetheless, differences between the two methods are coherent for all the parameters. This is evident when inspecting the distribution of their difference (see Figure 14), which are not normally distributed, clearly indicating a systematic difference between the two approaches. Hence, the proposed algorithm based on a Gaussian process regression yields robust parameters estimations also in the case of non-measured points.
TA B L E 3 u T FI after the correction of non-measured points according to the two considered methods

NM-points correction method / µm
Interpolation 2.475 1.570 GPR-based algorithm F I G U R E 1 3 Average and 95% confidence interval of topographical parameters evaluated on the EBM sample after the management of non-measured points according to the two considered management methods F I G U R E 1 4 (A) NPP and (B) histogram of NM-points corrected according to the GPR (red), the smooth interpolation (cyan) and their differences (black). Notice the strong deviation from normality of the difference, which allows concluding on the presence of a systematic factor between the two approaches

CONCLUSIONS
The measurement of surface topography is relevant in several industrial fields to assess the quality of manufactured surface components and to guarantee their functionality. Surface topography measurements are affected by some disturbances, namely spikes and non-measured points, that the literature has demonstrated to impact severely on the topographical characterisation. Consequently, they need to be identified and corrected properly. However, the current literature presents several different approaches with robustness and sensitivity criticalities, which disregard the presence of the spatial correlation in the topographical measurements. The present paper proposes a supervised machine learning approach based on a Gaussian process regression to identify and correct disturbances of surface topography measurements. Differently to other literature alternatives, this approach offers definitive advantages as it is univocally defined, it is capable of managing, at once, both spikes and non-measured points, and it inherently takes into account the spatial correlation. Such approach was tested on different types of quasi-flat technological surfaces showing better sensitivity in identifying the disturbances with respect to literature alternatives, and also overall improved performance in the disturbance correction, both in terms of robustness and precision, as far as the characterisation of surface topography parameters is concerned. The presented approach relied on random sampling to train the machine learning algorithm, future work will investigate the sensitivity of the method on the size of the training dataset and on the sampling criterion, for example, latin-hypercube and space filling designs, to improve the trade-off between the computation and the prediction performances. Future development will deal with the concurrent presence of spikes and non-measured points in quasi-free form surfaces.

A C K N O W L E D G E M E N T S
This work has been partially supported by "Ministero dell'Istruzione, dell'Università e della Ricerca". Award "TESUN-83486178370409 finanziamento dipartimenti di eccellenza CAP. 1694 TIT. 232 ART. 6".

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.