NeuroImage 56 (2011) 961–973

Contents lists available at ScienceDirect

NeuroImagej o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / y n i m g

Local cortical surface complexity maps from spherical harmonic reconstructions☆
Rachel A. Yotter a,⁎, Igor Nenadic a, Gabriel Ziegler a, Paul M. Thompson b, Christian Gaser aab

Department of Psychiatry, Friedrich-Schiller University, D-07743 Jena, Germany
Laboratory of Neuro Imaging, Department of Neurology, UCLA School of Medicine, Los Angeles, CA 90095, USAa r t i c l e

i n f oArticle history:
Received 8 October 2010
Revised 1 February 2011
Accepted 1 February 2011
Available online 17 February 2011
Keywords:
Complexity
Shape analysis
Gyriﬁcation
Morphology
Schizophrenia
MRI

a b s t r a c t
Altered cortical surface complexity and gyriﬁcation differences may be a potentially sensitive marker forseveral neurodevelopmental disorders. We propose to use spherical harmonic (SPH) constructions tomeasure cortical surface folding complexity. First, we demonstrate that the complexity measure is accurate,by applying our SPH approach and the more traditional box-counting method to von Koch fractal surfaces
with known fractal dimension (FD) values. The SPH approach is then applied to study complexity differences
between 87 patients with DSM-IV schizophrenia (with stable psychopathology and treated with
antipsychotic medication; 48 male/39 female; mean age = 35.5 years, SD = 11.0) and 108 matched healthy
controls (68 male/40 female; mean age = 32.1 years, SD = 10.0). The global FD for the right hemisphere in theschizophrenia group was signiﬁcantly reduced. Regionally, reduced complexity was also found in temporal,
frontal, and cingulate regions in the right hemisphere, and temporal and prefrontal regions in the left
hemisphere. These results are discussed in terms of previously published ﬁndings. Finally, the anatomical
implications of a reduced FD are highlighted through comparison of two subjects with vastly different
complexity maps.
© 2011 Elsevier Inc. All rights reserved.

Introduction
In clinical neuroscience, the brains of individuals with a neurological or psychiatric disorder are often found to be structurally
different from the brains of healthy subjects. Furthermore, patterns ofsystematic structural differences have been found to be consistent for
subjects with a given disorder, such that the pattern of differences
depends on the speciﬁc disease. Detection of characteristic neuroanatomical patterns in patients with psychiatric disorders is likely tolead to a variety of advances in diagnosis, treatment, and preventionof psychological illness. For instance, it is at least logically possiblethat early diagnosis may be made in a prodromal phase, before theonset of symptoms, allowing preventative measures to be taken sothat no symptoms ever appear. Since the introduction of magnetic
resonance imaging (MRI) and other imaging techniques that allow invivo measurement of brain structure, ﬁnding speciﬁc structuraldifferences that characterize different psychiatric and neurogenetic
disorders is a real possibility.
Already, several methods have been proposed to ﬁnd potentiallyvery subtle structural differences (Thompson et al., 2004). Approaches
to measure differences may either be based on 3D volumetricassessments or 3D statistical mapping of gray or white matter☆ Funding sources: BMBF 01EV0709, BMBF 01GW0740.
⁎ Corresponding author at: Friedrich-Schiller University, Department of Psychiatry,
3 Jahnstrasse, 07743 Jena, Germany. Fax: +49 3641 934755.
E-mail address: rachel.yotter@uni-jena.de (R.A. Yotter).
1053-8119/$ – see front matter © 2011 Elsevier Inc. All rights reserved.
doi:10.1016/j.neuroimage.2011.02.007

structures, or may involve making measurements on 3D surface
models of a structure, such as the cortex or hippocampus (MacDonald
et al., 2000; Mangin et al., 1995; Morra et al., 2009; Wang et al., 2010;
Xu et al., 1999). Volume-based approaches are most appropriate formeasuring differences in gray matter density or in regional volumes,
such as for subcortical structures. Surface-based methods canmeasure cortical thickness or cortical folding patterns, as well as
shape or curvature measures derived from the imposed surface
coordinate system (Gutman et al., 2009; Luders et al., 2006; Lui et al.,2010). Our focus here is on analyzing cortical folding patterns.
Previously, the most common approach to measuring corticalfolding complexity was to use a metric such as the gyriﬁcation index(GI), which is deﬁned as the ratio of the inner surface size to the outersurface size of an outer (usually convex) hull. The GI can be measured
in two dimensions by examining cortical slices (Zilles et al., 1988), orin 3D by using a reconstructed surface mesh. However, the GI metricdepends on how the outer hull is deﬁned and how brains arenormalized to reduce the effect of brain size, which may result in large
differences in reported values across studies. There may also be an
inﬂuence of slicing direction for measures derived from 2D sections of
a 3D object – if more convolutions are apparent in a given plane ofsection, the measured complexity may tend to be higher. In otherwords, the method is not intrinsic. Another potential confound isnoise in the surface reconstruction, which could artiﬁcially inﬂate
surface area without corresponding to the underlying anatomy.
These drawbacks can be circumvented by using the fractal
dimension (FD), which does not rely on deﬁning an explicit outer962

R.A. Yotter et al. / NeuroImage 56 (2011) 961–973

hull (for a review, see Lopes and Betrouni, 2009). It has been proposedthat the brain is approximately a fractal (Kiselev et al., 2003), at least
over a limited range of scales (Jiang et al., 2008; Lee et al., 2004;
Mandelbrot and Blumen, 1989). This property has been evaluatedusing voxel-based information about the shape of the white mattervolume (Bullmore et al., 1994; Liu et al., 2003; Zhang et al., 2007;
Zhang et al., 2006). The FD can also be applied to measure corticalfolding complexity, and previous studies have shown signiﬁcantdifferences in regional FD for psychiatric disorders such as schizophrenia (Casanova et al., 1989; Narr et al., 2004; Narr et al., 2001),
obsessive-compulsive disorder (Ha et al., 2005), epilepsy (Free et al.,1996), Alzheimer's disease (King et al., 2009, 2010), multiple sclerosis(Esteban et al., 2007s, 2009), cerebellar degeneration (Wu et al.,2010), stroke (Zhang et al., 2008a), autism (Awate et al., 2008), andWilliams syndrome (Thompson et al., 2005), as well as for normal
development (Awate et al., 2010; Blanton et al., 2001; Kalmanti andMaris, 2007; Shyu et al., 2010; Wu et al., 2009), early-life blindness(Zhang et al., 2008b), and IQ (Im et al., 2006). Furthermore, somestudies report local FD or folding differences for genetic factors such as
22q11 deletion syndrome (Schaer et al., 2008), gender (Awate et al.,2009; Luders et al., 2004), and brain size (Toro et al., 2008).
Most calculations of the cortical surface FD rely on the boxcounting method, which calculates regional areas for progressively
lower sampling resolutions. Since the number of vertices steadilydecreases, the position of these vertices can have a large impact on the
FD metric and can potentially overlook relevant cortical folding
information (Free et al., 1996; King et al., 2009). This concern can beaddressed by aligning sulci across subjects to approximate the same
cortical location for each vertex for all subjects (Thompson et al.,1996). However, alignment is a complicated endeavor that often
requires manual delineation of areas that correspond geometrically
across subjects, such as gyral landmarks or cortical sulci. Anotherpossibility is to use a series of dilated disks of different sizes to cover
the surface, reducing but not eliminating orientation dependencies(Free et al., 1996).
Rather than using box-counting approaches, we propose to extract
surface complexity information using spherical harmonic (SPH)
reconstructions (Yotter et al., 2010). By using SPH reconstructions,the number of vertices remains the same for all reconstructedsurfaces, thus reducing the inﬂuence of individual vertex alignment
and avoiding the need to progressively re-grid the surface, which
incurs error from interpolation. Furthermore, analysis of structuralcharacteristics in some clinical disorders may be improved byinvestigating the pattern of regional differences rather than a single
global metric. For instance, it has been shown that patternclassiﬁcation techniques can enhance the speciﬁcity of morphometricﬁndings (Davatzikos et al., 2005; Kawasaki et al., 2007; Soriano-Mas
et al., 2007; Sun et al., 2009; Yushkevich et al., 2005).
One advantage of this approach is the use of information frommultiple brain regions when developing a characteristic neuroanatomical signature. This may be helpful in studies of mental illnesses,
e.g., schizophrenia, which we examine here. Using SPH-derived
reconstructions, we will demonstrate that it is possible to calculate
a local FD for each vertex within the reconstruction, which is the ﬁrst
step towards obtaining accurate neuroanatomical signatures of
disease states.
This paper is divided into two parts. In the ﬁrst section, we
compare our SPH-derived FD values to those obtained using the more
traditional box-counting approach. We then compare both measures
to the ground truth. We accomplish this by analyzing “von Koch”
surfaces with known FD values. Furthermore, we analyze how robust
the two measures are, using a brain surface rotated around each of thethree axes. In the second section, we then apply our approach to datafrom schizophrenia and control subjects, as an illustrative application,and to assess evidence for regional and local differences in corticalcomplexity. The ﬁndings from this analysis are compared topreviously published results from other methods. Finally, we discuss
the implications of the cortical complexity value in clinical and structural terms.
Methods
Spherical harmonic approach to measuring complexity
Generally, FD is computed by ﬁnding the slope of a plot regressing
log(area) versus log(dimension), over a certain range of scales. In thisregression – also called a multi-fractal plot – the dimension indicatesthe scale of measurement, which is varied either by sub-sampling theobject, or by reducing the degree of the shape representation. Whenusing spherical harmonic reconstructions, the plot is modiﬁed to use
the maximum l-value (or degree, see Appendix) of the reconstruction,and the slope may be found by regressing log(area) versus log(max lvalue) (Fig. 1). Complexity can be calculated at different scales: global,
regional, and local. The global complexity is a single value for theentire hemisphere; regional or local values are a set of values forregions of interest or vertices, respectively.
To obtain these values, for each hemisphere we extracted thespherical harmonic coefﬁcients of the central surface up to a
maximum l-value of 1024.1 The spherically remapped points weretransformed into harmonic space using a modiﬁcation of the fast
Fourier transform (Kostelec et al., 2000). The Appendix gives details
on the spherical harmonic analysis. Finally, the surfaces were reparameterized using a low-distortion spherical mapping of the central
surface.
Fractal dimension is calculated by ﬁnding the slope of the linear
portion of the log-log plot of area versus maximum l-value. This linear
portion generally occurs when the total reconstructed surface area isbetween 40% and 80% of the full surface area. To reduce computation
time, for brain surfaces 10 separate reconstructions were chosen using
maximum l-values between 11 and 29. These l-values achieved
optimal results for brain surfaces, resulting in surface areas that varied
from approximately 40% to 75% of the maximum surface area (Yotter
et al., 2010). For the von Koch surfaces, the linear portion of the slope
occurred for maximum l-values between 5 and 8; for the results
shown, all four reconstructions were used to calculate the FD values.
Calculation of the global complexity value was obtained directlyfrom the total surface area of the reconstructions. For each subject, a
set of initial local values was also produced. First, complexity values
for each polygon were estimated by regressing the log–log plot of
normalized area of the polygon versus the maximum l-value. Second,
point-wise complexity values were obtained by averaging thecomplexity values of all neighboring polygons. This resulted in a setof unsmoothed raw local complexity values for each subject.
For both local and global measures, the area values for the
spherical harmonic reconstructions were normalized by the areavalues in the full-coefﬁcient reconstruction (l-value = 1024), whose
surface area deviated by an inﬁnitesimal amount from the original
surface area. Although the original surface area could be used, the fullreconstruction surface area was used for consistency across allcomplexity measures (local and global). For the local complexity
calculation, the surface area of each polygon was normalized using thearea of the matching polygon in this reconstruction. This normalization step makes the complexity value somewhat independent of thesurface area of the original surface; however, since smaller brains tend
to have higher levels of folding detail, it is not possible to completely
remove correlations between brain size and complexity. Furthermore,
the regression slope is independent of the original surface area,
1
The maximum l-value does not signiﬁcantly affect the complexity calculation if thefull reconstruction retains most of the details of the original surface mesh (i.e., greaterthan approximately 128). The only effect is a shift in the x-intercept due to a different
surface area for the fully reconstructed mesh.

R.A. Yotter et al. / NeuroImage 56 (2011) 961–973

963

Fig. 1. Fractal dimension is found by ﬁnding the slope of a logarithmic plot of surface area versus the maximum l-value of the reconstruction (a measure of the bandwidth of
frequencies used to reconstruct the surface shape; shown here in red font for a range of bandwidths). Surface areas are normalized by the original surface area. A linear
approximation is reasonable over this range of scales.

which is reﬂected simply as a shift in the x-intercept of the linear
approximation.
Before comparing values across subjects, local values must ﬁrst be
represented in a common coordinate system, which is achieved via reparameterization of the unsmoothed raw values. This was accomplished using the registered spherical mapping for each subject (xh.
sphere.reg), then re-parameterizing the local FD values using the
fsaverage spherical mesh included in FreeSurfer. These values werethen smoothed using a 30-mm Gaussian heat kernel (Chung et al.,2005), where the width of the smoothing ﬁlter was chosen according
to the matched ﬁlter theorem. That is, the spatial frequency of thesulcal/gyral pattern suggests a ﬁlter that optimally enhances featuresin the range of the distance between sulci and gyri, which is about 20–
30 mm.
These smoothed, registered local complexity values were thenused for all statistical analyses. Regional values were obtained byaveraging the complexity values associated with a particular region ofinterest. A ﬂowchart of the steps taken to obtain complexity values isshown in Fig. 2.
Box-counting approach to measuring complexity
The deﬁnition of fractal dimension is related to the number of selfsimilar shapes versus the characteristic dimension, as follows:FD =log NðlÞ
;log 1 = lð1Þwhere l is the characteristic dimension and N(l) is the number of selfsimilar shapes with dimension l required to cover the original object.
Under certain conditions, namely when all objects are of uniform sizeand shape, this fractal dimension value may be calculated directly;
however, in most applications, there are no uniform shapes to be
counted. Fortunately, the number of self-similar objects N(l) can be
related to the total volume or area of the object. Since our application
focuses on surface shape analysis, we are interested in area. Generally,
the relationship between the number of objects N(l) and surface area
A can be related as follows:ANðlÞ∝ 2 :l

ð2ÞThus a complexity value can be calculated using the surface arearather than the number of self-similar objects, but the values will
deviate from the true Hausdorff dimension of the fractal surface if thefractal dimension were computed directly, by a factor of 1/log(l2).
For the box-counting method, the surface was re-sampled at
progressively lower resolution, and the surface areas were normalized, so that a value of 1 indicates that the surface area of the reparameterized surface is equivalent to the original surface area. We
opted to use regular Platonic solids to avoid non-uniform re-sampling
as much as possible. The FD was found by ﬁnding the slope of thelinear portion of a log–log plot in which the x-axis is the dimension (orcharacteristic length) and the y-axis is the surface area. For all surfaces(von Koch and cortical), the re-sampling resolutions within the linear
portion of the regression ranged from 80 to 5120 points. To derive thedimension, the x-axis was set to the square root of the inverse of thetriangle count, which is proportional to the length of a triangle edge.
To obtain local FD values, the down-sampled reconstruction wasthen re-sampled back into the original coordinate space. The FD value
was then calculated for each polygon by regressing log(area) versuslog(dimension), then point-wise values were calculated by averaging

Fig. 2. A ﬂowchart of the calculation process illustrates at which stages each set of complexity values are calculated. First, spherical harmonic reconstructions are derived from acentral surface mesh. Area values from these reconstructions are used to directly calculate the global FD values. The area information is also re-parameterized into a common
coordinate system across subjects; this data is then smoothed to derive the local FD values. Finally, regional values are averaged to obtain regional FD values.

964

R.A. Yotter et al. / NeuroImage 56 (2011) 961–973

the values of neighboring polygons. Regional FD values were obtainedby averaging the local FD values for all points within the region.
All attempts were made to obtain high-quality results for the boxcounting approach, including careful examination of the log–log plots,
the choice of dimensions, and the inclusion of a constant offset for all
values.

von Koch surfaces
To determine whether the FD values obtained from SPH-derived
reconstructions were accurate, we generated a variety of von Kochlike fractal surfaces of either tetrahedral or cubic structure. Thesesurfaces have an FD that is computable mathematically, as well as by
the SPH and box-counting algorithms. Because true von Koch surfaces
result in self-intersections, the surfaces needed to be modiﬁed to
avoid self-intersections, so that the Euler characteristic remained
equal to 2 and the surfaces could be projected to a sphere. This was
accomplished by simply reducing the length of the surface normal
used to calculate the projecting structure. Speciﬁcally, the surface
normal is l*sqrt(2/3) for tetrahedral and l for cubic von Koch surfaces,
while the modiﬁed surfaces had normal lengths of l*sqrt(2/3)/s and l/s,
respectively, where l is the characteristic dimension and s is what wewill term the projection scale.
Despite this modiﬁcation, the von Koch surfaces still had
characteristic FD values determined by measuring the slope of alog–log plot of total surface area versus characteristic dimension(Fig. 3). The accuracy of these measurements was validated bygenerating a true von Koch surface and conﬁrming that its FD value
matched the theoretical value (results not shown). These FD values
are the “true” values by which other approximation methods can beevaluated.
To test the ability to accurately discern local FD differences, a set of
surfaces was constructed such that each of the faces used a different
projection scale and thus had a unique fractal dimension.
In summary, the tetrahedral structures used six different projection scales (1.5, 2.0, 2.5, 3.0, 3.5, and 4.0) and three different characteristic dimensions or complexity levels, such that there were 18
surfaces plus three combination surfaces, which used four of the six
projection scales; the cubic structures used six different projection
scales and three complexity levels, resulting in 18 surfaces plus thethree combination surfaces. The point resolution of all surfaces waschosen so that the number of mesh points was approximately
equivalent to the number of points typically used for surface
representations of cortical surface structure. For both approaches,
the set of re-sampling resolutions or bandwidths was chosen so thatthe resultant global fractal dimension values were as closely matchedto the true fractal dimension values. Fitting of these values also
included a constant offset for each value in the set.

Statistical analysis
Statistical differences in global, regional, and local FD values wereassessed using a two-tailed t-test. Regression of age and sex in theschizophrenia study did not affect the results. The resulting t values
were thresholded at a p-value of p b 0.05, and local and regional values
were either corrected for multiple comparisons using the False Discovery Rate (FDR) method (Benjamini and Hochberg, 1995) or left
uncorrected.

Subjects
We studied 87 patients with a DSM-IV diagnosis of schizophrenia
and 108 healthy controls. The patients (n = 87; 48 male/39 female;mean age = 35.5 years, SD = 11.0) were recruited from the Department of Psychiatry in Jena and ﬁrst screened with a semi-structured
interview before being assessed by two psychiatrists establishing the
DSM-IV diagnosis. None of the patients had a second psychiatric,
neurological or major medical condition. All patients were inpatients.
None of them was in an acute episode of illness; all were in remitted
clinical state, showed stable psychopathology, and were on stableantipsychotic medication. Current psychopathology in patients wasassessed with the Scales for Assessment of Positive Symptoms (SAPS),
and Scales for Assessment of Negative Symptoms (SANS), which wasadministered by an experienced and trained clinical psychiatrist(Andreasen et al., 1995). Healthy controls (n = 108; 68 male/40
female; mean age = 32.1 years, SD = 10.0) were recruited from thecity of Jena and surrounding counties and matched to the patientswith regard to sex, age, and overall school/academic achievement.
They were screened to exclude a concurrent or past history ofpsychiatric, neurological, or major medical conditions using a semistructured interview. Further exclusion criteria for both samples
included a history of head trauma, concurrent or previous substancedependence or alcoholism, and learning disability (mental retardation). All subjects were right-handed, as scored from the short versionof the Edinburgh Handedness Scale (Oldﬁeld, 1971). All participants
gave written informed consent to a study protocol approved by the
Ethics Committee of the Friedrich-Schiller University of Jena. Detailsof this patient group may be found in Nenadic et al. (2010); a subset of
this sample (n = 12) were excluded due to poor surface reconstruction quality, which was based on visual inspection and rating by anexpert (G. Z.).

Fig. 3. Modiﬁed von Koch fractal surfaces with either a cubic (top row, left) or tetrahedral (bottom row, left) topology. Each set includes surfaces with progressively higher levels ofdetail, obtained by inserting self-similar shapes at lower dimensions. A log–log plot of surface area versus characteristic dimension (or width of the smallest feature) for three fractal
levels results in a linear ﬁt (right); the slope of this line is the characteristic FD.

R.A. Yotter et al. / NeuroImage 56 (2011) 961–973

Imaging protocol
We obtained a high-resolution structural MRI for each subject on a
1.5-T Phillips Gyroscan ASCII system using a T1-weighted sequence,
including 256 sagittal slices covering the entire brain (TR = 13 ms,
TE = 5 ms, 25° ﬂip angle, ﬁeld of view [FOV] = 256 mm, voxel
dimensions = 1 × 1 × 1 mm3) for all subjects. Foam pads were used,
where appropriate, to further restrict head movement. Prior to image
processing, each image was checked visually for artifacts. In addition,
we used an automated tool to detect outliers implemented in the
VBM8 package. All scans passed both manual and automated quality
checks.
Cortical surface extraction
Image volumes passed through a number of preprocessing steps
using mostly automated procedures included in the FreeSurfersoftware suite, version 4.5 (http://surfer.nmr.mgh.harvard.edu).
Details of the processing steps have been described previously (Daleet al., 1999; Fischl et al., 1999a, 1999b). In brief, images were preprocessed using intensity normalization and skull stripping, followed bynormalizing the head position along the commissural axis andlabeling cortical and subcortical regions. Images were then segmented, as follows (Fischl and Dale, 2000). First, images were rigid-bodyregistered to a probabilistic brain atlas, followed by non-linear
morphing to the atlas. Depending on the probability that a given
location is of a particular tissue class, the intensity of the image at thelocation, and the local spatial conﬁguration of the location in relation
to the labels, each voxel is then assigned to one speciﬁc tissue class(gray matter, white matter, CSF, or background).
A white matter surface was derived from the white matter tissue
segmentation map using the marching cubes algorithm, followed by
topology correction. By outwardly deforming the white mattersurface, the pial surface was also generated (Dale et al., 1999; Fischl
and Dale, 2000). As a ﬁnal step, the pial and white matter surfaces
were averaged together, vertex-by-vertex, to construct a central
surface. The central surface meshes were then used as the input forcomplexity analysis using our spherical harmonic approach. Before
comparing local and regional values across subjects, the local FD
values were re-sampled such that all values were registered to the
FreeSurfer average template using the default registration algorithm(Fischl et al., 1999b).
Validation
Results
Fractal dimension measurement accuracy for methods using boxcounting and spherical harmonics. The fractal dimension for the vonKoch surfaces was modiﬁed by adjusting the scale used for the
projecting structures, such that a projection scale value of 1.0 is a true
von Koch surface and values larger than that is the ratio by which the
projection scale has been reduced. This eliminates self-intersectionsin the structures and offers an approach to modify the fractal dimension values without large topological changes. Note that if the numberof triangles or squares is plotted against the characteristic dimension,
this process yields a Hausdorff dimension of 2.5849 for the tetrahedralsurfaces and 2.3347 for the cubic surfaces, independent of the projection scale.
When comparing reconstructions of level-3 surfaces from the boxcounting method (Fig. 4) to the SPH approach (Fig. 5), there are
striking similarities. Both lose ﬁne detail as the re-sampling resolution
or the maximum l-value decreases. Because the SPH approach does
not down-sample the mesh resolution, the reconstructions are moreregular and of higher quality than those created from the box-

965

counting method. However, the plots of surface area versus characteristic dimension or maximum l-value are similar.
For all of the modiﬁed von Koch surfaces, the FD measuresobtained from the SPH approach were more accurate than thosederived from the box-counting method, at both the global (or overall)and regional levels (Fig. 6). For tetrahedral topologies, the resultant
FD values are virtually identical to the “true” FD values, and themeasurement variability over different levels of complexity is low.
Local variations in fractal dimension were simulated using vonKoch fractal surfaces that had a different projection scale for eachsurface face. At the regional level, the SPH approach effectively measures the differences between the complexity values of each face,
while the box-counting approach produces large errors, to the extent
that some faces having lower complexity appear to have higher
complexity. Because the spherical harmonic approach re-samples thespherical mesh at very high resolution (rather than down-sampling aswith box-counting), the local fractal dimension maps are much moreregular than those obtained from the box-counting approach (Fig. 7).
Measures on a rotated cortical surface. Validation of the measures on
cortical surfaces is difﬁcult since the approach is new and there are no
cortical surfaces with precisely known fractal dimension values.
However, it is possible to manipulate a cortical surface in a prescribedway, measure the change in the resulting values, and assess whether
the effects coincide with expectation. One of the simplest manipulations is a rotation of the surface around rotational axes, e.g., pitch,yaw, and roll.
In theory, the complexity measure for a brain surface should
remain constant, even if the surface is rotated. However, the global
fractal dimension value for the box-counting approach varies widelywith respect to the angle by which the surface is rotated (Fig. 8). Thisis clearly undesirable, as the summary measure of the properties of an
object should not heavily depend upon its orientation. Furthermore,
the correlation between the local FD values of rotated and un-rotatedsurfaces for the box-counting method only approaches a value of 1 fora rotation of 180° (Fig. 9). In contrast, the spherical harmonic derived
values for global FD remain more constant, and the local FD
correlation between all rotated surfaces and the un-rotated surface
was never lower than 1.0.
Discussion
Validation of SPH approach using von Koch surfaces. Analysis of the vonKoch surfaces, whose fractal dimensions are known, demonstrated
that the complexity values obtained from spherical harmonicreconstructions are generally more accurate than those obtained
from the box-counting approach. Because they are less severely
affected by re-sampling, interpolation, and object alignment, theresultant local complexity maps from SPH reconstructions tend to be
more regular than those obtained from box-counting. This presumably makes measures more accurate and less prone to error at both
the regional and global levels.
Further evidence of the improved performance of the SPH
approach was obtained by simply rotating a cortical brain surface
around the three cardinal axes. The large measurement differences
imply that box-counting will only give measures that are comparable
across subjects if brain surfaces are ﬁrst carefully registered to eachother, such that all re-sampled points in the down-sampled meshes
roughly correspond to each other for all brain surfaces. Due to thedependence of the results on orientation of the brain, it can bepresumed that box-counting would also be highly dependent upon
the accuracy of the registration method employed. Conversely, thespherical harmonic approach is less dependent upon the orientationof the brain surfaces, and thus does not require careful registration inorder to extract valid complexity values.

966

R.A. Yotter et al. / NeuroImage 56 (2011) 961–973

Fig. 4. Box-counting can measure the fractal dimension. Here we show a level-3 cubic von Koch surface with a projection scale of 2.0 re-sampled at various resolutions; the numbers
indicate the triangle count in the re-sampled mesh. Below, we plot the normalized surface area versus dimension on a log–log scale, for detail levels 0–3. Fractal dimension can be
estimated by measuring the slope over a range of scales that results in an approximately linear function. The area of the re-sampled surface is divided by the surface area of theoriginal mesh.

Application to schizophrenia
Results
We detected no signiﬁcant group differences between patients andcontrols for the global FD values in the left hemisphere (control:
2.5820 ± 0.0002; schizophrenia: 2.5789 ± 0.0003; p = 0.4252), butthere was a signiﬁcant difference in the right hemisphere (control:
2.5871 ± 0.0003; schizophrenia: 2.5729 ± 0.0003; p = 0.0007). Interestingly, there were also signiﬁcant differences between groups forsurface area in both hemispheres (left hemisphere: control: 923.70
cm2; schizophrenia: 895.40 cm2; p = 0.0234; right hemisphere:
control: 925.58 cm2; schizophrenia: 894.91 cm2; p = 0.0147), but
not between right and left hemispheres within groups (control:
p = 0.8776; schizophrenia: p = 0.9689).
Regionally, the FD values in the schizophrenia group weresigniﬁcantly lower than those for controls, bilaterally in the precentral
gyrus; in the caudal anterior cingulate gyrus and frontal pole in theleft hemisphere; and the corpus callosum, lingual gyrus, and superiorparietal lobe in the right hemisphere (Table 1). Patients had
signiﬁcantly higher regional surface complexity in the cuneus,
isthmus cingulate gyrus, lateral orbitofrontal gyrus, paracentral
gyrus, posterior cingulate cortex, and transverse temporal gyrus in
the left hemisphere. These regions are highlighted in Fig. 10.
An unusual feature of the values in Table 1 is that two values are
not within the range of 2.0–3.0. Most deﬁnitions of fractal dimension
imply that the value must be lower than the embedding dimension,
which is 3.0 for a 3-dimensional surface (Mehaute, 1991). However,the deﬁnition of fractal dimension used here is not the exact Hausdorff
dimension, thus can result in values outside this range. In the righthemisphere, the frontal pole region deﬁes this rule. It is a highly
folded gyral region that contains a relatively small number ofpolygons, and this deviation may be due to loss of resolution in thatarea.
The local FD maps revealed that the lower overall complexity inthe right hemisphere was attributable to differences in the parietal
and temporal lobes (Fig. 11). Local complexity was also lower for the
prefrontal and temporal cortex in the left hemisphere, but only
prefrontal ROIs in the left hemisphere showed signiﬁcantly lower FD
in the patient group.
Discussion
Surface complexity in schizophrenia. For the whole hemisphere,
complexity analysis of the schizophrenia patient group revealed a
global reduction in surface complexity in the right hemisphere only.
This lateralization corresponds with some previous studies that
detected lateralized gray matter reductions in the right hemisphere(Pantelis et al., 2003). Even so, some studies detected gray matter
deﬁcits only in the left hemisphere (Honea et al., 2005; Hulshoff Pol
et al., 2006). Furthermore, prior studies of cortical folding found either
bilateral or left-hemisphere reductions in GI for the schizophrenia
group (Sallet et al., 2003) or no differences (Noga et al., 1995).
Measuring the FD of a skeleton representing the cortical folding
pattern revealed signiﬁcantly lower FD values for both hemispheres inR.A. Yotter et al. / NeuroImage 56 (2011) 961–973

967

Fig. 5. Fractal dimension may be computed from a series of SPH low-pass ﬁltered reconstructed surfaces. Top two rows: we show a series of level-3 cubic surfaces with projection
scale of 2.0 that have been low-pass ﬁltered; the numbers indicate the maximum l-value used in the ﬁlter. Below, we plot the normalized surface area versus the maximum l-value on
a log–log scale, for detail levels 0–3. Fractal dimension can be estimated by measuring the slope over a range of l-values that results in an approximately linear function. The surface
area of the reconstructed surface is divided by the surface area of the original mesh.

Fig. 6. Global and regional fractal dimension measurements obtained from spherical harmonic reconstructions are more accurate than those obtained from the box-counting
approach. The average FD values – shown with standard errors of the means (SEM) – are averaged from three levels of complexity for each projection scale. Regional (ROI) values are
obtained from averaged local FD values from a combined surface in which each face uses a different projection scale.

968

R.A. Yotter et al. / NeuroImage 56 (2011) 961–973

Fig. 7. A map of local fractal dimension based on the box-counting (top row) and spherical harmonics (bottom row) methods was calculated for a von Koch surface with varying
projection scales on each surface. The level-3 surface is shown here, and the numbers represent the projection scale of each face, such that face 1.5 has the largest features. The
spherical harmonic mapping clearly marks structures with higher complexity with a larger fractal dimension value, while the box-counting approach depends strongly upon the resampled points.

the patient group (Ha et al., 2005). This was similar to the result foundusing box-counting estimation of FD for the boundary between GM
and WM (Bullmore et al., 1994). However, another study found anincrease in FD for both hemispheres (Sandu et al., 2008). Since thesecomplexity measures were of the GM/WM boundary, however, they
may not account for any possible gray matter abnormalities and may
have a different etiology.
The pattern of regional structural differences in schizophrenia
versus controls is even more heterogeneous across studies, both interms of gray matter reduction (Weinberger and McClure, 2002) andin gyriﬁcation index (Wheeler and Harper, 2007). However, the mostconsistently reported ﬁndings include gray matter reduction in thetemporal (Goldman et al., 2008; Honea et al., 2005; McDonald et al.,2006; Wright et al., 2000) and frontal lobes (Honea et al., 2008;
Thompson et al., 2009), along with an increase in the size of the lateral
ventricles (Steen et al., 2006). Regional complexity analysis revealed
signiﬁcant decreases in complexity for the frontal lobe areas in bothhemispheres, which is similar to previous studies reporting bilateralcomplexity differences in the frontal cortex (Narr et al., 2001);
however, another study showed only right hemisphere GI increases in
the prefrontal cortex (Vogeley et al., 2000; Vogeley et al., 2001).

Fig. 8. The global fractal dimension value can depend upon the orientation of a brain
surface. Rotating one brain surface around the three axes results in relatively large
errors in global fractal dimension values for the box-counting approach, while theerrors remain vanishingly small for the proposed spherical harmonic approach.

The evaluation of the clinical results is further complicated by thefact that the underlying pathogenetic causes for schizophrenia areheterogeneous. For example, disturbances in neurodevelopment
might be inﬂuenced by genetic factors or obstetric complications,which might lead to different patterns of FD change compared tothose patients who have little genetic risk but suffer a “second hit,”e.g., during adolescence. It is therefore necessary to conceive ofmultiple subsyndromes within schizophrenia. Hence, when all
patients are analyzed as a single group, any structural abnormalities
potentially related to speciﬁc subsyndromes are lost or greatly
diluted. Furthermore, our local FD analyses were not corrected for
multiple comparisons, and were provided solely to demonstrate that
this methodological approach can calculate highly localized complexity values. To further complicate the analysis, direct comparison to GI
or gray matter volume reduction cannot be made, since thecomplexity is a different type of measure (Thompson et al., 2005).
This last point is discussed below in more detail.
Dividing patients according to the subtypes deﬁned in DBM-IIIR or
DSM-IV diagnostic criteria may be limited by the fact that the crosssectional psychopathology could vary from one psychotic episode toanother, hence obscuring the identiﬁcation of the causal pathogenetic
mechanism (Gaser et al., 2004). However, using voxel-based
morphometry, we have previously demonstrated that there isconsiderable brain structural variation related to symptom proﬁles,
especially for the prefrontal cortex, thalamus, and temporal lobe
(Nenadic et al., 2010).
The heterogeneity of these ﬁndings suggest that there is further
improvements to be made in both deﬁning schizophrenia as a
psychiatric disease and in developing a standardized method tomeasure potentially small structural differences. It could be argued
that the contradictory results are entirely due to methodologicaldifferences, diverse patient groups, or small sample sizes. It is also
possible, as in any study, that unmodeled confounds mediate thedifferences. The algorithm proposed in this paper may overcome
some of the methodological difﬁculties for measuring cortical surface
complexity. It computes localized maps of complexity values that may
be used to generate surface patterns. When combined with machinelearning approaches, these patterns may be able to accurately classify
certain disease states (Sun et al., 2009).
Interpretation of surface complexity measures. In order to clinicallyinterpret a difference in measured complexity, it is important to beable to relate the value to a more conventionally accepted structuraldifference. Unlike gray matter differences, cortical thickness, orR.A. Yotter et al. / NeuroImage 56 (2011) 961–973

969

Fig. 9. Despite rotation along axes, the correlation between the derived local fractal dimension values and the original values from the un-rotated surface should remain (nearly)
perfect. For the spherical harmonic approach, the correlation is unaffected by the rotation and is never lower than 1. However, only a rotation of 180° results in a correlation
approaching a value of 1 for the box-counting approach.

cortical folding differences, surface complexity, especially at a local
level, measures a structural aspect that is perhaps somewhat lessintuitive or less commonly considered. It helps to bear in mind that acomplexity measure does not directly measure the intuitive meaningof complexity (e.g., more detail), but is a characteristic of the surface
shape. This is true at the global level, and an assumption made in theabove analyses is that these concepts apply to the local level as well.
Table 1
Mean annotated FD for schizophrenia and control subjects, over regions deﬁned using
registration algorithms included in the FreeSurfer software suite.
Region of interest

Bankssts
Caudal anterior cingulate
Caudal middle frontal
Corpus callosum
CuneusEntorhinal
Fusiform
Inferior parietal
Inferior temporal
Isthmus cingulate
Lateral occipital
Lateral orbitofrontal
LingualMedial orbitofrontal
Middle temporal
Parahippocampal
Paracentral
Parsopercularis
Parsorbitalis
Parstriangularis
Pericalcarine
Postcentral
Posterior cingulate
Precentral
Precuneus
Rostral anterior cingulate
Rostral middle frontal
Superior frontal
Superior parietal
Superior temporal
Supramarginal
Frontal pole
Temporal pole
Transverse temporal⁎ p b 0.05, uncorrected.
⁎⁎ p b 0.05, FDR.

Left hemisphere

Right hemisphere

Control

Schizo.

Control

Schizo.

2.81532.29522.7336
2.17582.65122.90312.52982.69152.4693
2.18482.50522.26682.59982.42082.67832.90802.47772.76052.96482.77672.58022.77102.43592.73022.5641
2.11292.58542.38072.66692.78732.62622.79542.81072.64222.77922.2453 ⁎2.73152.16342.6982 ⁎2.92132.50952.70352.46592.2214 ⁎2.50742.2838 ⁎2.58252.44762.68692.91342.5018 ⁎2.74872.95952.74542.55232.75352.4832 ⁎2.7100 ⁎2.87982.06252.67582.54802.72092.33852.47282.64452.52672.74892.47962.35582.68442.39822.87992.53542.90152.18452.29022.53942.85692.49792.46822.58442.71352.31202.58502.75302.73262.44792.4442
3.18062.55001.95812.84692.04152.66122.5117 ⁎2.75432.34772.46542.62692.51522.74482.48272.35142.6466 ⁎2.5463
2.09762.59432.3631 ⁎2.65602.76942.62212.7510 ⁎2.82252.7777 ⁎2.39252.88922.50182.90112.19652.28392.54022.84362.51042.45622.5687 ⁎2.68812.31512.57762.74412.6679 ⁎⁎2.43892.43973.21772.57801.9695

Indeed, quantitative support for this assumption is provided by the
relative accuracy of regional fractal dimension differences for the
combined von Koch surfaces.
A fractal is a structure that is self-similar across a range of scales,
and the complexity analysis roughly corresponds to how “spaceﬁlling” the fractal surface is. However, a plot of local complexity
measures for two brains – one with a large FD and one with a small FD– demonstrates that the relationship between self-similarity, spaceﬁlling, and cortical folding pattern is not particularly straightforward
at the local level (Fig. 12). Generally, regions with a high FD value
appear to be more periodically spaced, such as a sine wave withregular peaks and troughs. This may be because the more periodicallyspaced structures also tend to ﬁll more space over the range of scales
examined for derivation of complexity values. As support for thisconcept, the Pearson correlation coefﬁcient between the summed
fractal dimension for both hemispheres and the ratio of surface area
to volume for the entire brain is signiﬁcantly positive (r = 0.4906,
n = 195, p = 1.67 × 10−13). Further evidence comes from examining
the change in fractal dimension for each surface. In the above analysis,Fig. 10. The highlighted regions have signiﬁcantly different FD values, on average,
between schizophrenia and control subjects. The FD values were signiﬁcantly different
in the precentral gyrus, caudal anterior cingulate gyrus, frontal pole, cuneus, isthmus
cingulate gyrus, lateral orbitofrontal gyrus, paracentral gyrus, posterior cingulate
cortex, and transverse temporal gyrus in the left hemisphere; and the precentral gyrus,
corpus callosum, lingual gyrus, and superior parietal lobe in the right hemisphere. Bluecolors indicate that the mean values for the schizophrenia group are lower than thosefor the control group.

970

R.A. Yotter et al. / NeuroImage 56 (2011) 961–973

Fig. 11. Local average FD values are mapped onto a representative cortical surface model for control and schizophrenia groups. Left and right hemispheres are shown. A map of mean
complexity differences (3rd column) between groups and p-values showing the local signiﬁcance of the inter-group difference are shown (4th column). Vertices are highlighted if(uncorrected) p b 0.05. For both the mean group difference and p-value maps, blue colors indicate that the mean values for the schizophrenia group are lower than those for thecontrol group.

the central surface was used rather than the white matter or pial
surfaces. The fractal dimensions of these surfaces were relatively
stably offset from the fractal dimension value of the central surface,such that the pial surfaces had higher fractal dimension values and thewhite matter surfaces had lower values (Fig. S2). Again, this may be
due to the more (or less) space-ﬁlling structure of the pial (or white)surfaces as compared with the central surface.
The measured FD value for an object depends, to some extent, on
the range of scales examined, and may depend upon the resolution of
the re-sampling, or the largest l-value in the case of sphericalharmonic reconstructions. A very high resolution re-sampling could
retain high-frequency noise that may be due more to noise in the scanor the reconstruction protocol than biological aspects of brain
structure. For this reason, we deliberately used fairly low-frequency
reconstructions such that the underlying base shape determines
the resulting FD value. This avoidance of very high frequency information should make our approach less dependent on scanner or noise.
When this was tested, the FD value was not completely independentof scanner, as intensity differences are reﬂected to some extent
in structural differences in the reconstructed surface mesh. ScannerFig. 12. These two example surfaces have very different cortical complexity values. To visualize what a complexity value means in terms of structure, two low-pass ﬁltered
reconstructions are shown, as well as a local complexity map. In regions of high complexity, the underlying brain structure is more regular, similar to a sine wave, which is mostapparent in the temporal and frontal lobes.

R.A. Yotter et al. / NeuroImage 56 (2011) 961–973

biases may inﬂuence pre-processing steps such as segmentationand surface reconstruction. If that is the case, these biases are theninevitably reﬂected in the resulting FD value, but should not induce
systematic differences between groups of patients and controls.
However, the pattern of local FD values for the same brain scannedwith different scanners appears to be highly reproducible (Fig. S3).
Conclusion
In this paper, we proposed a new approach to measure corticalsurface complexity from brain surface meshes. The approach relies onhighly smoothed spherical harmonic reconstructions to derivecomplexity values, using a new deﬁnition of fractal dimension as
the regression of log(surface area) versus log(max l-value). Despite thenew deﬁnition, the fractal dimension values can accurately measuredifferences in global (or overall) fractal dimension, and is demonstrably more accurate than the box-counting method. The new
deﬁnition is also consistent for this method, thus allowing potential
comparison of results from multiple studies.
Because the number of vertices in each construction remains thesame, calculation of local (vertex-wise) complexity values is straightforward. These local values can be used to obtain regional complexity
values by averaging the values from vertices contained within eachregion of interest. Using artiﬁcial von Koch surfaces, we demonstrated
that the regional values followed the trend for differences in fractal
dimension (i.e., a plot of projection scale versus fractal dimension
remained monotonic). Visual inspection of local FD maps on the vonKoch surfaces revealed a highly regular pattern.
As demonstrated in the rotated cortical surface, the box-counting
method must be preceded by careful feature alignment of the corticalsurfaces to attempt to map re-parameterized points to equivalentlocations across subjects (Thompson et al., 1996). Due to the highresolution re-sampling in the spherical harmonic approach, we could
bypass this step without loss of accuracy in global and local complexity measures. However, to minimize variability, subjects wereregistered to a common template before analysis.
Cortical surface complexity may be a valid clinical marker forschizophrenia (Narr et al., 2004), Williams syndrome (Gaser et al.,2006; Thompson et al., 2005), bipolar disorder (McIntosh et al., 2009),
and obsessive-compulsive disorder (Ha et al., 2005). When applied toa large sample of schizophrenia and control subjects, the schizophrenia subjects had signiﬁcantly lower cortical surface complexity in theright hemisphere. Regional and local analyses suggest that there aredifferences in the frontal, parietal, and temporal lobes, which should
be examined in more detail. These preliminary ﬁndings demonstratethat our approach can detect relatively small, local differences incomplexity.
Supplementary materials related to this article can be found online
at doi:10.1016/j.neuroimage.2011.02.007.
Acknowledgments
R. A. Y., G. Z., and C. G. were supported by grants BMBF 01EV0709
and BMBF 01GW0740 from the German Bundesministerium für
Bildung und Forschung. P.T. was supported by NIH grants EB008432,
EB008281, EB007813, HD050735. I.N. was partly supported throughthe EU grant: RTN, FP6; "EUTwinS," MRTN-CT-2006-035987.
Appendix. Spherical harmonic analysis
To analyze the harmonic content of a surface mesh, the ﬁrst
required step is to change the parameter space from Cartesian tospherical coordinates θ and ϕ, where θ is the co-latitude and ϕ is theazimuthal coordinate. The original spherical mapping is from thestandard FreeSurfer pipeline. Then, grid points are interpolated forequally sampled values of θ and ϕ for all members in the sets, suchthat there are 2B points per set, where B is the degree, or maximum l-971

value. For each regularly sampled point, the closest polygon on thespherical mapping is found. Within the closest polygon, a spatial
location for the interpolated vertex is approximated using barycentric
coordinates. The result is a regularly sampled spherical map in which
every point is associated with a coordinate that gives its location on
the original surface.
Once the surface mesh is re-parameterized, the harmonic content
of a spherical mesh may be obtained using normalized sphericalharmonics Ylm (θ, ϕ):
m

m

imϕ

Yl ðθ; ϕÞ = Pl ðcos θÞe;ð3Þ

where l and m are integers with |m| ≤ l, and Pm
l is the associatedLegendre function deﬁned by:mPl ðxÞ = l+m 
l
1 
2 m d
2
1−x 2 l + m x −1 :
l
dx
2 l!

ð4Þ

A square-integrable function f(θ,ϕ) on the sphere may be
expanded in the spherical harmonic basis such that:Bf ðθ; ϕÞ = ∑l∑

l = 0 m = −ljjYlm jj2 f̂ ðl; mÞ⋅Ylm ;−2

ð5Þ

where the coefﬁcients f̂(l,m) are deﬁned by f̂(l,m) = hf, Ylm i and theL2-norm of Ylm is given by:ðl + mÞ!
⋅:jjYlm jj2 = 2l 4π+ 1 ðl−mÞ!
−2

ð6Þ

It is possible to solve this system directly by ﬁnding the coefﬁcients
ﬁrst, but a more efﬁcient approach is to use a divide-and-conquer
scheme (Healy et al., 1996).
These coefﬁcients can then be low-pass ﬁltered, such that only
lower frequencies have non-zero values, and passed through an inverse
transform to produce a surface reconstruction. For FD calculations,ten reconstructions are produced using an upper l-value between 11
and 29.
References
Andreasen, N.C., Arndt, S., Alliger, R., Miller, D., Flaum, M., 1995. Symptoms ofschizophrenia: methods, meanings, and mechanisms. Arch. Gen. Psychiatry 52,
341–351.
Awate, S., Win, L., Yushkevich, P., Schultz, R., Gee, J., 2008. 3D cerebral corticalmorphometry in autism: increased folding in children and adolescents in frontal,
parietal, and temporal lobes. In: Metaxas, D., Axel, L., Fichtinger, G., Székely, G.
(Eds.), Medical Image Computing and Computer-Assisted Intervention – MICCAI2008. Springer Berlin, Heidelberg, pp. 559–567.
Awate, S., Yushkevich, P., Licht, D., Gee, J., 2009. Gender differences in cerebral corticalfolding: multivariate complexity-shape analysis with insights into handling brainvolume differences. In: Yang, G.-Z., Hawkes, D., Rueckert, D., Noble, A., Taylor, C.
(Eds.), Medical Image Computing and Computer-Assisted Intervention – MICCAI2009. Springer Berlin, Heidelberg, pp. 200–207.
Awate, S.P., Yushkevich, P.A., Song, Z., Licht, D.J., Gee, J.C., 2010. Cerebral cortical folding
analysis with multivariate modeling and testing: studies on gender differences andneonatal development. Neuroimage 53, 450–459.
Benjamini, Y., Hochberg, Y., 1995. Controlling the false discovery rate: a practical andpowerful approach to multiple testing. J. R. Stat. Soc. B Methodol. 57, 289–300.
Blanton, R.E., Levitt, J.G., Thompson, P.M., Narr, K.L., Capetillo-Cunliffe, L., Nobel, A.,
Singerman, J.D., McCracken, J.T., Toga, A.W., 2001. Mapping cortical asymmetry andcomplexity patterns in normal children. Psychiatry Res. Neuroimaging 107, 29–43.
Bullmore, E., Brammer, M., Harvey, I., Persaud, R., Murray, R., Ron, M., 1994. Fractal
analysis of the boundary between white matter and cerebral cortex in magnetic
resonance images: a controlled study of schizophrenic and manic-depressive
patients. Psychol. Med. 24, 771–781.
Casanova, M.F., Daniel, D.G., Goldberg, T.E., Suddath, R.L., Weinberger, D.R., 1989. Shapeanalysis of the middle cranial fossa of schizophrenic patients: a computerized
tomographic study. Schizophr. Res. 2, 333–338.
Chung, M.K., Robbins, S.M., Dalton, K.M., Davidson, R.J., Alexander, A.L., Evans, A.C.,
2005. Cortical thickness analysis in autism with heat kernel smoothing. Neuroimage 25, 1256–1265.
Dale, A.M., Fischl, B., Sereno, M.I., 1999. Cortical surface-based analysis: I. Segmentationand surface reconstruction. NeuroImage 9, 179–194.

972

R.A. Yotter et al. / NeuroImage 56 (2011) 961–973

Davatzikos, C., Shen, D., Gur, R.C., Wu, X., Liu, D., Fan, Y., Hughett, P., Turetsky, B.I., Gur, R.
E., 2005. Whole-brain morphometric study of schizophrenia revealing a spatiallycomplex set of focal abnormalities. Arch. Gen. Psychiatry 62, 1218–1227.
Esteban, F.J., Sepulcre, J., de Mendiz·bal, N.V., GoÒi, J., Navas, J., de Miras, J.R., Bejarano,B., Masdeu, J.C., Villoslada, P., 2007. Fractal dimension and white matter changes in
multiple sclerosis. Neuroimage 36, 543–549.
Esteban, F.J., Sepulcre, J., Ruiz de Miras, J., Navas, J., Velez de Mendizabal, N., Goni, J.,
Bejarano, B., Villoslada, P., 2009. Fractal dimension analysis of grey matter inmultiple sclerosis. J. Neurol. Sci. 282, 67–71.
Fischl, B., Dale, A.M., 2000. Measuring the thickness of the human cerebral cortex frommagnetic resonance images. Proc. Natl Acad. Sci. USA 97, 11050–11055.
Fischl, B., Sereno, M.I., Dale, A.M., 1999a. Cortical surface-based analysis: II. Inﬂation,
ﬂattening, and a surface-based coordinate system. Neuroimage 9, 195–207.
Fischl, B., Sereno, M.I., Tootell, R.B.H., Dale, A.M., 1999b. High-resolution intersubject
averaging and a coordinate system for the cortical surface. HBM 8, 272–284.
Free, S.L., Sisodiya, S.M., Cook, M.J., Fish, D.R., Shorvon, S.D., 1996. Three-dimensional
fractal analysis of the white matter surface from magnetic resonance images of thehuman brain. Cereb. Cortex 6, 830–836.
Gaser, C., Luders, E., Thompson, P.M., Lee, A.D., Dutton, R.A., Geaga, J.A., Hayashi, K.M.,
Bellugi, U., Galaburda, A.M., Korenberg, J.R., Mills, D.L., Toga, A.W., Reiss, A.L., 2006.
Increased local gyriﬁcation mapped in Williams syndrome. Neuroimage 33, 46–54.
Gaser, C., Nenadic, I., Volz, H.-P., Buchel, C., Sauer, H., 2004. Neuroanatomy of ‘hearing
voices’: a frontotemporal brain structural abnormality associated with auditory
hallucinations in schizophrenia. Cereb. Cortex 14, 91–96.
Goldman, A.L., Pezawas, L., Mattay, V.S., Fischl, B., Verchinski, B.A., Zoltick, B.,
Weinberger, D.R., Meyer-Lindenberg, A., 2008. Heritability of brain morphology
related to schizophrenia: a large-scale automated magnetic resonance imaging
segmentation study. Biol. Psychiatry 63, 475–483.
Gutman, B., Wang, Y., Morra, J., Toga, A.W., Thompson, P.M., 2009. Disease classiﬁcation
with hippocampal shape invariants. Hippocampus 19, 572–578.
Ha, T.H., Yoon, U., Lee, K.J., Shin, Y.W., Lee, J.-M., Kim, I.Y., Ha, K.S., Kim, S.I., Kwon, J.S.,
2005. Fractal dimension of cerebral cortical surface in schizophrenia and obsessivecompulsive disorder. Neurosci. Lett. 384, 172–176.
Healy, D.M., Rockmore, D.N., Moore, S.S.B., 1996. FFTs for the 2-Sphere-Improvements
and Variations. Dartmouth College.
Honea, R., Crow, T.J., Passingham, D., Mackay, C.E., 2005. Regional deﬁcits in brain
volume in schizophrenia: a meta-analysis of voxel-based morphometry studies.
Am. J. Psychiatry 162, 2233–2245.
Honea, R.A., Meyer-Lindenberg, A., Hobbs, K.B., Pezawas, L., Mattay, V.S., Egan, M.F.,
Verchinski, B., Passingham, R.E., Weinberger, D.R., Callicott, J.H., 2008. Is graymatter volume an intermediate phenotype for schizophrenia? A voxel-based
morphometry study of patients with schizophrenia and their healthy siblings. Biol.
Psychiatry 63, 465–474.
Hulshoff Pol, H.E., Schnack, H.G., Mandl, R.C.W., Brans, R.G.H., van Haren, N.E.M., BaarÈ,W.F.C., van Oel, C.J., Collins, D.L., Evans, A.C., Kahn, R.S., 2006. Gray and white matter
density changes in monozygotic and same-sex dizygotic twins discordant for
schizophrenia using voxel-based morphometry. Neuroimage 31, 482–488.
Im, K., Lee, J.-M., Yoon, U., Shin, Y.-W., Hong, S.B., Kim, I.Y., Kwon, J.S., Kim, S.I., 2006.
Fractal dimension in human cortical surface: multiple regression analysis withcortical thickness, sulcal depth, and folding area. Hum. Brain Mapp. 27, 994–1003.
Jiang, J., Zhu, W., Shi, F., Zhang, Y., Lin, L., Jiang, T., 2008. A robust and accurate algorithm
for estimating the complexity of the cortical surface. J. Neurosci. Meth. 172,
122–130.
Kalmanti, E., Maris, T.G., 2007. Fractal dimension as an index of brain cortical changes
throughout life. In Vivo 21, 641–646.
Kawasaki, Y., Suzuki, M., Kherif, F., Takahashi, T., Zhou, S.-Y., Nakamura, K., Matsui, M.,
Sumiyoshi, T., Seto, H., Kurachi, M., 2007. Multivariate voxel-based morphometry
successfully differentiates schizophrenia patients from healthy controls. Neuroimage 34, 235–242.
King, R., George, A., Jeon, T., Hynan, L., Youn, T., Kennedy, D., Dickerson, B., 2009.
Characterization of atrophic changes in the cerebral cortex using fractal
dimensional analysis. Brain Imaging Behav. 3, 154–166.
King, R.D., Brown, B., Hwang, M., Jeon, T., George, A.T., 2010. Fractal dimension analysisof the cortical ribbon in mild Alzheimer's disease. Neuroimage 53, 471–479.
Kiselev, V.G., Hahn, K.R., Auer, D.P., 2003. Is the brain cortex a fractal? Neuroimage 20,
1765–1774.
Kostelec, P.J., Maslen, D.K., Dennis, M., Healy, J., Rockmore, D.N., 2000. Computational
harmonic analysis for tensor ﬁelds on the two-sphere. J. Comput. Phys. 162,
514–535.
Lee, J.-M., Yoon, U., Kim, J.-J., Kim, I.Y., Lee, D.S., Kwon, J.S., Kim, S.I., 2004. Analysis of thehemispheric asymmetry using fractal dimension of a skeletonized cerebral surface.
IEEE Trans. Biomed. Eng. 51, 1494–1498.
Liu, J.Z., Zhang, L.D., Yue, G.H., 2003. Fractal dimension in human cerebellum measured
by magnetic resonance imaging. Biophys. J. 85, 4041–4046.
Lopes, R., Betrouni, N., 2009. Fractal and multifractal analysis: a review. Med. Image
Anal. 13, 634–649.
Luders, E., Narr, K.L., Thompson, P.M., Rex, D.E., Jancke, L., Steinmetz, H., Toga, A.W.,
2004. Gender differences in cortical complexity. Nat. Neurosci. 7, 799–800.
Luders, E., Thompson, P.M., Narr, K.L., Toga, A.W., Jancke, L., Gaser, C., 2006. A curvaturebased approach to estimate local gyriﬁcation on the cortical surface. Neuroimage
29, 1224–1230.
Lui, L.M., Wong, T.W., Zeng, W., Gu, X., Thompson, P.M., Chan, T.F., Yau, S.T., 2010.
Detection of shape deformities using Yamabe ﬂow and Beltrami coefﬁcients.
Inverse Probl. Imaging 4, 311–333.
MacDonald, D., Kabani, N., Avis, D., Evans, A.C., 2000. Automated 3-D extraction of innerand outer surfaces of cerebral cortex from MRI. Neuroimage 12, 340–356.

Mandelbrot, B.B., Blumen, A., 1989. Fractal geometry: what is it, and what does it do?
[and Discussion]. Proceedings of the Royal Society of London. Series A. Math. Phys.
Sci. 423, 3–16.
Mangin, J.-F., Frouin, V., Bloch, I., Régis, J., López-Krahe, J., 1995. From 3D magnetic
resonance images to structural representations of the cortex topography using
topology preserving deformations. J. Math. Imaging Vis. 5, 297–318.
McDonald, C., Marshall, N., Sham, P.C., Bullmore, E.T., Schulze, K., Chapple, B., Bramon,E., Filbey, F., Quraishi, S., Walshe, M., Murray, R.M., 2006. Regional brain
morphometry in patients with schizophrenia or bipolar disorder and theirunaffected relatives. Am. J. Psychiatry 163, 478–487.
McIntosh, A.M., Moorhead, T.W.J., McKirdy, J., Hall, J., Sussmann, J.E.D., Stanﬁeld, A.C.,
Harris, J.M., Johnstone, E.C., Lawrie, S.M., 2009. Prefrontal gyral folding and its
cognitive correlates in bipolar disorder and schizophrenia. Acta Psychiatr. Scand.
119, 192–198.
Mehaute, A.l., 1991. Fractal geometries: theory and applications. CRC Press, Boca Raton,FL.
Morra, J.H., Tu, Z., Apostolova, L.G., Green, A.E., Avedissian, C., Madsen, S.K., Parikshak,N., Hua, X., Toga, A.W., Jack, C.R., Schuff, N., Weiner, M.W., Thompson, P.M., 2009.
Automated 3D mapping of hippocampal atrophy and its clinical correlates in 400
subjects with Alzheimer's disease, mild cognitive impairment, and elderly controls.
Hum. Brain Mapp. 30, 2766–2788.
Narr, K.L., Bilder, R.M., Kim, S., Thompson, P.M., Szeszko, P., Robinson, D., Luders, E.,
Toga, A.W., 2004. Abnormal gyral complexity in ﬁrst-episode schizophrenia. Biol.
Psychiatry 55, 859–867.
Narr, K.L., Thompson, P.M., Sharma, T., Moussai, J., Zoumalan, C., Rayman, J., Toga, A.W.,
2001. Three-dimensional mapping of gyral shape and cortical surface asymmetries
in schizophrenia: gender effects. Am. J. Psychiatry 158, 244–255.
Nenadic, I., Sauer, H., Gaser, C., 2010. Distinct pattern of brain structural deﬁcits insubsyndromes of schizophrenia delineated by psychopathology. Neuroimage 49,
1153–1160.
Noga, J.T., Aylward, E., Barta, P.E., Pearlson, G.D., 1995. Cingulate gyrus in schizophrenic
patients and normal volunteers. Psychiatry Res. Neuroimaging 61, 201–208.
Oldﬁeld, R.C., 1971. The assessment and analysis of handedness: the Edinburgh
inventory. Neuropsychologia 9, 97–113.
Pantelis, C., Velakoulis, D., McGorry, P., Wood, S., Suckling, J., Phillips, L., Yung, A.,
Bullmore, E., Brewer, W., Soulsby, B., Desmond, P., McGuire, P., 2003. Neuroanatomical abnormalities before and after onset of psychosis: a cross-sectional andlongitudinal MRI comparison. Lancet 361, 281–288.
Sallet, P.C., Elkis, H., Alves, T.M., Oliveira, J.R., Sassi, E., de Castro, C.C., Busatto, G.F.,
Gattaz, W.F., 2003. Reduced cortical folding in schizophrenia: an MRI morphometric study. Am. J. Psychiatry 160, 1606–1613.
Sandu, A.-L., Rasmussen, Lundervold, A., Kreuder, F., Neckelmann, G., Hugdahl, K.,
Specht, K., 2008. Fractal dimension analysis of MR images reveals grey matterstructure irregularities in schizophrenia. Comput. Med. Imaging Graph. 32,
150–158.
Schaer, M., Cuadra, M.B., Tamarit, L., Lazeyras, F., Eliez, S., Thiran, J.P., 2008. A surfacebased approach to quantify local cortical gyriﬁcation. IEEE Trans. Med. Imaging 27,
161–170.
Shyu, K.K., Wu, Y.T., Chen, T.R., Chen, H.Y., Hu, H.H., Guo, W.Y., 2010. Measuring
complexity of fetal cortical surface From MR images using 3-D modiﬁed boxcounting method. IEEE Trans. Instrum. Meas. 1–10.
Soriano-Mas, C., Pujol, J.S., Alonso, P., Cardoner, N., MenchÛn, J.M., Harrison, B.J., Deus, J.,
Vallejo, J., Gaser, C., 2007. Identifying patients with obsessive-compulsive disorderusing whole-brain anatomy. Neuroimage 35, 1028–1037.
Steen, R.G., Mull, C., McClure, R., Hamer, R.M., Lieberman, J.A., 2006. Brain volume inﬁrst-episode schizophrenia: systematic review and meta-analysis of magnetic
resonance imaging studies. Br. J. Psychiatry 188, 510–518.
Sun, D., van Erp, T.G.M., Thompson, P.M., Bearden, C.E., Daley, M., Kushan, L., Hardt, M.E.,
Nuechterlein, K.H., Toga, A.W., Cannon, T.D., 2009. Elucidating a magnetic
resonance imaging-based neuroanatomic biomarker for psychosis: classiﬁcation
analysis using probabilistic brain atlas and machine learning algorithms. Biol.
Psychiatry 66, 1055–1060.
Thompson, P.M., Bartzokis, G., Hayashi, K.M., Klunder, A.D., Lu, P.H., Edwards, N., Hong,
M.S., Yu, M., Geaga, J.A., Toga, A.W., Charles, C., Perkins, D.O., McEvoy, J., Hamer, R.
M., Tohen, M., Tollefson, G.D., Lieberman, J.A., Group, T.H.S., 2009. Time-lapse
mapping of cortical changes in schizophrenia with different treatments. Cereb.
Cortex 19, 1107–1123.
Thompson, P.M., Hayashi, K.M., Sowell, E.R., Gogtay, N., Giedd, J.N., Rapoport, J.L., de
Zubicaray, G.I., Janke, A.L., Rose, S.E., Semple, J., Doddrell, D.M., Wang, Y., van Erp, T.
G.M., Cannon, T.D., Toga, A.W., 2004. Mapping cortical change in Alzheimer's
disease, brain development, and schizophrenia. Neuroimage 23, S2–S18.
Thompson, P.M., Lee, A.D., Dutton, R.A., Geaga, J.A., Hayashi, K.M., Eckert, M.A., Bellugi,U., Galaburda, A.M., Korenberg, J.R., Mills, D.L., Toga, A.W., Reiss, A.L., 2005.
Abnormal cortical complexity and thickness proﬁles mapped in Williams
syndrome. J. Neurosci. 25, 4146–4158.
Thompson, P.M., Schwartz, C., Lin, R.T., Khan, A.A., Toga, A.W., 1996. Three-dimensional
statistical analysis of sulcal variability in the human brain. J. Neurosci. 16,
4261–4274.
Toro, R., Perron, M., Pike, B., Richer, L., Veillette, S., Pausova, Z., Paus, T., 2008. Brain sizeand folding of the human cerebral cortex. Cereb. Cortex 18, 2352–2357.
Vogeley, K., Schneider-Axmann, T., Pfeiffer, U., Tepest, R., Bayer, T.A., Bogerts, B., Honer,
W.G., Falkai, P., 2000. Disturbed gyriﬁcation of the prefrontal region in maleschizophrenic patients: a morphometric postmortem study. Am. J. Psychiatry 157,
34–39.
Vogeley, K., Tepest, R., Pfeiffer, U., Schneider-Axmann, T., Maier, W., Honer, W.G., Falkai,P., 2001. Right frontal hypergyria differentiation in affected and unaffected siblings

R.A. Yotter et al. / NeuroImage 56 (2011) 961–973
from families multiply affected with schizophrenia: a morphometric MRI study.
Am. J. Psychiatry 158, 494–496.
Wang, Y., Zhang, J., Gutman, B., Chan, T.F., Becker, J.T., Aizenstein, H.J., Lopez, O.L.,
Tamburo, R.J., Toga, A.W., Thompson, P.M., 2010. Multivariate tensor-based
morphometry on surfaces: application to mapping ventricular abnormalities in
HIV/AIDS. Neuroimage 49, 2141–2157.
Weinberger, D., McClure, R., 2002. Neurotoxicity, neuroplasticity, and magnetic
resonance imaging morphometry: what is happening in the schizophrenic brain?
Arch. Gen. Psychiatry 59, 553–558.
Wheeler, D.G., Harper, C.G., 2007. Localised reductions in gyriﬁcation in the posterior
cingulate: schizophrenia and controls. Prog. Neuro Psychopharmacol. Biol.
Psychiatry 31, 319–327.
Wright, I.C., Rabe-Hesketh, S., Woodruff, P.W.R., David, A.S., Murray, R.M., Bullmore, E.
T., 2000. Meta-analysis of regional brain volumes in schizophrenia. Am. J.
Psychiatry 157, 16–25.
Wu, Y.-T., Shyu, K.-K., Chen, T.-R., Guo, W.-Y., 2009. Using three-dimensional fractal
dimension to analyze the complexity of fetal cortical surface from magnetic
resonance images. Nonlinear Dyn. 58, 745–752.
Wu, Y.-T., Shyu, K.-K., Jao, C.-W., Wang, Z.-Y., Soong, B.-W., Wu, H.-M., Wang, P.-S.,
2010. Fractal dimension analysis for quantifying cerebellar morphological changeof multiple system atrophy of the cerebellar type (MSA-C). Neuroimage 49,
539–551.

973

Xu, C., Pham, D.L., Rettmann, M.E., Yu, D.N., Prince, J.L., 1999. Reconstruction of thehuman cerebral cortex from magnetic resonance images. Med. Imaging IEEE Trans.
18, 467–480.
Yotter, R.A., Thompson, P.M., Gaser, C., 2010. Estimating local surface complexity maps
using spherical harmonic reconstructions. In: Jiang, T. (Ed.), MICCAI. Springer
Berlin, Beijing, China.
Yushkevich, P., Dubb, A., Xie, Z., Gur, R., Gur, R., Gee, J., 2005. Regional structuralcharacterization of the brain of schizophrenia patients. Acad. Radiol. 12, 1250–1261.
Zhang, L., Butler, A.J., Sun, C.-K., Sahgal, V., Wittenberg, G.F., Yue, G.H., 2008a. Fractal
dimension assessment of brain white matter structural complexity post stroke inrelation to upper-extremity motor function. Brain Res. 1228, 229–240.
Zhang, L., Dean, D., Liu, J.Z., Sahgal, V., Wang, X., Yue, G.H., 2007. Quantifying
degeneration of white matter in normal aging using fractal dimension. Neurobiol.
Aging 28, 1543–1555.
Zhang, L., Liu, J.Z., Dean, D., Sahgal, V., Yue, G.H., 2006. A three-dimensional fractal
analysis method for quantifying white matter structure in human brain. J. Neurosci.
Meth. 150, 242–253.
Zhang, Y., Jiang, J., Lin, L., Shi, F., Zhou, Y., Yu, C., Li, K., Jiang, T., 2008b. A surface-based
fractal information dimension method for cortical complexity analysis. Med.
Imaging Augmented Real. 133–141.
Zilles, K., Armstrong, E., Schleicher, A., Kretschmann, H.-J., 1988. The human pattern of
gyriﬁcation in the cerebral cortex. Anat. Embryol. 179, 173–179.