Statistical Study of the Bias and Precision for Six Estimation Methods for the Fractal Dimension of Randomly Rough Surfaces

: The simulation and characterisation of randomly rough surfaces is an important topic in surface science, tribology, geo-and planetary sciences, image analysis and optics. Extensions to general random processes with two continuous variables are straightforward. Several surface generation algorithms are available, and preference for one or another method often depends on the specific scientific field. The same holds for the methods to estimate the fractal dimension D . This work analyses six algorithms for the determination of D as a function of the size of the domain, variance, and the input value for D , using surfaces generated by Fourier filtering techniques and the random midpoint displacement algorithm. Several of the methods to determine fractal dimension are needlessly complex and severely biased, whereas simple and computationally efficient methods produce better results. A fine-tuned analysis of the power spectral density is very precise and shows how the different surface generation algorithms deviate from ideal fractal behaviour. For large datasets defined on equidistant two-dimensional grids, it is clearly the most sensitive and precise method to determine fractal dimension.


Introduction
The description of surface topography by the superposition of harmonic waves with randomly oriented wave vectors and phase angles was pioneered by Longuet-Higgins in relation to the moving sea surface [1,2].It uses the relationship between the magnitude of the wave vector and the power spectrum (PS).Whitehouse et al. [3] were among the first to analyse the roughness of engineering surfaces under the concept of random process theory.With increasing measurement size and resolution, it was found that roughness increases with the length (area) of the measured domain [4][5][6][7].
Fractals are often associated with graphics in a plane, such as the Julia and Mandelbrot sets [8] or with problems like the coastline paradox [9], explained in terms of a Hausdorff dimension in a classical paper by Mandelbrot [10].A geometrical analysis of fractal dimension follows naturally from the graphical nature of the object studied.For random fractal surfaces, or more generally, for the self-affine stochastic processes found in many scientific disciplines, the problem must be approached differently.A mathematical definition of the basic problem is the most direct way to start.
The increase in the variance with an increase in the size of the domain over which the random variable is measured follows from Plancherel's theorem [11] for self-affine surfaces with a power-law PS: let the surface topography be described as a function z(x, y), with an isotropic PS: Here, r = √ k 2 + l 2 , the magnitude of the wave vector k with components (k, l).r l and r u are the lower and upper cut-off radii in the frequency domain.C is a constant.The measurement length corresponds to the maximum wavelength λ u = 2π/r l and the resolution of the measurement λ l = 2π/r u , where H is the Hurst exponent and the fractal dimension D = 3 − H.The variance of the height is given by the following: The operator < • > refers to the expected value.For r l /r u ≪ 1, the second term in the numerator can be neglected, and the standard deviation σ = √ z 2 ∼ r −h l ∼ λ h u , i.e., the root mean square roughness increases with the measurement length λ u .
It has been shown that λ l can be of atomic size in diamond coatings [12,13], while λ u can be as large as the size of the observable universe [14].The cut-off radii are generally not imposed by the physics of the phenomenon but by instrumental limitations, although a combination of measurement techniques at overlapping observation scales has been used to successfully extend λ u /λ l [4,12,13].In geo-and planetary sciences, multifractal descriptions are often used, where the value of H is allowed to vary with λ u and/or λ l [15][16][17].In materials science and tribology, an upper cut-off λ u is universally observed, beyond which the roughness no longer increases with the measurement length [12,13,18].λ u can be associated with the intuitive concept of coarseness, while H can be connected to complexity [19].In numerical simulations, λ u /λ l is determined by the size of the dataset, and λ u can be set to 1 without loss of generality.
Different definitions of the fractal dimension can give different results for D [20].A distinction must be made between the mathematical definition of D and the numerical methods used to estimate D from a discrete set of measured or simulated data.From a mathematical viewpoint, it is sufficient to find a single example in which two definitions of D produce different results.Although there are many counterexamples, the Minkowski and Hausdorff dimensions coincide in almost all cases, and objective criteria exist to test this coincidence [21].
In numerical methods, if D M is the estimate of D by a given method M. If two methods produce distinct values D M , it must be checked if they correspond to established mathematical exceptions.If not, the difference is more likely due to the numerical procedure missing some essential details of the geometry or because of the finite size and resolution of the available dataset.All definitions of fractal dimensions involve a limit for λ l → 0 .Limitations on λ u /λ l may severely distort the relationship between D M and D.
Most methods for determining D are only applicable to a limited subset of geometries [20].This study assumes that a smooth reference surface exists on a simply connected domain on which a function z(x, y) can be defined.For engineering surfaces and many geological and geomorphological studies, z is the height.More generally, z can be any random variable measured along a line or surface, such as the grey scale or colour value of pixels in an image [22][23][24][25] or the temperature of the cosmic background radiation [26,27].Contours of constant z 0 = z(x, y) (or w 0 = w(x, y, z) in three dimensions [28]) cannot be analysed by the methods studied here, except for box counting methods (BCM) and a strongly modified power spectral density method [29].This work will provide clear definitions of the bias and precision of the different methods M for the determination of H, with the goal of determining which method provides the best estimate for H. Section 2 of this paper describes the methods used to simulate the random surfaces and calculate H. Section 2.1 explains how this combination defines a random process.Section 2.2 presents three surface generation algorithms.It justifies the selection of only two of them for further analysis.The following seven subsections describe the methods used for determining H. Section 3 describes the simulation scheme, identifies the parameters of interest, and provides a statistical analysis of the results based on linear regression and analysis of variance (ANOVA).The results are discussed in terms of bias, dispersion, precision, and computational and information efficiency.This will allow us to reach a clear conclusion with respect to the problem posed in the title of this paper.

Surface Simulation and Characterisation as a Random Process
Simulated surfaces will be used to analyse the methods for estimating H.The process is illustrated in Figure 1.For Gaussian surfaces with zero mean, the input values are the variance < z 2 > in , H in , λ u , and λ l .These are deterministic values.The surface generation algorithm S acts on a finite set {g i } of random values, producing a single representation z S (x, y).Because of the random nature of {g i }, the mean value of µ S ̸ = 0 and < z 2 > S ̸ =< z 2 > in in general.Shifting and rescaling the simulated surface can correct these random variations.A numerical procedure M then calculates the simulated Hurst exponent H SM (or H out ), which is a random variable, which will be characterised by its mean and standard deviation.The following seven subsections describe the methods used for determining H. Section 3 describes the simulation scheme, identifies the parameters of interest, and provides a statistical analysis of the results based on linear regression and analysis of variance (ANOVA).The results are discussed in terms of bias, dispersion, precision, and computational and information efficiency.This will allow us to reach a clear conclusion with respect to the problem posed in the title of this paper.

Surface Simulation and Characterisation as a Random Process
Simulated surfaces will be used to analyse the methods for estimating H.The process is illustrated in Figure 1.For Gaussian surfaces with zero mean, the input values are the variance 〈 〉 ,  ,  , and  .These are deterministic values.The surface generation algorithm  acts on a finite set  of random values, producing a single representation   (, ).Because of the random nature of  , the mean value of   ≠ 0 and 〈 〉  ≠ 〈 〉 in general.Shifting and rescaling the simulated surface can correct these random variations.A numerical procedure ℳ then calculates the simulated Hurst exponent  ℳ (or  ), which is a random variable, which will be characterised by its mean and standard deviation.The bias is defined as  ℳ ( ) = 〈 ℳ 〉 −  .The variance 〈 ℳ 〉 ≠ 〈  〉 ; the right-hand side is the variance of the unknown Hurst exponent of the surface simulated with .The total dispersion can be characterised by the standard deviation  ( ) = 〈( ℳ −  ) 〉. 〈  〉 cannot be estimated independently of ℳ, but the smallest 〈 ℳ 〉 obtained from a given set of methods will give an upper bound estimate for 〈  〉.

Surface Generation
Surfaces used in this study are defined on an equidistant set of points  ,  , with , : 1 → 2 and  = 8, 9, or 10.The side of the square domain can be arbitrarily set equal to 1.Then,  = √2 and  = 2 .Three surface generation methods  are commonly used to generate (, ).The generalised Weierstrass-Mandelbrot function [30,31] ( = ) has recently recovered some popularity [32][33][34][35][36][37][38].It is described by the following: G is called the fractal roughness, and  is the lacunarity.The random phase angles  are chosen from a uniform distribution on −π, π .The observation that the second cosine term in the sum is independent from the random phase angles allows speeding up the calculation of the WM surface during Monte Carlo simulations.Even so, the method is not computationally efficient compared to the following algorithms.The bias is defined as S >; the right-hand side is the variance of the unknown Hurst exponent of the surface simulated with S. The total dispersion can be characterised by the standard deviation S > cannot be estimated independently of M, but the smallest < H 2 SM > obtained from a given set of methods will give an upper bound estimate for < H 2 S >.

Surface Generation
Surfaces used in this study are defined on an equidistant set of points x i , y j , with i, j : 1 → 2 n 0 and n 0 = 8, 9, or 10.The side of the square domain can be arbitrarily set equal to 1.Then, λ u = √ 2 and λ l = 2 −n 0 .Three surface generation methods S are commonly used to generate z(x, y).The generalised Weierstrass-Mandelbrot function [30,31] (S = W M) has recently recovered some popularity [32][33][34][35][36][37][38].It is described by the following: (3) G is called the fractal roughness, and γ is the lacunarity.The random phase angles φ mn are chosen from a uniform distribution on [−π, π].The observation that the second cosine term in the sum is independent from the random phase angles allows speeding up the calculation of the WM surface during Monte Carlo simulations.Even so, the method is not computationally efficient compared to the following algorithms.
The midpoint displacement method (S = MP) has been widely used in computer graphics [39] and in early studies on fractal landscapes [40,41].The original midpoint algorithm suffers from the problem that the initial few values can dominate the entire simulated surface, creating clear cross-patterns in the result.This can be corrected by modifying the way in which random values are added during each simulation step [39].It is also feasible to ensure periodicity of the surface, which is important in some simulations, e.g., in contact mechanics, where a Fourier transform boundary element method is employed [42,43].
The Fourier transform method (S = FT) [43,44] uses normally distributed variables g ij on a regular grid with a size of 2 n 0 × 2 n 0 .This set is transformed to the frequency domain using a fast Fourier transform (FFT), defining ĝij = F g ij .By multiplying with the PS P(k, l) in the frequency domain and applying the inverse Fourier transform (IFT), ĝij is filtered: where i = √ −1; the linear combination of real and imaginary parts in Equation ( 4) generates a purely real z(x, y).
It shall be remembered that the FFT is only an efficient way of summing cosine and sine terms in a discrete Fourier transform (DFT).This is equivalent to the original approach by Longuet-Higgins [1], but Equation (4) allows for the efficient inclusion of 2 2n 0 harmonic terms, which would require unacceptable computer times if the surface were generated by naive algorithms.However, by using the FFT, aliasing is introduced, which results in a statistical height distribution that is not Gaussian.This aspect can be improved by eliminating the longest wavelengths from the analysis [45].
Figure 2 shows three topographies generated by WM, MP, and FT.Their differences are seen in P S and the autocorrelation R S (Section 2.7).P W M shows many individual peaks corresponding to the individual combinations (m, n) in Equation (3) and has ridges along the coordinate axes.The rotational symmetry inherent to Equation (3) is visible in R W M as small wrinkles in the conical part of the surface.Elaborating Equation (3) is computationally inefficient and limits the use of WM for large-scale simulations. and  show pronounced edges along the coordinate axes.The preferential 90°-45° orientations imposed by the algorithm are clearly visible in  . consists of a noisy signal around the centre of the spectrum that decays according to Equation (1). consists of a conical peak superposed on a background ridge.Only the conical part is important for the analysis (Section 2.8).P MP and R MP show pronounced edges along the coordinate axes.The preferential 90 • -45 • orientations imposed by the algorithm are clearly visible in R MP .P FT consists of a noisy signal around the centre of the spectrum that decays according to Equation (1).R FT consists of a conical peak superposed on a background ridge.Only the conical part is important for the analysis (Section 2.8).

Box Counting Methods
Box counting methods (BCM) are, in principle, numerical implementations of the theoretical definition of a Minkowski dimension [18].For fractal curves defined on a unit square domain, the method consists of dividing the image into a set of m × m identical squares with length a = 1/m and counting the number n m of squares which contain at least one point of the curve.The box counting dimension is determined as the following: If there is an algorithm for infinitely refining the scale of the curve, m is limited only by computational precision, and the method is very accurate [46][47][48].For experimental data, m is limited by the spatial resolution.The procedure is analogous in 3 dimensions [28].
For functions defined as z(x, y) on a finite grid, the method produces considerable bias [49].Several recent publications have focussed on the rule that n m should be the infimum of all possible covering configurations at a given size [22,50].This infimum can be approximated by shifting the covering boxes with respect to their ideal position [51].Differential box counting (DBC) adapts for differences in vertical and horizontal scale by calculating the range ∆z over the domain and using boxes of vertical size ∆z/m.A lack of vertical resolution is addressed by considering only the highest and lowest boxes that contain part of the surface in each column of boxes; all intermediate boxes will contain a part of the surface, even if the corresponding points do not appear in the dataset [49].This modifies Equation (5) as the following: where the sum is taken over all m 2 boxes, ⌈•⌉ represents the ceiling function, and M i and m i are the maximum and minimum values within square i.The method is closely related to the Hall-Wood estimator [52,53].Differences between the original definition of a box counting dimension, BCM and DBC, are illustrated in Figure 3.
Fractal Fract.2024, 8, x FOR PEER REVIEW 6 of 22 vertical resolution is addressed by considering only the highest and lowest boxes that contain part of the surface in each column of boxes; all intermediate boxes will contain a part of the surface, even if the corresponding points do not appear in the dataset [49].This modifies Equation (5) as the following: where the sum is taken over all  boxes, ⌈•⌉ represents the ceiling function, and  and  are the maximum and minimum values within square i.The method is closely related to the Hall-Wood estimator [52,53].Differences between the original definition of a box counting dimension, BCM and DBC, are illustrated in Figure 3. Notice that Figures 3-5 use time series to exemplify the basic concepts of the methods and are for illustrative purposes only.The curve shown in Figures 3a and 4a corresponds to an arbitrary contour line of a surface produced by FT.The data series represents an arbitrary vertical section through this surface in arbitrary units.All methods described in this paper refer to the fully three-dimensional functions, as shown in Figure 2.  3a and 4a corresponds to an arbitrary contour line of a surface produced by FT.The data series represents an arbitrary vertical section through this surface in arbitrary units.All methods described in this paper refer to the fully three-dimensional functions, as shown in Figure 2.
Notice that Figures 3-5 use time series to exemplify the basic concepts of the methods and are for illustrative purposes only.The curve shown in Figures 3a and 4a corresponds to an arbitrary contour line of a surface produced by FT.The data series represents an arbitrary vertical section through this surface in arbitrary units.All methods described in this paper refer to the fully three-dimensional functions, as shown in Figure 2.  [39] for the estimation of the fractal dimension of a time series.The length of the lines between the dots is approximated by the length of the blue lines over intervals of width a in successive refinements.By Pythagoras' theorem, the method is sensitive to vertical scaling, and the finite resolution of the data means that the limit for  → 0 can only be approached crudely.[39] for the estimation of the fractal dimension of a time series.The length of the lines between the dots is approximated by the length of the blue lines over intervals of width a in successive refinements.By Pythagoras' theorem, the method is sensitive to vertical scaling, and the finite resolution of the data means that the limit for a → 0 can only be approached crudely.
vertical shift and counts the total range in individual intervals.Box ** is still missed.
Notice that Figures 3-5 use time series to exemplify the basic concepts of the methods and are for illustrative purposes only.The curve shown in Figures 3a and 4a corresponds to an arbitrary contour line of a surface produced by FT.The data series represents an arbitrary vertical section through this surface in arbitrary units.All methods described in this paper refer to the fully three-dimensional functions, as shown in Figure 2.  [39] for the estimation of the fractal dimension of a time series.The length of the lines between the dots is approximated by the length of the blue lines over intervals of width a in successive refinements.By Pythagoras' theorem, the method is sensitive to vertical scaling, and the finite resolution of the data means that the limit for  → 0 can only be approached crudely.

Triangular Prism Method
The triangular prism method (TPM) [40,41,54] is an extension of the Hausdorff dimension as used by Mandelbrot [10] to explain the coastline paradox [9].The original method is applied to contour lines in the plane, where the number n i of segments (yardsticks) with a given length l required to cover the entire contour is determined by geometrical construction.TPM extends this approach by covering a surface with a mesh of triangles of a given size.The total surface area s Tot is the sum of all areas obtained by approximating the surface by triangles over a grid of square domains with side a and the following equation: Because surface area does not linearly scale with height, this method is sensitive to the vertical scale, which would limit its use to self-similar surfaces [55].However, it has been proven that Equation (7) always converges to the correct value of D, but very small values of a may be required, corresponding to very large datasets [46].The difference between the original definition of a Hausdorff dimension and its application to a time series is shown in Figure 4.For more details on the geometry in the case of two-dimensional surfaces, the reader is referred to Ref. [49].

Detrended Fluctuation
Detrended fluctuation (DTF) was originally developed for the study of single-variable processes that show nonstationarity [56].Instead of analysing the variance of the function over a time interval ∆t, a linear or polynomial fit is made to the accumulated function over ∆t, and the residual variance of this fit is analysed as a function of ∆t.The method was adapted to two-variable random processes by Gu and Zhou [57,58], with an emphasis on multifractal analysis.
Here, the cumulative function Z(x, y) is calculated over each square subdivision with side a.A fitting polynomial f (x, y) is adjusted to Z(x, y) on the individual squares.The residual variance v a , averaged over all squares, obeys v a ∼ a H , hence the following equation: Linear fits and second-degree polynomials are commonly used for f (x, y) [49].This work uses the following Equation:

Roughness-Length Method
The roughness-length method (RLM) is based on the log-log plot of the variance measured over square areas of length a, according to Equation (2).Early Monte Carlo simulations [59] showed that its reliability is comparable to the power spectral method for time series.Its use is popular in the field of rock mechanics [60].It was used to analyse the quality of surfaces generated by the FT and WM methods [32].The number of simulations in this study was limited.
The implementation and use of the RLM are simple.Its theoretical basis (Equation ( 2)) is straightforward.Therefore, extended studies on its theoretical background or numerical improvement have not been published, contrary to TPM, BCM, DBC, and DTF.Its use for the characterisation of random surfaces is nonetheless interesting, as will be shown in Section 3.An illustration of the difference between RLM and DTF is shown in Figure 5 for the simplified case of time series.

Power Spectral Density and Related Functions
In most cases of interest, the function z(x, y) can be described with respect to a fixed mean value, which can be set to 0 without loss of generality.This is certainly true for the numerically generated surfaces studied here.Three closely related functions are defined for discrete datasets on a regular grid x i , y j with i, j : 1 → 2 n 0 .The autocorrelation of a real-valued function z(x, y) is defined as the following: (10) or for discrete data: The convolution in Equation ( 10) can be easily performed as a multiplication in the frequency domain.Likewise, the discrete convolution in Equation ( 11) is most easily performed as a multiplication after FFT followed by IFT.The power spectral density (PSD) is the Fourier transform of the autocorrelation function, i.e., The variogram [61], height correlation function [62], or structure function [63] is defined as follows: Given that γ(ξ, ψ) = 1 − R(ξ, ψ), all three functions hold the same information [62,64].Their numerical implementation and fields of application may differ, which may induce differences in the estimated value of D M , with M referring to PSA (power spectral analysis) or SFA (structure-function analysis).
Variograms or structure functions are widely used in geosciences [15][16][17] but are seldom found in tribology and surface sciences.One reason is that the variogram can be applied without defining a reference level (mean value) and can be calculated for limited datasets obtained from irregularly spaced sample points, which are typical for geological exploration [65].The other reason is that surface profilers and atomic force microscopes usually come with an extensive postprocessing package, which includes background corrections and PSD calculation in the form of black-box procedures [66,67].

Tuning of the PSA
A common strategy to calculate the PSD from a measured or simulated surface is to apply an FFT to each individual line z x, y j and averaging the results.In this work,  An important aspect of PSA is illustrated in Figure 7.  () correctly reproduces the power-law PSD used as an input for FT.This is not a trivial result because the algorithm used to determine H by PSA does not "know" what kind of surface is used as an input, i.e., PSA does not use a trivial deconvolution of  and () because the input  are unknown to the algorithm. () flattens at high frequencies.Although the theoretical analysis of MP shows that this method produces a power-law PSD [39], its numerical implementation induces a deviation of this theoretical behaviour. () shows the strongest deviation from the ideal behaviour.The peaks in the spectrum are due to the finite values of N and M in Equation (3).Each peak corresponds to a circle of spikes in  (, ) (Figure 2).Using PSA without accounting for this behaviour leads to large errors [32].2.9.Tuning of the SFA The data in Figure 6b present two inconveniences for correct regression analysis.The residual variance is not constant, and the density of data points increases strongly with r. Figure 6c shows the results of applying a low-pass filter to the average data of Figure 6b, with threshold radii r Thr = 1, 2 −4 , 2 −8 .With the highest threshold, linearity is respected over the entire log 2 r range, but dispersion is considerable for large values of r.As r Thr is increased, dispersion decreases, but the curve flattens for small r Thr .This phenomenon is artificial and should not be confused with an upper cutoff.For the statistical simulations reported in Section 3.3, r Thr = 2 −5 was used, with linear regression performed for r : 2 −6 → 2 −1 .
An important aspect of PSA is illustrated in Figure 7. P FT (r) correctly reproduces the power-law PSD used as an input for FT.This is not a trivial result because the algorithm used to determine H by PSA does not "know" what kind of surface is used as an input, i.e., PSA does not use a trivial deconvolution of ĝij and P(r) because the input ĝij are unknown to the algorithm.P MP (r) flattens at high frequencies.Although the theoretical analysis of MP shows that this method produces a power-law PSD [39], its numerical implementation induces a deviation of this theoretical behaviour.P W M (r) shows the strongest deviation from the ideal behaviour.The peaks in the spectrum are due to the finite values of N and M in Equation (3).Each peak corresponds to a circle of spikes in P W M (k, l) (Figure 2).Using PSA without accounting for this behaviour leads to large errors [32].
rithm used to determine H by PSA does not "know" what kind of surface is used as an input, i.e., PSA does not use a trivial deconvolution of  and () because the input  are unknown to the algorithm. () flattens at high frequencies.Although the theoretical analysis of MP shows that this method produces a power-law PSD [39], its numerical implementation induces a deviation of this theoretical behaviour. () shows the strongest deviation from the ideal behaviour.The peaks in the spectrum are due to the finite values of N and M in Equation ( 3).Each peak corresponds to a circle of spikes in  (, ) (Figure 2).Using PSA without accounting for this behaviour leads to large errors [32].

Tuning of the SFA
A reference curve for the structure function () can be calculated by eliminating C from Equations ( 1) and (2) and taking the limit for  → ∞ for ().This produces the following: With the help of mathematical software, the IFT of Equation ( 14) can be found as the following:

Tuning of the SFA
A reference curve for the structure function γ(r) can be calculated by eliminating C from Equations ( 1) and (2) and taking the limit for r u → ∞ for P(r).This produces the following: with the help of mathematical software, the IFT of Equation ( 14) can be found as the following: where Γ(•) is the gamma function, and 1 F 2 (•) is the generalised hypergeometric function.As R(ρ) has zero crossings (Figure 8a), its representation in a log 2 − log 2 plot is meaningless.The log-log plot for γ(ρ) = 1 − R(ρ) shows asymptotic linear behaviour for log 2 ρ < 0 (Figure 8b). Figure 8c shows the average of γ(ρ) over all points with the same ρ for H = 0.25, 0.5, 0.75 and S = FT.
As () has zero crossings (Figure 8a), its representation in a log − log plot is meaningless.The log-log plot for () = 1 − () shows asymptotic linear behaviour for log  < 0 (Figure 8b). Figure 8c shows the average of () over all points with the same  for H = 0.25, 0.5, 0.75 and  = .Contrary to the radial average of the PSD, no filtering is needed for SFA, but the range where the theoretical prediction and linear approach are valid is limited compared with the filtered PSD.There is a visible difference between the slope of the simulated () and the theoretical curves.The horizontal offset of the simulated () is a random value with broad dispersion.For calculating H, this offset is irrelevant.

Simulation Scheme
Five factors were considered in this study: the surface generating algorithm, the size of the dataset, the vertical scale  = 〈 〉 , the method for calculating H, and the value of Hin.The domain is defined as a grid of 2 × 2 points, with  taking values of 8, 9, and 10.An earlier study [49] showed that TPM, and to a lesser extent, DTF, depend on vertical scaling; three roughness values  equal to 0.25, 0.5, and 1 are considered.The methods used are DBC, TPM, DTF, PSA, SFA, and RLM.
The value of  was varied from 1/128 to 127/128 in steps of 1/64.These values are chosen to exploit the perfectly parallel nature of the calculations on a processor with eight Contrary to the radial average of the PSD, no filtering is needed for SFA, but the range where the theoretical prediction and linear approach are valid is limited compared with the filtered PSD.There is a visible difference between the slope of the simulated γ(ρ) and the theoretical curves.The horizontal offset of the simulated γ(ρ) is a random value with broad dispersion.For calculating H, this offset is irrelevant.

Simulation Scheme
Five factors were considered in this study: the surface generating algorithm, the size of the dataset, the vertical scale σ 0 = < z 2 > in , the method for calculating H, and the value of H in .The domain is defined as a grid of 2 n 0 × 2 n 0 points, with n 0 taking values of 8, 9, and 10.An earlier study [49] showed that TPM, and to a lesser extent, DTF, depend on vertical scaling; three roughness values σ 0 equal to 0.25, 0.5, and 1 are considered.The methods used are DBC, TPM, DTF, PSA, SFA, and RLM.
The value of H in was varied from 1/128 to 127/128 in steps of 1/64.These values are chosen to exploit the perfectly parallel nature of the calculations on a processor with eight cores.Two surfaces were created for each combination of S, n 0 and H in .For each of these surfaces, the three values of σ 0 and all six methods were applied.The scheme allows the assessment of bias and dispersion by standard regression analysis and the comparison of results by ANOVA.

Data Analysis
For both S, each M and n 0 , H SM (or H out ) is plotted vs. H in , for each value of σ 0 .The linear regression is calculated with the following Equation: is determined by standard least squares.A fourth regression line is estimated for the pooled data, disregarding the value of σ 0 .The residual variance σ 2 Out of the pooled data is compared with the mean residual variance of the three individual regression lines.The resulting F-ratio is compared with the 90% quantile of the F-distribution with the corresponding degrees of freedom, and the p-value is calculated.If p > 0.1, it is concluded that M is insensitive to σ 0 in the range evaluated.
The value of a 1 characterises the bias of the methods.Although this value cannot be entirely independently of a 0 , a 1 affects the precision of a method.Having determined H SM (Equation ( 16)) and knowing a 0 and a 1 , a corrected value can be found as the following: σ out represents the dispersion of the method.The precision of M can be defined as the difference between H in and H crr for a measured H SM , quantified, for example, as the standard deviation between the predicted value H crr and H in is as follows: σ out interacts with a 1 in determining σ Pred .Small values of a 1 magnify σ out when projected on the horizontal axis of the H SM vs. H in -plot.σ Pred is obtained from the regression analysis on a H in vs. H out -plot, under the assumption that the residuals of this regression are independent of H SM .If the residuals follow a normal distribution, the confidence intervals can be calculated.The assumption of normality was assessed using the Kolmogorov-Smirnov test.

Results for S = FT
The complete results for n 0 = 10, σ 0 = 0.25, 0.5, and 1 and all M for S = FT are presented in Figure 9.As an example, Table 1 presents the complete analysis of the regression results for DTF and DBC.For all M, the regression results are independent of σ 0 , although the P F -values for TPM are close to 0.1.A summary of the pooled data for all M is shown in Tables 2 and 3, for FT and MP and all n 0 , together with P F and the 90% confidence interval for H crr .None of the P KM indicate deviations from normality, but note that for n 0 = 8, P F = 0.025 for TPM, confirming its dependence on σ 0 , as observed in an earlier study [49].All methods except RLM and PSA present a visible increase in the variance with H in .This effect is not studied in detail and is neglected in the statistical analysis.
hough the  -values for TPM are close to 0.1.A summary of the pooled data for all ℳ is shown in Tables 2 and 3, for FT and MP and all  , together with  and the 90% confidence interval for  .None of the  indicate deviations from normality, but note that for  = 8,  = 0.025 for TPM, confirming its dependence on  , as observed in an earlier study [49].All methods except RLM and PSA present a visible increase in the variance with  .This effect is not studied in detail and is neglected in the statistical analysis.Table 1.Results of the regression analysis for n 0 = 10, σ 0 = 0.25, 0.5, and 1 for S = FT and M = DTF and DBC.The first three lines present the results for individual σ 0 .R 2 is the coefficient of determination, and n is the number of simulations.The next line presents the results for the pooled values.Line 6 presents the ANOVA results for the effect of σ 0 ; F 90% and P σ 0 are the 90% confidence level and p-value for the Fisher test, respectively.Line 8 gives the residual standard deviation σ Pred , the 90% confidence interval for H crr , and the p-value for the Kolmogorov-Smirnov test (P KM ) applied to the residuals of the H in vs. H SM regression line.The results for S = MP differ considerably from S = FT.Figure 10 shows that SFA has a large bias, especially for high values of H (low D), where the predominance of the first few values in MP is not masked by the strong fluctuations associated with low H.In the statistical results (Table 3), the precision for H TPM and H DBC varies significantly with σ 0 , as found in an earlier study [49].This effect depends on both the surface simulation method and the algorithm used for the calculation of D. For DTF and TPM, the results depend significantly on  ; for DBC and SFA, the bias is much higher than in  ℳ .

Bias, Dispersion, and Precision
Most studies on the fractal dimension of measured data conclude that different methods yield different results, as is clearly confirmed here.However, if the exact input value  is a priori unknown, there is no way to decide which  ℳ is closer to this value.In the field of image recognition, there is a tradition of applying different methods to a standard library of bitmaps and evaluating whether a method can distinguish between similar images [23][24][25]51].Methods are then compared based on their consistency with earlier research, but such an approach does not provide information on bias or statistical dispersion.
Using simulated data with known  , two factors become clear at once.The first is the bias ( ), as captured by  .DBC is clearly underperforming in this regard.However, by Equation ( 17), this can be corrected by if  and  are calibrated through numerical simulation.Dispersion (σout) cannot be corrected.In this aspect, DBC outperforms the other methods, except for PSA.σout combines with B to determine  and the confidence interval for  .A low value of  will magnify the dispersion on  for a given  .The low dispersion found in DBC cancels out the effect of bias for  = .This effect has not received due attention in the literature; most studies comparing different methods for the determination of D do not define precision.

Computational Efficiency
In the context of this study, efficiency can be interpreted either in terms of computational speed or the optimal use of information.Computation speed is often a secondary concern if the effort involved in measuring (, ) is much greater than the work involved in the calculation of D. However, as the size of the datasets increases, the For DTF and TPM, the results depend significantly on σ 0 ; for DBC and SFA, the bias is much higher than in H FTM .

Bias, Dispersion, and Precision
Most studies on the fractal dimension of measured data conclude that different methods yield different results, as is clearly confirmed here.However, if the exact input value H in is a priori unknown, there is no way to decide which H SM is closer to this value.In the field of image recognition, there is a tradition of applying different methods to a standard library of bitmaps and evaluating whether a method can distinguish between similar images [23][24][25]51].Methods are then compared based on their consistency with earlier research, but such an approach does not provide information on bias or statistical dispersion.
Using simulated data with known H in , two factors become clear at once.The first is the bias B(H in ), as captured by a 1 .DBC is clearly underperforming in this regard.However, by Equation ( 17), this can be corrected by if a 0 and a 1 are calibrated through numerical simulation.Dispersion (σ out ) cannot be corrected.In this aspect, DBC outperforms the other methods, except for PSA.σ out combines with B to determine σ Pred and the confidence interval for H crr .A low value of a 1 will magnify the dispersion on H crr for a given H out .The low dispersion found in DBC cancels out the effect of bias for S = PS.This effect has not received due attention in the literature; most studies comparing different methods for the determination of D do not define precision.

Computational Efficiency
In the context of this study, efficiency can be interpreted either in terms of computational speed or the optimal use of information.Computation speed is often a secondary concern if the effort involved in measuring z(x, y) is much greater than the work involved in the calculation of D. However, as the size of the datasets increases, the computational cost may blow up.In applications such as real-time pattern recognition, computational speed is essential.PSA, DBC, and RLM will have an advantage here.
To compare computational efficiency, it must be ensured that the most efficient algorithm is used for each method.A straightforward example is the use of FFT in SFA and PSA, but efficient algorithms are also known for BCM [68] and direct calculation of R(ρ) and γ(ρ) [69].As the algorithms in this study have not been fully optimised for efficiency, only the most evident effects will be mentioned.As an example, the full PSA routine used here for n 0 = 10 takes 5 ′′ using a single core on an 11th Gen Intel Core i9-11950H @ 2.60 GHz.Additional software optimisation can certainly increase this speed.
Calculation of the cumulative function in DTF and of the variance in RLM can be significantly sped up using the results at length scale a to obtain the values for 2a.A naive calculation of z 2 was used in RLM because, even without optimisation, this method is very fast.In contrast, DTF blows up rapidly with n 0 .The use of the cumulative function instead of the original values, combined with a regression analysis for each subset in each step, generates computational overhead.TPM is not very efficient either and is not very precise.
Further development of DTF and TPM to analyse physical measurements in the form of z x i , y j shows little promise.Considering the algorithms used in this study, DBC, RLM, PSA, and SFA are at least ten times more efficient than TPM, which is significantly faster than DTF.SFA takes about twice the time of PSA simply because of the calculation of the IFT in addition to the FFT.Given the lower precision of SFA, PSA is preferred.

Information Efficiency
Information efficiency describes how much of the available information in the dataset is effectively used to estimate D, including the individual values of z x i , y j and their spatial correlation.Power spectral analysis and autocorrelation were originally developed for the analysis of time series.Contact surface profilers measure roughness along a line [3]; hence, estimating fractal dimensions by analysing data along a line is a deep-rooted tradition [61,63,70,71].
Most modern applications generate large datasets over a rectangular array of regularly spaced points, as is the case for images obtained by digital cameras, atomic force microscopes, optical roughness measurements, or photogrammetry data.To analyse single lines out of a full set of f x i , y j data induces a considerable loss of information and may even distort the conclusions of the analysis [72,73].The effect is clearly illustrated in Figure 6a, where the power spectrum of single lines shows extreme dispersion, but the average value of the 1025 lines appears as a well-defined straight line.Single-line analysis should never be performed if data are available over a 2D domain.
Even so, DTF, TPM, DBC, and RLM do not use all available information because they are typically implemented by subsequent division of intervals by two in each step.For the periodic surfaces generated by FT, there are 2 2n 0 possible positions to collocate the corner of the first interval.While the computational cost of DTF and TPM does not encourage exploration of this effect, it has been used to decrease the statistical dispersion in DBC [20,21,51] and can be used in RLM.According to the definition of the Minkowski dimension [20], in DBC, the infimum of n m should be determined, while in RLM, the mean value can be used.Calculating all possible values would significantly increase the calculation time for RLM and DBC.SFA and PSA are not subject to this observation.

Effect of Vertical Scaling and Resolution
Vertical scaling clearly affects TPM but was found to be insignificant in all other methods for FT.It also affects DTF for MP.It is always possible to apply vertical rescaling such that σ 0 is of the same order of λ u if the significant number of digits for x, y, and z is equal, i.e., image depth should be equal to image resolution [49].If not, scaling only increases the gaps between consecutive vertical values.Many scientific cameras obey this condition, and optical roughness measurements often have a higher vertical resolution than lateral resolution.In numerical simulations, the vertical resolution is limited by machine precision, while the lateral resolution is defined by n 0 .
In terms of cameras and measurement systems available in many laboratory settings, n 0 = 10 is rather modest.However, for time-consuming measurements, large sample sets, or computationally complex simulations, it may still be interesting to work at lower n 0 .It follows from Tables 2 and 3 that σ out increases with decreasing n 0 .The effect is small for all methods except for PSA, where the increase in σ out is proportional to the reduction in spatial resolution.Still, at n 0 = 8, PSA outperforms all other methods at n 0 = 10, i.e., better estimates can be obtained using smaller datasets.

Effect of the Surface Generation Algorithm
PSA is more precise than the other methods for FT by a factor of 10.P FT (r) is related to the input function P(r) in Equation ( 4) by ĝij , which randomises the result.It has been claimed that using PSA with S = FT is "unfair" [32] because one is presumably the inverse of the other.It was pointed out in Section 2.8 that this is not the case; the implementation of PSA used in this work is independent of the surface generating algorithm.The excellent results obtained by PSA from surfaces generated by FT show that both methods are consistent.
The radial PS of surfaces defined by MP has a reduced linear section in the log-log plot (Figure 7).DTF is slightly more successful here than for S = FT (Tables 2 and 3) because the principal characteristic of DTF is to filter out effects of nonstationarity, which is inherent to MP [37].Zhang et al. [32] detected a breakdown of PSA for WM.This can be associated with the discrete peaks in the PS (Figures 2b and 7c).These peaks are the reason why WM was not included in the analysis, as they do not appear in the experimental spectra available in the literature [12][13][14][15][16][17][18].
As a reference, a series of simulations was performed for WM with n 0 = 10 using RLM and PSA.The latter method does not produce useful results over the entire range of H. RLM gave H W M RLM = 0.082 + 0.81H in , with σ out = 0.02.This is comparable to RLM and SFA for FT.It follows that an approximate value for H can be obtained, even if the initial assumption (Equation (1)) is not fulfilled.Without the analysis of the PS, the peaks in the PS for WM would be overlooked.An extensive statistical analysis of the available surface generation algorithms, in terms of their PS, H, and the probability density distribution of surface heights, is an interesting topic for future investigation.
If the purpose of a study is to simulate a surface which conforms to Equation (1), the peaks in the PS created by WM induce a deviation from the basic hypothesis.Whether such preferred frequencies are important in the simulation depends on the application.For optical phenomena, they can be very important as they may interact with the frequency of the optical signal.In contact and fracture mechanics, the effect of the details of the PS has not been critically analysed.MP introduces a smaller deviation than WM, but its effect on specific applications has not been quantified either.
If the goal is to estimate the fractal dimension of some measured dataset z = z(x, y), it is important to investigate whether the hypothesis of a power-law PS is fulfilled.This can be carried out by plotting P √ k 2 + l 2 as in Figure 7 or by studying the residual standard deviation and coefficient of determination of the linear fit to the log 2 − log 2 plot for P(r).Failure to do so may lead to incorrect or incomplete conclusions about the character of the data.
The simulations made in this work are meant as a reference for the analysis of experimental data.In an experimental setting, D is a priory unknown, so it is not possible to compare and calibrate the methods objectively.For measurements on an irregular grid, the present version of PSA is not feasible, and other methods may be selected.For experimental results on a regular grid, which may show deviations from the exact power law (Equation ( 1)), PSA can detect such deviations, while techniques such as DBC and DTF will mask them, which may represent an unacceptable loss of information.
Refined methods for the characterisation of more complex random processes on a two-dimensional domain can be devised.Substituting the simple power law by more advanced fitting functions may be the first approach to this problem.For the simulation of surfaces which purposedly deviate from Equation (1), the FT technique can use more refined expressions for P(k, l) in Equation (4).

Conclusions
A basic rule in scientific research is that measurement techniques must be calibrated to reliable standards.Neither calibration techniques nor standards have been clearly defined for determining the Hurst exponent of random fractal surfaces defined as z = z(x, y) on a unit square.This study defines the bias, dispersion, and precision of the methods.Bias refers to the difference between the output and input values of the Hurst exponent.It is a function of the deterministic input value H in .Dispersion is the standard deviation of the random variable H out , and precision is the standard deviation after correction for bias, i.e., on H crr .
The consistency and high precision of the power spectral analysis for surfaces produced by the Fourier transform method establishes the latter as a useful standard for the other methods.The detrended fluctuation and triangular prism methods can be rejected because of their computational cost and low precision.Based on the bias, differential box counting should be rejected.Box counting methods were developed to characterise contours z 0 = z(x, y) in the plane but not for surfaces z = z(x, y) in three dimensions, where they underperform.
Differential box counting and detrended fluctuation can filter out deviations of the perfect power law behaviour of the power spectrum, but power spectral analysis accurately detects these deviations.Future research should, therefore, focus on exploring the unique capabilities of power spectral analysis, rather than trying to improve methods that are computationally inefficient or make suboptimal use of the information available in the measurements.
Author Contributions: J.L.F.A., Methodology, Software, Investigation, Formal analysis, Visualisation, and Writing-Original draft; C.G.F., Resources, Funding acquisition, Supervision, and Writing-Review and editing; V.H.J., Conceptualisation, Methodology, and Writing-Review and editing; F.V.V., Software, Supervision, Writing-Review and editing; R.S., Conceptualisation, Methodology, Software, Formal analysis, Supervision, and Writing-Final text.All authors have read and agreed to the published version of the manuscript.

Funding:
The PhD scholarship for JLFA was provided by CONACYT (Mexico).This research was partially funded by the DGAPA project PAPIIT-IA106720.

Figure 1 .
Figure 1.Representation of the process which creates  ℳ for a given input value  through a surface generating algorithm  and a method ℳ.By the random input  ,   (, ) is a random process, and  ℳ is a random value.

Figure 1 .
Figure 1.Representation of the process which creates H SM for a given input value H in through a surface generating algorithm S and a method M. By the random input {g i }, z S (x, y) is a random process, and H SM is a random value.

Figure 3 .
Figure 3. (a) Original concept of BCM as a method to estimate the Minkowski dimension of a fractal contour.The grey boxes contain part of the curve and are counted in  .(b) BCM with vertical scaling.The boxes marked by * and ** are missed because the curve passes through the box, but there are no data points inside the box due to a lack of vertical resolution.(c) DBC allows for a vertical shift and counts the total range in individual intervals.Box ** is still missed.

Figure 3 .
Figure 3. (a) Original concept of BCM as a method to estimate the Minkowski dimension of a fractal contour.The grey boxes contain part of the curve and are counted in n m .(b) BCM with vertical scaling.The boxes marked by * and ** are missed because the curve passes through the box, but there are no data points inside the box due to a lack of vertical resolution.(c) DBC allows for a vertical shift and counts the total range in individual intervals.Box ** is still missed.Notice that Figures3-5use time series to exemplify the basic concepts of the methods and are for illustrative purposes only.The curve shown in Figures3a and 4acorresponds to an arbitrary contour line of a surface produced by FT.The data series represents an arbitrary vertical section through this surface in arbitrary units.All methods described in this paper refer to the fully three-dimensional functions, as shown in Figure2.

Figure 4 .
Figure 4. (a) Original concept of the yardstick method to estimate the Hausdorff dimension of a fractal contour.The contour is approximated by lines of a fixed length.(b) Trapezium method[39] for the estimation of the fractal dimension of a time series.The length of the lines between the dots is approximated by the length of the blue lines over intervals of width a in successive refinements.By Pythagoras' theorem, the method is sensitive to vertical scaling, and the finite resolution of the data means that the limit for  → 0 can only be approached crudely.

Figure 4 .
Figure 4. (a) Original concept of the yardstick method to estimate the Hausdorff dimension of a fractal contour.The contour is approximated by lines of a fixed length.(b) Trapezium method[39] for the estimation of the fractal dimension of a time series.The length of the lines between the dots is approximated by the length of the blue lines over intervals of width a in successive refinements.By Pythagoras' theorem, the method is sensitive to vertical scaling, and the finite resolution of the data means that the limit for a → 0 can only be approached crudely.

Figure 4 .
Figure 4. (a) Original concept of the yardstick method to estimate the Hausdorff dimension of a fractal contour.The contour is approximated by lines of a fixed length.(b) Trapezium method[39] for the estimation of the fractal dimension of a time series.The length of the lines between the dots is approximated by the length of the blue lines over intervals of width a in successive refinements.By Pythagoras' theorem, the method is sensitive to vertical scaling, and the finite resolution of the data means that the limit for  → 0 can only be approached crudely.

Figure 5 .
Figure 5. (a) Illustration of RLM for a time series.The vertical variation (standard deviation) is marked as dashed lines relative to the mean (blue line).(b) Cumulative data (black) with superposed piecewise fitting curves (blue, Equation (9)), according to DBC.(c) Fitting residuals, showing the residual standard deviation (dashed line) for each interval.Standard deviations are shown for plotting purposes; RLM and DBC use the corresponding variance.
calculated and averaged over all values with identical radii.Both approaches are illustrated in Figure6.While there is a large dispersion on the individual P FT k, l j , the average over all j is restricted to a very narrow band.P FT √ k 2 + l 2 shows a much larger spread.Meanwhile, for j : 1 → 1025 for the line profiles, there are 82,798 different values for r = √ k 2 + l 2 in a 1025 × 1025 grid (1,050,625 points).The lower number of values for each value of r produces a larger variance of the sample mean.Fractal Fract.2024, 8, x FOR PEER REVIEW 9 of 22over the entire log  range, but dispersion is considerable for large values of r.As  is increased, dispersion decreases, but the curve flattens for small  .This phenomenon is artificial and should not be confused with an upper cutoff.For the statistical simulations reported in Section 3.3,  = 2 was used, with linear regression performed for : 2 → 2 .

Figure 6 .
Figure 6.(a) log − log plot of  ,  , : 1 → 1025 for a surface with 1025 × 1025 points, together with the average for all 1025-line profiles.(b) Individual values and averages over all points with the same radius.(c) Filtered data.As the low-pass threshold radius  decreases, dispersion is reduced, as is the linear part of the plot.

Figure 6 .
Figure 6.(a) log 2 − log 2 plot of P FT k, l j , j : 1 → 1025 for a surface with 1025 × 1025 points, together with the average for all 1025-line profiles.(b) Individual values and averages over all points with the same radius.(c) Filtered data.As the low-pass threshold radius r Thr decreases, dispersion is reduced, as is the linear part of the plot.

Figure 7 .
Figure 7. Filtered mean values of the log 2 − log 2 plot for (a) P FT (r); (b) P MP (r); and (c) P W M (r).

Figure 9 .
Figure 9. Plots of  ℳ vs.  for three values of  and all ℳ, with the ideal case  ℳ =  marked as a red line.The methods are listed in order the precision  .Figure 9. Plots of H FTM vs. H in for three values of σ 0 and all M, with the ideal case H FTM = H in marked as a red line.The methods are listed in order the precision σ Pred .

Figure 9 .
Figure 9. Plots of  ℳ vs.  for three values of  and all ℳ, with the ideal case  ℳ =  marked as a red line.The methods are listed in order the precision  .Figure 9. Plots of H FTM vs. H in for three values of σ 0 and all M, with the ideal case H FTM = H in marked as a red line.The methods are listed in order the precision σ Pred .

Figure 10 .
Figure 10.Plots of  ℳ vs.  for three values of  and all ℳ.For DTF and TPM, the results depend significantly on  ; for DBC and SFA, the bias is much higher than in  ℳ .

Figure 10 .
Figure 10.Plots of H MPM vs. H in for three values of σ 0 and all M.For DTF and TPM, the results depend significantly on σ 0 ; for DBC and SFA, the bias is much higher than in H FTM .
Section 2 of this paper describes the methods used to simulate the random surfaces and calculate H. Section 2.1 explains how this combination defines a random process.Section 2.2 presents three surface generation algorithms.It justifies the selection of only two of them for further analysis.
provides the best estimate for H.

Table 2 .
Summary of the statistical results for S = FT and all M. The results are pooled over all σ 0 .Notice how σ out and the 90% CI on H crr increase with decreasing n 0 .

Table 3 .
Summary of the statistical results for S = MP and all M. The results are pooled over all σ 0 .Generally, bias is higher, and precision is lower than in FT.
nm Random phase angles in WM M Method to determine H (generic) Range of z S x i , y j over entire domain M i Maximum of z S x i , y j in sub-domain i. m i Minimum of z S x i , y j in sub-domain i. f (x, y) Regression function for DTF v a Residual variance for sub-domain size a c 00 , c 10 , c 01 , c 20 , c 02 , c 11 , a 0 , a 1 parameters s Tot Total surface determined in TPM P S (k, l), P S (r) Power spectrum for method S R S (x, y), R S (ρ) Autocorrelation for method S γ S (ρ) Structure function for method S r Thr Treshold radius for low-pass filter H SM Estimated H for S, M, and H in H crr Corrected H for S, M, and H in B SM Bias for methods S and M σ out Standard error on H SM (dispersion)