# Generation and Reproduction of Random Rough Surfaces

Arthur Francisco, Noël Brunetière

<figure><img src="image_1.png" /></figure>

## To cite this version:

Arthur Francisco, Noël Brunetière. Generation and Reproduction of Random Rough Surfaces. Lubricants, 2025, 13 (11), pp.506. 10.3390/lubricants13110506. hal-05482212

HAL Id: hal-05482212

https://hal.science/hal-05482212v1

Submitted on 28 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Article

# Generation and Reproduction of Random Rough Surfaces

Arthur Francisco<fnref n="1" />,*,† D and Noël Brunetière <fnref n="2" />,† D

**Abstract:** In both mixed lubrication and dry contact studies, statistically meaningful conclusions often require a large set of numerical rough surfaces. While such surfaces can be acquired through metrological tools such as interferometers, they rarely exhibit the exact height or spatial parameters of interest, and the available datasets are typically limited. Although the numerical generation of rough surfaces is not a new subject, its relevance has grown and important challenges remain. Building on their earlier work, in which a new generation method was introduced, the authors extend its scope to produce surfaces with prescribed height and spatial parameters, under non-periodicity constraints and arbitrary orientation anisotropy. In addition, they propose the reproduction of existing rough surfaces for future AI training applications, and highlight topographic patterning as the next major challenge to address.

**Keywords:** Rough surface generation, statistical moments, autocorrelation function, non-Gaussian random surfaces, genetic algorithm

## 1. Introduction

Whatever the real surface under study, as long as it is examined at a sufficiently small scale, irregularities will appear. This provides an intuitive approach to roughness: microscopic deviations from a macroscopic mean line. Accordingly, to characterize the roughness level of a surface, several parameters have been defined, addressing the height distribution, the spatial arrangement of heights, the material ratio curve, and so forth. The reader may refer to Leach [1] for a comprehensive overview of rough surface characterization, particularly the ISO 25178 parameters, some of which are used in the present study. The most obvious source of surface generation is the mechanical industry. Depending on the requirements, the desired surface finish may involve more or less numerous and complex processes. These processes consist of successive material removals, followed by finishing operations. It is easy to understand that the final surface results from the cumulative effect of material removal, large at the beginning and small at the end of the process. The final surface, in the mathematical sense, is therefore a collection of heights, as a superposition of wavelengths of decreasing length.

The layered topography of common industrial surfaces can be numerically approximated by a simple mathematical tool, in order to recreate such surfaces. In this operation, fidelity to the mean line is crucial, since the long wavelengths that compose it are the most contributive to the load-carrying capacity of nearly conformal lubricated surfaces. The

**Citation:** Francisco, A.; Brunetière, N. Title. *Lubricants* **2025**, *1*, 0. https://doi.org/

Received:
Accepted:
Published:

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2025 by the authors. Submitted to *Lubricants* for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Institut Pprime, CNRS, Université de Poitiers, ISAE-ENSMA, 16021 Angoulême, France; arthur.francisco@univ-poitiers.fr</footnote>

<footnote n="2">² Institut Pprime, CNRS, Université de Poitiers, ISAE-ENSMA, 86962 Futuroscope Chasseneuil, France; noel.brunetiere@univ-poitiers.fr</footnote>

<footnote>* Correspondence: arthur.francisco@univ-poitiers.fr; Tel.: +33-545-673-249</footnote>
shortest wavelengths, i.e. roughness, are less involved in load support. However, when lubricated surfaces in relative motion approach each other, the occurrence of asperity contact on both sides of the film increases. As long as the event remains globally rare, the debris and thermal energy produced are absorbed: this is the situation of mixed lubrication. When surfaces approach further, occurrences multiply until degradation occurs, sometimes rendering the mechanism inoperative. It is therefore clear that the finishing quality of surfaces depends on their intended use, since systematically pushing it to a higher level would represent a prohibitive cost. It becomes evident that being able to generate a numerical surface with controlled topography represents a major challenge for accurate simulations.

To characterize a rough surface, we generally need two sets of information: one to describe the statistical profile of the heights themselves, and another to describe their spatial layout. The four statistical moments handle the profile – the mean $\mu$, the squared standard deviation $\sigma^2$, the skewness coefficient $Ssk$, and the kurtosis coefficient $Sku$ – while the autocorrelation function – acf – defines the layout.

In mixed lubrication, the parameters $Ssk$ and $Sku$ have a clear influence on load capacity and maximum pressure [2]. It is even advisable to increase $Sku$ in order to reduce friction [3], provided that the values of the roughness parameter $Ssq$ – the height standard deviation $\sigma$ – and the skewness coefficient $Ssk$ remain small. Sedlacek et al. [3] concluded that surfaces should be machined with these parameters in mind.

More recently, Ba et al. [4] confirmed the importance of these parameters as well as the influence of the manufacturing process, through the shape of textures imprinted on the surface. This conclusion – ultimately long integrated – led He et al. [5] to propose a type of topography for frictional surfaces. Finally, Pawlus et al. [6] suggested extending the range of parameters influencing lubrication to the category of spatial parameters. But this is not really new, since more than forty years have passed since the first proposals in this direction.

Since the mean $\mu$ and the standard deviation $\sigma$ of the heights are additive and multiplicative coefficients, respectively, the tribological behavior of a surface would thus be governed by four parameters: $Ssk$, $Sku$, $l_ζ$, and $l_η$ – respectively skewness, kurtosis, greater and smaller correlation lengths. Reality, however, is somewhat different, which is common sense: how could there exist a bijection between the space of heights and a four-dimensional space? Technically, the relation is surjective but not injective: for identical parameters, several possible surfaces exist. Even if one could establish a perfect correspondence between a surface and a reduced set of parameters, the processes themselves have intrinsic variability that affects both the distribution of heights and their spatial arrangement. Yang et al. [7] examine the relevance of rough surface parameters as defined in ISO 25178. Based on the correlation level among parameters from a dataset of one thousand surfaces, they identify a set of 13 independent parameters: $Sa$, $Ssk$, $Sku$, $Sp$, $Sv$, $Vvv$ from Abbott curve-related parameters, and $Str$, $Spk$, $Smr1$, $Sxp$, $Spd$, $Spc$, $Sdq$ from spatial height arrangement parameters. Among the latter, only $Str$ is controlled by convolution, and interestingly, $Sal$ appears to be linked to $Str$ although these two parameters are theoretically independent. In any case, the authors address a crucial issue, namely determining the minimal set of parameters capable of fully characterizing a rough surface. In their study Chen et al. [8] enhance a random switching system by introducing constraints for spatial features in the height set, enabling precise control of ISO 25178 roughness pa-
rameters (*Spd*, *Spc*, *Sdq*, *Sdr*) while preserving surface topography characteristics. This approach also compensates for one of the drawbacks of convolution-based methods, namely the absence of control over spatial parameters other than the correlation lengths. Chen *et al.* [9] provide a comprehensive overview of the context of rough surface generation. Building on the method of Francisco and Brunetière [10], with several improvements, they show that the statistical parameters of the generated surfaces adequately cover the requirements of industrial surfaces. They also emphasize – importantly – that convolution-based methods cannot directly reproduce special structures, such as those observed on worn surfaces, which has also been confirmed by Borodich *et al.* [11]. Both aspects will be addressed in the present work.

Another source of rough surfaces is biological in origin: dental surfaces. Indeed, nature often creates surface patterns to meet functional needs such as hydrophobicity, adhesion, abrasion, and so forth. In the present case, the patterns are related to the way the tooth has grown, and more specifically to the enamel. Finally, the different types of abrasion – attrition, scratches, chipping, chemical attacks, etc. – further add to the complexity of the surface. The study of dental microwear allows solid hypotheses to be made about the diet of the species concerned – at least its last meals. Many paleontologists use dental microwear patterns to classify herbivorous species into categories: grazers, browsers, frugivores, or mixed feeders. Dental surfaces exhibit such variability in terms of roughness that classifications are often difficult and only weakly reproducible from one group of individuals to another, see Francisco *et al.* [12] for details. Since the acquisition of dental surfaces is a rather time-consuming process, and machine learning – particularly convolutional methods – requires large volumes of data, the ability to generate realistic dental surfaces for each dietary group would make it possible to build an accurate classification tool.

The generation of numerical surfaces poses a significant challenge, driven by its importance in two distinct fields: the development of predictive models for lubrication and the creation of comprehensive dental datasets for training machine learning algorithms. We will therefore first focus on making our method for generating rough digital surfaces more robust and more general. In a second step, we will experiment with ways in which the algorithm can reproduce existing surfaces in terms of both height and spatial parameters.

## 2. Materials and Methods

### 2.1. Modified Hu and Tonder method

There exist several methods for generating rough surfaces, as discussed by Pawlus *et al.* [13] and Wang *et al.* [14], and we have chosen the one proposed by Hu and Tonder [15] as the working basis, which can be summarized in five steps, with steps 3 and 5 being modified.

#### 2.1.1. Step 1: Choice of parameters

Let us denote by $z_i$ a height, realization of a random variable $z$, and $n$ the sample size. The four statistical moments $\mu, va = \sigma^2, Ssk, \text{ and } Sku$ are considered representative of a height unimodal distribution. Since $Ssk$ and $Sku$ are obtained after centering and normalizing the distribution, we can choose $\mu = 0$ and $va = 1$; recall that:
$$ \mu = \frac{1}{n} \sum_{i=1}^{n} z_i ; va = \sigma^2 = \frac{1}{n} \sum_{i=1}^{n} (z_i - \mu)^2 ; Ssk = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{z_i - \mu}{\sigma} \right)^3 ; Sku = \frac{1}{n} \sum_{i=1}^{n} \left( \frac{z_i - \mu}{\sigma} \right)^4 \quad (1) $$

The reader will not be surprised by the use of biased statistical estimators, since for a numerical surface we typically have $n > 512 \times 512$, thus making the bias entirely negligible.

### 2.1.2. Step 2: Determination of the digital filter H
Analytical surface generation

Before addressing numerical surface generation, let us start with the analytical form of a surface $z(x,y)$ obtained by filtering $h(x,y)$ of a white noise $\eta(x,y)$:

$$ z(x,y) = h(x,y) * \eta(x,y), \text{ where } * \text{ denotes the convolution operator} $$

The autocorrelation function $acf(z)$ can be written, omitting spatial variables:

$$ \begin{aligned} acf(z) &= z \otimes z = (h * \eta) \otimes (h * \eta) \\ &= (h * \eta) \otimes (\eta * h), \text{ by commutativity of convolution} \end{aligned} $$

$z$ being real, $acf(z)$ is real and even.

Let $\mathcal{F}(f)(\omega_x, \omega_y)$ denote the Fourier transform of a function $f(x,y)$, and define:

* $Z(\omega_x, \omega_y) = \mathcal{F}(z)$

* $A(\omega_x, \omega_y) = \mathcal{F}(\eta)$

* $H(\omega_x, \omega_y) = \mathcal{F}(h)$

* $\text{PSD}(\eta) = |A|^2$

* $\text{PSD}(z) = |Z|^2$

Using Fourier transform properties regarding convolution and correlation products, we obtain:

$$ \text{PSD}(z) = Z\bar{Z} = H \times \text{PSD}(\eta) \times \bar{H} \quad \text{where } \bar{z} \text{ denotes complex conjugation} $$

Since $\eta$ is chosen as white noise, $\text{PSD}(\eta)$ is a constant $c$. An intuitive way to see this is through the Wiener-Khinchin theorem, which states that $\text{PSD}(z) = \mathcal{F}(acf(z))$. $\eta$ autocorrelation is a Dirac peak centered at $(0,0)$, implying the presence of all frequencies with equal contributions in its Fourier series expansion.

Consequently, there exists a simple relation between the Fourier transform of the autocorrelation of $z$ and the squared modulus of the digital filter $H$: $\mathcal{F}(acf(z)) = c|H|^2$. This relation is fundamental since the knowledge of $acf(z)$ directly yields $|H|^2$ because the constant $c$ is absorbed into the variance of $z$, which will be standardized anyway.

Thus,

$$ \mathcal{F}(acf(z)) = |H|^2 \implies H = \sqrt{\mathcal{F}(acf(z))} e^{i\phi(\omega_x, \omega_y)} $$
We may choose a nonzero phase function $\phi$. However, a linear function is of little practical interest, while a nonlinear function makes the convergence of the process described later, in Step 5, more difficult. Moreover, we do not fully control the effect of a nonlinear function in creating patterns through distortion. For all these reasons, $\phi$ is set to zero. Since $PSD(z)$ is real and even, $H$ is real and even as well, which guarantees that $h$ is real and even - symmetric about to the origin.

$$h(x, y) = \mathcal{F}^{-1}\left(\sqrt{\mathcal{F}(acf(z))}\right) \qquad (2)$$

### Numerical surface generation

Transposing the above derivation to the numerical case requires some adjustments, particularly concerning the quality of the random signal, which is not strictly white noise. More precisely, since solving the system $Z^2 = H A^2 \tilde{H}$ is not computationally feasible, we retain the expression of $H$ but correct $A$:

$$H = \sqrt{\mathcal{F}(acf(z))} \quad \text{and} \quad A \leftarrow \frac{A}{|A|} \qquad (3)$$

#### 2.1.3. Step 3: Determination of initial statistical parameters

After filtering the initial random signal $\eta$, the statistical moments are modified as follows:

$$\begin{align}
va_z &= \sigma_z^2 \approx \sum_{i=1}^{n} h_i^2 \nonumber \\
Ssk_z &\approx \frac{\sum_{i=1}^{n} h_i^3}{\left(\sum_{i=1}^{n} h_i^2\right)^{\frac{3}{2}}} Ssk_{\eta} \tag{4} \\
Sku_z - 3 &\approx \frac{\sum_{i=1}^{n} h_i^4}{\left(\sum_{i=1}^{n} h_i^2\right)^2} (Sku_{\eta} - 3) \nonumber
\end{align}$$

The formulas in Equation 4 are approximations that assume all multiple and weighted correlations of $\eta$ to be strictly zero. In practice, for large correlation lengths and a non-Gaussian surface, this assumption may even lead to a violation of Pearson's inequality. As an illustrative example, we consider the generation of a surface whose parameters are listed in Table 1.

<table><caption>**Table 1.** An example illustrating the limitations of the relations in Equation 4.</caption><thead><tr><th>width (µm)</th><th>height (µm)</th><th>w</th><th>h</th><th>Ssk<sub>z</sub></th><th>Sku<sub>z</sub></th><th>l<sub>ζ</sub> (µm)</th><th>l<sub>η</sub> (µm)</th></tr></thead><tbody><tr><td>200.0</td><td>200.0</td><td>1024</td><td>1024</td><td>-3.0</td><td>15.0</td><td>30.0</td><td>30.0</td></tr></tbody></table>

Table 2 shows the results obtained using the relations of system 4. Pearson's inequality is clearly incompatible with the values of $Ssk_{\eta}$ and $Sku_{\eta}$, and to proceed with the process, $Ssk_{\eta}$ and $Sku_{\eta}$ must be corrected according to $Sku_{\eta} = Ssk_{\eta}^2 + 1$.
<table><caption>**Table 2.** Resulting statistics Ssk and Sku with the parameters of the Table 1.</caption><thead><tr><th>Ssk<sub>h</sub></th><th>Sku<sub>h</sub></th><th>Ssk<sub>η</sub></th><th>Sku<sub>η</sub></th></tr></thead><tbody><tr><td>21.8</td><td>1537.3</td><td>-141.0</td><td>8187.9</td></tr></tbody></table>

The expressions, Equation 4, as used in the literature, notably by Bakolas [16], do not exploit the fact that $\mu_h$ can be set to 0:

$$z = h * \eta = (\tilde{h} + \mu_h) * \eta = \tilde{h} * \eta \quad (5)$$

because $\mu_\eta = 0$. As proved in Appendix [A.2-A.4], we have the following approximations:

$$ \begin{align} \sigma_z^2 &\approx n\sigma_h^2 \nonumber \\ Ssk_z &\approx \frac{1}{\sqrt{n}}Ssk_h Ssk_\eta \tag{6} \\ Sku_z - 3 &\approx \frac{1}{n}Sku_h (Sku_\eta - 3) \nonumber \end{align} $$

instead of the classical ones, Equation 4.

These relations allow us to determine the statistical parameters after filtering, given the initial ones. One can deduce from these relations that a Gaussian signal is transformed into a Gaussian-filtered signal: $Ssk_z$ remains zero, as does the excess kurtosis $Sku_z - 3$. An interesting effect is the division by $n$ or $\sqrt{n}$ of the moments of $z$: unless the initial moments are very large, the moments of $z$ will be close to 0. In other words, convolution makes the initial signal more Gaussian, which suggests numerical challenges. This convergence toward the normal distribution was predictable since, by virtue of the Central Limit Theorem (CLT), a sum of random variables tends toward a normal distribution. Since $z_k = \sum_i h_i \eta_{i+k}$, the CLT applies. One can even conjecture that, starting from a given signal $\eta$, repeated convolutions will yield a distribution arbitrarily close to the normal law.

Finally, the reader, perhaps a future user, must be warned that these are approximations: any nonzero autocorrelation of $\eta$, however small, will move these relations away from exact equality.

### 2.1.4. Step 4: Transformation of the initial random surface

Since the target parameters are Ssk and Sku, the distribution $\eta$ must be transformed so that it exhibits the desired statistical properties. The conventional approach is to use Johnsons system of transformations [17], originally designed to normalize a distribution. Since the transformations are strictly monotonic, they can be inverted.

As shown by Francisco and Brunetière [10], it is sometimes difficult, or even impossible, to reach certain values of Ssk and Sku using Johnsons system of transformations following Hills algorithm. The improvement proposed by the authors consists in using analytical functions – the tangent and exponential functions – to generate a series $\eta$ with the correct statistical moments.
2.1.5. Step 5: Computation of z by filtering

Once the initial distribution $\eta$ has been determined, it is filtered by $h$, as follows:

$$Z = HA \quad \text{hence} \quad z = \mathcal{F}^{-1}(HA) \qquad (7)$$

When applying the filter $h$ to the resulting random series $\eta$, by construction of $h$ the obtained autocorrelation function is strictly identical to the target one. The issue is that we did not start from a perfectly white noise, but only from an approximation of it. Consequently, the filtering process – which is a form of smoothing – does not yield the correct $Ssk$ and $Sku$ values for the series $z$.

The idea, presented Francisco and Brunetière is then to identify the permutation $O_z$ that transforms the obtained series $z$ into a sorted series: the heights may not be exactly correct, but their arrangement is. It then suffices to substitute these approximate heights $z$ with a series $z'$ of correct heights by applying the inverse permutation $O_z^{-1}$.

As a result, the autocorrelation function is somewhat modified, but the heights themselves are exact. For a non-Gaussian resulting series $z$, the larger the correlation lengths, the stronger the required smoothing, and the farther the initial series will deviate from normality. It is therefore sufficient to embed these actions within an iterative scheme: after replacing $z$ with $z'$, the filter is reapplied, and so on. When the initial series $\eta$ has been properly generated, $z$ is close enough to $z'$, and after only a few iterations the resulting autocorrelation function matches the expected one; Figure 1 illustrates the general structure of the algorithm.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Flowchart of the modified Hu and Tonder algorithm for rough surface generation. Block 1 can be replaced by the recovery of an existing surface. Block2 can be replaced by a prescribed texture.</figcaption></figure>

By following the modified algorithm proposed by Hu and Tonder, it is possible to generate a rough surface from virtually any set of height or spatial parameters. However, some realizations exhibit a lack of realism, as will be discussed in the results section, where corrective strategies will also be introduced. Moreover, the use of an antisymmetric exponential function imposes an additional constraint, which can in fact be removed in order to simplify the algorithm, as we shall see immediately.
## 2.2. No magic universal function

Francisco and Brunetière [10], have shown that two analytical expressions are sufficient to generate all pairs (Ssk, Sku), provided they do not fall below Pearson's limit. The functions used to generate n points are the tangent function Equation 8 and an antisymmetric exponential function, Equation 9.

$$y_i = \tan\left(-\frac{\pi}{2}(1-a) + (i-1)h\right), \text{ with } h = \frac{\pi}{2(n-1)}(2-a-b) \quad (8)$$

$$y_i = \operatorname{sgn}(x_i)(1 - \exp(-|x_i|)), \text{ with } x_i = 1 - \frac{1}{a} + (i-1)h, \text{ and } h = \frac{1}{n-1}\left(\frac{1}{a} + \frac{1}{b} - 2\right) \quad (9)$$

with $i \in [1, n]$ and $a, b \in [0, 1]$.

The dimensionless parameters a and b are determined using a genetic algorithm, as will be discussed later, so that {$y_i$}$_{i=1,n}$ exhibits the desired statistical moments, namely skewness (Ssk) and kurtosis (Sku). The function in Equation 8 is used in the case Sku ≥ 1.34Ssk² + 1, and the function in Equation 9 in the case Sku < 1.1Ssk² + 1. The authors have shown that these two functions must be combined within the limits 1.1Ssk² + 1 < Sku < 1.34Ssk² + 1 in order to obtain optimal results.

To avoid dealing with subdomains, one may ask whether a universal function exists that would allow generating a dataset of n points for any admissible (Ssk, Sku) pair. Unfortunately, to the best of our knowledge, no such parametric function exists. Indeed, as illustrated in Figure 2, this function would need to approximate both a binary dataset – dashed red line – and extreme cases of Ssk, Sku – dashed green line. It can be noted that the functions represented by solid red lines have their reciprocal functions in solid green lines. Thus, the universal function h(x) would be written as $h(x) = (1 - u)f(x) + uf^{-1}(x)$, with $u \in [0, 1]$ a partition parameter, but it significantly complicates the integration of powers 2, 3, and 4 required for the statistical moments.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The principle of a parametric universal function able to cover the whole Ssk, Sku domain. The red curves may approach the Pearson boundary, whereas the green curves can lead to high values of Ssk and Sku.</figcaption></figure>
Thus, a more general function than the antisymmetric exponential should be sought, so as to eliminate the need to handle the intermediate region.

## 2.3. Choice of a more versatile function than exp

Francisco and Brunetière have shown that the exponential function $f(x) = \operatorname{sgn}(x)(1 - \exp(-|x|))$ is a relevant choice for approximating Pearson's limit; however, it becomes less suitable when moving away from it – $Sku > 1.1Ssk^2 + 1$.

The idea is to complement it with a function $g(x)$ that can exhibit peaks at the integration bounds, while still ensuring that the resulting function $h(x) = f(x) + g(x)$ is analytically integrable up to the 4th power. The simplest function that meets this requirement is the cubic function $g(x) = x^3$. To provide flexibility to $h(x)$, a parameter $c$ must be added to $a$ and $b$. Thus:

$$h(x) = \operatorname{sgn}(x)(1 - \exp(-|x|)) + cx^3 \quad (10)$$

with $x \in [1 - \frac{1}{a}, -1 + \frac{1}{b}]$ and $a, b \in [0, 1]$. The issue is that the size of the interval will have a strong influence on $g(x)$, which must be compensated by the value of $c$. To avoid this phenomenon, $g(x)$ is restricted to the interval $[-1, 1]$, and, as with $a$ and $b$, $c \in [0, 1]$. We then can write:

$$h(x) = \operatorname{sgn}(x)(1 - \exp(-|x|)) + c \left( \frac{x}{\frac{1}{a} + \frac{1}{b} - 2} \right)^3 \quad (11)$$

and it is obvious that $(a, b) = (1, 1)$ must be prevented.

The next question is: why restrict the use of the modified exponential to the domain $Ssk^2 + 1 \le Sku \le 1.34 Ssk^2 + 1.8$? Based on the material curve of a rotor worn surface, Figure 3 shows that the tangent function is more realistic than the modified exponential, and to confirm this, it suffices to switch to logarithmic scales, Figure 4.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Best fit of a rotor surface material curve with tangent and a modified exponential functions; $Ssk = -4.3$ and $Sku = 27.0$.</figcaption></figure>

The rotor heights were measured using a white-light interferometric microscope. From the original surface containing approximately one million points, a representative subset of about 2,000 points preserving the minima and maxima was retained. The adjusted heights, on the other hand, are generated numerically.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Best fit of a rotor surface material curve with tangent and modified exponential functions. Logarithmic scale.</figcaption></figure>

Finally, it can be easily verified that the distributions obtained with modified exponentials are unimodal when moving away from Pearson's limit, as shown in Figure 5 for the previous case. In the previous version of the algorithm, a mixture of two distributions was used, which could introduce undesired bimodality, see Chen *et al.* [9].

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Height histogram of the rotor surface showing unimodality.</figcaption></figure>

## 2.4. On the relevance of the modified exponential function

To ensure that the tangent and modified exponential functions are capable of covering the entire range of *Ssk*, *Sku* values, it is necessary to determine the boundaries *Sku*<sub>min</sub>, *Sku*<sub>max</sub> for *Ssk* ∈ [0, √*n*].

To explain the maximum value $S_{skmax} = \sqrt{n}$, one must recall that $S_{sk}$ is a measure of the asymmetry of a distribution.

If we take *n* points such that $y_1 = \sqrt{n-1}^{-1}$ and $y_{2..n} = \sqrt{n-1}$, we automatically obtain $\mu = 1$ and $\sigma = 1$. Consequently:
$$ Ssk = \frac{1}{n}(-\sqrt{n-1}^3 + (n-1)\sqrt{n-1}^{-3}) \approx -\sqrt{n-3} \approx -\sqrt{n} \quad (12) $$

To determine the boundaries of the (Ssk, Sku) domain, one could scan the entire set of values by varying a, b, and c. However, most of the computations would be unnecessary since they fall within the interior of the domain. An interesting alternative is then to explore only the boundary points.

### 2.4.1. The Pareto Front

Suppose we have two conflicting objectives to optimize, for example the comfort of a vehicle – to be maximized – and its price – to be minimized. It is easy to imagine that when one improves one of the two objectives, the other deteriorates. One can then compile a list of pairs (comfort, price) that result from the best trade-offs.

Let the two objectives be $f_1 = \frac{1}{\text{comfort}}$ and $f_2 = \text{price}$, both of which must be minimized. In a graph ($f_1, f_2$), one obtains a convex curve, where each point is non-dominated. A point is said to be dominated when both objectives are simultaneously improved by another point as explained in Figure 6.

<figure><img src="image_8.png" /><figcaption>Figure 6. The Pareto front principle.</figcaption></figure>

How, then, can we directly determine the Pareto front?

The first intuitive approach is to fix a value of k, say k = 1, and minimize $f^0(x_i) = f_1(x_i) + k f_2(x_i)$. Once the $x_i$ are determined, they are reused as initial values to minimize $f^1(x_i)$ with $k = k \pm \Delta k$, and so on. The drawback of this method is that it does not guarantee a homogeneous distribution of solutions along the front: one may end up in a situation where particular points of the front aggregate most of the nearby solutions. We therefore need to choose a minimization algorithm that is robust against local minima and ensures a sufficient and homogeneous density of the front.

Among the major families of optimization methods are metaheuristics. These are generally iterative stochastic algorithms: the random component of the search helps to avoid pitfalls such as local optima, while successive iterations build solutions that rely heavily on the search history. A comprehensive presentation of multi-objective methods
– and more generally metaheuristics – can be found in Collette and Siarry [18]. Most metaheuristics implement natural paradigms such as:

* the evolution of species (evolutionary algorithms)

* the cooling of crystalline structures (simulated annealing)

* ant colonies

Unlike gradient-descent methods, metaheuristics are essentially combinatorial approaches that do not rely, for example, on derivative computations. Nevertheless, they all suffer from a high number of calculations and cost function evaluations, and in some cases from the need to tune a large number of parameters (as in simulated annealing). Evolutionary algorithms are not exempt from the computational burden; however, the analytical determination of statistical moments accommodates this very well. Moreover, there are relatively few parameters to tune. In fact, one could even say that no real tuning is required, since the adjustable parameters have default values that already yield very good results.

### 2.4.2. General principles of evolutionary algorithms

For a problem of minimizing a function $f(x_i)$, the idea is as follows:

* starting from a set of size $p$ of parameter sets $x_i$ – or equivalently a set of $p$ individuals – one generates a set of $p$ solutions $f$

* the $p$ solutions are ranked in ascending order, and the best individuals are set aside

* the $x_i$ thus identified are recombined with others, forming a new set

* randomly, some of the $x_i$ are perturbed

After iterating this process many times, the $x_i$ obtained are those for which $f$ is minimal. One may then ask: how are the $x_i$ "recombined" and how are they "perturbed"? Inspired by genetic coding, one works with the binary representation of the parameters. Thus, if an $x_i$ is encoded on 32 bits and the parameter set consists of $p$ such $x_i$, the chromosome of $f(x_i)$ will consist of 32 genes, i.e., $32p$ bits. To combine two chromosomes, one simply takes a portion of one and completes it with the other.

### 2.4.3. The NSGA-II algorithm

It would add little value here to provide a complete inventory of the various multi-objective evolutionary techniques. To convince the reader of the abundance of options, one may refer to the non-exhaustive list compiled by Coello *et al.* [19]. In this list, the NSGA-II technique of Deb *et al.* [20] is highlighted; this is the one we adopt. This choice is based on two pragmatic considerations. The first is its reputation for efficiency, as the Pareto front is determined "quickly" and is uniformly populated. The second is the Python code made freely available in the Pymoo library [21]. Although the program for computing statistical moments is written in Fortran, Pymoo can use it via the free library Fmodpy [22]. Finally, the computation can be easily parallelized using the free multiprocessing library [23].

### 2.4.4. Analytical expressions of the lower boundaries

The NSGA-II algorithm was used to determine the following boundaries:

* $f_{e,min}(Ssk) = \min(Sku)$, modified exponential basis, lower bound

* $f_{t,min}(Ssk) = \min(Sku)$, tangent basis, lower bound

* $f_{e,max}(Ssk) = \max(Sku)$, modified exponential basis, upper bound
* $f_{t,max}(Ssk) = \max(Sku)$, tangent basis, upper bound

343

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Analytical function boundaries in log scale. Note that for rough surfaces, *Ssk* is most of the time negative.</figcaption></figure>

The upper bounds were not fitted, since their analytical form is of no practical significance; the grey lines are only used to delimit the relevant region.

From the results of Figure 7, in logarithmic scales to visualize the entire domain, several essential elements can be identified:

344

345

1. The function that best approximates $f_{e,min}(Ssk)$ is $Sku = 1.017 Ssk^2 + 1$, which is very close to Pearson's limit. This also shows that the multi-objective optimization works very well.

346

347

348

2. The function that best approximates $f_{t,min}(Ssk)$ is $Sku = 1.24 Ssk^2 + 1.8$, at least up to $Ssk = \sqrt{\frac{n}{2}}$; we shall return to this point later, in Appendix. In Francisco and Brunetière [10], the boundary for using the tangent function was set to $Sku = 1.34 Ssk^2 + 1$, which is consistent since, near the boundary, the transition from analytic to numerical results suffers from some inaccuracies.

349

350

351

352

353

3. As a consequence of the previous point, the tangent/exponential transition is set at $Sku = 1.34 Ssk^2 + 1.8$ to ensure the value of $Sku$ at $Ssk = 0$.

354

355

4. The theoretical limit of $f_{t,max}$ is approximately $n$, as previously observed. However, it has also been shown that this limit depends on the size of the integration interval $[1 - \frac{1}{a'}, -1 + \frac{1}{b}]$, which explains why it does not form a horizontal line.

356

357

358

5. The upper bound of the exponential $f_{e,max}$ lies well above the transition curve, unlike in previous works [10], thanks to the introduction of the cubic function $x^3$. Thus, there is no longer a transition zone, and the weakness of the previous algorithm, pointed out by Chen *et al.* [9], no longer exists.

359

360

361

362

We can now define two working regions depending on the position of a point (*Ssk*, *Sku*) relative to the transition, see Figure 8. If the point lies above, the tangent function is used; otherwise, the exponential basis is applied.

363

364

365
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Analytical function boundaries.</figcaption></figure>

The pair of tangent and modified exponential functions adequately covers the required values (*Ssk*, *Sku*) compatible with the Pearson limit and the number *n* of heights.

## 2.5. Reference surfaces for reproduction

As mentioned in the introduction, dental surfaces pose a challenge for algorithms dedicated to the generation or reproduction of numerical surfaces due to the spatial arrangement of heights. The unworn surface exhibits smooth patterns and resembles skin. Additional features may include deep striations left by the friction of silica grains, or even pits caused by ingested hard particles.

The reference surface used throughout the reproduction process was selected from herbivore dental surfaces. It was extracted from an original surface and reshaped into a rectangular size, in order to demonstrate the versatility of the generation program, Figure 9.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Dental surface of a grazer-type herbivore. 1024 × 512 px surface, size = 132 × 66 µm. *Ssk* = -1.2, *Sku* = 6.0, *Stri* = 1.4, *Sal* = 5.2 µm. The white ellipsis is the autocorrelation function cut at *z* = 0.5. Note that *Stri* stands for *Str*⁻¹.</figcaption></figure>
The capability of the proposed algorithm is also tested on typical industrial surfaces. The objective is to reproduce them not only in terms of the imposed parameters but also with respect to realism. Five surfaces are therefore considered, Figure 10.

378

379

380

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Industrial rough surfaces used to test the algorithm. The white ellipse are the autocorrelation functions cut at $z = 0.5$. Note that *Stri* stands for *Str*⁻¹.</figcaption></figure>

The following table, Table 3, provides a qualitative overview of the characteristics of the six surfaces considered.

381

382
<table><caption>**Table 3.** Qualitative overview of the six surfaces, highlighting the difficulty of each characteristic.</caption><thead><tr><th>Surface</th><th>non-Gaussian</th><th>Anisotropy Stri</th><th>Correlation length l<sub>ζ</sub></th></tr></thead><tbody><tr><td>Dental</td><td>++</td><td>-</td><td>+++</td></tr><tr><td>Rotor</td><td>+++</td><td>-</td><td>-</td></tr><tr><td>Milled</td><td>-</td><td>++</td><td>++</td></tr><tr><td>Ground</td><td>-</td><td>+++</td><td>-</td></tr><tr><td>Polished</td><td>+</td><td>+</td><td>++</td></tr><tr><td>Mirror-polished</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

## 3. Results and discussion

The generation of homogeneous isotropic rough surfaces is no longer a challenge, even though the visual appearance is not always realistic. In the presence of anisotropy, with an arbitrary orientation, the situation becomes more complex since it departs from the assumption of periodicity required for the use of FFT. This departure from periodicity becomes increasingly significant as the correlation lengths increase. We will therefore first focus on generating oriented rough surfaces, whether periodic or not, with a non-uniform distribution of heights. In a second stage, we will attempt to reproduce existing surfaces while respecting not only their height and spatial parameters, but also their visual appearance.

### 3.1. Generating rough surfaces from scratch

Anisotropic surfaces are defined as those whose topography exhibits a preferred orientation. The most obvious case is that of unidirectionally striated surfaces. When the direction is aligned with one of the principal axes *x* or *y*, the periodicity condition inherent to the use of FFT is easily satisfied. Conversely, an orientation different from $0[k\pi/2]$ may lead to a break in periodicity along *x* or *y*.

It should be recalled that applying a digital filter to white noise – or approximately white noise – performs a convolution over the entire surface, which cannot yield localized, regular topographic patterns as shown in Figure 11. Achieving such patterns would instead require using a texture as input rather than a white-noise surface.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Localized scratches that can not result from a convolution operator applied to homogeneous noise.</figcaption></figure>

#### 3.1.1. Periodic surfaces

The rough surface to be generated has the following characteristics: *w* = 200 µm, *h* = 100 µm, *µ* = 0, *σ* = 1, *Ssk* = -3, *Sku* = 15, *α* = 40, *l*<sub>ζ</sub> = 30 µm, *l*<sub>η</sub> = 10 µm, *z* = 0.5. It is immediately noticeable that the correlation lengths are large compared to the surface size, which suggests potential issues with realism in the generation process.
As seen in Figure 12, even though all parameters are strictly satisfied, Figure 13, the result is unsatisfactory. The reason is that the theoretical autocorrelation function is not periodic – see Figure 14 – which induces spurious striations.

To overcome this difficulty, the autocorrelation function can be made periodic by applying a modified Tukey apodization function. The major advantage of the Tukey window lies in its flexibility: it provides an adjustable compromise between the frequency performance of the rectangular window – good resolution – and the strong attenuation of the Hann window – effective suppression of spurious high frequencies.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Generated rough surface showing the impact of the *acf* non-periodicity.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Autocorrelation functions along principal axes $\xi$ and $\eta$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Power Spectral Density of the dental surface that emphasizes the *acf* non-periodicity. Scales are irrevelant here, the information is qualitative.</figcaption></figure>

Here, the goal is to preserve both the surface orientation and the anisotropy factor in the windowing process. The modification of the Tukey window therefore consists of a product of cosines along the principal axes $(\xi, \eta)$ of the *acf*. More specifically, if we define an ellipse whose principal axes coincide with those of the *acf*, and whose semi-axes are $l_{\xi}$ and $l_{\eta}$, then for every point inside the ellipse the *acf* remains unchanged, while outside it is attenuated according to the Tukey window, as shown on Figure 15.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Autocorrelation function in 3D view, showing the effect of the modified Tukey window.</figcaption></figure>

Setting the boundaries of the *acf* to zero enforces periodicity, and because the transitions are smooth, the unwanted striations become barely perceptible. As a result, the generated surface appears more realistic and a way to further enhance the realism of the generated surface is to initiate the iterative convolution/height-replacement process with a very smooth surface containing very few valleys. The resulting surface exhibits fewer scattered pits, Figure 16. It is worth noting that the statistical parameters and the autocorrelation function – for distances shorter than the correlation lengths – remain strictly identical.
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Generated surface after *acf* apodization and a smooth iterative process initiation.</figcaption></figure>

### 3.1.2. Non-periodic surfaces

The use of FFT in the image generation process has the advantage of computational speed. The drawback, however, is the periodicity of the resulting surfaces. A relatively simple way to circumvent this issue is to generate a surface larger than the desired size – taking into account the correlation lengths – and then extract a sub-surface of the correct dimensions whose spatial parameters are closest to the imposed values.

In Figure 17, it can be observed that the realism of the surface has not been affected and that, clearly, the surface is no longer periodic. The statistical parameters strictly match the imposed values, and as shown in Figure 18, the autocorrelation function also remains consistent with the theoretical one. It should be noted, however, that discrepancies between the obtained *acf* and the imposed *acf* may appear when moving significantly away from the center –beyond the correlation lengths.

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Generated non-periodic surface.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Autocorrelation functions of the non-periodic surface.</figcaption></figure>

## 3.2. Reproducing real rough surfaces

At the scale of a few tens of microns, surfaces resulting from industrial material removal processes often appear rather homogeneous, even in the presence of anisotropy. Reproducing a fairly homogeneous surface, even an anisotropic one, is not a challenge in itself, although, as we will see, fine scratches require specific handling. Across all the test surfaces shown in the Figure 10, only the rotor surface is not homogeneous at the scale of correlation lengths, which represents a challenge when using convolutions. Indeed, the use of convolution in the surface generation process tends to homogenize the resulting topography, but progressive wear can create local flat areas.

The dental surface is preferentially used to test the capabilities of the developed algorithm, since it is the most difficult to faithfully reproduce. However, the other surfaces will also be used, as each has its own specific features. The steps of the algorithm are almost identical to those used for generating surfaces from the parameters *Ssk*, *Sku*, *l*<sub>ξ</sub>, and *l*<sub>η</sub>. The essential difference is that the height dataset is retrieved along with the autocorrelation function *acf*.

The direct consequence of restoring the original heights and reinjecting them into the final surface at the last step is that all height parameters are strictly identical between the original and the reproduced surfaces. Therefore, it is unnecessary to compare them before and after reconstruction. Regarding the *acf*, the situation is slightly different, as the reinjection of the original heights may alter it. Nevertheless, since the algorithm iteratively enforces compliance with the condition *acf*(*x*, *y*) ≤ 0.5, the resulting *acf* remains highly consistent with the original one. Other types of parameters are also worth measuring. The parameter *Sh* is a simple yet non-standard quantity, which represents the percentage of surface facets whose normals lie within a cone of 5° half-angle at the apex. It provides an intuitive indicator of the surface flatness. A comparison that appears particularly relevant concerns the fractal parameters – as designated in ISO 25178 – as they account for both the surface heights and their spatial arrangement. We have selected two such parameters: *Sdr* and *Sasfc*. The first corresponds to the developed interfacial area ratio, or more precisely, the excess relative area. The rougher the surface, the higher the *Sdr* value. The second parameter follows the same principle, but the calculation is performed across the different wavelength scales of the surface. The parameter *Sasfc* is commonly regarded as a measure of surface complexity, [1].
### 3.2.1. Basic surface reproduction

The idea is simply to follow, in broad terms, the approach illustrated in Figure 1. The difference is that the heights do not need to be generated, since those from the original surface are used, along with its autocorrelation function. Consequently, Block 1 is replaced by the retrieval of the heights and the autocorrelation function. Following this approach, all requirements are met: the statistical moments are obviously identical and the autocorrelation function is also preserved, as shown in 19.

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Autocorrelation functions along $\xi$ and $\eta$ for the dental surface.</figcaption></figure>

The main limitation of this method, in its raw form, is primarily visual, as it can be clearly seen in Figure 20. The unwanted folds parallel to $x$ and $y$ are caused by the autocorrelation function. Indeed, this function is obtained through FFT, which inherently assumes periodicity of the surface, whereas no such periodicity exists naturally. The consequence is a discontinuity at the edges of the surface, manifesting as these folds. This phenomenon is referred to as spectral leakage.

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Basic reproduction of the dental surface showing spectral leakage.</figcaption></figure>

A first improvement therefore consists in making the surface periodic for the computation of the autocorrelation function.
### 3.2.2. Windowing the surface to be reproduced

The simplest method is the addition of a windowing function – or apodization – which consists of attenuating the heights progressively as one approaches the edges. This is not without consequences on the height spectrum, but given the poor quality of the generated surfaces, the benefit/risk balance leans in favor of its use. It should be noted that apodization is applied only for the determination of the autocorrelation function.

Table 4 clearly illustrates the effect of windowing on the original *acf*, as the correlation lengths, the texture aspect ratio, and the main orientation are altered. By contrast, the original height distribution is preserved.

<table><caption>**Table 4.** Changes in spatial parameters with and without surface windowing. Sal: *acf* fastest decay rate - Str: Texture aspect ratio - Std: Texture direction.</caption><thead><tr><th>Surface</th><th>Windowing</th><th>Sal (μm)</th><th>Str</th><th>Std (°)</th></tr></thead><tbody><tr><td>Original</td><td>-</td><td>5.22</td><td>7.00 × 10<sup>-1</sup></td><td>-13.0</td></tr><tr><td>Reproduced, figure 20</td><td>no</td><td>5.23</td><td>7.00 × 10<sup>-1</sup></td><td>-13.0</td></tr><tr><td>Reproduced, figure 21</td><td>yes</td><td>5.87</td><td>5.58 × 10<sup>-1</sup></td><td>-19.0</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>**Figure 21.** Reproduction of the dental surface after apodization, showing spurious pits (white circles).</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure 22.** Autocorrelation functions along $\xi$ and $\eta$ after apodization, showing its limited impact.</figcaption></figure>

Figure 21 shows the result of the generation along with the new autocorrelation functions along $\xi$ and $\eta$, Figure 22. The correlation lengths are larger than previously, since the surface gradually decays to zero near the edges. All imposed parameters are satisfied, but the visual appearance of the generated surface is not yet satisfactory. In particular, one observes pits scattered across the entire surface – highlighted with white circles – that are not present in the original surface.

The idea to eliminate these spurious pits is to initialize the iterative process with a much smoother surface than the one generated randomly, Block2 in Figure 1.

### 3.2.3. Introducing a smooth topography

Depending on the type of surface to be reproduced, one may, for example, create an initial periodic surface composed of valleys and hills, Figure 23, distributed randomly.

<figure><img src="image_25.png" /><figcaption>**Figure 23.** Smooth topography for the iterative process initiation.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure 24.** Reproduction of the dental surface with the smoothing texture.</figcaption></figure>

The result is clear: the spurious pits disappear, Figure 24, while all imposed parameters are still perfectly satisfied. Typically, the random surface used before applying the digital filter should be white noise, so that after convolution the autocorrelation function is strictly identical to the imposed autocorrelation function. This statement, however, must be qualified: when the statistical parameters *Ssk* and *Sku* are large, the noise is no longer white. An example of random surfaces departing from the characteristics of white noise is provided in Appendix C. Consequently, a few iterations are required to obtain the correct autocorrelation function. Based on this principle, one can deliberately deviate from white noise for the surface to be convolved, in order to enforce a smoother final surface.

When the original and reproduced surfaces are visually compared, it can be observed that the numerically generated surface, although highly consistent with the imposed parameters, exhibits a higher degree of complexity than the original one. The *Sh* parameter indicates that the surface is globally less flat, while the increased value of *Sasfc* reflects a greater rise in surface area with increasing sampling resolution, Table 5.

<table><caption>**Table 5.** Comparison of hybrid and spatial parameters between the original and the reproduced surface – dental surface case.</caption><thead><tr><th>Surface</th><th>Sh (%)</th><th>Sdr</th><th>Sasfc</th><th>Sal (μm)</th><th>Str</th><th>Std (°)</th></tr></thead><tbody><tr><td>Original</td><td>4.30 × 10<sup>1</sup></td><td>1.70 × 10<sup>-2</sup></td><td>2.29</td><td>5.87</td><td>5.57 × 10<sup>-1</sup></td><td>19.0</td></tr><tr><td>Reproduced</td><td>2.39 × 10<sup>1</sup></td><td>2.08 × 10<sup>-2</sup></td><td>2.61</td><td>5.87</td><td>5.58 × 10<sup>-1</sup></td><td>19.0</td></tr></tbody></table>

Most machined surfaces can be reproduced following this principle. The milled surface, which exhibits strong anisotropy, is faithfully reproduced not only in terms of statistical and spatial parameters but also from a purely visual standpoint, Figure 25.
<figure><img src="image_27.png" /><figcaption>**Figure 25.** Reproduction of the milled surface.</figcaption></figure>

Nevertheless, a quantitative analysis of the results reveals the same discrepancies in overall flatness and complexity, Table 6 A closer match to the original surface could be achieved by performing local height permutations, as proposed by Chen *et al.* [8], which would increase *Sh* and decrease *Sasfc*.

<table><caption>**Table 6.** Comparison of hybrid and spatial parameters between the original and the reproduced surface - milled surface case.</caption><thead><tr><th>Surface</th><th>Sh (%)</th><th>Sdr</th><th>Sasfc</th><th>Sal (μm)</th><th>Str</th><th>Std (°)</th></tr></thead><tbody><tr><td>Original</td><td>2.16 × 10<sup>1</sup></td><td>1.12 × 10<sup>-1</sup></td><td>4.45 × 10<sup>1</sup></td><td>2.40 × 10<sup>1</sup></td><td>1.01 × 10<sup>-1</sup></td><td>6.0</td></tr><tr><td>Reproduced</td><td>1.25 × 10<sup>1</sup></td><td>1.46 × 10<sup>-1</sup></td><td>5.80 × 10<sup>1</sup></td><td>2.35 × 10<sup>1</sup></td><td>9.50 × 10<sup>-2</sup></td><td>6.0</td></tr></tbody></table>

The ground surface, which displays similar characteristics, also yields very good results, Figure 26.

<figure><img src="image_28.png" /><figcaption>**Figure 26.** Reproduction of the ground surface.</figcaption></figure>

In this case as well, the surface reproduction procedure exhibits the same behavior, characterized by reduced flatness and increased complexity, Table 7.
<table><caption>**Table 7.** Comparison of hybrid and spatial parameters between the original and the reproduced surface – ground surface case.</caption><thead><tr><th>Surface</th><th>Sh (%)</th><th>Sdr</th><th>Sasfc</th><th>Sal (μm)</th><th>Str</th><th>Std (°)</th></tr></thead><tbody><tr><td>Original</td><td>3.43 × 10<sup>1</sup></td><td>9.63 × 10<sup>-2</sup></td><td>1.75 × 10<sup>1</sup></td><td>5.40</td><td>4.99 × 10<sup>-2</sup></td><td>47.0</td></tr><tr><td>Reproduced</td><td>1.12 × 10<sup>1</sup></td><td>1.06 × 10<sup>-1</sup></td><td>2.00 × 10<sup>1</sup></td><td>5.38</td><td>4.95 × 10<sup>-2</sup></td><td>47.0</td></tr></tbody></table>

When a surface exhibits highly localized features, the qualitative appearance is somewhat degraded, as can be seen on the polished surface, Figure 27. Nevertheless, the outcome remains highly satisfactory.

<figure><img src="image_29.png" /><figcaption>**Figure 27.** Reproduction of the polished surface, showing a lack of high and flat areas, resulting from the polishing process.</figcaption></figure>

It is noteworthy that, in this instance, the situation differs slightly. While Sh provides little meaningful information due to the inherent flatness of the surface, the overall complexity appears to be well preserved, Table 8. A more relevant assessment in this context would involve considering the distribution of topographic features; however, such an analysis falls outside the scope of the present study.

<table><caption>**Table 8.** Comparison of hybrid and spatial parameters between the original and the reproduced surface – polished surface case.</caption><thead><tr><th>Surface</th><th>Sh (%)</th><th>Sdr</th><th>Sasfc</th><th>Sal (μm)</th><th>Str</th><th>Std (°)</th></tr></thead><tbody><tr><td>Original</td><td>1.00 × 10<sup>2</sup></td><td>2.12 × 10<sup>-6</sup></td><td>2.81 × 10<sup>-4</sup></td><td>2.66 × 10<sup>1</sup></td><td>4.89 × 10<sup>-1</sup></td><td>-10.0</td></tr><tr><td>Reproduced</td><td>1.00 × 10<sup>2</sup></td><td>2.15 × 10<sup>-6</sup></td><td>2.94 × 10<sup>-4</sup></td><td>2.63 × 10<sup>1</sup></td><td>4.90 × 10<sup>-1</sup></td><td>-10.0</td></tr></tbody></table>

Another case, comparable to the dental surface, concerns the rotor surface. Since the contact area is relatively smooth, it is now well known that applying a convolution cannot reproduce the surface with high fidelity: the resulting surface is clearly less smooth than the original, Figure 28, even though the prescribed parameters are strictly satisfied.
<figure><img src="image_30.png" /><figcaption>**Figure 28.** Reproduction of the rotor worn surface. Here again, flat areas should be observed.</figcaption></figure>

The situation here is broadly similar to that observed previously, with no significant deviations in complexity. What differentiates the generated surface from the original one is primarily a matter of homogeneity in the height distribution, Table 9.

<table><caption>**Table 9.** Comparison of hybrid and spatial parameters between the original and the reproduced surface – rotor surface case.</caption><thead><tr><th>Surface</th><th>Sh (%)</th><th>Sdr</th><th>Sasfc</th><th>Sal (μm)</th><th>Str</th><th>Std (°)</th></tr></thead><tbody><tr><td>Original</td><td>9.59 × 10<sup>1</sup></td><td>2.64 × 10<sup>-3</sup></td><td>4.37 × 10<sup>-1</sup></td><td>4.10</td><td>7.75 × 10<sup>-1</sup></td><td>3.0</td></tr><tr><td>Reproduced</td><td>9.31 × 10<sup>1</sup></td><td>2.92 × 10<sup>-3</sup></td><td>4.80 × 10<sup>-1</sup></td><td>4.05</td><td>7.74 × 10<sup>-1</sup></td><td>3.0</td></tr></tbody></table>

### 3.2.4. Using a low-pass filter

Another way to achieve a similar result is to apply a low-pass filter – such as a top-hat filter – immediately before each convolution step. This forces the process to limit high frequencies, which also leads to a smoother surface, Figure 29.

<figure><img src="image_31.png" /><figcaption>**Figure 29.** Reproduction of the dental surface, after top-hat low-pass filter.</figcaption></figure>
It should be emphasized that no filter is applied after the final convolution iteration, since this would degrade the autocorrelation function. In the example presented, the cutoff frequency is set to 0.01, i.e., 100 times the grid spacing, or 12.9 µm.

In order to obtain a final surface that is as smooth as possible, the top-hat filter must remove as many high frequencies as possible. However, it can be observed that in the iterative convolution/height-replacement procedure, removing too many high frequencies prior to convolution slows down the process. In the present case, when the cutoff falls below 0.01, the procedure no longer converges satisfactorily, and a larger deviation in the acf would therefore have to be accepted, without any apparent improvement in the final result.

As the iterations proceed, the effect of the top-hat filter remains minor, as can be seen in Table 10 when this filter is applied in the final step. In all cases, the procedure concludes with the convolution and the restoration of the original height values, which ensures that both the height-related and spatial parameters are correctly preserved. One may be surprised by the limited impact of the filter, given the visual appearance of the surface, as all wavelengths shorter than 100 mesh steps are absent. This observation is, in fact, simply a consequence of the surface characteristics, particularly its correlation lengths of approximately 50 and 80 times the mesh step.

<table><caption>**Table 10.** Changes in spatial and height parameters after top-hat smoothing.</caption><thead><tr><th>Windowing</th><th>Sv (µm)</th><th>Sp (µm)</th><th>Smd (µm)</th><th>Sa (µm)</th><th>Sq (µm)</th><th>Ssk</th><th>Sku</th><th>Sal (µm)</th><th>Str</th><th>Std (°)</th></tr></thead><tbody><tr><td>no</td><td>-2.31</td><td>1.25</td><td>6.84 × 10<sup>-2</sup></td><td>3.46 × 10<sup>-1</sup></td><td>4.85 × 10<sup>-1</sup></td><td>-1.20</td><td>6.02</td><td>5.86</td><td>5.58 × 10<sup>-1</sup></td><td>19</td></tr><tr><td>yes</td><td>-2.14</td><td>1.14</td><td>7.23 × 10<sup>-2</sup></td><td>3.45 × 10<sup>-1</sup></td><td>4.85 × 10<sup>-1</sup></td><td>-1.27</td><td>6.07</td><td>6.12</td><td>5.64 × 10<sup>-1</sup></td><td>19</td></tr></tbody></table>

Although the result is quite satisfactory, the generated surface, while matching the statistical parameters of the original one, still does not sufficiently "resemble" it in appearance.

### 3.2.5. Two-level reconstruction

To more faithfully reproduce the surface across its different frequency components, one can attempt to decompose it into two levels: a low-frequency level, Figure 30 and a high-frequency level, Figure 31.
<figure><img src="image_32.png" /><figcaption>**Figure 30.** Dental surface: low frequencies (LF).</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Figure 31.** Dental surface: high frequencies (HF).</figcaption></figure>

Note that there is no universally accepted definition to distinguish low from high frequencies. In this work, the cutoff wavelength was set to 200 grid steps, i.e., 25.8 µm.

The modified algorithm consists of reproducing as accurately as possible the two surfaces corresponding to the frequency-partitioned components of the original surface – while preserving the four statistical moments and the autocorrelation function – and then recombining them with a few final iterations so that the resulting surface is fully compliant. The resulting surfaces are shown in Figure 32 and Figure 33, and the final result in Figure 34. Once again, all imposed parameters are satisfied, but differences remain observable.
<figure><img src="image_34.png" /><figcaption>**Figure 32.** Reproduction of the LF dental surface.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>**Figure 33.** Reproduction of the HF dental surface.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>**Figure 34.** Final result after combination of BF and HF results.</figcaption></figure>

The decomposition did not bring any substantial improvement regarding the "skin-like" aspect of the original surface. Moreover, the parameters *Sh*, *Sdr*, and *Sasfc* are nearly
identical to those obtained for the surface generated by introducing a smooth surface. 588
Nevertheless, one could consider adding a few scratches to the numerical surface, which 589
would provide greater realism. 590

3.2.6. Introducing scratches

By introducing into the program a routine that models scratches of random lengths 591
and orientations, a scratched texture can be imposed at the beginning of the iterations, 592
causing them to appear in the final surface. It should be noted that the digital filter is 593
then applied to a scratched surface, which ultimately smooths the scratches while leaving 594
scattered artifacts, Figure 35. 595
596

<figure><img src="image_37.png" /><figcaption>**Figure 35.** Reproduction of the dental surface with added scratches.</figcaption></figure>

The use of a scratched surface can be particularly useful for reproducing polished 597
mirror-like industrial surfaces, since, as shown in Figure 36, the result is satisfactory. The 598
quality cannot truly be improved by reducing the thickness of the scratches, since even 599
when starting from scratches as narrow as a single pixel, the convolution process broadens 600
them, resulting in scratches that remain wider than desired. 601

<figure><img src="image_38.png" /><figcaption>**Figure 36.** Reproduction of a polished mirror-like industrial surfaces.</figcaption></figure>
In contrast to the other surfaces, the generated surface exhibits lower complexity, Table 11. This can be explained by the fact that the original scratches display more varied characteristics than those introduced during the generation process.

<table><caption>**Table 11.** Comparison of hybrid and spatial parameters between the original and the reproduced surface – polished mirror-like surface case.</caption><thead><tr><th>Surface</th><th>Sh (%)</th><th>Sdr</th><th>Sasfc</th><th>Sal (μm)</th><th>Str</th><th>Std (°)</th></tr></thead><tbody><tr><td>Original</td><td>1.00 × 10<sup>2</sup></td><td>3.10 × 10<sup>-6</sup></td><td>6.48 × 10<sup>-4</sup></td><td>9.38</td><td>6.67 × 10<sup>-1</sup></td><td>-34.0</td></tr><tr><td>Reproduced</td><td>1.00 × 10<sup>2</sup></td><td>2.92 × 10<sup>-6</sup></td><td>6.12 × 10<sup>-4</sup></td><td>9.07</td><td>6.67 × 10<sup>-1</sup></td><td>-34.0</td></tr></tbody></table>

## 4. Conclusions

With regard to the generation of random rough surfaces with prescribed parameters, the algorithm presented has demonstrated the ability to achieve the following objectives:

* preservation of strong directional anisotropies,

* accurate reproduction of large correlation lengths,

* enforcement of high third- and fourth-order statistical moments,

* near-realistic rendering of the resulting surfaces.

Concerning the reproduction of rough surfaces, the contributions are as follows:

* automatic compliance with both height and spatial parameters,

* preservation of realism for surfaces that remain homogeneous at the scale of the correlation lengths,

* approximate reproduction of localized topographic features.

The advantages of the proposed process over the REPR (Representative Elementary Pattern of Roughness) method developed by Borodich *et al.* [11] are as follows:

* the generated surfaces are, in theory, not limited in size;

* the statistical parameters of the original surface are preserved, whatever the chosen size;

* heterogeneity can be introduced into the surface;

* the spatial arrangement of heights through the PSD, although imperfect, allows the generated surface to approximate the original one without being a direct extraction of it.

This last point is crucial, particularly in lubrication-related problems, since the wavelengths of the surface heights do not contribute equally to load-carrying capacity. With the REPR method, the extracted surface patch may exhibit a wavelength distribution that differs significantly from that of the original surface, and thus fail to reproduce its tribological behavior accurately. Finally, by construction, the REPR method cannot provide diversity in the generated surfaces, which limits the statistical analysis of surface parameters with respect to a given quantity (such as load-carrying capacity or efficiency coefficient).

The attempt to decompose rough surfaces into low- and high-frequency components did not yield the expected results. Nevertheless, we believe this approach should be further explored through a multi-scale decomposition that is not limited to only two components. During the recomposition phase, scaling factors and offsets could then be adjusted for each surface component, so as to better match both the imposed parameters and the
visual rendering. Working with a multi-scale decomposition also has the advantage of providing explicit control over:

* low-frequency components, which are most relevant for lubrication-related issues, and

* high-frequency components, which capture micro-wear patterns associated with diet.

If this decomposition principle proves to be effective, it will then become possible to generate large quantities of rough surfaces for statistical studies, and even for training within the framework of machine learning.

**Author Contributions:** Conceptualization, methodology, software, A.Francisco; validation, N.Brunetière and A.Francisco; writing—original draft preparation, A.Francisco; writing—review and editing, N.Brunetière and A.Francisco. All authors have read and agreed to the published version of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Nomenclature

638

639

640

641

642

643

644

645

647

648

649

650

651
a, b, c Analytical function parameters

acf Autocorrelation function

h Filter (Hu and Tonder)

l_ζ, l_η Autocorrelation lengths along principal axes ζ and η, at z = 0.5.
The default value for z is z = 0.2, however for large acf lengths,
a higher z value should be selected.

n Number of heights

va Surface height variance

w, h Surface dimensions

x, y Horizontal and vertical axes of a surface

x_i z_i x-coordinate (analytical height function)

z Height random variable; acf plane cut

z_i Surface height number i

Oz z permutation operator for ordering

Sa ISO 25178 parameter, arithmetical mean height of the surface

Sasfc ISO 25178 parameter, area-scale fractal analysis complexity pa-
rameter

Sal ISO 25178 parameter corresponding to l_η

Sdr ISO 25178 parameter, excess relative area

Sh Percentage of surface facets whose normals lie within a cone of 5° half-angle at the apex

Smd Height median

Sp ISO 25178 parameter, maximum height of peaks

Sq ISO 25178 parameter, root mean square height of the surface

Ssk, Sku ISO 25178 parameter, height skewness and kurtosis

Std ISO 25178 parameter, texture direction of the surface

Stri Inverse of the ISO 25178 anisotropy factor Str

Sv ISO 25178 parameter, maximum height of valleys

α Texture orientation (deg)

μ Surface height mean

σ Surface height standard deviation

**Appendix A. Mathematical justification of the statistical moments of z = h * η**

652

The following results hold only for infinite series and for white noise η with zero mean and unit variance.

653

654

655

Appendix A.1. Mean $\mu_z$

656

$$
\frac{1}{n} \sum_{i=1}^{n} z_i = \sum_{k=1}^{n} h_k \left( \frac{1}{n} \sum_{i=1}^{n} \eta_{k+i} \right) = n \mu_h \mu_\eta
$$

as $\mu_{\eta} = 0, \mu_{z} = 0$.

657

Appendix A.2. Variance $\sigma_z^2$

658

z being zero-mean :

659
$$ \sigma_z^2 = \frac{1}{n} \sum_{i=1}^{n} z_i^2 = \sum_{k=1}^{n} h_k^2 \left( \frac{1}{n} \sum_{i=1}^{n} \eta_{k+i}^2 \right) + \frac{2}{n} \sum_{k=1}^{n-1} \sum_{p=k+1}^{n} h_k h_p \sum_{i=1}^{n} \eta_{k+i} \eta_{p+i} $$

$\eta$ being unit variance and considered as white noise:

$$ \sum_{i=1}^{n} \eta_{k+i} \eta_{p+i} = 0 \quad k \neq p $$

therefore:

$$ \sigma_z^2 = \sum_{k=1}^{n} h_k^2 \cdot \sigma_\eta^2 = \sum_{k=1}^{n} h_k^2 = n \sigma_h^2 $$

Appendix A.3. Skewness $Ssk_z$

$$ \frac{1}{n} \sum_{i=1}^{n} z_i^3 = \sum_{k=1}^{n} h_k^3 \left( \frac{1}{n} \sum_{i=1}^{n} \eta_{k+i}^3 \right) \\ + \frac{3}{n} \sum_{k=1}^{n-1} \sum_{p=k+1}^{n} h_k^2 h_p \sum_{i=1}^{n} \eta_{k+i}^2 \eta_{p+i} + \\ + \frac{6}{n} \sum_{k=1}^{n-2} \sum_{p=k+1}^{n-1} \sum_{q=p+1}^{n} h_k h_p h_q \sum_{i=1}^{n} \eta_i \eta_{(p-k)+i} \eta_{(q-k)+i} $$

Under the key assumption that the $\eta_i$ are independent, identically distributed:

$$ E[\eta_{k+i}^2 \eta_{p+i}] = E[\eta_{k+i}^2] \cdot E[\eta_{p+i}] = \sigma_\eta^2 \cdot 0 = 0 $$

For the same reason:

$$ \sum_{i=1}^{n} \eta_i \eta_{(p-k)+i} \eta_{(q-k)+i} = 0 $$

Therefore:

$$ Ssk_z = \frac{1}{n\sigma_z^3} \sum_{i=1}^{n} z_i^3 = \frac{1}{\sqrt{n}\sigma_h^3} \sum_{k=1}^{n} h_k^3 Ssk_\eta $$

and as $\mu_h = 0$,

$$ Ssk_z = \frac{1}{\sqrt{n}} Ssk_h Ssk_\eta $$
Appendix A.4. Kurtosis Sku<sub>z</sub>

667

$$
\begin{aligned}
\frac{1}{n} \sum_{i=1}^{n} z_{i}^{4} = & \sum_{k=1}^{n} h_{k}^{4} \left( \frac{1}{n} \sum_{i=1}^{n} \eta_{k+i}^{4} \right) \\
& + \frac{4}{n} \sum_{k=1}^{n-1} \sum_{p=k+1}^{n} h_{k}^{3} h_{p} \sum_{i=1}^{n} \eta_{k+i}^{3} \eta_{p+i} \\
& + \frac{6}{n} \sum_{k=1}^{n-1} \sum_{p=k+1}^{n} h_{k}^{2} h_{p}^{2} \sum_{i=1}^{n} \eta_{k+i}^{2} \eta_{p+i}^{2} \\
& + \frac{12}{n} \sum_{k=1}^{n-2} \sum_{p=k+1}^{n-1} \sum_{q=p+1}^{n} h_{k}^{2} h_{p} h_{q} \sum_{i=1}^{n} \eta_{k+i}^{2} \eta_{p+i} \eta_{q+i} \\
& + \frac{24}{n} \sum_{k=1}^{n-3} \sum_{p=k+1}^{n-2} \sum_{q=p+1}^{n-1} \sum_{r=q+1}^{n} h_{k} h_{p} h_{q} h_{r} \sum_{i=1}^{n} \eta_{i} \eta_{(p-k)+i} \eta_{(q-k)+i} \eta_{(r-k)+i}
\end{aligned}
$$

In the same way as we have eliminated certain products earlier by invoking that $\eta$ is white noise, we can eliminate several terms here as well. What remains is:

668

669

$$
\frac{1}{n} \sum_{i=1}^{n} z_{i}^{4} = \sum_{k=1}^{n} h_{k}^{4} \left( \frac{1}{n} \sum_{i=1}^{n} \eta_{k+i}^{4} \right) + \frac{6}{n} \sum_{k=1}^{n-1} \sum_{p=k+1}^{n} h_{k}^{2} h_{p}^{2} \sum_{i=1}^{n} \eta_{k+i}^{2} \eta_{p+i}^{2}
$$

Under the same assumptions as before about $\eta$, we can state ([24]):

670

$$
\frac{1}{n} \sum_{i=1}^{n} \eta_{k+i}^{2} \eta_{p+i}^{2} = \frac{1}{n} \sum_{i=1}^{n} (\eta_{k+i} \eta_{p+i})^{2} = \sigma_{\eta}^{2} \sigma_{\eta}^{2} = 1
$$

hence:

671

$$
\frac{1}{n} \sum_{i=1}^{n} z_{i}^{4} = \sum_{k=1}^{n} h_{k}^{4} \left( \frac{1}{n} \sum_{i=1}^{n} \eta_{k+i}^{4} \right) + 6 \sum_{k=1}^{n-1} \sum_{p=k+1}^{n} h_{k}^{2} h_{p}^{2}
$$

It can be easily established that:

672

$$
2 \sum_{k=1}^{n-1} \sum_{p=k+1}^{n} h_{k}^{2} h_{p}^{2} = \left( \sum_{k=1}^{n} h_{k}^{2} \right)^{2} - \sum_{k=1}^{n} h_{k}^{4}
$$

hence

673

$$
\begin{aligned}
\frac{1}{n} \sum_{i=1}^{n} z_{i}^{4} &= n\sigma_{h}^{4} Sku_{h} Sku_{\eta} + 3 \left[ (n\sigma_{h}^{2})^{2} - n\sigma_{h}^{4} Sku_{h} \right] \\
&= n^{2}\sigma_{h}^{4} \left[ \frac{1}{n} Sku_{h} Sku_{\eta} + 3 \left( 1 - \frac{1}{n} Sku_{h} \right) \right]
\end{aligned}
$$

Finally,

674
$$ Sku_z - 3 = \frac{1}{n} Sku_h (Sku_\eta - 3) $$

## Appendix B. Some analytical justifications regarding the boundaries

The maximum value of the pair $(Ssk, Sku)$ is approximately $(-\sqrt{n}, n)$, so the function $Sku = 1.24 Ssk^2 + 1.8$ is not compatible with this point. It turns out that for $Ssk < -\frac{\sqrt{n}}{2}$, the points align along a straight line, whereas for $Ssk > -\frac{\sqrt{n}}{2}$, the parabola $Sku = 1.24 Ssk^2 + 1.8$ fits perfectly.

This explains why the junction of the two functions occurs at $Ssk \approx -\frac{\sqrt{n}}{2}$, when determining the parameters of the straight line passing through the point $(-\sqrt{n}, n)$ and tangent to the parabola $Sku = 1.24 Ssk^2 + 1.8$.

Moreover, for the transition, the intersection point of the two curves is exactly $(-\frac{\sqrt{n}}{2}, \frac{\sqrt{n}}{3})$, to first order in $\frac{1}{\sqrt{n}}$, if one approximates 1.34 by $\frac{4}{3}$ and replaces 1.8 with $\frac{9}{5}$.

### Appendix B.1. Minimum of Sku when $Ssk = 0$

#### Appendix B.1.1. Tangent based distribution

For the tangent-based distribution, the lowest value of $Sku$ is reached around 0. Indeed, near $x = 0$, $\tan(x) \approx x$, which is the region where all dataset points $\tan(x_i)$ are closest to each other. We then have:

$$ \begin{aligned} \sigma^2 &= \frac{1}{2\epsilon} \int_{-\epsilon}^{+\epsilon} (x - \mu)^2 dx = \frac{1}{2\epsilon} \int_{-\epsilon}^{+\epsilon} x^2 dx \\ &= \frac{1}{3}\epsilon^2 \end{aligned} $$

$\mu = 0$ because the interval is symmetric and the function $f(x) = x$ is odd.

$$ \begin{aligned} Sku &= \frac{1}{2\epsilon} \int_{-\epsilon}^{+\epsilon} \left(\frac{x-\mu}{\sigma}\right)^4 dx \\ &= \frac{1}{2\epsilon\sigma^4} \int_{-\epsilon}^{+\epsilon} x^4 dx \\ &= \frac{9}{2\epsilon^5} \frac{2\epsilon^5}{5} = \frac{9}{5} \end{aligned} $$

#### Appendix B.1.2. Modified exponential based distribution

The smallest attainable $Sku$ lies on Pearson's bound, namely $Sku = Ssk^2 + 1$. This limit is reached for binary distributions. The exponential-based distribution can approach this limit when the $x^3$ component is canceled out and the integration bounds are "far away".
$$ \mu = \lim_{l \to \infty} \left( \frac{1}{2l} \int_{-l}^{0} (-1 + \exp(x)) dx + \frac{1}{2l} \int_{0}^{l} (+1 - \exp(-x)) dx \right) = 0 $$

$$ \sigma^2 = \lim_{l \to \infty} \left( \frac{1}{2l} \int_{-l}^{0} (-1 + \exp(x))^2 dx + \frac{1}{2l} \int_{0}^{l} (+1 - \exp(-x))^2 dx \right) \\ = \lim_{l \to \infty} \left( \frac{1}{2l} (-3 + 2l + 4 \exp(-l) - \exp(-2l)) \right) = 1 $$

$$ \begin{aligned} Sku &= \lim_{l \to \infty} \left( \frac{1}{2l} \int_{-l}^{0} (-1 + \exp(x))^4 dx + \frac{1}{2l} \int_{0}^{l} (+1 - \exp(-x))^4 dx \right) \\ &= \lim_{l \to \infty} \left( -\frac{25}{6} + 2l - 8 \exp(-l) - 6 \exp(-2l) + \frac{8}{3} \exp(-3l) - \frac{1}{2} \exp(-4l) \right) = 1 \end{aligned} $$

## Appendix B.2. Maximum of Sku when Ssk = 0

### Appendix B.2.1. Tangent based distribution

For the tangent-based distribution, the analysis is complex since it involves integrating the fourth power of the tangent function close to its vertical asymptotes located at $\pm \frac{\pi}{2}$, within $\epsilon$. The analytical expression can be obtained using a symbolic computation tool such as Maxima [25]

$$ Sku = - \frac{\pi(\epsilon + 1)(-3\pi(1-\epsilon)s^3 - 8c^3 + 6c)}{3s((\pi^2(1-\epsilon)^2 - 4)s^2 - 4\pi(1-\epsilon)cs + 4)} \sim \frac{1}{3\epsilon} \quad (A1) $$

$$ \text{with } s = \sin(\frac{\pi\epsilon}{2}) \text{ and } c = \cos(\frac{\pi\epsilon}{2}). $$

Using 128-bit floating-point numbers, one can numerically approach the theoretical value quite closely with $\epsilon = 10^{-3}$, namely $Sku \approx 333$. Obviously, this value is far above any physically encountered value for a surface, since it would correspond to a globally flat surface with isolated peaks and valleys of very large amplitude. However, during the generation of a random series prior to the application of a digital filter, such values may be required.

The key point here is that the $Sku_{max}$ curve cannot be easily plotted since it depends on the implementation of the distribution (floating-point encoding and the value of $\epsilon$). For readers aiming to approximate the theoretical curve as closely as possible, it may be necessary to decompose the statistical moments according to intervals:

$$ \begin{aligned} \left[-\frac{\pi}{2}(1-\epsilon_l), \frac{\pi}{2}(1-\epsilon_u)\right] &= \left[-\frac{\pi}{2}(1-\epsilon_l), -\frac{\pi}{2}(1-\alpha)\right] \\ &\cup \left[-\frac{\pi}{2}(1-\alpha), \frac{\pi}{2}(1-\alpha)\right] \\ &\cup \left[\frac{\pi}{2}(1-\alpha), \frac{\pi}{2}(1-\epsilon_u)\right] \end{aligned} $$

with $(a, b) = (\epsilon_l, \epsilon_u)$ a pair of very small values, e.g. $10^{-6}$, and $\alpha = 10^{-3}$ a value for which the integral can be computed numerically without accuracy issues. The idea is then to approximate the integrals near $\pm \frac{\pi}{2}$ by series expansions in $\epsilon_l$ and $\epsilon_u$. The point here is to test the overall robustness of the analytic-numeric approach.
## Appendix B.2.2. Modified exponential based distribution

For the exponential-based distribution, the number of parameters was limited to three: two parameters for the integration bounds and one parameter for the $x^3$ component. Moreover, the choice of normalizing $x^3$ over the integration interval has the following effects:

1. it favors the binary-distribution effect when the integration bounds move far from the origin,

2. it amplifies the influence of $x^3$ near the origin.

The maximum value of $Sku$ is reached in case 2, on $[x_a, x_b] = [-\epsilon, +\epsilon]$, therefore:

$$\mu = \lim_{l \to \infty} \frac{1}{2\epsilon} \int_{-\epsilon}^{+\epsilon} \left(\frac{x}{2\epsilon}\right)^3 dx = 0$$

$$\sigma^2 = \lim_{l \to \infty} \frac{1}{2\epsilon} \int_{-\epsilon}^{+\epsilon} \left(\frac{x}{2\epsilon}\right)^6 dx = \frac{1}{7.26}$$

$$Sku = \lim_{l \to \infty} \frac{1}{2\epsilon} \int_{-\epsilon}^{+\epsilon} \left(\frac{x}{2\epsilon}\right)^{12} .7^2 .2^{12} dx = \frac{49}{13} \approx 3.77$$

In practice, the NSGA-II algorithm in Pymoo does not succeed in properly approaching this limit. The reason lies in the infinitesimals induced by the integration interval and the presence of exponential and power functions. This theoretical limit replaces the numerical value obtained, which is lower.

## Appendix B.3. Maximum of both Sku and Ssk

$$\text{i.e. } (Ssk, Sku) = \max\{(Ssk, Sku)_{(a,b)}, (a,b) \in ]0, 1]^2\}$$

### Appendix B.3.1. Tangent based distribution

For the tangent based distribution, the situation is complex, as previously explained. To obtain large values for the pair $(Ssk, Sku)$, it suffices to take the integration interval $[\frac{\pi}{2}(1-\epsilon), 0]$ with $\epsilon$ very small. Once again, the values $(|Ssk|, Sku)$ are upper-bounded by $(\sqrt{n-3}, n-2)$.

### Appendix B.3.2. Modified exponential based distribution

For the exponential-based distribution, in order to decenter and flatten the distribution as much as possible, it suffices to select a lower integration bound of zero and an upper bound $l \to \infty$, and to cancel the $x^3$ function:

$$\mu = \frac{1}{l} \int_{0}^{l} (1 - \exp(-x)) dx \approx 1$$

$$\sigma^2 = \frac{1}{l} \int_{0}^{l} (1 - \exp(-x) - 1)^2 dx \approx \frac{1}{2l}$$

$$Ssk = \frac{1}{l} \int_{0}^{l} (-\exp(-x))^3 (2l)^{\frac{3}{2}} dx \approx -\frac{2^{\frac{3}{2}}}{3} \sqrt{l}$$

$$Sku = \frac{1}{l} \int_{0}^{l} (-\exp(-x))^4 4l^2 dx \approx l$$
This choice maximizes both Ssk and Sku. When $l$ is sufficiently large, the function $f(x) = 1 - \exp(-x)$ yields for most $x_i$ a value around 1, and for $x_0 = 0$ the value 0. Hence, the statistical moments can be deduced, computed over the corresponding dataset, by taking for the two values $-\sqrt{n-1}$ and $\sqrt{n-1}^{-1}$, which ensures $\mu = 0$ and $\sigma = 1$.

$$Ssk = \frac{1}{n} \left( -\sqrt{n-1}^3 + (n-1)\sqrt{n-1}^{-3} \right) \approx -\sqrt{n-3}$$

$$Sku = \frac{1}{n} \left( \sqrt{n-1}^4 + (n-1)\sqrt{n-1}^{-4} \right) \approx n-2$$

and Pearson's inequality is naturally satisfied, since the distribution is binary. Thus, for the exponential based distribution, ($|Ssk|, Sku) \le (\sqrt{n-3}, n-2)$, a result also confirmed numerically.

## Appendix C. Impact of the statistical moments Ssk and Sku on the noise autocorrelation

Table A1 details the number of acf values of a random surface that exceed a threshold $th$, as a function of skewness (Ssk) and kurtosis (Sku). The generated surfaces contain 1024 × 1024 data points. The statistical moments satisfy the relation $Sku = 1.1 \times Ssk^2 + 1.8$, ensuring that the values remain near the Pearson boundary. It is observed that, in contrast to white noise, secondary peaks can occur.

<table><caption>**Table A1.** Number of acf values, except the center, exceeding threshold $th$ as a function of Ssk and Sku for a random surface.</caption><thead><tr><th>Ssk</th><th>Sku</th><th>th > 0.05</th><th>th > 0.1</th><th>th > 0.2</th></tr></thead><tbody><tr><td>0</td><td>1.8</td><td>0</td><td>0</td><td>0</td></tr><tr><td>-20</td><td>441.8</td><td>0</td><td>0</td><td>0</td></tr><tr><td>-40</td><td>1761.8</td><td>0</td><td>0</td><td>0</td></tr><tr><td>-60</td><td>3961.8</td><td>0</td><td>0</td><td>0</td></tr><tr><td>-80</td><td>7041.8</td><td>0</td><td>0</td><td>0</td></tr><tr><td>-100</td><td>11001.8</td><td>0</td><td>0</td><td>0</td></tr><tr><td>-120</td><td>15841.8</td><td>0</td><td>0</td><td>0</td></tr><tr><td>-140</td><td>21561.8</td><td>0</td><td>0</td><td>0</td></tr><tr><td>-160</td><td>28161.8</td><td>0</td><td>0</td><td>0</td></tr><tr><td>-180</td><td>35641.8</td><td>18</td><td>0</td><td>0</td></tr><tr><td>-200</td><td>44001.8</td><td>86</td><td>0</td><td>0</td></tr><tr><td>-220</td><td>53241.8</td><td>132</td><td>0</td><td>0</td></tr><tr><td>-240</td><td>63361.8</td><td>132</td><td>0</td><td>0</td></tr><tr><td>-260</td><td>74361.8</td><td>116</td><td>8</td><td>0</td></tr><tr><td>-280</td><td>86241.8</td><td>114</td><td>20</td><td>0</td></tr><tr><td>-300</td><td>99001.8</td><td>94</td><td>26</td><td>0</td></tr><tr><td>-320</td><td>112641.8</td><td>84</td><td>30</td><td>0</td></tr><tr><td>-340</td><td>127161.8</td><td>66</td><td>42</td><td>0</td></tr><tr><td>-360</td><td>142561.8</td><td>60</td><td>28</td><td>0</td></tr><tr><td>-380</td><td>158841.8</td><td>44</td><td>30</td><td>0</td></tr><tr><td>-400</td><td>176001.8</td><td>40</td><td>26</td><td>6</td></tr><tr><td>-420</td><td>194041.8</td><td>30</td><td>20</td><td>10</td></tr><tr><td>-440</td><td>212961.8</td><td>26</td><td>20</td><td>12</td></tr><tr><td>-460</td><td>232761.8</td><td>26</td><td>18</td><td>8</td></tr><tr><td>-480</td><td>253441.8</td><td>24</td><td>16</td><td>6</td></tr><tr><td>-500</td><td>275001.8</td><td>18</td><td>12</td><td>6</td></tr></tbody></table>

## References

1. Leach, R. Characterisation of Areal Surface Texture; Springer-Verlag Berlin and Heidelberg GmbH & Co. K: Heidelberg, 2013.

754

755
2. Wang, W.z.; Chen, H.; Hu, Y.z.; Wang, H. Effect of Surface Roughness Parameters on Mixed Lubrication Characteristics. *Tribology International* **2006**, *39*, 522–527. https://doi.org/10.1016/j.triboint.2005.03.018.

756

757

3. Sedlaček, M.; Podgornik, B.; Vižintin, J. Correlation between Standard Roughness Parameters Skewness and Kurtosis and Tribological Behaviour of Contact Surfaces. *Tribology International* **2012**, *48*, 102–112. https://doi.org/10.1016/j.triboint.2011.11.008.

758

759

760

4. Ba, E.C.T.; Dumont, M.R.; Martins, P.S.; Drumond, R.M.; Martins da Cruz, M.P.; Vieira, V.F. Investigation of the Effects of Skewness Rsk and Kurtosis Rku on Tribological Behavior in a Pin-on-Disc Test of Surfaces Machined by Conventional Milling and Turning Processes. *Materials Research 2021*. https://doi.org/10.1590/1980-5373-mr-2020-0435.

761

762

763

5. He, P.; Lu, S.; Wang, Y.; Li, R.; Li, F. Analysis of the Best Roughness Surface Based on the Bearing Area Curve Theory. *Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology* **2021**, p. 13506501211018937. https://doi.org/10.177/13506501211018937.

764

765

766

6. Pawlus, P.; Reizer, R.; Wieczorowski, M. Functional Importance of Surface Texture Parameters. *Materials* **2021**, *14*, 5326. https://doi.org/10.3390/ma14185326.

767

768

7. Yang, D.; Tang, J.; Xia, F.; Zhou, W. Rough Surface Characterization Parameter Set and Redundant Parameter Set for Surface Modeling and Performance Research. *Materials* **2022**, *15*, 5971. https://doi.org/10.3390/ma15175971.

769

770

8. Chen, J.; Tang, J.; Shao, W.; Li, X.; Yang, D.; Zhao, B.; Dong, H. A New Numerical Simulation Method of 3D Rough Surface Topography with Coupling 3D Roughness Parameters Sdr, Sdq, Spd, Spc, and Characteristic Functions. *Tribology International* **2024**, *200*, 110117. https://doi.org/10.1016/j.triboint.2024.110117.

771

772

773

9. Chen, J.; Zang, F.; Zhao, X.; Li, H.; Tong, Z.; Yuan, K.; Zhu, L. Generating Non-Gaussian Rough Surfaces Using Analytical Functions and Spectral Representation Method with an Iterative Algorithm. *Applied Mathematical Modelling* **2025**, 137, 115665. https://doi.org/10.1016/j.apm.2024.115665.

774

775

776

10. Francisco, A.; Brunetière, N. A Hybrid Method for Fast and Efficient Rough Surface Generation. *Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology* **2016**, *230*, 747–768. https://doi.org/10.1177/1350650115612116.

777

778

11. Borodich, F.M.; Pepelyshev, A.; Jin, X. A Multiscale Statistical Analysis of Rough Surfaces and Applications to Tribology. *Mathematics* **2024**, *12*, 1804. https://doi.org/10.3390/math12121804.

779

780

12. Francisco, A.; Brunetière, N.; Merceron, G. Gathering and Analyzing Surface Parameters for Diet Identification Purposes. *Technologies* **2018**, *6*, 75. https://doi.org/10.3390/technologies6030075.

781

782

13. Pawlus, P.; Reizer, R.; Wieczorowski, M. A Review of Methods of Random Surface Topography Modeling. *Tribology International* **2020**, *152*, 106530. https://doi.org/10.1016/j.triboint.2020.106530.

783

784

14. Wang, Y.; Azam, A.; Wilson, M.C.T.; Neville, A.; Morina, A. A Comparative Study for Selecting and Using Simulation Methods of Gaussian Random Surfaces. *Tribology International* **2022**, *166*, 107347. https://doi.org/10.1016/j.triboint.2021.107347.

785

786

15. Hu, Y.Z.; Tonder, K. Simulation of 3-D Random Rough Surface by 2-D Digital Filter and Fourier Analysis. *International Journal of Machine Tools and Manufacture* **1992**, *32*, 83–90. https://doi.org/10.1016/0890-6955(92)90064-N.

787

788

16. Bakolas, V. Numerical Generation of Arbitrarily Oriented Non-Gaussian Three-Dimensional Rough Surfaces. *Wear* **2003**, *254*, 546–554. https://doi.org/10.1016/S0043-1648(03)00133-9.

789

790

17. Johnson, N.L. Systems of Frequency Curves Generated by Methods of Translation. *Biometrika* **1949**, *36*, 149–176, [2332539]. https://doi.org/10.2307/2332539.

791

792

18. Collette, Y.; Siarry, P. *Optimisation multiobjectif*; Eyrolles: Paris, 2002.

793

794

19. Coello, C.C.; Lamont, G.B.; van Veldhuizen, D.A. *Evolutionary Algorithms for Solving Multi-Objective Problems*, 2 ed.; Genetic and Evolutionary Computation, Springer US, 2007. https://doi.org/10.1007/978-0-387-36797-2.

795

20. Deb, K.; Pratap, A.; Agarwal, S.; Meyarivan, T. A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II. *IEEE Transactions on Evolutionary Computation* **2002**, *6*, 182–197. https://doi.org/10.1109/4235.996017.

796

797

21. Blank, J.; Deb, K. Pymoo: Multi-Objective Optimization in Python. *IEEE Access* **2020**, *8*, 89497–89509. https://doi.org/10.1109/ACCESS.2020.2990567.

798

799

22. Fmodpy: A Lightweight, Efficient, Highly Automated, Fortran Wrapper for Python.

800

23. Multiprocessing — Process-based Parallelism. https://docs.python.org/3/library/multiprocessing.html.

801

24. Goodman, L.A. On the Exact Variance of Products. *Journal of the American Statistical Association* **1960**, *55*, 708–713. https://doi.org/10.1080/01621459.1960.10483369.

802

803

25. Maxima. Maxima, 2025.

804