PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
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 2014 March 1.
Published in final edited form as:
PMCID: PMC3785796
NIHMSID: NIHMS444728

Cortical Surface Reconstruction via Unified Reeb Analysis of Geometric and Topological Outliers in Magnetic Resonance Images

Yonggang Shi, Rongjie Lai, Arthur W. Toga,corresponding author and Initiative the Alzheimer’s Disease Neuroimaging**

Abstract

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.

Index Terms: Cortical surface reconstruction, Reeb graph, topology correction, Laplace-Beltrami eigenfunctions, tissue classification

I. Introduction

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.

Fig. 1
An overview of our cortical reconstruction system.

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.

II. Enhancing Tissue Classification

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 Ψ : AR denote the skull-stripped MR image, where A is the lattice of voxels that we represent as A = {(i, j, k) [set membership] Z3|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 pGM, pWM, pCSF 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

equation M1
(1)

for all p = (i, j, k) [set membership] 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).

Fig. 2
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 [nabla]Ψ(p), and Hessian as:

equation M2
(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(p), and the corresponding eigenvector is denoted as vmax(p). At each voxel, we define the tissue enhancement feature as

equation M3
(3)

where the inner product between [nabla]Ψ(p) and vmax(p) is the gradient of Ψ along vmax(p) at the voxel. For voxels at the top of ridges or bottom of valleys, we have | < [nabla]Ψ(p), vmax(p) > | ≈ 0. Using the sign of λmax, we can determine whether it is a ridge or valley. With this tissue enhancement feature, we can help resolve ambiguities in intensity-based tissue classification. As an example, we show in Fig. 2(b) the tissue enhancement feature. Compared with the image in Fig. 2(a) and the original tissue map in Fig. 2(e), we can see that this feature successfully detects ridges and valleys, but there are also false positives that need to be cleaned.

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 equation M4 denote the warped GM atlas in the image space. As an illustration, the GM atlas equation M5 is shown in Fig. 2(c), and the warped anatomical labels are shown in Fig. 2(d). Using the anatomical labels, we denote RSUB = {p [set membership] A|p [set membership] (putamen, caudate, ventricle)} as the set of voxels belonging to the sub-cortical structures putamen, caudate and ventricles. Using the tissue enhancement feature and atlas labels, we can filter out outliers and detect CSF valleys and WM ridges.

Let RCSF = {p [set membership] A|Φ(p) ≤ 2, TE(p) > δCSF, equation M6, p [negated set membership] RSUB} be the set of candidate voxels on the CSF valley, where δCSF is a threshold we typically choose as 2. Because the CSF between touching gyri will occupy tiny fractions of voxels, we apply a thinning algorithm as listed in Table I to RCSF using −TE as the ridge/valley detection feature. Starting from the boundary voxel of RCSF with the least feature value TE, this thinning algorithm iteratively removes voxels if its feature TE is above the threshold δCSF and is a simple point [24], [55], i.e., peeling it off will not result in a topological change. As a result, we obtain a thin set of voxels RCSF with the maximal tissue enhancing feature.

TABLE I
Ridge/Valley Detection Algorithm

For partial WM voxels, we detect voxels on the top of ridges to follow the geometry of partial volume models. Let RWM = {p [set membership] A|2 < Φ(p) < 3, TE(p) < δWM, p [negated set membership] RSUB} be the set of candidate voxels on the ridge, where δWM is a threshold for ridge detection that we typically choose as −2. Applying the thinning algorithm to the set RWM with TE as the ridge/valley detection feature, we obtain a thin set of voxels RWM with partial WM.

Using these results, we can define an enhanced tissue map [Phi w/ hat] as

equation M7
(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:

equation M8
(5)

and the GM boundary is:

equation M9
(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.

III. Intrinsic Reeb Graphs on Tissue Boundaries

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.

A. Mesh representation

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 Ao and the background region Ab. To extract the continuous representation of the boundary between regions, we consider each grid point p = (i, j, k) [set membership] A as the center point of a rectangular cuboid An external file that holds a picture, illustration, etc.
Object name is nihms444728ig1.jpg(p) = {(x, y, z) [set membership] R3| |x−ihx| ≤ hx/2, |y − jhy| ≤ hy/2, |z − khz| ≤ hz/2}, where hx, hy, and hz are the spatial sampling resolutions in the x, y and z direction of the MR image. Under this cuboid representation, each voxel has six rectangular faces.

Because the well-composedness condition [30] is satisfied in our reconstruction [31], the boundary between Ao and Ab is a manifold and composed of a set of rectangular faces B = {B1, B2, · · ·, BNB}. Each face Bi is the intersection of the cuboids of two voxels: equation M10 and equation M11. This establishes a map from boundary faces to voxels in the object region equation M12 with equation M13 and the background region equation M14 with equation M15. By dividing each rectangular face Bi into two triangular faces An external file that holds a picture, illustration, etc.
Object name is nihms444728ig2.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms444728ig3.jpg, we have a triangular mesh representation of the boundary An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg = { An external file that holds a picture, illustration, etc.
Object name is nihms444728ig5.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig6.jpg}, where An external file that holds a picture, illustration, etc.
Object name is nihms444728ig5.jpg = { An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg|i = 1, · · ·, NV} and An external file that holds a picture, illustration, etc.
Object name is nihms444728ig6.jpg = { An external file that holds a picture, illustration, etc.
Object name is nihms444728ig8.jpg|i = 1, · · ·, 2NB} are the set of vertices and triangles. For each face An external file that holds a picture, illustration, etc.
Object name is nihms444728ig9.jpg, we know it is on the boundary face Bi with i = [left floor](j + 1)/2[right floor], where [left floor]x[right floor] denotes the greatest integer less than or equal to x. We can define the maps from each triangular face to voxels in Ao and Ab as: equation M16 and equation M17. These maps allow us to move freely between the surface representation and voxel representation of the tissue boundary under consideration.

B. Laplace-Beltrami eigen-system

With the triangular mesh representation of the boundary An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg = { An external file that holds a picture, illustration, etc.
Object name is nihms444728ig5.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig6.jpg}, we can compute its LB eigen-system by solving a matrix eigenvalue problem [31]–[33]:

equation M18
(7)

where λ is the eigenvalue, f : An external file that holds a picture, illustration, etc.
Object name is nihms444728ig5.jpgR 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:

equation M19

where An external file that holds a picture, illustration, etc.
Object name is nihms444728ig10.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg) is the set of vertices in the 1-ring neighborhood of An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig10.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig11.jpg) is the set of triangles sharing the edge ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig12.jpg), equation M20 is the angle in the triangle An external file that holds a picture, illustration, etc.
Object name is nihms444728ig13.jpg opposite to the edge ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig12.jpg), and Area( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig13.jpg) is the area of the l-th triangle An external file that holds a picture, illustration, etc.
Object name is nihms444728ig13.jpg.

The LB eigen-system is discrete and the eigenvalues can be ordered as 0 = λ0 ≤ λ1 ≤ λ2 ≤ · · ·. Correspondingly, the eigen-functions are f0, f1, f2, · · ·. 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 f1 here as the Morse function for Reeb graph construction [59], which is the solution of the minimization problem

equation M21
(8)

and can be viewed as the smoothest map from the manifold An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg to the real line R. As an illustration, the LB eigenfunction f1 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.

Fig. 5
Geometric and topological outlier detection with the Reeb graph of the LB eigenfunction. Zoomed views of geometric and topological outliers are plotted in (c) and (d), respectively. Yellow: Geometric outliers. Red: Topological outliers.

C. Reeb graph

Given a Morse function f on the mesh An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg, its Reeb graph is defined as follows [44].

Definition 1

Let f : An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpgR. The Reeb graph R(f) of f is the quotient space with its topology defined through the equivalent relation x [similar, equals] y if f(x) = f(y) for ∀x, y [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg.

Following this abstract definition, the Reeb graph is intuitively a graph of level contours of f on An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg. 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 An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms444728ig5.jpg, and its one-ring neighborhood An external file that holds a picture, illustration, etc.
Object name is nihms444728ig10.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg). Let An external file that holds a picture, illustration, etc.
Object name is nihms444728ig35.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg) = { An external file that holds a picture, illustration, etc.
Object name is nihms444728ig12.jpg [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms444728ig10.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg)|f ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig12.jpg) < f ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg)} and An external file that holds a picture, illustration, etc.
Object name is nihms444728ig34.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg) = { An external file that holds a picture, illustration, etc.
Object name is nihms444728ig12.jpg [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms444728ig10.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg)|f ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig12.jpg) > f( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg)} denote the lower and upper neighbors in the 1-ring neighborhood of a vertex An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg. Let # denote the number of connected components in a set. Using the number of connected components in An external file that holds a picture, illustration, etc.
Object name is nihms444728ig35.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg) and An external file that holds a picture, illustration, etc.
Object name is nihms444728ig34.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg), we can classify the vertices as follows:

equation M22
(9)

Let C = {C1, C2, · · ·, CN} be the set of critical points of f on the surface An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg. We assume all critical values are different and sort the critical points according to the critical values such that f(C1) < f (C2) < · · ·< f (CN). For anatomical shapes without perfect symmetry, we find this assumption always holds in our experience. For synthetic shapes with perfect symmetry, we can perturb the metric [37], [43] and make sure this assumption is valid. To accurately represent the partition of the surface by neighboring saddle points on the Reeb graph, which could have very subtle differences in the function values, we will augment the original mesh by splitting its triangles along the level contours during the Reeb graph construction process. We denote this augmented mesh as An external file that holds a picture, illustration, etc.
Object name is nihms444728ig14.jpg = ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig15.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig16.jpg), where An external file that holds a picture, illustration, etc.
Object name is nihms444728ig15.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms444728ig16.jpg are the augmented set of triangles and vertices of [M with circumflex]. The Reeb graph is represented as R(f) = (C, E), where C = {C1, C2, · · ·, CN} is the set of critical points and act as nodes of the graph, and E = {E1, E2, · · ·} is the set of arcs. Each arc in R(f) is represented as equation M23, where SCi, ECi [set membership] C are the start and end node, ETi [subset or is implied by] T is the set of triangles belonging to this arc in R(f), equation M24 and equation M25 are the level contours on the boundary of the arc. The set of triangles ETi represent the subset of An external file that holds a picture, illustration, etc.
Object name is nihms444728ig14.jpg enclosed by equation M26 and equation M27 and form a manifold with boundaries.

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 An external file that holds a picture, illustration, etc.
Object name is nihms444728ig15.jpg and vertices An external file that holds a picture, illustration, etc.
Object name is nihms444728ig16.jpg are initialized as the original triangle set An external file that holds a picture, illustration, etc.
Object name is nihms444728ig6.jpg and vertex set An external file that holds a picture, illustration, etc.
Object name is nihms444728ig5.jpg, respectively. We also assign a label function ArcLabel and initialize it to be −1 for all vertices in An external file that holds a picture, illustration, etc.
Object name is nihms444728ig16.jpg. We scan through these critical points sequentially as follows to build the Reeb graph.

If Ci is a minimum, we increase the current arc label j by one and create a new arc Ej. We set the start node of Ej as SCj = Ci, and its start level contour as equation M28. We also set ArcLabel(Ci) = j. Note that this arc is incomplete and we need to find the end node ECj, the end level contour equation M29, and the triangles ETj.

If Ci 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:

equation M30
(10)

For outgoing arcs, the isovalue is defined as:

equation M31
(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 An external file that holds a picture, illustration, etc.
Object name is nihms444728ig35.jpg (Ci) of the saddle point Ci, we trace a level contour at the value equation M32 on the mesh An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg. From the definition of equation M33, we know that this level contour crosses all edges between Ci and vertices in this component. This level contour is represented as a polyline P = (p1, p2, · · ·, pN) composed of points intersecting edges of An external file that holds a picture, illustration, etc.
Object name is nihms444728ig14.jpg at the level equation M34. We augment the mesh with this set of new vertices by adding edges. Let pk and pk+1 be two consecutive points that intersect the two edges of a triangle An external file that holds a picture, illustration, etc.
Object name is nihms444728ig17.jpg= ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig18.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig19.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig20.jpg) in An external file that holds a picture, illustration, etc.
Object name is nihms444728ig15.jpg. As shown in Fig. 3, we either keep this triangle intact or split it into two or three triangles and add them to An external file that holds a picture, illustration, etc.
Object name is nihms444728ig15.jpg. Repeating this process for all line segment in this level contour, we add all points to An external file that holds a picture, illustration, etc.
Object name is nihms444728ig16.jpg as new vertices in the mesh. Starting with these new vertices in P, we then grow backward using the algorithm in Table. II and complete an arc Ej in the Reeb graph R(f). We set ECj = Ci, equation M35, and ETj = Ξ, which are the set of triangles in An external file that holds a picture, illustration, etc.
Object name is nihms444728ig15.jpg in the augmented mesh An external file that holds a picture, illustration, etc.
Object name is nihms444728ig14.jpg that belong to this arc of the Reeb graph. We repeat this process until all arcs connecting Ci and An external file that holds a picture, illustration, etc.
Object name is nihms444728ig35.jpg (Ci) is crossed by a level contour.

Fig. 3
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.
TABLE II
Reeb Graph Arc Growing

For each component in the upper neighborhood An external file that holds a picture, illustration, etc.
Object name is nihms444728ig34.jpg (Ci) of the saddle point Ci, we trace a level contour at the value equation M36 that cross all edges connecting Ci and vertices in this component. Similarly, this level contour is represented as P = (p1, p2, · · ·, pN) and all intersecting points in this level contour are added to the augmented vertex set An external file that holds a picture, illustration, etc.
Object name is nihms444728ig16.jpg of the mesh An external file that holds a picture, illustration, etc.
Object name is nihms444728ig14.jpg. A new, but incomplete, arc in the Reeb graph is created and we increase the current arc label j by one. We set SCj = Ci, equation M37, and set ArcLabel(pk) = j for all newly added vertices pk [set membership] P. This process is repeated until all edges connecting An external file that holds a picture, illustration, etc.
Object name is nihms444728ig34.jpg (Ci) is crossed.

If Ci is a maximum, we grow backward with it as the starting point using the algorithm in Table. II to complete an arc Ej. We set ECj = Ci, equation M38, and the triangle set ETj = Ξ.

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 An external file that holds a picture, illustration, etc.
Object name is nihms444728ig14.jpg, which has non-uniform triangles generated from the triangle-splitting process, we can build a map μ: An external file that holds a picture, illustration, etc.
Object name is nihms444728ig15.jpgAn external file that holds a picture, illustration, etc.
Object name is nihms444728ig6.jpg to relate the properties derived from R(f) to the original mesh An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg that has regular triangles. For each triangle An external file that holds a picture, illustration, etc.
Object name is nihms444728ig21.jpg [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms444728ig15.jpg, the map μ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig21.jpg) denotes the triangle in An external file that holds a picture, illustration, etc.
Object name is nihms444728ig6.jpg such that An external file that holds a picture, illustration, etc.
Object name is nihms444728ig21.jpg is a subset. This is a many-to-one map since the triangles in An external file that holds a picture, illustration, etc.
Object name is nihms444728ig15.jpg is obtained by splitting triangles in An external file that holds a picture, illustration, etc.
Object name is nihms444728ig6.jpg. By establishing a map between the triangles of An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms444728ig14.jpg, we build a bridge between the Reeb graph on the augmented surface An external file that holds a picture, illustration, etc.
Object name is nihms444728ig14.jpg 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 An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg 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 Ci is illustrated. One level contour is generated that crosses edges connecting the lower neighborhood An external file that holds a picture, illustration, etc.
Object name is nihms444728ig35.jpg (Ci) and this saddle point. Two level contours are generated crossing the edges connecting the upper neighborhood An external file that holds a picture, illustration, etc.
Object name is nihms444728ig34.jpg (Ci) and Ci. The arc construction process is illustrated in Fig. 4 (b), where all arcs are plotted in different colors. For each arc, we also plot its two nodes as black dots. In the first step, an arc between a saddle point and a minimum is constructed. In the 2nd, 3rd, and 4th step, arcs between saddle points are added sequentially to the Reeb graph. In the final step, an arc between a maximum and a saddle point is constructed. The final Reeb graph is visualized in Fig. 4 (c) by plotting triangles on each arc of R(f) with corresponding colors used in Fig. 4 (b), where the graph structure is evident from the neighboring relation of the arcs.

Fig. 4
Reeb graph construction of a double torus using a LB eigenfunction. (a) An illustration of the level contours in the neighborhood of a saddle point. (b) An illustration of the arc construction process in building the Reeb graph. (c) The final partition ...

IV. Unified Analysis of Geometric and Topological Outliers

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.

A. Paired boundary estimates

Given the WM evolution speed FWM, 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 equation M39 and the background region is denoted as equation M40. The triangular mesh representation of the boundary is An external file that holds a picture, illustration, etc.
Object name is nihms444728ig22.jpg = ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig23.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig24.jpg). For topological analysis, we turn off the topology-preserving constraint and continue the boundary evolution process under the same speed FWM to obtain the second boundary estimate, which can have arbitrary topology. The object and background region determined by this boundary are denoted as equation M41 and equation M42, and the boundary mesh is denoted as An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg = ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig26.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig27.jpg).

By comparing the paired boundary estimates An external file that holds a picture, illustration, etc.
Object name is nihms444728ig22.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg, we can locate paired patches on An external file that holds a picture, illustration, etc.
Object name is nihms444728ig22.jpg that together fill a handle or tunnel in An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg. For all faces on An external file that holds a picture, illustration, etc.
Object name is nihms444728ig22.jpg with their interior voxel in the background region of An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg, i.e., equation M43, where the map equation M44 was defined in section III-A, we perform a connected component labeling and denote the set of connected components as equation M45. For faces in the k-th component equation M46, we define a connected component (CC) label equation M47. To find paired patches on An external file that holds a picture, illustration, etc.
Object name is nihms444728ig22.jpg that jointly fill a tunnel or handle, we want to find the set of filling voxels that connect different components in CCG and fill topological handles or tunnels. For every voxel p in equation M48, which is the set of all interior voxels of faces in CCG, we define a set PL(p) as the paired labels of connected components it bridges. For each face equation M49, we first assign its CC label to equation M50. Note that there can be at most two components sharing an interior voxel, so there cannot be more than two labels in the paired label set PL(p) for any voxel p. Because a subset of the voxels were included in equation M51 to satisfy the well-composedness constraint, they do not directly touch two opposing patches in An external file that holds a picture, illustration, etc.
Object name is nihms444728ig22.jpg. For each equation M52 with #PL(p)=1, we search its six-neighborhood An external file that holds a picture, illustration, etc.
Object name is nihms444728ig28.jpg (p). If there is a voxel q [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms444728ig28.jpg (p) that satisfies #(PL(q) [union or logical sum] PL(p)) = 2, which means they together form a bridge to connect two patches, we set PL(p) and PL(q) as the union of these two sets. For each face equation M53, we assign it the paired label equation M54 of its interior voxel if equation M55. Using paired labels, we group triangles with the same paired labels into paired components equation M56 and denote equation M57 as the set of filling voxels for topological artifacts. As an illustration, we show in Fig. 7 (b) and (d), the filling voxels for a handle and tunnel, respectively. By analyzing the underlying tissue properties on filling voxels, cut or fill decisions could be made in topological analysis with the inclusion of information from tissue classification.

Fig. 7
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 An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg, we compute its LB eigenfunctions and construct the Reeb graph. Let f = f1 be the first LB eigenfunction of An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg, as illustrated for a WM boundary in Fig. 5 (a), and R(f) = (C, E) the Reeb graph of f on An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg, 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 [M with circumflex]F = ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig29.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig30.jpg). The map from the triangles on An external file that holds a picture, illustration, etc.
Object name is nihms444728ig31.jpg to An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg is denoted as μ: An external file that holds a picture, illustration, etc.
Object name is nihms444728ig29.jpgAn external file that holds a picture, illustration, etc.
Object name is nihms444728ig26.jpg. Each arc of the Reeb graph equation M58, where SCi, ECi [set membership] C are the start and end node, and ETi [subset or is implied by] An external file that holds a picture, illustration, etc.
Object name is nihms444728ig29.jpg is the set of triangles on this arc in R(f), equation M59 and equation M60 are the level contours at the start and end node. The degree of a node in the Reeb graph is the number of arcs that it is either the start or the end node. A node is called a leaf node if its degree is one, and an internal node if its degree is greater than one. The arc connected to a leaf node is called a leaf arc.

B. Geometric outlier correction

We detect geometric outliers from the set of leaf arcs in the Reeb graph. Left Ei be a leaf arc with equation M61 and equation M62 as the start and end level contour. Using the level contour equation M63 and equation M64, we define an arc feature H(Ei) as follows:

equation M65
(12)

where n is the outward normal on the surface An external file that holds a picture, illustration, etc.
Object name is nihms444728ig31.jpg, equation M66 and equation M67 are the mean coordinates of points on equation M68 and equation M69. If H(Ei) > 0, it means Ei is an outward leaf pointing toward the exterior of the cortex; otherwise, it is an inward leaf that points toward the interior of the cortex.

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

equation M70
(13)

where Area(ETi) is the sum of area of all triangles in ETi, the parameters ρ1 and ρ2 are thresholds selected to identify sharp and small outliers. To further localize geometric outliers, we project the original mesh An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg onto a subset of its LB eigenfunctions S = {f0, f1, f2, f3} and calculate the area distortion of triangles. For each triangle in An external file that holds a picture, illustration, etc.
Object name is nihms444728ig26.jpg, its area distortion factor (ADF) is defined as

equation M71
(14)

where equation M72 and equation M73 are the area of the triangle equation M74 in the original mesh An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg and the projected mesh An external file that holds a picture, illustration, etc.
Object name is nihms444728ig31.jpg. For an outlier leaf arc, we define the set of outlier triangles as Outlier(ETi) = { An external file that holds a picture, illustration, etc.
Object name is nihms444728ig8.jpg [set membership] μ(ETi)|ADF ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig8.jpg) > γ}, where μ(ETi) is the set of triangles obtained by mapping the triangles in ETi to the original mesh An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg, and the parameter γ is a threshold used for further localization of geometric outliers to triangles that exhibit large area distortions during the projection onto the subspace of LB eigenfunctions [31]. To localize outliers on sub-cortical surfaces, the method in [31] relies exclusively on eigen-projection and needs to use hundreds of eigen-functions to form the subspace. For the much more complicated cortical surfaces, we observe in our experiments that a much larger number of eigen-functions are needed to achieve similar localization of outliers. This is computationally very expensive and clearly infeasible for efficient processing because the number of vertices on a cortical surface is typically two orders of magnitude larger than many sub-cortical structures. If we use the same subspace S of the first four eigen-functions, the eigen-projection method in [31] would produce a large number of false positives, especially on the frontal lobe, as shown in Fig. 5(b). In contrast, we show in Fig. 5(c) and (d) that we can achieve much better performance with the Reeb analysis approach proposed here using the same number of eigen-functions. A zoomed view of detected geometric outliers is plotted in Fig. 5(c), which shows that the proposed method achieves accurate outlier localization with very few eigen-functions.

To remove the outlier, we modify the evolution speed FWM as follows. If Ei is an outward leaf, we set equation M75 for all An external file that holds a picture, illustration, etc.
Object name is nihms444728ig8.jpg [set membership] Outlier(ETi) so that its interior will be removed. If Ei is an inward leaf, we set equation M76 for all equation M77 such that its exterior will be filled and the outlier can be removed. In contrast to previous methods that rely on global smoothness regularization to avoid geometric outliers, our method only modifies the evolution speed locally while leaving other parts of the boundary intact. This greatly decreases shrinkage effect and help improve the accuracy in surface reconstruction.

C. Topological outlier correction

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 Ei [set membership] R(f), we set W (i1, i2) = Area(ETi), where i1 and i2 are the indices of the start and end critical points, respectively, and Area(ETi) is the area of this arc.

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 Ci in the directed graph W, if there are more than one out-going arcs, we use breadth-first-search (BFS) to test if an outlier exists. Let j1, · · ·, jK denote the set of neighboring nodes which satisfies W (i, j1) > 0, · · ·, W (i, jK) > 0. For every outgoing arc, we set W (i, jk) = 0 and perform BFS with Ci as the starting node. This generates a spanning tree of nodes reachable from Ci after the removal of the arc (Ci, Cjk). The spanning tree is recorded as an array PARk. For each node Cj, PARk(Cj ) is the parent node in this spanning tree if it is reached by Ci from PARk(Cj). By repeating this step for all outgoing arcs at Ci, we generate a set of BFS trees PAR1, …, PARK. The intersection node CJmin of these trees is defined as the node with the smallest index j such that PARk(Cj) > 0 for all k = 1, · · ·, K. Starting from the intersection node CJmin, we trace backward to the current node Ci on each BFS tree PARk and record the path as PATHk = (Ci, · · ·, PARk(CJmin), CJmin). The cost of this PATHk is the total area of all arcs on this path. Among all paths we pick the one with the smallest cost and add all arcs on this path to the outlier set TO as topological outliers. As an illustration, the topological outliers detected with Reeb analysis on the surface in Fig. 5(a) are plotted in red in Fig. 5(c) and (d). A zoomed view of a topological outlier is plotted in Fig. 5(d).

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(Ei) as defined in (12) for each arc Ei [set membership] TO. If H(Ei) > 0, which means the centroids of the start and end level contours are enclosed by the surface as in the case plotted in Fig. 6(a), we classify Ei as a handle. If H(Ei) < 0, which means the centroids of the start and end level contours fall outside the surface as shown in Fig. 6(b), we classify Ei as a tunnel. For a handle or tunnel, we use a different method to find a cutting path with minimal length if a cut operation is needed. If Ei is a handle, we uniformly sample a set of level contours of the eigenfunction f between the start contour equation M78 and end contour equation M79, and pick the contour with the least length as the cutting path. For example, a handle on a WM boundary is plotted in red in Fig. 7(a) and the cutting path as the shortest level contour is plotted in Fig. 7(b). For a tunnel, we combine two geodesic paths to form a cutting path. The first geodesic goes from the start to the end node within the triangles ETi. The second geodesic goes from the start to the end node without passing through the triangles in ETi. For a tunnel on a WM boundary as shown in Fig. 7(c), the cutting path that is composed of two geodesics is plotted in Fig. 7(d).

Fig. 6
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 An external file that holds a picture, illustration, etc.
Object name is nihms444728ig22.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg, we have paired patches equation M80 on the genus-zero boundary estimate An external file that holds a picture, illustration, etc.
Object name is nihms444728ig22.jpg. These patches fill the handles and tunnels on An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg. For each paired patches, its interior voxels are equation M81 and they are filled inside equation M82 to satisfy the genus-zero constraint, which are illustrated in Fig. 7 for a handle and tunnel. At a paired patch equation M83, a cut decision should be made if either of the following two conditions is met

  • If the number of voxels equation M84 with [Phi w/ hat](p) ≤ 1 is greater than THDCSF.
  • If the number of triangles equation M85 with ADF ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig9.jpg) > γ is greater than THDGEO.

Both parameters THDCSF and THDGEO are thresholds we choose empirically to identify topological outliers that are not consistent with the underlying image intensity distributions and geometric regularity of cortical surfaces. The first condition checks if filling a handle or tunnel needs more than THDCSF voxels classified as CSF in the enhanced tissue map. The second condition measures the geometric regularity of the filling patches by calculating if there are more than THDGEO triangles with its ADF, which is defined in (14), greater than the threshold γ used for geometric outlier detection.

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 An external file that holds a picture, illustration, etc.
Object name is nihms444728ig22.jpg, we pick a cutting path on a handle or tunnel that is connected to this paired component. For the cutting path CPj of an arc Ej in R(f), we denote An external file that holds a picture, illustration, etc.
Object name is nihms444728ig26.jpg (CPj) as the set of triangles it passes on An external file that holds a picture, illustration, etc.
Object name is nihms444728ig25.jpg. We consider a paired component equation M86 connected to a cutting path CPj if the set equation M87 is not empty, which means the exterior voxels of faces on the cutting path intersects with the interior voxels of the paired patches. Among all cutting paths intersecting the paired component, we choose the one with the shortest length and denote it as CPJmin. To cut it open, we modify the evolution speed as

equation M88
(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.

D. WM and GM boundary estimate

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 An external file that holds a picture, illustration, etc.
Object name is nihms444728ig22.jpg 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).

Fig. 8
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).

V. Sub-voxel Accuracy with Adaptive Interpolation

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 Φl, we denote B = (B1, B2, · · ·, BNB) as the set of rectangular faces and An external file that holds a picture, illustration, etc.
Object name is nihms444728ig4.jpg = ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig5.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms444728ig6.jpg) the triangular mesh representation of the boundary surface as defined in section III-A. The pair (Φh, Φl) could be (3, 2) for the WM boundary or (2, 1) for the GM boundary. For a face Bi, its interior voxel neighbor is equation M89 and exterior voxel neighbor is equation M90. We denote BVin and BVout as the set of interior and exterior neighboring voxels of all rectangular faces on the boundary, respectively.

For each voxel p [set membership] BVin [union or logical sum] BVout, its cuboid An external file that holds a picture, illustration, etc.
Object name is nihms444728ig1.jpg (p) has two faces in each of the x, y, and z direction. To calculate the shift of faces for partial volume modeling, we count the number of faces that the intersection of An external file that holds a picture, illustration, etc.
Object name is nihms444728ig1.jpg (p) and B, i.e., An external file that holds a picture, illustration, etc.
Object name is nihms444728ig1.jpg (p) ∩ B, have in the x, y, z direction and denote them as Nx(p), Ny(p), and Nz(p), respectively. We also estimate locally the image intensity Ψh(p) for the high intensity tissue Φh as the mean intensity of voxels classified as Φh in the 5 × 5 × 5 neighborhood centered at p. Similarly, the image intensity Ψl(p) of the low intensity tissue Φl is estimated locally as the mean intensity of voxels classified as Φl in the same neighborhood. The fraction of high intensity tissue Φh contained in this voxel p can then be computed with linear interpolation as:

equation M91
(16)

To account for this partial volume effect, we will shift faces on the boundary. Let β(p) = min(1, Nx(p)) + min(1, Ny(p)) + min(1, Nz(p)) be the total number of directions that we will shift faces. The shrink factor in each direction is then α(p)1(p).

For a face Bi, if the tissue map of its interior voxel equation M92, it needs to move inward according to the shift factor α′ (p*) = 1 − α(p*)1(p*) computed with equation M93. If equation M94, it needs to move outward according to the partial volume shift factor α′ (p*) = α(p*)1(p*) computed from equation M95. Let equation M96 denote the coordinate of the point p*. We also denote equation M97 as the projection of p* onto the face Bi, which is the nearest point on the face Bi. The shift for this face is

equation M98
(17)

where the shift in each direction is multiplied by the corresponding spatial resolution, and divided by the number of faces that An external file that holds a picture, illustration, etc.
Object name is nihms444728ig1.jpg (p*) ∩ B has in that direction.

Given the shift of all faces in B, we calculate the shift for a vertex An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms444728ig5.jpg as the average shift of its neighboring faces:

equation M99
(18)

where An external file that holds a picture, illustration, etc.
Object name is nihms444728ig32.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg) denote the set of rectangular faces in the neighborhood of the vertex An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg, and # An external file that holds a picture, illustration, etc.
Object name is nihms444728ig32.jpg ( An external file that holds a picture, illustration, etc.
Object name is nihms444728ig7.jpg) is the number of rectangular faces in this set.

Let x denote the vector of coordinates of all vertices, and An external file that holds a picture, illustration, etc.
Object name is nihms444728ig36.jpg the shift of all vertices. We compute the final vertex coordinates with sub-voxel accuracy by minimizing the following energy function:

equation M100
(19)

where 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:

equation M101
(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.

VI. Experimental Results

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 = 100mm, ρ2 = 5, γ = 100, THDCSF = 5, THDGEO = 5, η = 10. For FreeSurfer, the default settings were used. All computations are carried out on the grid of the LONI pipeline. In general, our system takes around 2–4 hours to reconstruct surfaces from each subject. For FreeSurfer, the surface reconstruction process typically takes around 10 to 20 hours, and the whole workflow, including surface labeling, can take 20 to 30 hours.

A. Qualitative comparisons

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.

Fig. 9
The iterative process of outlier removal with Reeb analysis.

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.

Fig. 10
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 ...
Fig. 11
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 ...
Fig. 12
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.

Fig. 13
Comparison of cortical reconstruction results on four ADNI and four ICBM cases. For each case, the result from our method is plotted on the top, and the result from FreeSurfer is plotted on the bottom. Different views are selected for different subjects ...

B. Atrophy detection on simulated MR images

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.

Fig. 14
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.

Fig. 15
Projection of FreeSurfer gyral labels onto the GM surface reconstructed by our method.
Fig. 16
Gyrus-based map of p-values from the testing of longitudinal atrophy on LH and RH surfaces from our results (a, b) and FreeSurfer results (c, d).

C. Quantitative comparisons on ADNI data

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.

Fig. 17
Mean thickness of gyral regions on the left and right hemisphere of the NC and AD group.
Fig. 18
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 equation M102 [69], where MTNC and MTAD were the mean thickness of the NC and AD group, respectively, and σ was the pooled standard deviation. The power was computed with the software tool PS [70]. For both our method and FreeSurfer, we counted the number of gyri in each hemisphere with power greater than four thresholds 0.99, 0.95, 0.90, 0.85. The results are listed in Table III. We can see from the power analysis that our method is able to generate more sensitive thickness features on both hemispheres for all thresholds. From a study design point of view, this shows the thickness features derived from our method has the potential of reducing sample sizes in experiments for AD research.

TABLE III
Power analysis of the ADNI experiment: number of thickness features from our method and FreeSurfer with power above thresholds.

D. Quantitative comparisons on ICBM data

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.

Fig. 19
Gyrus-based thickness map of the young and elderly group.
Fig. 20
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 equation M103, where MTY and MTO were the mean thickness of the young and elderly groups, and σ was the pooled standard deviation. For each method, we counted the number of gyrus on each hemisphere with power exceeding the thresholds. The results were listed in Table IV. On this data set, we can see our method is able to generate a much larger number of sensitive thickness features than FreeSurfer. This demonstrates that our method achieves better performance on both the ADNI and ICBM datasets with very different demographic distributions.

TABLE IV
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.

TABLE V
Power analysis on repeat scans of ICBM dataset: number of thickness features from our method with power above thresholds.

VII. Conclusions

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

Acknowledgments

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.

Footnotes

1Part 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.

Contributor Information

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.

References

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]