r Human Brain Mapping 00:000–000 (2010) r

Topological Correction of Brain Surface Meshes
Using Spherical Harmonics
Rachel Aine Yotter,1* Robert Dahnke,1 Paul M. Thompson,2and Christian Gaser1
1

2

Department of Psychiatry, Friedrich-Schiller University, Jena, Germany
Laboratory of Neuro Imaging, Division of Brain Mapping, Department of Neurology,
UCLA School of Medicine, Los Angeles, Californiarr

Abstract: Surface reconstruction methods allow advanced analysis of structural and functional brain
data beyond what can be achieved using volumetric images alone. Automated generation of corticalsurface meshes from 3D brain MRI often leads to topological defects and geometrical artifacts that must
be corrected to permit subsequent analysis. Here, we propose a novel method to repair topologicaldefects using a surface reconstruction that relies on spherical harmonics. First, during reparameterization of the surface using a tiled platonic solid, the original MRI intensity values are used as a basis toselect either a ‘‘fill’’ or ‘‘cut’’ operation for each topological defect. We modify the spherical map of theuncorrected brain surface mesh, such that certain triangles are favored while searching for the bounding
triangle during reparameterization. Then, a low-pass filtered alternative reconstruction based on spherical harmonics is patched into the reconstructed surface in areas that previously contained defects. Selfintersections are repaired using a local smoothing algorithm that limits the number of affected points toless than 0.1% of the total, and as a last step, all modified points are adjusted based on the T1 intensity.
We found that the corrected reconstructions have reduced distance error metrics compared with a‘‘gold standard’’ surface created by averaging 12 scans of the same brain. Ninety-three percent of the
topological defects in a set of 10 scans of control subjects were accurately corrected. The entire process
takes 6–8 min of computation time. Further improvements are discussed, especially regarding the useof the T1-weighted image to make corrections. Hum Brain Mapp 00:000–000, 2010. VC 2010 Wiley-Liss, Inc.
Key words: topology correction; spherical harmonics; surface reconstruction; topological defects; noise;
self-intersections; MRIrr

INTRODUCTION

Contract grant sponsor: German Bundesministerium für Bildung
und Forschung; Contract grant numbers: BMBF 01EV0709, BMBF
01GW0740; Contract grant sponsor: NIH; Contract grant numbers:
EB008432, EB008281, EB007813, HD050735.
*Correspondence to: Rachel Aine Yotter, Department of Psychiatry, Friedrich-Schiller University, Jahnstrasse 3, 07745 Jena,
Germany. E-mail: rachel.yotter@uni-jena.deReceived for publication 28 December 2009; Revised 23 March
2010; Accepted 19 April 2010
DOI: 10.1002/hbm.21095
Published online in Wiley InterScience (www.interscience.wiley.
com).
C 2010 Wiley-Liss, Inc.
V

In many analyses of structural and functional brain
mapping data, it is often desirable to deform the gray matter (GM) sheet into an inflated brain surface or a sphere.
For example, the functional subdivisions of the cerebral
cortex are the target of most functional brain mapping
studies, but they lie within a convoluted sheet that is sohighly folded, at least in humans, that it is challenging toalign the cortical anatomy across subjects without performing further modeling of the cortex. As such, 2 decades of
surface-based cortical mapping research have been
devoted to finding convenient methods to generate corticalsurface reconstructions as well as cortical feature identification, alignment, and surface-based statistical analyses forr Yotter et al. r

Figure 1.
Common problems include topological defects and artifacts. Topological defects such as handles (a)and holes (b) prevent the surface from being homeomorphic with a sphere. By contrast, artifacts (c,highlighted in red) may permit a correct (spherical, genus-zero) topology but nonetheless are anatomically incorrect and should ideally be eliminated during topology correction. [Color figure can beviewed in the online issue, which is available at www.interscience.wiley.com.]comparing cortical signals across subjects and groups[Fischl et al., 1999b; MacDonald et al., 2000; Shattuck andLeahy, 2001; Thompson and Toga, 1996; Thompson et al.,
2004; Van Essen, 2004, 2005].
As well as improving the power and localization of functional imaging analyses [Formisano et al., 2004], surfacebased analyses have also discovered patterns of cortical
thickness alterations in numerous diseases, such as Alzheimer’s disease [Lerch et al., 2005; Thompson et al., 2004],
schizophrenia [Bearden et al., 2007a; Narr et al., 2005;
Thompson et al., 2009], epilepsy [Lin et al., 2007], major
depression [Ballmaier et al., 2004], and genetic syndromes
such as Williams syndrome and 22q deletion syndrome[Bearden et al., 2007b; Thompson et al., 2005]. Reconstructedcortical surface models can also provide the basis to understand subtle and distributed morphometric features that
may be altered in disease such as gyrification [Gaser et al.,2006], surface curvature [Tosun et al., 2006], surface metrictensors [Wang et al., 2010], and other surface abnormalities.
When analyzing cortical surface data from multiple subjects, a common coordinate system must be imposed to
extract meaningful comparisons between the subjects. Several coordinate systems have been proposed and used,
including flattened representations [Qiu et al., 2007; Van
Essen and Drury, 1997], hyperbolic spaces generated using
circle packing methods [Hurdal and Stephenson, 2009], slit
maps [Wang et al., 2008], and punctured disks [Wanget al., 2009], but the most commonly used is a sphericalcoordinate system, as the cortical surface is roughly homeomorphic to a sphere. Besides intersubject analysis, parametric models of the cortical surfaces can be used toenable a more precise intersubject registration by enforcing
a higher order matching of sulcal features or geometriclandmarks lying in the cortex [Desai et al., 2005; Hinds
et al., 2008; Lui et al., 2006, 2007; Thompson et al., 2004;
Wang et al., 2005].

r

However, the surface to be analyzed should be genus zero(e.g., contain no topological defects) before being inflated or
mapped to a sphere. It is fairly straightforward to cut andflatten a surface that has a spherical topology, but this task
becomes nearly impossible when there are topologicaldefects, making it necessary to use some approach for topological correction. Some methods extract the cortex by gradually deforming a spherical surface into the configuration of
the cortex, so a spherical topology is maintained and guaranteed [MacDonald et al., 2000]. However, many methods identify the cortex using a voxel-level segmentation by applying
a classification function to each voxel independently, andbecause of segmentation errors, there is no guarantee that theinitially created surfaces from the tissue segmentation maps
are genus zero (i.e., homotopic to a sphere). Correcting these
topological defects is a necessary prerequisite for intersubject
analysis. As such, several topology correction approaches
have been advocated using local ‘‘manifold surgery’’ or Reeb
graph methods [Abrams et al., 2002, 2004; Fischl et al., 2001].
There are three possible representations of the corticalsurface: the interface between GM and white matter (WM);
the interface between GM and cerebrospinal fluid(CSF), also called the pial surface; and the central surface(CS), which is approximately the midpoint between the
GM/WM and GM/CSF interfaces. Compared with the
WM/GM or GM/CSF interfaces, the CS provides an inherently less distorted representation of the cortical surface[Van Essen and Maunsell, 1980]. However, the intensity
differential across the WM/GM interface in the T1 image
is a source of useful information for labeling topologicaldefects for either fill or cut operations and improving themesh quality after topological correction.
Because of noise, partial volume effects, and other artifacts
during the MRI data acquisition process, a brain surface
mesh reconstructed from volumetric data often contains
topological defects and geometric artifacts. Figure 1 shows

2 r

r Cortical Surface Topology Correction r

some examples of topological defects and geometric artifacts for a surface generated using the standard FreeSurfer
pipeline, i.e., the intensity gradient in the T1 image is used
to classify tissue and a marching cubes algorithm generates the surface at the boundary between the gray andwhite tissue segmentation maps [Fischl et al., 2001].
Topological defects can include handles (Fig. 1a) andholes (Fig. 1b) that prevent the surface from being homeomorphic with a sphere. In terms of topology, handles andholes are equivalent; however, the two may be differentiated in that during the topology correction process, handles
should be cut but the holes should be filled. The number ofdefects in a given surface can be calculated using the Euler
characteristic, where an Euler characteristic of 2 represents
a surface that is homeomorphic with a sphere.1
Artifacts (Fig. 1c) are topologically correct sharp peaks(usually due to noise) that do not reflect true features ofcortical anatomy. All of these errors need to be repaired
before the surface can be inflated to a sphere. Because artifacts usually contain sharper edges than the rest of the cortical surface, they can be minimized via a smoothing
process. However, accurately correcting topology defects ismore difficult. Compared with previously proposed methods, we propose the use of regional spherical harmonic representations as a simple, straightforward method to correctthese topological errors and geometrical artifacts.

Prior Approaches to Correct Topological Defects
There are two general approaches used to address topology defects. The first approach is to start with a surface
with the desired topology and deform it to match thebrain surface (a top-down approach). The surface representation may either be discrete (such as a mesh) or continuous, and there are a number of approaches todeforming this initial surface. In this category are included
active contours or deformable models. There is a wealth ofliterature on variations of this method applied to brain
data; for a review, see [McInerney and Terzopoulos, 1996;
Montagnat et al., 2001]. Methods include fast marchingalgorithms [Bazin and Pham, 2005a,b, 2007a,b], level-set
methods [Bischoff and Kobbelt, 2002, 2003, 2004; Xiao
et al., 2003], minimization of an energy functional [Caselles
et al., 1997; Dale et al., 1999; Karaçali and Davatzikos,
2003; Kass et al., 1988; MacDonald et al., 2000], competitive
forces [Delingette and Montagnat, 2001; McInerney andTerzopoulos, 2000], a system of constraints [Lachaud andMontanvert, 1999], simplex meshes [Delingette, 1999],
homotopically deformable regions [Mangin et al., 1995],
and diffeomorphic transformations of a digitized anatomical template [Christensen et al., 1994]. As the cortical surface is highly convoluted, a disadvantage of theseapproaches is that the fitting function is highly nonlinear
and may require extensive computations. Furthermore, theresulting surface may not represent all of the deep sulci
accurately [Fischl et al., 2001]. This is because typicallyhigh curvature regions are penalized during the evolutionof the template, to make the extractions more robust tohigh frequency noise, and the resulting enforcement ofregularity may restrain the evolving surface from correctly
flowing into the deep sulci.
The second approach is to retrospectively correct the
topology after the brain data have been segmented into WM
and GM (a bottom-up approach). These corrections may beachieved either in volume space or on a mesh surface generated using an algorithm such as marching cubes [Lorensen
and Cline, 1987]. A common volume-based solution is to use
graph-based methods to detect topology defects and find theminimal modification required to correct the defect [Abrams
et al., 2002, 2004; Chen and Wagenknecht, 2006; Han et al.,2002, 2004; Jaume et al., 2005; Shattuck and Leahy, 2001].
It is also possible to find the minimal modification
directly on the surface [Guskov and Wood, 2001]. By
selecting just the minimal set, it is possible that the resulting correction is not identical to the changes that would bemade by an experienced operator. Another possible solution is to filter the WM membership until the largest triangular mesh has a geometry homeomorphic to a sphere [Xuet al., 1999]. A drawback of this approach is that the filtering process essentially smoothes the WM and some anatomical information can be lost. In volume space, it is also
possible to construct a skeleton to describe the shape of thevolume, correct cycles within the skeleton, and then regrow
the skeleton to form a topologically correct object [Zhouet al., 2007]. Finally, another volume-based solution is basedon region-growing algorithms to correct the initial segmentation [Kriegeskorte and Goebel, 2001; Ségonne et al., 2003].
Surface-based methods often do not correct large handles appropriately but are appropriate for small topological defects. One method to circumvent this problem is tomeasure the defect size and then correct small defects inthe surface and large defects in volume space [Woodet al., 2004]. The original MRI data may also inform retrospective topological corrections on the surface. The correction can then be carried out using purely surface-based
methods such as genetic algorithms [Ségonne et al., 2005],
retessellation of the defect [Fischl et al., 2001], or nonseparating loops with opening operators [Segonne et al., 2007].
These methods vary widely in speed and accuracy. For
deformable surfaces, the speed is usually related to the complexity of the fitting function. For retrospective correction
methods, pure surface-based methods tend to be slower,
while it is possible to quickly repair topological defects using
voxel-based methods but possibly at a cost to accuracy.

1

The Euler characteristic is defined for discrete surfaces (such as
polyhedra) according to the formula v 2 e 1 f, where v, e, and f are
the number of vertices (corners), edges, and faces, respectively, in agiven polyhedron.

r

Spherical Harmonics
We propose to use spherical harmonics for the first timeto accurately correct topological defects directly on the

3 r

r Yotter et al. r

Figure 2.
The full processing pipeline is illustrated using a simple shape the search for nearest intersecting triangle. These regularly
with a handle. The original surface mesh (a) is mapped onto a sampled points are forward- then reverse-transformed using
sphere (b). The radii of points in the handle are modified so spherical harmonics. High-frequency (e) and low-pass filtered (f)that they are no longer on the surface of the sphere (c). A regu- reconstructed surfaces are derived from the spherical harmonic
larly sampled grid (d, green) is then overlaid on top of this coefficients. In the regions containing the defect, vertices fromsphere to create uniform sampling in the parameter space. For the low-pass filtered reconstruction are patched into the higheach point in the regularly sampled grid (d, green), the intersect- bandwidth reconstruction (g). Finally, a T1 postcorrection step
ing triangle in the spherical mapping (d, blue) is found, and bary- removes any remaining artifacts from the removal of the handle
centric coordinates within this triangle are used to interpolate a (h). [Color figure can be viewed in the online issue, which is
spatial coordinate for this vertex lying on the original tessellated available at www.interscience.wiley.com.]
surface mesh. Points with modified radii are not favored during
surface. Spherical harmonics have been recently applied tobrain surface meshes, usually in the realm of shape analysis. First applied to brain imaging in 1995, as a method todetect shape differences between closed objects [Brechbühler et al., 1995], and for representing the cortex and flows
between surfaces [Thompson and Toga, 1996], sphericalharmonic analysis has been subsequently applied to quantify structural differences in subcortical structures [Geriget al., 2001; Gutman et al., 2006; Kelemen et al., 1999; Shenton et al., 2002] and full cortical surfaces [Chung et al.,2006; Shen and Chung, 2006].
Spherical harmonic analysis can be described as a Fourier transform on a sphere that decomposes the brain surface data into frequency components; the sphericalharmonics are the eigenfunctions of the Laplacian operator
on the sphere, and they are often used a basis to computeand estimate solutions to diffusion equations or flows in
spherical coordinates, as they provide a complete orthogonal sequence of functions that can approximate any(square-integrable) smooth function defined on the sphere(including, by analogy, the 3D spatial coordinates of aspherically parameterized surfaces, as they can be considered as functions on the sphere). The surface can also bereconstructed from the coefficients of the spherical har-r

monics. If spherical harmonics are used to approximate areconstructed surface, the RMS error of the reconstructedsurface decreases as the number of coefficients increases
and can become arbitrary small (i.e., uniform approximation). Furthermore, the coefficients of the spherical harmonics may be weighted to achieve a smoothing effect[Chung et al., 2007].
A drawback to spherical harmonics is the computation
time required to calculate the coefficients. However, amodification of the fast Fourier transform compatible withspherical coordinates dramatically decreases the computation time, as it is no longer necessary to directly calculate
the coefficients [Healy et al., 1996; Kostelec et al., 2000].
In this article, we present a method to accurately correct
topological defects and artifacts using spherical harmonics.
Spherical harmonic functions retrospectively correct topological errors directly on the brain surface mesh. In brief,
we resample a spherical mapping of the cortical surface
mesh and calculate 3D spatial coordinates for each spherically sampled point based on its approximate location in
the cortical surface mesh (Fig. 2b–d). The coordinates oflocation mapped on the regularly sampled sphere are forward transformed and then inversely transformed via
spherical harmonics to reconstruct a cortical surface

4 r

r Cortical Surface Topology Correction r

without topological defects (Fig. 2e,f). In a high-frequency
reconstruction, the topological defects are replaced with aspiked surface (high sharpness), which is repaired byreplacing local patches with points from a low-pass filtered reconstructed surface (Fig. 2g). The result is a moreaccurately reconstructed cortical surface free from topological defects and large geometrical artifacts (Fig. 2h).

METHODS
Sample Data Sets
Three data sets were used to evaluate the validity of thisapproach: (a) an artificial structure containing topologicaldefects and geometrical artifacts; (b) 12 MRI scans of thesame brain, including the creation of a ‘‘gold standard’’ byaveraging all 12 scans and manually correcting any
remaining topology defects; and (c) scans of 10 healthy
control subjects.
The topological phantom is essentially a cube that has
three holes, two handles, and two ‘‘artifacts.’’ The 60-mm3
phantom was made by simulating a 1-mm3 T1 intensity
image using routines available from the Statistical Parametric Mapping software package (http://www.fil.ion.ucl.
ac.uk/spm/software/) and in-house MatLab code andthen constructing a surface mesh with !50 K triangles
using a marching cubes algorithm. This surface mesh isthen corrected using our method. Another simulated ‘‘gold
standard’’ T1-weighted image was also generated as an
input for T1 image postcorrection. By comparing to the‘‘gold standard’’ mesh surface created directly from the‘‘gold standard’’ T1 intensity image, it is possible to quantitatively assess the validity of the spherical harmonicreconstruction.
The second approach was to create a ‘‘gold standard’’
brain by averaging 12 scans of the same subject. The averaging reduces Gaussian noise and the size and number of
the topological defects but cannot remove all topologicaldefects entirely. The remaining defects were extremely
small and corrected manually using FreeSurfer [Fischl
et al., 1999a]. The sample data set included 12 brain scansof the same subject performed on two different 1.5 T Siemens Vision scanners. Both scanners used 3D magnetization-prepared gradient echo (MP-RAGE) T1-weighted
sequences of 160 sagittal slices with voxel dimensions 1
mm " 1 mm " 1 mm and FOV ¼ 256 mm. Scanner 1
parameters were TR/TE/FA ¼ 11.4 ms/4.4 ms/15$ , andScanner 2 parameters were TR/TE/FA ¼ 15 ms/5 ms/30$ .
To verify that the spherical harmonic-based correction
approach was not optimized for a single brain, a third
dataset of 10 healthy control subjects’ brain scans was also
included. For the 10 control subjects (4 females/6 males,mean age 35.3 years, %11.3 SD), T1-weighted images wereobtained on a 1.5 T Philips Gyroscan ACSII. There were
256 sagittal slices per scan (1-mm thickness, TR ¼ 13 ms,
TE ¼ 5 ms, flip angle ¼ 25$ , field-of-view ¼ 256 mm) witha matrix size of 256 " 256, resulting in an isotropic voxelrsize of 1 mm3. As no ‘‘gold standard’’ is available for comparison, the correction for these brains was validated by
an expert who scored the validity of topological defect corrections in brains corrected either using our method or by
FreeSurfer. The brains were randomly labeled such thatthe scorer was not informed about which correction
method had been used to correct each brain.
For both datasets (12 scans of the same brain and 10
control subjects), each brain scan was processed using thedefault FreeSurfer processing pipeline to produce a WM
surface representation for each hemisphere [Fischl et al.,1999a]. The surface used as input into our topological correction method was the xh.smoothwm surface before topological correction. The spherical mapping was additionally
postprocessed using in-house tools to improve the accuracy of the reparameterization [Yotter et al., in-press]. All
individual surface meshes included a unique set of topological defects.

Software Implementation
The topology correction program was implemented in
the C programming language. The code related to spherical harmonics was based on a publicly available sphericalharmonic processing library [Kostelec and Rockmore,2004], and the final correction using the T1 intensity
images used functions extracted from the BICPL software
library (McConnell Brain Imaging Centre Programming
Library, http://packages.bic.mni.mcgill.ca/). These functions are described in more detail below. All other code
was developed in-house.

Cut/Fill Operations and the Initial
Spherical Mapping
A difficult problem to solve in brain surface mesh reconstruction is how to fill holes and cut handles. We propose
a simple and fast new method compatible with the spherical harmonics approach to accomplish this.
To apply spherical harmonic analysis to the brain surface mesh, a spherical mapping of the uncorrected surface
mesh must first be generated. This spherical mapping is
not homeomorphic with a sphere, because it still contains
defects. It is also beneficial to reduce area distortion in thespherical mapping, because this improves the reparameterization process, during which a coordinate of location
based on the estimated cortical surface location is generated for each regularly spaced point in spherical coordinates [Yotter et al., in-press]. Additionally, this sphericalmapping may also be used to identify topological defects
by locating self-intersections in the spherical mapping[Fischl et al., 1999a]. Points associated with topologicaldefects will be referred to as defect points in the following
discussion.
During the reparameterization process (detailed below),
regularly sampled spherical points are associated with a

5 r

r Yotter et al. r

Figure 3.
Defects are labeled as either holes for filling or handles for cut- mapping to have a slightly smaller radius (c). When searching
ting before the surface is reparameterized. In this hole, the point for the closest triangle during the reparameterization process,set is first bisected (a). The original spherical mapping with the outer triangles will be strongly favored. [Color figure canthe bisected points marked is shown in (b). Inner points (pre- be viewed in the online issue, which is available atviously marked as green) are then modified on the spherical www.interscience.wiley.com.]

cortical location by finding the closest triangle in thespherical mapping and then using barycentric coordinates
to interpolate its location on that triangular tile. By varying
the radius of the defect points within the defects, it is possible to favor a subset of defect points, such that the outer
points for holes (or the inner points for handles) are more
likely to be chosen. These points are favored because
they are closer to the surface of the reparameterized
sphere, whereas the other points with nonunity radius
are further away than they would have been had theirradii not been modified. This affects the search for theclosest point when mapping the original set of points to
the reparameterized mesh.
First, the defects must be labeled for either cut or fill
operations. This is accomplished by examining the intensity of the T1 image in the space central to the defect. If it
is above a T1 intensity threshold (set to 99% of the average
T1 threshold of nondefect points, a threshold which empirically performed best for correctly labeling holes and handles), then it is labeled as a hole to be filled; otherwise,
it is a handle to be cut. For large defect areas with many
topological errors (e.g., when the Euler characteristic isless than &10), this approach is not applicable. These large
defect clusters invariably occur either at the ventricles
or in the nonbrain area around the orbitofrontal cortex (thecorrection of which is shown in Fig. 11). For the ventricular
defects, they are identified using their 3D position within thebrain and marked as a hole to be filled; for the defects
around the orbitofrontal cortex, they are always cut.
Once a defect has been labeled as a target for either a‘‘cut’’ or a ‘‘fill’’ operation, the defect points are thenbisected into two groups of points, based on how closerthey are to the cortical centroid, along the axis of the cortical defect.2 For defects marked for filling, the defect points
closer to the cortical centroid are adjusted in the sphericalmapping such that their radii are slightly smaller; for cutting,defect points farther from the cortical center have adjusted
radii that are larger. Thus, during the reparameterization process, the closest polygon to the reparameterized point willmore likely be the outer (for filling) or inner (for cutting)points within the defect region that contains overlapping triangles in the spherical mapping, as nonfavored points havebeen moved away from the surface of the reparameterized
sphere. Figure 3 demonstrates this method for a hole.

Spherical Parameterization and Generationof Coordinates of Location
To analyze the harmonic content of a spherical surface,
the spherical surface needs to have regularly sampled
points with respect to y and f, where y is the colatitudeand f is the azimuthal coordinate. As the original corticalsurface does not generally map to a regularly sampled
sphere, the first step in using spherical harmonics is to
reparameterize the surface.
To do this, points are generated from equally sampledvalues of y and f for all members in the sets, such that there
are 2B points per set, where B is the bandwidth. For eachregularly sampled spherical point, the closest polygon on
2

The defect axis is determined by summing the normals of defectpoints and then normalizing. This method takes advantage of thefact that the normals inside of the defect (as well as along the sides)
cancel out, while the outer points have no equivalently matchedpoints on the opposite side. Thus, the result is a directional vector
along the defect that can be used along with the cortical centroid tocorrectly bisect the defect points.

6 r

r Cortical Surface Topology Correction r

Figure 4.
Using a Butterworth filter to filter the coefficients results in a ringing artifacts can also be reduced by using a higher order spherismoother reconstructed surface. Both a boxcar filter (a) and cal harmonic approximation. The sharpness of a mesh point is thea Butterworth filter (b) leave some ringing artifacts, but using a maximum angle between the normals of nearest neighbor polyButterworth filter reduces the magnitude of the ringing artifacts. gons, in degrees. [Color figure can be viewed in the online issue,
The filtered reconstructions were obtained by initially calculating which is available at www.interscience.wiley.com.]coefficients for l ¼ 1,024, then filtering such that ll ¼ 64. Thesethe cortical spherical mapping is found. Because of thechange in the radius of defect points, these are slightly disfavored when searching for the closest polygon.
Within the closest polygon, a spatial location for theinterpolated vertex in 3D is approximated using barycentric coordinates. The result is a regularly sampled spherical map in which every point is associated with acoordinate that gives its location on the original corticalsurface (Fig. 2b–d).

A square-integrable function f(y,f) on the sphere can be
expanded in the spherical harmonic basis such that:f ðh;/Þ ¼

m
im/
Ym
;
l ðh; /Þ ¼ Pl ðcos hÞe(1)where l and m are integers with |m| ) l, and Pm
l is theassociated Legendre polynomial function defined by:Pml ðxÞ ¼"m lþm ! 2
"l
1 !
2 2 d
x &1 :
1
&
x
2l l!
dxlþm

(2)

r

(3)

l¼0 m¼&l

where the coefficients f̂(l,m) are defined by the inner prod2
uct f̂(l,m) hf,Yml i, and the L norm of each basis functionm
Yl is given by:

Harmonic Analysis
The spherical harmonic series representation of a spherical mesh can be obtained using normalized spherical harmonics Yml (y,f):

B X
l #
X
#
#Ym #&2^f ðl; mÞ " Ym ;l 2
l# m #&2#Y # ¼l

24p ðl þ mÞ!
+:
2l þ 1 ðl & mÞ!

(4)This system can be solved directly by finding the bases
first, but a more efficient approach is to use a divide-andconquer scheme as described in [Healy et al., 1996].
Using the calculated coefficients, two surfaces are reconstructed. The first surface is a high-frequency surface thatuses all coefficients as is (Fig. 2e). The second surface is a
low-bandwidth (smoothed) surface that is reconstructedfrom filtered coefficients using a 128-order Butterworth
low-pass filter, thus excluding contributions from higher
frequency coefficients (Fig. 2f). A Butterworth filter

7 r

r Yotter et al. r

Figure 5.
Some detail is lost as more coefficients are filtered. Filtering coefficients such that ll ¼ 32 or lower
result in a surface that no longer resembles the original surface. The numbers shown are the lower
value of l, and all surfaces were constructed such that l ¼ 1,024 and then filtered with a Butterworth
filter. An artifact (red arrow) is no longer visible in the surfaces reconstructed using spherical harmonics.
[Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]
reduces ringing artifacts and results in a smoother reconstructed surface (see Fig. 4). To reconstruct the surfaces,
the harmonic content is processed through an inverse Fourier transform to produce coordinates associated with eachspherical point. These coordinates are used to reconstruct
the cortical mesh without topological defects using an
icosahedron as the base mesh.

Patching of Reconstructed Cortical Meshes
After reconstruction via spherical harmonics, the resulting cortical mesh is homeomorphic to a sphere. However,two problems need to be addressed in the reconstructedsurface. First, adjacent to regions formerly containing topological defects, the spherical harmonic-based reconstructedsurface generally replaces the defect with a spiked topology that does not correspond to actual brain anatomy. The
spiked appearance is only present if the bandwidth is highenough to admit higher frequencies; if a lower bandwidthis used, the surface near former topological defects issmooth and well reconstructed, but anatomical detail islost. Second, the reconstructed surface often contains selfintersections that should be repaired; this problem isaddressed in the next section.
To correct the first problem, these surfaces are combinedsuch that the regions that exhibit spiked topology in thehigher frequency (l ¼ 1,024) surface are replaced with patchesfrom the low-pass filtered (ll ¼ 64) cortical reconstruction, but
only in regions previously marked as defect areas (Fig. 2g).
The lower limit was chosen such that the spiked regions weresmooth but the surface retained the approximate shape of theoriginal cortical surface, such that the union of the two surfaces does not result in large discontinuities (see Fig. 5).

r

Regions with spiked topology are identified by calculating a sharpness value for each point. The sharpness is the
maximum angle between the normals of nearest neighbor
polygons. If the sharpness is above a threshold of ts ¼ 60$ ,and if this point is labeled as being within a region containing a topological defect, then this point and its set of nearest
neighbors are replaced with the corresponding points in the
low-pass filtered cortical surface. This angular threshold
was empirically determined to make sure that all brain
surfaces were appropriately patched; however, as the patching only occurs within regions marked as areas formerly
containing topological defects, lowering the threshold further had little effect on the final reconstructed surface.

Correcting Self-Intersections
In all cases studied here, the reconstructed cortical surface contained self-intersections, i.e., the cortical surface
intersects itself. These self-intersections generally occur inregions formerly containing topological defects and arefound by testing for edge-triangle intersections. Self-intersections are undesirable, as it would mean that the same
imaging data would be represented in at least two nonadjacent locations on the unfolded cortical surface. This problem should be addressed as part of the topologicalcorrection process. There are two approaches capable ofresolving self-intersections: further patching using the lowpass filtered surface and localized smoothing.
As the low-pass filtered reconstruction is smoother thanthe high-frequency surface, it tends to contain fewer selfintersections, as the lower frequency functions cannot easily generate a self-intersection (as an extreme case, consider that low-pass filtering such that ll ¼ 1 results in a
8 r

r Cortical Surface Topology Correction r

Figure 6.
Self-intersections are corrected via localized smoothing. A region containing self-intersections(a) is smoothed until the cortical surface no longer intersects itself (b). The algorithm is implemented to affect as few points as possible. [Color figure can be viewed in the online issue, whichis available at www.interscience.wiley.com.]
sphere). By patching self-intersection regions with theequivalent points from the low-pass filtered surface, thenumber of self-intersections can be greatly reduced, but
usually not eliminated. This is the first step taken when
repairing self-intersections.
For remaining self-intersections, the mesh is smoothedusing a localized smoothing algorithm until the region nolonger contains self-intersections [Yotter et al., in-press). Figure 6 shows a typical example of self-intersection smoothing.
The localized smoothing algorithm was implemented such that
it would not affect more than 0.1% of the total numberof surface points. For details of this algorithm, see Appendix.

T1 Image-Based Postcorrection
It could be argued that the reparameterization, patching,and localized smoothing operations may result in a loss offidelity in the spherical mapping. Furthermore, there may
be some mislabeled defects that were cut/filled when the
opposite operation would have been more appropriate.
These issues can be counteracted by including a small
amount of surface postprocessing based on the intensitiesof the associated T1-weighted image (Fig. 2h). In our
implementation, all points that were modified via patchingor via self-intersection repair are later adjusted using theoriginal T1-weighted intensity image and previously
approximated T1 intensity threshold using the publicly
available software package BICPL (http://packages.bic.mni.mcgill.ca/). The approach uses a surface deformation
approach to correct the surface mesh based on the T1 intensity image [MacDonald et al., 1994]. This correction process uses only the T1 intensity image to reposition thesurface mesh to locations in the T1 intensity image that are
at a particular intensity threshold and tends to smooth thersurface. By limiting this operation to only previously modified points and setting the maximum number of iterations
to 40, the amount of smoothing of the surface is minimized.

Quantitative Analysis of Results
For the first two data sets (the phantom and 12 scans of thesame brain), corrected surfaces were compared to the ‘‘gold
standard’’ surface using the Hausdorff distance, the mean
distance error, distance error histograms (for the second setonly), and percentage of outlier reduction (defined below).
The mean distance error de is the average minimum distance
between a set of points X and a surface S, whereas the Hausdorff distance is the maximum value within the minimum
distance set. The minimum distance function d(p,S) between
a point p[X and the surface S can be defined as:
dðp; XÞ ¼ min||p & p0 ||;
0
p 2S

(5)

where p0 is a point on the surface S. The mean distanceerror is then defined as follows
de ¼

1 X
dðp; SÞ;Np p2X(6)

where Np is the number of points in the set of points X.
The outlier reduction percent represents the fraction ofpoints that remain above a distance error threshold set in
the uncorrected brain surface, such that:
$
o%
Nt Np
" 100;(7)
OP ¼ 1 & o +
Nt Np

9 r

r Yotter et al. r

Figure 7.
Topological correction using spherical harmonics relies on a union the low-pass filtered surface (c). The union of the two surfaces
between a high-frequency (l ¼ 1,024) and low-pass filtered (ll ¼ (d) optimizes detail retention and topological error correction. T164) surface. The original surface (a) contains three holes (red postcorrection smooths out the holes and handles (e), so that thearrows), two handles (yellow arrows), and two geometrical arti- corrected surface closely matches the ideal surface (f). [Color figfacts (blue arrows). The high-frequency reconstruction (b) replaces ure can be viewed in the online issue, which is available attopological defects with sharp spikes but retains more detail than www.interscience.wiley.com.]
where N$t and Nt are the number of points whose distanceerrors are above threshold in the original brain surface
and the corrected surface, respectively, and N$p and Np are
the total numbers of mesh points in the original brain surface and the corrected surface, respectively. The thresholdis set to include the top 5% points with the largest distanceerrors in the uncorrected brain surface. A value of 100%
indicates that all outliers have been removed, whereas avalue of 0% indicates no improvement.
For the third dataset (10 scans of healthy control subjects), each topological defect was examined by an expert(R. D.) and given a score (þ/0/& for accurate/partiallyaccurate/inaccurate) based on the quality of the final surface. Here, a score of ‘‘0’’ or ‘‘partially accurate’’ meansthat the cut (or fill) was correct, but the corrected surface
did not completely match the contours of the inner (orouter) surface in the original mesh.

r

RESULTS AND DISCUSSION
Topological Phantom
Although a topological phantom does not resemble a
brain surface, topological defects in such a synthetic dataset should still be accurately repaired if the approach isfundamentally valid. Such datasets provide a way toassess accuracy when ground truth is known because it isspecified analytically.
As shown in Figure 7, the topology is accurately corrected in a 60-mm3 cubic phantom. Before T1 postcorrection, the topological correction algorithm shows goodresults filling holes and cutting handles. Furthermore,
topologically correct but ‘‘inaccurate’’ geometrical artifacts
are greatly reduced. The forward Hausdorff distance is
1.71 mm (due to the remaining geometrical artifacts (bluearrows) that originally deviated by 12 mm from the ideal
surface), and the forward mean distance error is 4.7 lm.

10 r

r Cortical Surface Topology Correction r

TABLE I. Topological correction using spherical harmonics reduces the mean distance error and the proportionof the top 5% outliers in comparison with the averaged ‘‘gold standard’’ brain surface
Mean distance error (mm)Hausdorff distance (mm)Correction

Forward

ReverseForward

ReverseOutlier reduction(top 5%)Total
self-intersections

None
FreeSurfer
SPH

0.749 (%0.007)
0.641 (%0.007)
0.537 (%0.007)

0.509 (%0.008)
0.563 (%0.007)
0.539 (%0.008)

10.949 (%0.022)
7.910 (%0.056)
6.855 (%0.044)

4.332 (%0.032)
5.270 (%0.041)
5.922 (%0.034)

—
65.4%
85.1%

13
21
3

Generally, the forward Hausdorff distance and mean distance error metrics measure how well the reconstructed surface matches the‘‘gold standard’’ mesh; the reverse Hausdorff distance and mean distance error indicate how much information in the ‘‘gold standard’’mesh was retained in the reconstruction.

After T1 postcorrection, there are still some roundededges in the phantom that deviate from the target surface.
These highlight a key feature of the spherical harmonics
approach, in that sharp edges are not preserved. As sharpedges do not exist in the folded surface of the brain, thismight be considered to be an advantage of this approach,because it would correct these edges that must be due tonoise. For surfaces with sharp edges, however, it would be
important to choose a T1 postcorrection algorithm that
does not contain a curvature limitation, as does the software package that we chose to use.

Twelve Scans of the Same Brain
The corrected surfaces were all homeomorphic to asphere, verified by calculating the Euler characteristic.
Quantitatively, our spherical harmonic method generated
a corrected surface that has a lower distance error, an
improved distance error histogram, and fewer outliers
compared with the uncorrected surface, when both surfaces
are compared with the ideal surface (Table I and Fig. 8). The
corrected surfaces also contained fewer self-intersections.
Visual inspection reveals that most corrections werecloser to the ‘‘gold standard’’ averaged brain surface compared with the original uncorrected surface. Typical corrections are shown in Figure 9. By mapping the distanceerror across the surface of a single brain, the distance error
located in nondefect regions closely matched the distanceerror in the uncorrected surface (see Fig. 10). In previously
published algorithms, large defects were sometimes problematic for the correction processes. However, ourapproach accurately corrects large defects. Figure 11 shows
typical corrections in the ventricular region and the nonbrain area around the orbitofrontal cortex.
Because the averaged ‘‘gold standard’’ brain may besmoothed during the averaging process, it could be argued
that the smoothing due to spherical harmonics is the primary cause of the reduced distance error metrics. However, by measuring the distance error from the original
uncorrected surface, it is clear that the spherical harmonics
reconstruction is highly representative of the original uncorrected surface, despite the fact that the surface has been
reparameterized (Table II). This is further shown by the
low mean distance errors. As expected, the reverse Haus-rdorff distance is large because of the removal of large
defects and artifacts. Changes in these areas are also
responsible for an increase in the reverse mean distance error.
The full topology correction process using spherical harmonics (including identification of defects, smoothing,
patching, and so on) requires !6–8 min on a 2.4-GHz
iMac for a mesh that contains 150,000 vertices. This
excludes the initial spherical mapping generation and T1
postcorrection; the latter step requires !8 min of additional processing time. The resulting mesh contains
approximately the same number of vertices. As a basis ofcomparison, the topology correction process in the standard FreeSurfer pipeline takes approximately an hour; however, as no source code is available, it is unknown
whether other operations besides topological correction are
performed during this processing step.

Figure 8.
The forward distance error histogram is improved for surfaces
corrected using the spherical harmonic method. The histograms
are an average compiled for 24 hemispheres. Compared with theuncorrected surfaces, there are more points with low or zero
distance error for the spherical harmonic corrected surfaces.
[Color figure can be viewed in the online issue, which is availableat www.interscience.wiley.com.]

11 r

r Yotter et al. r

Figure 9.
A typical correction of topological defects by spherical harmonics is shown. Holes (a) arefilled (b), and handles (c) are cut (d). [Color figure can be viewed in the online issue,
which is available at www.interscience.wiley.com.]

Figure 10.
Away from topological defects, surface information remains unchanged. The uncorrected surfaces(a and b) are almost identical to the spherical harmonics corrected surfaces (c and d), except inregions near topological defects. The distance in mm is the error between these surfaces andthe ‘‘gold standard’’ averaged surface. [Color figure can be viewed in the online issue, which isavailable at www.interscience.wiley.com.]

r

12 r

r Cortical Surface Topology Correction r

Figure 11.
Large defect regions in the ventricular region are filled, whereas large defects in the nonbrain
area around the orbitofrontal cortex are cut. (a) In the original uncorrected surface, the ventricular defect is marked in red and the defect around the orbitofrontal cortex is marked in yellow.
(b) Spherical harmonics correction fills the ventricular defect and cuts the defect around theorbitofrontal cortex. [Color figure can be viewed in the online issue, which is available atwww.interscience.wiley.com.]the errors arose from incorrect labeling for fill/cut operations or incorrect bisecting of the defect points due to an
unusual geometry of the defect. As a point of comparison,
the standard FreeSurfer pipeline accurately repaired 94%
of the defects, of which 3% were only partially correct.

Ten Scans of Healthy Controls
Ten additional brain surfaces were processed, and eachdefect correction was rated by an expert. The expert wasprovided with two corrections per hemisphere that had
been randomly grouped with the goal of obtaining an
unbiased correction score; however, when questioned
afterward, it was clear to the expert which brain surface
was corrected using which method, due primarily to thecorrections performed for the ventricular areas. Sphericalharmonic correction accurately repaired topological defects
in 93% of the cases (score of þ or 0), of which 9% were
only partially correct. For the remaining 7% (score of &),CONCLUSION
Topology correction using spherical harmonics successfully removes topology defects and reduces geometricalartifacts. The approach is also valid for large defects that
contain multiple holes and handles (Euler characteristic

TABLE II. The mean distance errors for surfaces corrected using topological correction are small whencompared with the original uncorrected surface
Mean distance error (mm)Hausdorff distance (mm)Correction

Forward

ReverseForward

ReverseFreeSurfer
SPH

0.227 (%0.006)
0.011 (%0.000)

0.408 (%0.005)
0.309 (%0.002)

4.186 (%0.037)
3.031 (%0.018)

10.592 (%0.075)
11.025 (%0.037)

However, the reverse Hausdorff distance is large because of the removal of defects.

r

13 r

r Yotter et al. rless than &10). It accurately corrects most of the topological defects, and considering that this is a new approach, it
is expected that, with future refinements, the correction
accuracy may increase in the future.
One source of error is in misclassification of holes andhandles, resulting in holes being cut and handles being filled.
The approach presented here relies on the T1-weighted
image intensity in the center of the defect to differentiate
between the two cases. A more complicated algorithm thattakes into account the location of the defect on the surfaceof the brain, the T1-weighted intensity image surrounding
the defect, and even curvature values may improve classification accuracy. One important tradeoff is the computation time required for a more complicated approach,which may offer limited improvement in accuracy.
Another source of improvement would be in the T1based postcorrection process. In this implementation, we
chose to use a publicly available software package that
does not account for the method with which the surfaces
were originally corrected. A global T1 postcorrection algorithm that is compatible with spherical harmonics (forexample, by making assumptions about the initial smoothness of the reconstructed surface or by examining the original uncorrected surface) could increase the accuracy ofthis correction step without a noticeable sacrifice in computation speed. Furthermore, the T1-based postcorrection
was only applied to the topological defects; it is expectedthat a T1-based postcorrection applied globally to theentire surface would further increase the accuracy of thereconstructed surface.
The topological errors (especially the large defect in thearea around the orbitofrontal cortex) arose initially fromincorrect segmentation of tissue types. The accuracy oftissue segmentation could be improved by using tissue
registration [Ashburner and Friston, 2005]. However,depending on the surface construction algorithm, theimproved segmentation accuracy may not help to reduce
the size and number of handles across gyri, because thespaces between gyri would still be relatively narrow.
In conclusion, we found that topological correction
using spherical harmonics offers a fast, accurate method tocorrect topological defects in cortical surfaces used in brain
mapping. We verified this both qualitatively and quantitatively. As the approach is fairly straightforward, it lends
itself to inclusion in a full processing pipeline for surface
data. The topology of surfaces generated using thismethod are more accurate than the original uncorrectedsurface, are homeomorphic with a sphere, and are suitablefor intersubject analysis.

REFERENCES
Abrams L, Fishkind DE, Priebe CE (2002): A proof of the sphericalhomeomorphism conjecture for surfaces. IEEE Trans Med
Imaging 21:1564–1566.

r

Abrams L, Fishkind DE, Priebe CE (2004): The generalized spherical homeomorphism theorem for digital images. IEEE TransMed Imaging 23:655–657.
Ashburner J, Friston KJ (2005): Unified segmentation. Neuroimage
26:839–851.
Ballmaier M, Kumar A, Thompson PM, Narr KL, Lavretsky H,
Estanol L, DeLuca H, Toga AW (2004): Localizing gray matter
deficits in late-onset depression using computational cortical
pattern matching methods. Am J Psychiatry 161:2091–2099.
Bazin P-L, Pham DL (2005a): Topology correction using fast
marching methods and its application to brain segmentation.
In: Duncan J, Gerig G, editors. MICCAI. Palm Springs, CA:
Springer. pp 484–491.
Bazin P-L, Pham DL (2005b): Topology preserving tissue classification with fast marching and topology templates. In: Christensen GE, Sonka M, editors. Information Processing in Medical
Imaging. Glenwood Springs, CO: Springer. pp 234–245.
Bazin P-L, Pham DL (2007a): Topology correction of segmented
medical images using a fast marching algorithm. Comput
Methods Programs Biomed 88:182–190.
Bazin P-L, Pham DL (2007b): Topology-preserving tissue classification of magnetic resonance brain images. IEEE Trans Med
Imaging 26:487–496.
Bearden CE, Erp TGMv, Thompson PM, Toga AW, Cannon TD(2007a): Cortical mapping of genotype-phenotype relationships
in schizophrenia. Hum Brain Mapp 28:519–532.
Bearden CE, van Erp TGM, Dutton RA, Tran H, Zimmermann L,
Sun D, Geaga JA, Simon TJ, Glahn DC, Cannon TD, EmanuelBS, Toga AW, Thompson PM (2007b): Mapping cortical thickness in children with 22q11.2 deletions. Cereb Cortex 17:1889–
1898.
Bischoff S, Kobbelt LP (2002): Isosurface reconstruction with topology control. In: Coquillart S, Shum H-Y, Hu S-M, editors. Proceedings of the 10th Pacific Conference on Computer Graphics
and Applications. Beijing: IEEE Computer Society. pp 246–255.
Bischoff S, Kobbelt L (2003): Sub-voxel topology control for levelset surfaces. Comput Graph Forum 22:273–280.
Bischoff S, Kobbelt L (2004): Topologically correct extraction of
the cortical surface of a brain using level-set methods. In: Tolxdorff T, Braun J, Handels H, Horsch A, Meinzer H-P, editors.
Proceedings of BVM. Berlin: Springer-Verlag. pp 50–54.
Brechbühler C, Gerig G, Kübler O (1995): Parametrization ofclosed surfaces for 3-D shape description. Comput Vis Image
Understand 61:154–170.
Caselles V, Kimmel R, Sapiro G (1997): Geodesic active contours.
Int J Comput Vis 22:61–79.
Chen L, Wagenknecht G (2006): Automated topology correction forhuman brain segmentation. In: Larsen R, Nielsen M, Sporring J,
editors. MICCAI. Copenhagen, Denmark: Springer. pp 316–323.
Christensen GE, Rabbitt RD, Miller MI (1994): 3D brain mapping
using a deformable neuroanatomy. Phys Med Biol 39:609–618.
Chung M, Shen L, Dalton K, Davidson R (2006): Multi-scale
voxel-based morphometry via weighted spherical harmonic
representation. In: Yang G-Z, Jiang T, Shen D, Gu L, Yang J,
editors. Medical Imaging and Augmented Reality. Shanghai:
Springer. pp 36–43.
Chung MK, Dalton KM, Li S, Evans AC, Davidson RJ (2007):
Weighted Fourier series representation and its application toquantifying the amount of gray matter. IEEE Trans Med Imaging 26:566–581.

14 r

r Cortical Surface Topology Correction r

Dale AM, Fischl B, Sereno MI (1999): Cortical surface-based analysis. I. Segmentation and surface reconstruction. Neuroimage
9:179–194.
Delingette H (1999): General object reconstruction based on simplex meshes. Int J Comput Vis 32:111–146.
Delingette H, Montagnat J (2001): Shape and topology constraints
on parametric active contours. Comput Vis Image Understand
83:140–171.
Desai R, Liebenthal E, Possing ET, Waldron E, Binder JR (2005):
Volumetric vs. surface-based alignment for localization of auditory cortex activation. Neuroimage 26:1019–1029.
Fischl B, Sereno MI, Dale AM (1999a): Cortical surface-based analysis. II. Inflation, flattening, and a surface-based coordinate
system. Neuroimage 9:195–207.
Fischl B, Sereno MI, Tootell RBH, Dale AM (1999b): High-resolution intersubject averaging and a coordinate system for thecortical surface. Hum Brain Mapp 8:272–284.
Fischl B, Liu A, Dale AM (2001): Automated manifold surgery: Constructing geometrically accurate and topologically correct modelsof the human cerebral cortex. IEEE Trans Med Imaging 20:70–80.
Formisano E, Esposito F, Di Salle F, Goebel R (2004): Cortex-based
independent component analysis of fMRI time series. Magn
Reson Imaging 22:1493–1504.
Gaser C, Luders E, Thompson PM, Lee AD, Dutton RA, Geaga
JA, Hayashi KM, Bellugi U, Galaburda AM, Korenberg JR,
Mills DL, Toga AW, Reiss AL (2006): Increased local
gyrification mapped in Williams syndrome. Neuroimage 33:
46–54.
Gerig G, Styner M, Jones D, Weinberger D, Lieberman J (2001):
Shape analysis of brain ventricles using SPHARM. In: Mathematical Methods in Biomedical Image Analysis (MMBIA).
Kauai, HI: IEEE. pp 171–178.
Guskov I, Wood ZJ (2001): Topological noise removal. In: Proceedings of Graphics Interface. Ottawa, Ontario, Canada: Canadian
Information Processing Society. pp 19–26.
Gutman B, Wang Y, Lui LM, Chan TF, Thompson PM (2006): Hippocampal surface analysis using spherical harmonic function
applied to surface conformal mapping. In: International Conference on Pattern Recognition (ICPR). Hong Kong, China:
IEEE Computer Society. pp 964–967.
Han X, Xu C, Braga-Neto U, Prince JL (2002): Topology correctionin brain cortex segmentation using a multiscale, graph-basedalgorithm. IEEE Trans Med Imaging 21:109–121.
Han X, Pham DL, Tosun D, Rettmann ME, Xu C, Prince JL (2004):
3RUISE: Cortical reconstruction using implicit surface evolution. Neuroimage 23:997–1012.
Healy DM, Rockmore DN, Moore SSB (1996): FFTs for the
2-Sphere—Improvements and Variations. Technical Report
PCS-TR96–292, Dartmouth College.
Hinds OP, Rajendran N, Polimeni JR, Augustinack JC, Wiggins G,
Wald LL, Diana Rosas H, Potthast A, Schwartz EL, Fischl B(2008): Accurate prediction of V1 location from cortical folds ina surface coordinate system. Neuroimage 39:1585–1599.
Hurdal MK, Stephenson K (2009): Discrete conformal methods forcortical brain flattening. Neuroimage 45:S86–S98.
Jaume S, Rondao P, Macq B (2005): Open topology: A toolkit forbrain isosurface correction. In: Duncan J, Gerig G, editors.
MICCAI. Palm Springs, CA: Springer.
Karaçali B, Davatzikos C (2003): Topology preservation and regularity in estimated deformation fields. In: Taylor C, Noble JA,
editors. Information Processing in Medical Imaging. Ambleside, UK: Springer. pp 426–437.

rKass M, Witkin A, Terzopoulos D (1988): Snakes: Active contour
models. Int J Comput Vis 1:321–331.
Kelemen A, Szekely G, Gerig G (1999): Elastic model-based segmentation of 3-D neuroradiological data sets. IEEE Trans Med
Imaging 18:828–839.
Kostelec PJ, Rockmore DN (2004): S2kit: A Lite Version of SpharmonicKit. Hanover, NH: Department of Mathematics, Dartmouth College.
Kostelec PJ, Maslen DK, Dennis M, Healy J, Rockmore DN (2000):
Computational harmonic analysis for tensor fields on the twosphere. J Comput Phys 162:514–535.
Kriegeskorte N, Goebel R (2001): An efficient algorithm for topologically correct segmentation of the cortical sheet in anatomical MR volumes. Neuroimage 14:329–346.
Lachaud J-O, Montanvert A (1999): Deformable meshes with automated topology changes for coarse-to-fine three-dimensionalsurface extraction. Med Image Anal 3:187–207.
Lerch JP, Pruessner JC, Zijdenbos A, Hampel H, Teipel SJ, Evans
AC (2005): Focal decline of cortical thickness in Alzheimer’s
disease identified by computational neuroanatomy. Cereb Cortex 15:995–1001.
Lin JJ, Salamon N, Lee AD, Dutton RA, Geaga JA, Hayashi KM,
Luders E, Toga AW, Engel J Jr, Thompson PM (2007): Reduced
neocortical thickness and complexity mapped in mesial temporal lobe epilepsy with hippocampal sclerosis. Cereb Cortex
17:2007–2018.
Lorensen WE, Cline HE (1987): Marching cubes: A high resolution 3D surface construction algorithm. In: Stone MC, editor.
SIGGRAPH: 14th Annual Conference on Computer Graphics
and Interactive Techniques. Anaheim, CA: ACM. pp 163–
169.
Lui L, Wang Y, Chan T, Thompson P (2006): A landmark-based
brain conformal parameterization with automatic landmark
tracking technique. In: Larsen R, Nielsen M, Sporring J,
editors. Medical Image Computing and Computer-Assisted
Intervention—MICCAI 2006. Copenhagen, Denmark: SpringerVerlag, New York. pp 308–315.
Lui LM, Wang Y, Chan TF, Thompson P (2007): Landmark constrained genus zero surface conformal mapping and its application to brain mapping research. Appl Numer Math 57:847–858.
MacDonald D, Avis D, Evans AC (1994): Multiple surface identification and matching in magnetic resonance images. In: Robb
RA, editor. Visualization in Biomedical Computing 1994. Rochester, MN: SPIE. pp 160–169.
MacDonald D, Kabani N, Avis D, Evans AC (2000): Automated
3-D extraction of inner and outer surfaces of cerebral cortex
from MRI. Neuroimage 12:340–356.
Mangin J-F, Frouin V, Bloch I, Régis J, Ló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.
McInerney T, Terzopoulos D (1996): Deformable models in medical image analysis: A survey. Med Image Anal 1:91–108.
McInerney T, Terzopoulos D (2000): T-snakes: Topology adaptive
snakes. Med Image Anal 4:73–91.
Montagnat J, Delingette H, Ayache N (2001): A review of deformable surfaces: Topology, geometry and deformation. Image Vis
Comput 19:1023–1040.
Narr KL, Bilder RM, Toga AW, Woods RP, Rex DE, Szeszko PR,
Robinson D, Sevy S, Gunduz-Bruce H, Wang Y-P, DeLuca H,
Thompson PM (2005): Mapping cortical thickness and gray15 r

r Yotter et al. r

matter concentration in first episode schizophrenia. Cereb
Cortex 15:708–719.
Qiu A, Younes L, Wang L, Ratnanather JT, Gillepsie SK, Kaplan
G, Csernansky J, Miller MI (2007): Combining anatomical
manifold information via diffeomorphic metric mappings forstudying cortical thinning of the cingulate gyrus in schizophrenia. Neuroimage 37:821–833.
Ségonne F, Grimson E, Fischl B (2003): Topological correction of
subcortical segmentation. In: Ellis RE, Peters TM, editors. MICCAI. Montreal, Canada: Springer. pp 695–702.
Ségonne F, Grimson E, Fischl B (2005): A genetic algorithm for the
topology correction of cortical surfaces. In: Christensen GE,
Sonka M, editors. Proceedings of Information Processing in
Medical Imaging. Glenwood Springs, CO: Springer. pp 393–405.
Segonne F, Pacheco J, Fischl B (2007): Geometrically accurate topology-correction of cortical surfaces using nonseparating
loops. IEEE Trans Med Imaging 26:518–529.
Shattuck DW, Leahy RM (2001): Automated graph-based analysis
and correction of cortical volume topology. IEEE Trans Med
Imaging 20:1167–1177.
Shen L, Chung MK (2006): Large-scale modeling of parametricsurfaces using spherical harmonics. In: 3DPVT. Chapel Hill,
NC: IEEE Computer Society. pp 294–301.
Shenton ME, Gerig G, McCarley RW, SzÈkely Gb, Kikinis R(2002): Amygdala-hippocampal shape differences in schizophrenia: The application of 3D shape models to volumetric
MR data. Psychiatry Res 115:15–35.
Thompson PM, Toga AW (1996): A surface-based technique for
warping three-dimensional images of the brain. IEEE TransMed Imaging 15:402–417.
Thompson PM, Hayashi KM, Sowell ER, Gogtay N, Giedd JN,
Rapoport JL, de Zubicaray GI, Janke AL, Rose SE, Semple J,
Doddrell DM, Wang Y, van Erp TGM, Cannon TD, Toga AW(2004): Mapping cortical change in Alzheimer’s disease, brain
development, and schizophrenia. Neuroimage 23:S2–S18.
Thompson PM, Lee AD, Dutton RA, Geaga JA, Hayashi KM, Eckert
MA, Bellugi U, Galaburda AM, Korenberg JR, Mills DL, Toga AW,
Reiss AL (2005): Abnormal cortical complexity and thickness profiles mapped in Williams syndrome. J Neurosci 25:4146–4158.
Thompson PM, Bartzokis G, Hayashi KM, Klunder AD, Lu PH,
Edwards N, Hong MS, Yu M, Geaga JA, Toga AW, Charles C, Perkins DO, McEvoy J, Hamer RM, Tohen M, Tollefson GD, Lieberman JA, Group THS (2009): Time-lapse mapping of corticalchanges in schizophrenia with different treatments. Cereb Cortex
19:1107–1123.
Tosun D, Reiss AL, Lee AD, Dutton RA, Hayashi KM, Bellugi U,
Galaburda AM, Korenberg JR, Mills DL, Toga AW, Thompson
PM (2006): Use of 3-D cortical morphometry for mapping
increased cortical gyrification and complexity in Williams syndrome. In: IEEE International Symposium on Biomedical Imaging (ISBI2006), Arlington, VA, USA.
Van Essen DC (2004): Surface-based approaches to spatial localization and registration in primate cerebral cortex. Neuroimage
23:S97–S107.
Van Essen DC (2005): A population-average, landmark- and surface-based (PALS) atlas of human cerebral cortex. Neuroimage
28:635–662.
Van Essen DC, Drury HA (1997): Structural and functional analyses of human cerebral cortex using a surface-based atlas.
J Neurosci 17:7079–7102.
Van Essen DC, Maunsell JHR (1980): Two-dimensional maps ofthe cerebral cortex. J Comp Neurol 191:255–281.

rWang Y, Chiang M-C, Thompson PM (2005): Automated surface
matching using mutual information applied to Riemannsurface structures. In: Duncan J, Gerig G, editors. Medical
Image Computing and Computer-Assisted Intervention—MICCAI 2005. Palm Springs, CA: Springer Berlin/Heidelberg.
pp 666–674.
Wang Y, Gu X, Chan TF, Thompson PM, Yau S-T (2008): Conformal slit mapping and its applications to brain surface parameterization. In: Metaxas DN, Axel L, Fichtinger G, Szekely G,
editors. International Conference on Medical Image Computing
and Computer Assisted Intervention—MICCAI 2008. New
York City, NY: Springer Berlin/Heidelberg. pp 585–593.
Wang Y, Dai W, Gu X, Chan TF, Yau S-T, Toga AW, Thompson
PM (2009): Teichmüller shape space theory and its application
to brain morphometry. In: Yang G-Z, Hawkes DJ, Rueckert D,
Noble JA, Taylor CJ, editors. Medical Image Computing andComputer Assisted Intervention (MICCAI2009). London, UK:
Springer Berlin/Heidelberg. pp 131–140.
Wang Y, Zhang J, Gutman B, Chan TF, Becker JT, Aizenstein HJ,
Lopez OL, Tamburo RJ, Toga AW, Thompson PM (2010): Multivariate tensor-based morphometry on surfaces: Application
to mapping ventricular abnormalities in HIV/AIDS. Neuroimage 49:2141–2157.
Wood Z, Hoppe H, Desbrun M, Schröder P (2004): Removing excess
topology from isosurfaces. ACM Trans Graph 23:190–208.
Xiao H, Chenyang X, Prince JL (2003): A topology preserving levelset method for geometric deformable models. IEEE Trans Pattern Anal Mach Intell 25:755–768.
Xu C, Pham DL, Rettmann ME, Yu DN, Prince JL (1999): Reconstruction of the human cerebral cortex from magnetic resonance images. IEEE Trans Med Imaging 18:467–480.
Yotter RA, Thompson PM, Gaser C: Algorithms to improve there-parameterization of spherical mappings of brain surface
meshes. J Neuroimag, in-press.
Zhou Q-Y, Ju T, Hu S-M (2007): Topology repair of solid models
using skeletons. IEEE Trans Visual Comput Graph 13:675–685.

APPENDIX: AREA SMOOTHING ALGORITHM
Self-intersections are resolved by smoothing the meshlocally. First, the self-intersection region is labeled. Edge points
are then ‘‘anchored’’ such that the points within the edge
points are the only ones that are smoothed. Given an inner
point p, each triangle adjacent to the point is assigned a weightthat is used to adjust the position of the point. The center of
each triangle tj containing point p is found using the formula:&x þ x þ x y þ y þ y z þ z þ z '
p
r
s
p
r
s
p
r
s
cj ¼
;
;
;
(8)
3
3
3
where r and s are the other two points in the triangle. The
weight wj of triangle tj is set such that:
aj
wj ¼
;
(9)
Ap
where Ap is the total area of all triangles containing point
p. The position of point p can then be updated as follows:

16 r