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

**|**HHS Author Manuscripts**|**PMC3818082

Formats

Article sections

- Abstract
- 1 Introduction
- 2 Intrinsic Reeb Graph
- 3 Unified Geometry and Topology Correction
- 4 Experimental Results
- 5 Conclusions
- References

Authors

Related links

Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2013 November 5.

Published in final edited form as:

Med Image Comput Comput Assist Interv. 2012; 15(0 1): 601–608.

PMCID: PMC3818082

NIHMSID: NIHMS444724

See other articles in PMC that cite the published article.

A key challenge in the accurate reconstruction of cortical surfaces is the automated correction of geometric and topological outliers in tissue boundaries. Conventionally these two types of errors are handled separately. In this work, we propose a unified analysis framework for the joint correction of geometric and topological outliers in cortical reconstruction. Using the Reeb graph of intrinsically defined Laplace-Beltrami eigenfunctions, our method automatically locates spurious branches, handles and holes on tissue boundaries and corrects them with image information and geometric regularity derived from paired boundary evolutions. In our experiments, we demonstrate on 200 MR images from two datasets that our method is much faster and achieves better performance than FreeSurfer in population studies.

Automated reconstruction of cortical models from MR images is fundamental for large scale brain mapping. While many approaches were proposed [1–5], critical challenges remain in improving the accuracy, robustness, and speed of reconstruction algorithms. In this paper, we propose a unified analysis framework based on intrinsic Reeb graphs that jointly tackles the correction of geometric and topological outliers in cortical reconstruction. We demonstrate that our method can efficiently reconstruct accurate cortical surfaces and achieve better performance in large scale population studies [6, 7].

Cortical reconstruction is a complicated system development problem that requires the successful integration of various preprocessing steps including non-uniformity correction, skull stripping, non-linear registration, the labeling of sub-cortical structures, and tissue classification that classifies image intensities into white matter(WM), gray matter(GM) and cerebrospinal fluid(CSF). At the core of this system, though, is the generation of smooth and topologically correct mesh models of tissue boundaries. In previous works, smoothness is typically achieved with the incorporation of a global regularization in surface evolution. Topology correction, on the other hand, was handled separately with graph-based methods [3,4] or spherical mapping [2]. Using Laplace-Beltrami (LB) eigenfunctions [8–10] and their Reeb graphs, accurate surface reconstruction methods were proposed recently [10, 12], but topological outliers were not considered. Based on the theory of Morse functions and Reeb graphs, we propose in this work a unified analysis framework that corrects both geometric and topological outliers in cortical reconstruction. Using the Reeb graph of LB eigenfunctions, our method can localize geometric and topological outliers and make decisions about correction strategies with information from tissue classification and geometric regularity. After the correction, accurate surface models are generated with adaptive interpolation.

The rest of the paper is organized as follows. In section 2, we develop a novel approach for the construction of Reeb graph of LB eigenfunctions on triangular meshes. The unified approach for topology and geometry correction is developed in section 3. Experimental results are presented in section 4, where we demonstrate that our method can achieve better performance than FreeSurfer [2]. Finally conclusions are made in section 5.

For intrinsic shape analysis, the Reeb graphs of LB eigenfunctions have been successfully applied for shape modeling and geometric outlier detection on genus-zero surfaces [12]. In this section, we develop a novel and general method for Reeb graph construction on surfaces of arbitrary topology. Instead of scanning through all level sets on surfaces [11], our method is efficient and uses only level sets in the neighborhood of saddle points.

Let
= (
,
) be a triangular mesh, where
and
are the set of vertices and triangles. Given a function *f*:
→ , the Reeb graph *R*(*f*) can be viewed as a graph of critical points on meshes [11]. In this work, we choose the function *f* as the intrinsically defined eigenfunctions of the Laplace-Beltrami operator on surfaces, which have the advantage of being invariant to scale and pose variations. Let
= {*c*_{1}, *c*_{2}, · · ·, *c _{N}*} be the set of critical points of

To construct the Reeb graph, we scan through all critical points sequentially. If *c _{i}* is a minimum, we add a new, but incomplete, edge to

$$\begin{array}{l}{f}_{\mathit{iso}}^{-}({c}_{i})=(f({c}_{i})+max(f({{c}_{i}}_{-1}),\underset{{\mathcal{V}}_{j}\in {N}_{-}({c}_{i})}{max}f({\mathcal{V}}_{j})))/2\\ {f}_{\mathit{iso}}^{+}({c}_{i})=(f({c}_{i})+min(f({{c}_{i}}_{+1}),\underset{{\mathcal{V}}_{j}\in {N}_{+}({c}_{i})}{max}f({\mathcal{V}}_{j})))/2.\end{array}$$

(1)

For each component in *N*_{+}(*c _{i}*), we trace a level contour

As an example, we show in Fig. 1 the construction of the Reeb graph on a double torus. In Fig. 1(a), we can see one level contour crosses the lower neighborhood, and two level contours cross the upper neighborhood of the saddle point. The edges of the Reeb graph are shown in Fig. 1(b), where the graph structure is evident from the neighboring relation of the edges and we can see the loops in the graph captures the topology of the surface as suggested in Morse theory.

For the analysis of a 3D MR image, we denote the lattice of voxels as *A* = {(*i*, *j*, *k*) ^{3}|0 ≤ *i* < *I* − 1, 0 ≤ *j* < *J* − 1, 0 ≤ *k* < *K* − 1} and the 3D image as a function *A* → . Following the preprocessing steps in [12], we can define an evolution speed *F* on the lattice *A* based on tissue classification of MR images, and use the fast evolution algorithm in [10, 13] to find the boundary between two tissue types. Let *A _{o}* and

For topological analysis, we compute a pair of boundary estimates with the evolution speed *F*. The first boundary satisfies the genus-zero condition and is denoted as
= (
,
). The second boundary can have arbitrary topology and is denoted as
= (
,
), which is obtained by turning off the topology-preserving condition in the evolution algorithm. Let *f* be the first non-constant LB eigenfunction on
, and *R*(*f*) = (*C*, *E*) the Reeb graph of *f* on
that we build with the algorithm in section 2. Here *C* denotes the set of critical points, and *E* the set of edges. Each edge is represented as
${E}_{i}=({SC}_{i},{EC}_{i},{EV}_{i},{LC}_{i}^{s},{LC}_{i}^{e})$, where *SC _{i}*,

The degree of a node in the Reeb graph is the number of edges 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 edge connected to a leaf node is called a leaf edge. We detect geometric outliers from the set of leaf edges in the Reeb graph. Let *E _{i}* be a leaf edge,
${LC}_{i}^{s}$ and
${LC}_{i}^{e}$ denote the level contours at the start and end node that are used to generate the edge. Let
(

To remove the outlier, we modify the evolution speed *F* as follows. For each outlier edge *E _{i}*, we use the level contour
${LC}_{i}^{s}$ and
${LC}_{i}^{e}$ to compute:

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

(2)

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

For topological analysis, we first remove duplicated edges in *R*(*f*). For any two edges with the same start and end node, we remove the one with smaller size from *R*(*f*) and add it to the set of topological outlier which we denote as *TO*. After that, we represent the Reeb graph *R*(*f*) as a matrix *W*. For any edge
${E}_{i}=({C}_{i1},{C}_{i2},{EV}_{i},{LC}_{i}^{s},{LC}_{i}^{e})$, we set *W*(*i*1, *i*2) = #(*EV _{i}*), which is the number of vertices on this edge.

For a node *C _{i}* in the directed graph

For an edge *E _{i}*

To make cut or fill decisions about topological outliers, we find paired patches on
that jointly fill a hole or handle in
. For all faces on
with their interior voxels in the background region of
, we perform a connected component labeling and denote the set of patches as
${CC}^{G}=\{{CC}_{1}^{G},{CC}_{2}^{G},\cdots \}$. Triangles in different patches are paired if they share a common interior voxel. These triangle pairs are then grouped into paired components
${PC}^{G}=\{{PC}_{1}^{G},{PC}_{2}^{G},\cdots \}$ according to their paired patch labels in *CC ^{G}*. For each paired patches, its interior voxels are
${\mathbf{x}}_{o}({PC}_{i}^{G})$ and they are filled inside
to satisfy the genus-zero constraint. A cut decision should be made if either conditions is met:

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

The first condition checks if filling a handle/hole needs voxels classified as CSF. The second condition measures the geometric regularity of the filling patches. To cut open a paired component in
, we search for a handle or hole with the shortest cutting path that is connected to this paired component. Let *CP _{i}* denote the cutting path of this hole or handle and the set of triangles it passes on
as
(

As an illustration, we plotted the detected outliers on a WM boundary with Reeb analysis in Fig. 2. By iteratively applying the steps in section 3.1, 3.2, 3.3, we can remove both geometric and topological outliers in a unified framework.

Let *M ^{G}* represent the genus-zero, triangular mesh representation of the boundary between high intensity tissue

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

(3)

where **x** is the vector of coordinates of all vertices, *Δ* is the discrete Laplacian matrix of the mesh, and *η* is a regularization parameter. The solution of this quadratic problem gives the coordinates of vertices on the final cortical surface:

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

(4)

In this section, we demonstrate our method and compare with FreeSurfer [2], which is widely used in brain imaging research, on two large datasets. The first dataset includes skull-stripped MR images of 50 normal controls(NC) and 50 Alzheimer’s disease (AD) patients from ADNI [6]. The second dataset includes 100 skull-stripped MR images from ICBM [7] and it has a wide age range from 19 to 80. For all images, the same set of parameters are used in our method: *α* = 100*mm*, *β* = 5, *γ* = 100*, THD _{CSF}* = 5,

In the first experiment, we present detailed comparisons using one MR scan of an AD patient in ADNI. 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 more than 10 hours. In Fig. 3, the left hemispherical WM surfaces reconstructed by both methods are plotted. In regions highlighted by the dashed curves, 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. 3(c). The left hemispherical GM surfaces reconstructed with our method and FreeSurfer are plotted in Fig. 4(a). The intersections of the surfaces with two sagittal slices shown in Fig. 4(b) illustrate that our surface can better capture deep sulcal regions.

A comparison of WM surfaces. (a) Our result. (b) FreeSurfer result. (c) Intersection of image slices with highlighted regions on surfaces. Red: our method. Blue: FreeSurer.

A comparison of GM surfaces. (a) Left: Our result. Right: FreeSurfer result. (b) Intersection of sagittal slices with surfaces. Red: our method. Blue: FreeSurfer.

In the second experiment, we applied both our method and FreeSurfer to the MR images of 50 NC and 50 AD from ADNI. Using a general surface mapping algorithm [14], we projected the gyral labels of the LPBA40 atlas [15] to all left hemispherical surfaces. For each gyrus, we computed the average gray matter thickness and tested group differences between NC and AD using two tailed t-tests. The p-values from our method and FreeSurfer are plotted in Fig. 5. Our method can achieve more significant p-values on 18 of the 24 gyral regions.

P-value maps of gyrus-based thickness differences between NC and AD. (a) Our results. (b) FreeSurfer results.

In the third experiment, we applied our method and FreeSurfer to the 100 MR images from the ICBM database and used the results to model the decrease of gray matter thickness in normal aging. The reconstructed surfaces were first automatically labeled into gyral regions with the same approach in the second experiment. Regression analysis was then applied to the mean thickness of each gyrus against subject age. The results are shown in Fig. 6, where the rates of decrease are plotted in Fig. 6(a) and (b), and the p-values of the regression analysis are plotted in Fig. 6(c) and (d). We can see that our method can achieve statistically much more significant results on all gyral regions.

In this paper we developed a novel approach for the unified correction of geometric and topological outliers in cortical surface reconstruction. Comparisons with a state-of-the-art tool showed that our method is computationally more efficient and can achieve better performance in population studies. The Reeb analysis framework developed here is general and can be valuable for general surface reconstruction and analysis problems.

1. Mangin JF, Frouin V, Bloch I, Regis J, Lopez-Krahe J. 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. Dale AM, Fischl B, Sereno MI. Cortical surface-based analysis i: segmentation and surface reconstruction. NeuroImage. 1999;9:179–194. [PubMed]

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

4. Han X, Pham DL, Tosun D, Rettmann ME, Xu C, Prince JL. CRUISE: Cortical reconstruction using implicit surface evolution. NeuroImage. 2004;23:997–1012. [PubMed]

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

6. Mueller S, Weiner M, Thal L, Petersen RC, Jack C, Jagust W, Trojanowski J, Toga A, Beckett L. The Alzheimer’s disease neuroimaging initiative. Clin North Am. 2005;15:869–877. xi–xii. [PMC free article] [PubMed]

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

8. Reuter M, Wolter FE, Shenton M, Niethammer M. Laplace-beltrami eigenvalues and topological features of eigenfunctions for statistical shape analysis. Computer-Aided Design. 2009;41(10):739–55. [PMC free article] [PubMed]

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

10. 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–22. [PMC free article] [PubMed]

11. 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–44.

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

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

14. 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–34. [PubMed]

15. Shattuck D, Mirza M, Adisetiyo V, Hojatkashani C, Salamon G, Narr K, Poldrack R, Bilder R, Toga A. Construction of a 3D probabilistic atlas of human brain structures. NeuroImage. 2008;39(3):1064–1080. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |