© The Author(s), published by EDP Sciences, 2026

1 Introduction

The progress in high-resolution and efficient optical metrology technologies –such as optical coherence tomography or intensity-based wavefront sensing [1–3] – applied to measuring either ocular interfaces or, indirectly, the wavefront surface has permeated ophthalmology practice. The highly dense sampling of these surfaces consequently leads, at least for some applications, to local surface representations such as splines [4, 5] or mesh surfaces [6], i.e., those constructed by a set of polygonal faces, typically triangles. Furthermore, these high-spatial-sampled surfaces encourage the application of spectral analysis to them. For instance, classical Fourier techniques have been proposed for analysing keratoconus eyes – i.e. an eye condition characterized by a thinning and bulging of the cornea degrading vision – through wavefront phase imaging [3]. Besides Fourier analysis, Zernike decomposition of optical surfaces is undoubtedly the most typical approach to analyze, or characterize, optical wavefront or surfaces. However, both Zernike and Fourier coefficients do not fully capture the shape information, especially the effect of curvature on the spectrum itself [7]. Indeed, this intrinsic limitation of conventional Zernike polynomials led to the proposal of some alternatives, such as the so-called Zernike curvature polynomials [8]. However, these are only valid when the surface gradient is small, so the wavefront curvature reduces to its Laplacian; an approximation that is not sustained in many practical cases.

By analogy to classical Fourier analysis, spectral methods over surfaces defined as meshes were introduced at the end of the last century for discrete geometry processing [7]. The goal of these methods is to capture the geometrical shape properties of the surface, which is referred to as shape analysis. Mesh surfaces are especially suitable for shape analysis thanks to the fact that discrete differential geometry techniques can be directly applied to them [9]. A pioneering work in spectral methods over meshes was that of Taubin [10], which defined the mesh vertex coordinates as a 3D signal and applied Laplacian operators. Another option is to define a scalar function over the mesh surface.

Surface spectral methods comprise three steps: 1) selecting an appropriate linear operator; 2) the eigendecomposition of that operator; and 3) employing the spectral eigenvalues and eigenvectors in a specific manner according to the processing purpose. In other words, spectral methods inspect a surface by examining the eigenvalue decomposition of a meaningful selected linear operator. Therefore, the application at hand determines the selection of the linear operator and how to employ the spectrum.

Surface spectral methods have been used intensively in a diversity of techniques [11], such as mesh compression [12, 13], segmentation [14], smoothing or denoising through filtering [15], symmetry detection [16, 17], surface reconstruction from point clouds [18], or watermarking [19], all of them relevant in image processing applications.

As mentioned before, the first and most crucial step in spectral methods is the selection of an appropriate linear operator. In optics, the curvature – in its various mathematical forms – is the most relevant geometrical property because it linearly depends on the optical power, the fundamental property in all imaging systems. In any form, the curvature involves second-order derivatives, so a mandatory choice for such an operator is that it should include them. Then, the Laplace–Beltrami operator emerges as the best option because it involves second-order derivatives, but it is linear within the surface geometry coordinates. Here, I use the so-called geometric mesh Laplacian, because it explicitly encodes the geometric information of the mesh [7]. Furthermore, selecting a curvature-based local shape index as the scalar function defined over the mesh leads to an optically meaningful Laplace–Beltrami spectrum [7, 20]. One of the potentials of the Laplace–Beltrami, for our purposes, stems from the fact that it characterizes the smoothness of a scalar function over a surface, grosso modo, the difference between the scalar at a point P and a weighted average in a small neighborhood of P (Sect. 6.2 [7]).

For certain applications, shape analysis focuses on capturing morphologically relevant shape descriptors, also called signatures or fingerprints. Such a signature has been called shape-DNA; so, from now on, I will follow this notation [11]. As it will be seen, employing a local curvature shape index scalar (for instance, the Gaussian curvature) notably increases the discrimination capacities of the Laplace–Beltrami spectrum.

To date, I am not aware that geometric mesh analysis has been applied in the visual optics literature, except for a paper in retina imaging [21]. I believe that there is a call to introduce these techniques, especially for analysing anterior chamber depth interfaces, ocular wavefront surfaces, or even ophthalmic lens surfaces. In this paper, I present a systematic approach to these techniques and provide some applications. A major hypothesis is that surface second-order differential features, captured by the Laplace–Beltrami spectrum, can be employed to detect, compare, and study ocular surface shape features. The paper is organized as follows. First, in Section 2, I introduce all the basic mathematical aspects of the Laplace–Beltrami spectra decomposition theory. Next, in Section 3, I introduce a particular, but meaningful, geometric eye surface-thickness model that will be used in subsequent sections, particularly a keratoconus model. Subsequently, in Sections 4 and 5, I present two relevant applications of the theory: the spectral decomposition and the shape-DNA analysis, respectively. Finally, in Section 6, I summarize results and discuss their implications.

2 Laplace–Beltrami spectral decomposition theory

The Laplacian operator (Δ), i.e., the divergence of the gradient $\Delta =\frac{{\partial }^2}{{\partial }^{2}x}+\frac{{\partial }^2}{{\partial }^{2}y}$, measures how a function differs from its average in a local neighborhood. The Laplace–Beltrami operator is a generalization of the Laplacian to functions defined over surfaces or manifolds in higher dimensions.

Let f ∈ C ² be a real-valued function (continuous second derivatives) defined over a smooth bounded surface S, locally represented by a parametrization: r(x, y) ≡ (x, y, z(x, y)): S ² → R ³, being S ² a two-dimensional (circular or not) Cartesian parametric domain. Denoting ${\mathbf{r}}_x=\left(\mathrm{1,0},\frac{\partial z}{\partial x}\right)\equiv \left(\mathrm{1,0},z_x\right)$ and ${\mathbf{r}}_x=\left(\mathrm{0,1},\frac{\partial z}{\partial y}\right)\equiv \left(\mathrm{0,1},z_y\right)$, the Riemannian metric tensor or first fundamental form matrix (symmetric) is:$$\left[\begin{array}{cc}{g}_{11}& g_{12}\\ g_{21}& g_{22}\end{array}\right]\equiv \left[\begin{array}{cc}{\mathbf{r}}_x{\mathbf{r}}_x& {\mathbf{r}}_x{\mathbf{r}}_y\\ {\mathbf{r}}_y{\mathbf{r}}_x& {\mathbf{r}}_y{\mathbf{r}}_y\end{array}\right]=\left[\begin{array}{cc}(1+z_x{)}^2& z_x{z}_y\\ z_x{z}_y& (1+z_y{)}^2\end{array}\right],$$

which determinant is $g=g_{11}{g}_{22}-g_{12}^2$.

Now, f_x, f_y, f_xx, f_yy, f_xy, denoting partial derivatives of f and considering relations between the elements of the Riemannian metric and those of the inverse matrix (p. 153 [22]), the Laplace–Beltrami operator over f is defined as [23]:$${\Delta }_g\left[f\right]=\frac{g_{12}\left(f_{\mathrm{yx}}-{\mathrm{\Gamma }}_{21}^1{f}_x-{\mathrm{\Gamma }}_{21}^2{f}_y\right)-g_{11}\left(f_{\mathrm{yy}}-{\mathrm{\Gamma }}_{22}^1{f}_x-{\mathrm{\Gamma }}_{22}^2{f}_y\right)}{g}-\frac{g_{22}\left(f_{\mathrm{xx}}-{\mathrm{\Gamma }}_{11}^1{f}_x-{\mathrm{\Gamma }}_{11}^2{f}_y\right)-g_{12}\left(f_{\mathrm{xy}}-{\mathrm{\Gamma }}_{12}^1{f}_x-{\mathrm{\Gamma }}_{12}^2{f}_y\right)}{g}.$$(1)

Here, ${\mathrm{\Gamma }}_{\mathrm{ij}1}^k$ are the so-called Christoffel symbols; their values as a function of partial derivatives of z have been derived elsewhere [24]:$${\mathrm{\Gamma }}_{11}^1=\frac{z_{\mathrm{xx}}{z}_x}{(1+z_x^2+z_y^2)}, \hspace{1em}{\mathrm{\Gamma }}_{12}^1=\frac{z_{\mathrm{xy}}{z}_x}{(1+z_x^2+z_y^2)},$$ $${\mathrm{\Gamma }}_{22}^1=\frac{z_{\mathrm{yy}}{z}_x}{(1+z_x^2+z_y^2)},\hspace{1em} {\mathrm{\Gamma }}_{11}^2=\frac{z_{\mathrm{xx}}{z}_y}{(1+z_x^2+h_y^2)},$$ $${\mathrm{\Gamma }}_{12}^2=\frac{z_{\mathrm{xy}}{z}_y}{(1+z_x^2+z_y^2)},\hspace{1em} {\mathrm{\Gamma }}_{22}^2=\frac{z_{\mathrm{yy}}{z}_y}{(1+z_x^2+z_y^2)}.$$

I note that the Laplace–Beltrami operator encodes the intrinsic geometry, i.e., that part of it not depending on the surface parameterization, so any other parametrization (e.g., polar parametrization) instead of the Cartesian one would provide the same geometrical structure. Also, it is worth mentioning that when the local curvature of the surface is quasi-flat (z_x, z_y ≪ 1), the Laplace–Beltrami operator is approximately equal to the Laplacian.

In discretized surfaces, the discrete Laplace operator replaces the continuous one. The surface is described by a triangular mesh M, described by a set of geometric vertices coordinates, a matrix V = (x₁, …, x_n, y₁, …, y_n, z₁, ..., z_n), and an ordering collection of those vertices through edges. Then, the discrete Laplace–Beltrami operator (Δ_M) of a scalar function f at a vertex i can be estimated with the cotangent formula [7]:$${\Delta }_M[f{]}_i=\frac{1}{2A_i}\sum _j\mathrm{}(\cot {\alpha }_{\mathrm{ij}}+\cot {\beta }_{\mathrm{ij}}\left)\right(f_j-f_i),$$(2)where the sum is over all adjacent vertices j of i, α_ij and β_ij are the two angles subtended by the edge joining i and j, and A_i is the vertex area of i.

The computation of equation (2) is facilitated if the operator is encoded in matrix form that is obtained by the product of a diagonal matrix and a symmetric matrix (Sect. 6.2 [7]):$$\mathcal{L}={\mathcal{B}}^{-1}\mathcal{C},$$(3)where $\mathcal{B}$ is the diagonal matrix whose entries are 1/2A_i and $\mathcal{C}$ is a symmetric matrix – the so-called cotangent – encoding the sum of the cotangent angles. The Laplacian matrix is largely sparse, with an average of seven non-zero values per row [7]. However, the matrix $\mathcal{L}$ is, in general, not symmetric, which makes the eigenvalue problem more demanding. Luckily, $\mathcal{C}$ is symmetric, so a good approach (followed here) is to obtain the spectral decomposition by solving the generalized eigenvalue problem (Sect. 9.1 [7]):$$\mathcal{C}f=\lambda \mathcal{B}f,$$(4)which provides the same spectral decomposition as finding the eigenvalues of $\mathcal{L}$ directly.

I solved the generalized eigenvalue problem of equation (4) using Matlab(R) built-in function eigs.m.

The eigenvectors obtained are sorted in ascending eigenvalues: $\mathcal{T}=[{\varphi }_1,...,{\varphi }_k]$, and excluding the smallest eigenvalue ϕ₀ that is zero, so such a truncated decomposition can be interpreted as a low-pass filter because the smallest eigenvalues correspond to the less ‘curved’ eigenvectors.

Then, the Laplace–Beltrami spectral reconstruction of the surface mesh is obtained by applying the matrix transformation [25]:$$\mathrm{V}=\mathcal{T}\left({\mathcal{T}}^T\mathcal{V}\right),$$(5)where, recall, that $\mathcal{V}$ is the mesh vertex matrix.

2.1 Laplace–Beltrami spectra of scalar functions defined over the mesh

Given the scalar function f defined over the mesh M, the set of Laplace–Beltrami eigenfunctions makes it possible to perform a spectral decomposition of f on S, provided that the eigenfunctions are normalized such that (ϕ_i , ϕ_j) = δ_i,j. We denote this set of normalized eigenvectors: $T_n=[{\varphi }_{n1},...,{\varphi }_{\mathrm{nk}}]$. Then, the Laplace–Beltrami spectrum is obtained by a projective matrix operation [7]:$$\widehat{f}=T_n^{T}f.$$(6)

As will be presented later, I have evaluated the Laplace–Beltrami spectra of two scalar functions: thickness and Gaussian curvature, and for a specific simulation the mean curvature.

Evaluating the Gaussian curvature (or any other differential quantity on the mesh) involves a trade-off in the neighbourhood size used for numerical estimations. Although for data free of noise, smaller sizes increase accuracy, in the presence of noise, larger neighbourhoods may be more robust to the effect of noise. In Section 3, this trade-off will be illustrated.

For Gaussian curvature computation, the first step is to estimate a normal vector n_i associated with each vertex of the mesh surface. This was done by taking the weighted average of the normals associated with each face touching the vertex [26]. The weights were chosen by following the procedure proposed at [27]. Then, selecting a direction determined by two, p_i and p_j, neighbourhood vertices points, the normal curvature at p_i in that direction is given by [26]:$$k_{\mathrm{ij}}^n=\frac{2{\mathbf{n}}_i({\mathbf{p}}_i-{\mathbf{p}}_j)}{|{\mathbf{p}}_i-{\mathbf{p}}_j{|}^2}.$$(7)

Now, similar to how the first fundamental form is defined, the second fundamental form, denoting n the surface normal, is:$$\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]\equiv \left[\begin{array}{cc}{\mathbf{r}}_{\mathrm{xx}}\mathbf{n}& {\mathbf{r}}_{\mathrm{xy}}\mathbf{n}\\ {\mathbf{r}}_{\mathrm{yx}}\mathbf{n}& {\mathbf{r}}_{\mathrm{yy}}\mathbf{n}\end{array}\right]=\left[\begin{array}{cc}\frac{z_{\mathrm{xx}}}{\sqrt{1+z_x^2+z_y^2}}& \frac{z_{\mathrm{xy}}}{\sqrt{1+z_x^2+z_y^2}}\\ \frac{z_{\mathrm{yx}}}{\sqrt{1+z_x^2+z_y^2}}& \frac{z_{\mathrm{yy}}}{\sqrt{1+z_x^2+z_y^2}}\end{array}\right].$$

From the set of normal and normal curvatures, an estimation of the second fundamental matrix above is obtained using least squares [28]. Finally, I recall that the principal curvatures, k₁ and k₂, are the eigenvalues of the second fundamental form, the Gaussian curvature is its product: $\mathcal{K}=k_1{k}_2$, and the mean curvature its average $\mathcal{H}=\frac{k_1+k_2}{2}$. In practice, for curvature computations, I used Matlab code developed by Itzik Ben, freely available at: https://es.mathworks.com/matlabcentral/fileexchange/47134-curvature-estimationl-on-triangle-mesh.

3 Example of application: artificial eye surface model

On the one hand, anterior segment ocular surface interfaces (crystalline lenses and cornea) can be modelled to some extent as a biconic function [29] plus some perturbation terms. On the other hand, a certain type of surface perturbation can be generically described with a Gaussian function; for instance, a keratoconus [30]. Considering the above, as an example of application, I modelled an anterior cornea surface topography at each (x, y) location with the following equation:$$z(x,y)=\mathrm{piston}+\frac{\frac{x^2}{R_x}+\frac{y^2}{R_y}}{1+\sqrt{1-\frac{(1+Q_x)x^2}{R_x^2}-\frac{(1+Q_y)y^2}{R_y^2}}}+h_0{e}^{-\frac{(x-x_0{)}^2}{2{\sigma }_x^2}-\frac{(y-y_0{)}^2}{2{\sigma }_y^2}},$$(8)and the thickness at each point with:$$z(x,y)=0.5 \mathrm{mm}-h_0{e}^{-\frac{(x-x_0{)}^2}{2{\sigma }_x^2}-\frac{(y-y_0{)}^2}{2{\sigma }_y^2}}.$$(9)

The piston term is just an arbitrary number – particularly, 10 – selected to avoid negative numbers in z, but, of course, without changing the shape. Following equation (8), I computed the Gaussian curvature following the procedure described in Section. 2.1. These functions are defined within a circular domain $\sqrt{x^2+y^2}<5\mathrm{ }\mathrm{mm}$. Considering that real biometric data is affected by measurement noise, I added white Gaussian noise to the z values with the Matlab function awgn. The level of Gaussian noise was simulated through the Signal-to-Noise Ratio (SNR), whose units are dBs.

Figure 1 plots the surface mesh, the Gaussian curvature, and the thickness corresponding to these equations with the following parameters: R_x = 7.79 mm, R_y = 7.61 mm, Q_x = −0.27, Q_y = −0.27, h₀ = 0.2 mm, σ_x = 0.5, σ_y = 0.5, x₀ = 1.5 mm, and y₀ = −1.5 mm and with a SNR = 90 dB. From now on, this cornea model will be referred to as test A.

Figure 1 Mesh, Gaussian curvature, and thickness (mm) over the mesh of reference cornea A.

As mentioned in Section 2.1, there is a trade-off between mesh density and noise level. To exemplify this fact, I computed the Gaussian curvature map for two noise levels: low noise (SNR = 90 dB) and moderate noise (SNR = 70 dB), and three levels of mesh sampling density determined by the number of mesh vertices: 5.025, 20.081, and 45.225 points. Figure 2 shows the results. Although for low noise, the more accurate Gaussian curvature estimation is achieved with high sampling density (45.225 mesh points), for moderate noise, the best estimation requires fewer mesh sampling points (between 20.081 and 45.225). Therefore, the choice a priori of the optimal mesh density, to improve the SNR ratio, is not an easy task because it also depends on the geometry shape. Nevertheless, a rule of thumb is that the stronger the local curvature changes, the denser the mesh required.

Figure 2 Gaussian curvature over mesh of surface test A for different mesh densities and two levels of noise: SNR = 70 dB and SNR = 90 dB.

4 Spectral decomposition

A first potential of the Laplace–Beltrami operator comes from its spectral decomposition obtained with equation (5). As mentioned in the introduction, it can be used for mesh compression or smoothing through filtering. In both cases, the Laplace–Beltrami operator serves to separate the curvature shape features of the surface in question. The first eigenvectors of the spectral decomposition (associated with low eigenvalues) capture the low-curvature shape signal, whereas the following ones apprehend the high-curvature part of the shape.

Figure 3 shows the mesh reconstruction of surface test A with an SNR = 80 dB for different numbers of eigenvectors in the spectral decomposition, particularly: 5, 10, 25, 50, 75, 100, 200, and 300. The figure shows that only when more than approximately 50 eigenvectors are included does the “bump” part of the keratoconus start to be revealed. Note that, as with any spectral method, some artifacts are introduced at the periphery. This problem, related to the boundary conditions in the Laplace–Beltrami will be discussed in Section 6.

Figure 3 Original mesh corresponding to surface model A and spectral mesh reconstruction using different numbers of eigenvectors in the spectral decomposition.

Generally, the problem of choosing the optimal number of eigenvectors would depend on the shape structure and the level of measurement noise. Therefore, preliminary modelling of the shape features of interest would surely help in that selection.

5 Laplace–Beltrami for shape-DNA analysis of ocular surfaces

A substantial benefit of the Laplace–Beltrami spectra is their invariance to isometric changes. These are deformations that preserve the length of every arc on the surface, so they conserve the shape. The most common isometric changes are rigid motions (rotation or global translation) and similarity scaling [11]. On the one hand, rotation or global translation is a typical issue due to eye fixation changes among measurements, for instance, in optical coherence tomography [31]. On the other hand, different eyes may differ in size but not in ocular shape; in this case, the Laplace–Beltrami spectra would be simply scaled, as will be shown.

The isometric invariance of the Laplace–Beltrami spectrum enables its use to extract geometrical fingerprints characterizing ocular surfaces, which has been coined a shape-DNA [11]. This signature is derived from the Laplace–Beltrami spectra as computed through equation (6).

The shape-DNA emerges akin to Fourier analysis as a spectral tool analysis, but contrary to the latter, a tool that captures better the shape geometry. Also, contrary to Fourier analysis, the surface mesh geometry determines different eigenbases of the Laplace–Beltrami spectrum, so the corresponding eigenvalues of different meshes cannot be directly compared; some type of normalization is required so that the shape-DNA is invariant to the surface’s size [7, 11]. Particularly, I propose scaling the Laplace–Beltrami spectra coefficients to the first coefficient: $\widehat{f}=\frac{\widehat{f}}{\widehat{f}\left(1\right)}$, a procedure that is appropriate when trying to identify a surface shape pattern [11].

Furthermore, to apply the shape-DNA for quantitative comparison purposes, or to identify surface-mesh shape patterns, a distance metric must be constructed [11]. Given two arbitrary shape-DNAs (normalized set of eigenvalues): λ = (λ₁, λ₂, … λ_n) and μ = (μ₁, μ₂, … μ_n), from now on it will be used the 2-norm distance, which is unitless because the above normalization, expressed as:$$d(\lambda ,\mu )={\left(\sum _{i=1}^n\mathrm{}({\lambda }_i-{\mu }_i{)}^2\right)}^{1/2}.$$(10)

5.1 Shape-DNA similarity: continuous deformation

The shape-DNA offers a way to quantitatively track, with a single metric, the timing of ocular surface shape changes. This could be useful, for instance, as a clinical tool to analyze the time evolution of pathologies manifested by surface shape abnormalities. The underlying mathematical principle is that of similarity; namely, the shape-DNA fingerprint should depend continuously on shape deformations [11].

As an example, I propose an artificial model of keratoconus variation more or less compatible with the reported literature [32]. Starting with surface model A, I modelled a surface elevation change and thinning of the cornea with time manifested by a change of h_o ranging from 0.2 to 0.4; and, synchronous with this process, a spatial expansion of the local abnormality modelled by a change of σ_x and σ_y values from 0.5 to 2. Three different shape-DNA were generated: associated with the Gaussian curvature, the mean curvature, and with the cornea thickness.

I computed the shape-DNA distance (Eq. (10)) between all these values and the initial reference value (h_o = 0.2, σ_x = 1, and σ_y = 1), sampling the ranges in 100 steps. Results are plotted in Figure 4, revealing that the similarity principle is fulfilled, because the shape-DNA fingerprint depends smoothly on the continuous variation of the h_o,σ_x and σ_y. However, it should be observed that this dependence is not necessarily linear or even monotonically increasing, as revealed by the Gaussian and mean curvature curvature fitting curves. Also, the relation depends strongly on the scalar function defined over the mesh, as revealed by the difference in the fitting curves of the Gaussian curvature and thickness. To complement this figure, Figure 5 shows the mesh, Gaussian curvature, and thickness for three representative stages of the previous timing model.

Figure 4 Shape-DNA distance to reference cornea model A for different scalar functions – Gaussian curvature, mean curvature and thickness – as function of timing surface deformation: explicitly the ho parameter, but implicitly also σx and σx.

Figure 5 Surface meshes and Gaussian curvature as function of shape timing deformation: A) ho = 0.264, σx = 0.98, and σy = 0.98; B) ho = 0.332, σx = 1.49, and σy = 1.49; C) ho = 0.4, σx = 2, and σy = 2.

5.2 Shape-DNA similarity: stability

As mentioned before, the isometry invariance of the Laplace–Beltrami operator implies that the fingerprint should not differ between surfaces that can be isometrically mapped onto each other. The most common isometrical mappings are rigid motions and scalings. Additionally, as it will be shown in Figure 6, the Gaussian curvature has the same value at corresponding points of isometric surfaces (p. 177 [33]).

Figure 6 Surface, Gaussian curvature and thickness over meshes: A) Reference mesh; B) Reference mesh transformed by an eye fixation change. C) Reference mesh transformed through a uniform scaling.

Starting again with reference surface A, I generated two isometric surfaces to it, exemplifying two typical scenarios found in visual optics biometry. In the first case, I modelled a change in eye fixation between two measurements of the same ocular surface as a direct rigid transformation between points of the original mesh [x_i, y_i, z_i] onto the transformed one [x_f, y_f, z_f] based on the matrix operator given by the equation:$$\left[\begin{array}{l}{x}_i\\ y_i\\ z_i\\ 1\end{array}\right]=\left[\begin{array}{llll}{R}_{11}& R_{12}& R_{13}& t_x\\ R_{21}& R_{22}& R_{23}& t_y\\ R_{31}& R_{32}& R_{33}& t_z\\ 0& 0& 0& 1\\ & & & \end{array}\right]\times \left[\begin{array}{l}{x}_f\\ y_f\\ z_f\\ 1\end{array}\right],$$(11)where t_x, t_y and t_z are translations along directions x, y and z; and R_ij are entries of the Euler rotation matrix [34]. I implemented this with the help of the Matlab function rigidtform3d.m. Particularly, I modelled an eye fixation change provided by an Euler angle α = 45° and a translation defined by [t_x, t_y, t_z] = [0.4, 0, 0]. Applying this transformation to the reference surface A, surface B is obtained.

In a second example, I modelled a uniform scaling transformation that isotropically enlarges or shrinks the ocular surface by a scale factor equally in all directions. This is implemented through a scaling matrix operator.$$\left[\begin{array}{l}{x}_i\\ y_i\\ z_i\end{array}\right]=\left[\begin{array}{lll}\upsilon & 0& 0\\ 0& \upsilon & 0\\ 0& 0& \upsilon \\ & & \end{array}\right]\times \left[\begin{array}{l}{x}_f\\ y_f\\ z_f\end{array}\right],$$(12)where υ is the scalar factor. Using a scaling dilatation υ = 1.005 onto the surface A, generate surface C.

The meshes and Gaussian curvatures over them of surfaces A, B, and C are depicted in Figure 6. It illustrates the change in eye fixation between surface B and A, as well as the uniform dilatation of mesh C with respect to A. Observing the Gaussian curvature map, its invariance with scaling shows the potential of this map with respect to plain meshes.

The shape-DNA distances between reference surface A, and surfaces B, C, are 8E–13, and 3E–13, respectively; quantities which are clearly negligible.

6 Conclusion

The large increase, in recent times, of the biometry capabilities applied to the anterior segment of the eye has led to a synchronous development of data processing techniques, which include compact and statistical surface representation [35, 36], image segmentation [37], or data interpolation/extrapolation [38, 39]. A similar process has occurred in reconstructing, modelling, or analyzing wavefronts [3, 40].

However, all these techniques, although indirectly dealing with curved surfaces, do not fully address the curvature shape information. In this sense, the Laplace–Beltrami spectral decomposition may emerge as a new tool offering this advantage. Of course, further exploration of the potential of the methodologies proposed in this paper requires testing them with real clinical data, such as optical coherence tomography measurements or highly dense wavefront sensor data.

As shown, the spectra are especially suitable for performing shape analysis through shape-DNA fingerprints. However, these fingerprints can be spatially localized, occupying different scales. When the objective is to perform a global shape comparison, but at the same time compare local geometrical features, a possible approach is to pick up the spectra associated with selected eigenfunctions [20]; for instance, in the discriminative representation of brain morphology [41]. This is something that could be explored in the future.

In the mathematical part, Laplace–Beltrami spectral theory still offers challenging open questions. One of which is how to address boundary conditions in open bounded surfaces. Although in the present approach, I have not complied with boundary conditions, because the discrete Laplace–Beltrami operator based on the cotangent formula can be implemented without them. However, boundary conditions affect both the eigenvalues and their corresponding eigenfunctions [11]. Classical boundary conditions, namely homogeneous Dirichlet – constant values of the solution along the boundary – or Neumann – zero derivative of the solution with respect to the outward normal –, are not strictly valid on the human eye anatomical or wavefront surfaces. Therefore, I hypothesize that some improvements in the solutions – related to edge reconstruction errors as shown in Figure 3 – are expected if proper numerical considerations for the boundaries are included. I plan to research rigorous options related to more realistic boundary conditions – the more reasonable are non-homogeneous Dirichlet conditions –, as explored in the mathematical literature [42].

Finally, it is worth mentioning that machine learning techniques –which are also increasingly being used in visual optics biometry processing [37, 43]– are boosting spectral shape analysis in what has been named geometric deep learning [44].

Funding

This work was supported by grant PID2023-150166NB-I00 funded by MCIN/AEI/10.13039/501100011033.

Conflicts of interest

The author declares that he has no competing interests to report.

Data availability statement

The research data associated with this article are included within the article.

References

All Figures

	Figure 1 Mesh, Gaussian curvature, and thickness (mm) over the mesh of reference cornea A.
In the text

	Figure 2 Gaussian curvature over mesh of surface test A for different mesh densities and two levels of noise: SNR = 70 dB and SNR = 90 dB.
In the text

	Figure 3 Original mesh corresponding to surface model A and spectral mesh reconstruction using different numbers of eigenvectors in the spectral decomposition.
In the text

	Figure 4 Shape-DNA distance to reference cornea model A for different scalar functions – Gaussian curvature, mean curvature and thickness – as function of timing surface deformation: explicitly the ho parameter, but implicitly also σx and σx.
In the text

	Figure 5 Surface meshes and Gaussian curvature as function of shape timing deformation: A) ho = 0.264, σx = 0.98, and σy = 0.98; B) ho = 0.332, σx = 1.49, and σy = 1.49; C) ho = 0.4, σx = 2, and σy = 2.
In the text

	Figure 6 Surface, Gaussian curvature and thickness over meshes: A) Reference mesh; B) Reference mesh transformed by an eye fixation change. C) Reference mesh transformed through a uniform scaling.
In the text