Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE Trans Med Imaging. Author manuscript; available in PMC 2011 December 1.
Published in final edited form as:
PMCID: PMC2995840

Robust Surface Reconstruction via Laplace-Beltrami Eigen-Projection and Boundary Deformation

Yonggang Shi, Member, IEEE
Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095, USA (ude.alcu.inol@ihsy)
Rongjie Lai
Department of Mathematics, UCLA, Los Angeles, CA 90095, USA (ude.alcu.htam@jrial)
Jonathan H. Morra, Ivo Dinov, and Paul M. Thompson
Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095, USA (ude.alcu.inol@arrom.nahtanoj; ude.alcu.inol@vonid.ovi; ude.alcu.inol@nospmoht)
Arthur W. Toga, Member, IEEE*


In medical shape analysis, a critical problem is reconstructing a smooth surface of correct topology from a binary mask that typically has spurious features due to segmentation artifacts. The challenge is the robust removal of these outliers without affecting the accuracy of other parts of the boundary. In this paper, we propose a novel approach for this problem based on the Laplace-Beltrami (LB) eigen-projection and properly designed boundary deformations. Using the metric distortion during the LB eigen-projection, our method automatically detects the location of outliers and feeds this information to a well-composed and topology-preserving deformation. By iterating between these two steps of outlier detection and boundary deformation, we can robustly filter out the outliers without moving the smooth part of the boundary. The final surface is the eigen-projection of the filtered mask boundary that has the correct topology, desired accuracy and smoothness. In our experiments, we illustrate the robustness of our method on different input masks of the same structure, and compare with the popular SPHARM tool and the topology preserving level set method to show that our method can reconstruct accurate surface representations without introducing artificial oscillations. We also successfully validate our method on a large data set of more than 900 hippocampal masks and demonstrate that the reconstructed surfaces retain volume information accurately.

Keywords: Surface reconstruction, mask, outlier, Laplace-Beltrami eigen-function, eigen-projection, deformation, topology

I. Introduction

To perform 3D shape analysis of anatomical structures, such as the mapping and analysis of sub-cortical regions in large scale brain imaging studies [1]–[8], a critical problem is the robust reconstruction of a smooth and triangulated surface from segmented volume masks. While this is a well-studied area, conventional solutions are often not satisfactory in the presence of spurious features due to segmentation errors. At the core of the difficulty is the localized filtering of boundary geometry without shrinking other parts of the mask or altering the topology. In this paper, we propose a novel solution to this challenge using iterated Laplace-Beltrami eigen-projection and boundary deformation for localized outlier detection and removal. The surface generated by our method provides an accurate and smooth representation of the boundary geometry and is guaranteed to have the correct topology.

While binary masks are usually sufficient for volume-based studies, smooth surface representations are important for many shape-based analyses of anatomical structures. A smooth curvature map will allow the robust detection of landmarks using the differential geometry of surfaces [9]. The smoothness of the normal directions over the surfaces is also critical for establishing correspondences with currents [10]. Removing spurious outliers from the surface models also improves regularity across population and can be useful for shape prior model construction [3]. It is important to point out that deformable models such as level-set methods [11] can be used to compute smooth surfaces directly from images for many medical problems [12]. The main motivation for our work is the reconstruction of surface models for neuroanatomical structures such as the hippocampus, caudate, and putamen in the sub-cortical region of the human brain, which are typically assumed to have the genus-zero topology. For sub-cortical segmentation, the most successful tools in large scale studies to date have been those based on voxel labeling and their outputs are binary masks [13]–[15]. While spatial smoothness is usually incorporated in automated segmentation, leakage into a neighboring structure can occur and as a result spurious spikes or branches can form on the mask boundary. The binary masks may also come from manual tracing and the difficulty of human tracers in enforcing 3D smoothness can lead to large discontinuities in neighboring slices. To ensure accuracy in the reconstructed surface, the ideal solution should thus eliminate these high frequency outliers while keeping other parts of the boundary intact.

The generation of smooth surfaces from segmentation is a well-studied problem and various solutions exist for different applications. By viewing the boundary of the mask as a series of 2D slices, interpolation-based approaches were popular for the generation of smoothed surfaces for visualization purposes [16]–[18], but the balance between outlier removal and surface volume shrinkage was not considered in these methods and they typically do not guarantee a specific topology. To smooth the mask in the 3D domain, interactive segmentation techniques can be applied [19], [20], but the need of human interaction makes it difficult to use them in large scale studies [21]. Morphological operators such as opening and closing were applied sequentially to modify the boundary [22], [23], but their effect is global and can perturb smooth regions and thus affect the accuracy at those locations. To guarantee the genus-zero topology, topology-preserving deformations [24]–[27] can be applied, but the smoothing force in deformable models such as the mean curvature flow [11] can lead to volume shrinkage. While volume-preserving deformations were proposed [28], such remedies are global and cannot account for the local differences in shrinkage. Diffeomorphic registration was also proposed to generate genus-zero surfaces from segmented masks [29], but it requires a pre-processed template and lacks explicit handling of spurious features since the regularization is global. The smoothing can be applied to a mesh representation of the mask boundary after topology correction is applied to the mask [30]–[33]. Even though popular approaches such as mesh filtering [34] or curvature flows [35] can help smooth the surface, these methods were designed to suppress small scale noise and cannot cleanly remove large outliers. By first mapping the mask boundary to a unit sphere, the spherical harmonics (SPHARM) were applied to provide a least square approximation to the boundary without shrinkage [1], [36]. The SPHARM are most appropriate when low order approximation is satisfactory and become less effective in preserving surface details as artificial oscillations start to appear when higher order basis functions are incorporated.

In this paper, we propose a novel approach that can locally detect and remove spurious features on the mask boundary and reconstruct a smooth surface representation of the anatomical structure with the genus-zero topology. The main contribution of our work is the development of the iterative mask filtering process, as illustrated in Fig. 1, that performs localized outlier detection and removal by going back and forth between the mask and mesh representations of the anatomical region. By projecting the boundary mesh onto the subspace of its Laplace-Beltrami eigen-functions [8], [37]–[44] in the first step, our method automatically locates the position of spurious features by computing the metric distortion in eigen-projection. Using this information, the second step is a mask deformation process that only removes the spurious features while keeping the rest of the mask intact, thus preventing unintended volume shrinkage. This deformation is topology-preserving and well-composed such that the boundary surface of the mask is a manifold. These two steps iterate until convergence and our method generates the final surface as the eigen-projection of the mask boundary, which is a smooth surface with genus-zero topology.

Fig. 1
An overview of our surface reconstruction method.

The rest of the paper is organized as follows. In section II, we introduce the background and necessary concepts for the problem of surface reconstruction from binary masks. In section III, we first introduce the numerical computation of Laplace-Beltrami eigen-functions on triangulated surfaces and then develop the metric for localized outlier detection. In section IV, we develop the well-composed and topology-preserving boundary deformation algorithm for outlier removal. Experimental results are presented in section V to demonstrate the effectiveness of our method and its application to a large data set of over 900 hippocampal masks. Finally, conclusions and future research plans are presented in section VI.

II. Continuous representation of masks

In order to robustly analyze the geometry of the mask boundary with its Laplace-Beltrami eigen-functions, we first need to construct a continuous mesh representation of the boundary. After that, the finite element method can be applied to compute the eigen-functions. In this section, we follow the notations in [45] and develop a continuous representation of the object region boundary, which allows us to move freely between the continuous surface representation and discrete mask representation. This enables us to extract outliers in the mesh representation and feed this information back to the discrete domain for deformation-based outlier removal.

Let B denote a 3D binary mask defined over a lattice of grid points A={(i,j,k)Z30i<I1,0j<J1,0k<K1}. Each point is usually referred to as a voxel in 3D image. The object region is the set of voxels Ao = {(i, j, k) [set membership] A|B(i, j, k) = 1}, and the background region is Ac = A\Ao.

To reconstruct a smooth surface representation of the object boundary, we consider each grid point p = (i, j, k) as the center point of a rectangular cuboid C(p)={(x,y,z)R3xihxhx2,yjhyhy2,zkhzhz2}, where hx, hy and hz are the spatial sampling resolutions in the x, y and z direction determined by the imaging process. The continuous representation of the object and background region are


The boundary of the object and background is then


For numerical processing, we use a triangular mesh representation of the boundary surface. For a cuboid, we can represent its boundary as twelve triangular faces by dividing each of its six rectangular faces diagonally into two triangles. A triangular face F of a cuboid is called a boundary face if it satisfies FC(p)C(q) where p [set membership] Ao and q [set membership] Ac, and we denote p [set membership] Ao as the object neighbor of F and q [set membership] Ac as the background neighbor of F in the lattice. The union of all boundary faces form a triangular mesh representation of the boundary surface and we denote it as M={V,T} where V={Vii=1,,Nv} and T={Tii=1,,Nt} are the set of vertices and triangles. Given the three vertices of any triangle Ti, we can use the definition of the cuboids to determine its object and background neighbor in the lattice uniquely. We denote these relations as two maps ON:TAo and CN:TAc. Because each cuboid has twelve triangular faces, both maps are many to one.

In order to use concepts from differential geometry, our method requires the boundary surface to be a manifold. For arbitrary masks, however, this is not necessarily the case. A binary mask is called well-composed [45] if its boundary M is a manifold, i.e., the surface M locally looks like an open set in the 2D plane. To check if a binary mask is well-composed, we can use the following proposition from Latecki [45].

Proposition 1: The mask B is well-composed iff the following two configurations do not occur for the continuous representation C(Ao) and C(Ac). (1) Four cuboids share an edge. Two of them that do not share a face are in C(Ao) and the other two are in C(Ac). (2) Eight cuboids share a corner. Two of them are in C(Ao)(C(Ac)) and they are corner adjacent but not edge/face adjacent, and the other six cuboids are in C(Ac)(C(Ao)).

For input masks that do not satisfy this requirement, perturbations can be introduced to make them well-composed [46]. In this work, we develop in section IV-B an efficient boundary deformation algorithm to construct an accurate and well-composed approximation to the segmented mask, which enables us to assume the well-composedness of the mask in the next section such that shape analysis techniques based on the eigen-functions of manifolds can be developed.

III. Laplace-Beltrami Eigen-Projection

Based on the manifold representation of the mask boundary, we use spectral analysis to study its geometry and robustly detect spurious outliers to be removed by boundary deformation. In this section, we first introduce the eigen-functions of the Laplace-Beltrami operator on the boundary surface M. After that, a novel outlier detection metric is developed based on the eigen-projection of this surface.

A. Laplace-Beltrami eigen-functions

Given a manifold M, the eigen-functions of its Laplace-Beltrami (LB) operator ΔM is defined as [47]


where λ is the eigenvalue and f is the corresponding eigen-function. The spectra of ΔM is discrete and we can order the eigenvalues according to their magnitude as 0 = λ0 ≤ λ1 ≤ λ2 ≤…. For λi, we denote the corresponding eigen-function as Fi.

There have been increasing interests in using the LB spectra to study 3D shapes in computer vision [37], [48], medical shape analysis [8], [38], [40], [42]–[44], and graphics [39], [41]. For the problem of surface reconstruction from masks, there are several reasons that we choose to use the LB eigen-functions. First of all, these functions are intrinsically defined and can be easily computed from the boundary surface with no need of any parameterizations. Secondly, the LB eigen-functions are isometry invariant, and thus, are robust to the jagged nature of the boundary surface. Last, but not least, the magnitude of the eigenvalues of the LB operator intuitively corresponds to the frequency in Fourier analysis, thus it provides a convenient mechanism to control the smoothness of the reconstructed surface.

To numerically compute the eigen-functions, we use the triangular mesh representation of M. In the discrete form, the eigen-function is defined on the vertices and we denote f as a vector of length Nv. We can then compute f by solving a generalized matrix eigenvalue problem:


where the two matrices Q and U are formed using the finite element method [49]. More specifically, the matrices are defined as:



where N(Vi) is the set of vertices in the 1-ring neighborhood of Vi, N(Vi,Vk) is the set of triangles sharing the edge (Vi,Vj), θli,j is the angle in the triangle Tl opposite to the edge (Vi,Vj), and Al is the area of the l-th triangle Tl2. To solve the eigenvalue problem, we use the software package ARPACK [50] and SuperLU [51].

B. Outlier detection via eigen-projection

Given a function g:MR, its eigen-projection onto the space of the first K eigen-functions is defined as:


In the discrete form, both the eigen-functions and the signal g are defined on the vertices and we treat them as vectors of size NV×1. By representing the first K eigen-functions as a matrix


we can write the eigen-projection of g as


where the matrix U takes into account the integral on the triangular mesh.

By considering each coordinate of the vertices as a function on M, we can define the eigen-projection of the boundary surface. Let x, y, and z denote the vector of the x, y, and z coordinates of all vertices in M. The i-th element of the vectors represent the coordinate of the i-th vertex, i.e., Vi=(xi,yi,zi). The eigen-projection of the coordinate vectors are denoted as x^,y^ and z^, respectively. By replacing the coordinate of each vertex with these projected coordinates, we denote the eigen-projection of the surface M as the mesh M^=(V^,T^) where V^={V^ii=1,,Nv} with V^i=(x^i,y^i,z^i) and T^=T.

To illustrate the effect of eigen-projection, we show in Fig. 2 the projection of the boundary surface of a manually segmented caudate nucleus onto the first 300, 100, and 50 eigen-functions. With the decrease of the number of eigen-functions, we can see the eigen-projection generates a surface that becomes smoother in most places. This is because the effect of eigen-projection on every coordination function is analogous to the low-pass filtering of a one-dimensional signal defined on the surface. However, sharp outliers such as the cusp and wedge highlighted in Fig. 2(a) are still clearly visible in Fig. 2(d) even when only 50 eigen-functions are used. This difficulty in removing sharp outliers via eigen-projections results from the challenge that the smoothness of the coordinate functions does not guarantee the smoothness of the geometry. In the graphics literature [41], such effects can be observed when eigen-projections are used to process animated human and animal shapes. This difficulty is also general for other basis functions such as spherical harmonics since it is well-known that high curvature features such as cusps can be formed even with low order spherical harmonics [52].

Fig. 2
The eigen-projection of a caudate mask.

While eigen-projection itself does not smooth out all the outliers, it offers an effective approach of detecting the outliers on the jagged mask boundary. From the result in Fig. 2, we can see clearly that the cusp and wedge become sharper after projection as the projection tries to squeeze them into a relatively smaller area on the projected surface M^. This motivates us to define the following area distortion factor(ADF) on each triangle of M to locate sharp features on the boundary surface:


where Al and A^l are the area of the corresponding triangle Tl and T^i. Because the ADF is determined by M and M^, it is a function of K, which is the number of eigen-functions used to generate the projected surface. One nice property of the ADF is that it is invariant to rotation, translation, reflection and scale differences. This makes it easy to choose a robust threshold that is applicable for the detection of sharp outliers on a wide variety of shapes. As an illustration, we show in Fig. 3 the ADF between the boundary surface in Fig. 2(a) and its eigen-projection in Fig. 2(b). Except for the sharp features, we can see the ADFs in most parts of the surface are very small. This shows that the ADF can successfully localize the sharp features and differentiate them from the rest of the surface. Because the ADF is derived from the eigen-functions, it is determined by the global geometry of anatomical structures that are typically stable across population and disease groups. Thus, it gives us a robust way of outlier filtering suitable for surface reconstruction in large brain mapping studies.

Fig. 3
The ADF of the boundary surface after the projection onto the first 300 eigen-functions.

Using the mapping function ONa and CN, we can define the ADF for lattice points on the boundary of the object and background region. Given an object region Ao and background region Ac, we first define two sets:


Here N6(·) denotes six-connected neighbors of a voxel, and the set Lin and Lout represent the voxels on the interior and exterior boundary of Ao, respectively. The six-connectedness is used in defining the neighborhood because it is a necessary condition for the well-composedness of the object and background region. For points in these two sets, we define their ADF as


Thus the ADF at a boundary voxel is the maximal area distortion experienced by any of its boundary faces during the eigen-projection. This feeds the outlier information computed with the mesh representation of the boundary back to the lattice domain. This allows us to design localized speed functions to drive the boundary deformation algorithm we develop in the next section such that outliers can be filtered out in the lattice domain.

IV. Outlier Removal with Boundary Deformation

Using information derived from the eigen-projection of the mask boundary and the relation between its continuous and discrete representation, we can filter out the outliers via boundary deformation without affecting other parts of the boundary. In this section, we first extend a fast algorithm for boundary deformation by incorporating well-composedness and topological constraints. For an arbitrary input mask, we use this algorithm to build a well-composed approximation such that eigen-projection can be applied to its boundary surface. After that, the deformation algorithm can be successfully applied to remove outliers and compute a smooth surface with genus-zero topology using an ADF-based speed function.

A. Well-composed and topology-preserving deformation

The deformation algorithm we develop here is an extension of the fast algorithm that approximates level-set-based curve evolution proposed in [53]. To ensure the boundary surface is a manifold, we incorporate well-composedness into the deformation algorithm according to Proposition 1. For many medical shape analysis problems, especially the study of brain structures, the reconstructed surface must have the genus-zero topology, thus we also include the topology-preserving ability into the deformation algorithm.

Starting from an initial mask, our method iteratively adds or removes voxels from the mask according to a speed function and their effect on the well-composedness and topology. Let A^o denote the initial guess of the object region on the lattice A that is well-composed and meets the topological constraint and its complement as A^c=AA^o. To evolve A^o, we construct the two lists Lin and Lout as defined in (9) and the following function for fast access of the regional information about the voxels:


Given a speed function F on the lattice, we can then run the fast evolution algorithm in Table I that we adapt from [53] to iteratively deform the mask until the stopping condition is satisfied. During each iteration, this algorithm updates the two lists according to the speed if the change neither violates the well-composedness nor changes the topology. To check whether the addition or removal of a voxel p from the two lists will change the well-composedness of A^o, we use its 26-connected neighborhood N26(p) in the lattice A and check if the configurations in Proposition 1 will occur. To ensure there is no topological change during the deformation process, we follow the work of topology-preserving level-set method by only removing or adding simple points from the two lists [27].

Well-Composed and Topology Preserving Deformation Algorithm

The design of the speed function is application dependent. Next we design two different speed functions to first reconstruct a well-composed approximation of an arbitrary input mask and then iteratively remove outliers detected via eigen-projections.

B. Well-composed mask with genus-zero topology

To apply the eigen-projection and outlier detection algorithm, we need the mask boundary to be a manifold. For sub-cortical structures, genus-zero topology is also required. Here we apply the deformation algorithm in Table I to construct a well-composed and genus-zero approximation to an arbitrary input mask.

At the start of the deformation algorithm, we choose the initial mask A^o as the bounding box of Ao such that it meets the requirement of being well-composed and having no holes and handles. To evolve the initial region A^o toward the boundary of Ao, we define the speed function as


which tells the boundary evolution algorithm to move the boundary inward if it is outside Ao and outward if it is inside Ao. As described in Table I, the evolution of a point on the boundary is determined by both the speed function in (12) and the condition of being a well-composed and simple point, which ensures the genus-zero topology [27]. The stopping condition used here is that no voxels change membership during one iteration. Once the algorithm converges, we obtain a well-composed mask with genus-zero topology as the input to the outlier detection and removal process.

As an illustration, we show in Fig. 4 the deformation process that starts from a bounding box and converges to an accurate approximation of the caudate nucleus shown in Fig. 2. To measure the difference of the approximation A^o and the input Ao, we computed the symmetric difference of the two sets A^o and Ao and it only accounts for 0.00023% of the volume of Ao, which clearly shows the accuracy of the approximation.

Fig. 4
The deformation from the bounding box to the well-composed approximation with genus-zero topology.

C. Outlier removal

Using the well-composed mask A^o as the input, we next develop a new speed function to iteratively remove outliers from the mask and build a smooth surface representation of the object boundary. Let M denote the mesh representation of the boundary of A^o. To define the speed function at each iteration, we first apply the eigen-projection process to M and compute the ADF for the two lists Lin and Lout of A^o. Using this information, we define the speed function below to remove the outliers via boundary deformation.


The threshold α is chosen to differentiate outliers from smooth parts of the surface. The smaller α is, more regions on the surface will be considered outliers. For sub-cortical or similar structures, the threshold α is typically chosen between 5 and 10 and this gives robust performances in our experience. The function H is defined as


to determine whether locally the boundary is convex or concave at the current voxel, where N27(p) denotes the 3 × 3 × 3 neighborhood centered at p, and H is the Heaviside function that is zero for negative argument and one otherwise. This simple rule is designed according to the diffusion-based approximation of the mean curvature flow [54] which showed that the mean curvature evolution of binary masks could be approximated by convolving the mask with symmetric kernels. The definition of H can be viewed as the convolution of the mask with a constant kernel of size 3 × 3 × 3. Intuitively this is how this speed can remove outliers. If H(p)12 for a point p [set membership] Lin, we decide the outlier at p is convex and it should be removed from the object region, thus a negative speed is assigned. On the contrary, if H(p)<12 for a point in Lout, we determine the outlier at p is concave and should be removed by moving the boundary outward with a positive speed.

By plugging the speed into the deformation algorithm in Table I, we can iteratively detect and remove outliers as illustrated in Fig. 1. At step 2 of the algorithm, the eigen-projection is applied to the mesh representation of the boundary surface and the ADF is computed for points in Lin and Lout. Using this information, the speed F is also computed. The boundary deforms according to the speed at step 3 to remove outliers from the mask. These two steps of eigen-projection and boundary deformation are repeated iteratively until convergence. As in section IV-B, the same stopping condition that no voxels change membership during step 3 is used at step 4. Once we have the final mask, we generate the reconstructed surface as the eigen-projection M^ of the converged mask boundary.

V. Experimental Results

In this section, we present experimental results to demonstrate our surface reconstruction algorithm. In the first experiment, we apply our method to the manually segmented masks of three sub-cortical structures and illustrate the impact of parameters on reconstruction results. In the second experiment, we test our method on different input masks of the same structure and measure the improvement of consistency between the filtered masks. In the third experiment, we compare with the SPHARM method and demonstrate that our approach can generate smoother and more naturally looking surfaces. In the fourth experiment, we compare our method with the SPHARM tool and the topology preserving level set(TPLS) method in reconstructing hippocampal surfaces from a large data set of more than 900 automatically segmented masks to validate its robustness and accuracy.

A. Reconstruction of sub-cortical surfaces

The input data for the first experiment are the manually segmented masks for a caudate nucleus, putamen, and hippocampus as shown in Fig. 5, where the caudate nucleus is the same as in Fig. 2(a). From the data visualized in Fig. 5, we can clearly see the jagged nature of the masks and the existence of outliers. Our goal is to reconstruct smooth surface representations of these structures and demonstrate the robustness of our method to outliers. We will present results obtained with various combinations of parameters and illustrate their effects on the reconstructed surfaces.

Fig. 5
The manually segmented masks of three sub-cortical structures.

There are two parameters in our algorithm that affect the quality and speed of the reconstruction: the number of eigen-functions K in eigen-projection and the threshold α for outlier detection. For all three masks, we tested 14 set of parameters with the combination of K = 100, 150, 200, 250, 300, 350, 400 and α = 5, 10. For all parameter selections, our method successfully computed the reconstructed surfaces. In general the computational cost is higher with the increase of K and the decrease of α. For the 14 parameter sets, the computational time of our method, which was implemented in C++, ranged from 2 to 10 minutes on a PC with a 1.6 GHz CPU.

To validate the outlier removal process only affects the mask boundary locally around the detected outliers, we computed the fraction of the grid points on the boundary of the input mask, i.e., the points in Lin as defined in (9), that stayed static during the deformation process. For all parameter selections, the percentage of static boundary points on the three masks are plotted as a function of K for each α in Fig. 6. Overall, we can see less boundary points are moved for larger α and K. Even for the smallest parameter set tested here (K = 100, α = 5), we can see around 85% of the boundary points are not moved. This clearly demonstrates the localizing effect of the outlier filtering process in our method.

Fig. 6
The percentage of static boundary points of the input mask after outlier removal.

In Fig. 7, ,88 and and9,9, we visualize eight representative results from the 14 parameter sets for the three structures. By comparing the results in Fig. 7, ,88 and and99 to the input masks in Fig. 5, we can see that our method is able to remove the outliers, such as the cusp and wedge highlighted in Fig. 2(a) for the caudate and the large inter-slice discontinuity in the head region of the putamen and hippocampus, and reconstruct smooth surfaces. Given a specific threshold α, we can see our method gradually adds more detail to the reconstructed surface with the increase of K. For a fixed K, more outliers were removed with the decrease of α.

Fig. 7
The reconstructed surface of the caudate nucleus with eight sets of different parameters.
Fig. 8
The reconstructed surface of the putamen with eight sets of different parameters.
Fig. 9
The reconstructed surface of the hippocampus with eight sets of different parameters.

In summary, we demonstrated that our method was able to successfully remove outliers and reconstruct smooth surfaces from jagged input masks. By choosing the set of parameters properly, our results showed that we can achieve a good balance between outlier removal and fidelity to the input data.

B. Validation with different masks of the same structure

As shown in the first experiment, our method can reconstruct smooth surfaces from jagged input masks. To check the accuracy of the reconstructed surface models, however, is an open problem due to the lack of ground truth. To partially overcome this difficulty, we apply our method to different input masks of the same structure in the second experiment and test the validity of our method by investigating the change of consistency between masks after the removal of random outliers. We used two different data sets in this experiment. In the first data set, two different manual segmentations are available for the left caudate and putamen of five subjects. In the second data set, only one set of manual segmentations are available and we compare them with automatically segmented results.

We performed two tests in the second experiment. In the first test, the input data are the left caudate and putamen of five subjects that were segmented manually by two different tracers. For each structure, we computed the Dice coefficient [55] of the segmented masks from the two tracers and the result is shown in Fig. 10. The iterative outlier removal algorithm of our method was applied to the two segmented masks of each structure. Two sets of parameters (K = 300, α = 10) and (K = 300, α = 5) were used to test the impact of parameters. Given the two filtered masks of each structure, we computed their Dice coefficient and the result is plotted in Fig. 10. From the results in Fig. 10(a) and (b), we can see that the filtered masks have higher Dice coefficients than the original masks from the two tracers. This shows that the removal of outliers help improve the consistency between the segmented masks from these two tracers. We showed in the first experiment that the decrease of the parameter α helped remove more outliers. By comparing the Dice coefficients obtained by the two sets of parameters, we observe a similar trend and the masks filtered by the second set of parameters have higher inter-tracer consistency than the results from the first set of parameters.

Fig. 10
Dice coefficients of masks from two different tracers.

In the second test, the input data are the manually and automatically segmented [15] right hippocampi from five subjects. We assume in this test the manually segmented masks are the gold standard and want to test if the application of the mask filtering process would help improve the agreement between automatically segmented results and the assumed ground truth. For each pair of automatically and manually segmented masks, the Dice coefficient was computed and plotted in Fig. 11. We applied our method to the automatically segmented masks with two sets of parameters: (K = 300, α = 5) and (K = 300, α = 10). The Dice coefficients of the filtered and manually segmented masks were also computed and plotted in Fig. 11. The results from the two sets of parameters are almost the same for this example. While the degree of improvement varies across subjects, we can see that the outlier removal process in our method helps improve the consistency between these automatically segmented masks and manual results.

Fig. 11
Dice coefficients of automatically and manually segmented masks.

By generating smooth approximations of segmented masks, we can see from these test results that our method improves the overall consistency of the segmentation results of sub-cortical structures. This demonstrates the value of removing outliers such as spikes and wedges for high quality surface reconstruction. On the other hand, the results should also be interpreted cautiously because the size of the data sets used in the tests is small. The encouraging results, however, makes it an important direction in future work to further quantify the impact of surface reconstruction on modeling accuracy with large data sets.

We only used the filtered masks generated by our method in this experiment. To compare with existing surface reconstruction methods, we will focus on the triangular mesh representation of the reconstructed surfaces in the next two experiments. This will allow us to compare surfaces in terms of their geometric features. We will also study the fidelity of the surfaces in retaining volume information.

C. Comparison with SPHARM

In this experiment, we compare our method with the SPHARM tool [1], [36] that is publicly available and well-known in the neuro-imaging community. Starting from a binary mask, the SPHARM tool first maps the boundary surface to a unit sphere and generates a smooth surface with spherical harmonics, which is represented as a triangular mesh of genus-zero topology.

To control the smoothness of the reconstructed surface, the SPHARM tool allows the selection of the highest order of spherical harmonics l used in the reconstruction, which means the number of basis functions is K = (l + 1)2. For the three masks shown in Fig. 5, we applied the SPHARM tool with l = 9 to 19 such that the number of basis functions used in the SPHARM reconstruction matches the range of the number of eigen-functions used in the first experiment. The spherical mesh used here for the SPHARM reconstruction is the icosahedron subdivision of the unit sphere at the division rate 20 with 8000 triangles. For each structure, we show four representative SPHARM results with the order l = 9, 13, 16 and 19 in Fig. 12, ,13,13, ,14.14. From these results we can see the SPHARM tool is able to generate smooth surfaces with low order spherical harmonics such as l = 9. But if we try to recover more details in the anatomical boundary by increasing l, artificial oscillations begin to appear in some portions of the reconstructed surface while the surface in general still appears overly smooth. The degree of artificial oscillations becomes more severe as more spherical harmonics are used in the reconstruction. By comparing the SPHARM results with the surfaces in Fig. 7, ,88 and and9,9, we can see our method is able to generate more naturally looking surfaces with the use of high order basis functions. Next we develop a metric to quantify this difference.

Fig. 12
The SPHARM reconstruction results for the caudate nucleus.
Fig. 13
The SPHARM reconstruction results for the putamen.
Fig. 14
The SPHARM reconstruction results for the hippocampus.

Let k:MR denote the mean curvature of a surface M. We define the mean curvature nodal length (MCNL), which is the length of the zero level-set contours of κ on the surface, as a metric to quantify how oscillatory a surface is and use it to compare the performance of different surface reconstruction results. Conceptually this metric can be viewed as an extension of the frequency of 1D sinusoids. The degree of oscillation for sinusoids is defined via counting the number of zero-crossings over a unit interval. Similarly, the oscillation on the surface is characterized by the transition between regions with opposite signs of mean curvature. The difference is that the zero-crossing of the function κ is a set of contours and its size is the length of these contours. As an illustration, we plotted the zero level-set contours on four surfaces in Fig. 15. The surface in Fig. 15(a) was reconstructed with our method using the parameter set K = 250 and α = 5, and its MCNL is 537.5mm. The surfaces in Fig. 15(b), (c) and (d) were reconstructed using SPHARM with l = 15, i.e., K = 256 basis functions. Because the SPHARM result is also affected by the spherical mesh used for reconstruction, three different spherical meshes were used to reconstruct the surfaces in Fig. 15(b), (c) and (d). The first mesh is the icosahedron subdivision of the unit sphere at the division rate 20 with 8000 triangles and the reconstruction result is shown in Fig. 15(b), whose MCNL equals 774.2mm. The second spherical mesh was also obtained by the icosahedron subdivision, but we increased the devision rate to 30 so it has 18000 triangles. The reconstruction result with this denser mesh is shown in Fig. 15(c) and its MCNL is 826.7mm. For the third spherical mesh, we used the spherical parameterization of the boundary mesh of the binary mask, which has 6560 triangles. The reconstruction result of this mesh is shown in Fig. 15(d) and its MCNL is 814.2mm. This shows the MCNL gives an intuitively correct measure that the surface in Fig. 15(b), (c) and (d) have a higher degree of artificial oscillations than that of Fig. 15(a).

Fig. 15
The mean curvature map on two surfaces and their zero level-set contours. (a) The caudate surface reconstructed with our method using the parameter (K = 250, α = 5). (b) The caudate surface reconstructed with SPHARM using the parameter l = 15( ...

For each mask in Fig. 5, we computed the MCNL for surfaces reconstructed with our method and the SPHARM tool. We plotted the MCNL with respect to the number of basis functions K in Fig. 16. For our method, K is the number of LB eigen-functions and we plotted the MCNL with respect to K for each different selection of the parameter α. For the SPHARM tool, K is the number of spherical harmonics determined by the parameter l. Three spherical meshes were used in computing the SPHARM reconstruction: the icosahedron subdivision at the division rate 20 and 30, and the spherical parameterization of the boundary mesh. While the results in Fig. 16 (b) show that spherical parameterization of the boundary mesh gives lower MCNL than the other two spherical meshes when l is small, overall their oscillation level as measured by the MCNL are comparable for the three structures. From all the plots in Fig. 16, we can see the MCNL of surfaces reconstructed with our method using the parameter set (K = 300, α = 5) are lower than all the SPHARM results. This shows that our method is able to reconstruct less oscillatory surfaces when a comparable number of basis functions are used, and the difference is especially significant when high order reconstructions are desired.

Fig. 16
The MCNL of surfaces with respect to the number of basis functions used in the reconstruction.

To further explore the reason for this difference, we investigated the impact of the different basis functions used in our method and the SPHARM tool, i.e., the LB eigen-functions and the spherical harmonics, by fixing the input masks to both the LB eigen-projection and the SPHARM tool. For each shape in Fig. 5, we chose the input mask for both the SPHARM tool and LB eigen-projection as the filtered mask generated by our outlier removal process using the parameter set K = 300 and α = 5. Given the same input mask, the SPHARM tool was applied with the set of parameters l = 9 to 19. Three spherical meshes were used as above to test the effect of different meshes. Similarly, the LB eigen-projection was applied to the boundary surface with K = 100, 150, 200, 250, 300, 350, 400. For all the results reconstructed with our method and the SPHARM tool, we computed their MCNL and plotted them in Fig. 17(a), (b) and (c) for each anatomical structure. These plots clearly show that spherical harmonics generate more oscillatory surfaces than the LB eigen-functions.

Fig. 17
The MCNL of surfaces with respect to the number of basis functions used in the reconstruction.

In summary, our method is able to produce smooth and more detailed surfaces than the SPHARM tool. Our experiment also suggests that the use of the intrinsically defined LB eigen-functions can help avoid the artificial oscillation in the SPHARM results and produce more naturally looking surfaces.

D. Application to a large data set

In this experiment, we apply our algorithm to a large data set and demonstrate its robustness. The input data are the masks of 926 automatically segmented right hippocampi [15] from the baseline and follow-up MRI scans of 145 normal controls, 230 patients of mild cognitive impairment (MCI), and 88 patients with Alzheimer's disease (AD) in the ADNI data [21].

For all the 926 masks, we used the same parameter K = 300 and α = 5 and our method successfully reconstructed smooth surface representations for all the masks. As an illustration, we plotted three examples in Fig. 18 by overlaying the surfaces with the input masks. We can see that our method is able to remove the outliers while producing a smooth surface that accurately represents other parts of the mask boundary.

Fig. 18
The reconstructed right hippocampal surface overlaid with the corresponding mask of three subjects.

To demonstrate that the smooth surfaces reconstructed by our method accurately approximate input masks, we first computed the Dice coefficient of each input mask and its filtered mask generated by our method. The cumulative distribution function (CDF) of the Dice coefficients is plotted in Fig. 19. This clearly shows the filtered masks overlay very well with the input masks after the outlier removal. Furthermore, we have computed the ratio between the volume difference of the reconstructed surface and the input mask and the mask volume. The CDF of this volume difference ratio is plotted in Fig. 20(a). The mean and standard deviation of the volume difference ratio is 0.00098 and 0.00086, and the maximum is 0.0052. For each surface, we computed its MCNL and the CDF of the MCNL is plotted in Fig. 20(b). As a comparison, we also applied the SPHARM tool to reconstruct all surfaces from the input masks. To use a comparable number of basis functions, the order was selected as l = 16. According to the result in Fig. 16(c), we have chosen here the icosahedron subdivision of the unit sphere at the division rate 20 as the spherical mesh because it achieved the least oscillation for the hippocampus in the first experiment. The CDF of the volume difference ratio and MCNL of the SPHARM results are plotted in Fig. 20(a) and (b). We can see from Fig. 20 that the reconstructions generated by our method are not only less oscillatory but also have less volume distortion.

Fig. 19
The CDF of the Dice coefficients between the input masks and the filtered masks generated by our method.
Fig. 20
The CDF of volume difference ratio and MCNL of results generated by our method and the SPHARM tool.

We have also applied the topology preserving level set (TPLS) method [27] for surface reconstruction that is available in the LONI pipeline [56], [57], where a curvature weight is used for smoothness regularization. The mesh representation of the reconstructed surface is extracted from the level-set function with a topologically consistent marching cube algorithm [27]. As an example, we show in Fig. 21(a) the overlay of an input mask with the reconstructed surface from the TPLS method. The curvature weight in this example is 0.3, which was chosen to be just high enough to remove the outlier. The nodal set of the mean curvature of the surface is plotted in Fig. 21(b), which clearly shows the extensive presence of oscillatory regions. For the 926 hippocampal masks, we applied the TPLS method with three curvature weights 0.3, 0.5 and 1.0. The CDF of the volume difference ratio and MCNL for results generated by the three weights are plotted in Fig. 22(a) and (b). Because the curvature weight is applied globally regardless of the presence of outliers, the TPLS results tend to have a shrinkage effect. From the result in Fig. 22(a), we see the volume difference ratio becomes larger with the increase of the curvature weight. Comparing the TPLS results in Fig. 22 with Fig. 20, we can see the surface models generated by both our method and the SPHARM tool are less oscillatory and have smaller volume distortion.

Fig. 21
An example of TPLS reconstruction. (a) The overlay of the input mask and reconstructed surface. (b) The nodal set of the mean curvature.
Fig. 22
The CDF of volume difference ratio and MCNL for TPLS results.

To demonstrate the impact of different surface reconstruction methods on brain mapping studies, we performed a simple test for shape differences between MCI and AD groups by using the volume-normalized MCNL of the reconstructed surfaces as a signature. By factoring out volume differences, we can study the shape differences of the hippocampi across population. Because the MCNL is computed from the mean curvature of a surface, we also demonstrate in this example the impact of surface reconstruction on computing the differential geometry of the anatomical region. The surfaces reconstructed in this experiment with our method, the SPHARM tool, and the TPLS method with the curvature weight chosen as 0.3 were used for comparison. For each method, a p-value was computed by applying the t-test to the MCNL of the surfaces in the MCI and AD group. The p-values of the three methods are listed in Table II. We can see that only the p-value from our method is below 0.05, which is the typical threshold for significance. The result shows that the differential geometry of surfaces with less artificial oscillations gives better separation of the two groups. This suggests the potential of the surface models reconstructed by our method in generating useful geometric signatures for population studies.

P-values of volume-normalized MCNL between the MCI and AD group.

VI. Discussion and Conclusion

We developed a novel approach to reconstruct genus-zero surfaces from segmented masks. Using the metric distortion in LB eigen-projection, our method detects spurious features on the mask boundary and performs localized outlier removal via well-composed and topology-preserving deformation. As a result, our method can remove the outliers without moving other parts of the boundary and generate a smooth, while faithful, surface representation of the input mask. Compared to the SPHARM tool, we showed that our method was able to preserve more geometric details in the reconstructed surfaces without introducing artificial oscillations. We also validated its robustness on a data set of more than 900 masks. In our current research, we are applying this method to a variety of structures in the brain. Using the reconstructed surfaces, we can compute detailed mappings and perform large scale studies that can be valuable for various clinical applications such as the generation of more sensitive bio-markers for the early detection of the Alzheimer's disease.

In our current work, we have used visual observation to select the parameters K and α in our method. For typical sub-cortical regions, this is not difficult since this process only needs to be performed once for each type of structure. On the other hand, the parameter set can be customized automatically for a specific segmentation algorithm if high quality manual segmentation data are available as training data. By minimizing the approximation error on this training set, the optimal parameter set can be determined.

Even though we developed our method for the reconstruction of anatomical structures with spherical topology, it can be easily extended to reconstruct high genus surfaces. Because the outlier detection algorithm is derived from the LB eigen-functions, it is valid for general manifold. Thus, the only change required is removing the simple point constraint in Table I, such that topological changes can occur freely in the boundary deformation process.

Because our method is based on the computation of eigen-functions on the mask boundary, the computational cost can be high for large structures. To overcome this difficulty, we will investigate multi-scale representations in our future research. The spectral shift technique can also be incorporated to parallelize the numerical computation of a large number of eigen-functions [41].

In our future work, we will also study the extension of our method to reconstruct a surface from multiple masks segmented with different algorithms. This can be viewed as a problem of computing the weighted average of multiple masks. The ADF of each mask naturally provides a confidence measure at each boundary point. By using the fast marching algorithm [58], [59], we can compute the distance transform of each mask and extend the confidence measure from the boundary to the whole domain. An evolution speed can then be designed to incorporate the weights from different masks and drive an initial mask toward a weighted average of multiple segmented masks, which can then be used to generate a smooth surface via the LB eigen-projection.


This work was funded by the National Institutes of Health through the NIH Roadmap for Medical Research, Grant U54 RR021813 entitled Center for Computational Biology (CCB). Information on the National Centers for Biomedical Computing can be obtained from


[1] Brechbühler C, Gerig G, Kübler O. Parameterization of closed surfaces for 3-D shape description. CVGIP: Image Understanding. 1995;61(2):154–170.
[2] Gerig G, Styner M, Jones D, Weinberger D, Lieberman J. Shape analysis of brain ventricles using SPHARM. Proc. Workshop on Mathematical Methods in Biomedical Image Analysis. 2001:171–178.
[3] Davies RH, Twining CJ, Allen PD, Cootes TF, Taylor CJ. Shape discrimination in the hippocampus using an MDL model. Proc. IPMI. 2003:38–50. [PubMed]
[4] Thompson PM, Hayashi KM, de Zubicaray GI, Janke AL, Rose SE, Semple J, Hong MS, Herman DH, Gravano D, Doddrell DM, Toga AW. Mapping hippocampal and ventricular change in Alzheimer disease. NeuroImage. 2004;22(4):1754–1766. [PubMed]
[5] Wang L, Miller JP, Gado MH, McKeel DW, Rothermich M, Miller MI, Morris JC, Csernansky JG. Abnormalities of hippocampal surface structure in very mild dementia of the alzheimer type. NeuroImage. 2006;30(1):52–60. [PMC free article] [PubMed]
[6] Shi Y, Thompson PM, de Zubicaray G, Rose SE, Tu Z, Dinov I, Toga AW. Direct mapping of hippocampal surfaces with intrinsic shape context. NeuroImage. 2007;37(3):792–807. [PMC free article] [PubMed]
[7] Yeo B, Sabuncu M, Vercauteren T, Ayache N, Fischl B, Golland P. Spherical demons: Fast surface registration. Proc MICCAI. 2008:745–753. [PMC free article] [PubMed]
[8] Shi Y, Morra J, Thompson PM, Toga AW. Inverse-consistent surface mapping with Laplace-Beltrami eigen-features. Proc. IPMI. 2009:467–478. [PMC free article] [PubMed]
[9] Thirion JP. The extremal mesh and the understanding of 3D surfaces. Int'l Journal of Computer Vision. 1996;19(2):115–128.
[10] Vaillant M, Glaunes J. Surface matching via currents. Proc. IPMI. 2005:381–392. [PubMed]
[11] Osher S, Sethian J. Fronts propagation with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. Journal of computational physics. 1988;79(1):12–49.
[12] Tejos C, Irarrazaval P, Cárdenas-Blanco A. Simplex mesh diffusion snakes: Integrating 2D and 3D deformable models and statistical shape knowledge in a variational framework. Int'l Journal of Computer Vision. 2009;85(1):19–34.
[13] Fischl B, Salat DH, Busa E, Albert M, Dieterich M, Haselgrove C, van der Kouwe A, Killiany R, Kennedy D, Klaveness S, Montillo A, Makris N, Rosen B, Dale AM. Whole brain segmentation: Automated labeling of neuroanatomical structures in the human brain. Neuron. 2002;33(3):341–355. [PubMed]
[14] Tu Z, Narr K, Dollar P, Dinov I, Thompson P, Toga A. Brain anatomical structure segmentation by hybrid discriminative/generative models. IEEE Trans. Med. Imag. 2008;27(4):495–508. [PMC free article] [PubMed]
[15] Morra J, Tu Z, Apostolova L, Green A, Avedissian C, Madsen S, Parikshak N, Hua N, Toga A, Jack C, Weiner M, Weiner M, Thompson P, the ADNI Validation of a fully automated 3D hippocampal segmentation method using subjects with Alzheimer's disease, mild cognitive impairment, and elderly controls. NeuroImage. 43(1):59–68. [PMC free article] [PubMed]
[16] Raya S, Udupa J. Shape-based interpolation of multidimensional objects. IEEE Trans. Med. Imag. 1990;9(1):32–42. [PubMed]
[17] Herman GT, Zheng J, Bucholtz CA. Shape-based interpolation. IEEE Computer Graphics and Applications. 1992;12(3):69–79.
[18] Treece GM, Prager RW, Gee AH, Berman L. Surface interpolation from sparse cross-sections using region correspondence. IEEE Trans. Med. Imag. 2000;19(11):1106–1114. [PubMed]
[19] Grady L. Random walks for image segmentation. IEEE Trans. Pattern Anal. Machine Intell. 2006;28(11):1768–1783. [PubMed]
[20] Bai X, Sapiro G. A geodesic framework for fast interactive image and video segmentation and matting. Proc. ICCV. 2007:1–8.
[21] Mueller S, Weiner M, Thal L, Petersen R, Jack C, Jagust W, Trojanowski J, Toga A, Beckett L. The Alzheimer's disease neuroimaging initiative. Clin. North Am. 2005;15:869–877. xi–xii. [PMC free article] [PubMed]
[22] Serra J. Introduction to mathematical morphology. Comput. Vision Graph. Image Process. 1986;35(3):283–305.
[23] Couprie M, Bertrand G. Topology preserving alternating sequential filter for smoothing two-dimensional and three-dimensional objects. J. Electron. Imaging. 2004;13(4):720–730.
[24] Mangin J-F, Frouin V, Bloch I, Régis J, López-Krahe J. From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations. Journal of Mathematical Imaging and Vision. 1995;5(4):297–318.
[25] Davatzikos C, Bryan RN. Using a deformable surface model to obtain a shape representation of the cortex. IEEE Trans. Med. Imag. 1996;15:785–795. [PubMed]
[26] MacDonald D. Ph.D. dissertation. McGill Univ.; Canada: 1998. A method for identifying geometrically simple surfaces from threee dimensional images.
[27] Han X, Xu C, Prince J. A topology preserving level set method for geometric deformable models. IEEE Trans. Pattern Anal. Machine Intell. 2003 Jun;25(6):755–768.
[28] Ruuth SJ, Wetton BTR. A simple scheme for volume-preserving motion by mean curvature. J. Sci. Comput. 2003;19:1–3. 373–384.
[29] Qiu A, Fennema-Notestine C, Dale AM, Miller MI. Regional shape abnormalities in mild cognitive impairment and alzheimer's disease. NeuroImage. 2009;45(3):656–661. [PMC free article] [PubMed]
[30] Shattuck D, Leahy R. Automated graph-based analysis and correction of cortical volume topology. IEEE Trans. Med. Imag. 2001;20(11):1167–1177. [PubMed]
[31] Han X, Xu C, Braga-Neto U, Prince JL. Topology correction in brain cortex segmentation using a multiscale, graph-based algorithm. IEEE Trans. Med. Imag. 2002;21(2):109–121. [PubMed]
[32] Ségonne F, Grimson WEL, Fischl B. Topological correction of subcortical segmentation. Proc. MICCAI (2) 2003:695–702.
[33] Bazin P-L, Pham DL. Topology correction of segmented medical images using a fast marching algorithm. Computer Methods and Programs in Biomedicine. 2007;88(2):182–190. [PMC free article] [PubMed]
[34] Taubin G. A signal processing approach to fair surface design. Proc. SIGGRAPH. 1995:351–358.
[35] Desbrun M, Meyer M, Schröder P, Barr AH. Implicit fairing of irregular meshes using diffusion and curvature flow. Proc. SIGGRAPH. 1999:317–324.
[36] Styner M, Oguz I, Xu S, Brechbhler C, Pantazis D, Levitt J, Shenton M, Gerig G. Framework for the statistical shape analysis of brain structures using spharm-pdm. Open Science Workshop at MICCAI. The Insight Journal. 2006 [Online]. Available: [PMC free article] [PubMed]
[37] Reuter M, Wolter F, Peinecke N. Laplace-Beltrami spectra as Shape-DNA of surfaces and solids. Computer-Aided Design. 2006;38:342–366.
[38] Qiu A, Bitouk D, Miller MI. Smooth functional and structural maps on the neocortex via orthonormal bases of the Laplace-Beltrami operator. IEEE Trans. Med. Imag. 2006;25(10):1296–1306. [PubMed]
[39] Rustamov RM. Laplace-beltrami eigenfunctions for deformation invariant shape representation. Eurographics Symposium on Geometry Processing (SGP) 2007:225–233.
[40] Niethammer M, Reuter M, Wolter F-E, Bouix S, Koo NPM-S, Shenton M. Global medical shape analysis using the Laplace-Beltrami spectrum. Proc. MICCAI. 2007;1:850–857. [PMC free article] [PubMed]
[41] Vallet B, Lvy B. Spectral geometry processing with manifold harmonics. Computer Graphics Forum. 2008;27(2):251–260.
[42] Shi Y, Lai R, Krishna S, Sicotte N, Dinov I, Toga AW. Anisotropic Laplace-Beltrami eigenmaps: Bridging Reeb graphs and skeletons. Proc. MMBIA. 2008:1–7. [PMC free article] [PubMed]
[43] Shi Y, Lai R, Kern K, Sicotte N, Dinov I, Toga AW. Harmonic surface mapping with Laplace-Beltrami eigenmaps. Proc. MICCAI. 2008;2:147–154. [PMC free article] [PubMed]
[44] Lai R, Shi Y, Dinov I, Chan TF, Toga AW. Laplace-Beltrami nodal counts: a new signature for 3D shape analysis. Proc. ISBI. 2009:694–697. [PMC free article] [PubMed]
[45] Latecki LJ. 3D well-composed pictures. Graph. Models Image Process. 1997;59(3):164–172.
[46] Siqueira M, Latecki LJ, Gallier J. Making 3D binary digital images well-composed. Vision Geometry XIII, Proc. SPIE. 2005;5675(1):150–163.
[47] Jost J. Riemannian Geometry and Geometric Analysis. 3rd ed. Springer; 2001.
[48] Reuter M. Hierarchical shape segmentation and registration via topological features of Laplace-Beltrami eigenfunctions. Int'l Journal of Computer Vision. 2009 in press,doi:10.1007/s11263-009-0278-1.
[49] Silvester PP, Ferrari RL. Finite Elements for Electrical Engineers. Cambridge University Press; 1996.
[50] Lehoucq RB, Sorensen DC, Yang C. ARPACK Users Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. SIAM; 1998.
[51] Demmel JW, Eisenstat SC, Gilbert JR, Li XS, Liu JWH. A supernodal approach to sparse partial pivoting. SIAM J. Matrix Analysis and Applications. 1999;20(3):720–755.
[52] Khairy K, Howard J. Spherical harmonics-based parametric de-convolution of 3D surface images using bending energy minimization. Med. Image. Anal. 2008;12(2):217–227. [PubMed]
[53] Shi Y, Karl WC. A real-time algorithm for the approximation of level-set-based curve evolution. IEEE Trans. Image Processing. 2008;17(5):645–657. [PMC free article] [PubMed]
[54] Merriman B, Bence J, Osher S. Diffusion generated motion by mean curvature. In: Taylor J, editor. Computational Crystal Growers Workshop. American Mathematical Society; Providence, Rhode Island: 1992. pp. 73–83.
[55] Dice L. Measures of the amount of ecologic association between species. Ecology. 1945;26:297–302.
[56] Rex DE, Ma JQ, Toga AW. The LONI Pipeline processing environment. NeuroImage. 2003;19(3):1033–1048. [Online]. Available: [PubMed]
[58] Sethian J. A fast marching level set method for monotonically advancing fronts. Proc. Nat. Acad. Sci. 1996;93(4):1591–1595. [PubMed]
[59] Tsitsiklis JN. Efficient algorithms for globally optimal trajectories. IEEE Trans. Automat. Contr. 1995 Sep;40(9):1528–1538.