J. Eur. Opt. Society-Rapid Publ. 22, 25( 2026) 261
where the sum is over all adjacent vertices j of i, a ij and b 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 ]):
L ¼ B �1 C;
where B is the diagonal matrix whose entries are 1 / 2A i and 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 L is, in general, not symmetric, which makes the eigenvalue problem more demanding. Luckily, 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 ]):
Cf ¼ kBf; ð4Þ
ð3Þ
which provides the same spectral decomposition as finding the eigenvalues of 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: T ¼½ / 1;:::; / k Š, and excluding the smallest eigenvalue / 0 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 ]:
V ¼ T ðT T VÞ; where, recall, that 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) = d i, j. We denote this set of normalized eigenvectors: T n ¼½ / n1;:::; / nk Š. Then, the Laplace – Beltrami spectrum is obtained by a projective matrix operation [ 7 ]:
^f ¼ T
T n 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
ð5Þ 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 n ij ¼ 2n iðp i � p j Þ: ð7Þ jp i � p j j 2
Now, similar to how the first fundamental form is defined, the second fundamental form, denoting n the surface normal, is:
2
b 11 b 12
r z xx z pffiffiffiffiffiffiffiffiffiffiffiffi xy
3 pffiffiffiffiffiffiffiffiffiffiffiffi xxn r xy n 1þz
¼ 2 x þz 2 y 1þz 2 x þz 2 y
4 5 z b 21 b 22 r yx n r yy n pffiffiffiffiffiffiffiffiffiffiffiffi yx z p ffiffiffiffiffiffiffiffiffiffiffiffi yy
:;
1þz 2 x þz 2 y
1þz 2 xþz 2 y
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 1 and k 2, are the eigenvalues of the second fundamental form, the Gaussian curvature is its product: K ¼ k 1 k 2, and the mean curvature its average H ¼ 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-curvatureestimationl-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: x 2
R x þ y2
R y zðx; yÞ ¼piston þ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 1 � ð1þQ xÞx 2
� ð1þQ yÞy 2
R 2 x
R 2 y þ h 0 e �ðx�x0Þ2 2r 2 � ðy�y 0 Þ2 x 2r 2 y
; ð8Þ and the thickness at each point with:
zðx; yÞ ¼0:5mm� h 0 e �ðx�x0Þ2 2r 2 � ðy�y 0 Þ2 x 2r 2 y
: ð9Þ