Skeletons of 3D Surfaces Based on the Laplace-Beltrami Operator Eigenfunctions

In this work we describe the algorithms to construct the skeletons, simplified 1D representations for a 3D surface depicted by a mesh of points, given the respective eigenfunctions of the Discrete Laplace-Beltrami Operator (LBO). These functions are isometry invariant, so they are independent of the object’s representation including parameterization, spatial position and orientation. Several works have shown that these eigenfunctions provide topological and geometrical information of the surfaces of interest [1] [2]. We propose to make use of that information for the construction of a set of skeletons, associated to each eigenfunction, which can be used as a fingerprint for the surface of interest. The main goal is to develop a classification system based on these skeletons, instead of the surfaces, for the analysis of medical images, for instance.


Skeletons
A curve-skeleton is a 1D model of a 3D object that captures the general characteristics of the original object, and it is also known as centerline. They are useful for visualization and virtual navigation. Another application is registration of 3D objects: given a query object, the task is to find similar objects in a database by using the curveskeleton as a fingerprint. A great variety of algorithms for the generation of skeletons have been developed in recent years [3] [4].
Shape intrinsic information should not depend on the given representation of the object. However, many of the current methods of skeleton construction have the weakness of being sensitive to changes on scale factors, changes in the surface's triangulation, orientation, etcetera. It is desirable that the curve-skeletons have certain properties in order to be used as fingerprints [1]: • Topology preserving. Two objects have the same topology if they have the same number of connected components and cavities. Though a 1D curve cannot have cavities, skeletons must be able to grasp objects characteristics related to its genus. • Scaling invariant. It is necessary for the skeletons to be measurement unit independent; i.e. that it does not depend on the way in which the object could be measured. • Isometry invariant. The skeleton of an object should be independent of the object's given depiction and location. • Rotation invariant. Therefore, checking if two objects are similar needs no prior alignment.
• Similarity. Similar objects should have similar fingerprints.
In particular, it is known, that the eigenfunctions of the Laplace-Beltrami Operator satisfy several of the properties required for a fingerprint [1], for instance: The eigenfunctions depend only on the gradient and divergence which are dependent on the Riemannian structure of the manifold, so they are clearly isometry invariant.
The eigenfunctions are normalizable, therefore, there is no need to concern about scale factors.
Recently, several methods have been developed making use of these eigenfunctions to construct the curveskeletons of objects of interest [5].

Surfaces Representation
Object File Format (.off) files are used to represent the geometry of a model by specifying a triangulation of the model's surface. The OFF files in the Princeton Shape Benchmark [6] conform to the following standard: OFF files are text files. It has a header line with the string OFF. The second line states the number of vertices, the number of faces, and the number of edges; however the number of edges can be ignored for our purpose.
The next lines describe the Cartesian coordinates of each vertex, written one per line. The enumeration of the vertices is given by the order they occurred in the file, starting with 0.
After the list of vertices, the faces are listed; starting with the number of sides and followed by the oriented list of vertices included. All the faces are oriented in the same direction.
The faces can have any number of vertices, although they usually are triangles. For example, Table 1 shows the description of a unitary cube in this format.

Case of Study
Our present case of study are surfaces of rat-hippocampus, obtained from MRI images. On reported works, a relation between morphological changes in the hippocampus and Alzheimer disease in early stages has been found; nowadays there are many studies in image analysis of this and other different brain structures [7].
As many works have shows, the first eigenfunction of the Laplace-Beltrami operator clearly identifies a principal direction of the surfaces of interest, so it is very useful to build a skeleton. Besides that, we noticed that the second eigenfunction reveals additional geometrical information, such as localization of protuberance. Thus, we decided to build a skeleton based on the second eigenfunction also. The construction of both skeletons will be described in detail in section 4. We expect that the properties of these skeletons will provide important information regarding the geometry of objects of interest in order to classify them in a more detailed manner.

Laplace-Beltrami Operator
be a real-valued function defined on a Riemannian manifold M: (1) The Laplace-Beltrami operator is defined as This operator appears in several equations in physics: The method of separation of variables allows us to isolate the spatial dependence of u from the temporal de- , substituting this into the wave equation produces The substitution into diffusion equation leads to the same result. The solution of Equation (6) is a set of eigenvalues 1 2 , , λ λ  and a set of eigenfunctions 1 2 , , f f  which can be thought as fundamental vibration modes.

Finite Element Method
We do not have an equation describing a differentiable variety M, thus we work with a discrete representation of a triangulated surface S described by an OFF file. We need a numerical integration method to solve Equation (6) on S. We have opted to use the Finite Element Method [2].
First of all, we choose N linearly independent form functions 1 2 , , , N F F F  as a basis of a vector space. These base functions : i F S → R are chosen to simplify the calculations, so they are constructed as linear functions that "sample" function f at each vertex of the triangulated surface: The function f then can be written as a linear combination of these base functions The substitution of this approximation in (6) reduces the equation to a generalized eigenvalue problem.
where the entry U i is the contribution of f at vertex i on S. The components of matrix A are given by: and the components of matrix B are: Figure 1(a) and Figure 1(b) show the first and second eigenfunctions, respectively, of Equation (9) applied on a triangulated surface of a rat-hippocampus. The lower values are colored in blue and the higher ones in red. We can see the monotonous behavior of the first function whereas the second function grows from the middle of the surface through its extremes.

The Construction of the Skeletons
The first step is to transform the OFF format into a graph, in this way we can use standard Graph Theory methods. Let be { } , = G V E , the vertices V are indexed as they are in the OFF file, and the edges E can be easily obtained from the list of faces.
We chose the adjacency lists representation for simplicity and efficiency on a triangulated surface with thousands of vertices, usually each one is adjacent only up to five or six vertices; although that number depends on the triangulation and the surfaces. Table 2 shows the adjacency lists of the unitary cube above.
Although the graph is undirected, we allow redundancy of edges in order to simplify the search for local maxima and minima of the eigenfunctions.   We do not want to modify the OFF files, so we use an auxiliary file with the values of the first eigenfunctions calculated for each vertex, this is done with the method described in the previous section. On the surfaces studied, we have observed that each eigenfunction gives different information: • The first eigenfunction f 1 has one maximum and one minimum at opposite points of the surface, these extreme points give us a principal axis, and a main direction (Figure 1(a)). So we decide to use this function to build a polygonal skeleton, which can give us information about curvature and torsion. • The second eigenfunction f 2 has several local maxima and minima at the prominent protuberances of the surface (Figure 1(b)). So we decide to use this function to build a tree-based skeleton that captures the arborescent structure of the surface. So we need two variations of the same algorithm, one for each eigenvalue.

First Eigenfunction Skeleton
Let M and m be the vertices with the absolute maximum and minimum values of the first eigenfunction: The vertices fall between these energy levels, however there are some edges with one vertex in one level and the other in the following, so we define a boundary as a subgraph These boundaries generate a partition of the original graph. We perform a deep-first search to find one vertex in the boundary B i , and a second deep-first search to get all the vertices in it. It is possible that for some energy levels these are so close to each other that the algorithm cannot find a definite boundary; in this case we skip to the next level.
Each boundary B i is a "ring" of vertices, which must be reduced into a centroid c i that can be calculated as a "center of mass": 1 , for 1, 2, , 1; where r j represents the coordinates of the vector of vertex j. Obviously, the boundaries for e m and e M have just one vertex, so 0 m c r = and n M c r = . Finally we connect these centroids to build the skeleton; this last step is performed by the Prim's algorithm [8]. This algorithm finds the minimum cost-spanning tree of the set of centroids c i , using their Euclidean distances as the costs function. Although the skeleton for the first eigenfunction is a polygonal, it can be seen as a degenerated tree (a one-degree tree).

The Second Eigenfunction Skeleton
As we exposed for the first eigenfunction, we search for the absolute maximum and minimum points The search can be easily done in the adjacency lists of E.
We use the coordinates vectors of these vertices as the base of the centroids set:

419
We define the energy levels e i as in (12) and the boundaries B i as in (13). However each one of these boundaries could have k connected components 1 2 , , , , each one has its own centroid 1 2 , , , , as defined in (14). The algorithm for the second eigenfunction process the surface as a tree structure: it performs a deep-first search for a set of boundaries with their respective branch of centroids, and then it performs a second search to find a second branch, an so on. The partition of the surface defined by the first search avoids the algorithm to fall into the same boundaries.
In this case, the Prim's algorithm shows all its performance in the construction of the skeleton, connecting the set of centroids as a minimum-cost spanning tree. In order to avoid that the skeleton cuts the surface, we add directions to the distances between centroids. We associate the information of the energy level to each centroid: The algorithm connects the centroids in increasing order of ( ) h c . Figure 3(a) and Figure 3(b) show the second eigenfunction skeletons of the same surfaces as Figure 2 and

Conclusions and Further Work
We have developed the algorithms for building the skeletons for the first and second eigefunctions of the Laplace-Beltrami operator calculated on an rat-hippocampus surface depicted by an OFF file. The skeleton for the first eigenfunction is built as a polygonal structure along the main axis of the surface. The skeleton for the second eigenfunction has a tree-structure with two main branches, each one with a similar structure to the first skeleton, and small branches for some local maxima at prominent protuberances. This is a first step for the developing of a classification system. The next steps are to generalize these results for different anatomical structures, and define characterization methods for these skeletons based on their geometrical and topological properties, such as critical points of curvature and torsion, bifurcation points, number of branches, etcetera.