Laplace Beltrami Filter on QuadEdge Meshes

This document describes a contribution to the Insight Toolkit intended to support the process of statistical analysis in Computational Anatomy. The methods included here are restricted to closed, triagulated surfaces (represented by a QuadEdgeMesh) with a Spherical geometry and topology. The ﬁlter assigns basis function values as Point Data on each vertex of the Mesh. This paper is accompanied with the source code, input data, parameters and output data that we used for validating the algorithm described in this paper. This adheres to the fundamental principle that scientiﬁc publications must facilitate reproducibility of the reported results.


Introduction
In the field of Computational Anatomy the shape of biological structures are compared mathematically between subjects. The difference or similarity in shape between subject structures is quantified and correlations between shape and disease can be produced and examined for statistical significance.
It is desirable to reduce the dimensionality of statistical comparisons between subjects, but maintain information that might be valuable in finding correlations between shape and pathology. A surface mesh can be represented in a way that its N most significant surface harmonics can be compared with other subjects.
The purpose of this filter is to use the Laplace-Beltrami operator to determine surface harmonics in terms of PointData at each vertex. In the same way that a sound signal can be approximately characterized by a combination of its most significant frequency components, so can a surface be expressed as a combination of its surface harmonics. This filter determines the requested N most significant harmonics.

Overview
The filter described in this report operates on itk::QuadEdgeMesh data, and provides the requested basis function as PointData in a copy of the Input Mesh.

Algorithm For Approximation of the Laplacian on Triangulated Surfaces
The basic algorithm implemented in this filter is based on [Qiu2006] and is similar to [Levy2009]. The algorithm visits all the faces of the triangulated mesh, determining each face's area and the areas of the associated vertices, then computes the laplacian operator as a sparse matrix over the vertices.
Let S be a triangulated surface with faces f ∈ F and vertices v ∈ V . Let ψ be a function defined on vertices. We want to define the laplacian of ψ.
First define the gradient, defined as a function indexed by faces. Let e 1 , e 2 , e 3 be three edges forming a face f (with the correct orientation). Let (v 1 , v 2 , v 3 ) be its vertices so that e 1 = v 3 − v 2 , e 2 = v 1 − v 3 and e 3 = v 2 − v 1 . Let c = (e 1 + e 2 + e 3 )/3 be the center of the face.
With this notation, the previous system is Mψ f = G f α; this implies |u| 2 = α T G f α = ψ T f M T G −1 f Mψ f . First note that det G f = |e 1 | 2 |e 2 | 2 − (e 1 · e 2 ) 2 = (|e 1 ||e 2 | sin θ 3 ) 2 where θ 3 is the angle at v 3 . It is therefore |e 1 | 2 e 1 · e 2 e 1 · e 3 e 1 · e 2 |e 2 | 2 e 2 · e 3 e 1 · e 3 e 2 · e 3 |e 3 | 2   after computation and using e 3 = −e 1 − e 2 . Let Σ f denote this last matrix. We can write the approximation: We want to identify the operator ψ → ∆ψ (the discrete Laplacian on S) such that is the edge opposed to v in face f and v and v are the other two vertices in f . This means that one should define One can rewrite this discrete Laplacian in terms of angles. For a vertex v and a face f such that )a( f ) and, since the sum of edges is 0, One can therefore write These expressions allow us to write the Laplacian operator as a sparse matrix indexed over the vertices of the triangulated surface. Letting V = (v 1 , . . . , v N ), define A to be a sparse matrix such that (A is sparse since A(k, l) = 0 if k = l and v k and v l are not connected by an edge.) Besides the zeros, A has the following entries. We first assume that S has no boundary, so that each edge belongs to exactly two faces. • If (k, l) is an edge, we can define β kl and β kl to be the two angles opposed to the edge in the two faces that contain it. Then where l ∼ k indicates that l and k are linked by an edge.
These definitions must be modified when surfaces have boundaries. One usually works with two types of boundary conditions: ψ = 0 on the boundary of S (Dirichlet boundary condition) or ∇ψ tangent to the boundary of S (von Neumann boundary condition).
For triangulated surfaces, Dirichlet boundary condition is directly implemented by restricting to nonboundary vertices. For such vertices, the formulae for A(k, k) and A(k, l) remain unchanged.
Von Neumann boundary condition can be discretized by adding dummy faces to the surfaces symmetrically to boundary edges. This results in an extended surfaceS and a function ψ on S can be extended toS by symmetry too.
The resulting new definition of A(k, l) is identical to the one with closed surfaces if (k, l) is an interior edge. If (k, l) is a boundary edge, then there is only one β kl and A(k, l) = cot(β kl )/a(v k ) (so that the angle is counted twice). A(k, k) is defined also in this case by (Notice that this is not necessarily satisfied with the Dirichlet boundary condition.) Letting D be the diagonal matrix with coefficients D(k, k) = a(v k ), then B = DA is a symmetric matrix (B is just defined like A without the normalizations nu a(v k )). Finding the eigenvalues and eigenvectors of A is equivalent to solving the generalized eigenvalue problem DefiningB(k, l) = B(k, l)/ a(v k )a(v l ) so thatBD −1/2 BD −1/2 , this is equivalent to solving the standard eigenvalue problemBψ = λψ and to let ψ(k) =ψ(k)/ a(v k ).
Finally, notice that the definition can be modified as where a f (v) is the area of the Voronoï¡F9¿ cell of v in f (region in the face for which points are closer to v than to the other two vertices).
The expression of a f (v) depends on whether f is obtuse (has an angle larger than π/2) or not.
In the non-obtuse case, one has (with the same notation as before) In the obtuse case, one must take a f (v) = a( f )/2 if v is the obtuse angle and a f (v) = a( f )/4 otherwise.

Implementation
This filter derives from itk::QuadEdgeMeshToQuadEdgeMeshFilter. The input consists of a triangulated itk::QuadEdgeMesh. Filter Data Generation code calls CopyInputMeshToOutputMesh so the output of the filter is structurally a copy of the input filter, of type TOutputMesh. The filter uses vnl sparse matrix to contain the Laplacian operator over the vertices. To compute eigenvalues and eigenvectors, the filter uses vnl sparse symmetric eigensystem.

Usage
This filter derives from itk::QuadEdgeMeshToQuadEdgeMeshFilter, so the user needs to specify TInputMesh and TOutputMesh to instantiate the class. The input should be a triangulated itk::QuadEdgeMesh based object. Meshes with polys of more than 3 sides will generate an exception in Update(). The user should specify the number of surface harmonics (lbFilter->SetEigenValueCount(eCount)) to generate before updating the filter. Once the filter data is generated (lbFilter->Update()), the user has access to: • GetLBOperator -the sparse matrix of the Laplace operator over the vertices

Results
The results in this example can be obtained by executing the test program itkLaplaceBeltramiFilterTest2: The program will determine the Laplacian Operator values for each matrix. If if the user specifies a harmonic count N with -e N or --eigenvalueCount N, the program will produce N .vtk files containing the original mesh with the PointData set to the first N surface harmonics. The following example computes the first nine surface harmonics on the input mesh fvo.vtk, a hippocampus surface.