## Abstract

We present a technique for automatically assigning a neuroanatomical label to each location on a cortical surface model based on probabilistic information estimated from a manually labeled training set. This procedure incorporates both geometric information derived from the cortical model, and neuroanatomical convention, as found in the training set. The result is a complete labeling of cortical sulci and gyri. Examples are given from two different training sets generated using different neuroanatomical conventions, illustrating the flexibility of the algorithm. The technique is shown to be comparable in accuracy to manual labeling.

## Introduction

Techniques for labeling geometric features of the cerebral cortex are useful for analyzing a variety of functional and structural neuroimaging data (Rademacher *et al.*, 1992; Caviness *et al.*, 1996; Paus *et al.*, 1996; Seidman *et al.*, 1997; Goldstein *et al.*, 1999, 2001, 2002; Lohmann *et al.*, 1999). In many cases, what is desired is not a labeling of a discrete set of cortical features, but rather the identification of every point in the entire cortex. This type of labeling is known as a *parcellation*. Unfortunately, despite their potential utility, cortical parcellations are not commonly used in the neuroimaging community due to the difficult and time-consuming nature of the task of manually parcellating the entire cortex from high-resolution MRI images. While many techniques exist for labeling parts of the cortex (Mangin *et al.*, 1995; Thompson *et al.*, 1996b; Valiant *et al.*, 1996; Sandor and Leahy, 1997; Lohmann, 1998; Lohmann and von Cramon, 1998, 2000; Goualher *et al.*, 1999; Rettmann *et al.*, 2002), none of them provides an automated procedure for generating a complete and detailed labeling of the entire cortex that can incorporate prior information as well as cortical geometry into the labeling.

The inclusion of prior information is a critical feature of a cortical parcellation algorithm. The reason why this is the case is that the divisions that are useful in a cortical parcellation scheme can be based on properties of the brain other than cortical geometry. Many disparate pieces of information, such as knowledge of structure–function relationship, and cytoarchitectonic or receptor labeling properties of regions, may be used by a neuroanatomist in generating a cortical parcellation — information that is not directly available to an automated parcellation procedure from magnetic resonance imaging (MRI) data. A trivial example of this is the fact that many unbroken sulci change names as lobar boundaries are crossed. In addition, functional heterogeneity can make it desirable for adjacent areas to be assigned different labels despite the absence of macroscopic cortical features to distinguish them. Conversely, there are situations in which the secondary folding structure of the cortex (i.e. folds within folds) carries information, such as the location of the hand area in primary motor cortex, which is well predicted by a posterior-pointing secondary fold in the precentral gyrus (Boling *et al.*, 1999).

Various approaches have been taken to the problem of labeling of cortical features. For example, Sandor and Leahy (1997) use a manually labeled atlas brain, which is then warped into correspondence with an individual subject’s surface model. The individual subject surface then inherits the labels from the atlas. Similar surface-fitting approaches have been developed by Thompson *et al.* (1996b), Valiant *et al.* (1996), Lohmann (1998) and Lohmann and von Cramon (1998, 2000). A graph-based approach is taken in Goualher *et al.* (1999), in which sulci are represented by vertices in the graph and arcs connecting them represent their relationship to one another. A manual training set is then used in order to generate a semi-automated procedure for assigning neuroanatomical labels to the detected sulci. Another graph-based approach was taken in Mangin *et al.* (1995), in which a 3-D skeletonization of cortical sulci that is robust to topological errors is generated and used to construct a graph for sulcal labeling. A somewhat different technique was developed by Lohmann (1998). In this approach, sulci are not represented as abstract quantities such as vertices in a graph, but instead are derived using a watershed transform on an MR volume. This work was extended by Rettmann *et al.* (2002) to the use of surface models. The segmented regions identified by the watershed transform are then manually labeled by a trained neuroanatomist.

Here we present a technique for using manually labeled data as the basis for an automated parcellation procedure, using an intrinsically cortical coordinate system to store prior statistics and class-conditional densities (Fischl *et al.*, 1999a,b). The use of the manual parcellation as a training set allows the technique to incorporate neuroanatomical convention into the parcellation in regions in which geometry alone is not predictive of a parcellation label. The procedure models the parcellation labels as a first order anisotropic non-stationary Markov random field (MRF), allowing it to capture the spatial relationships between parcellation units that are present in the training set. The anisotropy separates label probabilities in the first and second principal curvature directions, in order to encode the increased probability of differing labels in the direction of maximal curvature. This type of modeling allows the procedure to automatically encode information such as ‘precentral gyrus frequently neighbors central sulcus in the direction of high curvature (moving across the sulcus), but not in the direction of low curvature (along the banks of the sulcus)’. This type of parcellation has commonly been generated by a having a trained anatomist or technician manually label some or all of the structures in the cortex, a procedure that can take up to a week for high-resolution images. Here, we use the results of the manual labeling using the validated techniques of the Center for Morphometric Analysis (CMA) (Rademacher *et al.*, 1992; Caviness *et al.*, 1996) to automatically extract the information required for automating the parcellation procedure. In addition, in order to demonstrate that the technique is sufficiently flexible to encode a different set of neuroanatomical conventions, we show comparable results using an entirely different surface-based (SB) parcellation scheme (Destrieux *et al.*, 1998).

## Materials and Methods

Local spatial relationships between labeled structures have been encoded by modeling the labeled image using *Markov random fields* (MRFs) in a variety of image processing contexts, dating back to their introduction (Geman and Geman, 1984). In the MRF approach, the probability of a label at a given position is computed not just in terms of the data and prior probabilities at that position, but also as a function of the labels in a neighborhood around the point in question. In the context of segmenting MR images, *isotropic* (all directions are equal) and *stationary* (the probabilities are the same for all spatial locations) MRFs have been used to provide smoothness constraints on a given segmentation (Held *et al.*, 1997; Kapur *et al.*, 1998; Zhang *et al.*, 2001). In this way, the prior probability of a label is computed by examining how likely the label is given the labels of its neighbors, regardless of the direction of the neighbor, or the position within the brain. While this type of approach can obviate the need for prefiltering of the images, it does not provide for the use of information regarding the spatial relationships of classes to one another. Anisotropic MRFs have been used in the context of segmenting MRI volumes in order to disambiguate the problem of labeling subcortical structures (Fischl *et al.*, 2002), work which is extended here to the labeling of the cortical surface. [Note that the derivations in the following sections closely follow those in Fischl *et al.* (2002).]

Thus, we base the cortical parcellation procedure on a number of pieces of information incorporated into a space-varying classification procedure. That is, class statistics (e.g. means and covariance matrices) are tabulated regionally throughout an atlas space, using a registration procedure that optimally aligns cortical folding patterns (Fischl *et al.*, 1999b). Prior probabilities are computed via a frequency histogram in the atlas space, allowing the calculation of the probability that each parcellation label occurs at every atlas location. Finally, the prior probability of a given spatial arrangement of parcellation labels is incorporated into the final parcellation procedure. These priors are also computed from a training set, examples of which are given in Figure 1, for each point in the atlas by modeling the parcellation as an anisotropic non-stationary Markov random field, resulting in a procedure that is comparable in terms of accuracy to manual labeling.

### Problem Statement

The problem of automatically labeling cortical gyri and sulci can be naturally phrased within the framework of Bayesian parameter estimation theory. In this approach, one can relate the probability of a parcellation *P* given the observed surface model *S* to the probability *p*(*S*|*P*) of the geometry of the surface model occurring given a certain parcellation, together with the prior probability of the segmentation:

The primary advantage of the Bayesian approach is that it allows for the explicit incorporation of prior information via the *p*(*P*) term in equation (1). In order to render the problem more tractable in the face of the large degree of variability in cortical folding patterns, both the priors on *P* and the conditional probability of observing the surface given the classification *p*(*S*|*P*) can be expressed within an atlas space, allowing them to vary as a function of position on the cortical surface (hence making them non-stationary). The advantage of using an atlas space is that coordinates in the atlas have more anatomical meaning than the native coordinate system of the image (Bajcsy *et al.*, 1983; Bookstein, 1989; Miller *et al.*, 1993; Gee *et al.*, 1994; Vannier *et al.*, 1994; Christensen *et al.*, 1996; Collins *et al.*, 1996; Ashburner *et al.*, 1997; Collins and Evans, 1997; Woods *et al.*, 1998; Fischl *et al.*, 1999b; Thompson *et al.*, 2000). Classifiers can then be distributed throughout the atlas, allowing each one to focus on only the small number of classes that may occur within the region for which the classifier is responsible. The number of classes that occur within a region of space is then directly related to the accuracy of the atlas coordinate system. That is, *P*(*P*(r) = *c*) will be 0 for all but a few values of *c* at each atlas location r. In practice, if the classifiers are reasonably dense in the atlas space then the number of classes at each location is typically relatively small. In this way, the intractable problem of classifying each surface node into one of 80 or so labels is decomposed into a set of tractable problems of classifying the nodes in each region of the surface model into only a small number of labels.

The definition of the atlas requires the calculation of a function *f*(r), which takes native image coordinate as input, and returns the coordinate of the corresponding point in the atlas. For *f* to be useful in this context, the coordinates it returns should be related to the anatomical location of **r**. This type of mapping therefore provides the ability to meaningfully relate coordinates across subjects. In the most general case, we wish to maximize the joint probability of both the parcellation *P* and the atlas function *f*:

The terms *p*(*S*|*P*,*f*) and *p*(*P*|*f*) in equation (2) provide a natural means for incorporating atlas information into the parcellation procedure. The first term encodes the relationship between the class label at each atlas location and the predicted surface geometry. Using the atlas space, we can allow the class statistics to vary as a function of location, allowing the within-class variations in surface geometry, such as secondary and tertiary cortical folds, to be incorporated into the classification procedure. The second term allows the expression of prior information regarding the spatial structure of the parcellation labels. Finally, the term *p*(*f*) provides a means for constraining the space of allowable atlas functions (e.g. continuity, differentiability and invertibility). In the following, we will use the atlas function *f* described in Fischl *et al.* (1999b), which establishes a spherical surface-based coordinate system for each cortical hemisphere.

### Atlas Construction

In general, two different approaches have been taken to the construction of anatomical atlases from neuroimaging data. The first is to use an individual as a template, and either manually or automatically estimate a transformation that aligns each new dataset with the individual template (e.g. Talairach *et al.*, 1967; Talairach and Tournoux, 1988; Van Essen and Drury, 1997; Van Essen *et al.*, 1998). The alternative is to compile a probabilistic atlas based on the anatomy of a large number of subjects (e.g. Collins *et al.*, 1994; Fox *et al.*, 1994; Mazziotta *et al.*, 1995; Thompson *et al.*, 1997, 2000). Each of these approaches has strengths and weaknesses. The former allows one to represent anatomical structures at as fine a scale as the neuroimaging technology allows, but the atlas is then biased by the idiosyncrasies of the individual anatomy chosen as the template. The latter technique resolves this problem by averaging across anatomies, thus only retaining that which is common in the majority of subjects. Nevertheless, the cross-subject averaging removes potentially useful information in the atlas. Here, we preserve the advantages of each technique by using a group of subjects to construct an atlas that retains information about each parcellation unit at each point in space. Given the atlas function *f*, and a group of *N* manually parcellated subjects, we first estimate the prior probability of parcellation label *c* occurring at each atlas location, independent of all other locations:

The geometry of the observed surface *S* at each point is summarized by two quantities encoded as a vector **G**(**r**) at each point in the surface (note that here **r** refers to spherical atlas coordinates, and *f*(**r**) refers to individual subject spherical coordinates). The first of these is the average convexity as defined in Fischl *et al.* (1999a), while the second is the mean curvature of the surface (do Carmo, 1976). Additional quantities such as MRI derived variables (e.g. *T*_{1} and *T*_{2}), Gaussian curvature and cortical thickness can easily be incorporated into this framework. The likelihood of observing the surface geometry **G** given the parcellation label *P*(**r**) is modeled as a Gaussian, the parameters of which are computed from the manually labeled training set:

where **G*** _{i}* are geometric parameters extracted from each of the set of

*M*surfaces for which label

*c*occurs at location

*f*(

**r**) in the corresponding manually labeled surface models

*S*

*[i.e.*

_{i}*S*

*(*

_{i}*f*(

**r**)) =

*c*]. The covariance matrix for class

*c*at location

**r**in the atlas is then given by:

Information about surface geometry is thus maintained separately for each parcellation label at every location in the atlas, obviating the need to average geometric information across classes. Finally, we also estimate the pairwise probability that parcellation label *c*_{2} is the neighbor at r* _{i}* when parcellation label

*c*

_{1}is the label at

**r**, for r

*ϵ*

_{i}*N*(

**r**), where

*N*(r) is a neighborhood function of r.

As before, this information is stored separately for each atlas location. It is important to note that the probabilities are stored separately for each pair of classes as well as for each neighborhood location r* _{i}*. While it may seem that this would lead to combinatorial explosion and intractable memory requirements, in practice the space is sparse as relatively few configurations of parcellation labels actually occur. Two examples of the manual labeling that is the basis for the atlas are given in Figure 1. The spherical atlas space is made up of a pair of super-sampled icosahedra. The priors are stored at essentially the same resolution as the individual surface models, using an icosahedral model with 163 842 vertices at ∼1 mm spacing, while the class conditional densities are stored 2562 vertices at ∼4 mm spacing, reflecting the slower variation of surface geometry for each label.

### Parcellation

As mentioned above, the Bayesian approach allows us to incorporate prior information that is necessary for the parcellation procedure. This prior information takes two forms. The first makes use of the global spatial information provided by *f* and the atlas in order to express the probability that a given parcellation label occurs at a particular location in the atlas, independent of the local surface geometry. The second encodes the local spatial relationship between parcellation labels, allowing information such ‘precentral gyrus can be anterior, superior or inferior to central sulcus, but never posterior to it’ to be automatically detected and incorporated into the segmentation.

Formally, we compute the maximum *a posteriori* (MAP) estimate of the parcellation *P* given an input surface geometry **G**, and the nonlinear spherical transform *f*. The MAP estimate can be expressed as maximizing *p*(*P|***G**,*f*), the probability distribution of the parcellation given the observed surface geometry and the atlas function*.* Using Bayes rule we can relate this to the product of the probability of observing the surface geometry with the prior probability of a given spatial configuration of labels *p*(*P*):

Assuming the noise at each vertex is independent from all other vertices in the surface model, we can rewrite *p*(**G***|P*,*f*) as the product of the distribution at each vertex over the surface:

Note that in the more general case in which the spatial correlation structure of the noise is constant across space, the equality in equation (8) should be replaced by proportionality (Worsley *et al.*, 1992; Thompson *et al.*, 1996a), which does not effect any of the subsequent derivations. The distribution of the geometric properties of each class at each location in the atlas is modeled as a Gaussian, the mean vector µ* _{c}*(

**r**) and covariance matrix Σ

*(*

_{c}**) of which are computed using equations 4 and 5. The probability of observing the geometry at**

**r****G**(

*f*(

**r**)) is then expressed as:

All that remains is to find an expression for the prior probability of a given parcellation *P*. Here we assume that the spatial distribution of labels can be well approximated by an anisotropic non-stationary Markov random field. This allows one to encode prior information about the relationship between labels as a function of location within the cortex (i.e. non-stationary), as well as with local direction (i.e. anisotropic). Formally, the Markov assumption can be expressed as:

That is, the prior probability of a label at a given vertex **r** is only influenced by the labels within some neighborhood of ** r**. The locality restriction imposed by the Markov model permits the probability of the entire parcellation to be written in terms of neighborhood or

*clique*potentials

*V*

*(*

_{c}*P*) via the Hammersley–Clifford theorem (Besag, 1974). That is, the probability

*p*(

*P*) can be equivalently characterized by a Gibbs distribution:

where *Z* is a normalizing constant and will be dropped in the following, and *U*(*P*) is an energy function that can be written in the form:

The clique potentials *V** _{c}*(

*P*) encode the energy associated with a certain configuration of labels within the

*c*th clique. Choosing

*V*

*(*

_{c}*P*) to be

*–*log(

*p*(

*P*(

**)**

**r***|P*(

**r**_{1}),

*P*(

**r**

_{2})

*… P*(

**r**

*)), where*

_{K}**is the central vertex of the**

**r***c*th clique, allows us write the probability of the entire parcellation as the product of the probability of the label at each vertex, given its neighborhood:

Using Bayes rule, we can rewrite this as:

Equation (14) allows the probability of a given label to be modulated by any configuration of neighboring labels. While this would be extremely useful, it is unfortunately not computationally tractable to implement, as one would need to compute separate prior probabilities for every combination of neighboring labels that occur. Instead, we make the simplifying assumption that only the first order conditional dependence is important. That is, that the dependence of a label on its neighbors can be expressed as the product of the probability given each of the neighbors:

where again we have explicitly included the dependence on the neighbor location **r*** _{i}* to emphasize that the probability densities are maintained separately for each neighbor position in

*N*(

**r**). Using this assumption, we arrive at an expression for the prior probability of the full parcellation:

Equation (16) allows two types of prior information to be incorporated into the segmentation procedure. The approximate location a surface feature may occupy within the cortex is given by *p*(*P*(**r**)), which is computed and stored in the atlas using equation (3). The local relationship between parcellation labels is encoded in *p*(*P*(**r*** _{i}*)

*|P*(

**r**),

**r**

*) using equation (6). We currently let the neighborhood function*

_{i}*N*(

**r**) include the two principal directions on the surface at each location in the atlas space. This allows the parcellation procedure to store information separately for the direction of highest curvature from that of the direction of lowest curvature, as the majority of borders occur in the former.

Directly computing the global MAP estimate of *P* in equation (7) using the Markov model of equation (16) is computationally intractable. Instead, we employ the iterated conditional modes (ICM) algorithm proposed by Besag (Besag, 1986). In this approach, the parcellation is initialized with the MAP estimate assuming *p*(*P*(r* _{i}*)

*|P*(r),r

*) is uniform, as no label has yet been assigned to each vertex. The parcellation is then sequentially updated at each location by computing the label*

_{i}*P*(r), that maximizes the conditional posterior probability

*p*(

*P*(

**r**)

*|P*(r

*),*

_{i}**G**,r

*):*

_{i}Equation (17) is then iteratively applied until no vertices are changed. Snapshots of the evolution of this procedure are given in Figure 2. As can be seen convergence is rapid, resulting in a procedure that requires ∼3 min on a 1.5 GHz Pentium III. A small amount of and post-processing is carried out on the results of the labeling to relabel small isolated patches of label to the most likely of the neighboring labels.

## Results and Discussion

Here we present the results of applying the surface-based parcellation technique to two different parcellation schemes. The first is the volumetric one in use at the MGH CMA as described in (Rademacher *et al.*, 1992; Caviness *et al.*, 1996), which defines ∼58 separate labels (see Appendix for a complete list), while the second is an intrinsically surface-based (SB) parcellation (Destrieux *et al.*, 1998) based on the conventions established in Duvernoy (1991) with 85 separate parcellation units (see Appendix). For the CMA parcellation we first sample the volumetric labeling onto the reconstructed cortical surface of each subject (Dale *et al.*, 1999; Fischl *et al.*, 1999a, 2001; Fischl and Dale, 2000), the subsequent procedures for the two parcellations are identical.

The data used for the CMA parcellation are part of an ongoing study of schizophrenia (Seidman *et al.*, 1997; Goldstein *et al.*, 1999; Seidman *et al.*, 1999; Goldstein *et al.*, 2002). As part of this study, 36 MRI volumes (2 MP-RAGE scans per subject, motion corrected and averaged) have been manually parcellated by trained technicians. (Note that as the study is not yet complete, we are blinded to diagnosis.) Cortical models were reconstructed for each of the subjects using previously presented techniques, including non-rigid surface-based alignment to a previously constructed spherical atlas (Dale and Sereno, 1993; Dale *et al.*, 1999; Fischl *et al.*, 1999a,b; Fischl and Dale, 2000; Fischl *et al.*, 2001). The manual parcellations were then sampled onto the cortical models at the midpoint of the cortical ribbon using the thickness estimates described in Fischl and Dale (2000). For the SB parcellation, cortical surface models were reconstructed for 12 subjects. Surface-based drawing tools that are part of the FreeSurfer software package (http://surfer.nmr.mgh.harvard.edu) were then used to generate a complete parcellation for all 24 cortical hemispheres. In order to validate the accuracy of the automated parcellation algorithm, we employed a jackknife/leave-one-out technique in which one dataset is removed from each training set, an atlas is constructed, and used to parcellate the left-out set. The error is then quantified point-by-point as the percentage of points at which the manual and automated parcellation are in disagreement. Figure 3 presents the results of this study on an average inflated surface (dark gray are sulcal regions, and light gray are gyral ones). Note that the vast majority of the surface is below threshold (i.e. accuracy >75%) for both the CMA parcellation on the left as well as the SB parcellation on the right. With the median accuracy at ∼81% for the left hemisphere and 83% for the right for the CMA parcellation, and 80% (left) and 79% (right) for the SB (the slightly lower accuracy for the SB parcellation is no doubt due to the greater number of parcellation units in it, some of which are quite small). In addition, the low-accuracy regions are clustered at the boundaries between labels, the precise positions of which are somewhat arbitrary in the manual training set.

The distribution of accuracy collapsed across each cortical hemisphere is given in Figure 4 for the CMA parcellation and Figure 5 for the SB parcellation. The top rows display histograms of the percent correct for the left and right hemispheres, while the bottom row is the cumulative histogram of the same data. As can be seen, almost 40% of the surface was correct 100% of the time for both hemispheres for both of the parcellation schemes.

A scatter plot of the accuracy of the individual subjects is given in Figure 6 CMA parcellation (left) and the SB parcellation (right). Interestingly, a single outlier can be seen in the left hemisphere of both parcellation sets. In both cases this is a subject with a rare folding pattern — a split central sulcus, and even in these case, the accuracy is close to 70%. As our database of manually labeled examples grows, we expect these types of rare neuroanatomies to be handled with better accuracy (the training set that was used for this subjects actually contained *no* examples of this type of folding pattern, as it was of course left out during construction of the classifiers for labeling it).

In order to directly compare the automated and manual parcellations we computed the mean and standard errors of the surface area of a set of cortical features across all the subjects in each dataset, using the two techniques for each of the parcellation schemes. This comparison is given in Figure 7 for the CMA parcellation, and Figure 8 for the SB parcellation. As can be seen, the surface area for a variety of cortical regions as assessed by the manual (light bars) and automated (dark bars) parcellations is essentially indistinguishable from a statistical standpoint, although it is worth pointing out the uniformly smaller error bars of the automated technique, indicating it is a potentially more powerful technique for detecting variations in cortical surface area.

Finally, a visual comparison of the automated and manual parcellation for the subject with the median accuracy in each the dataset is given in Figure 9. Note that many of the small protuberances in the manual CMA parcellation (on the left) are almost certainly errors in the manual labeling. We are currently developing an intrinsically surface-based adaptation of this parcellation that will minimize this type of error. In addition, the automated parcellation can be used to correct errors in the manual labelings, although of course these datasets could not be used for validation of the automated technique after correction.

## Conclusion

Labeling of cortical features is a technique that has been shown to have many neuroscience applications dating back to classic studies on functional lateralization (Geschwind and Levitsky, 1968). More recently, Schlaug and colleagues have found differences in musicians and non-musicians in perirolandic cortex (Schlaug *et al.*, 1995), while Boling *et al.* have shown folding patterns to be predictive of the location of the motor representation of the hand (Boling *et al.*, 1999). Witelson and Kigar demonstrated asymmetries in the Sylvian fissure that relate to handedness (Witelson and Kigar, 1992), while Goldstein and colleagues carried out an extensive study of sexual dimorphism both cortically and subcortically (Goldstein *et al.*, 2001), as well as the disruption of normal sexual dimorphism in schizophrenia (Goldstein *et al.*, 2002). Other applications include assessments of normal anatomical asymmetry in the sylvian fissure (Rubens *et al.*, 1976), planum temporale changes in neurofibromatosis (Billingsley *et al.*, 2002) and autism (Rojas *et al.*, 2002). Schizophrenia researchers have found effects in many cortical regions including the planum temporale (Barta *et al.*, 1995; Hirayasu *et al.*, 2000), the angular gyrus (Niznikiewicz *et al.*, 2000), the hippocampus (Seidman *et al.*, 1997), as well as the middle frontal gyrus, orbitofrontal cortex, the anterior and paracingulate gyri and the supramarginal gyrus (Goldstein *et al.*, 1999).

While manual methods exist for assessing this type of change, the process of manually parcellating an entire high-resolution structural MR volume requires on the order of a week for a trained neuroanatomist or technician. This makes the routine analysis of large patient and control populations untenable. Finally, manual labeling procedures do not generalize well to the use of multidimensional inputs such as cortical thickness measures, intrinsic MRI tissue parameters and cortical geometry, which can be easily incorporated into the automated parcellation procedure.

The automated method described in this paper for assigning a neuroanatomical label to every point in the cortex has been shown to be comparable in terms of accuracy to a previously validated method of manual parcellation. The accurate labeling of a large number of cortical structures is enabled through the use of both global and local spatial information. The global information is encoded by distributing classifiers throughout an atlas in a surface-based coordinate system and maintaining class statistics on a per-class, per-location basis, allowing the classifiers to be robust to variations cortical folding patterns. Local information is incorporated into the classification procedure by modeling the parcellation as a non-stationary anisotropic Markov random field. The introduction of anisotropy and non-stationarity into the parcellation model allows the spatial relationships of parcellation units to one another to be incorporated into the procedure in a principled fashion.

The automatic parcellation procedure can also be used to automatically define regions of interest (ROIs) for use in functional imaging studies. Specifically, this will allow one to generate average time courses by cortical region, or even parts of regions, facilitating for example the comparison of the response of the motor cortex to that of somatosensory cortex, or anterior superior temporal sulcus (STS) to the posterior STS.

While the results presented here are comparable to those achievable by current manual labeling techniques, we intend to extend them in a number of ways. Primarily, we will investigate the relationship between the manual and automated labelings of each brain. In particular, the CMA parcellation technique is defined in the volume using slice-based viewers, and is currently being updated to use surface-based drawing and visualization tools. When complete, we expect the surface-based extension of the CMA parcellation to be significantly more accurate than the labeling based on the volumetric parcellation. As part of this process, we will perform structure-by-structure analysis of variability both across manual labelings, and between manual and automated labelings for each of the manual training sets. This analysis should allow us to identify regions of low reliability, which will require either updating the definition of the parcellation labels in the case that neither the automated nor the manual labelings is reliable, or the inclusion of additional information into the automated parcellation procedure. This may include the relationship of the cortical regions to subcortical structures, as well as additional cortical information (e.g. regional thickness measures, intrinsic tissue parameters).

The automated nature of the methods described here, in contrast with existing manual or semi-automated techniques, allows for their routine application in large-scale studies. Having access to this type of detailed morphometric information for large populations including various disorders as well as a spectrum of normal controls should facilitate the characterization of the anatomical signatures associated with specific disorders. Ultimately, this may provide a more accurate and sensitive tool for early diagnosis of brain disorders.

## Appendix

### Labels for CMA cortical parcellation

AG angular gyrus

BFsbcmp basal forebrain

CALC intracalcarine cortex

Cga cingulate cortex, anterior part

CGp cingulate cortex, posterior part

CN cuneus

CO central operculum

F1 superior frontal gyrus

F2 middle frontal gyrus

F3o inferior frontal gyrus, pars opercularis

F3t inferior frontal gyrus, pars triangularis

FMC frontomedial cortex

FO frontal operculum

FOC fronto-orbital cortex

FP frontal pole

H1 Heschl’s gyrus

INS insula

LG lingual gyrus

MW medial wall

OF occipital fusiform gyrus

Oli occipital lateral gyri, inferior part

Ols occipital lateral gyri, superior part

OP occipital pole

PAC paracingulate cortex

PCN precuneus

Pha parahippocampal gyrus, anterior part

PHp parahippocampal gyrus, posterior part

PO parietal operculum

POG postcentral gyrus

PP planum polare

PRG precentral gyrus

PT planum temporale

SC subcallosal cortex

SCLC supracalcarine cortex

Sga supramarginal gyrus, anterior part

SGp supramarginal gyrus, posterior part

SMC supplementary motor cortex or juxtapositional lobule cortex

SPL superior parietal lobule

T1a superior temporal gyrus, anterior part

T1p superior temporal gyrus, posterior part

T2a middle temporal gyrus, anterior part

T2p middle temporal gyrus, posterior part

T3a inferior temporal gyrus, anterior part

T3p inferior temporal gyrus, posterior part

Tfa temporal fusiform gyrus, anterior part

TFp temporal fusiform gyrus, posterior part

TO2 middle temporal gyrus, temporo-occipital part

TO3 inferior temporal gyrus, temporo-occipital part

TOF temporal fusiform gyrus, temporo-occipital part

TP temporal pole

UN unknown

### Labels for surface-based parcellation

Corpus_callosum

G_and_S_Insula_ONLY_AVERAGE

G_cingulate-Main_part

G_frontal_inf-Opercular_part

G_frontal_inf-Triangular_part

G_frontal_superior

G_insular_long

G_occipital_inferior

G_occipital_superior

G_occipit-temp_med-Lingual_part

G_orbital

G_parietal_inferior-Angular_part

G_parietal_superior

G_precentral

G_rectus

G_subcentral

G_temp_sup-Lateral_aspect

G_temp_sup-Planum_tempolale

G_temporal_middle

Lat_Fissure-ant_sgt-ramus_horizontal

Lat_Fissure-post_sgt

Pole_occipital

S_calcarine

S_central_insula

S_cingulate-Marginalis_part

S_circular_insula_inferior

S_collateral_transverse_ant

S_frontal_inferior

S_frontal_superior

S_intermedius_primus-Jensen

S_intraparietal–and_Parietal_transverse

S_occipital_inferior

S_occipital_superior_and_transversalis

S_occipito-temporal_medial_and_S_Lingual

S_orbital_medial-Or_olfactory

S_paracentral

S_pericallosal

S_precentral-Inferior-part

S_subcentral_ant

S_suborbital

S_supracingulate

S_temporal_superior

S_transverse_frontopolar

Unknown

This research is funded jointly by the National Center for Research Resources (R01 RR16594-01A1), the National Institute of Neurological Disorders and Stroke, the National Institute of Mental Health and the National Cancer Institute (R01-NS39581, R01-NS34189). Further support was provided by the National Center for Research Resources (P41-RR14075 and R01-RR13609), the National Institute of Mental Health (NIMH MH56956), the Theodore and Vada Stanley Foundation, the Amytotrophic Lateral Sclerosis Association, Faculté de Médecine/CHRU, Tours, France, the Mental Illness and Neuroscience Discovery (MIND) Institute and the National Association for Research in Schizophrenia and Depression. Thanks to Steve Hodge and Jennifer Koch for providing much-needed support and organization. Thanks also to Andrea Boehland, Bobbie Holland and Nicole Cullen for providing manual parcellations for use in this paper.

## References

*et al.*(

*in vivo*magnetic resonance imaging.

*et al*. (

*In vivo*morphometry of the intrasulcal gray matter in the human cingulate, paracingulate, and superior-rostral sulci: hemispheric asymmetries, gender differences and probability maps.

*In vivo*evidence of structural brain asymmetry in musicians.

*et al.*(

*et al.*(