Spherical Cap Harmonic Analysis (SCHA) for Characterising the Morphology of Rough Surface Patches

We use spherical cap harmonic (SCH) basis functions to analyse and reconstruct the morphology of scanned genus-0 rough surface patches with open edges. We first develop a novel one-to-one conformal mapping algorithm with minimal area distortion for parameterising a surface onto a polar spherical cap with a prescribed half angle. We then show that as a generalisation of the hemispherical harmonic analysis, the SCH analysis provides the most added value for small half angles, i.e., for nominally flat surfaces where the distortion introduced by the parameterisation algorithm is smaller when the surface is projected onto a spherical cap with a small half angle than onto a hemisphere. From the power spectral analysis of the expanded SCH coefficients, we estimate a direction-independent Hurst exponent. We also estimate the wavelengths associated with the orders of the SCH basis functions from the dimensions of the first degree ellipsoidal cap. By windowing the spectral domain, we limit the bandwidth of wavelengths included in a reconstructed surface geometry. This bandlimiting can be used for modifying surfaces, such as for generating finite or discrete element meshes for contact problems. The codes and data developed in this paper are made available under the GNU LGPLv2.1 license.


Introduction
The mechanics of granular materials, such as soil mechanics, investigates the interaction of particles at various loading conditions where the shape, size and material properties of the particles can vary widely.Numerical models are therefore often used to study how the mechanical behaviour of a material is influenced by the various particle morphology and microstructures created by these particles.Similar problems exist in structural engineering when modelling, for example, the microstructure of concrete aggregates or the various stones in stone masonry elements, as the shape and roughness of the particles affect the matrices of these materials.
Addressing these problems numerically requires the modelling of randomly shaped closed objects and their interactions.In such models, the shape of the particle is typically represented by a finite or discrete element mesh, while the roughness of the surface is modelled through laws that describe the interaction of two surfaces, such as friction models.Efficiently separating "shape" and "roughness" of a randomly shaped particle is a challenging problem, though work in this area helps us better understand the associated constitutive laws in numerical modelling.In recent decades, the concept of spherical harmonics (SH) analysis has been frequently used to describe the shape and roughness of stones or aggregates.However, while this approach is transformative for closed objects, it is difficult to apply to the morphology of nominally flat and open surfaces (topological disks with free edges).In this paper, we therefore propose the use of spherical cap harmonic analysis (SCHA) for describing the geometry of open shapes.The difference here is that while a traditional spherical harmonics analysis (SHA) describes a function mapped onto an entire unit sphere, the SCHA describes a function mapped onto a spherical cap via a half-angle θ c (see Figure 1).SCHs are widely used in geophysics to describe field data on planetary caps, but to our knowledge have never been used to describe the geometry of a patch with free edges.Here we give example problems where SCHA can be useful for describing the topological characterisation of a patch (i.e. an open shape with a free edge): • For studying the topology of a rough surface patch, such as a rock surface, and quantifying its fractal dimension (FD).
• For characterising different roughness regions on rocks or faults, such as sedimentary rocks made of several strata.The traditional SH approach, which assumes a closed object, only exploits basis functions to fit all the surface variations and does not distinguish regions of different roughness (multifractal surfaces).
• For matching objects based on surface characteristics.This can be useful for medical and forensic anthropology applications, such as matching fractured bones or objects (see [1] and [2] as examples), or for mechanical applications, such as characterising the surfaces of fragmented rocks or soil particles [3,4] by directly comparing and matching the scanned fractured regions.
In this paper, we propose a new method for characterising the morphology of open surface patches using SCH, focusing specifically on nominally flat patches, such as when the global radius of the curvature of the targeted patches approaches infinity.The herein proposed SCHA method allows us to study the structure of any arbitrary surface patch with faster convergence than the traditional SHA as it scales to any level of detail with reasonable computational complexity.It can also be used for reconstructing surfaces and modifying the morphology of digital twins of real surfaces, which is useful when numerically studying contact problems such as friction.Additionally, the chosen orthogonal basis functions allow us to conduct power spectral analyses of the coefficients obtained from SCHA to estimate the fractal dimensions of the surfaces.
To adapt the SCHA for these applications, a new conformal parameterisation algorithm with minimal area distortion is herein proposed for parameterising the surface patches over a unit spherical cap.The parameterisation assures a unique assignment for each vertex and a uniform distribution of the morphological features over the spherical cap.Initially, the θ≤θc with half-angle θ c , where the spherical cap harmonics analysis (SCHA) is defined over.algorithm parameterises the surface patch to a planar disk using the conformal mapping algorithm proposed by Choi and Lui (2015) [5] followed by a suitable rescaling.Then, using the south-pole inverse stereographic projection, we conformally map the disk to a polar spherical cap.Finally, we find an optimal Möbius transformation to minimise the area distortion on the conformal spherical cap.
This paper is structured as follows: We first provide a brief literature review in Section 2. We then introduce the SCHA basis functions and their solutions in Section 3. In Section 4, we describe the proposed conformal mapping algorithm and then demonstrate the overall SCHA process in Section 5. Next, we explain the shape descriptors and the derived fractal dimensions (FD) in Section 6.We solve the problem of finding the optimal θ c for the analysis in Section 7. In Section 8, we explain how to modify the morphology by a proposed roughness projection methodology for the microstructures.The complete results are depicted in Section 9. Finally, we summarise the main results of the paper and discuss possible future works in Sections 10 and 11, respectively.

Literature review
We herein briefly present key publications on SHA and its application in mechanics, generation and reconstruction methods for both closed objects and those with free edges as well as SCHA and its current applications.
method that provides a unique mapping of each vertex on a surface onto a unit sphere S 2 through the use of equivalent coordinates x(θ, φ), y(θ, φ) and z(θ, φ).In other words, their method can describe 3D, closed and randomly-shaped objects with a relatively small number of shape descriptors derived from 2D signals mapped onto a unit sphere.The SHA method can therefore be considered a 3D generalisation of the traditional Elliptic Fourier Descriptors (EFDs) method (see [7,8,9]), which was developed to model 2D shapes (1D signals).The SHA method can also be used to accurately reconstruct 3D shapes from their shape descriptors.Recently, Su et al. [10] found a correlation between the SHA and EFD descriptors and used that for predicting the size of the irregular particles from 2D images.
Since the original SHA work was published, many contributions generalised it for a broader range of applications.While the original SHA dealt with voxel-based and simply connected meshes, the Control of Area and Length Distortions (CALD) algorithm [11] extended the parameterisation method to triangulated meshes, such as Standard Tessellation Language (STL) meshes.The original SPHARM determines the coefficients of the matrix by a least-squares algorithm that requires the calculation of the inverse.Modelling complex shapes with small details requires such a vast number of vertices that calculating the inverse is burdensome with the traditional SHA approach.Shen and Chung (2006) [12] used an iterative approach to expand the coefficients for the complex and large-scale modelling of shapes.Other notable works include the weighted SPHARM method by Chung et al. (2007) [13], the SPHARM-PDM statistical toolbox for shape descriptors by Styner et al. (2006) [14], and the 3D analytic framework with improved spherical parameterisation and SPHARM registration algorithms developed by Shen et al. (2009) [15].
In biomedical imaging, shape descriptors based on SHA were used to reconstruct surfaces and compare the morphology of organs.In engineering mechanics, two further questions were investigated: (i) How can SHA shape descriptors be linked to the well-known, traditional shape descriptors?(ii) How can the SHA shape descriptors be linked to mechanical parameters, such as friction coefficients, and traditional morphology measures, such as fractal dimensions (FD)?Here, we review studies in mechanics that used SHA for generating particles as well as for simulating particle mechanics.
Zhou et al. (2015) [16] scanned real particles, computed the distribution of the SHA shape descriptors and then used these descriptors as inputs for a particle generator.The particle generator used the shape descriptors to reconstruct particle shapes, which became inputs for Discrete Element Method (DEM) simulations.Wei et al. (2018a) [17] established links between the coefficients expanded from the SHA and traditional shape descriptors, namely, the form, roundness and compactness, and used these as inputs for a particle generator.They applied their method to Leighton Buzzard sand particles, which were scanned by Zhao et al. (2017) [18] to reveal an exponential decay of the rotation-invariant descriptors against the degree of expansion on a log-log scale.In Wei et al. (2018b) [19], the authors estimated the Hurst exponent (H) from the exponential decay of the power spectrum computed from SHA coefficients; thus the SHA is able to capture a self-similar phenomenon akin to the traditional Power Spectral Density (PSD) analysis associated with Fourier expansion for 1D signals.Apart from methods based on power spectral analysis, the authors in [20,21] correlated the FD in SHA by investigating the surface area evolution at an increasing reconstruction degree in a log-log scale.
More recently, Wei et al. (2020) [22] studied the contact between rough spheres and a smooth base using finite element analysis.For this, they constructed roughened spheres using spherical harmonics and considering more than 2000 degrees.The locally roughened patch on the sphere was generated by assuming a certain FD and relative roughness measure.
In contrast with the traditional concept of computing Fourier coefficients along certain lines (see Russ [1994] [23] for details), the advantage of SHA-based fractal dimensions is that they are not direction-dependent because they are intrinsically defined everywhere along the sphere and do not expand sampled sections (1D signals).Furthermore, the SH basis functions constitute a complete orthogonal set integrated over the entire surface of the sphere in Hilbert space S 2 , which allows the power spectral analysis.Lately, Feng (2021) [24] proposed a new energy-conserving contact model for triangulated star-shaped particles reconstructed with the SH.The model reduces the computational complexity of the discrete element simulations with controllable mesh size obtained via the golden spiral lattice algorithm.
The SHA is therefore a very powerful tool for contact analysis.However, it is currently only valid for closed objects.In this paper, we address this by offering an alternative method for analysing the surface morphology-instead of analysing the surface of an entire closed object using SHA, we study the morphology of a patch on this object using SCHA.Depending on the size of the patch on a sphere, this usually requires significantly less degrees.In addition, this method makes it possible to analyse regions of different roughness on the same closed object.Here, we provide a general procedure suitable for nominally flat and open surface patches and further show how the fractal dimension concept is an intrinsic property in such basis functions defined over spheres, hemispheres and spherical caps.

Modelling arbitrarily shaped objects with free edges
For modelling arbitrarily shaped objects using SHA, most of the literature relates exclusively to genus-0 closed objects.SHA requires an appropriate parameterisation algorithm (bijective function) over a unit sphere S 2 , but does not cover open objects, which have a free boundary along their edge.This problem was previously addressed using a hemispherical harmonic (HSH) basis instead of SH basis.The HSH basis uses shifted associate Legendre polynomials to describe a 2D function projected onto a hemisphere S 2 θc≤π/2 with a free edge.Huang et al. (2006) [25] was the first to apply HSH functions to extract shape descriptors for anatomical structures with free edges via a modified parameterisation algorithm of the original one proposed by [6].By applying the proposed method to images of ventricles, they concluded that it is more robust and natural than the SH basis for representing free edges, which requires more effort to satisfactorily converge along the discontinuities.In addition to having the naturally defined free edge in HSH, the convergence is also faster than with the ordinary SH because of the distortion introduced by the parameterisation algorithm and the maximum wavelength defined on the regionally-selected patch on the targeted object, which will be discussed and demonstrated in later sections.As a similar strategy, Giri et al. (2020) [26] recently developed a hemispherical area-preserving parameterisation algorithm and analysed 60 hemispherical-like anatomical surfaces.They compared the results with ordinary SH and also found significant improvements from 50 to 80% on the reconstruction quality.They also proposed variants of HSH for other anatomical shapes [27,28].

Regional modelling of spatial frequencies on a sphere
This paper uses SCH basis functions to study the local roughness and morphological features of rough surface patches that are regionally distributed over, for example, natural rocks, stones, bones, faults, etc.The idea of applying an SCHA to a regional subset of data distributed over a sphere was previously exploited in geophysics to study physical quantities such as geomagnetic and gravitational fields (refer to [29] for a comprehensive review of other applications).Many other methods can also be found in the literature for dealing with regional data.
Multitaper spectral analysis is one common method for localising data using global basis functions (integrated over the whole domain) such as the spherical harmonics.In this method, the data outside the region of interest is tapered to zero.When the region of interest is a patch with a free edge, it can be parameterised over an appropriate cap with an assumption of a grid-of-zeros over the rest of the sphere.This technique can be implemented efficiently using sparse matrices; however, the data tapering will affect the analysis results [30], and special relationships will be needed to link the localised analysis to the global one [31].
Another common approach is wavelet analysis, where the signals are multiplied (inner product) with a wavelet function and then integrated over the entire domain.However, this analysis becomes a function of the wavelet size (scale) and not the degrees of the SH expansion [30].Although it is tempting to use these approaches, and their computation is fast and stable, they alter the actual signals and properties of the basis functions.Instead, we herein consider a direct and analytical approach that defines global functions over the entire regional domain without altering the actual data or basis functions.

Spherical cap harmonic (SCH) basis functions
The SCH basis functions were first proposed by Haines (1985) [32], and his codes with the implementation details were published shortly thereafter [33].In his paper, he solved the SturmLiouville (self-adjoint) boundary value problem to derive two mutually orthogonal sets over a constrained spherical cap that satisfy the Laplacian.By applying one of the Neumann or Dirichlet boundary conditions on the free edge of the cap, we can get a general form of the associated Legendre function of real fractional (non-integer) degrees and integer orders.To satisfy the applied boundary conditions, we first solve for eigenvalues of the SturmLiouville boundary value problem, and these real and non-integer degrees of the associated Legendre functions can be used to define the basis functions.The resulting basis set depends on the applied boundary conditions; two types of bases are widely used and are named the "even" and the "odd" basis functions (for more details, see Section 3).
Haines's basis functions have been studied and revised several times.In general, they were derived to analyse signals on a solid spherical cap with a finite thickness, which is described by an internal (r = a) and external (r = b) radii.These basis functions were used to expand physical quantities through the Lithosphere.However, in this work, we use the basis for a spherical cap with zero thickness, i.e. a → b.This simplifies the problem because the normalisation of the spherical cone is ignored and the associated extra boundary conditions do not need to be applied through the thickness b − a. Haines's derivation of the normalisation factor was indeed an approximation (Schmidt-normalised), and he claims that the normalisation itself is not crucial as soon as the same normalisation factor is used for the reconstruction.An analytical derivation for the normalisation factor was proposed by Hwang et al. (1997) [34], and they named it the fully normalised SCH.
To capture smaller minimum wavelengths at the surface of the earth than with the original SCH, De Santis (1991) [35] proposed a modified version developed by shifting the origin of the spherical cap from the origin of the earth towards the surface.In 1992, De Santis [36] proposed using the ordinary SH defined over a hemisphere (using integer degrees) by scaling the spherical cap data to fit over a hemisphere.Although it is tempting to use such approximations to relax the complexity of the proposed functions in the original SCH, they are only valid for small half-angles (shallow spherical caps).More recently, Thébault and his co-authors [37,38] proposed a new revised version of the SCH to deal with problems mainly associated with the solid spherical cone.These papers are considered out-of-scope here because we only deal with the surfaces of spherical caps.More details about the history and the evolution of the method and its applications in geophysics can be found in [29,39].

Parameterisation algorithms
Before expanding spatial data using the SCH basis, the surface patch needs to be mapped onto the spherical cap.This mapping algorithm assigns each vertex coordinate of the surface of the original object a unique coordinate on the spherical cap.A proper parameterisation algorithm should distribute the morphological features uniformly over the surface without clustering [6].The mapping algorithms used here can be categorised into angle-preserving (conformal) [40,41] or area-preserving [42,43].The former preserves the angles between the mesh elements such that that the local morphological properties are preserved, but this could introduce undesirable distortion in the element areas.The latter algorithm preserves the area of the elements but can introduce sufficient angular distortion to jeopardise the representation of the local morphological features.In general, angle and area preservation cannot be achieved at the same time [44], hence the choice of parameterisation method depends on the application.In mesh fitting, it is often desirable to balance angle and area distortions [45].

Spherical cap harmonic (SCH) basis functions
In this section, we describe the set of orthogonal functions on the surface of the spherical cap that we use for studying the geometry of a nominally flat and open-edged patch projected onto the cap.To do this, we propose using the basis functions of the SCHs that satisfy the Laplace equation in 3D; the derivations for these functions can be found in [32,34].Here, we introduce the derivation only to the extent that the definition and meaning of the various components of the functions are clear.It serves also as documentation for the accompanying code, which is shared openly and in which these equations have been implemented numerically.
Here, we are considering only polar spherical caps, meaning the local polar axis of the cap coincides with the global polar axis of the unit sphere.The coordinates of any point on the unit spherical cap are defined in terms of the angles θ and φ.The angle θ represents the latitude (elevation angle), measured from the positive z-axis ∈ [0, θ c ], and the angle φ represents the azimuth angle (longitude), measured counter-clockwise from the positive xaxis ∈ [0, 2π].The size of the spherical cap S 2 θ≤θc is therefore characterised by the half-angle of the spherical cap θ c , which is the latitude of the rim of the spherical cap as illustrated in Figure 1.
A 2D signal f (θ, φ), for which each value is paired with (θ, φ), can be written using the SCH basis as (see [32]): where θc C m k (θ, φ) are the SCHs and θc q m k are specific constants (SCH coefficients).The complex-valued SCHs of degree l and order m for a cap of size θ c can be expressed as: One way of solving the Laplace equation is by separating the independent variables, θ and φ in our problem, into separate ordinary differential equations and then assuming the overall solution to be the product of the separated solutions.The SCHs can therefore be written as a product of the term that describes the variation with regard to the latitude θ and a term that describes the variation with regard to the longitude φ.In polar caps, the longitudinal direction (φ) of the spherical cap covers, by definition, complete circles (Parallels of Latitude), so we apply the traditional Fourier method to capture the circular variations.
The variations along the latitudes, represented by (θ), are captured using, for example, associated Legendre functions (ALF) as a mutually orthogonal set of basis functions.In the following subsections, we describe the ALFs, how to stably compute them and how the degrees ls of these ALFs are found.Then we explain one possible normalisation factor that we used in this paper.

Associated Legendre function (ALF)
In classical SH, the degree l and the order m are integers, and the associated Legendre polynomials (ALPs) are defined for the interval x ∈ [−1, 1].In this interval, they form an orthogonal set of equations on a unit sphere (for x = cos θ, θ ∈ [0, π]).In the HSH basis functions, a linear transformation shifts the orthogonality interval to either [0, 1] or [−1, 0] depending on whether the lower S 2 θ≥π/2 or upper S 2 θ≤π/2 hemisphere is considered (see the works in [25,26]).The hemisphere introduces a new boundary (free edge) at the equator that must be defined with an additional appropriate boundary condition.However, note that the tangent at the edge of the equator (the gradient with respect to θ) is zero.As will be seen next, this is one of the Neumann boundary conditions for the boundary value problem of ALFs, which makes it a special case of the SCH basis with integer degrees and orders.The main challenges with SCHs are shifting the ALFs to be defined over the interval [cos θ c , 1] and assuring their orthogonality.
Using SH or HSH functions to describe a spherical cap (sub-region over the whole valid domain) will produce a set of ill-conditioned equations [34].Although these sets could be used to some extent to reconstruct the input surfaces, the power spectral analysis will be chaotic and without the clear attenuating trends normally produced when the equations are well-conditioned.To find a suitable set of basis functions, Haines [32] proposed treating the problem as an inverse problem and determining the degrees l for the ALPs that produce an orthogonal basis over the spherical cap while satisfying the boundary conditions.These degrees are not constrained to integers, and in our case, the degrees will be real-valued numbers.The orders m are still integers because, for polar caps, φ always varies over complete circles.
The basis functions proposed by Haines (1985) [32] satisfy three boundary conditions: (i) continuity along the longitudes, (ii) regularity at the pole of the cap, and (iii) an arbitrary signal value at the free-edge (see [46,29]).The first boundary condition is assured by the use of polar spherical caps (notice the complete circles in Figure 1), as they force the periodicity and the continuity along the longitudes to make the order m a positive or negative integer: The second boundary condition, the regularity at the cap's pole θ = 0, is satisfied through the use of the first kind of ALFs, such that: ∂f Note that even with the two above boundary conditions, the solution functions are still not unique.To obtain unique solution functions, we need to impose a proper boundary condition at the edge of the spherical cap.The third boundary condition is applied to the free edge at θ c (x c = cos θ c ): Here, we use the degree notation l(m) k instead of l because the degrees are now back calculated for a specific order m and index k, which will be explained in detail below.In Eq. ( 7), the case of A c = 0 and B c = 0 corresponds to the application of the Neumann boundary condition where the gradient with respect to θ at the cap's rim at x = cos θ is zero.It implies that f (θ c , φ) = 0 ∀φ ∈ [0, 2π] and can take any arbitrary value at the edge of the spherical cap.Applying this boundary condition results in basis functions referred to by Haines (1985) [32] as the "even" (k − m = even) basis or by Hwang et al. (1997) [34] as set 2. The case of A c = 0 and B c = 0 corresponds to the Dirichlet boundary conditions f (θ c , φ) = 0 ∀φ ∈ [0, 2π], and the resulting basis functions are known as the "odd" (k −m = odd) set or set 1 by Haines [32] or Hwang et al. [34], respectively.The boundary conditions for the two cases can therefore be rewritten as: The latter boundary condition, Eq. ( 8) or ( 9), is met by back calculating the roots l(m) k at different orders m and indices k and then using them for constructing the basis functions for different θ c .Although it is tempting to use a combination of both sets to expand a signal (e.g.see Haines [32]), for the uniform convergence of the series, it is noteworthy that these basis functions are mutually orthogonal and cannot be used together for the power spectral analysis [47].In this paper, we will only consider the even set, as it allows arbitrary values for the signal at the cap's rim (f (θ c , φ) = 0 ∀φ ∈ [0, 2π]), which contrasts with the odd set that equates the signal at the rim to zero.This consideration therefore allows us to correctly reconstruct any surface geometry.
Unlike the traditional ALPs, here the degrees l(m) k are determined such that the third boundary condition is satisfied; thus, l(m) k are real numbers, though not necessarily integers; hence, we use the index k to rank the orders.k is an integer index that arranges m from 0 to k, and the degree l(m) k is the degree at the k th index and order m.The k index is useful for finding the degrees on the applied boundary conditions.For instance, when k = 0, we find the first degree l(m) k ≥ 0 where the boundary condition equation first crosses zero (first eigenvalue).
For constructing the ALFs with integer orders and real degrees, Hobson (1965) [48] proposed the use of hypergeometric functions.Note that in the following, we refer to associated Legendre functions rather than polynomials because these basis functions are truncated numerically at a predefined numerical threshold.These functions can also be referred to as fractional ALFs, because the associated degrees are not integers.The ALFs can then be written as follows (see [32,34,48]): where Γ(•) is the gamma function, 2 F 1 (a, b; c; z) is the hypergeometric function defined over a unit disk |Z| < 1 and x = cos θ.
To find the zeros (eigenvalues) that satisfy the third boundary condition, we rewrite Legendre's differential equation in a Sturm-Liouville (S-L) form [32,34]: and subject it to the boundary conditions explained in Eq. ( 7).The two equations that result from the application of the Neumann and Dirichlet boundary conditions can then be written in a simpler factorised form, as proposed by Hwang et al. (1997) [34] for hypergeometric functions: where Equation ( 12) is the equivalent of applying the Neumann boundary conditions (even set), and Eq. ( 13) corresponds to the Dirichlet boundary conditions (odd set).Both Eqs. ( 12) and ( 13) have an infinite number of roots (eigenvalues) for the abovementioned S-L problem.Here, we are interested in finding the first k + 1 eigenvalues l(m) k where k ∈ {0, . . ., |m|}, which is done using a simple stepping solver based on Mueller's method for finding zeros (roots) [49].Additional results about the roots for θ c = 5π/18 are shown in Appendix 13.4.
After finding the needed roots for the different m and k, we can substitute these roots directly into Eq.( 10) to obtain the needed ALF basis.Figure 2 shows the ALFs associated with the roots for different k, where m = 0 and θ c = 10 • = π/18 (i.e.cos(π/18) = 0.985).Graphically, the value k − m indicates how many times the basis function changes signs along the interval [cos θ c , 1]; the (k + 1) th eigenvalue changes signs once more than the k th root, as it divides the interval into regions with different signs.Together, k and m define the frequency classification (harmonics) of the basis functions as zonal (m = 0), sectoral (|m| = k) or tesseral (0 < |m| < k) regions on the SCH, as will be illustrated next.Figure  12); (B) the odd basis and the corresponding eigenvalues from Eq. ( 13).
2 also shows the difference between even and odd sets depending on the applied boundaries at x = cos θ c .Equations ( 10), ( 12) and ( 13) are dependent on the ordinary Gaussian hypergeometric function evaluation 2 F 1 (a, b; c; z).Unfortunately, the implementations given in [32,34] using the recursive series are not numerically stable for large degrees when k > 12 about the radius of convergence |Z| < 1. Indeed this is a well-known problem for such functions, and different formulations can be found in the literature with different implementations for different ranges of a, b, c and z (see Pearson et al. [2015] [50], who reviewed different hypergeometric functions and their different implementations and stability conditions).Therefore, we herein use a Taylor series expansion for the first degrees (up to k = 12) when a, b and c are relatively small.Otherwise, we use a 5 th order Runge-Kutta algorithm along with the Dormand-Prince method for solving the hypergeometric differential equation (see [50] for more details).

Normalising the SCH basis functions
Normalising the basis functions is important for ensuring their orthonormality and for keeping the calculation of the coefficients computationally manageable [34,46].Many normalisation methods have been proposed in the literature, such as the "Heiskanen and Moritz" method [34] and the "Schmidt" or "Neumann" methods [48], and any of these can be used as long as the same normalisation method is used for both the expansion and the reconstruction [32].In the following, we use the "Schmidt semi-normalised" method by Haines [32].
This normalisation method requires solving the integral for the ALFs over the spherical cap (mean square product), which can be retrieved from [51] with a corrected sign from [52]: However, this complex expression is difficult to solve (see [34] for one approach using recursive expressions).In this paper, we use the Schmidt semi-normalised harmonics introduced by Haines (1985) [32]: where with The final form of the normalised ALFs can be written as follows (see [32]): We use Eq. ( 21) as a part of the final SCH basis functions.

The SCH series
The complete form of the SCH can then be expressed by combining Legendre's functions with Fourier expansions: We can split this expression into a real part and an imaginary part: C m k (θ, φ) .Note that at the breathing mode when k = m = 0, the surface inflation is unity and is therefore not included in the colorbar.The complex part of the basis are similar but rotated π/(2|m|) about the z-axis.
The complex-valued equation ( 22) is only valid for non-negative orders |m|.To address this, the redundant negative-orders equation can be written in terms of the positive ones, as they reassemble π/2 rotations about the positive z-axis.By making use of the Condon-Shortly phase factor, we can write the SCH functions for m < 0 as: where C m * k (θ, φ) is the complex conjugate of the SCH functions with positive orders.Figure 3 shows the real SCH functions (with positive and negative orders) for the half-angle θ c = 5π/18 and up to k = 4.
We can write a bi-directional signal f (θ, φ) that is explicitly parameterised to unique θ ∈ [0, θ c ] and φ ∈ [0, 2π] over a unit spherical cap S 2 θ≤θc as a linear sum of the independent harmonics defined over the spherical cap: where θc q m k are the complex-valued SCH coefficients at order m and index k.In our case, following Brechbühler's work in [6], the signal f (θ, φ) is the 3D vector that comprises three orthogonal components in the Cartesian space extracted from the coordinates of the vertices: In practice, we truncate the series in Eq. ( 25) at k = K max , which reassembles the minimumaccumulated wavelength (equivalent to the cut-off frequency) ω min in the harmonic sum from the original surface patch.Then, f (θ, φ) can be approximated by: The coefficients θc q m k can be expanded as: However, this integral is valid only when the signal f (θ, φ) is defined everywhere over the spherical cap (continuous signal).For discrete signals, we can estimate the coefficients using least-square fitting algorithms, as will be discussed later in Section 5.

Spherical cap conformal parameterisation
Here we develop a conformal parameterisation method for mapping any surface S with disk topology to the spherical cap S 2 θ≤θc , where θ c is the prescribed half angle.The overall spherical cap parameterisation f : S → S 2 θ≤θc is desired to be not only conformal but also with low area distortion.
First, note that the south-pole stereographic projection τ given by maps the spherical cap S 2 θ≤θc to a planar disk D r centred at the origin with radius and the inverse south-pole stereographic projection τ −1 given by maps D r to S 2 θ≤θc .As both τ and τ −1 are conformal, we were motivated to find an optimal conformal mapping from the given surface S to the planar disk D r so that we could apply the inverse projection τ −1 to obtain the desired parameterisation.4.

PSS
To achieve this, we first apply the disk conformal mapping method [5] to map S onto the unit disk D in an angle-preserving manner (see [5] for the algorithmic details).Denote the conformal mapping by g : S → D. Now, we search for a conformal mapping h : D → D r that can further reduce the distortion in area of the parameterisation.Mathematically, a Möbius transformation is a function that maps a complex number z = X + Y i (where X, Y are real numbers and i is the imaginary number with i 2 = −1) to another complex number (az + b)/(cz + d), where a, b, c, d are complex coefficients with ad − bc = 0. Möbius transformations are conformal and hence are a good candidate for the mapping h in our problem.More specifically, here we consider a Möbius transformation h : where A, B are two real numbers to be determined.To find the optimal h, we minimise the following area distortion measure where F is the set of all triangular faces of S. The desired spherical cap parameterisation is then given by In particular, since every step described above is conformal, the overall mapping f is also conformal.
To explain the formulation of the area distortion measure d area in Eq. ( 34) more clearly, we first consider the mathematical definition of area distortion.For any two surfaces M, N in the Euclidean space R 3 with the same total surface area, a mapping f : M → N is said to be area-preserving if for every open set U on M, the surface area of U is equal to the surface area of f (U ).More generally, the area distortion of a mapping f can therefore be evaluated using the dimensionless quantity log Area(f (U )) Area(U ) , which equals 0 if and only if Area(f (U )) = Area(U ).In our problem, M corresponds to the input triangulated surface S and hence it is natural to evaluate the area distortion by considering every triangular face T .In other words, by minimising Eq. ( 34) we look for a Möbius transformation h such that the area change of every triangular face under the spherical cap conformal parameterisation is as small as possible.Moreover, note that the total surface area of the input surface S is not necessarily equal to that of the target spherical cap domain S 2 θ≤θc .Therefore, we use the two summation terms in Eq. ( 34) for normalising the total surface area of the input surface and that of the spherical cap in the area distortion measure.
In practice, the solution of this objective function can be found using a global search heuristic such as the Pareto-like Sequential Sampling (PSS) approach (see [54] for details).Using a population size of 30 and an acceptance rate of α = 0.97, only 20 iterations are needed for most of the problems to reach a good approximation.This approach is gradientfree and does not require an initial value that could converge locally.Table 1 shows the stability of the PSS approach, while Figure 4 shows the convergence and the final result of the parameterisation step.Alternatively, the minimisation problem can also be solved using MATLAB's fmincon, using the modulus and argument of the complex number A + Bi as variable.The solution obtained from this function was d area = 1.621523941729933 with X = 0.032055812166766 and Y = 2.416521766598053, which is consistent with the PSS result.Figure 5 shows several spherical cap parameterisation results with different prescribed θ c for the Stanford bunny model [55].To further assess the conformality of the parameterisation results, we define the angle distortion d angle as the average of the absolute difference between every angle [v i , v j , v k ] in a triangular face in the input surface and that in the spherical cap parameterisation f : The values of d angle for the parameterisations in Figure 5(C), (D), and (E) are 0.00828069, 0.00869725, and 0.00919996 respectively.The results show that the proposed parameterisation method is highly conformal.
5 Numerical implementation of the spherical cap harmonic analysis (SCHA) In this section, we describe the numerical implementation of the SCHA method and use it to reconstruct triangulated surface meshes.

Spherical cap harmonic analysis (SCHA)
SCHA comprises the computation of the weights θc q m k by the integral in Eq. ( 29).However, this integral cannot be calculated directly as the data is not defined everywhere over the spherical cap.Instead, it can be estimated for discrete points that are distributed randomly over the cap.For this reason, we estimate these coefficients using the least squares method (statistical regularisation).Following the indexing system proposed in [6], the coefficients can be estimated as follows: where B contains the basis functions evaluated at the coordinates V that define the geometry of the surface patch.If n v vertices describe the geometry of this surface patch in spatial coordinates, the matrix V can be written as: It should be noticed that the vertices defined in Eq. ( 38) are normalised against the mean, so the actual expanded signals will be identified about a mean of zero.The matrix B holds the SCH basis functions that can be obtained by evaluating Eqs. ( 22) and ( 24) for the positive and negative orders, respectively, and can be written as: where c i,j is the basis function evaluated at the i th vertex for index j.The index j is calculated as j(m, k) := k 2 + k + m and −k ≤ m ≤ k for the zero-based numbering programming convention.This means we need to fit 3(K max +1) 2 coefficients simultaneously for all the x(θ, φ), y(θ, φ) and z(θ, φ) signals in f (θ, φ).The final complex-valued θc q m k coefficient matrix can be written as: The least squares fitting is a statistical regularisation method that provides a good approximations for the coefficients θc q m k , as long as the number of equations on the spherical cap is significantly larger than the expanded coefficients.Here, this condition is always satisfied because our surfaces are defined by dense point clouds.Other regularisation methods that enhance the accuracy of the least squares fitting and obtain a better orthogonality among the columns in Eq. ( 39) can be found in the literature (e.g.see [56]) but are not applied here.

Spherical cap harmonic reconstruction
We can reconstruct the surface patch from Eq. ( 27) when the coefficients are known.To obtain a homogeneous representation of all the morphological features over the cap domain, this step requires uniformly sampled points (vertices).In ordinary SHA, convex icosahedrons (geodesic polyhedrons) with different mesh refinements are often used for the reconstruction because the vertices of these icosahedrons are equally spaced.Analogously, for SCHA, we use geodesic domes with uniformly distributed vertices over a spherical cap (see Figure 7).
To construct these domes, we first construct a uniformly meshed unit disk and then use the inverse Lambert azimuthal equal-area projection to project it to the desired spherical cap.Like the approach presented in Section 4, we need to rescale the unit disk by an appropriate scaling factor r l to yield a spherical cap with the correct area under the projection: For reconstructing the surface patch, we need to choose a suitable mesh refinement for practical and theoretical aspects.Practically, the complexity of the numerical models, e.g.finite and discrete element methods, depends on the number of elements used to discretise the continuum.Theoretically, we need a mesh size that represents the smallest required detail depending on the application.If we imagine that a simple 1D wave propagates along a discrete line, we need at least five knots (points) to capture the sign changes along a full wave.This means that the distance between two adjacent vertices at a certain mesh refinement must be at least one-fourth of the target minimum wavelength.The trade-off between Figure 7: Geodesic domes for θ c = π/3 with different refinement cycles (N ).The domes in A-C were obtained using the unit disk mesh generation method in [57] followed by the rescaling in Equation ( 41) and inverse Lambert's projection in Equation (42).The domes in D-E were obtained using the DistMesh method [58] followed by rescaling and the inverse Lambert's projection.A) A mesh grid resolution of 12 results in 36 vertices; B) A mesh grid resolution of 20 results in 100 vertices; C) A mesh grid resolution of 60 results in 900 vertices; D) An edge element size of 0.325 results in 35 vertices; E) An edge element size of 0.185 results 103 vertices; F) An edge element size of 0.635 results in 896 vertices.
capturing morphological features and the numerical complexity should be calibrated through several reconstruction trials.More details about the minimum and maximum wavelength associated with the SCHA will be discussed in Section 7.1.
We consider two approaches for constructing a disk with a uniform grid.The first approach, proposed by Roşca (2010) [57], builds the disk-grid by transforming a rectangular grid of squares to a quad-grid disk using an area-preserving map, which can be subsequently triangulated using standard Delaunay triangulation.This approach gives us direct control of the number of vertices on the spherical cap.Another approach is the DistMesh method proposed by Persson and Strang (2004) [58], in which one can construct a uniform mesh on a unit disk by a force-based smoothing approach.For this approach, the desired length of the edge elements on the unit disk is a required input that indirectly controls the number of vertices on the surface mesh.For both approaches, we need to rescale the resulting unit disk meshes using Eq. ( 41) (the rescaled disk is denoted by D l ) before applying the inverse Lambert's projection τ −1 l so that every (x, y) point on the rescaled disk is projected to the (X, Y, Z) on a spherical cap to S 2 θ≤θc as: Figure 6 visualises the process in [57], and Figure 7 shows the final geodesic domes obtained using the two approaches.

Shape descriptors and fractal dimension (FD)
Like the traditional Elliptical Fourier Descriptors (EFD) [48] and their 3D extension to SH [6], SCHA also yields coefficients that can describe and reconstruct surfaces.Because the SCH is a "local extension" of the ordinary SH, it matches most of its properties [47], such as the spectral power analysis for determining the surface fractality.In this section, we describe the interpretation of the shape descriptors derived from SCHA and their relationship to those derived from other traditional methods like HSH.We also show how to estimate the wavelengths associated with each degree and extract the fractal dimension (FD) from the attenuation of the shape descriptors.

Shape descriptors derived from SCHA
Shape descriptors tell us how similar or dissimilar surfaces are to each other and are particularly useful for matching or classifying complex topologies.These shape descriptors must be invariant to the location, rotation and size of the targeted surface.In this section, we compare the shape descriptors derived from SCHA to those derived from SHA or EFD and explain some additional properties of the SCH regional analysis.By expanding the signal f (θ, φ) = (x(θ, φ), y(θ, φ), z(θ, φ)) T , we simultaneously obtain the coefficients for the three orthogonal signals.We first consider the "breathing" mode (corresponding to k = 0), where we have a constant shift for all the points in the signal.The corresponding real-valued (k = m = 0) coefficients along the x, y and z axes are represented by θc q 0 0,x , θc q 0 0,y and θc q 0 0,z .These coefficients correspond to the geometric centroid of the surface, which is similar to the EFD and SHA methods.Practically, we ignore these coefficients and set them to zero to make the shape descriptors and the reconstruction invariant to the location.
For k = 1, EFD basis functions define an ellipse in 2D, and SH functions define an ellipsoid in 3D.However, in our case, the basis functions describe an ellipsoidal cap that is sized by the 3×3 matrix θc q m 1 , m ∈ {−1, 0, 1}, herein referred to as the first-degree ellipsoidal cap (FDEC).The FDEC determines the size and the orientation of the reconstructed surface, and it also shapes the main domain (spatial space) of the other harmonics (for all |m| > 1).Following Brechbühler et al. [6], the size of the FDEC can be determined by rearranging the coefficients from k = 1 to correspond to the ellipsoidal equation in the Cartesian coordinates: Solving the eigenvalue problem of AA T will provide three eigenvalues correspond to the amount of "stretch" along the three principal axes of the ellipsoidal cap.The directions of the "stretches" will follow three eigenvectors; two of them are in-plane orthogonal vectors and the third is an out-of-plane perpendicular vector.These three unit vectors describe the principal directions of the ellipsoidal cap. Figure 8 shows the FDEC in an arbitrarily rotated Cartesian space.When the size of the ellipsoidal cap is determined, the higher frequency basis functions will be overlaid onto the domain of the FDEC.Numerically, the basis functions of the SCH are similar for different θ c 's, but they are either "stretched" or "compressed" over S 2 θ≤θc .This constrains all the domains to the size of the FDEC, which makes the proposed method invariant to θ c .From the point of view of the basis functions, this theoretically implies that there is no difference between the SCHA and the HSHA, providing that we are using the same normalisation factor (e.g.Schmidt semi-normalised).However, the results of the analysis will depend on the introduced distortion between the real Cartesian space and the parameterisation space.This problem will be discussed in Section 7.
The FDEC size is also important for estimating the wavelengths associated with the harmonics expanded over the surfaces.In geophysics, the SCHA is used for a regional analysis over planets and moons to directly estimate the average circumference.As we do not have a closed surface that is approximately spherical, however, we use the FDEC to estimate the average circumference of a full "virtual" ellipse that lies at θ c ≡ π/2 (Greenwich line) to estimate the circumference (more details will be discussed in Section 7).

Shape descriptors for estimating the fractal dimension
Shape descriptors have been extensively studied in the SH literature.The most commonly used descriptors are rotation-invariant ones that do not require the registration of the surface.Based on the amplitudes calculated in the expansion stage θc q m k , we can express the shape descriptors (for instance see [18,59,60]) as follows: Then, the resultant (radial) rotation-invariant descriptors can be expressed as: These descriptors need to be normalised with the size of the FDEC to obtain size-invariant descriptors and to estimate the Hurst exponent, important for removing the dimensionality of the data [23].Thus, the final shape descriptors can be written as: For nominally flat and circular surface patches, where the roughness (noise) only shifts the z−coordinates, D2 k,x ≈ D2 k,y → 0 for k > 1.The power spectral analysis determines the power accumulated at a certain wavelength of the analysed surface, while the spectral attenuation computes the FD of the surface.Similar to the PSD analysis, here we found that the accumulated power at each degree correlates with the Hurst exponent H (see [23,61] for details) as follows: The FD can then be computed as (see [23]): For fractal surfaces, the Hurst coefficient is between 0 and 1 (0 < H < 1), and the FD is between 2 and 3 (2 < F D < 3).The higher the H exponent, the smoother the surface.The Hurst exponent can be estimated from the rate of power attenuation with an increasing k, as the slope of the fitted line on a log-log graph is −1/(2H).One common way to estimate the exponent is through the least squares fit of the equation D 2 k = ak −2H , where a is the proportionality constant in Eq. ( 45).This approach estimates a radially symmetric Hurst exponent for self-affine surfaces.Unlike the Hurst exponents computed by traditional methods that investigate profiles along different directions, the Hurst exponent herein computed on the basis of SCH coefficients is invariant to directions or rotations.We benchmarked this approach against artificial fractal surfaces, and the results are presented in Section 9.3.
7 Choosing the optimal half-angle θ c and the appropriate surface patch 7.1 Choosing the optimal half-angle θ c In Section 6, we explained that the SCH basis is invariant to θ c , raising the question: What is the optimal half-angle θ c for the analysis?To answer this, we need to consider two aspects: i) the numerical stability of the hypergeometric functions and ii) the distortion between the object's space and the parameterisation space.The stability of the hypergeometric functions was briefly discussed in Section 3. In general, the smaller the θ c , the faster the convergence (less terms needed), but θ c should not be too close to the singularity points θ c ∈ {0, π}.
If this restraint is considered, we can choose the half-angle θ c such that the distortion introduced by the parameterisation is minimised.As compared to HSH, SCH offers this additional parameter θ c to minimise the distortion.We here adopt the area distortion (d area ) used in Eq. ( 34) to determine the optimal θ c , though other measures can also be used.
To find the optimal θ c that minimises the area distortion, we used the PSS algorithm [54], which we used already for finding the optimal Möbius transformation for h(X, Y ) (see Section 4).Note that we used the area distortion measure d area as the objective function as the proposed parameterisation method is highly conformal and always yields a minimal angular distortion regardless of θ c .We searched for the optimal θ c in a domain limited by a lower bound θ LB and an upper bound θ U B .The optimisation problem can then be written as: By investigating several benchmarks, we observed that the algorithm chooses the θ c that best represents the global concavity of the open surface.For instance, for the nominally flat surfaces examined in Section 9, we find that θ c is always as low as θ LB , which we set to π/18 = 10 • .For the portion extracted from the scanned stone (see Section 13.2), we obtain an optimal half angle of 99.9158 • , which corresponds to an area distortion of 0.240406.For practical reasons, we can assume θ c ≈ 90 • (a hemisphere) with an area distortion of 0.253300.In addition to the natural free edge, this helps explain why using HSH led to better results than SH when studying skulls and ventricles in medical-imaging problems [25,26].

Wavelength analysis with θ c
The literature review in Section 2.1 highlighted that separating the shape of an object from its roughness is a challenging problem, with most prior works proposing the use of empirical rules based on the reconstruction quality.We here offer an alternative approach adopted from the geophysics literature on SCH that determines the required number of degrees based on the targeted wavelength range.Bullard (1967) [62] defined the wavelength in SH by ω = 2πr/n, where n is the SH degree and r is the average radius of the sphere.Haines (1988) [33] derived the relationship between the degree and the spherical cap half-angle θ c using an asymptotic approximation (for large k and small m): where θ c is in radians.Analogous to [62], Haines concluded that the minimum spherical cap wavelength associated with the k max is ω min = 2πr/l (m) k,max .Solving Eq. ( 48) for k max gives: where k max is the maximum index required to cover or represent a minimum wavelength on the expanded cap.
In geophysics, the term 2πr is the average circumference of the earth (or any other planet under consideration).As we mentioned earlier and unlike in geophysics, here we expand three orthogonal signals instead of the radial one.The size of the domain is described by the three parameters a, b and c, which are determined by the FDEC for k = 1.Bullard (1967) [62] derived this expression from the number of divisions along the circumference for a certain degree n; the maximum circle defined over the sphere (equator) occurs when m = n (sectoral harmonics).Furthermore, he proved that this is true at any point on the sphere for any arbitrary pole.Here, we apply the same concept by replacing n with k (the index of the degrees).At m = k, the largest ellipse defined over the spherical cap changes its sign k times.Thus, ω = ζ/k, where ζ is the circumference of the largest ellipse defined over the first degree ellipsoidal cap (FDEC) (Figure 8).Additionally, ζ can be estimated from the FDEC, which is described by the orthogonal vectors θc q −1 1 , θc q 0 1 and θc q 1 1 (Section 6).The wavelength at an index k can then be written as: Equation ( 50) can also be reached by assuming that the largest ellipse on an ellipsoidal cap corresponds to the equator (runs through the Greenwich line), and thus θ c → π/2.Then, by letting r ≈ (a 2 + b 2 )/2, the asymptotic formula in Eq. ( 49) will be the same as in Eq. ( 50).Note that the (≈) sign in the above expressions denote that these formulae have been derived asymptotically [29], the circumference of the largest ellipse is also approximated, and the circumference of the FDEC is not the final circumference of the expanded surface patch.This equation is more accurate for circular patches where the FDEC perimeter is more circular than elliptical.

Choosing the appropriate surface patch
Choosing the appropriate sampled surface patch eases the parameterisation process.It also avoids over-/under-representation of certain morphological features over the targeted surface, accordingly preventing potential problems in the subsequent analysis step.The following points are recommendations for sampling a surface patch: • Manifold edges and surfaces without opening (genus-0) meshes are required for the parameterisation algorithm proposed in this paper.
• The input mesh must not contain duplicated faces, edges, vertices or skewed faces with close-to-zero areas.
• A round surface patch avoids point clusters at sharp corners as a result of the parameterisation algorithm, which would lead to a poorly represented morphology.
• Circular patches with no sharp edges will cause the first ellipse to be mostly circular.As a result, the wavelength definition will be more accurate because the circumference estimate is more precise.
• The scanned surface patch should be sampled as uniformly as possible, which will make the parameterisation algorithm more accurate due to the use of the least-squares fitting method for computing the coefficients of the SCHA.
• When the surface patch does not have enough vertices to compute all degrees (results in badly conditioned matrix B in Eq. [37]), the surface should be subdivided as needed to increase the number of expanded equations and increase the quality of the coefficients.
• When the results are wavy when reconstructed with high degrees, we recommend subdividing the overall surface or the local wavy areas to increase the number of expanded equations, thus reducing the fitting error.

Windowing the spectral domain and roughness projection
In this section, we introduce a method for projecting the roughness on surfaces using the analysis results expanded from the SCH, HSH and ordinary SH methods.This is useful for studies of contact problems, as it allows for systematically modifying the bandwidth of wavelengths included in the reconstruction and for altering a surface geometry.The isolated or generated surface roughness can then be projected onto a donor mesh of the user's choice.
Let us assume that we have obtained the SCH coefficients up to the order K max either through expansion or through artificial generation (not included in this paper).The rough surface is to be projected onto a mesh that donates the general shape of the object, henceforth called the "donor mesh".The targeted vertices on the donor mesh can be determined and parameterised to obtain For the SH, we combined the parameterisation approaches by Choi et al. in 2015 [63] and 2020 [64] to obtain a conformal spherical parameterisation of the donor mesh with minimal area distortion.For the HSH and the SCH, we used the approach proposed in Section 4.However, instead of parameterising only the spherical cap, we parameterised the whole donor mesh over a unit sphere and then constrained the targeted domain of the vertices to a spherical cap with θ c .
After parameterising the targeted vertices over the donor mesh, we projected the roughness onto the surface using a simple bandwidth-limiting (windowing) function in the spectral domain.Let us assume that kmin ≥ 1 and kmax ≤ K max are the minimum and maximum indices, which correspond to ω max and ω min wavelengths computed from Eq. ( 50).The roughness projection can then be implemented by modifying Eq. ( 27) to: where C m k,d (θ, φ) are the basis functions evaluated for the targeted domain on the donor mesh.This expression is equivalent to multiplying the reconstruction expression in Eq. ( 27) by the traditional rectangular eigenfunction window H(k) in the spectral domain (a sharp cut-off in the spectral domain) where: The coefficients can be manipulated in the spectral domain by rotating and scaling as needed.In this paper, we will not demonstrate the coefficients manipulation methods or the random generation of fractal surfaces with a certain Hurst exponent.

Numerical examples and discussions for SCHA
In this section, we present and discuss the results of the proposed SCHA method.We begin with a visual benchmark that helps build the intuitive arguments we raise about the scalability and wavelengths.Second, we benchmark the SCHA against various random fractal surfaces to estimate the Hurst exponent.Then, we analyse and study a rough surface patch of a rock whose geometry was obtained by laser scanning.Finally, we demonstrate an example of a roughness projection from the scanned stone to alter the microstructure of a donor mesh.

Visual benchmark
We use a sculpture of a face (retrieved from [65]) for visualising the reconstruction evolution with k.The 3D-scanned geometry of the face was first cleaned and remeshed to a manifold surface with a reasonable number of vertices that sufficiently describe the details. Figure 9(A) shows the input geometry that was used.This benchmark was chosen for the following reasons: • It is a nominally flat surface that contains large and small wavelengths (facial features).
To capture all the small details, a high number of orders is required.
• It has a significant number of vertices (more than 20, 000), which increases the problem complexity.
In this example, we used θ c = π/18 to parameterise, expand and reconstruct the surface, as this half angle is low enough to minimise the area distortion but also high enough for the hypergeometric functions to be stable.The results in Figure 9(B)-(H) show the evolution of the reconstruction stages with k.
The method starts with fitting the large wavelengths (low k's).The large facial features, such as the nose bone, begin to take shape at k = 3.The eyebrow ridges begin to take shape at k = 10.Fine details like the mouth and the hair curls become visible only with higher k's.One way to measure the convergence of the spatial details is by simply measuring the distance between the sampled points (Monte Carlo) on the input surface and the closest equivalent ones on the reconstructed surface at a certain k.The root mean square error (RMSE) of these points is then calculated using MeshLab [66] as the Euclidean distance normalised by the diagonal of the bounding box of the surface.The results are shown in Figure 10; for k = 40, the RMSE averaged over all points is 0.158282.Another measurement that can compare the reconstruction of the original shape is the Hausdorff distance; see the inset included in the same figure.It should be mentioned that when k = 40, the largest captured feature varied at a wavelength of 184.589 mm, and the smallest varied at 4.858 mm.
Additionally, to investigate the influence of θ c on the quality of the reconstruction, we also expanded the input mesh up to k = 40 assuming θ c = 5π/18.The RMSE for an expansion and reconstruction with θ c = 5π/18 was only 0.1215% larger than the one with θ c = π/18.Practically, this means that the analysis is invariant to θ c , with the small difference due to the numerical accuracy of the hypergeometric functions and the area distortion introduced by the parameterisation.The area distortion after finding the optimal Möbius transformation was 0.377619 for θ c = 5π/18 and 0.341917 for θ c = π/18.

Wavelength and patch size
As the face sculpture example contains large and small wavelengths, high orders of SCH were needed to reconstruct the small details; compare Figure 9(G) and (H), where 20 additional orders were added to capture the hair curls, the mouth and the nasal openings.To   The inset shows the Hausdorff distance at k = 15 where we used 602, 087 points sampled on the original mesh and compared these with the reconstructed mesh.At θ c = π/18, the RMSE for the Hausdorff distance was 0.149861, the mean error was 0.087238 and the highest error value scored was 1.594269.The reported data was extracted with the open source package MeshLab [66].
demonstrate that the number of orders needed for reconstructing a certain detail depends on the size of the patch, we considered a regional patch from the hair details to filter out the larger wavelength (facial features, e.g. the nose).The considered local patch is the red-dotted region circled on Figure 10(A).
We expanded the local patch by assuming θ c = π/18 and using only k = 15.As we see from Figure 11, 15 degrees were enough to reach almost the same RMSE accuracy as the whole face with k = 40.This resembles about 15.23% of the overall computed SCH basis functions, though this ratio does not scale linearly with the computational time.The minimum RMSE for the simple distance between the surfaces at k = 15 was 0.104812.By changing the range of targeted wavelengths, we can clearly see the faster convergence of this method with relatively close scales of detail (wavelengths).At only k = 15, the largest captured feature varied at a wavelength of 59.761 mm and the smallest varied at 4.269 mm, in comparison with k = 40, with 4.858 mm for the full sculpture.Furthermore, we compared this with the HSHA, and we show that the error at least doubles for the different measures; see Section 13.3 for full details.

Analysis of rough fractal surfaces
Here, we analysed rough artificially generated fractal surfaces, computed the Hurst coefficient from the SCHA for these surfaces and compared the obtained coefficients to the input values used for the generation.Many libraries have been developed for generating such surfaces using the traditional power spectral density (PSD) method (reviewed in [67,68]), such as Tamaas Library [69].These benchmarks are important in contact mechanics, which usually depends on exact asperity shapes, distributions and their effects on the real contact area, as per Hyun et al. (2004) [70].
We generated four fractal surfaces with a root-mean-squared roughness of 0.7E−2 mm and Hurst exponents of H ∈ {0.4,0.5, 0.6, 0.9}; see Figure 12.To avoid biasing the first few k's on fitting the shape of the boundary lines and to expand the surfaces with fewer ks, notice that the generated samples are all disks (circular boundary lines).To ensure reproducibility, all the surfaces were generated with a fixed seed number of 10 and with no roll-off wave vector.In the same figure, we show the analysis results (up to k = 40) in terms of the shape descriptors (solid lines) in Eq. ( 44) and the fit of Eq. (45) (dashed lines).
The surface in Figure 12(A) was generated with a Hurst exponent of 0.4, and we estimated a Hurst exponent of H = 0.4175, which corresponds to a fractal dimension of 2.5825 and an error of +4.4%.For the generated surface with H = 0.5 (Figure 12(B)), we estimated H = 0.5197 (+3.9%) and F D = 2.4802.For the third example with H = 0.6 (Figure 12(C)), we estimated the H = 0.5922 (−1.3%) and F D = 2.4078.However, for the fourth sample generated with H = 0.9 (Figure 12(D)), we obtained H = 0.8370 (−7.0%) and F D = 2.1630.This error is typical for higher values of H, as the surface seems finer and most of the spatial details are part of smaller wavelengths, meaning that additional ks are needed to capture these finer details.Making use of the scalability of the SCHA together with the nature of the fractal surfaces (repeats itself at all scales), we cut out a smaller area at the centre of the patch in Figure 12(D).We analysed the latter with k = 8, and the estimated Hurst exponent was 0.8850 and F D = 2.1150.This agrees with the conclusions drawn by Florindo and Bruno (2011) [71], who found using the PSD that the attenuation slope of F D changes with the value of k and that most of the fractal information is contained in the microscopic wavelengths when the change in the slope is negligible.

Laser-scanned rough surface patch
Here, we scanned a real stone that had been previously used in an experimental campaign [72].This type of stone is common in historical stone masonry walls, so studying the stones' interfacial properties helps define the mechanical behaviour of the walls.The stone surface was acquired by laser scanning (Figure 13) with an accuracy of 0.01 mm (Figure 17).Additionally, we sampled (extracted) a small patch over the scanned surface of the stone, see Figure 17, to study the roughness of the surface.
The surface patch was expanded up to k = 40.The estimated dimensions of the FDEC (computed from Eq. ( 43)) were: a = 24.7879mm, b = 23.3502mm and c = 5.7899 mm.From FDEC, we can estimate the wavelengths of any index k from Eq. (50).For instance, with an estimated average FDEC circumference of 151.2980 mm, the maximum and minimum computed wavelengths ω 2 ≈ 75.6490 mm and ω 40 ≈ 1.9908 mm, respectively.Figure 14 summarises the results of the analysis of the descriptors.From the analysis, the decay of the shape descriptors gives us a Hurst exponent of H = 0.9464, which suggests that the surface is not very rough.Later, this roughness information will be used for modifying the microstructure of a selected object.

Example of roughness projection
As explained in Section 7, we start the roughness projection by parameterising the patch on the donor mesh over a sphere with a prescribed θ c . Figure 15 shows the parameterisation of a selected patch on a donor mesh with various half-angles θ c .For the projection, we used the coefficients extracted from the patch in Figure 14.To alter the microstructure of the selected face on the donor mesh, we only include wavelengths between w max = 15.12 mm (k min = 10) and w min = 1.99 mm (k max = 40).The resulting microstructure (see Eq. ( 51)) is shown in Figure 16(A).To obtain the final finite element mesh, we remeshed the STL surface using the OpenCASCADE library [73] with various refinements (see Figure 16(B) and (C)).The final volumetric finite element mesh was obtained using GMESH [74].

Conclusions
In this paper, we conducted morphological analyses of open and nominally flat rough surfaces using the spherical cap harmonics (SCH), whose basis functions are a generalisation of the hemispherical harmonics (HSH) basis.SCHs form an orthogonal basis and is suitable for analysing an open surface that is projected onto a spherical cap with any half-angle θ c (Figure 1), while HSHs form a basis for a hemisphere (θ c = π/2).SCHs have been widely used in geophysics for describing fields (e.g. the magnetic field over a continent), but to our knowledge, they have not yet been applied for characterising and reconstructing the morphology of surfaces.We exploited the solution for the even set of Legendre functions of the first kind together with the Fourier functions to constitute the final SCHs that satisfy the Laplace equation on a spherical cap.We chose the even set to perform surface reconstruction with free edges using the Neumann boundary conditions while simultaneously conducting a power spectral analysis using the orthogonality of the even basis functions.We also used SCH to analyse simply connected surfaces with vertices, edges and triangulated faces (the Standard Triangle Language STL files).The first step of such analyses is to couple each vertex on the surface with a unique pair of θ and φ (surface parameterisation).Thus, we proposed a parameterisation algorithm that provides conformal one-to-one mapping with minimal area distortion over a unit spherical cap with a prescribed half-angle θ c (S 2 θ≤θc ).We showed that the proposed analysis approach is invariant to θ c yet is sensitive to the distortion introduced in the parameterised surface at different half angles (θ c 's).The unavoidable distortion caused by the parameterisation algorithm can limit the analysis to certain half angles.Therefore, we used a global optimisation approach to identify the optimal half-angle θ c that minimises the area distortion between the input and parameterised surfaces.As expected, we found that the conventional HSH analysis is suitable for surfaces that can be parameterised over a unit hemisphere without introducing significant area and angular distortions.For nominally flat surfaces, however, the traditional HSH fails to analyse and reconstruct the surfaces correctly because of the distortion.
To show the robustness of the morphological analysis and reconstruction convergence, we applied the proposed method to a scanned face sculpture as a visual benchmark.This benchmark was chosen to contain relatively large and small wavelengths to test the convergence to the smallest details.The face was analysed for two selected half angles to demonstrate that the method is invariant to θ c and that the quality of the analysis depends only on the distortion introduced on the parameterised surface.We then demonstrated the convergence and the scalability of the method to expand and reconstruct complex surfaces of a selected local patch on the face.Because the patch was smaller than the entire face, fewer degrees were needed to capture the same details (faster convergence).We finally compared this method with the HSH and found that the reconstruction from HSH was wavy and showed larger error margins, as was expected due to the distortion in the parameterisation.
We derived an approximate formula that links the size of the first degree ellipsoidal cap (FDEC) with the wavelengths as a function of the index k.We then used the power spectral analysis of the rotation-invariant shape descriptors that are invariant to translation, rotation, scale and half-angle θ c to find the fractal dimension (FD) of the analysed open surfaces.To benchmark this approach, we generated fractal surfaces with the traditional PSD method and computed the Hurst exponent for these surfaces from the SCHA results.The method yielded good estimates of the Hurst exponent (errors between −7.0 and +4.4%).
As this method was proposed to study rough surfaces in contact, we also outline a method for projecting real or artificially generated roughness onto a selected donor mesh.This method can be used for the traditional SH, HSH and the SCH.We limited the bandwidth in the spectral domain to reconstruct and alter the roughness of surfaces between user-defined upper and lower limits.
The computational cost of this method was found especially expensive when evaluating the sequential hypergeometric function, the critical loop, for k > 12.This is mostly due to using a high-level programming language such as MATLAB and can be bested with using C or C++.

Future works
In this section, we outline future developments for SCHA: • This work hinges on the numerical stability and accuracy of the computed basis functions of SCH, which depends on the Gaussian hypergeometric function evaluation • The herein proposed method can be seen as a corner stone for regional wavelet analyses for shape morphology.Wavelet analyses conducted via traditional methods like SH could take advantage of the stable and well-studied evaluation techniques of the traditional basis function as well as any traditional recursive formulae (e.g.SH basis).
The SCH serves as an exact solution to compare with other potential regional analysis methods using traditional basis functions.
• Implementing a recursive approach for computing the SCH basis can improve the speed and efficiency.
• A physical regularisation method based on optimising the estimated SCH coefficients of the basis will improve the accuracy and noise levels of the analysis (see [12] as an example for SHA).
12 Reproducibility The codes and generated data and the surfaces are all made available in this paper and can be found on Zenodo (GitHub) [75]: • Zenodo: 10.5281/zenodo.4890809 • GitHub: git@github.com:eesd-epfl/spherical-cap-harmonicsThe reconstruction of the stone made using SH to analyse the whole stone (L max = 45), HSH for a regional analysis (L max = 45) and SCH for the regional analysis of nominally-flat patches (K max = 40).Figure 19: The optimal half-angle θ c for the spherical cap parameterisation obtained for the area extracted from the scanned stone in Figure 17; the map shows the radial distance from the centroid of the extracted part.

Additional results for the Sturm-Liouville eigenvalues
In this section, we show the Neumann boundary conditions and eigenvalues of the S-L problem for visualising the sequence of eigenvalues and their effects on the wavelengths.Figure 21 shows the boundary condition of the even set.The locations of the eigenvalues on this figure occur where we see the asymptotic points all over the surface that approach −∞ when the dP m l(m) k (xc) dx → 0 (eigenvalues), localising the points (asymptotics).These asymptotic lines are more obvious when the equations are plotted with greater accuracy (smaller step-size).Figure 22 shows the boundary conditions for three selected orders and explains how the size of the plateaus increase with m and thus how they affect the wavelength contributions at different levels when they come into effect in the expansion series.Figure

Figure 1 :
Figure 1: Representation of a polar spherical cap S 2θ≤θc with half-angle θ c , where the spherical cap harmonics analysis (SCHA) is defined over.

Figure 3 :
Figure 3: The normalised spherical cap harmonic (SCH) basis functions up to k = 4 at θ c = 5π/18.The figures show R

Figure 4 :
Figure 4: Parameterisation of the monkey head model.A) The monkey head (Suzanne) benchmark mesh [53], with the heat map (texture) illustrating the radial distance from the model's centroid; B) The back view of the model reveals the open boundary location; C)The conformal mapping onto the unit disk obtained using[5]; D) The convergence curve for finding the optimal Möbius transformation that minimises the area distortion on the spherical cap using the PSS algorithm f : S → S 2 θ≤θc for θ c = 40 • (i.e.r = 0.3640).

Figure 6 :
Figure6: Transformation of a square grid to a geodesic dome grid (before triangulating the surfaces) using the algorithm proposed by Roşca (2010)[57].A) A square grid with 60 horizontal and vertical lines resulting in 900 points of intersection (vertices); B) An areapreserving map of the square grid onto a quad-grid disk; C) The final geodesic dome with θ c = π/2 (a hemisphere); D) The final geodesic dome with θ c = π/3; E) The final geodesic dome with θ c = π/18.

FDECFigure 8 :
Figure 8: The first-degree ellipsoidal cap (FDEC) defined by k = 1; a is the half-major axis, b is the half-minor axis and c is the ellipsoidal cap depth, where |c| ≤ |b| ≤ |a|.

Figure 9 :
Figure 9: The reconstruction of the face sculpture with SCH.A) The input surface after remeshing with 21, 280 vertices.It was expanded up to k = 40 and reconstructed with an edge length of 0.032 on a unit disk; B) The reconstruction with k = 1 results in an FDEC with 2a = 123.381,2b = 111.336and c = 19.844;C) Reconstruction at k = 3; D) Reconstruction at k = 5; E) Reconstruction at k = 10; F) Reconstruction at k = 15; G) Reconstruction at k = 20; H) Reconstruction at k = 40.

Figure 10 :
Figure 10: The root mean square error (RMSE) for the face sculpture.The inset shows the Hausdorff distance at k = 40, where 1, 221, 276 points were sampled on the original mesh and were compared with the reconstructed mesh.For θ c = π/18, the RMSE for the Hausdorff distance was 0.229886, the mean error was 0.141412 and the highest error value scored was 2.799484.The reported data was extracted with the open source package MeshLab [66].

Figure 11 :
Figure 11: The root mean square error (RMSE) for the local patch on the visual benchmark.The inset shows the Hausdorff distance at k = 15 where we used 602, 087 points sampled on the original mesh and compared these with the reconstructed mesh.At θ c = π/18, the RMSE for the Hausdorff distance was 0.149861, the mean error was 0.087238 and the highest error value scored was 1.594269.The reported data was extracted with the open source package MeshLab[66].

FreqencyFigure 12 :
Figure 12: The analysis results for artificially generated fractal surfaces with different Hurst exponents.The solid lines represent the descriptors in Eq. (45) and the dashed lines represent the fit of Eq. (45).(A) a surface generated with H = 0.4; (B) a surface generated with H = 0.5; (C) surface generated with H = 0.6; (D) surface generated with H = 0.9.

Figure 13 :
Figure 13: The laser scanner setup using the Quantum FARO Arm ® ; the left side of the picture shows the laser scanner used to scan the stone, and the right side shows the stone pinned on a steel frame for the scan to minimise the loss of surface data.

Figure 14 :Figure 15 :Figure 16 :
Figure 14: The SCHA for a patch sampled over the stone in Figure 17.The analysis of the results for k = 40 and the best fit for the first k = 12 are shown, with H = 0.9464 and F D = 2.0537.The colour bar shows the height map of the patch in mm.The upper x-axis shows the corresponding wavelengths ω k .

Figure 17 :
Figure 17: The laser scanning data after down-sampling the point cloud to 444, 653 points instead of nearly 23, 000, 000 points; the colour map shows the radial distance measured from the centroid of the stone.

Figure 18 :
Figure 18:  The reconstruction of the stone made using SH to analyse the whole stone (L max = 45), HSH for a regional analysis (L max = 45) and SCH for the regional analysis of nominally-flat patches (K max = 40).

Figure 20 :
Figure20: The root mean square error (RMSE) for the local patch on the visual benchmark, as computer using HSH.The inset shows the Hausdorff distance at k = 15 where we used 1, 202, 087 points sampled on the original mesh and compared this value with the reconstructed mesh.At θ c = π/18, the RMSE for the Hausdorff distance was 0.310431, and the mean error was 0.087238, while the highest error value scored was 1.594269.The reported data was extracted using the open source package MeshLab[66].

Table 1 :
The results summary of 51 consecutive runs of the PSS algorithm for the benchmarking test shown in Figure