PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2011 January 1.
Published in final edited form as:
Med Image Comput Comput Assist Interv. 2010; 13(Pt 3): 49–56.
PMCID: PMC2970518
NIHMSID: NIHMS242596

Automated Sulci Identification via Intrinsic Modeling of Cortical Anatomy

Abstract

In this paper we propose a novel and robust system for the automated identification of major sulci on cortical surfaces. Using multiscale representation and intrinsic surface mapping, our system encodes anatomical priors in manually traced sulcal lines with an intrinsic atlas of major sulci. This allows the computation of both individual and joint likelihood of sulcal lines for their automatic identification on cortical surfaces. By modeling sulcal anatomy with intrinsic geometry, our system is invariant to pose differences and robust across populations and surface extraction methods. In our experiments, we present quantitative validations on twelve major sulci to show the excellent agreement of our results with manually traced curves. We also demonstrate the robustness of our system by successfully applying an atlas of Chinese population to identify sulci on Caucasian brains of different age groups, and surfaces extracted by three popular software tools.

1 Introduction

The automated identification of major sulci on the human cortex is a challenging problem with important applications in brain mapping[1]. While their form and location can vary quite significantly, there is no difficulty for an anatomist to observe the regularity of major sulci based simply on the geometry of cortical surfaces regardless of their size, orientation, and the software used to extract them from MRI images. From an engineering point of view, this simplicity is critical for a computational system to achieve the same level of robustness, which essentially lies in its ability of modeling sulcal anatomy with the geometry of cortical surfaces. To this end, we propose in this work a novel system for automated sulci identification by integrating local geometry, i.e., curvature, with global descriptors derived from the Laplace-Beltrami eigen-system [2,3,4] and intrinsic surface mapping [5].

At the core of our system is an intrinsic atlas of major sulci that represents prior knowledge in training data and enables the integration of curvature information into major sulci. This atlas-based approach is most related to learning-based methods for sulci identification in previous work [6,7,8,9,10]. Principal component analysis were used to model sulcal basins [6] and sulcal lines on the sphere [7]. Boosting methods were used in [9,10]. The Markovian relation of multiple sulci was incorporated with graphical models [8,10]. The main novelty in our system is that the modeling of cortical anatomy is based entirely on intrinsic geometry, which eliminates the use of coordinates in conventional MRI atlases such as the Talairach atlas. This makes our system robust to pose differences and variations across populations. We demonstrate this robustness by using an atlas constructed from a Chinese population to detect twelve major sulci on Caucasian brains of different age groups and surfaces extracted by different software tools.

The rest of the paper is organized as follows. In section 2, we describe the construction of the intrinsic atlas for the modeling of sulcal anatomy. The sulci identification algorithm based on this atlas is developed in section 3. Experimental results are presented in section 4. Finally conclusions are made in section 5.

2 Atlas Construction

As illustrated in Fig. 1, there are two main steps in the construction of the atlas. In the first step, a multiscale representation of the cortical surface is constructed. In the second step, surface maps at a selected scale are computed to bring manually traced sulcal lines to an atlas surface. We describe the two steps in detail next.

Fig. 1
Atlas construction

2.1 Multiscale Surface Representation

Let An external file that holds a picture, illustration, etc.
Object name is nihms242596ig1.jpg = (An external file that holds a picture, illustration, etc.
Object name is nihms242596ig2.jpg, An external file that holds a picture, illustration, etc.
Object name is nihms242596ig3.jpg) denote the triangular mesh representation of a cortical surface, where An external file that holds a picture, illustration, etc.
Object name is nihms242596ig2.jpg and An external file that holds a picture, illustration, etc.
Object name is nihms242596ig3.jpg are the set of vertices and triangles. We construct the multiscale representation of An external file that holds a picture, illustration, etc.
Object name is nihms242596ig1.jpg by using the eigen-system of its Laplace-Beltrami operator An external file that holds a picture, illustration, etc.
Object name is nihms242596ig4.jpg, which is defined as

ΔMf=λf.
(1)

The spectra of An external file that holds a picture, illustration, etc.
Object name is nihms242596ig4.jpg is discrete and we denote the eigenvalues as λ0λ1 ≤ · · · and the corresponding eigenfunctions as f0, f1, · · ·. To numerically compute the eigenfunctions, we use the finite element method and solve a generalized matrix eigenvalue problem[2,3,4].

Let X(·, 0) : An external file that holds a picture, illustration, etc.
Object name is nihms242596ig1.jpgR3 denote the coordinate function on An external file that holds a picture, illustration, etc.
Object name is nihms242596ig1.jpg, i.e., X(p, 0) = p for p [set membership] An external file that holds a picture, illustration, etc.
Object name is nihms242596ig1.jpg. Using the eigen-system, we can express the heat diffusion of the coordinate function as

X(p,t)=n=0eλntfn(p)Mfn(q)X(q,0)dM
(2)

By replacing the coordinates of the vertices on An external file that holds a picture, illustration, etc.
Object name is nihms242596ig1.jpg with X(·, t), we have a multiscale representation of An external file that holds a picture, illustration, etc.
Object name is nihms242596ig1.jpg. For numerical implementation, we approximate the diffusion with the first 300 eigenfunctions that are computed efficiently with the spectrum shift technique [11]. As an illustration, we show in Fig. 2 the multi-scale representation of a cortical surface at the scale t = 0, 10, 100, 1000. With the increase of the scale, we can see the surface exhibits more regularity that is common across population. However, sulcal landmarks on the original surface might be overly distorted in terms of length and angle if t is too large. In our work, we usually choose t = 100 as a tradeoff between surface regularity and landmark distortion.

Fig. 2
Multiscale representation of a cortical surface

2.2 Atlas Construction via Intrinsic Surface Mapping

To construct the atlas at a selected scale, we extend the intrinsic surface mapping technique developed for sub-cortical surfaces in [5] to cortical surfaces. Given a pair of surfaces M1t and M2t at a scale t, we compute two maps u1:M1tM2t and u2:M2tM1t by minimizing the following energy:

E=M1t[j=13αDj(ξ1jξ2j[composite function (small circle)]u1)2+αIC(Iu2[composite function (small circle)]u1)2+αH||Ju1||2]dM1t+M2t[j=13αDj(ξ2jξ1j[composite function (small circle)]u2)2+αIC(Iu1[composite function (small circle)]u2)2+αH||Ju2||2]dM2t
(3)

where ξij(i=1,2;j=1,2,3) are the feature functions defined on the surfaces that characterize their intrinsic geometry. In each integral, the first term penalizes the difference of feature functions, the second term encourages inverse consistency, and the third term uses the harmonic energy [12] for smoothness regularization, where Ju1 and Ju2 are the Jacobian of the maps. The regularization parameters αD, αIC, αH control the weight of different terms.

To model the sulcal anatomy intrinsically, we define three feature functions on cortical surfaces as shown in Fig. 3(a) by using the Reeb graph of the first and second eigen-function of its Laplace-Beltrami operator [5]. We can see these functions characterize the frontal-posterior, superior-inferior, and medial-lateral profile of the surface intrinsically. By projecting each point p[set membership]Mit to ( ξi1(p)ξi2(p)ξi3(p)), we construct an embedding of each surface in the feature space. As shown in Fig. 3(b), the two cortical surfaces align much better in the embedding space, which enables us to find a good initial map between them by simply looking for the nearest point on the other surface. Starting from the initial maps, we iteratively evolve them by solving a pair of PDEs on the surfaces in the gradient descent directions [5].

Fig. 3
(a) Intrinsic feature functions. (b) Embedding in the feature space.

By choosing one of the surface An external file that holds a picture, illustration, etc.
Object name is nihms242596ig5.jpg in the training set as the atlas surface, we compute the maps from all other surfaces in the training data to this atlas surface An external file that holds a picture, illustration, etc.
Object name is nihms242596ig5.jpg. All the manually traced sulcal lines can then be projected onto the atlas surface as shown in the third row of Fig. 1. Since both the multiscale representation and surface mapping are established via intrinsic geometry, we denote the collection of sulcal lines projected onto the atlas surface An external file that holds a picture, illustration, etc.
Object name is nihms242596ig5.jpg as the intrinsic atlas of major sulci.

3 Automated Identification of Major Sulci

Using the atlas of major sulci, we develop in this section an automated system for their identification on cortical surfaces. Given a new surface An external file that holds a picture, illustration, etc.
Object name is nihms242596ig1.jpg, the skeletal representation of its sulcal region is first computed based on the mean curvature of An external file that holds a picture, illustration, etc.
Object name is nihms242596ig1.jpg [13] as shown in Fig 4(a). The skeleton is decomposed into a set of branches B = {B1, B2, · · ·}, where each branch is a polyline on An external file that holds a picture, illustration, etc.
Object name is nihms242596ig1.jpg. The multiscale representation An external file that holds a picture, illustration, etc.
Object name is nihms242596ig6.jpg of An external file that holds a picture, illustration, etc.
Object name is nihms242596ig1.jpg and the map u : An external file that holds a picture, illustration, etc.
Object name is nihms242596ig6.jpgAn external file that holds a picture, illustration, etc.
Object name is nihms242596ig5.jpg is then computed. With the map u, the branches are projected onto the atlas surface and denoted as [B with circumflex] = {[B with circumflex]1, [B with circumflex]2, · · ·,}. Our goal is to extract anatomically consistent sulcal lines from these skeletal branches.

Fig. 4
(a) Skeletal branches of the sulcal region. (b) Candidate branches (black) together with the training curves of the superior frontal sulcus. (c) Sample paths (black) and the most likely path (red). (d) The detected curve on the original surface.

Let Si={Si1,Si2,(...),SiJ} denote the set of training curves for the i-th sulcus on the atlas surface. Given the prior model Si, the challenge is how to model the likelihood of skeletal branches in [B with circumflex]. The difficulty arises from the fact that B contains only partial observations of major sulci as they frequently cross gyral regions. To bridge the gap between the prior model of complete sulci and the partial observation in skeletal branches, we propose below a projection operator to model the likelihood of a branch on a major sulcus. Given a curve segment C on An external file that holds a picture, illustration, etc.
Object name is nihms242596ig5.jpg, its projection onto the training data Si is defined as:

PSi(C)={y(x)[set membership]Sij,x[set membership]C[mid ]||xy||=minz[set membership]Sij||xz||}
(4)

where Sij=argminjdH(C,Sij) and dH is the Hausdorff distance. Given a branch [B with circumflex]k and its projection PSi ([B with circumflex]k), we calculate the Hausdorff distance Dk = dH ([B with circumflex]k, PSi ([B with circumflex]k), and the projection ratio Rk = ||PSi ([B with circumflex]k)||/||[B with circumflex]k||, where ||·|| denotes the length of a curve. Using the distance Dk and projection ratio Rk, we can pick a set of candidate branches as N = {[B with circumflex]k [set membership] [B with circumflex]| DkT HD1, RkT HD2} where T HD1 and T HD2 are thresholds which we choose as 15mm and 0.5 in our practice. The candidate branches are plotted together with the training set over the atlas surface in Fig. 4(b).

Using the candidate branches in N, we construct a directed graph to generate a set of sample paths on An external file that holds a picture, illustration, etc.
Object name is nihms242596ig5.jpg as candidate curves for the sulcus. Given two nodes [B with circumflex]p and [B with circumflex]q in N, whose points are ordered according to the indices of their projection on Si, we form a new curve Cp,q = ([B with circumflex]p, [B with circumflex]q) by connecting the end point B^pe of [B with circumflex]p with the start point B^qs of [B with circumflex]q. Once again we compute the projection of Cp,q onto Si to evaluate the likelihood of Bp and Bq belonging to the same sulcus. Let Dp,q and Rp,q denote the distance and projection ratio of Cp,q. If Rp,qT HD2, we add an edge from [B with circumflex]p to [B with circumflex]q and define the weight as 1Dp,q||B^peB^qs||, where the distance ||B^peB^qs|| is included to encourage the connection of close branches. Starting from any branch without parents, we perform random walks on the graph to generate sample paths on the atlas surface. The probability of taking an edge during any walk is in proportion to its weight.

Let Ci={Ci1,Ci2,(...)} be the set of sample paths for the i-th sulcus. For a candidate curve, its distance to the training data is defined as:

d(Cik;Si)=minj(d¯(Cik,Sij)+d¯(Sij,Cik))
(5)

where d(·, ·) is the average distance from points on a curve to the other curve. We also define the “sulcality” of each curve as:

Sulcality(Cik)=B^p[subset or is implied by]Cik||B^p||/||Cik||
(6)

which measures how good the path follows the sulcal regions. We then define the likelihood of each curve as

L(Cik[mid ]Si)=ed(Cik;Si)/Sulcality(Cik)
(7)

and choose the detection result as the sample curve in Ci with the maximum likelihood. As an illustration, we show on An external file that holds a picture, illustration, etc.
Object name is nihms242596ig5.jpg the set of sample paths and the path with the maximum likelihood for the superior frontal sulcus in Fig. 4(c). The branches in this path are then connected via a curvature-weighted geodesic on the original surface to obtain the detected sulcus as plotted in Fig. 4(d).

Let Ci1 and Ci2 denote the candidate curves of the i1-th and i2-th sulcus, we define the joint distance from a pair of curves Ci1p[set membership]Ci1 and Ci2q[set membership]Ci2 to the training curves Si1 and Si2 as:

d(Ci1p,Ci2q;Si1,Si2)=minj(d¯(Ci1p,Si1j)+d¯(Si1j,Ci1p)+d¯(Ci2q,Si2j)+d¯(Si2j,Ci2q))
(8)

The joint likelihood of the two curves is then defined as:

L(Ci1p,Ci2q[mid ]Si1,Si2)=ed(Ci1p,Ci2q;Si1,Si2)/(Sulcality(Ci1p)*Sulcality(Ci2q))
(9)

The joint detection results are the pair of curves in the sample space Ci1 × Ci2 achieving the maximum likelihood.

4 Experimental Results

In this section, we present experimental results on four different datasets, including surfaces extracted from 3 software tools, to demonstrate our sulci detection system.

4.1 Atlas Construction and Quantitative Validation

In the first experiment, we applied our method to pial surfaces extracted from the 3T MRI images of 65 Chinese young subjects (18 ~ 27 years) by Freesurfer [14]. The left hemispheres were used in this work. Twelve major sulci as listed in Table 1 were manually traced on each surface. We used 40 surfaces as training data to construct the intrinsic atlas of major sulci, which is shown on the third row of Fig. 1, and the algorithm developed in section 3 was used to identify the twelve major sulci on the other 25 surfaces for testing and quantitative validation. For robustness, the central and post-central sulcus were detected jointly by maximizing the joint likelihood in (9). The other 10 sulci were detected separately by maximizing the likelihood in (7).

Table 1
Quantile statistics in testing data. (S1:central; S2:post-central; S3:pre-central; S4:superior-temporal; S5:intraparietal; S6:inferior-frontal; S7:superior-frontal; S8:olfactory; S9:collateral; S10:parietal-occipital; S11:cingulate; S12:calcarine.)

As an illustration, the results from 4 subjects in the testing data are plotted in Fig. 5. For better visualization, we plotted the detected curves on the surfaces at the scale t = 10 in Fig. 5(a) and (b). As shown in Fig. 5(c), the automatically detected sulcal lines align very well with manually traced sulci plotted in black. Quantile statistics of two distances dam and dma of each sulcus were listed in Table 1, where dam is the distance from points on a detected curve to the corresponding manually traced curve, and dma denotes the distance from points on each manual curve to the detected sulcus. For example, the first number in the column of S1 means that 70% of the points on the automatically identified central sulcus are within a distance of 3.2mm to the manually traced curve. While there is variability across different sulci, we can see the quantile statistics show the automated results accurately capture the main body of the major sulci.

Fig. 5
Detection results on testing data. (a) Lateral view. (b) Medial view. (c) Overlay with manually traced curves (black).

4.2 Robustness

In the second experiment, we applied the atlas built from the Chinese population in the first experiment to detect sulci on three different datasets of Caucasian brains. The first dataset consists of the left pial surfaces of eight elderly subjects (63 ~ 85 years) extracted by Freesurfer. The second dataset is composed of the left pial surfaces of two young adults extracted by BrainSuite [15]. The third dataset consists of white matter surfaces of two young adults extracted by BrainVisa [16]. For better visualization, we also plot the detected curves on all surfaces at the scale t = 10 in Fig. 6. These results demonstrate the robustness of our method across ethnic and age groups, and different software tools for surface extraction.

Fig. 6
Row 1–4: results on elderly subjects. Row 5: results on BrainSuite surfaces. Row 6: results on BrainVisa surfaces.

5 Conclusion

In summary, we have developed a novel system for the automated detection of major sulci on cortical surfaces. Quantitative validations on twelve major sulci showed excellent agreement between our system and manual tracing. We also demonstrated the robustness of our system across different populations and surface extraction tools. For future work, we will augment our system with the intrinsic modeling of gyral landmarks to further improve its performance.

References

1. Thompson PM, et al. Mapping cortical change in alzheimers disease, brain development, and schizophrenia. NeuroImage. 2004;23:S2–S18. [PubMed]
2. 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]
3. Reuter M, Wolter F, Peinecke N. Laplace-Beltrami spectra as Shape-DNA of surfaces and solids. Computer-Aided Design. 2006;38:342–366.
4. Shi Y, et al. Anisotropic Laplace-Beltrami eigenmaps: Bridging Reeb graphs and skeletons. Proc MMBIA. 2008:1–7. [PMC free article] [PubMed]
5. Shi Y, et al. Inverse-consistent surface mapping with Laplace-Beltrami eigen-features. Proc IPMI. 2009:467–478. [PMC free article] [PubMed]
6. Lohmann G, Cramon D. Automatic labelling of the human cortical surface using sulcal basins. Med Image Anal. 2000;4:179–188. [PubMed]
7. Tao X, Prince J, Davatzikos C. Using a statistical shape model to extract sulcal curves on the outer cortex of the human brain. IEEE Trans Med Imag. 2002;21(5):513–524. [PubMed]
8. Rivière D, et al. Automatic recognition of cortical sulci of the human brain using a congregation of neural networks. Med Image Anal. 2002;6:77–92. [PubMed]
9. Tu Z, et al. Automated extraction of the cortical sulci based on a supervised learning approach. IEEE Trans Med Imag. 2007;26:541–552. [PubMed]
10. Shi Y, et al. Joint sulcal detection on cortical surfaces with graphical models and boosted priors. IEEE Trans Med Imag. 2009;28(3):361–373. [PMC free article] [PubMed]
11. Vallet B, Lèvy B. Spectral geometry processing with manifold harmonics. Computer Graphics Forum. 2008;27(2):251–260.
12. Mèmoli F, Sapiro G, Osher S. Solving variational problems and partial differ-ential equations mapping into general target manifolds. Journal of Computational Physics. 2004;195(1):263–292.
13. Shi Y, et al. Hamilton-Jacobi skeleton on cortical surfaces. IEEE Trans Med Imag. 2008;27(5):664–673. [PMC free article] [PubMed]
14. Dale AM, Fischl B, Sereno MI. Cortical surface-based analysis i: segmentation and surface reconstruction. NeuroImage. 1999;9:179–194. [PubMed]
15. Shattuck D, Leahy R. BrainSuite: An automated cortical surface identification tool. Med Image Anal. 2002;8(2):129–142. [PubMed]
16. Mangin JF, et al. From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations. Journal of Mathematical Imaging and Vision. 1995;5(4):297–318.