Technology Update
Algorithms to Improve the Reparameterization of Spherical
Mappings of Brain Surface Meshes
Rachel A. Yotter, PhD, Paul M. Thompson, PhD, Christian Gaser, PhD
From the Department of Psychiatry, Friedrich-Schiller University, Jena, Germany (RAY, CG); Laboratory of Neuro Imaging, Department of Neurology, Division of Brain Mapping,
UCLA School of Medicine, Los Angeles, California (PMT).

ABSTRACTA spherical map of a cortical surface is often used for improved brain registration, foradvanced morphometric analysis (eg, of brain shape), and for surface-based analysis offunctional signals recorded from the cortex. Furthermore, for intersubject analysis, it isusually necessary to reparameterize the surface mesh into a common coordinate system.
An isometric map conserves all angle and area information in the original cortical mesh;however, in practice, spherical maps contain some distortion. Here, we propose fast newalgorithms to reduce the distortion of initial spherical mappings generated using one ofthree common spherical mapping methods. The algorithms iteratively solve a nonlinear
optimization problem to reduce distortion. Our results demonstrate that our correction
process is computationally inexpensive and the resulting spherical maps have improved
distortion metrics. We show that our corrected spherical maps improve reparameterizationof the cortical surface mesh, such that the distance error measures between the originaland reparameterized surface are significantly decreased.

Keywords: Spherical mapping, surface
mesh, distortion, intersubject analysis,
parameterization.
Acceptance: Received November 12,
2009, and in revised form February 23,
2010. Accepted for publication February
27, 2010.
Correspondence: Address correspondence to Rachel A. Yotter, Department of
Psychiatry, Friedrich-Schiller University,
3 Jahnstrasse, 07743 Jena, Germany.
E-mail: rachel.yotter@uni-jena.de
Conflict of Interest: The authors declareno conflicts of interests.
J Neuroimaging 2010;XX:1-14.
DOI: 10.1111/j.1552-6569.2010.00484.x

Introduction
Most brain MRI scanning protocols acquire volumetric data
about the anatomy of the subject. However, it is sometimes
desirable to conduct analyses that focus exclusively on the geometry of the cortical surface. This type of analysis may be
conducted directly in volume space or by first generating asurface mesh from the volumetric data.1-3 A surface mesh is
useful for computing average surfaces, and the surface coordinates are useful for making comparisons of data across subjects that take into account the cortical folding pattern. Furthermore, a surface mesh makes it easier to perform some
types of shape analyses, for example, spherical harmonic analysis, Laplace–Beltrami (LB) eigenmaps (another form of shapeanalysis), gyrification indices that measure surface complexity in 3-dimensional (3D), and complexity analysis that regrids the cortex.4-6 To proceed with these analyses, a surface
mesh that accurately represents cortical anatomy must first be
generated.
Surface analyses have special advantages that are not present
using volumetric data alone. For instance, brain surface meshes
have been shown to increase the accuracy of brain registrationcompared with linear Talairach registration.7-10 Brain surface
meshes also permit new forms of analyses, such as gyrification
indices that measure surface complexity in 3D,6,11,12 cortical
thickness,13,14 and data compression and searching in mining

databases.15 Furthermore, inflation or spherical mapping of thecortical surface mesh raises the buried sulci to the surface sothat mapped functional activity in these regions can be easily visualized. However, the primary purpose behind sphericalmapping of the brain surface mesh is to generate a common
coordinate system for intersubject analysis. To place a subject
into a common coordinate system, the surface mesh usuallymust be reparameterized.
One valuable type of spherical map is a pseudo-isometricmap, which attempts to conserve both angle and area information from the cortical surface mesh. Any mapping may retain area information by computing the Jacobian; however,a pseudo-isometric map preserves the local areas as much as
possible, such that measuring relative areas in the sphere is essentially the same as doing so on the original surface, without
having to know the Jacobian of the mapping.
Because of the highly folded nature of the cortical surface,
however, it is not possible to obtain a perfectly isometric spherical map of the cortical surface mesh. In addition, depending on
the priorities of the application, either angle or area preservationmust be optimized at the expense of the other. However, it has
not been explored which metric most influences the accuracyof reparameterization. It is hypothesized that area distortion ismore important than angular distortion. Because there are an
infinite number of solutions for an equiareal mapping, it can◦

Copyright C 2010 by the American Society of Neuroimaging

1

be advantageous to choose a solution that also keeps angulardistortion at a minimum.16
Previously, most efforts to minimize distortion concurredwith the algorithms used to generate the initial spherical mapping. For instance, one can create a conformal, or anglepreserving, spherical mapping. There is some evidence that
conformal maps may exist in the V1 visual cortex for processing visual data.17 Also, it may be easier to solve partial differential equations on a grid if the grid is conformal, such as for
signal smoothing or for surface-to-surface registration applications. A conformal map of a surface mesh can be generated byeither solving a partial differential equation that involves the
LB operator of the surface coordinates,18,19 often using sulcal
features as explicit landmark constraints,20-22 or by minimizing
the harmonic energy of the mapping to the sphere using variational methods.23-26 Another quasi-conformal mapping procedure uses circle packing as a discrete approximation to theclassical continuous mapping theory.27 More recently, some researchers have used a branch of differential geometry known
as the exterior calculus to build conformal grids directly on surfaces, using holomorphic differential one-forms.28 Conformal
grids on surfaces may also be generated using diverse methods
such as the Ricci flow method,29,30 algebraic function theory,31slit mapping,32 and cohomology theory.28 The resulting conformal grids may be used for shape analysis using Teichmuller
space theory or tensor-based morphometry on the surface
metric.28,33
Alternatively, specific landmark features, such as sulcal
curves lying in the cortex, can be forced to map to designatedlocations on the sphere, using either covariate partial differential equations derived from linear elasticity,34 level set methods and implicit function theory,35 or approximately, by using
mutual information defined on surfaces.24,36 Although someapproaches have computed surface-to-surface registrations byintermediate mappings of both surfaces to the sphere, an important class of methods has performed direct surface-to-surface
registration.37 In addition, methods using covariant partial differential equations attempt to produce surface-to-surface registrations that, in the continuous case at least, are independentof the intermediate spherical mappings used to impose grids
on the surfaces.34,38 In these mappings, a flow vector field inthe surface coordinates is developed in which the differential
operators are made covariant to the surface metrics, leading tosurface registrations that are provably independent of the way
the surfaces are gridded, so long as the surfaces are sufficiently
finely sampled.
Within the grid generation field, there exist penalty functionals that will conserve a linear combination of length, orthogonality, and area so that semi-isometric flat maps of sectionsof the cortical surface can be generated. However, unless thesurface is developable, it is not possible to conserve all threefunctionals. The brain surface mesh may be iteratively inflated
by minimizing an energy functional that is related to metric distortion or geodesic lengths, creating a spherical mapping withlow metric distortion.39,40 Alternately, an initial simple mapping may be optimized by finding the minimum solution to anerror functional that accounts for both area and angular distortion41 or area and length distortion.42 Semi-isometric flat maps

2

Journal of Neuroimaging Vol XX No X XX 2010of sections of the cortical surface may also be generated using a constrained optimization problem that accounts for bothangles and area.43 Such approaches are analogous to grid generation problems for surfaces using the compound functional
method, and can be expanded to n dimensions in the general case.44 In fact, equiareal mappings could be considered
a subcase achievable by using only one of the three penalty
functionals.45
In this paper, we propose to apply three algorithms to optimize spherical maps with respect to area distortion (and, secondarily, local metric distortion). The algorithms are uniqueand the problem is solved as an optimization problem to arrive at the minimum (distortion) solution. The method starts
with initial spherical maps generated using a variety of previously developed methods. As a first step, the area distortionis minimized by finding the minimum solution to a nonlinear optimization problem based on an algorithm (the distortalgorithm, Eq 8) that is directly related to area distortion. After the minimum solution is found, the spherical map is further optimized using more selective variants of all three algorithms. The resulting spherical map has improved distortion
metrics compared to the original map. With respect to reparameterization of the spherical map into a spherical coordinate
system, the optimized spherical maps result in a reparameterization that has significantly lower distance error values compared to the original map. The 3D distance error may be important in the field of morphometry, because if the original
surface is not accurately approximated during reparameterization, some statistical power may be lost for shape and area
analysis.

Methods
Minimizing Distortion
The definition of distortion is by no means established. Becauseof the nonheterogeneity of the distortion metrics used, it is feasible to select one metric to optimize. However, the optimizationof one metric usually involves the degradation of another distortion metric, and this information is lost if other distortion
metric values are not reported. Reporting a diverse set of distortion metrics for all resulting spherical maps can circumvent
this.
Three standard distortion metrics include metric distortion,area distortion, and angle distortion.44 The metric distortion
metric is a universal measure that incorporates both area andangular distortion, whereas the area and angular distortion metrics selectively measure distortion in areas or angles, respectively. In the following discussion, the surface mesh is assumed
to consist of triangles, but the metrics can be generalized to
non-triangular meshes as well.
Metric distortion is a measure based on differences in distances between vertices. Because two triangles with three congruent sides must by definition be isometric, any deviation corresponds to distortion. Assume that there is a vertex vo,i in theoriginal brain surface and its matching vertex on the unit spherevi . Let do,ij and dij be the distances between the ith and jth vertices in the original brain and the spherical map, respectively.

The metric distortion is then defined as follows:


Nv
! $
$%
1
1 !
$do,ij − dij $ dij ,!met =
Nv
n(i)
i=1

(1)

j ∈n(i)where Nv is the total number of vertices in the mesh and n
are the neighboring vertices.39,46 The metric distortion metric
could be applied locally such that the set of neighboring vertices
n(i) includes only the 1-neighbors of the ith vertex, or it can beapplied globally such that the neighboring vertices n(i) include
all other vertices in the mesh; despite whether the measure is
global or local, the optimal value signifying true isometry is
zero. For global metric distortion calculations, it is necessary to
select a random set of 1,000 vertices to reduce the computation
time. For is subset of vertices, the same calculation shown inEquation 1 is applied.
Area distortion examines the areal difference for each triangle in the mesh. Assume that there is a triangle to,i in theoriginal brain surface and its spherical map equivalent ti . We
define area distortion of triangle ti as follows:()
Ao a i.
(2)εi,area = log10
A ao,iHere, Ao and A are the total surface area of the brain meshand the unit sphere mesh, whereas ai and ao,i are the areas of
the triangles to,i in the brain mesh and ti in the unit sphere mesh,respectively. The areas are normalized with respect to the total
surface area of the mesh, allowing meshes of different physicaldimensions to be compared. The area distortion metric ε i has
some interesting properties, in that it will create a distributionthat is centered around zero, and a value of +1 means that atriangle is scaled by a factor of 10 compared to its equivalent
triangle in the original mesh. The average area distortion of theentire mesh is defined as follows:N!area =t$$1 !
$εi,area $,
Nt

(3)

i=1

where Nt is the total number of triangles in the mesh.
To measure angle distortion, each of the three angles ∠ai ,∠bi , and ∠ci in triangle ti need to be calculated and comparedto the original three angles ∠ao,i , ∠bo,i , and ∠co,i in triangle to,i .
We then define the angular distortion for triangle ti as follows:
$
$ $
$ $
$
(4)
εi,angle = $∠ai − ∠ao,i $ + $∠b i − ∠b o,i $ + $∠c i − ∠c o,i $ .

For this metric, the ideal value (eg, no angular distortion)is zero. The average angle distortion can then be calculated asfollows:!angle =

Nt
$
$
1 !
$εi,angle $.
3Nt

(5)

i=1

Weighted Algorithms for Area Distortion Correction
Weighted algorithms are applied directly upon the initial spherical map. Given a vertex v in the spherical map, each triangle
adjacent to the vertex can be assigned a weight that can be used
to adjust the position of the vertex. The center of each triangle

tj containing vertex v is found using the formula:)(
xv + xr + xs y v + yr + y s z v + z r + z s
,
,
,
cj =
3
3
3

(6)

where r and s are the other two vertices in the triangle. One
possible weighting scheme adjusts the vertex v so that the areas
are more evenly distributed, such that the weight wj of triangle
tj is
wj =

aj
,
Av

(7)

where Av is the total area of all triangles containing vertex v .
The position of vertex v can then be updated as follows:$=
v areaN!

wjc j,(8)

j =1

where N is the number of triangles containing vertex v . Thisweighting scheme could be considered as a form of smoothing
in that it distributes area more equally among triangles, and it
is beneficial only if the original mesh contains triangles of fairly
uniform area or the spherical map has large area distortions.
Applied indiscriminately to a mesh with nonuniform triangle
sizes, the area distortion will eventually increase. However, thealgorithm can be limited by calculating the local area distortion
post hoc, and only updating the vertex v if local area distortion is
reduced. The change in a patch of the spherical mesh after one
iteration of this algorithm on the LB conformal map is shownin Figure 1C.
A modification of the above algorithm uses a weighting
scheme based on the area distortion, as follows:
χj =

dj
.
Dv

(9)

Here, dj is the area distortion of triangle tj that includes
vertex v , and Dv is the total distortion of all of the triangles that
include vertex v . The area distortion dj of triangle tj can befound as follows:
aj
,
(10)
dj =
ao, j
where aj is the area of the triangle in the spherical map andao,j is the area of the triangle in the original mesh. The updated
vertex v $distort can be found as follows:$=
v distort

N!

χjc j.

(11)

j =1

The advantage of this algorithm over the previous one isthat it directly corrects for area distortion. It is also possible tocalculate the local area distortion post hoc, updating the locationof vertex v only if v $distort is found to reduce local area distortion.
The change in a patch of spherical mesh resulting from one
iteration of this algorithm on the LB conformal map is shownin Figure 1D.

Ray-Tracing Algorithms for Metric Distortion Correction
An inherent property of triangles is that two triangles with threecongruent sides must by default be isometric. Thus, by adjusting the length of each edge in the spherical map to be equal

Yotte et al: Algorithms to Improve Reparameterization of Spherical Mappings3

Table 1.
Correction

CoarseFine

Iterative Optimization Pipeline
Algorithm

Options

Max iterations

Distort (Eq 11)
Area (Eq 8)
Ray-trace (Eq 13)
Ray-trace (Eq 13)
Ray-trace (Eq 13)
Distort (Eq 11)
Ray-trace (Eq 13)

Selective
Selective
Large triangles only
Selective
Selective
Selective

100,000
1,000
1,000
1,000
1,000
1,000
1,000

The algorithm traces all line segments adjacent to vertex v
to update the position of the vertex such that metric distortion is
reduced. If the directional vector vj is the vector representationof a line segment in the spherical mesh that contains vertex v ,then the ray-traced vector v $ j can be represented mathematically as follows:+,*+ + +
v $j = +v j + − +v o, j + v̂ j .
(12)Here, vo,j is the vector representation of the corresponding
line segment in the original brain mesh, and v̂ j is the unit vector
along the direction of vj . The updated position of the vertex v
is then
N

$
=v+
v ray

Fig 1. This figure shows a patch from the spherical mapping processed with a variety of algorithms for one iteration only. The entire
map was processed and the patch was subsequently extracted. The
original patch (A) is first processed using the Laplace-Beltrami operator (B) to generate a conformal map. The conformal map is thenprocessed with only one iteration of the area algorithm (C), the distortion algorithm (D), the ray-tracing algorithm (E), or the ray-tracing
algorithm selectively applied to large triangles (F). Without testing fortriangle flips, the ray-tracing algorithm results in the mesh shown in(G). As a reference, a representative final patch is shown in (H).
to its corresponding edge in the original cortical surface map,the distortion in the spherical map should be reduced. Thiscan be accomplished by tracing line segments in the sphericalmesh along their vector direction such that the magnitude of the
traced vector is equivalent to the length of the corresponding
line segment in the original brain mesh.

4

Journal of Neuroimaging Vol XX No X XX 2010

1 ! $
v j,
N

(13)

j =1

where N is the total number of line segments that are adjacent
to vertex v . The ray-tracing function does not guarantee thatthe updated vertex positions lie upon the surface of the sphere,
so a final step is to project the vertex onto the surface of thesphere. The change in a patch of the spherical mesh after one
iteration of this algorithm on the conformal LB map is shownin Figures 1E–G.
Two major limitations of this algorithm need to be addressedto achieve a satisfactory solution. First, the updated vertex position can result in triangle flips (Fig 1G). This problem can beeasily circumvented by testing whether the angles between thenormals of neighboring triangles in the pre- and post-updated
spherical mesh are larger than 1◦ , and excluding the vertex
from the set of updated vertices if true (Fig 1E).
Second, the contribution to v $ray from large triangles canoverwhelm contributions from small triangles and can thus
generate large metric distortions in small triangles. This is especially valid if the spherical map has large area distortion, suchas the distortion that can be caused by the LB conformal map(which ignores area distortions to enforce angle conservation).
To compensate for this limitation, the algorithm may be selectively applied only to vertices where the local area distortion islarge (Fig 1F).
Finally, as earlier, it is also possible to calculate the local area
distortion post hoc and update the location of the vertex only if
local area distortion is reduced.

Optimization in Relation to Area Distortion
For meshes with large area distortion, empirical results demonstrate that the most effective algorithm for reducing area distortion is the distortion correction algorithm (Eq 11) (see “Results

Fig 2.

A surface is mapped to a sphere, then resampled onto regularly spaced grid. The original surface mesh (A) is mapped onto a sphere(B, blue). A regularly sampled sphere (B, red) is overlaid on top of this sphere. For each point in the regularly sampled grid (C, red), theintersecting triangle in the spherical mapping (C, blue) is found, and the barycentric coordinates within this triangle are calculated to obtain acoordinate of location within the original surface mesh.

and Discussion” section for the effect of each algorithm on
distortion). If all vertices are included, this equation is essentially a nonlinear optimization problem whose ideal solutionis achieved when the weights are equal to the inverse of thenumber of neighboring triangles. To arrive at the minimum solution, the problem can be solved iteratively, such that the surface mesh is only updated if the total area distortion of the meshis significantly reduced. It is presumed that an initial spherical
map that has more area distortion will require more iterations
to arrive at the minimum solution than a spherical mappingwith low initial area distortion.
After optimization using the distortion correction algorithm,
the other algorithms are applied to further reduce area distorTable 2.

LB
SI1
SI2

Spherical Coordinate Reparameterization and Error
Quantification
Using the original and optimized spherical maps, the surfaces
are reparameterized using a spherical coordinate system, such

Average Metric, Area, and Angular Distortion Metrics and Real Computation Time for 20 Central Surface Meshes
Local metricdistortion

CStion in the order outlined in Table 1. This order was optimized
empirically for computation speed and reduced distortion metrics. Generally, the “coarser” algorithms are applied first; thesealgorithms are “coarse” in that they result in relatively largechanges in the distortion metrics. The mesh is then optimized
with increasingly “finer” algorithms that make smaller changes
in the distortion metrics. All initial spherical maps were processed using an identical algorithm pipeline.

IC

√
√
√

Global metricdistortion

Area distortion

Angular distortion

Time (s)

Mean

!

Mean

!

Mean

!

Mean

!

Mean

!

.3762
.2372
.3856
.2372
.2615
.2235

–
−.1390–−.1484–−.0380

.3155
.2949
.2963
.2875
.2878
.2874

–
−.0206–−.0088–−.0004

.3879
.1677
.2373
.1299
.1632
.1345

–
−.2202–−.1074–−.02872.326
13.30023.17016.18616.50114.417–+10.974–−6.984–−2.084129.98780.944.4
899.3
4523.55115.2–+8651.0–+854.9–+591.7#: Change from uncorrected surface; IC: isometric correction; LB: Laplace–Beltrami; SI1: surface inflation method #1; SI2: surface inflation method #2.

Table 3.

Average Metric, Area, and Angular Distortion Metrics and Real Computation Time for 20 White Matter Surface Meshes
Local metricdistortion

WM

LB
SI1
SI2

IC

√
√
√

Global metricdistortion

Area distortion

Angular distortion

Time (s)

Mean

!

Mean

!

Mean

!

Mean

!

Mean

!

.3932
.2230
.3506
.2239
.1988
.1983

–
−.1702–−.1267–−.0005

.3285
.3022
.3079
.2993
.3066
.3056

–
−.0263–−.0086–−.0010

.4235
.1569
.2199
.1198
.0769
.0760

–
−.2466–−.1001–−.00092.044
13.00322.52315.91515.90315.919–+10.959–−6.608
–
+.016

104.2
9243.544.4
874.0
4797.94960.7–+9139.3–+829.6–+162.8#: Change from uncorrected surface; IC: isometric correction; LB: Laplace–Beltrami; SI1: surface inflation method #1; SI2: surface inflation method #2.

Yotte et al: Algorithms to Improve Reparameterization of Spherical Mappings

5

Fig 3.

The average metric, area, and angular distortion histograms are shown here for 20 central surface meshes. The isometry correction
process improves most profiles, except for the angular distortion for the Laplace-Beltrami initial mapping. LB: Laplace–Beltrami; SI1: surface
inflation method #1; SI2: surface inflation method #2.

that the surfaces have regularly sampled points with respect toθ and φ, where θ is the colatitude and φ is the azimuthal coordinate. To accomplish this, points are generated from equally
sampled values of θ and φ for all members in the sets, such that
there are a varying number of points per set (128, 256, 512, or
1,024). For each regularly sampled spherical point, the bounding triangle on the cortical spherical map is found. Within thatbounding triangle, a coordinate of location is approximated using barycentric coordinates, for example, the location of theregularly sampled point within the bounding triangle on thespherical map determines a certain set of barycentric coordinates, and these barycentric coordinates are then used within
the corresponding triangle in the original cortical surface tofind the cortical location. The result is a regularly sampled
spherical map in which every point is associated with a coordinate related to the location on the original cortical surface (Fig 2). Near the poles (located along the north–southaxis), the mesh has a higher resolution than at the equator.
This increased point density at the poles may be avoided byregridding using refinement of Platonic solids; however, thisadditional step was not implemented for the results shown
later.47
The reparameterized surfaces were then compared to theoriginal brain surface mesh using the mean distance error and$the Hausdorff distance. The minimum distance function d(p,S )between a point p∈X on the reparameterized mesh and the6

Journal of Neuroimaging Vol XX No X XX 2010

original brain surface S $ can be defined as++d ( p, S$ ) = min + p − S$ + .

(14)

The mean distance error de is the average minimum distance
between a set of points X and a surface S $
de =

1 !
d( p, S$ ),
Np

(15)

p∈X

where Np is the number of points in the set of points X .
The Hausdorff distance is simply the maximum distance error within the set of minimum distances:d H = max{d( p, S)},

(16)

There are two surfaces that are being compared, S and S $ . It
is key to note that the minimum distance from the set of points
on surface S to the other surface S $ is not equivalent to theminimum distance from the set of points on surface S $ to surface
S . When measuring the distance error from the reparameterized
surface to the original brain surface, this will be referred to as
the forward Hausdorff distance and mean distance error; in theother direction, it is the reverse Hausdorff distance and mean
distance error.

Fig 4. Shown is the distribution of the metric distortion across a sample central surface mesh. The isometry correction process (rightcolumn) reduces metric distortion for all initial spherical mapping methods (left column). For SI2, only the fine correction was used. P: isometry
correction processed; LB: Laplace–Beltrami; SI1: surface inflation method #1; SI2: surface inflation method #2.

Sample Data
All results were generated using 10 healthy control subject
brains. For the 10 control subjects (4 females and 6 males,mean age 35.3 years, SD ±11.3), T1-weighted images wereobtained on a 1.5T Philips Gyroscan ACSII. There were 256
sagittal slices per scan (1 mm thickness, T R = 13 ms, T E =5 ms, flip angle = 25◦ , field-of-view = 256 mm) with a matrix size of 256 × 256, resulting in an isotropic voxel sizeof 1 mm3 . All images were preprocessed using the processing pipeline (recon-all-all) included in FreeSurfer v4.0.5(http://surfer.nmr.mgh.harvard.edu/) to obtain the initial whitematter (WM) and pial surface meshes. Because there is a direct correspondence between the vertices in these two surfaces, these surfaces were averaged to obtain the central
surface.
The initial spherical mappings were generated using one ofthree previously published methods: solving the LB operator
with a Möbius transformation optimization19 ; inflating the brain
surfaces using a smoothing operation (SI1)48 ; and inflating thebrain surfaces while optimizing for metric distortion (SI2).39
The isometry correction process was applied to both the WM
and the central surface spherical maps for both the right andleft hemispheres.

For all methods, the time to calculate the resulting spherical
map mesh was obtained by running all analyses on an Apple
Xserve with 2 Quad-Core Intel Xeon with 2.6 GHz running
Mac OS X 10.5. No other programs were run at the same time
while running the analysis.

Results and Discussion
Distortion Metrics Improve via Algorithmic Correction
The isometry correction process based on the three distortion correction algorithms outlined above creates a pseudoisometric map of the cortical surface. Applied to the three
initial spherical maps, the result is a spherical map that has
significantly lower area and local metric distortion, with varying results on angular distortion and minimal cost in regards
to computational time. There is low variability in the distortion
results across subjects, indicating that the isometry correction
process is fairly robust independent of the input spherical map.
Computation times depend on both the mesh size (which forthe sample data is an average of 126,000 points per mesh) andthe average area distortion of the input spherical map.
The average area, angular, and local metric distortion values, as well as the real computation time, for the central surface

Yotte et al: Algorithms to Improve Reparameterization of Spherical Mappings

7

Fig 5. Shown is the distribution of the area distortion across a sample central surface mesh. The isometry correction process (right column)
reduces area distortion for all initial spherical mapping methods (left column). For SI2, only the fine correction was used. P: isometry correction
processed; LB: Laplace–Beltrami; SI1: surface inflation method #1; SI2: surface inflation method #2.
and WM surface are given in Tables 2 and 3, respectively. The
initial spherical map that has fairly low metric distortion is the
one generated using the SI2 method. In this case, the isometry correction process makes only minimal modification to thestarting map, especially for the WM surface. However, for boththe central and WM surfaces, the area and metric distortions are
reduced. Conversely, the large initial area distortions of the LB
conformal maps result in a significant increase in computation
time compared to the other two initial maps.
The change in the distortion distribution due to the isometry correction process for the central surface is shown in
Figure 3. In most cases, the isometry correction process results in an improved distortion profile. However, the LB conformal map had initially low angular distortion that was lost
during the isometry correction process. The histogram profile
is highly consistent for the mesh passed through the isometry
correction process, independent of the initial spherical mapping
method.
Figures 4–6 show the distribution of the metric, area, and angular distortion of the spherical map, mapped across the original
central surface mesh. For the LB and SI1 mapping processes,
both the coarse and fine corrections were used; the SI2 map
was only processed with the fine correction. Under ideal circumstances, the distortion would be evenly distributed across

8

Journal of Neuroimaging Vol XX No X XX 2010

the brain surface; however, it is not. With respect to angular distortion, the LB conformal map has very low angular distortion,but other methods create a map that has relatively large angulardistortion along the gyri and not the sulci. The LB conformal
map also compresses triangles at the rostral and caudal poles,
as reflected in the large negative area distortion values locatedin these regions.
Most notably, the distortion in the final spherical map appears to be independent of the starting spherical map, as reflected numerically in the distortion metric values and visuallyin the histograms. Although there are an infinite number of solutions for mapping points onto a sphere, the minimum-distortion
solution is a small subset of possible solutions. The similarity between the final mappings indicates that the isometry correction
process converges toward a common minimum. The rate ofconvergence is initially fast, and then decreases as the sphericalmapping solution approaches convergence.
We must also confirm that the cortical mapping onto thesphere does not have any large distortions. As Figure 7 shows,the isometry correction process dramatically reduces the distortion evident in the LB and SI1 mapping, while there is no
significant change from the SI2 mapping. The LB conformal
map devotes significantly more surface area to gyral regions
than was originally present in the brain surface.

Fig 6. Shown is the distribution of the angular distortion across a sample central surface mesh. The isometry correction process (rightcolumn) reduces angular distortion for all initial spherical mapping methods (left column) except for LB. For SI2, only the fine correction wasused. P: isometry correction processed; LB: Laplace–Beltrami; SI1: surface inflation method #1; SI2: surface inflation method #2.
Although the maps were of the central surface of one brain,
the results generalize to the other central surface and WM maps.
The isometry correction process generates a spherical map withlow overall distortion, independent of the initial spherical mapping method.

Effects of Disparate Algorithms on Area and Angle Distortion
The comprehensive effects of each algorithm on the sphericalmapping cannot be fully understood due to the complex natureof the mesh being acted upon. However, it is worthwhile toattempt to tease out the peculiar effects of each algorithm on thespherical mappings by passing a map through 1,000 iterationsof each algorithm alone and measuring the resulting area andangular distortion (Fig 8). The LB conformal map was chosen as
the initial map due to its relatively high initial metric distortion.
The ray-tracing algorithm (Eq 13) and its large-triangle-only
variant produced maps that had approximately the same area
distortion but much larger angular distortion compared to theoriginal LB conformal map. An examination of the maps reveals
that this is due precisely to the complex nature of the meshbeing acted upon. Some triangles are corrected at a high costto neighboring triangles. The ray-tracing algorithm was indeed

found to perform well only when the area and angular distortion
were already fairly low.
The two weighted algorithms (Eqs 8 and 11) generated mapsthat had significantly lower area distortion than the original LB
conformal map, with the distortion variant (Eq 11) outperforming the area-only algorithm (Eq 8) in both metrics. In general,the distortion variant (Eq 11) was found to be the best solution
for generating a spherical map from the LB conformal map interms of computation speed and distortion metrics, comparedto other algorithms.

Reparameterization Accuracy as a Function of Distortion
and Resolution
Generally, the forward Hausdorff distance and mean distanceerror metrics measure how well the reconstructed surface
matches the original surface mesh; the reverse Hausdorff distance and mean distance error indicate how much informationfrom the original surface mesh was retained in the reconstruction. It is therefore expected that the forward metrics may increase as a function of resolution, since there are simply morepoints in the mesh that will have potentially larger distanceerrors. However, these measures should be extremely smallvalues, which was experimentally confirmed; the two metrics

Yotte et al: Algorithms to Improve Reparameterization of Spherical Mappings

9Fig 7. Cortical mappings onto the spherical mesh are shown on the sphere for the lateral (left) and medial (right) cortical central surfaces
for one subject. The top row shows the inflated surface as a reference. The surface achieved via the isometry correction process is similar,
independent of the initial spherical mapping. For the SI2 method, there is almost no change, since the original mapping had low distortion. IC:
isometry correction; LB: Laplace–Beltrami; SI1: surface inflation method #1; SI2: surface inflation method #2.
are on the order of nanometers rather than millimeters (Figs 9
and 10).
For the reverse metrics, it is expected that a higher resolution
mesh will decrease these values, which was found to be the case
(Figs 9 and 10). The mean distance error can be considered to bea more critical measure than the Hausdorff distance, because
the Hausdorff distance is the absolute maximum rather thanan overall measure for the entire surface mesh. A large mean
distance error indicates that geometric information from theoriginal cortical surface mesh is not adequately represented in
the reparameterized mesh.
For both the central and WM surfaces, the reverse mean
distance error decreased if the isometry correction process wasapplied; for the SI2 method, there was only a slight decreasein this measure (Tables S1 and S2). The reverse Hausdorff distance also decreased for all cases, with the one exception of
SI2 for the CS at low resolutions. It can be assumed that thisoccurred due to the way that the SI2 mesh was created. The
surface was inflated while minimizing a squared edge length
error. Because the value is squared, outliers are punished moreso than in the isometry correction process, which uses an error
function that is directly proportional to the length error. Thiscan be translated into fewer distance error outliers during the10

Journal of Neuroimaging Vol XX No X XX 2010

Fig 8. A scatter plot of average area distortion metrics versus
average angular distortion metrics is shown here, for central surface
spherical maps using just a single algorithm over 1,000 iterations.
The average area and angular distortion metrics (with standard errorof the mean) are included. LB: Laplace–Beltrami; A: area; D: distort;
R: ray-trace; LR: ray-trace for large triangles only.

Fig 9. For the central surface, the average Hausdorff distance and mean distance error change as a function of mesh resolution. For

the reverse Hausdorff distance, the PLB and PSI2 methods are almost equivalent; for the mean distance error, the PLB, PSI1, SI2, and
PSI2 methods overlap for most mesh resolutions. Error bars show the standard error of the mean for the 20 meshes. P: isometry correction
processed; LB: Laplace–Beltrami; SI1: surface inflation method #1; SI2: surface inflation method #2.

reparameterization process, which will have an effect on thereverse Hausdorff distance (a maximization term) but not on
the reverse mean distance error. Instead, the mean distanceerror is minimized if the distortion is minimized, which was
found to be the case, because the isometry correction process reduces both the distortion metrics and the reverse mean distanceerrors.

Conclusion
Our isometry correction process significantly decreases the distortion of spherical maps using a combination of three originalalgorithms. The starting spherical map was generated using oneof three common methods: a conformal map obtained by solving a partial differential equation involving the LB operator of

the surface coordinates and optimized with a Möbius transformation, iteratively smoothing and inflating the brain surface
mesh until a spherical shape is obtained, and inflating the brain
surface mesh while minimizing a metric distortion energy functional. Independent of the method used to generate the starting
spherical map, the map that resulted after processing with theisometry correction algorithms exhibited lower metric and area
distortion, at relatively minimal cost to angular distortion.
The analyses were conducted on the central and WM surfaces. Conceivably, there are three potential surfaces that canbe analyzed: the WM surface, the pial surface, and the central
surface. The central surface can be constructed by averaging
the pial and WM surfaces, or by extracting it directly fromthe volumetric data. Generally, the pial surface tends to havemore pronounced folds, whereas the WM surface tends to have

Yotte et al: Algorithms to Improve Reparameterization of Spherical Mappings

11

Fig 10. For the WM surface, the average Hausdorff distance and mean distance error change as a function of mesh resolution. For thereverse mean distance error, the PLB, PSI1, SI2, and PSI2 methods overlap for most mesh resolutions. Error bars show the standard error of
the mean for the 20 meshes. P: isometry correction processed; LB: Laplace–Beltrami; SI1: surface inflation method #1; SI2: surface inflation
method #2.

higher average curvature, or a rougher surface. We chose tofocus on the central surface to have a midway representation
between these two extremes. The central surface might represent the cortical ribbon more accurately than surfaces along
either the inner or outer boundary,15 because both sulci andgyri are better defined. However, to confirm the generality of
the isometry correction process, the numerical results for theWM surfaces were also included. These results demonstrate amore significant correction via the isometry correction process
than was observed for the central surface.
Each algorithm used in the processing pipeline has an antecedent, and elucidating these relationships may help future
research. For instance, the weighted algorithm that updates thevertices based solely on area is related to smoothing functions,because by default it creates a distribution of triangle areas that

12

Journal of Neuroimaging Vol XX No X XX 2010

is more uniform. An improved variant of this algorithm would
be directly related to area distortion, and it was experimentally confirmed that choosing a weight based on area distortion
results in an improved outcome.
A class of algorithms that could be considered antecedent
to the ray-tracing algorithm is those that calculate energy functions that minimize metric distortion.39,49 However, metric distortion is an indefinite term. Some previous definitions include
the variance of the distance error extremes, the square-mean
error of distance, or the direct distance error. Furthermore,
the distance error can be local or global, in which the latterincorporates the geodesic distance error of each vertex to allother vertices in the mesh. It may be of interest to elucidate
one algorithm from this class in detail. The FreeSurfer suite reduces distortion by calculating an energy functional that must beminimized. The energy functional includes a value representing
a form of local metric distortion as well as a spring force to drive
surface inflation, where the metric distortion term is the meansquared local distance error. Minimizing this energy functional
during surface inflation produces a surface with low distortion
but at relatively high computational cost. In contrast, the raytracing algorithm is applied iteratively and potentially independently of the mesh distortion, thus reducing the computation
time.
Certainly, it is optimal to reduce the global metric distortion, which would minimize the change in geodesic distance
between prominent surface features.50 The computational costof computing the global metric distortion naturally restricts the
direct use of this metric when inflating the surface or mapping
to a sphere. Theoretically, the local overall reduction of metric distortion inherently leads to low global metric distortion,
yet forthcoming developments in computational power may
permit a more direct application of global metric distortion in
spherical mapping.
Reparameterization of the spherical mapping was improved
with respect to distance error values if the area distortion of thespherical mapping was reduced. A related problem in spline
theory optimizes this criterion (eg, distance error) directly byadjusting the location of the reparameterized points. If the goal
is to maximize the surface fidelity using a fixed number ofpoints, the best mapping will be one that is adaptive to the
local curvature of the surface, for example, flatter areas of thesurface will be sampled using fewer points, and the point density
will increase in high-curvature regions. An equiareal mapping
achieves this approximately, which is reflected in the decreased
distance errors during re-parameterization.
This work was supported by the following grants from the German
Bundesministerium für Bildung und Forschung: BMBF 01EV0709 and
BMBF 01GW0740.
Funding sources: BMBF 01EV0709, BMBF 01GW0740.

References
1. Memoli F, Sapiro G, Thompson P. Implicit brain imaging. NeuroImage 2004;23(Suppl 1):S179-S188.
2. Thompson PM, MacDonald D, Mega MS, et al. Detection andmapping of abnormal brain structure with a probabilistic atlas ofcortical surfaces. J Comput Assist Tomogr 1997;21(4):567-581.
3. Thompson PM, Schwartz C, Toga AW. High-resolution randommesh algorithms for creating a probabilistic 3D surface atlas of thehuman brain. NeuroImage 1996;3:19-34.
4. Narr KL, Bilder RM, Kim S, et al. Abnormal gyral complexity infirst-episode schizophrenia. Biol Psychiatry 2004;55(8):859-867.
5. Luders E, Narr KL, Thompson PM, et al. Gender differences incortical complexity. Nat Neurosci 2004;7(8):799-800.
6. Thompson PM, Schwartz C, Lin RT, et al. Three-dimensional statistical analysis of sulcal variability in the human brain. J Neurosci
1996;16(13):4261-4274.
7. Desai R, Liebenthal E, Possing ET, et al. Volumetric vs. surfacebased alignment for localization of auditory cortex activation. NeuroImage 2005;26(4):1019-1029.
8. Hinds OP, Rajendran N, Polimeni JR, et al. Accurate predictionof V1 location from cortical folds in a surface coordinate system.
NeuroImage 2008;39(4):1585-1599.

9. Fischl B, Sereno MI, Tootell RBH, et al. High-resolution intersubject averaging and a coordinate system for the cortical surface. Hum
Brain Mapp 1999;8(4):272-284.
10. Thompson PM, Toga AW. A surface-based technique for warping
three-dimensional images of the brain. IEEE Trans Med Imaging
1996;15(4):402-417.
11. Luders E, Thompson PM, Narr KL, et al. A curvature-based approach to estimate local gyrification on the cortical surface. NeuroImage 2006;29(4):1224-1230.
12. Schaer M, Cuadra MB, Tamarit L, et al. A surface-based approach to quantify local cortical gyrification. IEEE Trans Med Imaging 2008;27(2):161-170.
13. Han X, Jovicich J, Salat D, et al. Reliability of MRI-derived
measurements of human cerebral cortical thickness: the effectsof field strength, scanner upgrade and manufacturer. NeuroImage
2006;32(1):180-194.
14. Fischl B, Dale AM. Measuring the thickness of the human cerebral
cortex from magnetic resonance images. Proc Natl Acad Sci USA
2000;97(20):11050-11055.
15. Van Essen DC, Drury HA, Dickson J, et al. An integrated software
suite for surface-based analyses of cerebral cortex. J Am Med Inform
Assoc 2001;8(5):443-459.
16. Floater MS, Hormann K. Surface parameterization: a tutorial andsurvey. In: Dodgson N, Floater MS, Sabin M, eds. Advances inMultiresolution for Geometric Modelling . Vol IV. Berlin Heidelberg:
Springer; 2005:157-186.
17. Schwartz EL. Spatial mapping in the primate sensory projection: analytic structure and relevance to perception. Biol Cybern
1977;25(4):181-194.
18. Angenent S, Haker S, Tannenbaum A, et al. On the LaplaceBeltrami operator and brain surface flattening. IEEE Trans Med
Imaging 1999;18(8):700-711.
19. Tosun D, Rettmann ME, Prince JL. Mapping techniques for aligning sulci across multiple brains. Med Image Anal 2004;8(3):295-309.
20. Lui LM, Wang Y, Chan TF, et al. Automatic landmark tracking
applied to optimize brain conformal mapping. Proceedings of theIEEE International Symposium on Biomed Imaging . Arlington, VA,
2006.
21. Lui LM, Wang YL, Chan TF, et al. Brain anatomical feature detection by solving partial differential equations on a general manifold.
J Discrete Cont Dyn Syst, Ser B 2007;7(3):605-618.
22. Lui LM, Wang Y, Chan TF, et al. Landmark constrained genus zerosurface conformal mapping and its application to brain mapping
research. Appl Numer Math 2007;57(5-7):847-858.
23. Gu X, Wang Y, Chan TF, et al. Genus zero surface conformal
mapping and its application to brain surface mapping. IEEE TransMed Imaging 2004;23(8):949-958.
24. Wang Y, Chiang M-C, Thompson PM. 3D surface matching withmutual information and Riemann surface structures. Comput Graph
Imaging . Honolulu, HI; 2005:94-99.
25. Joshi AA, Shattuck DW, Thompson PM, et al. Cortical surface
parameterization by p-harmonic energy minimization. Proceedingsof the IEEE International Symposium of Biomed Imaging . pp. 428-431.
Arlington, VA, 2004.
26. Joshi AA, Shattuck DW, Thompson PM, et al. Surface-constrained
volumetric brain registration using harmonic mappings. IEEE TransMed Imaging 2007;26(12):1657-1669.
27. Hurdal MK, Stephenson K. Cortical cartography using the discrete
conformal approach of circle packings. NeuroImage 2004;23(Suppl
1):S119-S128.
28. Wang Y, Gu X, Chan TF, et al. Multivariate tensor-based brain
anatomical surface morphometry via holomorphic one-forms. Med
Image Comput Comput Assist Interv 2009;5761:337-344.
29. Wang Y, Gu X, Chan TF, et al. Brain surface conformal parameterization with the Ricci flow. Proc IEEE Int Symp Biomed Imaging .
pp. 1312-1315. Arlington, VA, 2007.

Yotte et al: Algorithms to Improve Reparameterization of Spherical Mappings1330. Wang Y, Gu X, Chan TF, et al. Brain mapping with the Ricci flow
conformal parameterization and multivariate statistics on deformation tensors. Proceedings of the 2nd MICCAI Workshop on Mathematical
Foundations of Computational Anatomy. pp. 36-47. New York, NY,
2008.
31. Wang Y, Gu X, Chan TF, et al. Brain surface conformal parameterization with algebraic functions. Med Image Comput Comput Assist
Interv 2006;9(Pt 2):946-954.
32. Wang Y, Gu X, Chan TF, et al. Brain surface conformal parameterization with the slit mapping. Proceedings of the IEEE International
Symposium on Biomed Imaging . pp. 448-451. Paris, 2008.
33. Wang Y, Gu X, Chan TF, et al. Shape analysis with conformal
invariants for multiply connected domains and its application to
analyzing brain morphology. NeuroImage 2009;47(Suppl 1):S39S41.
34. Thompson PM, Hayashi KM, Sowell ER, et al. Mapping corticalchange in Alzheimer’s disease, brain development, and schizophrenia. NeuroImage 2004;23(Suppl 1):S2-S18.
35. Leow A, Yu CL, Lee SJ, et al. Brain structural mapping using anovel hybrid implicit/explicit framework based on the level-set
method. NeuroImage 2005;24(3):910-927.
36. Wang Y, Lui LM, Chan TF, et al. Optimization of brain conformal
mapping with landmarks. Med Image Comput Comput Assist Interv
2005;8(Pt 2):675-683.
37. Shi Y, Thompson PM, Dinov I, et al. Direct cortical mapping via
solving partial differential equations on implicit surfaces. Med Image
Anal 2007;11(3):207-223.
38. Thompson PM, Hayashi KM, de Zubicaray GI, et al. Detecting dynamic and genetic effects on brain structure using high-dimensional
cortical pattern matching. Proceedings of the IEEE International Symposium on Biomed Imaging . pp. 473-476. Washington, DC, 2002.
39. Fischl B, Sereno MI, Dale AM. Cortical surface-based analysis. II.
Inflation, flattening, and a surface-based coordinate system. NeuroImage 1999;9(2):195-207.
40. Carman GJ, Drury HA, Van Essen DC. Computational methods
for reconstructing and unfolding the cerebral cortex. Cereb Cortex
1995;5(6):506-517.
41. Kruggel F. Robust parametrization of brain surface meshes. Med
Image Anal 2008;12(3):291-299.
42. Shen L, Makedon F. Spherical parameterization for 3D surface
analysis in volumetric images. Proceedings of the International Conference on Information Technology: Coding and Computing . Vol 2, pp.
643-649. Las Vegas, NV, 2004.

14

Journal of Neuroimaging Vol XX No X XX 2010

43. Timsari B, Leahy RM. Optimization method for creating semiisometric flat maps of the cerebral cortex. SPIE Symposium on Medical
Imaging 2000: Image Processing . Vol. 3979, pp. 698-708. San Diego,CA, 2000.
44. Liseikin VD. Grid Generation Methods. Heidelberg: Springer-Verlag;
1999
45. Tinoco-Ruiz J-G, Barrera-Sanchez P, Cortes-Medina A. Some
properties of area functionals in numerical grid generation. Proceedings of the 10th International Meshing Roundtable. Newport, pp.
43-54. Beach, CA, 2001.
46. Ju L, Stern J, Rehm K, et al. Cortical surface flattening using leastsquare conformal mapping with minimal metric distortion. Proceedings of the IEEE International Symposium on Biomed Imaging . pp. 77-80.
Arlington, VA, 2004.
47. MacDonald D. A method for identifying geometrically simple surfaces
from three dimensional images. Montreal: School of Computer Science, McGill Unversity, 1998.
48. Drury HA, Van Essen DC, Anderson CH, et al. Computerized mappings of the cerebral cortex: a multiresolution flattening method and a surface-based coordinate system. J Cogn Neurosci
1996;8(1):1-28.
49. Praun E, Hoppe H. Spherical parametrization and remeshing. ACM
Trans Graph 2003;22(3):340-349.
50. Schwartz EL, Shaw A, Wolfson E. A numerical solution to the generalized mapmaker’s problem: flattening nonconvex polyhedralsurfaces. IEEE Trans Pattern Anal Machine Intell 1989;11(9):10051008.

Supporting Information
Additional supporting information may be found in the online
version of this article:
Table S1. Mean distance error and Hausdorff distances for 20
central surface meshes
Table S2. Mean distance error and Hausdorff distances for 20
white matter surface meshes
Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting materials supplied by theauthors. Any queries (other than missing material) should be directed to the corresponding author for the article.