Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE Trans Med Imaging. Author manuscript; available in PMC 2008 May 5.
Published in final edited form as:
PMCID: PMC2366175

Statistical Analyses of Brain Surfaces Using Gaussian Random Fields on 2-D Manifolds


Interest in the morphometric analysis of the brain and its subregions has recently intensified because growth or degeneration of the brain in health or illness affects not only the volume but also the shape of cortical and subcortical brain regions, and new image processing techniques permit detection of small and highly localized perturbations in shape or localized volume, with remarkable precision. An appropriate statistical representation of the shape of a brain region is essential, however, for detecting, localizing, and interpreting variability in its surface contour and for identifying differences in volume of the underlying tissue that produce that variability across individuals and groups of individuals. Our statistical representation of the shape of a brain region is defined by a reference region for that region and by a Gaussian random field (GRF) that is defined across the entire surface of the region. We first select a reference region from a set of segmented brain images of healthy individuals. The GRF is then estimated as the signed Euclidean distances between points on the surface of the reference region and the corresponding points on the corresponding region in images of brains that have been coregistered to the reference. Correspondences between points on these surfaces are defined through deformations of each region of a brain into the coordinate space of the reference region using the principles of fluid dynamics. The warped, coregistered region of each subject is then unwarped into its native space, simultaneously bringing into that space the map of corresponding points that was established when the surfaces of the subject and reference regions were tightly coregistered. The proposed statistical description of the shape of surface contours makes no assumptions, other than smoothness, about the shape of the region or its GRF. The description also allows for the detection and localization of statistically significant differences in the shapes of the surfaces across groups of subjects at both a fine and coarse scale. We demonstrate the effectiveness of these statistical methods by applying them to study differences in shape of the amygdala and hippocampus in a large sample of normal subjects and in subjects with attention deficit/hyperactivity disorder (ADHD).

Index Terms: Conformal mapping, Gaussian random field, manifold, Riemannian curvature tensor, statistical shape analyses, students’ t-statistics

I. Introduction

Maturational, neuroplastic, and degenerative processes in the central nervous system (CNS) are regionally specific, causing localized variability in brain morphology, within and across groups of healthy individuals and persons with various diseases, [1], [2] that can be measured in terms of differences in regional volumes, shapes, and asymmetries of brain regions [3]. Understanding normal and pathological morphology of the brain, therefore, requires detailed statistical analyses of the shapes of cortical and subcortical regions [4]-[13]. Statistical analyses of complex shapes depend critically, however, on the choice of method for formal geometric representation of the morphology of these surfaces.

Geometric representations provide a quantitative, summary descriptor for a morphological surface that can then be used for formal statistical analyses across surfaces in individuals or groups of individuals [14]-[17]. A number of geometric representations have been proposed to describe and analyze a shape. For example, after selecting points along the boundary of a region as a descriptor of the surface, principal components analysis (PCA) can quantify its global variability [18]. Alternatively, a linear sum of the Fourier basis functions can be used to describe the surface, thereby allowing a statistical analysis of an overall shape that is based on the weights of the Fourier basis functions [19], [20]. Statistical analyses using discrete, multiscale medial descriptors (m-reps) [21]-[23] of the surface can capture both coarse- and fine-scale variability in shape [14], [24], [25]. A representation of shape based on biological landmarks and their warps along thin-plate splines can be used to analyze the global characteristics of shape [26]. A parametric description of a surface [27] using spherical harmonics as basis functions can model the surface of a simply connected structure [28]. Finally, the statistical learning of shapes using distance maps has also been proposed [29]. All of these methods aim to reduce the complexity of a surface by retaining only a quantitative, summary descriptor that best describes the shape of that surface; hence, each local feature of the surface influences to some extent the value of this global descriptor. Whereas multiscale descriptors can be used to capture fine-scale variability in surface morphology [30], [24], [14], the statistical dependences of these various features are either poorly described or unknown, thus limiting the utility or validity of these descriptors for analyzing surface morphologies.

To overcome the limitations of multiscale descriptors, we use the entire set of points on a discretized surface as its descriptor, thereby allowing the detection of small local differences in surface morphology across individuals or groups of individuals. The statistical variability in morphology at a particular point on the surface can be studied by statistically comparing sets of corresponding points across subjects; using these descriptors for statistical comparisons, therefore, requires the accurate determination of points of correspondence for the surface across individuals. These correspondences across surfaces can be identified by first using affine transformations to coarsely register the surfaces within a common coordinate space [15]. Point correspondences across the surfaces can then be determined using a number of methods including geodesics [31], nearest neighbors [8], warpings based on fluid dynamics or elastic deformations [3], [32], measures of curvature [33], [34], robust point matching [35], [36], or hierarchical attribute matching [37].

These approaches, however, require the selection of a reference brain with points already identified on the surface of interest to serve as markers for labeling points of correspondence on the surfaces of analogous regions in brains that have been registered to this reference. After determining these points of correspondence, Euclidean distance is calculated between each point on the surface of the reference region and the corresponding point on the surface of each coregistered surface; for a group of subjects, this process generates for the group a set of Euclidean distances calculated for each point on the surface of the reference region. The reference region combined with the sets of distances at every point on its surface yields a complete statistical description of the surface’s morphology for the subjects who belong to that group. To compare surface morphologies across groups of subjects, this same statistical descriptor is calculated for subjects in the second group, which can then be used to compare the sets of distances at each point on the reference surface.

Students’ t-statistics [38] can be used for statistical comparisons of surface morphology across individuals and groups of subjects. An image depicting the t-statistics computed at each point on the surface of a reference region is called a t-map. Usually the null hypothesis to be tested is that no between-group differences will be detected in the signed Euclidean distances (distances for the points on the undeformed surface of each subject’s region that are positioned inside the boundary of the reference region are labeled as negative, whereas distances for points positioned outside of it are labeled as positive) at each point on the reference surface. However, because the number of points on the surface of a reference region will be large (ranging typically from several thousand to several hundred thousand), the false positive (or Type I) error rate quickly becomes unacceptably high [39]. Although Bonferroni correction [40], [41] can be employed to correct for the number of statistical comparisons performed, Bonferroni correction is not the ideal statistical adjustment [42] because the sets of Euclidean distances are not statistically independent of one another. Bonferroni correction therefore tends to be overly conservative, leading to an unacceptable false negative (or Type II) error rate. On the other hand, nonparametric statistical methods [43] do not account for the intrinsic curvature of the two-dimensional (2-D) surfaces embedded in R3 Euclidean space.

Our statistical model for morphological analyses describes a correlated, smooth random field defined on the surface of the reference region. Comparison of these random fields is desired for valid statistical analyses of morphological surfaces across groups of subjects. Under the null hypothesis of no shape difference across the entire surface, and for a given significance level (α) [38], if a correct threshold can be computed for the random field, then the t-map can be thresholded to yield statistically significant differences in surface morphology. The smallest significance level at which a voxel is statistically significant is called the probability value, or p-value, for that point. A p-value for a zero-mean and homogeneous Gaussian random field (GRF), f(t), that is defined on N-dimensional Euclidean space, RN, has been developed previously [44]; this work has been used extensively to detect statistically significant voxels in Functional Magnetic Resonance Imaging (FMRI) datasets [45], [46] and has been extended to non-Gaussian random fields [46], [47]. These results have also been used to detect statistically significant differences in the shapes of brain regions between groups. This prior work, however, does not account for the intrinsic curvature of surfaces. Additionally, nonparametric statistical methods [43] have been used for the analyses of surface morphometry and although they make minimal assumptions about the distribution of the samples, how nonparametric methods account for correlations across the surface in Euclidean distances is unclear. Additionally, current nonparametric methods also ignore intrinsic curvature of the surfaces, and therefore, would not correctly address the issue of multiple comparisons. The degree of improvement gained through analyses that account for the curvature of surfaces, as compared to the analyses that ignore it, will depend upon the surface morphology of regions under investigation. For example, because amygdala is a small region, and hence has considerable curvature across its surface, analyses that account for curvature will be more accurate because the random fields are correctly modeled. Because our statistical model for the comparison of morphological surfaces defines a random field on the surface of the reference region, i.e., on a 2-D manifold in three-dimensional (3-D) space, we use previously described, detailed methods [48] to derive expressions for the statistical analysis of differences in surface morphology across groups of individuals.

II. Methods

Although the concept of using sets of signed Euclidean distances at each point on a surface as a statistical description of that surface is relatively simple, the actual statistical analysis based on this model is challenging because the model comprise of a general random field across the surface. To simplify the analysis, we assume that the sets of signed Euclidean distances at each point of a reference surface form a Gaussian random field (GRF) f on the surface of the reference region; furthermore, we assume that the GRF is smooth and suitably regular [44] and that the surface is a 2-D manifold embedded in three-dimensional (3-D) Euclidean space, R3.

Statistical analyses using the expected Euler characteristics of a GRF defined in Rn Euclidean space have been described previously [44]. For our analyses, an excursion set A of f is defined as A(f, u) = f−1[u, +∞), i.e., A is a set of all points in a space in which the values of a realization of the random field exceeds a specified value u. The Euler characteristics, χ(A(f, u)), of the excursion set is then defined by extending the concept of the number of up-crossings of a one-parameter random process. Others [44] have shown that the expected value E[χ(A(f, u))] for a smooth, isotropic GRF, f on Rn, and an open domain T [subset or is implied by] Rn bounded by a C2 compact hyper-surface, is asymptotically equal to the probability that the value of the GRF f will be greater than the value u, i.e.,


where f|T is the restriction of f to T. Thus, the expected Euler characteristics are shown to be good approximations of Gaussian excursion probabilities [44].

A. Evaluation of the Expected Euler Characteristics on A 2-D Manifold

In our statistical model, a random field is defined on the surface of the reference region, which is a 2-D manifold in 3-D space. Thus, to compute P[supy[set membership]T f(y) > u] we need to calculate the expected Euler characteristics for the random fields defined on this 2-D manifold. The expected Euler characteristics, E[χ(A(f, u))]—of a suitably regular and centered Gaussian random field of unit variance on a C3 n-dimensional compact manifold M—have been derived previously [48]. Let g denote the Riemannian metric [49] induced by f, and let R denote its Riemannian curvature tensor. Then the expected Euler characteristics can be written as




are the Minkowski functionals, Hj is the jth Hermite polynomial, Φ is the cumulative distribution of a standard Gaussian random variable, gij = g(Ei, Ej), the Riemannian metric defined on the manifold, and Volg is defined as Volg=det(gij)(i=1ndxi). x21330(M) integrates the trace of the Riemannian curvature tensor over the entire surface, and x21332(M) computes the total surface area of the manifold. Thus, from (1) and (2), we calculate


where Volg is the total surface area of the manifold M, and TrM (R) is the trace of the Riemannian curvature tensor at a point on the surface.

The Trace of the curvature tensor TrM (R) for an orthonormal measurable section (Xi)1≤i≤2 of O(M) with dual frames (θi)1≤i≤2 is derived as


where S(2) is the symmetric group of permutations of 2 letters and εσ is the sign of a permutation [48]. Using the symmetric and skew symmetric properties of the Riemannian curvature tensor [49], the trace TrM (R) can be further simplified to be TrM(R)=R1212X. Thus, the expected value of the Euler characteristics in (3) can now be written as


for any X in the orthonormal frame section O(M) on the manifold. Thus, to evaluate the expected Euler characteristics, E[χ(A(f, u))], only R1212X needs to be evaluated.

To evaluate the Riemannian curvature tensor RX, consider (Ei)1≤i≤2, the coordinate vector fields on the surface manifold M, and gij = g(Ei, Ej), the Riemannian metric defined on the manifold. Let (dei)1≤i≤2 denote the dual basis of (Ei)1≤i≤2, i.e., g = gij dei [multiply sign in circle] dej. Then, the dual frames (dei)1≤i≤2 and (θi)1≤i≤2 are related as


where gij1/2 are defined as gij1/2=g(Ei,Xj). In matrix notation, the above equation can be written as


where we defined


If components of the Riemannian curvature tensor are RijklE in the coordinate frames, then RijklX are evaluated as [48], [49]


which can be further simplified using the symmetric and antisymmetric properties of the Riemannian curvature tensor [49] to be


From (5), it can be shown easily that (g11(1/2)g22(1/2)g12(1/2)g21(1/2))=1/det(gij1/2); hence, R1212X can be written as


For a given coordinate frame field, (Ei)1≤i≤2, an orthonormal frame field, (Xi)1≤i≤2, can be chosen as


Thus, gij1/2, s at a point p on the surface M, can be evaluated as


Thus, to evaluate the expected Euler characteristics in (4), we need to evaluate gp(Ei, Ej) and RijklE. Following others [48], gp(Ei, Ej) is defined as


and the components of the Riemannian curvature tensor RijklE are given as [49]


where, Γikl are the Christoffel’s symbol of the second kind [49], and (gst) = (gst)−1 is the inverse of the Riemannian metric, such that k=1ngikgkj=δij. Explicitly, (gst) are evaluated as


Thus, R1212E in (7) can be written as


where the Christoffel’s symbols Γijk are evaluated as [49]


Thus, we need to evaluate [partial differential]2f/[partial differential]u2, [partial differential]2f/[partial differential]v2, [partial differential]2f/[partial differential]u[partial differential]v, [partial differential]f/[partial differential]u, and [partial differential]f/[partial differential]v, and the expected values of their products. The first two terms, E1Γ212 and E2Γ112, in (9) can be written as


Therefore, their difference


involves only second-order partial derivatives of the random field f. Finally, (9), (8), and (7) are substituted in (4) to evaluate the expected Euler characteristics of the GRF, f, defined on the 2-D manifold.

B. Estimating Coordinate Vectors and Their Derivatives

Evaluation of the Christoffel’s symbols in (9) and the Riemannian metric in (8) requires determination of the coordinate vectors in 3-D Euclidean space and their derivatives as well as the first- and second-order derivatives of a realization of the random field with respect to the coordinate vectors. We compute derivatives of a realization of the random field by first locally fitting a quadratic surface at a point. Fitted quadratic surface is locally smoothed approximation of the realization. Then derivatives of the quadratic surface yield smoothed approximations of the derivatives of the realization of the random field. Finally, the expected values of the derivatives are estimated by calculating the numerical averages of appropriate terms across all realizations of the random field.

Let κ be a diffeomorphism from a patch in 2-D to a patch on the surface in 3-D, i.e., (x) = κ(u, v), where (u, v) and (x) are the coordinates of the points in the 2-D patch and on the surface, respectively. The coordinate vectors E1, E2, and E12 are defined as E1 = xu, E2 = xv, and E12 = xuv, That is


The coordinate vectors can be evaluated directly if the diffeomorphic mapping κ can be written in an analytic form. However, because the analytic form generally is unknown, the mapping is approximated locally by a quadratic equation at each point of the surface. Thus, once the mapping κ has been determined, the (x, y, z) coordinates of a point x on the surface are written in the quadratic form as a function of (u, v) as


For six neighboring points on a surface, unknowns ax, bx, cx, dx, ex, fx can be evaluated by solving the linear equations for six neighboring points. Then the coordinate vectors and their derivatives are computed as


Note that to evaluate the partial derivatives Ei f of a realization of the GRF f, we can write


Similarly, the partial derivatives of a realization of the random field f are evaluated by approximating the f(u, v) as a quadratic function of (u, v), i.e.,


Using this approximation and solving for the unknowns af, bf, cf, df, ef, and ff using six neighboring points, the partial derivatives of f are evaluated as


For a point on the surface with fewer than six neighbors, a linear fit of the realization of the random field is used to evaluate the partial derivatives as


Finally, the expected values of the partial derivatives and their products are approximated by calculating the numerical averages across each realization of the random field.

C. Least Squares Conformal Maps

Evaluation of the coordinate vectors and the partial derivatives of the GRF in (10) and (11) requires a diffeomorphism κ between a patch on the surface in the 3-D Euclidean space and a patch in the 2-D Euclidean space. Least squares conformal maps [50] are an appropriate mapping and are estimated using previously described algorithms [51]. While the conformal mapping we used required cutting the surface to define the κ for each patch, cutting of surfaces itself is not necessary for analysis of the surface morphology. For example, a spherical parameterization of the surfaces will avoid cutting the surface arbitrarily.

Let, κ denote a mapping from a subset U of R2 into a subset W of a surface in R3, (u, v) denote the coordinates in the subset U, and let (x, y, z) denote the coordinates in R3; then, ([partial differential]κ/[partial differential]u) (u, v) and ([partial differential]κ/[partial differential]v) (u, v) denote the coordinate vectors at the point κ(u, v) in W. The homeomorphic mapping κ, along with the subset W, that maps W to a disc in R2 is called a chart. If N(u, v) denotes the unit normal vector at the point κ(u, v), then the mapping κ is conformal if the relation


is satisfied—i.e., if the vectors ([partial differential]κ/[partial differential]u) (u, v) and ([partial differential]κ/[partial differential]v) (u, v) are orthogonal and have the same norm. Furthermore, if κ/u2=1, then the conformal mapping is locally isotropic, i.e., it preserves angles and distances.

A parameterization of a surface satisfying the conformality relation [(12)] can always be found because any surface, S, is homeomorphic to a disc (according to Riemann’s theorem [49]). The conformality condition may not be completely satisfied, however, if additional constraints are imposed on the mapping—for example, if the edges of the triangulations are mapped to straight lines, or if the mappings vary linearly in each triangle. Thus, a mapping is estimated that satisfies the above constraints and is close to being conformal in the least squares sense.

Consider a triangulation (G) of a surface (S) given by G = {[1,…,n]n≥3, T, (pj)1≤jn}, where [1,…,n]n≥3 denotes the set of vertices, (pj)1≤jn is the location of a vertex j [set membership] R3, and T denotes the set of triangles represented by triples of vertices in the triangulation. In a local orthonormal basis, let the coordinates of a vertex j be denoted by (x1, y1), (x2, y2), (x3, y3). Then the mapping κ maps (u, v) to (x, y); i.e., κ(u, v) = (x(u, v) y(u, v))T. Using complex number notations, we can denote the mapping as κ(u, v) = x(u, v) + iy(u, v). The conformality equation can then be written as


Let u denote the inverse mapping on a triangle given by u(x, y) = (u(x, y) v(x, y))T; then, by the theorem of derivatives of the inverse mapping, the conformality equation can be described as


Because the above equation cannot be satisfied with the conformality constraints, the cost C(T) for a triangle T with area AT is defined as


Hence, the total cost C(T), for the triangulation G, minimized to determine the least-squares conformality mapping, is given as


The optimal mapping is then estimated by inverting the appropriate matrix [50].

D. Example of Application in a Real-World Dataset

As an example of analyzing differences in shape in real world dataset based on GRF-theory, we compared the shapes of the amygdala and hippocampus across a group of 62 healthy individuals and a group of 51 subjects with attention deficit/hyperactivity disorder (ADHD). The amygdala and hippocampus were manually segmented in a randomized order by an expert in neuroanatomy using a well-developed and highly reliable algorithm [52] while blind to subject characteristics.

1) Image Acquisition

Head positioning was standardized using cantho-meatal landmarks. High-resolution T1-weighted MRIs were acquired from a single 1.5-T scanner (GE Signa; General Electric, Milwaukee, Wisconsin) using a sagittal 3-D volume spoiled gradient echo sequence with repetition time = 24 msec, echo time = 5 msec, 45° flip angle, frequency encoding superior/inferior, no wrap, 256 × 192 matrix, field of view = 30 cm, 2 excitations, slice thickness = 1.2 cm, and 124 contiguous slices encoded for sagittal slice reconstruction, voxel dimensions 1.17 × 1.17 × 1.2 mm.

2) Subjects

We acquired anatomical MR images in 51 subjects who met DSM-IV criteria [53] for the combined type of ADHD and in 62 healthy individuals. Subjects were 7 to 18 years of age. Healthy individuals had no lifetime history of tic disorder, obsessive compulsive disorder (OCD), ADHD, or current DSM-IV Axis I disorder. Subjects were predominantly right-handed.

3) Selection of the Reference Region

Detection, localization, and interpretation of the statistically significant differences across groups of subjects conceivably could depend on the choice of subject in whom the reference region is defined because the estimated registration parameters and the established correspondences of points on the surfaces of the region will depend upon the degree to which the selected reference brain is representative of the population of subjects being studied. Use of a synthetic, average reference region can also be considered. Creating a synthetic average, however, may not always be possible, particularly for complex regions. For example, the pattern of gyri and sulci on the cortical surface vary significantly across individuals even within the healthy population [54]. Thus, when not all subjects have a particular gyrus or sulcus, how to generate a synthetic, average cortical surface is not at all clear. Moreover, interpreting any statistically significant differences across surfaces, regardless of how the synthetic average is generated, will be difficult.

To overcome the limitations associated with use of a synthetic reference, we instead selected as a reference a brain that was as representative as possible of the sample of subjects sample being studied. First, one subject (a 23.5-year-old right-handed, Caucasian male) whose demographic characteristics were nearest the group average was selected preliminarily as the reference brain. The brains for all remaining subjects in the sample were registered to that preliminary reference. The point correspondences across their surfaces were determined (described below), and then the distances of those points from the corresponding points on the reference surface were calculated. The brain for which all points across its surface were closest, in the least squares sense to the average of the distances across those points for the entire sample was selected to be the final reference. The registration process, determination of point correspondences, and calculation of distances across surfaces was then repeated for all subjects in the sample against this final reference brain. Distances between surfaces were then compared across individuals or groups.

4) Smoothing the Random Field

To use the theory of GRFs, each realization of a random field should be smooth [55]. Effectively, we use two methods to smooth each realization of the random field. First, an isotropic Gaussian kernel that is five pixels wide and has a full-width at half-maximum (FWHM) of three pixels is used to smooth each realization. Second, the coordinate vectors and their derivatives are estimated by locally fitting a quadratic form to each realization. This fitting of the quadratic form can be regarded as a further smoothing of the random field.

5) Controlling for Covariates

Surface morphology, and hence computed signed Euclidean distances, of brain regions could depend upon Sex and Age of subjects, in addition to their clinical diagnosis. Multiple regression analysis allows assessment of relationship between an independent variable and the dependent variable after taking into account the remaining independent variables [38]. Therefore, to control for the effects of covariates (age, sex) on surface morphology, we performed a multiple linear regression analysis [38] at each point on the reference surface


where di was the set of signed Euclidean distances for the ith subject. We computed the regression coefficient β3 between the distances and clinical diagnosis (having ADHD or not) and then computed the p-value of this coefficient using the Student’s t-test. We did not covary for Whole Brain Volumes (WBVs) because each subject’s cerebrum was coregistered to the reference subject’s cerebrum, thereby accounting for scaling effects within individual subjects.

For the statistical analysis based on the theory of GRFs, we first computed the residuals of the distances using multiple linear regression analysis without including the ADHD term in the model, i.e.,


Then a statistical analysis on the residual of the distances localized significant differences between the diagnostic groups.

6) Surface Analysis of the Amygdala and Hippocampus

The cerebrums of all subjects were first coregistered to the cerebrum of the subject whose brain was selected as the reference (described above in Sections II-D–III). Registration was performed using a similarity transformation, the parameters of which (three translation, three rotations, and a global scale) were estimated such that the mutual information [56] across the entire image for each subject was maximized. While different methods to register images will scale brains differently, thus affecting the group comparisons of brain subregions, our prior work [57] has shown that our method accurately and consistently estimates scaling parameters over a wide range of values. The amygdala and hippocampus in the coregistered brains were then additionally rigidly registered to the reference to further improve their coregistration across subjects. This two-step strategy for registration preserved differences in the size of the regions that were not addressed by coregistration and rescaling of the cerebrums.

Correspondences between the points on the surfaces of the regions were then identified by deforming these regions within the brains of all subjects into the reference region using an algorithm that is based on fluid dynamics [58], [57]. This method has been shown to perform well for real and synthetic data [57]. Signed Euclidean distances between the corresponding points of the regions of the subject and reference brains were then calculated. Thus, a group of 62 healthy individuals determined a set of 62 signed Euclidean distances at each point on the surface of the reference region. These sets of signed Euclidean distances across the entire surface of the reference determined the random field on its surface.

The morphological descriptions of the surfaces of the regions were compared across diagnostic groups by comparing the sets of distances at each point on the reference surface. Student’s t-statistics were calculated for each point on the surface of the reference region and were converted to Gaussian statistics. The corresponding z-map was then constructed.

7) Converting the T-Statistic to a Gaussian Statistic

Equation (4) was used to evaluate the expected Euler characteristics and hence the significance level for a specified threshold, u. However, in (4) the expected Euler characteristics were derived only for a GRF defined on a 2-D manifold. Thus, the Student’s t-statistics t-statistics t(x) at a location (x) was converted to the Gaussian-distributed random variable z(x) using the equation, z(x) = p−1{pn[t(x)]} [45]; here pn(t) was the probability that a Student’s-t distributed random variable of (n−1) degrees of freedom exceeded t, and p(z) was the probability that a standard Gaussian random variable exceeded a value of z. The computed statistics g(x) were then thresholded using the significance level estimated with (4) to determine the locations where the morphologies of the surfaces of the sample N of subjects deviated from the mean of the control population to a statistically significant degree.

For a specified significance level (α = 0.05), statistically significant differences were then estimated by thresholding the z-map at the threshold uα = 4.61. This threshold was computed using (4) for the normalized GRF, which that was estimated from the signed Euclidean distances for the entire healthy individuals.

8) Expected Euler Characteristic of the T-Field

We can write the expected Euler characteristic of a random t-field with ν degrees-of-freedom defined on 2-D manifold as






are the EC densities for the t-field [46], [59]. Using (15) obviates the need to convert the Student t-statistic to the Gaussian statistic. However, for many degrees-of-freedom ν the two analyses, generated using (1) and (15), will be similar.

III. Results

A. Estimated Random Field and Euler Characteristics

Evaluation of the expected Euler characteristics in (4) assumes that the estimated random field is Gaussian distributed. Skewness (a measure of asymmetry) and kurtosis (a measure of peakedness relative to a Gaussian distribution) [60] were evaluated for two sets of signed Euclidean distances at two different points on the surface of the left amygdala chosen as the reference region (Fig. 2). The sets of distances were obtained by matching the left amygdala of the 62 healthy individuals to the reference region. Histograms in Fig. 2(a) and (b) show that the distributions of the distances are close to Gaussian, an assertion supported by the quantile-quantile (QQ) plots of these distances in Fig. 2(c) and (d).

Fig. 2
Distribution of signed euclidean distances for a group of 62 Healthy Individuals. Histogram (a) and QQ plot (b) of the estimated signed Euclidean distances for a group of 62 healthy individuals. Plots of random variables that are zero-mean, Gaussian-distributed ...

The expected Euler characteristics, E[χ(A(f, u))] in (4), is valid asymptotically only for large values of the thresholds u; for small values of u, E[χ(A(f, u))] can be greater than 1, and they can even be negative (see Fig. 3). Thus, (4) is used to evaluate thresholds for small significance levels only.

Fig. 3
Expected Euler characteristics, E[χ(A(f, u))], as a function of the threshold U. The expected Euler characteristic was computed using the normalized GRF, over the left amygdala, for the group of 62 healthy individuals. Here, E[χ(A(f, u ...

B. Analyses of the Amygdala and Hippocampus

The flattened amygdala and hippocampus are shown in rows A and D of Fig. 1. The color patterns on the flattened surfaces demonstrate that angles of the surface triangulation have been preserved locally by the conformal mapping.

Fig. 1
Surface analyses of the amygdala and hippocampus. Rows A-C: left amygdala and left hippocampus. Rows D-F: right amygdala and right hippocampus. Columns 1–4: group comparison using only the Student’s t-test computed on the β3 regression ...

Surface analyses of the group of 62 healthy individuals demonstrate (Fig. 1) that the analyses using either the uncorrected Student’s t-test alone or in combination with analyses based on GRF theory captured similar large-scale differences between the surface of the reference regions and the corresponding surfaces of the healthy individuals. Use of GRF theory, however, produced far fewer statistically significant differences, as expected. The significant differences between the surfaces in the reference subject and those in the controls therefore represent significant individual variability of even the reference subject relative to the group average in this large sample of healthy individuals.

Comparing the group of 62 healthy individuals to the group of 51 subjects with ADHD, the Student’s t-test alone detected local decreases in volume of the amygdala and local increases in volume of the hippocampus in the subjects with ADHD. GRF analyses detected similar large regions of differences between the groups of subjects; however, because GRFs address the multiple comparisons using Student’s t-tests, these analyses were conservative and detected fewer small regions of differences as being significant.

We studied the effects of differing methods of segmentation of the 3-D surface on the results of our statistical analyses (Fig. 4). This produced findings highly similar to the original (Fig. 1, row C). From these limited analyses, we conclude that the particular methods used to segment and parameterize the surface seem not to affect the results of surface analyses using our method for GRF correction.

Fig. 4
Effect of cut on the surface analysis of the left amygdala. We obtained a differing segmentation of the 3-D surface for the left amygdala using a set of parameters that segmented the surface into many pieces. Images from left to right: The segmentation ...

To study the effect of segmenting the 3-D surface on the statistical analysis, we parameterized the surface of the left amygdala with a segmentation than differs from those in Fig. 1, row A. Even with different segmentation and parameterization of the surface of the left amygdala, the two analysis in Figs. 4 and and1,1, row C matched well. Therefore, our method for segmenting and parameterizing the surface did not affect the surface analyses.

Analyses that used the expected EC where the curvature term [the first term on the right in (4)] was set equal to zero were more conservative for this dataset than the analyses that accounted for the intrinsic curvature of surfaces (Fig. 5). Ignoring the intrinsic curvature of surfaces may therefore increase the rate of false negatives in the analyses and thereby increase the difficulty of detecting features of interest in surface morphology.

Fig. 5
Surface analysis without the curvature term to compute the expected euler characteristic (EC). In these analyses of the surface of the amygdala, the signed Euclidean distances were modeled as Gaussian random fields (GRFs). However, to study the effect ...

Analyses using expected ECs computed for the t-field (Fig. 6) matched well with the analyses computed for the GRF in which we first converted the t-statistic to the Gaussian statistic (Fig. 6, rows C and F, columns 5–8). Because analyses based on a t-field that has many degrees-of-freedom (as is the true for our dataset) converge to the analyses based on GRFs, their results here are similar.

Fig. 6
Analyses with expected euler characteristic (EC) Computed for T-Field. Analyses of the surfaces of the amygdala and hippocampus for a group of 62 healthy individuals and 51 subjects with ADHD. For these analyses, we modeled the random field of signed ...

IV. Discussion

To compare surface morphologies across individuals and groups, we devised a statistical assessment in which a random field of distances is defined between points on the surface of a reference region and the corresponding points on the surfaces of the regions in the brains of other individuals in the study sample. The random field was estimated by first using a similarity transformation to coregister the surface of the region in each subject to the surface of the reference region. We then identified correspondences between the points on the surface of subject regions and the points on the surface of the reference region. We calculated the sets of signed Euclidean distances between each point on the surface of the reference region and the corresponding point on the subject surfaces. These sets of distances defined a random field across the entire surface of the reference region. To compare groups, we calculated sets of distances from the same reference to each group. The corresponding random fields were then estimated and compared at each point.

In processing medical images, a surface historically has been represented by features such as landmark points [26], basis functions [19], [27], m-reps [21], or a discrete subset of points along the boundary of an object [18]. The complexity of the representation was thereby reduced, and the representation retained only those few morphological features that were deemed to capture the most important variance in surface topography. A statistical description of a shape was then estimated based on the probability distribution of each feature using a sample of example shapes, which allowed a comparison of the shapes between groups. However, because each feature tends to affect the shape globally, studying fine-scale shape differences in shape using these shape representations has proved difficult. Moreover, fine-scale m-reps have been proposed to detect local shape differences [14], yet localization of these shape differences has also proved challenging.

A statistical description of surface morphologies that is based on a random field defined over the surface of a reference region allows a comparison of surfaces at both fine and coarse scales of description, while simultaneously accounting for the nonindependence and intercorrelations of topological features across the surface. Other investigators have described methods for calculating the probability of a GRF greater than a specified value; this probability has been shown to be approximately equal to the expected value of the Euler characteristics [48]. We used this calculation to derive the mathematical expressions for the statistical comparison of surfaces that are 2-D manifolds embedded in a 3-D Euclidean space.

Application of our method to the comparison of surfaces of the amygdala and hippocampus in a large, real-world dataset confirmed that analyses based on the theory of GRF produce results that are conservative compared with those based on the Student’s t-test for which the signed Euclidean distances at each point of the surface are assumed to be independent statistically from one another. Both maps detected similar large-scale differences across diagnostic groups. The maps generated the Student’s t-test, however, showed a number of smaller-scale differences across surfaces compared with the maps generated using GRFs.

To use the theory of GRFs, the random fields need to be optimally smoothed. We used an isotropic Gaussian kernel of three pixels FWHM to smooth the GRF defined on the 2-D manifold. Further smoothing of the GRF was implicit in local fitting of quadratic forms to estimate derivates of the GRF. Because we estimated the expected EC from the smoothed GRF, results of the anaylses could be affected by the specific methods of smoothing employed. While our smoothing parameters determined empirically performed well, we will systematically investigate the effects of smoothing parameters, such as isotropic or anisotropic Gaussian kernel, FWHM of the kernel, different shapes of kernel, on the analysis to determine optimal smoothing parameters.

When comparing surface morphologies of brain regions across individuals or groups, we must correct for scaling effects in the brain that typically contribute enormously to interindividual variability in size of brain subregions [61]. We correct for differences in scaling across subjects by scaling each individual brain to the reference brain and then scaling the manually defined subregion in each brain by the same amount. However, if the average whole brain volume (WBV) differs across groups of subjects but the absolute size of the subregion does not, global scaling may obfuscate the analysis of size of the subregion. Rescaling the smaller WBVs in a patient group that has subregions of normal absolute size, for example, would cause those subregions to appear enlarged relative to those of a control group. In this instance, two separate analyses—one with and another without global scaling—would help to correctly identify the nature of differences in surface morphology across groups. In our example, the difference between the average global scaling parameters for the two groups (1.05751 for the ADHD group and 1.05547 for the healthy individuals) was not statistically significant (p = 0.758). Therefore, our analyses that included rescaling of volumes of the subregions was likely valid, and reanalysis without scaling was unnecessary.

Cutting the surface arbitrarily, as we do in this method for GRF-based correction for multiple statistical comparisons on a 2-D manifold in 3-D space, may affect the results of the statistical analyses in some instances, especially at the boundary of a cut. Therefore, alternative methods for parameterizing the surface that avoid cutting—for example, spherical parameterization [19]– may sometimes be desirable. Practical implementations of both spherical harmonics and superquadratics, however, inherently lead to smoothing of the surface, thereby potentially smoothing out fine-scale differences. Our method of cutting the surface, however, preserves surface contours by not smoothing them in an arbitrary fashion, and it therefore is able to detect small-scale differences in surface morphologies, which generally is a decided advantage, give that detection of statistical effects across a small spatial scale is the primary goal of these statistical analyses. In addition, our method for parameterizing surfaces satisfies the constraint that the angles between points of the surface are preserved, thereby allowing for an undistorted visualization of the 3-D surface in 2-D space; this constraint is not explicitly imposed in current implementations of spherical harmonic representations of surfaces. Our method can therefore be used to unfold surfaces with complex geometry—for example, the surface of the cortical mantel—and it will assist in the statistical analysis and interpretation of morphology deep inside sulci by providing to an investigator an undistorted representation of that complex surface. Furthermore, particular instantiations of our arbitrary cutting of the surface contour will have minimal influence on our localized estimations of the curvature of the surface and of derivatives of the GRF because we smooth the GRF data both on the 3-D surface and on the flattened surface, and because we estimate the local curvature of the surface by fitting second-order polynomial to the surface.

Arbitrarily cutting the surface may affect analyses in some instances, especially at the boundary of a cut. Therefore, parameterizations of the surfaces, for example spherical parameterization [19], that avoid cutting may sometimes be desirable. However, by not smoothing a surface, our method preserves surface contours. This is important because our goal is to detect small-scale differences in the surface morphologies, and by preserving surface contours, our method allows for the detection of these fine-scale differences. On the other hand, practical implementations of both spherical harmonics and superquadratics inherently lead to smoothing of the surface, thereby potentially smoothing out fine-scale differences. In addition, our method for parameterizing the surface satisfies the constraint that the angles be preserved, thereby allowing for an undistorted visualization of the 3-D surface in 2-D space; this constraint is not explicitly imposed in the current implementations of spherical harmonic representations of surfaces. Our method can therefore be used to unfold surfaces with complex geometry—for example, the surface of the cortical mantel—and will assist in the interpretation of the differences detected deep inside sulci by presenting an undistorted view to an investigator. Furthermore, because we smooth the GRF data both on the 3-D surface and on the flattened surface, and we estimate the local curvature of the surface by fitting second-order polynomial to the surface, localized operations to estimate the curvature of the surface and derivatives of the GRF will not be affected by cutting. Additionally, cutting the surfaces affects equally the GRF for the two groups of subjects. Thus, in this instance, arbitrarily cutting surfaces is not likely to lead to erroneous detection of between-group differences.


The authors would like to thank J. Royal for his editorial assistance.

This work was supported in part by the National Institute of Mental Health (NIMH) under Grants MH068318, MH36197, and K02-74677, in part by the National Institute on Drug Abuse (NIDA) under Grant DA017820, in part by The Mental Health Research Association NARSAD Independent Investigator Award, in part by the Thomas D. Klingenstein and Nancy D. Perlman Family Fund, and in part by the Suzanne Crosby Murphy Endowment at Columbia University.


1. Van Essen D. A tension–based theory of morphogenesis and compact wiring in the central nernous system. Nature. 1997;385(6614):313–318. [PubMed]
2. Van Essen D, Drury H. Structural and functional analyses of human cerebral cortex using surface–based atlas. J Neurosci. 1997;17:7079–7102. [PubMed]
3. Thompson P, Giedd J, Woods R, MacDonald D, Evans A, Toga A. Growth patterns in the developing human brain detected using continuum–mechanical tensor mapping. Nature. 2000;404:190–193. [PubMed]
4. Sowell E, Peterson B, Thompson P, Welcome S, Henkenius A, Toga A. Mapping cortical change across the human life span. Nature Neurosci. 2003:309–315. [PubMed]
5. Shenton ME, Gerig G, McCarley RW, Szekely G, Kikinis R. Amygdala–hippocampal shape differences in schizophrenia: The application of 3-D shape models to volumetric MR data. Psychiatric Res Neuroimag. 2002;115:15–35. [PMC free article] [PubMed]
6. Wang L, Joshi S, Miller M, Csernansky J. Statistical analysis of hippocampal asymmetry in schizophrenia. NeuroImage. 2001;14:531–545. [PubMed]
7. Thompson P, Mega M, Vidal C, Rapoport J, Toga A. Detecting disease-specific patterns of brain structure using cortical pattern matching and a population-based probabilistic brain atlas. Proc 17th Int Conf Inf Process Med Imag (IPMI) 2001:488–501. [PMC free article] [PubMed]
8. Zilles K, Kawashima R, Dabringhaus A, Fukuda H, Thorsten S. Hemispheric shape of European and Japanese brains: 3-D MRI analysis of intersubject variability, ethical, and gender differences. NeuroImage. 2001;13:262–271. [PubMed]
9. Thompson P, Moussai J, Zohoori S, Goldkorn A, Khan A, Mega GCJ, Small MSa, Toga A. Cortical variability and asymmetry in normal aging and Alzheimer’s disease. Cereb Cortex. 1998;8:492–509. [PubMed]
10. Salat D, Buckner R, Snyder A, Greve D, Desikan R, Busa E, Morris J, Dale A, Fischl B. Thinning of the cerebral cortex in aging. Cereb Cortex. 2004 [PubMed]
11. Shen L, Ford J, Makedon F, Flashma L, Saykin A. Surface-based morphometric analysis for hippocampal shape in schizophrenia. Human Brain Mapp. 2003
12. Martin J, Pentland A, Sclaroff S, Kikinis R. Characterization of neuropathological shape deformations. IEEE Trans Pattern Anal Mach Intell. 1998 Feb;20(2):97–112.
13. Csernansky JG, et al. Hippocampal morphometry in schizophrenia by high dimensional brain mapping. Proc Nat Acad Sci. 1998;95(19):11406–11411. [PubMed]
14. Yushkevich P, et al. Proc IPMI’01. Vol. 2082. LNCS; 2001. Intuitive, localized analysis of shape variability; pp. 402–408.
15. Dryden I, Mardia K. Statistical Shape Analysis. New York: Wiley; 1998.
16. Davatzikos C, et al. A computerized method for morphological analysis of the corpus callosum. J Comp Assis Tomogr. 1996;20:88–97. [PubMed]
17. Golland P, Fischl B, Spiridon M, Kanwisher N, Buckner RL, Shenton ME, Kikinis R, Dale A, Grimson WEL. Discriminative analysis for image-based studies. MICCAI’02. 2002;2488:508–515.
18. Cootes T, et al. Active shape models–their training and application. Comput Vision Image Understand. 1995;61(1):38–59.
19. Staib L, Duncan J. Model–based deformable surface finding for medical images. IEEE Trans Med Imag. 1996 May;15(5):720–731. [PubMed]
20. Szekely G, Kelemen A, Brechbuhler C, Gerig G. Segmentation of 2-D and 3-D objects from MRI volume data using constrained elastic deformations of flexible Fourier contour and surface models. Med Image Anal. 1996;(1):19–34. [PubMed]
21. Pizer S, Fletcher T, Fridman Y, Fritsch D, Gash A, Glotzer J, Joshi S, Thall A, Tracton G, Yushkevich P, Chaney E. Deformable M–reps for 3-D medical image segmentation. Int J Comput Vision. 2003;55(2):85–106. [PMC free article] [PubMed]
22. Styner M, Gerig G, Lieberman J, Jones D, Weinberger D. Statistical shape analysis of neuroanatomical structures based on medial models. Med Image Anal. 2003;7:207–220. [PubMed]
23. Styner M, Lieberman JA, Pantazis D, Gerig G. Boundary and medial shape analysis of the hippocampus in schizophrenia. Med Image Anal. 2004;8(3) [PubMed]
24. Joshi S, Pizer SM, Fletcher PT, Yushkevich P, Thall A, Marron JS. Multi–scale deformable model segmentation and statistical shape analysis using medial descriptions. IEEE Trans Med Imag. 2002 May;21(5):538–550. [PubMed]
25. Fletcher P, Joshi S, Lu C, Pizer SM. Gaussian distribution on lie groups and their applications to statistical shape analysis. Inf Process Med Imag (IPMI’03) 2003:450–462. [PubMed]
26. Bookstein F. Morphometric Tools for Landmark Data. Cambridge, U.K.: Cambridge Univ.; 1991.
27. Brechbuhler C, Gerig G, Kubler O. Parameterization of closed surfaces for 3-D shape description. Comput Vision Image Understand. 1995;61:154–170.
28. Shen L, Ford J, Makedon F, Wang Y, Steinberg T, Ye S, Saykin A. Med Image Comput Comput-Assisted Intervention (MICCAI’03) Vol. 2879. LNCS; 2003. Morphometric analysis of brain structures for improved discrimination; pp. 513–520.
29. Gollard P, Grimson WEL, Shenton ME, Kikinis R. MICCAI’00. Vol. 1935. LNCS; 2000. Small sample size learining for shape analysis of anatomical structures; pp. 72–82.
30. Bengtsson A, Eklundh J-O. Shape representation by multiscale contour approximation. IEEE Trans Pattern Anal Mach Intell. 1991 Jan;13(1):85–93.
31. Wang Y, Peterson BS, Staib LH. 3-D brain surface matching based on geodesics and local geometry. Comput Vision Image Understand. 2003;89:252–271.
32. Davatzikos C. Spatial transformation and registration of brain images using elastically deformable models. Comput Vision Image Understand. 1997;66(2):207–222. [PubMed]
33. Tagare H. Shape-based non-rigid correspondence with applications to heart motion analysis. IEEE Trans Med Imag. 1999 Jul;18(7):570–579. [PubMed]
34. Shi P, Robinson G, Chakraborty A, Stain L, Constable R, Sinusas A, Duncan J. Computer Vision, Virtual Reality, and Robotics in Medicine. Vol. 905. New York: Springer Verlag; 1995. A unified framework to assess myocardial function from 4D images; pp. 327–340. Lecture Notes Computer Science.
35. Rangarajan A, Chui H, Bookstein F. Information Processing in Medical Imaging (IPMI’97) Vol. 1230. LNCS; 1997. The softassign procustes matching algorithm; pp. 29–42.
36. Papademetris X, Jackowski AP, Schultz RT, Staibl LH, Duncan JS. MICCAI’03. Vol. 2879. LNCS; 2003. Computing 3-D non-rigid brain registration using extended robust point matching for composite multisubject fmri analysis; pp. 788–795.
37. Shen D, Davatzikos C. HAMMER: Hierarchical attribute matching mechanism for elastic registration. Proc IEEE Workshop Math Methods Biomed Image Anal. 2001:29–36.
38. Rosner B. Fundamentals of Biostatistics. New York: Duxbury; 1995.
39. Shaffer JP. Multiple hypothesis testing. Ann Rev Psych. 1995;46:561–584.
40. Bonferroni C. Studi in Onore del Professore Salvatore Ortu Carboni. Rome, Itlay: 1935. Il calcolo delle assicurazioni su gruppi di teste; pp. 13–60.
41. Bonferroni C. Teoria statistica delle classi e calcolo delle probabilitá Pubblicazioni del R Istituto Superiore di Scienze Economiche e Commerciali di Firenze. 1936;8:3–62.
42. Perneger T. What’s wrong with bonferroni adjustments. Brit Med J. 1998;316:1236–1238. [PMC free article] [PubMed]
43. Pantazis D, Leahy RM, Nichols TE, Styner M. Surface-based morphometry using a non-parametric approach. IEEE Int Symp Biom Imag (ISBI) 2004:1283–1286.
44. Adler RJ. The Geometry of Random Fields. New York: Wiley; 1981.
45. Worsley K, Evans A, Marrett S, Neelin P. A three–dimensional statistical analysis of CBF activation studies in human brain. J Cereb Blood Flow and Meta. 1992;12:900–918. [PubMed]
46. Worsley K, Marrett S, Neelin P, Vandal A, Friston K, Evans A. A unified statistical approach for determining significant signals in images of cerebral activation. Human Brain Mapp. 1996;4:58–73. [PubMed]
47. Worsley K, Andermann M, Koulis T, MacDonald D, Evans A. Detecting changes in nonisotropic images. Human Brain Mapp. 1999;8:98–101. [PubMed]
48. Taylor J, Adler R. Euler characteristics for Gaussian fields on manifolds. Ann Probab. 2003;31(2):533–563.
49. Boothby WM. An Introduction to Differential Manifolds and Riemannian Geometry. New York: Academic; 1986.
50. Levy B, Petitjean S, Ray N, Maillot J. Least squares conformal maps for automatic texture atlas generation. Special Interest Group on Computer Graphics—SIGGRAPH’02. 2002 Jul;:362–371.
51. Graphite. 2003. Online Available:
52. Blumberg HP, Kaufman J, Martin A, Whiteman R, Gore J, Chamey D, Krystal J, Peterson B. Amygdala and hippocampus volumes in adolescents and adults with bipolar disorder. Arch Gen Psychiatry. 2003;60:1201–1208. [PubMed]
53. A. P. Association. Diagnostical and Statistical Manual of Mental Disorders. Fourth. Washington, DC: Amer. Psychiatric Assoc.; 1994.
54. Ono M, Kubik S, Abernathey CD. Atlas of the Cerebral Sulci. New York: Thieme Med.; 1990.
55. Lerch J, Evans A. Cortical thickness analysis examined through power analysis and a population simulation. NeuroImage. 2005;24:163–173. [PubMed]
56. Viola P, Wells WM. Alignment by maximization of mutual information. 5th Int Conf Comput Vision. 1995:16–23.
57. Bansal R, Staib LH, Wang Y, Peterson BS. Roc-based assessments of 3-D cortical surface-matching algorithms. NeuroImage. 2005 [PubMed]
58. Christensen G, et al. 3-D Brain Mapping Using a Deformable Neuroanatomy. 1994;39:609–618. [PubMed]
59. Worsley KJ, Taylor JE, Tomaiuolo F, Lerch J. Unified univariate and multivariate random field theory. NeuroImage. 2004;23:S189–195. [PubMed]
60. Papoulis A. Probability, Random Variable, and Stochastic Processes. 3. New York: McGraw–Hill; 1991.
61. Thompson PM, Hayashi KM, Simon SL, Geaga JA, Hong MS, Sui Y, Lee JY, Toga AW, Ling W, London ED. Structural abnormalities in the brains of human subjects who use methamphetamine. J Neurosci. 2004;24(26):6028–6036. [PubMed]