Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3785796

Formats

Article sections

- Abstract
- I. Introduction
- II. Enhancing Tissue Classification
- III. Intrinsic Reeb Graphs on Tissue Boundaries
- IV. Unified Analysis of Geometric and Topological Outliers
- V. Sub-voxel Accuracy with Adaptive Interpolation
- VI. Experimental Results
- VII. Conclusions
- References

Authors

Related links

IEEE Trans Med Imaging. Author manuscript; available in PMC 2014 March 1.

Published in final edited form as:

PMCID: PMC3785796

NIHMSID: NIHMS444728

Yonggang Shi, Rongjie Lai, Arthur W. Toga,^{} and Initiative the Alzheimer’s Disease Neuroimaging^{**}

Yonggang Shi, Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095, USA;

The publisher's final edited version of this article is available at IEEE Trans Med Imaging

See other articles in PMC that cite the published article.

In this paper we present a novel system for the automated reconstruction of cortical surfaces from T1-weighted magnetic resonance images. At the core of our system is a unified Reeb analysis framework for the detection and removal of geometric and topological outliers on tissue boundaries. Using intrinsic Reeb analysis, our system can pinpoint the location of spurious branches and topological outliers, and correct them with localized filtering using information from both image intensity distributions and geometric regularity. In this system, we have also developed enhanced tissue classification with Hessian features for improved robustness to image inhomogeneity, and adaptive interpolation to achieve sub-voxel accuracy in reconstructed surfaces. By integrating these novel developments, we have a system that can automatically reconstruct cortical surfaces with improved quality and dramatically reduced computational cost as compared with the popular FreeSurfer software. In our experiments, we demonstrate on 40 simulated MR images and the MR images of 200 subjects from two databases: the Alzheimer’s Disease Neuroimaging Initiative (ADNI) and International Consortium of Brain Mapping (ICBM), the robustness of our method in large scale studies. In comparisons with FreeSurfer, we show that our system is able to generate surfaces that better represent cortical anatomy and produce thickness features with higher statistical power in population studies.

Cortical surface reconstruction from magnetic resonance (MR) images is a critical problem in brain mapping which provides the geometric foundation for measuring cortical morphometry and tissue integrity. While many sophisticated algorithms were developed for its solution [1]–[11], significant challenges remain in improving the accuracy, robustness, and speed of cortical reconstruction. In this work, we develop a novel system for the automated reconstruction of cortical surfaces from T1-weighted MR images based on intrinsic analysis of geometry and topology using the Reeb graph of Laplace-Beltrami eigenfunctions. We demonstrate that our system can robustly reconstruct high quality cortical surfaces on large scale data sets.

The successful reconstruction of cortical surfaces requires the development of a complicated image analysis system that includes preprocessing steps such as inhomogeneity correction [12], nonlinear registration [13]–[16], skull stripping [17]–[19], and tissue classification [20]–[23]. Even though different combinations of preprocessing methods were made in various systems, a popular choice for surface formation is that a smooth white matter (WM) surface, which represents the boundary between WM and gray matter (GM), with correct topology is first reconstructed, and then deformed to find the GM surface that represents the boundary of the GM and cerebrospinal fluid (CSF). Due to limited image resolution, partial failure of preprocessing steps, or anatomical variability across population, however, *geometric* outliers, such as spurious spikes, and *topological* outliers, such as handles and tunnels, occur frequently during the surface reconstruction process. To build high quality surface representations of tissue boundaries, the challenge is the removal of both types of outliers without sacrificing accuracy.

Geometric and topological outliers were handled separately in previous works. To avoid geometric outliers in the reconstructed surface, smoothness regularization was applied *globally* together with data terms in surface evolution [4], [5], [7], [8], [11]. This regularization-based method, however, has two problems. First, it only helps but does not guarantee the occurrence of geometric outliers can be prevented, which could be due to leakage into non-cortical areas not removed by skull-stripping. Second, shrinkage can arise, especially in deep sulcal regions, because the regularization is enforced everywhere on the surface, even at places it is unnecessary. This could affect the accuracy in cortical thickness measurements and lead to decreased power in statistical analysis.

To ensure genus-zero topology in reconstructed surfaces, topology-preserving deformations can be used [1], [3], [4], [24], [25] after topological outliers are removed. There are mainly two different approaches for topological outlier detection and correction. The first approach works in the voxel space and uses graph analysis and morphological operations for topology analysis [22], [26], [27], however cut decisions are made according to the geometric size of outlier branches only. The second approach uses the triangular mesh representation of the surface and detects topological outliers from overlapping triangles after mapping the surface to the sphere [28], [29]. Cut/fill decisions are then made according to image intensity and curvature distributions. The drawback of this approach is that it relies on the spherical mapping for outlier detection, which is computationally expensive and not suitable for complicated surfaces with large number of handles and tunnels.

In this work, we present a novel system for automated reconstruction of cortical surfaces from MR images. At the core of our system is a unified approach for the correction of geometric and topological outliers based on intrinsic geometry. With well-composed boundary deformation driven by evolution speeds derived from the MR image [30], [31], our system builds manifold representations of tissue boundaries and computes the Laplace-Beltrami (LB) eigenfunctions [32]–[43] as Morse functions for Reeb analysis [44]–[49], where we develop a novel algorithm for the efficient construction of Reeb graphs to capture manifold geometry and topology. By analyzing the loops and branches of the intrinsically defined Reeb graph, we develop a unified approach for the detection and correction of geometric and topological outliers. For geometric outliers, intrinsic Reeb analysis naturally leads us to perform localized outlier filtering without causing unintended shrinkage to other parts of the cortex. By computing paired boundary estimates with different topology constraints, our system carefully takes into account information from tissue analysis and geometric regularity for the cut or fill decisions in the correction of topological outliers.

Besides the novel method for geometric and topological outlier correction, there are two other key components in our system. The first is the enhancement of tissue classification with geometric features to help overcome ambiguities in tissue classification. By using features derived from the Hessian of images [50], [51], this improvement enables us to rely on information in the data instead of artificial assumptions to separate touching gyri and alleviate the impact of partial volume effects on misclassification. The second key element is an adaptive interpolation algorithm that achieves sub-voxel accuracy for the final surface. With locally estimated tissue property, we position the surface at sub-voxel locations for the accurate representation of cortical anatomy.

By integrating these novel developments, we build a robust and efficient system for automated cortical reconstruction. An overview of our system is shown in Fig. 1. Given a skull-stripped MR image, we perform enhanced tissue classification and nonlinear registration as described in section II to design initial evolution speeds for the WM and GM surface. Using methods developed in section III and IV, geometric and topological outliers in the WM evolution speed are then corrected with our unified approach, which produces an evolution speed for the estimation of a clean WM boundary. To obtain the GM boundary, we deform the WM boundary outward with the GM evolution speed. Geometric outliers in the GM boundary are removed with intrinsic Reeb analysis. Finally adaptive interpolation developed in section V is applied to both the WM and GM boundary to generate the WM and GM cortical surfaces with sub-voxel accuracy.

The rest of the paper is organized as follows. In section II, we present the tissue enhancement algorithm and derive the evolution speeds for WM and GM boundary estimation. The construction of the manifold representation of tissue boundaries and their Reeb graphs are developed in section III. Using paired boundary estimates, we develop in section IV the unified approach for geometric and topological outlier correction with intrinsic Reeb analysis and localized filtering. In section V, the adaptive interpolation algorithm is developed for sub-voxel accuracy. Experimental results are presented in section VI, where we present detailed comparisons with FreeSurfer [5], [52] and demonstrate that our system can achieve better performance on both simulated and real MR images from two databases. Finally conclusions are made in section VII.

Tissue classification typically uses statistical models to map intensities to three tissue types in the brain: gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF) [20]–[23], but inhomogeneity and partial volume effects usually compromise the effectiveness of such models in correctly identifying thin structures. To enhance fractions of CSF between touching gyri, skeletons of weighted distance transforms from WM were proposed to open up CSF regions [7], but the challenge is that the WM map by itself has inherent uncertainty. In this section, we propose to use features directly derived from the data, which is the eigen-structure of the Hessian, to enhance tissue maps and capture thin structures. Instead of just enhancing the CSF maps, our method applies to thin structures in both WM and CSF maps, which leads to more integral reconstructions of both WM and GM surfaces. Based on the enhanced tissue maps and atlas labels, we then define evolution speeds for boundary estimation.

Let Ψ : *A* → denote the skull-stripped MR image, where *A* is the lattice of voxels that we represent as *A* = {(*i, j, k*) ^{3}|0 *≤ i < I −* 1, 0 *≤ j < J −* 1, 0 *≤ k < K −* 1}. With partial volume models, a tissue classifier typically uses the intensity at each voxel to assign it as either one of three brain tissue types: GM, WM, CSF, or a weighted combination of WM and GM at the WM boundary, and GM and CSF at the GM boundary. Let *p _{GM}, p_{WM}, p_{CSF}* denote the individual tissue maps of GM, WM, and CSF, which we compute in this work with the FAST tool in the FSL software [21]. Tissue classification is a critical step in our system as it provides the basis for our enhancement algorithm. We designed our system with the flexibility of taking inputs from any software that generates tissue maps with partial volume models. Here we choose the FAST tool because of its robustness across various datasets. It is also relatively efficient and can compute results in less than 10 minutes. At each voxel, the individual tissue map represents the fraction of the voxel belonging to the corresponding tissue. Voxels where all tissue maps are zero are labeled as background (BG). Using the tissue maps, we can define a composite tissue map as

$$\mathrm{\Phi}(\mathbf{p})=\{\begin{array}{cc}0& \text{BG}\\ {p}_{\mathit{CSF}}(\mathbf{p})& \text{if}\phantom{\rule{0.16667em}{0ex}}{p}_{GM}(\mathbf{p})=0,{p}_{WM}(\mathbf{p})=0\\ 1+{p}_{GM}(\mathbf{p})& \text{if}\phantom{\rule{0.16667em}{0ex}}{p}_{\mathit{CSF}}(\mathbf{p})>0,{p}_{GM}(\mathbf{p})>0\\ 2& \text{if}\phantom{\rule{0.16667em}{0ex}}{p}_{GM}(\mathbf{p})=1\\ 2+{p}_{WM}(\mathbf{p})& \text{if}\phantom{\rule{0.16667em}{0ex}}{p}_{GM}(\mathbf{p})>0,{p}_{WM}(\mathbf{p})>0\\ 3& \text{if}\phantom{\rule{0.16667em}{0ex}}{p}_{WM}(\mathbf{p})=1.\end{array}$$

(1)

for all **p** = (*i, j, k*) *A*. This composite map provides a simple way of representing the partial volume tissue models. As an example, the composite tissue map of an MR image in Fig. 2 (a), which is from a patient with Alzheimer’s disease, is plotted in Fig. 2 (e).

Tissue enhancement with Hessian features and atlas labels. Using enhanced tissue maps and atlas labels, evolutions speeds are designed to reconstruct the WM (red) and GM (blue) boundaries of the left hemisphere (LH) and right hemisphere as plotted in **...**

Because tissue classification relies on intensity information to assign tissue types, it does not necessarily follow the geometric assumptions of partial volume models. For example, voxels with small factions of WM tissue should be either on the boundary of definite WM voxels or *ridges* of the image. When neighboring gyri touch each other, voxels in between with tiny fractions of CSF may exhibit as *valleys* in image intensities. Even though such ridges and valleys have distinctive geometric characteristics, partial volume modeling with intensity distributions may not be sensitive enough to detect them. To overcome this difficulty, we will compute an enhanced tissue map based on Hessian features of the image, which is well-known for ridge detection in image analysis [50], [51].

At a voxel **p** = (*i, j, k*) in the MR image Ψ, we denote its gradient as Ψ(**p**), and Hessian as:

$$\text{Hessian}(\mathbf{p})=\left[\begin{array}{ccc}{D}_{xx}(\mathbf{p})& {D}_{xy}(\mathbf{p})& {D}_{xz}(\mathbf{p})\\ {D}_{yx}(\mathbf{p})& {D}_{yy}(\mathbf{p})& {D}_{yz}(\mathbf{p})\\ {D}_{zx}(\mathbf{p})& {D}_{zy}(\mathbf{p})& {D}_{zz}(\mathbf{p})\end{array}\right]$$

(2)

where the entries are the second order gradients of Ψ at this voxel. The eigenvalue of Hessian(**p**) with the maximal magnitude is denoted as *λ _{max}*(

$$TE(\mathbf{p})=\frac{{\lambda}_{\mathit{max}}(\mathbf{p})}{\mid <\nabla \mathrm{\Psi}(\mathbf{p}),{v}_{\mathit{max}}(\mathbf{p})>\mid}$$

(3)

where the inner product between Ψ(**p**) and *v _{max}*(

To filter out these false positives, we use a nonlinear registration software ANTS [16] to register the image to the LPBA40 atlas [53] and warp the probabilistic GM atlas, which assigns to each voxel a probability of belonging to the GM, and anatomical labels to the image space. With the help of nonlinear registration, we can take advantage of atlas information to improve the robustness of our tissue enhancement algorithm and guide the design of evolution speeds for surface reconstruction. Here we choose the ANTS tool for nonlinear warping because of its efficiency in obtaining high quality MR image registration [54]. Typically we can obtain very accurate warped labels and GM atlas in around 30 minutes with the ANTS software. Let
${p}_{GM}^{\mathit{Atlas}}$ denote the warped GM atlas in the image space. As an illustration, the GM atlas
${p}_{GM}^{\mathit{Atlas}}$ is shown in Fig. 2(c), and the warped anatomical labels are shown in Fig. 2(d). Using the anatomical labels, we denote *R _{SUB}* = {

Let *R _{CSF}* = {

For partial WM voxels, we detect voxels on the top of ridges to follow the geometry of partial volume models. Let *R _{WM}* = {

Using these results, we can define an enhanced tissue map as

$$\widehat{\mathrm{\Phi}}(\mathbf{p})=\{\begin{array}{ll}\hfill 0\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}\mathbf{p}\in {\widehat{R}}_{\mathit{CSF}}\hfill \\ \hfill 3\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}\mathbf{p}\in {\widehat{R}}_{WM}\hfill \\ \hfill \mathrm{\Phi}(\mathbf{p})\hfill & \text{else}.\hfill \end{array}$$

(4)

For the MR image shown in Fig. 2(a), the enhanced tissue map is shown in Fig. 2(f). Compared with the original tissue map in Fig. 2(e), we can see that it successfully opens up deep sulci in buried gyral regions, and enhances thin WM regions such as the top of the superior frontal gyrus.

Using the enhanced tissue map, we can design evolution speeds for surface generation based on the fast evolution algorithm [31], [56], [57]. By dividing the image into left hemisphere (LH) and right hemisphere (RH) using the warped anatomical labels as shown in Fig. 2(d), we reconstruct the surfaces of each hemisphere separately. In each hemisphere, the evolution speed defined over the grid *A* for the WM boundary is:

$${F}_{WM}(\mathbf{p})=\{\begin{array}{ll}\hfill 1\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}\widehat{\mathrm{\Phi}}(\mathbf{p})>2.5\phantom{\rule{0.16667em}{0ex}}\text{or}\phantom{\rule{0.16667em}{0ex}}\mathbf{p}\in {R}_{\mathit{SUB}}\hfill \\ \hfill -1\hfill & \text{otherwise}\hfill \end{array}$$

(5)

and the GM boundary is:

$${F}_{GM}(\mathbf{p})=\{\begin{array}{ll}\hfill 1\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}\widehat{\mathrm{\Phi}}(\mathbf{p})>1.5\phantom{\rule{0.16667em}{0ex}}\text{or}\phantom{\rule{0.16667em}{0ex}}\mathbf{p}\in {R}_{\mathit{SUB}}\hfill \\ \hfill -1\hfill & \text{otherwise}\hfill \end{array}$$

(6)

The thresholds in (5) and (6) are chosen for the robustness to remaining inhomogeneities and minimization of mesh distortion in reconstructed surfaces. While the correction of inhomogeneity is applied before tissue classification, some remaining inhomogeneities can still exist in various regions of the brain. For example, in regions with hyper-intensities, large blocks of CSF voxels could be included to the interior of the GM boundary if the threshold is selected toward one in (6). On the other hand, in regions with relatively low intensities, large areas of GM voxels could be left out if the threshold is chosen toward two in (6). To achieve a balance in both situations, we choose the thresholds in (5) and (6) such that only voxels containing at least 50% of the WM or GM tissue are included in the WM or GM speed. As illustrated in Fig. 2 (g) and (h), the final reconstructed boundaries are optimally determined with the adaptive interpolation method in section V, which shifts the position of voxel boundaries with locally estimated tissue properties. By putting voxels with more than 50% of the WM or GM tissue to the interior of the boundary and voxels with less than 50% of the WM or GM tissue to the exterior of the boundary, we restrict the movement of all boundary points to be around half of the voxel resolution. This avoids large distortion of triangles and leads to improved mesh quality in reconstructed WM and GM surfaces.

In the next two sections, we will develop the unified approach for the removal of both geometric and topological outliers in the evolution speeds. After that, sub-voxel accuracy will be achieved with the adaptive interpolation scheme in section V to generate the final cortical surfaces.

In this section, we develop numerical algorithms for the construction of intrinsic Reeb graphs on tissue boundaries that correspond to the interface of different tissue types. We first build a mesh representation of tissue boundaries using the evolution speeds computed in section II. With the LB eigenfunction as the Morse function, a novel numerical algorithm is then developed to construct Reeb graphs that intrinsically capture the geometric and topological structure of the tissue boundary.

Starting from an initial mask, we run a fast boundary evolution algorithm, which we developed in our previous work [31], [56], [57], following the evolution speed in (5) or (6). The boundary evolution process separates the lattice of voxels *A* into two regions: the object region *A _{o}* and the background region

Because the well-composedness condition [30] is satisfied in our reconstruction [31], the boundary between *A _{o}* and

With the triangular mesh representation of the boundary = { , }, we can compute its LB eigen-system by solving a matrix eigenvalue problem [31]–[33]:

$$Qf=\lambda U\phantom{\rule{0.16667em}{0ex}}f$$

(7)

where *λ* is the eigenvalue, *f* :
→ is the eigenfunction, and the two matrices *Q* and *U* are formed using the finite element method [58]. More specifically, the matrices are defined as:

$$\begin{array}{c}{Q}_{ij}=\{\begin{array}{ll}{\scriptstyle \frac{1}{2}}\sum _{{\mathcal{V}}_{j}\in \mathcal{N}({\mathcal{V}}_{i})}\sum _{{\rfloor T}_{l}\in \mathcal{N}({\mathcal{V}}_{i},{\mathcal{V}}_{j})}cot{\theta}_{l}^{i,j},\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}i=j;\hfill \\ -{\scriptstyle \frac{1}{2}}\sum _{{\mathcal{T}}_{l}\in \mathcal{N}({\mathcal{V}}_{i},{\mathcal{V}}_{j})}cot{\theta}_{l}^{i,j},\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}{\mathcal{V}}_{j}\in \mathcal{N}({\mathcal{V}}_{i});\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}\\ {U}_{ij}=\{\begin{array}{ll}{\scriptstyle \frac{1}{12}}\sum _{{\mathcal{V}}_{j}\in \mathcal{N}({\mathcal{V}}_{i})}\sum _{{\mathcal{T}}_{l}\in \mathcal{N}({\mathcal{V}}_{i},{\mathcal{V}}_{j})}\mathit{Area}({\mathcal{T}}_{l}),\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}i=j;\hfill \\ {\scriptstyle \frac{1}{12}}\sum _{{\mathcal{T}}_{l}\in \mathcal{N}({\mathcal{V}}_{i},{\mathcal{V}}_{j})}\mathit{Area}({\mathcal{T}}_{l}),\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}{\mathcal{V}}_{j}\in \mathcal{N}({\mathcal{V}}_{i});\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}\end{array}$$

where
(
) is the set of vertices in the 1-ring neighborhood of
,
(
,
) is the set of triangles sharing the edge (
,
),
${\theta}_{l}^{i,j}$ is the angle in the triangle
opposite to the edge (
,
), and *Area*(
) is the area of the *l*-th triangle
.

The LB eigen-system is discrete and the eigenvalues can be ordered as 0 = *λ*_{0}
*≤ λ*_{1}
*≤ λ*_{2}
*≤ · · ·*. Correspondingly, the eigen-functions are *f*_{0}*, f*_{1}*, f*_{2}, · · ·. The LB eigen-system is intrinsically defined on manifolds and has the nice property of being isometric invariant. It has been applied successfully for various shape analysis tasks in computer vision and medical imaging [32]–[43]. In particular, we will use the first non-constant LB eigenfunction *f*_{1} here as the Morse function for Reeb graph construction [59], which is the solution of the minimization problem

$$\underset{{\int}_{\mathcal{M}}fd\mathcal{M}=0}{min}{\int}_{\mathcal{M}}{\left|\right|{\nabla}_{\mathcal{M}}f\left|\right|}^{2}d\mathcal{M}$$

(8)

and can be viewed as the smoothest map from the manifold
to the real line . As an illustration, the LB eigenfunction *f*_{1} on a WM boundary surface is plotted in Fig. 5(a). Similar to our previous work in hippocampal modeling with LB eigenfunctions [39], we can see the LB eigenfunction intrinsically models the front-to-posterior trend of the cortical surface.

Given a Morse function *f* on the mesh
, its Reeb graph is defined as follows [44].

Let *f* :
→ . The Reeb graph *R*(*f*) of *f* is the quotient space with its topology defined through the equivalent relation *x* *y* if *f*(*x*) = *f*(*y*) for ∀*x*, *y*
.

Following this abstract definition, the Reeb graph is intuitively a graph of level contours of *f* on
. For the construction of Reeb graphs on surfaces with general topology, previous methods typically need to scan through the whole mesh to detect topological changes of level contours [46]–[48]. Because level contours change topology only at critical points of *f*, Reeb graph is essentially a graph of critical points. By using level contours at saddle points, an efficient Reeb graph construction algorithm was proposed in [49], but it uses triangle-based region growing and cannot handle densely distributed saddle points. The partition generated by the method in [49] produces non-manifold regions which are not suitable for further analysis with intrinsic geometry. To overcome these drawbacks, we develop a novel method for Reeb graph construction on high genus meshes. The key idea here is that we augment the mesh with level contours crossing the upper and lower neighborhood of saddle points to process close saddle points. In our method, the surface partition generated by the Reeb graph ensures each surface patch is a manifold, so that further analysis of geometry and topology, such as the computation of geodesics, can be performed.

The critical points of *f* can be classified into maximum, minimum, and saddle points. For a vertex
, and its one-ring neighborhood
(
). Let
(
) = {
(
)|*f* (
) *< f* (
)} and
(
) = {
(
)|*f* (
) *> f*(
)} denote the lower and upper neighbors in the 1-ring neighborhood of a vertex
. Let # denote the number of connected components in a set. Using the number of connected components in
(
) and
(
), we can classify the vertices as follows:

$$\mathit{CLASS}({\mathcal{V}}_{i})=\{\begin{array}{ll}\hfill max\hfill & \#{\mathcal{N}}_{-}({\mathcal{V}}_{i})=1,\#{\mathcal{N}}_{+}({\mathcal{V}}_{i})=0\hfill \\ \hfill min\hfill & \#{\mathcal{N}}_{-}({\mathcal{V}}_{i})=0,\#{\mathcal{N}}_{+}({\mathcal{V}}_{i})=1\hfill \\ \hfill \text{saddle}\hfill & \#{\mathcal{N}}_{-}({\mathcal{V}}_{i})\ge 2,\#{\mathcal{N}}_{+}({\mathcal{V}}_{i})\ge 2\hfill \\ \hfill \text{regular}\hfill & \text{otherwise}.\hfill \end{array}$$

(9)

Let *C* = {*C*_{1}, *C*_{2}, · · ·, *C _{N}*} be the set of critical points of

To construct the Reeb graph, we need to find the arcs connecting these critical points. Let *j* denote the current arc label that is initialized as −1. The augmented set of triangles
and vertices
are initialized as the original triangle set
and vertex set
, respectively. We also assign a label function *ArcLabel* and initialize it to be −1 for all vertices in
. We scan through these critical points sequentially as follows to build the Reeb graph.

If *C _{i}* is a minimum, we increase the current arc label

If *C _{i}* is a saddle point, we first define the isovalues of level contours for arcs entering and leaving this node. For the incoming arcs, the isovalue is defined as:

$${f}_{\mathit{iso}}^{-}({C}_{i})=\frac{f({C}_{i})+max(f({C}_{i-1}),{max}_{{\mathcal{V}}_{j}\in {\mathcal{N}}_{-}({C}_{i})}f({\mathcal{V}}_{j}))}{2}.$$

(10)

For outgoing arcs, the isovalue is defined as:

$${f}_{\mathit{iso}}^{+}({C}_{i})=\frac{f({C}_{i})+min(f({C}_{i+1}),{min}_{{\mathcal{V}}_{j}\in {\mathcal{N}}_{+}({C}_{i})}f({\mathcal{V}}_{j}))}{2}.$$

(11)

The definition ensures that there is no interference from other critical points and that every arc in the Reeb graph would have a manifold structure. This enables us to use existing algorithms for surface analysis such as finding geodesics [60].

For each component in the lower neighborhood
(*C _{i}*) of the saddle point

Mesh augmentation by splitting triangles. Original edge: black. Level contour: red. Inserted edge: green. (a) Keep the triangle. (b) Split into two triangles. (c) Split into three triangles.

For each component in the upper neighborhood
(*C _{i}*) of the saddle point

If *C _{i}* is a maximum, we grow backward with it as the starting point using the algorithm in Table. II to complete an arc

By repeating the above process for all critical points, we complete the construction of the Reeb graph *R*(*f*). While this Reeb graph is defined on the augmented mesh
, which has non-uniform triangles generated from the triangle-splitting process, we can build a map *μ*:
→
to relate the properties derived from *R*(*f*) to the original mesh
that has regular triangles. For each triangle
, the map *μ*(
) denotes the triangle in
such that
is a subset. This is a many-to-one map since the triangles in
is obtained by splitting triangles in
. By establishing a map between the triangles of
and
, we build a bridge between the Reeb graph on the augmented surface
and the evolution speed defined in the voxel space, which we want to modify with intrinsic Reeb analysis. Once the outliers detected by Reeb analysis are corrected, we can generate the final WM and GM surfaces with the *regular* mesh structure of
to produce high quality mesh representations.

As an example, we show in Fig. 4 the construction of the Reeb graph on a double torus. The Morse function, which is the first non-constant LB eigen-function of the double torus, is plotted in Fig. 4 (a), where the level contours at a saddle point *C _{i}* is illustrated. One level contour is generated that crosses edges connecting the lower neighborhood
(

Similar to previous works [5]–[7], our system first reconstructs a clean WM surface with the correct topology and then deform it to obtain the GM boundary. In this section, we develop the unified approach for the removal of geometric and topological outliers on the WM boundary. With topology-preserving evolution, the cleaned WM boundary with genus-zero topology is deformed to generate the GM boundary. While we develop this unified approach in the context of WM surface reconstruction, the method is general and applicable to the analysis of surfaces of arbitrary topology. In particular, we apply it to the GM boundary with spherical topology to remove geometric outliers.

To remove outliers on the WM boundary, we iteratively modify the WM evolution speed with intrinsic Reeb analysis. Using the WM evolution speed, we first compute paired estimates of the WM boundary with and without topology constraints. By comparing the paired boundary estimates, we derive filling voxels for topological artifacts, which provide the basis for cut or fill decisions in topological analysis if the operation is consistent with underlying tissue properties. By analyzing the intrinsically defined Reeb graph on the WM boundary, we detect geometric and topological outliers and design modifications to the WM evolution speed for their removal. With the modified evolution speed, the above steps are repeated until no change is necessary. Next we describe the details of each step.

Given the WM evolution speed *F _{WM}*, we use the fast evolution algorithm described in [31] to compute a pair of estimates of the WM boundary for topological analysis. Starting from a genus-zero mask, we evolve the mask toward the boundary under the topology-preserving constraint, which is ensured by only updating points satisfying the simple and well-composedness condition. This generates a genus-zero estimate of the object boundary. The region enclosed by the boundary is denoted as
${A}_{o}^{G}$ and the background region is denoted as
${A}_{b}^{G}$. The triangular mesh representation of the boundary is
= (
,
). For topological analysis, we turn off the topology-preserving constraint and continue the boundary evolution process under the same speed

By comparing the paired boundary estimates
and
, we can locate paired patches on
that together fill a handle or tunnel in
. For all faces on
with their interior voxel in the background region of
, i.e.,
${\mathbf{x}}_{o}^{T}({\mathcal{T}}_{i}^{G})\in {A}_{b}^{F}$, where the map
${\mathbf{x}}_{0}^{T}$ was defined in section III-A, we perform a connected component labeling and denote the set of connected components as
${CC}^{G}=\{{CC}_{1}^{G},{CC}_{2}^{G},\cdots \}$. For faces in the *k*-th component
${CC}_{k}^{G}\subset {\mathcal{T}}^{G}$, we define a connected component (CC) label
$L({\mathcal{T}}_{i}^{G})=k$. To find paired patches on
that jointly fill a tunnel or handle, we want to find the set of *filling voxels* that connect different components in *CC ^{G}* and fill topological handles or tunnels. For every voxel

Topological analysis of handles and tunnels detected with Reeb analysis. (a) A handle detected with Reeb analysis. (b) Red: arc of the Reeb graph. Green: filling voxels. (c) A tunnel detected by Reeb analysis. (d) Red: arc of the Reeb graph. Green: filling **...**

For the unified analysis of geometric and topological outliers on the WM boundary
, we compute its LB eigenfunctions and construct the Reeb graph. Let *f* = *f*_{1} be the first LB eigenfunction of
, as illustrated for a WM boundary in Fig. 5 (a), and *R*(*f*) = (*C, E*) the Reeb graph of *f* on
, where *C* denotes the set of critical points, and *E* the set of arcs. The augmented mesh from the Reeb graph construction is denoted as * ^{F}* = (
,
). The map from the triangles on
to
is denoted as

We detect geometric outliers from the set of leaf arcs in the Reeb graph. Left *E _{i}* be a leaf arc with
${LC}_{i}^{s}$ and
${LC}_{i}^{e}$ as the start and end level contour. Using the level contour
${LC}_{i}^{s}$ and
${LC}_{i}^{e}$, we define an arc feature

$$H({E}_{i})={\int}_{{LC}_{i}^{s}}<\mathbf{p}(s)-{\overline{\mathbf{q}}}_{i}^{s},\mathbf{n}(\mathbf{p}(s))>ds+{\int}_{{LC}_{i}^{e}}<\mathbf{p}(s)-{\overline{\mathbf{q}}}_{i}^{e},\mathbf{n}(\mathbf{p}(s))>ds$$

(12)

where **n** is the outward normal on the surface
,
${\mathbf{q}}_{i}^{s}$ and
${\mathbf{q}}_{i}^{e}$ are the mean coordinates of points on
${LC}_{i}^{s}$ and
${LC}_{i}^{e}$. If *H*(*E _{i}*)

We consider this leaf arc as an outlier if it satisfies two conditions:

$$\begin{array}{l}\mathit{length}({LC}_{i}^{s})+\mathit{length}({LC}_{i}^{e})<{\rho}_{1}\\ \frac{\mathit{Area}({ET}_{i})}{\mathit{length}({LC}_{i}^{s})+\mathit{length}({LC}_{i}^{e})}>{\rho}_{2}\end{array}$$

(13)

where *Area*(*ET _{i}*) is the sum of area of all triangles in

$$\mathit{ADF}({\mathcal{T}}_{i}^{F})=\frac{\mathit{Area}({\mathcal{T}}_{i}^{F})}{\mathit{Area}({\widehat{\mathcal{T}}}_{i}^{F})}$$

(14)

where
$\mathit{Area}({\mathcal{T}}_{i}^{F})$ and
$\mathit{Area}({\widehat{\mathcal{T}}}_{i}^{F})$ are the area of the triangle
${\mathcal{T}}_{i}^{F}$ in the original mesh
and the projected mesh
. For an outlier leaf arc, we define the set of outlier triangles as *Outlier*(*ET _{i}*) = {

To remove the outlier, we modify the evolution speed *F _{WM}* as follows. If

For topological analysis, we first remove duplicated arcs in *R*(*f*). For any two arcs with the same start and end node, we remove the one with smaller size from *R*(*f*) and add it as a topological outlier to a set we denote as *TO*. After that, we can represent the Reeb graph *R*(*f*) as a matrix *W* and apply standard graph search algorithms to locate the remaining topological outliers. For any arc *E _{i}*

We detect topological outliers from redundant paths between saddle points in the Reeb graph. Using the matrix representation of the Reeb graph, we analyze paths starting at each node to locate topological outliers. For a node *C _{i}* in the directed graph

To provide further localization of topological operations, we compute cutting paths on the handles and tunnels detected by Reeb analysis. For each arc in the outlier set *TO*, we first decide if this is a handle or tunnel on the boundary and then apply different analysis strategies to calculate a cutting path. For better illustration, we show in Fig. 6 the two possible cases that a topological outlier can be represented in the Reeb graph. In Fig. 6(a), the arc of the Reeb graph, which is composed of the triangles plotted in red, captures this outlier attached to a box as a handle. On the other hand, this outlier could also be captured by an arc that forms a tunnel as plotted in Fig. 6(b). To differentiate these two cases, we compute the arc feature *H*(*E _{i}*) as defined in (12) for each arc

A Reeb graph can capture a topological outlier with an arc in the form of a handle (a) or tunnel (b).

We use information from tissue classification and geometric regularity to make cut/fill decisions about topological outliers. Based on tissue maps, we enforce the anatomical knowledge that gyri on WM surfaces should not enclose CSF. Using the area distortion during eigen-projection, we apply the assumption that the WM surfaces should have similar geometric regularity as a folded sheet. From the analysis of the paired boundary estimates and , we have paired patches ${PC}_{1}^{G},{PC}_{2}^{G},\cdots $ on the genus-zero boundary estimate . These patches fill the handles and tunnels on . For each paired patches, its interior voxels are ${\mathbf{x}}_{o}^{T}({PC}_{i}^{G})$ and they are filled inside ${A}_{o}^{G}$ to satisfy the genus-zero constraint, which are illustrated in Fig. 7 for a handle and tunnel. At a paired patch ${PC}_{i}^{G}$, a cut decision should be made if either of the following two conditions is met

- If the number of voxels $\mathbf{p}\in {\mathbf{x}}_{o}^{T}({PC}_{i}^{G})$ with (
**p**) ≤ 1 is greater than*THD*._{CSF} - If the number of triangles ${\mathcal{T}}_{j}\in {PC}_{i}^{G}$ with
*ADF*( )*> γ*is greater than*THD*._{GEO}

Both parameters *THD _{CSF}* and

To implement the cut/fill decision from the analysis of paired patches, we use the handles and tunnels detected by the Reeb analysis process. To cut open a paired component in
, we pick a cutting path on a handle or tunnel that is connected to this paired component. For the cutting path *CP _{j}* of an arc

$${F}_{WM}(\mathbf{p})=-1\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\forall \mathbf{p}\in {\mathbf{x}}_{o}^{T}({\mathcal{T}}^{F}(C{P}_{{J}_{\mathit{min}}})).$$

(15)

Similar to the modification for the removal of geometric outliers, this change to the evolution speed is *local* and peels off a layer of voxels along the cutting path. To completely cut open a large outlier, this process might be applied to the same outlier during each iteration of our unified outlier correction algorithm. The number of iterations it takes for the topology correction algorithm to converge depends on the complexity of the topological outliers in the mask boundary. In our experience from experiments on various MR datasets, this process typically takes less than five iterations to converge.

With the modified speed obtained from the geometric and topological analysis process in section IV-B and IV-C, we can repeat the above steps in section IV-A, IV-B and IV-C until no changes are made to the WM speed. The genus zero estimate then gives a clean estimate of the WM boundary. For the example in Fig. 2(a), the cleaned WM boundary after the removal of outliers shown in Fig. 5 is plotted in Fig. 8(a).

The WM and GM surface before and after sub-voxel calculations. (a) WM boundary after removal of outliers. (b) WM surface with sub-voxel accuracy. (c) GM boundary after removal of outliers. (d) GM surface with sub-voxel accuracy.

In order to reconstruct the GM boundary, we evolve the cleaned WM boundary outward using the GM speed defined in (6) under the topology-preserving constraint. To remove geometric outliers on the GM boundary estimate, the Reeb analysis is applied to this genus zero surface and the geometric outlier correction process described in section IV-B is applied iteratively. By removing the spurious outliers, we obtain a clean estimate of the GM boundary. As an illustration, the cleaned GM boundary corresponding to the example in Fig. 2 is plotted in Fig. 8(c).

After removing geometric and topological outliers, we obtain a clean reconstruction of the boundary between two tissue types. For the WM surface, it is the boundary between WM and GM. For the GM surface, it is the boundary between GM and CSF. To achieve sub-voxel accuracy in the final surface reconstruction, we develop an adaptive interpolation method in this section.

For a tissue boundary between a high intensity tissue Φ* _{h}* and low intensity tissue Φ

For each voxel **p** *BV _{in}*

$$\alpha (\mathbf{p})=\{\begin{array}{ll}\hfill 1\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}\mathrm{\Psi}(\mathbf{p})\ge {\mathrm{\Psi}}_{h}(\mathbf{p});\hfill \\ \hfill 0\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}\mathrm{\Psi}(\mathbf{p})\le {\mathrm{\Psi}}_{l}(\mathbf{p});\hfill \\ \hfill {\scriptstyle \frac{\mathrm{\Psi}(\mathbf{p})-{\mathrm{\Psi}}_{l}(\mathbf{p})}{{\mathrm{\Psi}}_{h}(\mathbf{p})-{\mathrm{\Psi}}_{l}(\mathbf{p})}}\hfill & \text{otherwise}.\hfill \end{array}$$

(16)

To account for this partial volume effect, we will shift faces on the boundary. Let *β*(**p**) = min(1, *N _{x}*(

For a face *B _{i}*, if the tissue map of its interior voxel
$\mathrm{\Phi}({\mathbf{x}}_{o}^{B}({B}_{i}))<{\mathrm{\Phi}}_{h}$, it needs to move inward according to the shift factor

$${S}_{B}(i)=\{\begin{array}{cc}(-{\scriptstyle \frac{{\alpha}^{\prime}({\mathbf{p}}^{\ast})}{{N}_{x}({\mathbf{p}}^{\ast})}}{h}_{x}\phantom{\rule{0.16667em}{0ex}}0\phantom{\rule{0.16667em}{0ex}}0)& \text{if}\phantom{\rule{0.16667em}{0ex}}{\mathbf{p}}_{x}^{\ast}-{\widehat{\mathbf{p}}}_{x}^{\ast}<0\\ ({\scriptstyle \frac{{\alpha}^{\prime}({\mathbf{p}}^{\ast})}{{N}_{x}({\mathbf{p}}^{\ast})}}{h}_{x}\phantom{\rule{0.16667em}{0ex}}0\phantom{\rule{0.16667em}{0ex}}0)& \text{if}\phantom{\rule{0.16667em}{0ex}}{\mathbf{p}}_{x}^{\ast}-{\widehat{\mathbf{p}}}_{x}^{\ast}>0\\ (0\phantom{\rule{0.16667em}{0ex}}-{\scriptstyle \frac{{\alpha}^{\prime}({\mathbf{p}}^{\ast})}{{N}_{y}({\mathbf{p}}^{\ast})}}{h}_{y}\phantom{\rule{0.16667em}{0ex}}0)& \text{if}\phantom{\rule{0.16667em}{0ex}}{\mathbf{p}}_{y}^{\ast}-{\widehat{\mathbf{p}}}_{y}^{\ast}<0\\ (0\phantom{\rule{0.16667em}{0ex}}\frac{{\alpha}^{\prime}({\mathbf{p}}^{\ast})}{{N}_{y}({\mathbf{p}}^{\ast})}{h}_{y}\phantom{\rule{0.16667em}{0ex}}0)& \text{if}\phantom{\rule{0.16667em}{0ex}}{\mathbf{p}}_{y}^{\ast}-{\widehat{\mathbf{p}}}_{y}^{\ast}>0\\ (0\phantom{\rule{0.16667em}{0ex}}0-{\scriptstyle \frac{{\alpha}^{\prime}({\mathbf{p}}^{\ast})}{{N}_{z}({\mathbf{p}}^{\ast})}}{h}_{z})& \text{if}\phantom{\rule{0.16667em}{0ex}}{\mathbf{p}}_{z}^{\ast}-{\widehat{\mathbf{p}}}_{z}^{\ast}<0\\ (0\phantom{\rule{0.16667em}{0ex}}0\phantom{\rule{0.16667em}{0ex}}{\scriptstyle \frac{{\alpha}^{\prime}({\mathbf{p}}^{\ast})}{{N}_{z}({\mathbf{p}}^{\ast})}}{h}_{z})& \text{if}\phantom{\rule{0.16667em}{0ex}}{\mathbf{p}}_{z}^{\ast}-{\widehat{\mathbf{p}}}_{z}^{\ast}>0\end{array}$$

(17)

where the shift in each direction is multiplied by the corresponding spatial resolution, and divided by the number of faces that
(**p**^{*}) ∩ *B* has in that direction.

Given the shift of all faces in *B*, we calculate the shift for a vertex
as the average shift of its neighboring faces:

$${S}_{\mathcal{V}}(i)=\frac{{\sum}_{{B}_{j}\in {\mathcal{N}}_{R}({\mathcal{V}}_{i})}{S}_{B}(j)}{\#{\mathcal{N}}_{R}({\mathcal{V}}_{i})}$$

(18)

where ( ) denote the set of rectangular faces in the neighborhood of the vertex , and # ( ) is the number of rectangular faces in this set.

Let **x** denote the vector of coordinates of all vertices, and
the shift of all vertices. We compute the final vertex coordinates with sub-voxel accuracy by minimizing the following energy function:

$$\mathcal{E}(\widehat{\mathbf{x}})={\left|\right|\mathbf{x}+{S}_{\mathcal{V}}-\widehat{\mathbf{x}}\left|\right|}^{2}+\eta {\left|\right|\mathrm{\Delta}\widehat{\mathbf{x}}\left|\right|}^{2}$$

(19)

where **$\widehat{x}$** is the optimized vertex coordinate, Δ is the discrete Laplacian matrix of the mesh, and *η* is a regularization parameter. Note that we have chosen the Laplacian instead of the gradient operator in the regularization term to avoid large shrinkage effects from the gradient operator. The solution of this quadratic problem gives us the coordinates of vertices on the smoothed surface:

$$\widehat{\mathbf{x}}={(I+\eta {\mathrm{\Delta}}^{\prime}\mathrm{\Delta})}^{-1}(\mathbf{x}+{S}_{\mathcal{V}}).$$

(20)

By applying the sub-voxel calculation to both the WM and GM surfaces, we obtain our final solution for cortical surface reconstruction. As an illustration, we have plotted the surfaces before and after sub-voxel interpolation in Fig. 8. Clearly we can see the improved quality in the reconstructed surfaces after the application of the adaptive interpolation scheme.

In this section, we present experimental results on 40 simulated MR images and T1-weighted real MR images of 200 subjects to demonstrate that our method can reconstruct high quality cortical surfaces on large scale data sets. We also compare our results with surfaces reconstructed with FreeSurfer [5], [52], a widely used software for cortical surface reconstruction and analysis in brain imaging research, to demonstrate that our method can obtain better results with significantly less computational cost.

The simulated dataset includes a set of 20 simulated MR images from a publicly available dataset [61], and a set of 20 images we simulated with known gray matter atrophy. The real T1-weighted MR images used in our experiments are from two publicly available databases. The first dataset has 100 baseline MR images, from 50 normal controls (NC) and 50 Alzheimer’s disease (AD) patients, from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database [62]^{1}. The second dataset includes MR images of 100 subjects from the International Consortium of Brain Mapping (ICBM) [63]. The 100 subjects have a wide age range from 19 to 80.

All MR images were first automatically skull-stripped with a meta algorithm [64] in the LONI pipeline [65]. Using the skull-stripped images, our system and FreeSurfer automatically reconstructed the WM and GM surfaces on the left and right hemispheres in all subjects. For our method, the same set of parameters are used: *ρ*_{1} = 100*mm, ρ*_{2} = 5*, γ* = 100, *THD _{CSF}* = 5,

As an illustration, we first present in this experiment a detailed comparison between our method and FreeSurfer using the same MR scan plotted in Fig. 2. This is an MR image of an AD patient from the ADNI dataset. Representative examples from both the ADNI and ICBM datasets will then be presented to demonstrate that our method can reconstruct surfaces that better preserve the integrity of cortical anatomy.

For the MR image of an AD patient, we applied our method and FreeSurfer to reconstruct the WM and GM surfaces on both hemispheres. Computationally our method took around 3 hours and is much more efficient than FreeSurfer, which took over 10 hours. Both our method and FreeSurfer successfully reconstructed cortical surfaces with genus-zero topology. The unified outlier removal process of our method took around 20 minutes for each hemisphere. The topological correction process in FreeSurfer took more than one and half hour for each hemisphere. As a more detailed illustration of the Reeb analysis process for outlier removal, we plotted in Fig. 9 the number of nodes in the Reeb graph, the number of geometric outliers, and the number of total handles and tunnels at each iteration of our algorithm. With the increase of iteration, we can see the complexity of the Reeb graph decreases with the number of outliers in the boundary. For this image, it took 3 iterations for the unified outlier correction algorithm to converge when all geometric outliers were removed and no more cuts to be made. After that, a topology-preserving evolution was applied to reconstruct a genus-zero boundary. The adaptive interpolation scheme developed in section V was finally applied to reconstruct a smooth and accurate representation of the cortical surface. As an illustration, we plotted the intersection of the image slice shown in Fig. 2(a) with the WM and GM surfaces reconstructed by our method in Fig. 2(g) and (h). From the results we can see that our method is able to robustly reconstruct the WM and GM boundaries even in the presence of tissue degeneration such as the white matter degeneration on the right hemisphere of this AD patient.

To compare the results from our method and FreeSurfer, the left hemispherical WM surfaces reconstructed by both methods are plotted in Fig. 10. The two WM surfaces have similar number of vertices because they are both derived from the WM mask boundary. The FreeSurfer WM surface has 137309 vertices, and our WM surface has 138130 vertices. In regions highlighted in the dashed circles in Fig. 10, we can see that our method produces more complete reconstruction of the boundary, which is better illustrated with the intersections of surfaces and two axial slices shown in Fig. 10(c) and (f). The left hemispherical GM surfaces reconstructed with our method and FreeSurfer are plotted in Fig. 11. Because the FreeSurfer GM surface is obtained by deforming the WM mesh to the GM boundary, it has the same number of vertices as its WM surface and typically has irregular triangles. On the contrary, our method follows the boundary between GM and CSF to generate a more uniform mesh representation of the GM surface. The GM surface reconstructed with our method here has 145136 vertices, which has slightly more vertices than the WM surface. This is natural in order for the GM surface to have the same level of vertex density and accuracy as the WM surface because it has a larger area than the WM surface. From the highlighted mesh structures of four ROIs in Fig. 11, we can see that FreeSurfer produced *broken* gyri in these regions while our method is able to generate more complete reconstruction at these gyri. The intersections of the GM surfaces reconstructed by our method and FreeSurfer with three sagittal slices are shown in Fig. 12, which illustrate that our surface can better capture deep sulcal regions. This shows that the localized outlier detection and filtering approach in our method is able to avoid shrinkage and produce a more accurate surface representation in these regions.

Comparison of WM surfaces reconstructed by our method and FreeSurfer. Surfaces are colored by their mean curvature and plotted in both lateral and medial views. (a)(d) Our result. (b)(e) FreeSurfer result. (c)(f) Intersection of an image slice with circled **...**

Comparison of the GM surface reconstructed by our method and FreeSurfer. (a)(b) Lateral and medial views of our result. (c)(d) Lateral and medial views of FreeSurfer result. The mesh structures of four highlighted ROIs are plotted to demonstrate differences **...**

Intersection of image slices with GM surfaces reconstructed with our method (red contours) and FreeSurfer (blue contours). Green arrows point to locations where our surface can better capture deep sulcal regions.

To further demonstrate that surfaces generated by our method are able to better represent the integral anatomy of the cortex, we have plotted in Fig. 13 more representative reconstruction results from ADNI and ICBM data. For each subject, the left hemispherical GM surfaces reconstructed by our method and FreeSurfer are plotted with differences highlighted in dashed circles. Compared with results from our method, we can see that FreeSurfer results failed to reconstruct a complete gyrus in the regions enclosed in the dashed circles. This shows that our method is able to produce a more complete representation of cortical anatomy.

In this experiment, we will compare our method and FreeSurfer on the detection of simulated atrophy in MR images where the ground truth of WM, GM, and CSF tissue maps are known. By matching reconstructed cortical surfaces with the underlying tissue boundaries, we can analyze the performance of different topology correction strategies in our method and FreeSurfer. Using the reconstructed WM and GM surfaces, we can calculate cortical thickness and compare quantitatively the ability of these two methods in detecting *sub-voxel* tissue atrophy.

To simulate longitudinal brain atrophy, we first downloaded from BrainWeb [66] a set of 20 simulated MR images of normal subjects with known WM, GM, and CSF tissue maps [61] and used them as the baseline images. All images have an isotropic spatial resolution of 1mm. Using these tissue maps, we simulated a 0.1 mm GM atrophy everywhere on the cortex by modifying intensities in the MR images. For each voxel with more than 10 percent GM tissue and on the boundary of GM and CSF, we introduced a sub-voxel atrophy by subtracting its intensity by 10 percent of the difference between the average GM intensity and CSF intensity in the image because we assumed the lost GM tissue in this voxel will be replaced by CSF. By applying this procedure to all baseline images, we obtained a set of 20 images with simulated GM atrophy, which we denoted as follow-up scans of these 20 subjects. Both our method and FreeSurfer were applied to the 40 simulated images to reconstruct the WM and GM cortical surfaces.

From the examples shown in Fig. 10 and Fig. 13, we see various cases that our method is able to generate a more integral reconstruction of the superior temporal cortex than FreeSurfer. With simulated images, we have the opportunity to demonstrate this on images with known tissue boundaries. In Fig. 14(a) and (d), we plotted the left WM surfaces reconstructed from one of the baseline images with our method and FreeSurfer, where we can clearly see that our result has a more integral representation of the superior temporal cortex. For both methods, we extracted the WM boundary before topology correction to investigate the cause of this difference. In Fig. 14 (b) and (d), we plotted the corresponding WM boundary without topology correction in the region highlighted by the dashed ellipses in Fig. 14 (a) and (d), respectively. In both cases, we can see a hole is clearly present on the superior temporal gyrus. Different correction strategies were adopted by these two methods: our method chose to fill the hole, while FreeSurfer decided to cut it open. As a result, different surface reconstruction results were obtained as shown in Fig. 14(c) and (f). The intersection of the reconstructed surfaces with two sagittal slices of the underlying true WM tissue map were plotted in Fig. 14 (g) and (h), where our result was plotted in red and the FreeSurfer result was plotted in blue. As highlighted by the green arrows in these two pictures, we can see that our result provides a more faithful representation of the underlying tissue boundary at the superior temporal cortex.

Comparison of our method and FreeSurfer on a simulated image with known WM tissue map. (a) The WM surface reconstructed by our method. (b)(c) A zoomed view of the circled region in (a) before topology correction and after final reconstruction. (d) The **...**

To compare the performance of our method and FreeSurfer in the detection of simulated longitudinal atrophy, we parcellated all cortical surfaces into gyral regions and tested for GM thickness changes in each region. In this work, we used the gyral labels on GM surfaces generated by FreeSurfer for a localized comparison. In order to project the FreeSurfer label onto our surface, we aligned both surfaces in the image space. For each vertex on our surface, we found the nearest vertex on the FreeSurfer surface such that the normal directions of both vertices have positive inner products and pulled back the label on this vertex. As an illustration, we plotted in Fig. 15 the labels on the GM surfaces generated by FreeSurfer and our method as shown in Fig. 11. We can see that we have successfully mapped the gyral labels onto the GM surface generated by our method. At each point of a GM surface, we used its distance to the closest point on the WM surface as the measure of thickness [7], [67]. Numerically, we first computed the signed distance function of the WM surface with the fast marching method [68], and then calculated the thickness at each vertex of the GM surface with sub-millimeter accuracy by interpolating the value of the signed distance function. For surfaces generated by both our method and FreeSurfer, the same approach of thickness calculation was applied for consistency. On results generated by both our method and FreeSurfer, we calculated the change of mean thickness on each gyrus between the follow-up and baseline scan of each subject. On each gyrus, a two-tailed t-test was applied to the thickness changes of the 20 subjects and tested if we could reject the null-hypothesis that there is no significant change of thickness between the baseline and follow-up scan. For cortical surfaces generated by our method and FreeSurfer, the p-values obtained from the t-tests at each gyrus were mapped onto the corresponding gyrus of a subject and plotted in Fig. 16, where p-values larger than 0.01 were considered insignificant and plotted in light gray. From the results, it is clear that our method was able to detect significant atrophy in more regions and generate more significant p-values than FreeSurfer. More specifically, on the left hemisphere, our results produced more significant p-values on 20 of the 34 gyral regions. On the right hemisphere, our results produced more significant p-values on 25 of the 34 gyral regions. Next we will demonstrate that our method can achieve more statistically significant results than FreeSurfer on the real data from ADNI and ICBM.

In the third experiment, we applied both our method and FreeSurfer to the 100 MR images from ADNI. To compare the performance of our method and FreeSurfer in detecting group differences, we parcellated cortical surfaces into gyral regions and use the mean thickness of each gyrus as the variable for statistical tests. As described in the previous experiment, we used the gyral labels on GM surfaces generated by FreeSurfer to parcellate our surfaces via nearest neighbor projection.

For surfaces generated by both methods, the thickness of each vertex on the GM surface was computed as its distance to the corresponding WM surface. For each gyrus, we computed the average gray matter thickness and tested group differences between NC and AD using two tailed t-tests. As an illustration of the thickness differences across groups, we plotted in Fig. 17 the mean thickness of gyral regions from the NC and AD group. By comparing the thickness of corresponding gyrus from the two groups, we can clearly see that the AD group have thinner cortex throughout the brain. For the results produced by our method and FreeSurfer, the p-values from t-tests were mapped onto the corresponding gyrus of the left hemisphere (LH) and right hemisphere (RH), and plotted in Fig. 18. Note that p-values larger than 0.01 were considered insignificant and plotted in light gray. We can see that both methods are very effective in detecting differences between the NC and AD groups on both hemispheres, whereas the degrees of difference vary on different gyri. For the left hemisphere, results from our method produces more significant p-values on 33 of the 34 regions. For the right hemisphere, our results generate more significant p-values on 30 of the 34 regions.

Gyrus-based map of p-values from the testing of NC versus AD group differences on LH and RH surfaces from our results (a,b), and FreeSurfer results (c, d).

For the two-tailed t-test on the mean thickness of each gyrus, we also performed power analysis [69] to compare the power of thickness features generated by our method and FreeSurfer in detecting group differences. The power of the t-test is the probability that the null hypothesis will be rejected if it is false. The type I error, i.e., the p-value, was set to 0.01. For each method, the effective size at a gyrus was computed as
${\scriptstyle \frac{{MT}_{NC}-{MT}_{AD}}{\sigma}}$ [69], where *MT _{NC}* and

In the fourth experiment, we tested on MR images of 100 subjects from the ICBM database. We split the 100 subjects into young and elderly groups using an age threshold of 48 years such that each group had 50 subjects. For each subject, we downloaded two scans for our experiments. For the first scan of each subject, we applied our method and FreeSurfer for comparison. The repeat scan of each subject will be used to demonstrate the reproducibility of our method.

Using surfaces generated by both methods on the first scan of each subject, we tested for group differences between the young and elderly groups and compared the statistical power of our method and FreeSurfer. Using the same method described in VI-B, we projected FreeSurfer labels onto the surfaces generated by our method and applied two-tailed t-tests to the mean thickness of each gyrus from the young and elderly groups. The mean thickness of gyral regions obtained from our method were plotted in Fig. 19, where we can see that the elderly group have lower thickness throughout the cortex. The p-value maps generated by our method and FreeSurfer were plotted in Fig. 20, where the same colormap as in Fig. 18 was used and p-values over 0.01 were plotted in light gray. For each hemisphere, t-tests with thickness measures generated by our method have more significant p-values than FreeSurfer on 33 of the 34 gyri.

Gyrus-based map of p-values from the testing of young and elderly group differences on LH and RH surfaces from our results (a, b) and FreeSurfer results (c, d).

Similar to the experiment on ADNI data, we also performed power analysis for the mean thickness measure generated by our method and FreeSurfer at each gyrus. For both methods, the p-value was set to 0.01, and the effective size was computed as
${\scriptstyle \frac{{MT}_{Y}-{MT}_{O}}{\sigma}}$, where *MT _{Y}* and

Power analysis on first scans of ICBM dataset: number of thickness features from our method and FreeSurfer with power above thresholds.

To test the reproducibility of our method, we also applied our method to the repeat scan from each of the 100 subjects in the ICBM dataset. To parcellate each GM surface reconstructed from a repeat scan into gyral regions, we rigidly aligned it with the corresponding GM surface from the first scan and used the nearest neighbor projection as described in VI-B to generate the gyral labels. For the results obtained from the repeat scans, the same power analysis was applied to the mean thickness of each gyrus from the young and elderly group. The results were listed in Table V. Compared with the results from the first scans, we can see that our method is able to successfully generate very similar and large numbers of sensitive thickness features in the repeat scans, which demonstrates the reproducibility of our method.

In this paper we have developed a novel system for automated cortical reconstruction from T1-weighted MR images. We have demonstrated that our system is able to generate high quality surfaces with much reduced computational cost when compared with FreeSurfer. The unified approach presented in this work for the analysis of geometric and topological structures, while developed in the context of cortical reconstruction, is general and applicable for other medical image analysis problems.

In this work, we focused on the reconstruction of the WM and GM surfaces. There are also considerable interests in obtaining a central cortical surface from MR images as a balanced representation of the WM and GM surfaces [7]. There are two different approaches we can follow to achieve this goal. Using the reconstructed WM and GM surfaces, a surface evolution algorithm can be developed as in [7] to evolve the WM surface to the middle of the GM tissue. We can also take a shape analysis approach to find the central cortical surface. In our recent work, we developed an intrinsic surface mapping method that can evolve the WM and GM surface toward each other in a high dimensional embedding space by optimizing their conformal metrics [43]. By establishing one-to-one correspondences between the WM and GM surface with this intrinsic mapping algorithm, we can compute the central cortical surface as the mean shape of the WM and GM surface. Using the one-to-one map between the WM and GM surfaces, we can also develop a more robust method to compute cortical thickness than using the distance transform of the WM surface, which can result in possible underestimation of thickness in the extreme case of neighboring sulcal banks with highly uneven thickness.

In future work, we will also perform more extensive validations on data from a wide range of neuroimaging studies, and extend the system to the consistent reconstruction of surfaces on data from longitudinal studies [11], [71].

Data collection and sharing from ADNI was funded by National Institutes of Health Grant U01 AG024904. ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: Abbott; Alzheimers Association; Alzheimers Drug Discovery Foundation; Amorfix Life Sciences Ltd.; AstraZeneca; Bayer Health-Care; BioClinica, Inc.; Biogen Idec Inc.; Bristol-Myers Squibb Company; Eisai Inc.; Elan Pharmaceuticals Inc.; Eli Lilly and Company; F. Hoffmann-La Roche Ltd and its affiliated company Genentech, Inc.; GE Healthcare; Innogenetics, N.V.; Janssen Alzheimer Immunotherapy Research & Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development LLC.; Medpace, Inc.; Merck & Co., Inc.; Meso Scale Diagnostics, LLC.; Novartis Pharmaceuticals Corporation; Pfizer Inc.; Servier; Synarc Inc.; and Takeda Pharmaceutical Company. The Canadian Institutes of Health Research is providing funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health (www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer’s Disease Cooperative Study at the University of California, San Diego. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of California, Los Angeles. Data collection and sharing from ADNI was also supported by NIH grants P30 AG010129 and K01 AG030514.

This work was in part supported by the National Institute of Health (NIH) under Grant K01EB013633 and Grant 5P41RR013642.

^{1}Part of the data used in the preparation of this article were obtained from the Alzheimers Disease Neuroimaging Initiative (ADNI) database (adni.loni.ucla.edu). ADNI was launched in 2003 by the National Institute on Aging (NIA), the National Institute of Biomedical Imaging and Bioengineering (NIBIB), the Food and Drug Administration (FDA), private pharmaceutical companies and non-profit organizations, as a $60 million, 5- year public-private partnership. The primary goal of ADNI has been to test whether serial magnetic resonance imaging (MRI), positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and early Alzheimers disease (AD). Determination of sensitive and specific markers of very early AD progression is intended to aid researchers and clinicians to develop new treatments and monitor their effectiveness, as well as lessen the time and cost of clinical trials.

Yonggang Shi, Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095, USA.

Rongjie Lai, Department of Mathematics, University of Southern California, Los Angeles, CA 90089, USA.

Arthur W. Toga, Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095, USA.

1. Mangin JF, Frouin V, Bloch I, et al. From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations. J MATH IMAGING VIS. 1995;5(4):297–318.

2. Carman GJ, Drury HA, Van Essen DC. Computational methods for reconstructing and unfolding the cerebral cortex. Ceberal cortex. 1995;5:506–517. [PubMed]

3. 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]

4. MacDonald D. PhD dissertation. McGill Univ; Canada: 1998. A method for identifying geometrically simple surfaces from threee dimensional images.

5. Dale AM, Fischl B, Sereno MI. Cortical surface-based analysis i: segmentation and surface reconstruction. NeuroImage. 1999;9:179–194. [PubMed]

6. Shattuck D, Leahy R. BrainSuite: An automated cortical surface identification tool. Med Image Anal. 2002;8(2):129–142. [PubMed]

7. Han X, Pham DL, Tosun D, et al. CRUISE: Cortical reconstruction using implicit surface evolution. NeuroImage. 2004;23:997–1012. [PubMed]

8. Kim JS, Singh V, Lee JK, Lerch J, Ad-Dab’bagh Y, MacDonald D, Lee JM, Kim SI, Evans AC. Automated 3-D extraction and evaluation of the inner and outer cortical surfaces using a laplacian map and partial volume effect classification. NeuroImage. 2005;27:210–221. [PubMed]

9. Xue H, Srinivasan L, Jiang S, Rutherford M, Edwards AD, Rueckert D, Hajnal JV. Automatic segmentation and reconstruction of the cortex from neonatal MRI. NeuroImage. 2007;38(3):461– 477. [PubMed]

10. Liu T, Nie J, Tarokh A, Guo L, STW Reconstruction of central cortical surface from brain MRI images: method and application. NeuroImage. 2008;40(3):991–1002. [PMC free article] [PubMed]

11. Li G, Nie J, Wu G, Wang Y, Shen D. Consistent reconstruction of cortical surfaces from longitudinal brain mr images. NeuroImage. 2012;59(4):3805– 3820. [PMC free article] [PubMed]

12. Sled J, Zijdenbos A, Evans A. A nonparametric method for automatic correction of intensity nonuniformity in mri data. IEEE Trans Med Imag. 1998;17(1):87–97. [PubMed]

13. Christensen GE, Rabbitt RD, Miller MI. Deformable templates using large deformation kinematics. IEEE Trans Imag Process. 1996;5(10):1435–1447. [PubMed]

14. Joshi S, Miller MI. Landmark matching via large deformation diffeomorphisms. IEEE Trans Imag Process. 2000;9(8):1357–1370. [PubMed]

15. Shen D, Davatzikos C. HAMMER: hierarchical attribute matching mechanism for elastic registration. IEEE Trans Med Imag. 2002;21(11):1421–39. [PubMed]

16. Avants BB, Epstein CL, Grossman M, Gee JC. Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain. Med Image Anal. 2008;12(1):26–42. [PMC free article] [PubMed]

17. Smith S. Fast robust automated brain extraction. Hum Brain Mapp. 2002;17(3):143–55. [PubMed]

18. Ségonne F, Dale A, Busa E, Glessner M, Salat D, Hahn H, Fischl B. Robust brain extraction across datasets and comparison with publicly available methods. NeuroImage. 2004;22(3):1660–75.

19. Iglesias J, Liu CY, Thompson P, Tu Z. Robust brain extraction across datasets and comparison with publicly available methods. IEEE Trans Med Imag. 2011;30(9):1617–34. [PubMed]

20. Leemput KV, Maes F, Vandermeulen D, Suetens P. Automated model-based tissue classification of MR images of the brain. IEEE Trans Med Imag. 1999;18(10):897–908. [PubMed]

21. Zhang Y, Brady M, Smith S. Segmentation of brain MR images through a hidden Markov random field model and the expectation maximization algorithm. IEEE Trans Med Imag. 2001;20(1):45–57. [PubMed]

22. Shattuck D, Leahy R. Automated graph-based analysis and correction of cortical volume topology. IEEE Trans Med Imag. 2001;20(11):1167–1177. [PubMed]

23. Bazin P, Pham D. Topology-preserving tissue classification of magnetic resonance brain images. IEEE Trans Med Imag. 2007;26(4):487–96. [PubMed]

24. 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.

25. Bazin P, Pham D. Topology correction using fast marching methods and its application to brain segmentation. Proc MICCAI. 2:484–91. [PubMed]

26. 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]

27. Chen L, Wagenknecht G. Topological correction of volumetric binary brain segmentation using a multiscale algorithm. Proc ISBI. 2007:1308–11.

28. Ségonne F, Grimson E, Fischl B. A genetic algorithm for the topology correction of cortical surfaces. Proc IPMI. 2005;19:393–405. [PubMed]

29. Yotter R, Dahnke R, Thompson P, Gaser C. Topological correction of brain surface meshes using spherical harmonics. Hum Brain Mapp. 2011;32(7):1109– 24. [PubMed]

30. Latecki LJ. 3D well-composed pictures. Graph Models Image Process. 1997;59(3):164–172.

31. Shi Y, Lai R, Morra J, Dinov I, Thompson P, Toga A. Robust surface reconstruction via Laplace-Beltrami eigen-projection and boundary deformation. IEEE Trans Med Imag. 2010;29(12):2009–2022. [PMC free article] [PubMed]

32. Reuter M, Wolter F, Peinecke N. Laplace-Beltrami spectra as Shape-DNA of surfaces and solids. Computer-Aided Design. 2006;38:342–366.

33. 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]

34. Rustamov RM. Laplace-beltrami eigenfunctions for deformation invariant shape representation. Eurographics Symposium on Geometry Processing (SGP) 2007:225–233.

35. Niethammer M, Reuter M, Wolter F-E, et al. Global medical shape analysis using the Laplace-Beltrami spectrum. Proc MICCAI. 2007:850–857. [PMC free article] [PubMed]

36. Vallet B, Lvy B. Spectral geometry processing with manifold harmonics. Computer Graphics Forum. 2008;27(2):251–260.

37. Shi Y, Lai R, Krishna S, et al. Anisotropic Laplace-Beltrami eigenmaps: Bridging Reeb graphs and skeletons. Proc MMBIA. 2008:1–7. [PMC free article] [PubMed]

38. 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]

39. 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]

40. 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]

41. Lai R, Shi Y, Scheibel K, Fears S, Woods R, Toga A, Chan T. Metric-induced optimal embedding for intrinsic 3D shape analysis. Proc CVPR. 2010:2871–2878.

42. Shi Y, Lai R, Toga A. CoRPORATE: Cortical reconstruction by pruning outliers with Reeb analysis and topology-preserving evolution. Proc IPMI. 2011:233–244. [PubMed]

43. Shi Y, Lai R, Gill R, Pelletier D, Mohr D, Sicotte N, Toga A. Conformal metric optimization on surface (CMOS) for deformation and mapping in Laplace-Beltrami embedding space. Proc MICCAI. 2011;2:327–334. [PubMed]

44. Reeb G. Sur les points singuliers d’une forme de Pfaff completement integrable ou d’une fonction nemérique. Comptes Rendus Acad Sciences. 1946;222:847–849.

45. Shinagawa Y, Kunii TL. Constructing a Reeb graph automatically from cross sections. IEEE Computer Graphics & Applications. 1991;11(6):44–51.

46. Cole-McLaughlin K, Edelsbrunner H, Harer J, Natarajan V, Pascucci V. Loops in Reeb graphs of 2-manifolds. Discrete Comput Geom. 2004;32(2):231–244.

47. Pascucci V, Scorzelli G, Bremer PT, Mascarenhas A. Robust on-line computation of Reeb graphs: simplicity and speed. ACM Trans Graph. 2007;26(3):58, 1–9.

48. Wood Z, Hoppe H, Desbrun M, Schroder P. Removing excess topology from isosurfaces. ACM Trans on Graph. 2004;23(2):190–208.

49. Patanè G, Spagnuolo M, Falcidieno B. A minimal contouring approach to the computation of the reeb graph. IEEE Trans Vis and Comp Graph. 2009;15(4):583–595. [PubMed]

50. Frangi A, Niessen W, Hoogeveen R, van Walsum T, Viergever M. Model-based quantitation of 3-D magnetic resonance angiographic images. IEEE Trans Med Imag. 1999;18(10):946–956. [PubMed]

51. Sato Y, Westin CF, Bhalerao A, Nakajima S, Shiraga N, Tamura S, Kikinis R. Tissue classification based on 3D local intensity structure for volume rendering. IEEE Trans Vis and Comp Graph. 2000;6(2):160–180.

52. Fischl B, Sereno MI, Dale AM. Cortical surface-based analysis ii: Inflation, flattening, and a surface-based coordinate system. NeuroImage. 1999;9(2):195–207. [PubMed]

53. Shattuck D, Mirza M, Adisetiyo V, et al. Construction of a 3D probabilistic atlas of human brain structures. NeuroImage. 2008;39(3):1064–1080. [PMC free article] [PubMed]

54. Klein A, Andersson J, Ardekani BA, Ashburner J, Avants B, Chiang MC, Christensen GE, Collins DL, Gee J, Hellier P, Song JH, Jenkinson M, Lepage C, Rueckert D, Thompson P, Vercauteren T, Woods RP, Mann JJ, Parsey RV. Evaluation of 14 nonlinear deformation algorithms applied to human brain mri registration. NeuroImage. 2009;46(3):786– 802. [PMC free article] [PubMed]

55. Siddiqi K, Bouix S, Tannebaum A, Zuker S. Hamilton-Jacobi skeletons. Int’l Journal of Computer Vision. 2002;48(3):215–231.

56. Shi Y, Karl WC. Real-time tracking using level sets. Proc CVPR. 2005;2:34–41.

57. 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]

58. Silvester PP, Ferrari RL. Finite Elements for Electrical Engineers. Cambridge University Press; 1996.

59. Uhlenbeck K. Generic properties of eigenfunctions. Amer J of Math. 1976;98(4):1059–1078.

60. Kimmel R, Sethian JA. Computing geodesic paths on manifolds. Proc Natl Acad Sci USA. 1998;95(15):8431–8435. [PubMed]

61. Aubert-Broche B, Griffin M, Pike G, Evans A, Collins D. Twenty new digital brain phantoms for creation of validation image data bases. IEEE Trans Med Imag. 2006;25(11):1410–6. [PubMed]

62. Mueller S, Weiner M, Thal L, et al. The Alzheimer’s disease neuroimaging initiative. Clin North Am. 2005;15:869–877. xi–xii. [PMC free article] [PubMed]

63. Mazziotta JC, Toga AW, Evans AC, et al. A probabilistic atlas and reference system for the human brain: international consortium for brain mapping. Philos Trans R Soc Lond B Biol Sci. 2001;356:1293–1322. [PMC free article] [PubMed]

64. Leung K, Parker D, Cunha A, et al. IRMA: An image registration meta-algorithm evaluating alternative algorithms with multiple metrics. Proc Int Conf on Scientific and Statistical Data Base Management. 2008:612–617.

65. Dinov I, Van Horn J, Lozev K, et al. Efficient, distributed and interactive neuroimaging data analysis using the LONI Pipeline. Front Neuroinform. 2009;3 [PMC free article] [PubMed]

66. Brainweb: Simulated Brain Database. [online] Available: http://www.bic.mni.mcgill.ca/brainweb.

67. Miller MI, Massie AB, Ratnanather J, Botteron KN, Csernansky JG. Bayesian construction of geometrically based cortical thickness metrics. NeuroImage. 2000;12(6):676– 687. [PubMed]

68. Sethian J. A fast marching level set method for monotonically advancing fronts. Proc Nat Acad Sci. 1996;93(4):1591–1595. [PubMed]

69. Cohen J. Quantitative methods in psychology: a power primer. Psychological Bulletin. 1992;112(1):155–59. [PubMed]

70. Dupont W, Plummer W. Power and sample size calculations: A review and computer program. Controlled Clinical Trials. 1990;11:116–28. [PubMed]

71. Reuter M, Schmansky NJ, Rosas HD, Fischl B. Within-subject template estimation for unbiased longitudinal image analysis. NeuroImage. 2012 doi: 10.1016/j.neuroimage.2012.02.084. in print. [PMC free article] [PubMed] [Cross Ref]