Author manuscript, published in "MICCAI Workshop- Manifolds in Medical Imaging: Metrics, Learning and Beyond, New-York: États- Unis d'Amérique (2008)" The Labeling of Cortical Sulci using Multidimensional Scaling

The task of classifying or labeling cortical sulci is made difficult by the fact that individual sulci may not have unique distinguishing features and usually need to be identified by a multivariate feature set that takes the relative spatial arrangement into account. In this paper, classical multidimensional scaling (MDS), which gives a geometric interpretation to input dissimilarity data, is used to classify 180 sulci drawn from the ten major classes of sulci. Using a leave-one-out validation strategy, we acheive a success rate of 100% in the best case and 78% in the worst case. For these more difficult cases, we propose a second stage of classification using shape based features. One of these features is the geodesic distance between sulcal curves obtained from a new open curve representation in a geometric framework. With MDS, we oﬀer a simple and intuitive approach to a challenging problem. Not only can we easily separate left and right brain sulci, but we also narrow the classification problem from, in this case, a 10-class to a 2-class problem. More generally, we can identify a region-of-interest (ROI) within which one can carry out further classification.


Introduction
The large variability present in the human brain cortex makes it one of the most challenging of brain anatomies to segment, label, classify or model.The sulci, which are fissures or grooves in the cortical surface (Figure 1), vary not just in their shape and placement, but even in their number.Even so, the underlying architectonic and functional organization of the brain influences to some degree the size, shape and placement of the folding patterns [1].The length, depth, the pure shape (i.e. after removing scaling, rotation and translation), the spatial arrangement and orientation of some of these individual sulci are characteristic features shared across a population.Collectively, these features can form the basis for classification, labeling and the study of selective differentiation of tissue that results from neurodegenrative pathologies, the action of genes or even as a result of cognitive activity.This study fills an important need and can lead to such wide-ranging applications as landmark or guidance-based neurosurgical procedures or the differential diagnosis of degenerative brain disease.
Accurate sulcal identification can be a challenge even for expert neuroanatomists.A strategy that is used in manual annotation is to identify sulci in relation to the major sulci for which we can show consistent anatomical correspondence between subjects.The major sulci provide a context and coupling this information with shape attributes, which are usually not unique to a class of sulci, helps in identification [2].The automated and semi-automated labeling methods that incorporate this idea typically involve using a large feature set and computationally intensive algorithms [3].We propose a simpler approach, not considered before, which uses only distance or dissimilarity data to individually identify unlabeled sulci against a reference graph (Figure 1(b)).
Multidimensional scaling (MDS) is a technique that gives a geometric interpretation to dissimilarity data and is a natural tool for such an application.By minimizing a loss function calculated for different possible configurations, we can assign a set of coordinates to objects.The resulting map, or embedding, places objects that have similar attributes close to each other.MDS can be used as a classification tool by introducing unlabeled sulci to an existing map.Membership is then assigned to the nearest node using a minimum distance classifier.
There is a growing recognition, in medical image analysis, that to take advantage of the full range of information presented in an image, one has to consider shape attributes.Sulcal shape analysis has typically used a multivariate feature set to describe shape attributes such as angle, curvature, depth, length [4].To date, little work has been done using a formal definition of shape spaces of sulcal curves.In this paper, we use a new geometric representation for studying shapes of three-dimensional sulcal curves using metrics invariant to the standard shape-preserving transformation.
In Section 2 we briefly explain the MDS method and how it is used for labeling the major sulci.Next, we present a new shape representation for open curves in Section 3. Though the formalism itself has a wider scope, the focus in this paper is its application to sulcal curves.In Section 4, we describe the other features used in our dissimilarity calculations.Finally, in Section 5, we present the preliminary results.In the discussion that follows, Section 6, we outline our plans to develop this classification method further.

Metric Multidimensional Scaling
Multidimensional scaling is a method that allows us to compute an optimal configuration, X, for a set of n observations using only their interpoint distances.If ∆ = [δ ij ] is a non-negative, symmetric matrix of dissimilarities computed from a set of (hidden) coordinates, X r of rank r, with MDS, we can represent X r as X d , d ≤ r, such that the distance matrix D computed from it is close to ∆.An intermediate step would be to derive the scalar product matrix, B = XX T , from the squared dissimilarities by double centering.An eigen decomposition, B = Γ ΛΓ T , gives us Λ, a diagonal matrix of eigenvalues, and Γ, a matrix of orthonormal eigenvectors.X, the square root of XX T , is then Γ Λ 1 2 .Since B is of rank d, X is also of rank d.The optimal solution for X d is found by minimizing the stress function: (1)

Out-of-sample Embedding
For classification, our objective is to introduce k additional points, y 1 ,. . .,y k , without disturbing an existing representation.Using Eqn. 1 to embed n + k points, however, gives an entirely new map.Intead, if we optimize where A = [a ij ] is the augmented (n + k) × (n + k) dissimilarity matrix, we can get y i , the coordinates for the k points, while x i , the original points remain fixed.This treatment, detailed in Trosset et al. [5], is an extension of the optimization solved by Eqn. 1.

Elastic Shape Analysis of Sulcal Curves
We use a new square-root velocity description for shapes of three-dimensional sulcal curves.This representation is related to the velocity function used in [6] and to the complex square-root representation for planar curves in [7].The advantage of this representation, first proposed by Joshi et al. [8] for closed curves, is that the resulting pre-shape becomes a unit hypersphere space and one can use the geometry of the sphere for computations.

Curve Representation and Pre-Shape Space
Let β : [0, 1] → R 3 be a unit-length, parameterized curve such that its speed β(t) is non-zero everywhere.We define its square-root velocity function as: The space of all square-root velocity representations of curves in R 3 , B, is a unit hypersphere and submanifold of L 2 ([0, 1], R 3 ).Since B is a Riemannian manifold with the inner product metric w 1 , w 2 , we can derive an exponential and an inverse exponential map.We can also compute geodesic paths between curves using where θ = cos −1 ( q 1 , q 2 ) is the geodesic distance, denoted by d c (q 1 , q 2 ), between two points.

Shape Space
The elements of B do not represent the shape of a curve uniquely.Any reparametrization or rigid rotation of β, results in a different square-root velocity function, q, while preserving its shape.We will call B the pre-shape space of open curves.This variability can be represented using the actions of the reparameterization group, Γ , and the rotation group, SO(3), on B. These actions commute and both groups act on B in an isometry.We can thus define the shape space as: S = B/(Γ × SO(3).
Let [q] denote an orbit resulting from all possible re-parameterizations and rotations of a curve denoted by q ∈ B. S is the set of all such orbits and it inherits a Riemannian structure from B, under which, the geodesic distance between two orbits [q 1 ] and [q 2 ] is: The geodesic is given by [ψ t ], where ψ t is the geodesic in B between q 1 and γ * O * q 2 (γ * ).Here (O * , γ * ) are the optimal transformations of q 2 that minimize d c .A closer look at that distance function reveals: argmin γ∈Γ,O∈SO(3) This equality says that to minimize the arc length on a unit sphere is the same as minimizing the chord length.

Karcher Means of Open Curves
To get the mean shape for a given set of curves, we compute a Karcher mean, μn = argmin A gradient-based approach to find the minimum is commonly used for finding means on nonlinear manifolds.The extrinsic mean, µ → µ/ µ , is used to initialize the gradient search.We update µ by computing the exponential map, exp µ (ǫv), to get a local minimizer for the cost function in Eqn. 5.

Feature Set
For the MDS evaluation, we used four descriptors, position, depth, length and orientation, in addition to shape.These contain physical information not captured in the shape of a 1D curve alone.For each of these descriptors, our features were the distances between sulci or the class averages of sulci.

Spatial Position
The mid-point of each sulcal curve was used for the position coordinates.
Mean Depth The depth profile for each sulcus was first obtained by resampling 100 points along the length of the sulcus and computing the distance between points on the top sulcal curve with corresponding points on the bottom curve.The mean depth was then calculated by averaging these.
Length The length of a sulcus was computed using 1 0 β(t) dt, the length of the curve β.Orientation Using singular value decomposition (SVD), the relative orientations, O, between sulcal curves were computed for each subject.
The class average for each of these descriptors was calculated by averaging across the 18 subjects.For the relative orientation, the average was obtained by computing a SVD of i O i .The Karcher mean gave us the mean shape.
The distance matrix was computed using the feature distances between two sulci.For position we computed a euclidean distance.For depth and length, we used an absolute distance.For orientation, we used , where . is the Frobenius norm.For shapes, the geodesic distance was used.The composite distance between any two sulci was also computed by using a combination of the scaled distances, where the w i could be chosen by trial and error.Fig. 2. Evaluating the number of dimensions required to represent the data.The scree plot (left), shows that 3 dimensions is adequate for the MDS map.The Shepard plot for 2D data (middle) shows more spread than the corresponding 3D plot (right).

Data and Initial Processing
We used a database of 18 T1-MR 3D SPGR brain images of healthy subjects, matched for sex (male), handedness (right), and age (35 ± 10 years).Five major sulci, the central, postcentral, superior frontal, sylvian fissure and superior temporal, were extracted from each hemisphere using the method described in [2,9].This gave us 18 × 10 = 180 sulci for classification.These primary sulci were chosen because they can consistently be identified in all normal individuals.The left and right precentral sulci were later included to test and extend the classifier.This expanded the dataset to 18 × 12 = 216 sulci.For analysis, the top and bottom sulcal lines extracted from each sulcus (see Figure 1(c)), were used.

Evaluating Dimensionality, Assessing Fit
MDS reproduces a high dimensional input in a lower dimensional space.A scree test was used to determine the dimension that best represents the sulcal data.In Figure 2, an elbow in the plot at d = 3, represents the optimal dimension beyond which only minor reductions in stress are acheived.The Shepard diagram, which plots the scaled MDS distances against the original dissimilarity data, was then used to assess fit for 2 and 3D data.A good fit is indicated by a minimal spread in the scatter plot and evidence of a monotonically ascending relation.A comparison between the two Shepard plots, shows that the output distances are highly correlated to the input dissimilarities for the 3D data.As a final test, the MDS map we obtain, reflects the actual arrangement of sulci in the brain.

Classification
Our goal was to label pre-segmented sulci making no assumptions about predefined relationships or any other information.For this we designed a 10-class MDS classifier.We evaluated the classification using leave-one-out (LOO) crossvalidation.The 18 sulci for each of the 10 sulcal classes were split 17:1 into a training set and a testing set.As expected, the best results for the 18 LOO iterations were obtained when spatial distances were used for the dissimilarities.We got an overall true positive rate of 90% and 90.6% for the top and bottom sulcal classifiers respectively (Table 1), and were able to correctly label all the superior frontal sulci (100%).The sylvian fissure, the most difficult class, had a success rate of 78%.The sylvian fissure is composed of a number of small folds that are difficult to extract completely.As a result, the features used to describe the sylvian fissure have a noisy multimodal distribution.
Classification was also done for shape, using geodesic paths for distances, and for scaled combinations of shape, position, length and depth.The shape classifier labeled all the central sulci, which are very regular, but was less successful with the other classes.Combining the features did not improve the success rate.Rather, the results achieved by the spatial distance classifier were compromised (Figure 4).The depth and length classifiers could not discriminate right and left brain sulci, and while the orientation information improved results in cases where the orientations of the sulci differed, these were offset by the errors introduced.An MDS spatial distance classifier offers several advantages.First, if n is the number of classes, we need only make n(n−1) 2 distance measurements, off-line, to construct a reference map.For classification, data from unlabeled sulci supplement this distance matrix.Second, our results suggest that we can very quickly assign a given sulcus to a small region of the left or right hemisphere.All 180 sulci were correctly identified when we relaxed the classification criteria to include the second nearest sulcus.Third, the classifier is robust to normal population variation.The results we obtained, were for data that had not been spatially normalized.Both the top and more stable bottom sulcal classifiers, moreover, gave comparable results (Table 1).We were even able to correctly identify a displaced and distended superior frontal sulcus, taken from a tumor patient.

Discussion
Central to our approach to sulcal classification is the use of a reference graph.This imposes several constraints which narrow down possible solutions so that a high average classification rate is achieved.The classifier works best when only a small number of well separated sulcal classes are used in the training set.We saw the performance degrade from 90.6% to 84.3% when two additional classes, the left and right precentral sulcus were added (Figure 5).Based on preliminary results and characterization, we propose a two stage classification scheme where the spatial distance classifier will serve as a front end to identify a region of interest.Specialized feature classifiers can then be applied in the second stage.A depth classifier, for instance, can discriminate between the sylvian fissure(SF) and the superior temporal sulcus(ST) (Figure 3).These two sulci were sometimes misclassified, one for the other, with the position classifier.

Fig. 1 .
Fig. 1.(a) The five major sulci used in this study.(b) Top: Sulcal curves, in blue, form a complex and widely varying pattern on the cortex.Bottom: Graph representation where the distances are the edges and the sulci are the nodes.(c) An individual sulcus.The top or superficial curve, which points in the direction of the outer cortex, is in blue.The bottom curve is in green.