YNIMG-09817; No. of pages: 14; 4C:
NeuroImage xxx (2012) xxx–xxx

Contents lists available at SciVerse ScienceDirect

NeuroImage
journal homepage: www.elsevier.com/locate/ynimg

Q1 2
3

Cortical thickness and central surface estimation
Robert Dahnke ⁎, Rachel Aine Yotter 1, Christian Gaser 1

F

1

Department of Psychiatry, University of Jena, Jahnstrasse 3, D-07743 Jena, Germany

i n f oa b s t r a c t

RO

a r t i c l e

Several properties of the human brain cortex, e.g., cortical thickness and gyriﬁcation, have been found to cor- 20
relate with the progress of neuropsychiatric disorders. The relationship between brain structure and function 21
harbors a broad range of potential uses, particularly in clinical contexts, provided that robust methods for the 22
extraction of suitable representations of the brain cortex from neuroimaging data are available. One such rep- 23
resentation is the computationally deﬁned central surface (CS) of the brain cortex. Previous approaches to 24
semi-automated reconstruction of this surface relied on image segmentation procedures that required man- 25
ual interaction, thereby rendering them error-prone and complicating the analysis of brains that were not 26
from healthy human adults. Validation of these approaches and thickness measures is often done only for 27
simple artiﬁcial phantoms that cover just a few standard cases. Here, we present a new fully automated 28
method that allows for measurement of cortical thickness and reconstructions of the CS in one step. It uses 29
a tissue segmentation to estimate the WM distance, then projects the local maxima (which is equal to the 30
cortical thickness) to other GM voxels by using a neighbor relationship described by the WM distance. This 31
projection-based thickness (PBT) allows the handling of partial volume information, sulcal blurring, and 32
sulcal asymmetries without explicit sulcus reconstruction via skeleton or thinning methods. Furthermore, 33
we introduce a validation framework using spherical and brain phantoms that conﬁrms accurate CS construc- 34
tion and cortical thickness measurement under a wide set of parameters for several thickness levels. The re- 35
sults indicate that both the quality and computational cost of our method are comparable, and may be 36
superior in certain respects, to existing approaches.
37© 2012 Published by Elsevier Inc. 38

Article history:
Accepted 20 September 2012
Available online xxxx

E

D

P

Keywords:
MRI
Cortical thickness
Central surface
Surface reconstruction
Brain
Phantom
Validation

E

C

T

5
6
7
8
10
9
11
12
13
14
15
16
17
18
19

O

4

42

R

41

Introduction

44

The cerebral cortex is a highly folded sheet of gray matter (GM) thatlies inside the cerebrospinal ﬂuid (CSF) and surrounds a core of whitematter (WM). Besides the separation into two hemispheres, the cortex
is macroscopically structured into outwardly folded gyri and inwardly
folded sulci (Fig. 1). The cortex can be described by the outer surface(or boundary) between GM and CSF, the inner surface (or boundary)
between GM and WM, and the central surface (CS) (Fig. 1). Corticalstructure and thickness were found to be an important biomarker for
normal development and aging (Fjell et al., 2006; Sowell et al., 2004,
2007) and pathological changes (Kuperberg et al., 2003; Rosas et al.,2008; Sailer et al., 2003; Thompson et al., 2004) in not only humans,
but also other mammals (Hofman, 1989; Zhang and Sejnowski, 2000).
Although MR images allow in vivo measurements of the humanbrain, data is often limited by its sampling resolution that is usuallyaround 1 mm3. At this resolution, the CSF is often hard to detect insulcal areas due to the partial volume effect (PVE). The PVE comes
into effect for voxels that contain more than one tissue type and have

49
50
51
52
53
54
55
56
57
58
59
60

N
C
O

47
48

U

45
46

R

43

⁎ Corresponding author. Fax: +49 3641 934755.
E-mail addresses: robert.dahnke@uni-jena.de (R. Dahnke),rachel.yotter@uni-jena.de (R.A. Yotter), christian.gaser@uni-jena.de (C. Gaser).
1
Fax: +49 3641 934755.

40
39

an intensity gradient that lies somewhere between that of the pure tis- 61
sue classes. Normally, the PVE describes the boundary with a sub-voxel 62
accuracy, but within a sulcus the CSF volume is small and affected by 63
noise, rendering it difﬁcult to describe the outer boundary in this region 64(blurred sulcus, Fig. 2). Thus, to obtain an accurate thickness measure- 65
ment, an explicit reconstruction of the outer boundary based on the 66
inner boundary is necessary. This can be done by skeleton (or thinning) 67
methods or alternatively by model-based deformation of the inner sur- 68
face. Skeleton-based reconstruction of the outer boundary is used by 69
CLASP (Kim et al., 2005; Lee et al., 2006a,2006b; Lerch and Evans, 702005), CRUISE (Han et al., 2004; Tosun et al., 2004; Xu et al., 1999), 71
Caret (Van Essen et al., 2001), the Laplacian approach (Acosta et al., 72
2009; Haidar and Soul, 2006; Hutton et al., 2008; Jones et al., 2000; 73
Rocha et al., 2007; Yezzi and Prince, 2003), and other volumetric 74
methods (Eskildsen and Ostergaard, 2006, 2007; Hutton et al., 2008; 75
Lohmann et al., 2003). Methods without sulcal modeling will tend to 76
overestimate thickness in blurred regions (Jones et al., 2000; 77
Lohmann et al., 2003) or must concentrate exclusively on non-blurred 78
gyral regions (Sowell et al., 2004). Alternatively, cortical thickness 79
may be estimated via deformation of the inner surface (FreeSurfer 80(Dale et al., 1999; Fischl and Dale, 2000), DiReCT (Das et al., 2009), 81
Brainvoyager (Kriegeskorte and Goebel, 2001), Brainsuite (Shattuck 82
and Leahy, 2001; Zeng et al., 1999) or coupled surfaces (ASP (Kabani 83
et al., 2001; MacDonald et al., 2000). Considering that the accuracy of 841053-8119/$ – see front matter © 2012 Published by Elsevier Inc.
http://dx.doi.org/10.1016/j.neuroimage.2012.09.050

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

F

2

E

D

P

achieved by FreeSurfer and a signiﬁcant improvement over the 96
Laplacian approach, it may be concluded that PBT is a highly accurate 97
volume-based method for measuring cortical thickness. For situations 98
in which extensive surface analysis is not required, PBT would allow 99
the exclusion of cortical surface reconstruction steps with no loss of 100
accuracy for cortical thickness measurements.
101
We also propose a suite of test cases using a variety of phantoms 102
with different parameters as a suggestion for how a cortical thickness 103
measurement approach could be rigorously tested for validity and 104
stability. Previously published validation approaches that used a 105
spherical phantom (Acosta et al., 2009; Das et al., 2009) often 106
addressed only one thickness and curvature (radii) of the inner and 107
outer boundary. The problem is that the measure may work well for 108
this special combination of parameters, but performance can change 109
for different radii. Another limitation is that this phantom describes 110
only areas where the CSF intensity is high enough, but most sulcal 111
areas (that comprise over half of the human cortex) are blurred. 112
Our test suite directly addresses these concerns.
113
The cortical thickness map may also be subsequently used to gen- 114
erate a reconstruction of the CS. Compared to the inner or outer sur- 115
face, the CS allows a better representation of the cortical sheet (Van 116
Essen et al., 2001), since neither sulcal or gyral regions are over- or 117
underestimated (Scott and Thacker, 2005). As the average of two 118
boundaries, it is less error-prone to noise and it allows a better map- 119
ping of volumetric data (Liu et al., 2008; Van Essen et al., 2001). Gen- 120
erally, a surface reconstruction allows surface-based analysis that is 121
not restricted to the grid and allows metrics, such as the gyriﬁcation 122
index (Schaer et al., 2008) or other convolution measurements 123(Luders et al., 2006; Mietchen and Gaser, 2009; Rodriguez-Carranza 124
et al., 2008; Toro et al., 2008), that can only be measured using sur- 125
face meshes (Dale et al., 1999). It provides surface-based smoothing 126
that gives results superior to that obtained from volumetric smooth- 127
ing (Lerch and Evans, 2005). Furthermore, surface meshes allow a 128
better visualization of structural and functional data, especially 129
when they are inﬂated (Fischl et al., 1999) or ﬂattened (Van Essen 130
and Drury, 1997). Due to these considerations, we have explored 131
the quality of the cortical surface reconstructions.
132

T

C
E
R

94
95

R

92
93

O

90
91

C

88
89

the measurement depends strongly upon the precision of cortical surface
reconstruction at the inner and outer boundaries, and that the computation time is often related to the anatomical accuracy of the reconstruction,such measurements may require intensive computational resources inorder to achieve the ﬁnal measurement.
Here, we present a new volume-based algorithm, PBT (Projection
Based Thickness), that uses a projection scheme which considers
blurred sulci to create a correct cortical thickness map. For validation,
we compare PBT to the volumetric Laplacian approach and thesurface-based approach included in the FreeSurfer (v 4.5) software
package. If the results from PBT are approximately the same as that

N

86
87

U

85

R
O

O

Fig. 1. The cortex: Shown is an illustration of the cortical macro- and microstructure. The cerebral cortex is a highly folded sheet of gray matter (GM) that lies inside the cerebrospinal ﬂuid (CSF) and surrounds a core of white matter (WM). Inwardly folded regions are called sulci whereas outwardly folded areas are denoted as gyri. There are three common
surfaces to describe this sheet: the outer surface, the inner surface, and the central surface (CS). The CS allows a better representation of the cortical GM sheet and improved accuracy of cortical surface measurements. Cortical thickness describes the distance between the inner surface and the outer surface and is related to cortical development and diseases such as Alzheimer's.

Fig. 2. Main ﬂow diagram: Shown is a ﬂow diagram of the pre-processing steps of the
CS and thickness estimation. A tissue segmentation algorithm (from VBM8) is used tocreate a segmentation image SEG from an anatomical image. This segmentation image
is used for (manual) separation of the cortex into two hemispheres and removal of thecerebellum with hindbrain, resulting in a map SEP. This map creates the map SEGPF’, amasked version of SEG with ﬁlled ventricular and subcortical regions. Both approaches
used an interpolated version of the map SEGPF’ to create a CS with a cortical thickness
value of each vertex. The red subﬁgure shows blurred sulcal regions, where CSF voxels
were detected as GM due to noise removal included in the segmentation algorithm.
These blurred regions need an explicit reconstruction of the outer surface for the
Laplacian approach (Fig. 4), whereas PBT uses an inherent scheme to account forthese regions (Fig. 3). (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)

Material and methods

133

We start with a short overview about the main steps of our meth- 134
od and the Laplacian approach; algorithmic details are separately de- 135
scribed in the following subchapters.
136

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

157
158
159

Q3 160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180

185
186
187
188
189
190
191

192
193
194

C

Q2 156E2

http://dbm.neuro.uni-jena.de/vbm/
http://mouldy.bic.mni.mcgill.ca/brainweb/
4
http://dbm.neuro.uni-jena.de/vbm/5
http://www.ﬁl.ion.ucl.ac.uk/spm/
3

FDistance measure

209

O

154
155

R

152
153

R

150
151

N
C
O

148
149

U

146
147

To take into account the asymmetrical structures, we used the 210
Eikonal equation with a non-uniform speed function F(x) to ﬁnd the 211
closest boundary voxel B(x) of a GM voxel x without passing a differ- 212
ent boundary. To allow sub-voxel accuracy, the normalized vector be- 213
tween B(x) and x is used to ﬁnd a point G(x) between x and B(x). The 214
intensity gradient between B(x) and G(x) allows a precise estimation 215
of the boundary point P(x), which is used to estimate the distance of x 216
to the boundary.
217
In a more formal way, we solved the following Eikonal equation: 218

RO

To achieve exact and stable results for thickness measures, the segmentation plays an important role. In principle, any segmentation for
GM, WM, and CSF can be used. The segmentation could be binarymaps, but to achieve more stable and exact results, it is important touse probability maps that are able to describe the boundary positions
with sub-voxel accuracy (Hutton et al., 2008). Furthermore, inclusionof an additional noise removal step increases the accuracy and stabilityof the thickness measurements (Coupe et al., 2008). We used the VBM8
4
toolbox (revision 388) for SPM8 5(Ashburner and Friston, 2005) (revision 4290) for segmentation of all T1 images, which includes anoptimized Rician non-local mean (ORNLM) (Coupe et al., 2008) anda Gaussian Hidden Markov Random Field (GHMRF) (Cuadra et al.,

144
145

P

183
184

143

F ðxÞt∇DEi ðxÞt ¼ 1;

D

Segmentation

141
142

DEi ðxÞ ¼ 0;

f or

f or

x∈Ω;

x∈Γ;

ð1Þ

E

182

139
140

2005) ﬁlter for noise reduction (NR). The probability tissue maps 195
CSF, GM, and WM are combined in one probability image SEG 196(Tohka et al., 2004). Pure tissue voxels are coded with integers 197(background = 0, CSF = 1, GM = 2, WM = 3), whereas values be- 198
tween integers describe the percentile relation between the tissues. 199
For example, a voxel with an intensity of 2.56 contains 44% GM and 200
56% WM and a value of 1.92 contains 92% GM and 8% CSF. Hence, tis- 201
sue boundaries are at 0.5 between background and CSF, 1.5 between 202
CSF and GM and 2.5 between GM and WM. Note that this map is only 203
able to describe two tissue classes per voxel. However, this does not 204
degrade our analyses, because most anatomical images do not pro- 205
vide more information for the segmentation. Furthermore, most re- 206
gions with no GM layer, such as the brainstem or the near the 207
ventricles, are cut or ﬁlled and thus are not included in the analysis. 208

T

181

MRI images are ﬁrst segmented into different tissue classes using
VBM8 2(Fig. 2; see Segmentation). This segmentation is used for(manual) separation of the hemispheres and removal of the cerebellum with hindbrain, resulting in a map SEP. This map creates the map
SEGPF, a masked version of SEG with ﬁlled ventricular and subcorticalregions. To take into account the small sulci with thicknesses ofaround 1 mm, SEGPF was linearly interpolated to 0.5 × 0.5 × 0.5 mm 3
(Hutton et al., 2008; Jones et al., 2000).
For each GM voxel, the distance from the inner boundary was estimated within the GM using a voxel-based distance method (seeDistance measure). The result is a WM distance map WMD, whose
values at the outer GM boundary represent the GM thickness. These
values at the outer boundary were then projected back to the innerboundary, resulting in a GM thickness map GMT. The relation between
the WMD and GMT maps creates the percentage position map PP that isused to create the CS at the 50% level (see Projection-based thickness).
As a basis of comparison, we constructed another CS using the
Laplacian-based thickness measure (Jones et al., 2000) on the ﬁlled
tissue segmentation map to create another set of GMT and PP maps.
This method requires an explicit sulcal reconstruction step (Bouix
and Kaleem, 2000) (see Laplacian-based thickness).
A topology correction based on spherical harmonics was used tocorrect the topology of the surfaces generated with the PBT and theLaplacian approach (Yotter et al., in press).
For validation, a set of spherical (SP; see Spherical phantoms) andbrain phantoms (BP; see Brain phantoms) with uniform thickness
were used to simulate different curvature, thickness, noise, and resolution levels. Since thickness and the location of the cortical surfaces were
known, the two data sets could be directly compared. For thickness RMS
error, the measured thickness was reduced by the expected thickness.
In addition to the spherical phantoms with equal thickness, we used
the Collins brain phantom with different noise levels 3(Collins et al.,1998) and a real data set of 12 scans of the same subject of our database(see Real data) to compare our results to FreeSurfer 4.5. Because the real
thickness of both data sets is unknown, we compare the results of eachtested surface to the results of a surface that was generated on an averaged scan. RMS error was calculated for all vertices of a surface, including vertices of the ﬁlled subcortical regions and the corpus callosum. For
these data sets, we evaluated the number of topological errors using
Caret. To count the number of defects, the uncorrected CS was used
for PBT and Laplacian, whereas for FreeSurfer the uncorrected WM
surface was used. The CS of FreeSurfer was generated via Caret, where
the positions of CS vertices were given by the mean positions of corresponding vertices of the inner and outer surface. Thickness RMS error
was estimated based on the original FreeSurfer thickness results.

137
138

3

where x is a voxel, Ω is given by the GM, Γ is the object (the WM or 219
220
the CSF and background), DEi is the Eikonal distance map, and F(x) 221
is the non-uniform speed map (FWM(x) for the WM distance and 222
FCSF(x) for the CSF distance) that is given by the image intensity of 223
SEGPF:
224
F WM ðxÞ ¼ minð1; maxð0; SEGPF ðxÞ−1ÞÞ;F CSF ðxÞ ¼ minð1; maxð0; 3−SEGPF ðxÞÞÞ:

ð2Þ
226
225

In GM areas, FWM(x) has a high “speed” which results in shorter 227
distances, whereas in CSF areas the “speed” is very low and thus re- 228
sults in longer distances, whereas FCSF(x) allows high speeds in GM 229
and CSF areas, but not in WM regions. Because the distance map DEi 230
contains distortions, it is only used to ﬁnd the closest object voxel 231
for each GM voxel x ∈ Ω:232
BEi ðx; Ω; Γ; F Þ;

ð3Þand to calculate the Euclidean distance DEu between the GM voxel x 233
234
and its nearest WM voxel BEi(x, Ω, Γ, F):235
DEu ðx; Ω; Γ; F Þ ¼ tx; BEi ðx; Ω; Γ; F Þt2 :

ð4Þ
237
236

We solve the above equations as follows: By solving the Eikonal equa- 238
tion within Ω, we also note the closest WM voxel BEi. To allow sub-voxel 239 Q4
accuracy, the normalized vector between x and BEi(x, Ω, Γ, F) is used to es- 240
timate a point G(BEi(x, Ω, Γ, F)) within one voxel distance to BEi(x, Ω, Γ, F). 241
The intensity gradient between BEi(x, Ω, Γ, F) and G(BEi(x, Ω, Γ, F)) can 242
then be used to estimate the exact boundary of Γ.
243
Projection-based thickness244

For simpliﬁcation we will use the terms of the GM, WM, and CSF 245
probability maps for the operations, even though only the map 246
SEGPF is used. Cortical thickness can be described as the sum of the 247
inner (WMD, Fig. 3b2) and outer (CSFD, Fig. 3b3) boundary distance. 248

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

C

T

E

D

P

R
O

O

F

4

253
254
255
256
257
258
259
260
261
262
263
264
265
266

N

251
252

Blurring of the outer boundary in sulcal regions due to the PVE leads
to an overestimation of the CSFD. To avoid the explicit reconstructionof the outer boundary by a skeleton, we focus on the information
given by the WMD. At the outer boundary, and also within blurred regions, the GMT is fully described by the WMD, because the CSFD iszero (Lohmann et al., 2003; Sowell et al., 2004). In other words, thehighest local WMD within the GM is identical to the GMT of thisarea, and it is only necessary to project this information to otherGM voxels.
This can be done using the successor relationship of the WMD. A
neighbor voxel v2 of a voxel v1 is a successor of v1, if the WMD of v2
is around one voxel greater than the WMD of v1. Similarly, if the
WMD of v2 is around one voxel smaller than v1, v2 is labeled as theparent voxel. In this case, v1 gets the thickness value of v2. Neighborvoxels with a WMD similar to v1 that are too close to be either a parent or a successor are called siblings, and their thicknesses remain
unrelated to v1. If v1 has no successor, then it is a local maximumthat is located at the CSF boundary and its GMT is given by its WMD.

U

249
250

C

O

R

R

E

Fig. 3. PBT: Subﬁgure (a) shows a ﬂow diagram of the PBT approach, whereas subﬁgure (b1–b8 with simpliﬁed titles) shows 2D illustrations of the volume maps of (a). In subﬁgure(c), we illustrate the most relevant cases of our PBT method — a gyral and a blurred sulcal case with initialization and projection step. For distance calculations, the Eikonal equation
is solved to account for partial volume information. PBT starts with the (interpolated) masked segmentation image SEGPF shown in Fig. 2 and estimates the distance to the inner(b2) and outer (b3) boundary. The blurring of outer boundary in sulcal regions leads to strong overestimation of the real distance and ﬁnally to an overestimation of the cortical
thickness. To get the correct values in these regions, PBT uses a modiﬁed version GMTI (b4) of the WMD, in which the local maximum describes the position of the outer boundary
and the correct thickness. It now uses the successor relation succ(v) of a voxel v (Eq. (8)), given by the WM distance WMD (b2), to project thickness values from the outer boundary(b4) over the whole GM (b5). PBT additionally uses the direct GM thickness GMTD (b6) – which is overestimated in blurred areas, but helps to reduce artifacts such as blood vessels – tocreate a ﬁnal map GMTF (b7) of the minimum thickness from both thickness maps. After estimation of cortical thickness, a percentage position map PP is generated to create the CS andmap cortical thickness onto it. The projection scheme shown in subﬁgure (c) uses the WM distance map to project the maximum local WM distance that is equivalent to the local thickness
to other voxels. The WM distance map allows the deﬁnition of successors (neighbors of a voxel v with a slightly larger distance than v) and siblings (neighboring voxel with a similar
distance to v), and a voxel v gets the mean thickness of its successors. If a voxel has no successors, then it is located at the outer boundary and its WM distance is related to its size.

We now want to describe this process in a more formal way, 267
starting with the WMD:
268
WMDðvÞ ¼

DEu ðv; GM > 0; WM; F WM Þ ; if GM ðvÞ > 0;0; otherwise

ð5Þ

where DEu gives the Euclidean distance of a voxel v to the nearest WM 269
270
boundary that was found by solving the Eikonal equation for the 271
speed map FWM (Eq. (2)). The distance to the CSF boundary is now 272
given by:273
8< −DEu ðv; CSF & GM; CSF; & BG; 1Þ ; if GM ðvÞ > 0 & CSF ðvÞ > 0CSFDðvÞ ¼;
DEu ðv; GM > 0; CSF&BG; F CSF Þ; if GMðvÞ > 0:0; otherwise

ð6Þ
where BG (background) describes all voxels that contain no tissue. 274
275
The cortical thickness map GMTI is initialized as a modiﬁed version 276

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

279
280

PBT was used to reconstruct problematic regions in an additional 331
preprocessing step that estimates the cortical thickness in the GM 332
with ﬂipped boundaries. These problematic regions are those that 333
are highly susceptible to errors due to the PVE, which creates prob- 334
lems in both gyral and sulcal regions. In the gyral case, thin WM 335
structures are blurred rather than the CSF blurring that occurs in nar- 336
row sulcal regions. This occurs most frequently in the superior tem- 337
poral gyrus, the cingulate gyrus, and the insula, and may be 338
addressed similarly to the idea proposed in (Cardoso et al.) for seg- 339
mentation reﬁnement. If a voxel of the inverse thickness map has 340
lower thickness than the original thickness map and if the thickness 341
of both is larger than 2 mm while SSEG > 2.0, we expect that the in- 342verse thickness map has identiﬁed a gyrus that is blurred by the 343
PVE. For these blurred regions, the thickness and percentage position 344
of the inverse maps are used.
345

of the WMD, because the WMD describes the distance only to thecenter of a GM voxel. GM voxels with more than 50% CSF need additional correction by the CSFD, in which the minimum correction ishalf of the voxel resolution res:
GMT I ðvÞ ¼ WMDðvÞ þ minðCSFDðvÞ; res=2Þ:

ð7Þ

282
281
283
284
285

Let N26 be the 26-neighborhood of a voxel v, and D26 the associated distance of v to its neighbors. A voxel n ∈ N26(v) is a successor of
the voxel v if the WM distance of s meets the following conditional:ð Þðð ÞÞbWMDðnÞbðWMDðvÞ þ a2  D26 ðnÞÞ%succðv; nÞ¼ 1 ; if WMD v þ a1  D26 n; otherwise

0

297
298
299
300
301
302
303
304
305
306

pt ðvÞ ¼∑n∈N26 ðvÞ succðv; nÞGMT I ðnÞ∑n∈N26 ðvÞ succðv; nÞ:The initial thickness GMTI can now be used to estimate the ﬁnal
projection-based thickness map GMTP, by projecting the values overthe GM region:

313
312
314

Q5 315
316
317

N
C
O

GMT P ðvÞ ¼ maxðGMT I ðvÞ; pt ðvÞÞ;

R

311

322

326
327
328
329
330

U

ð12Þ

that is corrected for the voxel resolution res. The percentage position
map PP can now be described as:
PP ðvÞ ¼ ðGMT F ðvÞ−WMDðvÞ=resÞÞ=GMT F ðvÞ þ ðSEGPF ðvÞ≥2:5Þ:

323
324
325

ð11Þ

is used to create the ﬁnal thickness map:
GMT F ðvÞ ¼ minðGMT P ðvÞ; GMT D ðvÞÞ=res;

321
320

ð10Þ

This mapping can be done in O(n) time using the same principle
described for voxel-based distance calculation (Rosenfeld and John,
1966). To reduce overestimations in the GMTP map due to GM fragments such as blood vessels or dura mater, the direct thickness map:
GMT D ðvÞ ¼ CSFDðvÞ þ WMDðvÞ;

319
318

ð9Þ

R

308
307
309
310

R
O

P

295
296

ð13Þ

Finally, the CS is generated from the PP map and reduced toaround 300,000 nodes using standard Matlab functions. Each vertexof the mesh is assigned a thickness value via linear interpolation of
the closest GM thickness map values. Fig. 3 shows the ﬂow diagramof our method and illustrates the idea for most relevant examples in
2D.

The Laplacian approach requires an explicit sulcal reconstruction 347
step (Jones et al., 2000; Tosun et al., 2004) that uses a skeleton map 348
to reconstruct the outer boundary in blurred regions of the segment 349
image SEGPF (Fig. 4b2) by changing the tissue class of the re- 350
constructed boundary voxels from GM to CSF resulting in a map SEGPFS 351(Figs. 4b3, 4a). To create the skeleton map S, we ﬁrst generate WM and 352
CSF distance maps with the same distance measure used for PBT to 353
allow asymmetrical structures. We then ﬁnd areas with high divergence 354of the gradient ﬁeld, resulting in a map SR. This map is normalized with- 355
in a low and a high boundary slow = 0.5 and shigh = 1.0 resulting in the 356
skeleton map S (Bouix and Kaleem, 2000), with & as a logical AND op- 357 Q6
erator:
358

D

293
294

346

E

291
292

Laplacian-based thickness

T

289
290

where 0 b a1 ≤ 1 ≤ a2 b 2 are weights depending on the used distancemetric; these weights allow the inclusion of more thickness information from neighboring voxels to achieve a smoother thickness map. If
there are no successors, then v is a border voxel and the WM distancesets its thickness. The lower threshold a1 deﬁnes the boundary between siblings and successors, whereas the higher threshold a2 is alimit for direct successors. An a1 threshold that is too low will create
too many siblings and lead to smoother results, while an a1 thresholdthat is too high will lead to coarser results. Likewise, an a2 thresholdthat is too low will exclude more neighbors of v from the successor
relationship and lead to less smooth images and in the worst-caseto a breaking of the projection because all possible successors are excluded, whereas an a2 threshold that is too high will lead tooversmoothed results with overestimation in gyral regions. For aquasi-Euclidean metric, which is not useful for cortical thickness but
acceptable for the PP map, a1 and a2 are equal and can be set by the
distance of v to its neighbors. Good results with minimal smoothing
were achieved using a1 = 0.5 and a2 = 1.25. If there are no successors,then v is a border point and the WM distance sets its thickness, else it
uses the mean of all successors:

C

288

E

287
286

F

ð8Þ

O

277
278

5

SR ¼∇ΔðWMD
Þ
 




 
S ¼ SR  ðSR > slow Þ& SRbshigh −slow  shigh −slow þ SR≥shigh :
ð14Þ
360
359

The skeleton map accurately represents the sulci that have been 361
blurred in the tissue segmentation process. We correct all voxels of 362
SEG ≥ 1 by:
363
SEGPFS ¼ SEGPF − maxð1; 2−SÞðSEGPF ≥1ÞÞ

ð15Þ(Fig. 4b3). The changing of GM voxels to CSF voxels leads to an un- 364
365
derestimation of the GM volume and local thickness, which will be 366
considered later. The corrected segment map SEGPFS can now be 367
used to solve the Laplace equation between the GM/WM and GM/CSF 368
boundary:
369
2

∇ ψ¼

∂ψ
∂ψ ∂ψ
þ
þ
¼ 0:
∂x2 ∂y2 ∂z2

ð16Þ
371
370

The above equation is solved iteratively using an initial potential 372
image with Dirichlet boundary conditions. The WM (SEGPFS ≥ 2.5) 373
forms the higher potential boundary with values of 1, whereas the 374
CSF (SEGPFS ≤ 1.5) represents the lower potential boundary with 375
values of 0. To accelerate convergence, all GM voxels are initialized 376
with a potential of 0.5. Eq. (17) is applied only to GM voxels 377(SEGPFS > 1.5 and SEGPFS b 2.5) and simply describes the mean of the 378
six direct neighbors of a voxel:379
23
ψ ðx þ Δx; y; zÞ þ ψi ðx−Δx; y; zÞþ
1 4 i
ψi þ 1ðx; y; zÞ ¼  ψi ðx; y þ Δy; zÞ þ ψi ðx; y−Δy; zÞþ 5:6
ψi ðx; y; z þ ΔzÞ þ ψi ðx; y; z−ΔzÞ

ð17Þ
381
380

The solution has converged when the error ε = (ψi − 1 − ψi)/ψi − 1 is 382
below a threshold value of 10 −3. After generating the potential 383

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

P

R
O

O

F

6

R

389
388

397
398

O

C

395
396

U

siþ1 ¼ si þ hNx ðsi Þ þ hNy ðsi Þ þ hNz ðsi Þ:
400
399
401
402
403
404
405
406
407
408

410
409

ð19Þ

Three normalized potential difference maps are then created: Nx,
Ny, and Nz (Fig. 4b3 — blue vectors). From these maps, we calculate
gradient streamlines for every GM voxel. A streamline s is a vectorof points s1,…,sn that describes the path from the starting point s1 toa border. The following point, si + 1, of si is estimated by using theEuler's method, or by adding the weighted normalized gradientN(si) to si:

N

393
394

ð18Þ

R

Δψðx; y; zÞ=Δx ¼ ½ψðx þ Δx; y; zÞ−ψðx−Δx; y; zÞ=2:
390
391
392

We correct for errors introduced by the skeleton S using the vol- 411
ume difference between the uncorrected tissue segment SEGPF and 412
the corrected tissue segment SEGPFC.
413

T

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Nx ¼ ðΔψ=ΔxÞ= ðΔψ=ΔxÞ2 þ ðΔψ=ΔyÞ2 þ ðΔψ=ΔzÞ2 ;

C

386
387

image, we calculate the gradient ﬁeld N of the Laplace map as thesimple normalized two-point difference for each dimension. For example, along the x-direction the normalized potential difference Nx
is calculated as follows:

E

384
385

E

D

Fig. 4. Subﬁgure (a) shows a ﬂow diagram of the Laplacian approach, where subﬁgure (b1–b8 with simpliﬁed titles) shows 2D illustrations of the volume maps of subﬁgure (a). First, askeleton based on the WM distance map (see Fig. 3b2) is used to reconstruct blurred sulcal regions (b1–b3). Next, the Laplace equation is solved in the GM area and a vector ﬁeld N isgenerated (b3). This vector ﬁeld allows the creation of streamlines that follow the vectors to each boundary to measure the distance (b4–b6). To avoid an underestimation due to sulcus
reconstruction (b5), the CSF distance L(sCSF(v)) was corrected for changes from the sulcus reconstruction (Eq. (22)). The addition of both distance maps gives the cortical thickness map
GMT that allows the creation of the percentage position map PP, which in turn is used to create the CS and map cortical thickness onto the surface.

ð20Þ

The weight h describes the step size of the streamline calculation
and was set to 0.1 mm as a compromise between speed and quality.
For every GM voxel v, we calculate the streamline sWM(v) starting at
the position of v to the WM boundary and other streamline sCSF(v)from v to the CSF boundary. To calculate sCSF(v), it is necessary to
use the inverse gradient ﬁeld. The length of a streamline L(s) can befound by summing the Euclidean distance of all points si to their successor si + 1:

n−1LðsÞ ¼ ∑i¼1

rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2 
2 
2ﬃ
siþ1;x −si;x þ siþ1;y −si;y þ siþ1;z −si;z :

ð21Þ





LC ðsÞ ¼ LðsÞ þ SEGPF sn;x ; sn;y ; sn;z −SEGPFC sn;x ; sn;y ; sn;z :

ð22Þ
415
414

The summation of the length of both streamlines sWM(v) and 416
sCSF(v) gives the GM thickness at voxel v (Figs. 4b5 and b6). The 417
RPM can also be calculated using the values for the lengths of 418
sWM(v) and sCSF(v), with all WM voxels set to one:
419




GMT ðvÞ ¼ L sWMðvÞ þ LC sCSF ðvÞ ;

ð23Þ



PP ðvÞ ¼ LC sCSF ðvÞ =GMT ðvÞ þ ðSEGPFC ðvÞ > 2:5Þ:

ð24Þ

421
420

(Figs. 4b7 and b8). Finally, the CS surface is generated at a resolution 422
423
of 0.5 mm from the PP map and reduced to around 300,000 nodes 424
using standard Matlab functions. Each vertex of the mesh is assigned 425
a thickness value via a linear interpolation of the closest GMT map 426
values.
427
Spherical phantoms

428

A variety of spherical phantoms were used for validation. For the 429
standard gyral case, the spherical phantom consisted of a cortical 430
GM ribbon around a WM sphere in the center of the tissue map 431(Fig. 5). To explore the ability to reconstruct blurred sulcal regions, 432
a second spherical phantom was constructed such that it contained 433
a cortical GM ribbon sandwiched in between two WM regions: the 434
center sphere and an outer shell. Between the ribbon boundaries, a 435
small gap allows testing of the inﬂuence of the presence of CSF. To 436
simulate asymmetrical structures, the size of the second ribbon was 437
deﬁned as a ratio of the size of the ﬁrst GM ribbon. To realize this 438
phantom with PVE, a distance map SPD that measures the distance 439

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

7

P

R
O

O

F

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx:

1
r þ 0:5−SPDðvÞ
0; if SPDðvÞ≤ðr−0:5Þ
; if SPDðvÞ > ðr−0:5Þ & SPDðvÞbðr þ 0:5Þ;; if SPDðvÞ≥ðr þ 0:5Þ

ð25Þ

R

T PVE ðv; r Þ ¼

8
<

C

T

from the center of the volume at a resolution of 1 × 1 × 1 mm 3 is used
to create the tissue map TS:

E

440
441

E

D

Fig. 5. Spherical phantom validation matrix: Two spherical phantom types – one for the gyral case (CGW) and one for the blurred sulcal case (WGW) – were used with different
parameters for each method. The rows and columns are used to describe curvature (given by inner radius) and thickness values under varying conditions (see Table 1 and Fig. 7).

443
442

N
C
O

R

TSPVE ðv; r; t; sw; rspÞ ¼ inner  WM  sphere þ inner  GM  sphere þ CSF  sphereþ
outer  GM  sphere þ outer  WM  sphere¼ T PVE ðv; r Þ þ T PVE ðv; r þ t Þ þ 1
þð1−T PVE ðv; r þ t þ swÞÞþ
ð1−T PVE ðv; r þ t þ sw þ ððt=rspÞ  ð1−rspÞÞÞÞ

ð26Þ

453

Brain phantoms

To ensure an equal thickness for the brain phantom, it is necessary 454
to expand sulcal regions such that they are able to achieve full thick- 455
ness without intersections. To accomplish this, a CS of a healthy adult 456
test subject was generated with Caret (Van Essen et al., 2001) and 457
manually corrected for geometrical and topological errors (Fig. 6b). 458
Twenty iterations of weighted nearest neighbor surface-based 459
smoothing (included in the Caret package) were used to remove 460
high frequency structures that can lead to problems in later ma- 461
nipulation steps (Figs. 6c, a1). A graph-based distance measure 462
Dgb is used to create a distance map that describes the distance 463
with sub-voxel accuracy to a given iso-surface that was generated 464
via Matlab iso-surface functions. The initial surface was linearly in- 465
terpolated once to reduce missed measurements. This distance 466
map allows ﬁnding the new inner boundary at half distance 467(Fig. 6a2). From this new inner boundary, we estimated the new 468
outer boundary based on the distance map generated from the 469
inner boundary. If a sulcus is too small to allow increased thick- 470
ness without intersections, then it is blurred (Fig. 6a3). From this 471
new outer boundary, a new distance map allows the creation of 472
the ﬁnal central boundary (Fig. 6a4). The distance map BPD from 473
the central boundary now allows the creation of a segment image 474

t1:1
t1:2
t1:3
t1:4
t1:5

Table 1
Overview of parameters for the brain phantom test cases. For each test case, the anatomically expected range of each parameter was tested while all other parameters remained
constant. The bracketed values give the number of test cases for each parameter. In most cases, only one default value was used. Because the cortex contains both blurred (CGW)and non-blurred (WGW) regions, both cases were tested. When testing sulcal width and position, more than one default value was necessary due to high variance in the results. For
example, a symmetrical sulcus will produce better results than an asymmetric sulcus.

447
448
449
450

U

451
452

where v is voxel of TS, r gives the inner boundary radius, t describes the
thickness, sw is the sulcus width, and rsp is the relative sulcus position.
Thickness is only evaluated for the inner ribbon because for asymmetrical structure the outer ribbon has a different thickness and curvature.
Table 1 shows the values of each parameter to be tested. To test
one parameter, all other parameters were ﬁxed to standard values.
The range of values chosen for the parameters was based on anatomical and technical considerations.

444
445
446

t1:6

Parameter

CurvatureThickness

PVE vs no PVE

V vs. SType

Sulcus width

rel. sulcus pos

t1:7
t1:8

RangesDefaults

1.0:0.01:5.0 (401)
2.5 (1)

0.0:0.01:5.0 (501)
2.5 (1)0 1 (2)1 (1)V S (2)S (1)CGW WGW (2)CGW WGW (2)0.0:0.01:2.0 (200)0.0:0.50:1.0 (3)0.2:0.001:0.8 (601)
0.3:0.2:0.7 (3)Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

R
O

O

F

8

with WM, GM, and CSF (Fig. 5a5) (for a resolution of 0.5 × 0.5 ×
0.5 mm3):
8
3; if BPDðvÞ≥ð−t=2−0:25Þ
>>>>< −t=2−0:25 þ BPDðvÞ ; if BPDðvÞbð−t=2−0:25Þ & BPDðvÞ > ð−t=2 þ 0:25Þ
TBPVE ðv; r Þ ¼
2; if BPDðvÞ≤ð−t=2 þ 0:25Þ & BPDðvÞ≥ðt=2−0:25Þ :>>t=2 þ 0:25−BPDðvÞ ; if BPDðvÞbðt=2−0:25Þ & BPDðvÞ > ðt=2 þ 0:25Þ
>>:1; if BPDðvÞ≤ðt=2 þ 0:25Þ

T

ð27Þ
477
478
479

applicable, to FreeSurfer. The ﬁrst test consisted of the set of spherical 506
phantoms, which were used to test the approaches over a wide set of 507
parameters under simple but precise conditions. The second test, 508
consisting of the brain phantoms, was used to explore the performance 509
of the approaches under the more realistic condition of a highly convo- 510
luted surface with equal thickness. For the third test, we used the Collins 511
phantom with different noise levels, both to add more realism and to di- 512
rectly compare the results to the FreeSurfer software package. Finally, 513
we used real MR data of one subject for a test–retest of all three 514
methods.
515

D

476

E

475

P

Fig. 6. Brain phantom generation: Subﬁgure (a) illustrates the generation process for the brain phantom, the difference between the original and ﬁnal surface (a.I), and problems forhigher thickness levels (a.II), whereas (b) to (f) show the changes from the individual surface to the brain phantom. A smoothed individual surface of a healthy adult (a.1) istransformed by distance operations to a surface that allows the creation of a t-mm thick cortical ribbon (a.5). This process removes high-frequency WM structures (a.2) and enlarges sulcal
regions (a.3) to ensure an actual thickness level between 0.5 and 4.0 mm. Larger thicknesses destroy most of sulcal structures of a normally folded brain (a.II).

484

Collins phantoms

485
486

489

To test different thickness levels on one surface and stability forimages interferences, 6 BrainWeb T1-weighted phantom datasets(1-mm resolution) with 1%, 3%, 5%, 7%, 9% noise and 20% inhomogeneity
were compared to a dataset without noise and inhomogeneity (Collins
et al., 1998).

490

Real data

491
492

501
502

The sample data set included 12 brain scans of the same healthy
adult subject performed on two different 1.5 T Siemens Vision scanners
within one year. Both scanners used 3D magnetization prepared gradientecho (MP-RAGE) T1-weighted sequences of 160 sagittal slices with voxel
dimensions 1×1×1 mm and FOV=256 mm. Scanner 1 parameters
were TR/TE/FA=11.4 ms/4.4 ms/15°, and Scanner 2 parameters were
TR/TE/FA=15 ms/5 ms/30°.
Each reconstructed surface of the 12 scans was compared to an average scan to estimate the surface reconstruction and thickness errors,
similar to the analysis used for the Collins phantom. Ideally, all reconstructions should be identical, and they should produce identical thickness measurements.

Over all test parameters, PBT shows better results than the Laplacian 517
approach for both thickness estimation and surface generation (Fig. 7a). 518
As expected, both methods have higher RMS error for thickness estima- 519
tion than for surface generation, both produce better results with PVE, 520
and both perform better for the simpler gyral case compared to the sulcal 521
case. The voxel-based results of the Laplacian approach are much worse 522
than after projection to the surface, whereas PBT produced equally accu- 523
rate results due to the smoothness parameter of the projection. Com- 524
pared to gyral regions, sulcal regions show higher RMS error, which is 525
strongly related to the width of the sulcal gap and its relative position. 526
Predicted by the sampling theorem, both show a strong increase of 527
RMS error below sampling resolution for thickness measurements, 528
but not for surface generation (Figs. 7b and c).
529
Furthermore, the Laplacian approach had larger ﬂuctuations of error 530
across the test cases. Relatively small variations of the test parameters 531
led to vastly different error values (Fig. 7c above 2.5 mm). This strong 532
variation can only be found if the step size of the parameter is very 533
small — around 0.01 mm. Especially, asymmetrical structures (Fig. 7d) 534
and small sulcal gaps (Fig. 7e) vastly increase the RMS error of the 535
Laplacian approach.
536
For the Laplacian approach, we were able to produce good results 537
such as those published in the literature only for cases with relatively 538
large CSF regions and low asymmetry, whereas PBT produced more 539
exact and stable results over the full range of test parameters.
540

503

Results

Brain phantoms

504

Four different test matrices were used to validate PBT; these results were then compared to the Laplacian approach and, wherever

For the brain phantoms, if the two approaches are compared, the PBT 542
method has much lower RMS error for the thickness measurements and 543

493
494
495
496
497
498
499
500

505

E

R

R

O

C
N

487
488

U

481
482

C

483

The default parameters (2.5 mm thickness, 1 × 1 × 1 mm3 resolution, 0% noise) were modiﬁed individually, resulting in 14 thickness
levels, 9 noise degrees, and 7 isotropic and 7 anisotropic grid resolutions. The dataset is available under http://dbm.neuro.uni-jena.de/phantom/.

480

Spherical phantoms

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

516

541

9

N
C
O

R

R

E

C

T

E

D

P

R
O

O

F

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

Fig. 7. Spherical phantom: PBT results in lower RMS error for all test categories, compared to the Laplacian approach (a). Below sampling resolution, both methods show a predictable
increase of thickness measurement error due to the sampling theorem (c), whereas the position error stays stable (b). Most errors happen for sulcal cases with low sulcus
width (d) and higher asymmetries (e).

544
545

555

Collins phantom

556
557

The advantage of using an additional Collins phantom is that thetwo approaches described here (PBT and Laplacian) can be compared
to a commonly used approach for both reconstructing cortical surfaces

546
547
548
549
550
551
552

558

U

553
554

similar RMS error for the surface position compared to the Laplacian
method (Fig. 8). The errors occur mostly in sulcal regions where the sulcus reconstruction step cut strongly into the fundi of the sulci such thatthat a complete correction was not possible (Fig. 8a). However, using aweaker sulcus reconstruction step or stronger corrections led to thickness
overestimation, more defects, and greater RMS errors, thus it was impossible to circumvent this problem. Generally, the largest errors occurred
for anisotropic resolutions, thickness levels below the sampling resolution, and higher noise levels. It can be assumed that these factors would
apply to any cortical data set, and thus should be considered before applying any cortical reconstruction method.

and measuring thickness, e.g., FreeSurfer. To summarize the ﬁndings, 559
the PBT approach had comparable or lower RMS error compared to 560
both the Laplacian and FreeSurfer approaches (Fig. 7, Supplementary 561
Fig. A1). If the noise level is increased, all thickness measures also 562
had increasing error. A two-sample unpaired t-test showed no signiﬁ- 563
cant differences of the RMS position error between PBT and Laplacian 564(t= −0.048, df = 8, p = 0.963) and PBT and FreeSurfer (t= 1.348, 565
df=8, p = 0.215). A signiﬁcant difference in thickness between PBT 566
and Laplacian was found (t= −2.95, df= 8, p = 0.019), but not for 567
PBT vs. FreeSurfer (t= −0.944, df=8, p = 0.374). Furthermore, the 568
PBT method provides an advantage in terms of reduced numbers of to- 569
pological defects (an average of 15.1 for PBT, compared to 28.2 for 570
Laplacian and 18.5 for FreeSurfer). PBT had signiﬁcantly fewer defects 571
compared to the Laplacian approach (t= −8.656, df= 10, p b 0.001), 572
but not compared to FreeSurfer (t= −1.481, df = 10, p = 0.182). The 573
defects associated with the Laplacian and PBT approaches were mostly 574
bridges between two gyri and were removed by the topology 575

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

P

R
O

O

F

10

576

T

E

D

Fig. 8. Brain phantom: Subﬁgure (a) shows the resulting surfaces for a simulated thickness of 2.5 mm, with an isotropic resolution of 1 × 1 × 1 mm3 and no noise. PBT produced
overall good results (left), whereas the Laplacian approach showed strong underestimation in sulcal regions due to the sulcus reconstruction step (right), even if sulcus error correction was used (middle). The Laplacian approach (b — red) produced much higher thickness RMS errors than PBT (b — blue). Low sample resolution, anisotropic resolutions, andnoise may increase the RMS error for both thickness measurements as well as CS position (c–e). (For interpretation of the references to color in this ﬁgure legend, the reader isreferred to the web version of this article.)

581

Twelve scans of one subject

582

As a ﬁnal approach for quantifying the performance of three approaches (PBT, Laplacian, FreeSurfer), we analyzed twelve separatescans of a single brain, then compared the results to an averaged scanof the same brain. Since the elapsed time between scans was less than
one year, cortical thickness should be unchanged. Again, the PBT approach provided some advantages over the other methods (Fig. 10, Supplementary Fig. A2). First, the PBT approach is comparable to or better
than other approaches in terms of the RMS thickness measurement errors (Fig. 10c; PBT: 0.39±0.02 mm; Laplacian: 0.64±0.02 mm;
FreeSurfer: 0.53±0.05 mm), and the RMS position error of the CS reconstructions was similar to the other two approaches (Fig. 10d; PBT:
0.50±0.05 mm; Laplacian: 0.54±0.05 mm; FreeSurfer: 0.60±
0.23 mm). There was no signiﬁcant difference in the RMS position
error between PBT and the Laplacian approach (t=−1.922, df=22,
p=0.067) and PBT and FreeSurfer (t=−1.409, df=22, p=0.172),
whereas the difference of the RMS thickness error was signiﬁcant(t=−8.177, df=22, pb 0.001). A major difference between the PBTand Laplacian approaches compared to FreeSurfer is a general underestimation of thickness in the motor cortex (Fig. 10b). Finally, the PBT approach produced far fewer topological defects per hemispherecompared to Laplacian (t=−6.036, df=24, pb 0.001) and FreeSurfer(t=−4.030, df=24, pb 0.001) (Fig. 10a; PBT: 21.5; Laplacian: 34.6;
FreeSurfer: 54.6).
Inline supplementary Fig. A2 can be found online at http://dx.doi.
org/10.1016/j.neuroimage.2012.09.050.

589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606

E

R

R

O

C

587
588

N

585
586

U

583
584

C

579
580

correction. These results were highly dependent upon the quality of
the initial tissue segmentation, the implications of which are discussed
more fully in the Discussion section (Fig. 9).
Inline supplementary Fig. A1 can be found online at http://dx.doi.
org/10.1016/j.neuroimage.2012.09.050.

577
Q7 578All calculations were done on an iMac 3.4 GHz Intel Core i7 with 607
8 GB RAM and Matlab 7.12. For both hemispheres with a resolution 608
of 0.5 mm and with topology error correction, PBT needed around 609
20 min, whereas the Laplacian approach takes around 2 h. Although 610
the FreeSurfer processing pipeline is structured differently than the 611
PBT and Laplacian approaches, rendering comparison difﬁcult, an es- 612
timate of the time to perform cortical reconstruction and thickness 613
measurement is several hours. 6
614

Discussion

615

For nearly all test cases, PBT had much lower thickness and position 616
errors than the Laplacian approach, because PBT uses an inherent model 617
that detects sulci, whereas the Laplacian method requires an explicit 618
sulcus reconstruction step that changes the tissue class of sulcal voxels 619
and may lead to the introduction of additional errors, even if these tis- 620
sue class changes are compensated for within the algorithm. The differ- 621
ent tests of the spherical phantom clearly show that the strong errors of 622
the Laplacian approach only happen in asymmetric sulcal regions, al- 623
though both methods are based on the same Eikonal distance measure 624
that accounts for the sulcal gap. Because the real cortex also contains 625
asymmetrical structures, it is important that the thickness measure is 626
able to accurately evaluate these asymmetries (Das et al., 2009; Fischl 627
and Dale, 2000; Kim et al., 2005). In addition, the brain phantoms indi- 628
cate errors on the fundi of the sulci for the Laplacian method, whereas 629
the continuous model of PBT allows a stable estimation over the whole 630
cortex.
631
The correct reconstruction of blurred sulci is still a challenging 632
process, since the result depends strongly on the used method and 633
its parameters (Acosta et al., 2008, 2009; Cardoso et al.; Dale et al., 634
1999; Das et al., 2009; Han et al., 2004; Hutton et al., 2008; Kim et 635
6

https://surfer.nmr.mgh.harvard.edu/fswiki/ReconAllRunTimes

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

11

E

D

P

R
O

O

F

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664

R

R

639
640

al., 2005; MacDonald et al., 2000; Scott and Thacker, 2005). The
shown results of the different methods allow only a rough impression
about the quality of the sulcus reconstruction step, by showing thatthe modeling of sulcal blurring leads to results that are closer to thesimulated cortical thickness. Most differences in these methods arevisible especially on the fundi of the sulci, which is where some approaches have thickness over- or underestimation. Although otherauthors, i.e. (Das et al., 2009; Hutton et al., 2008; Kim et al., 2005), illustrate the reconstruction of blurred regions for principle examples, this is
the ﬁrst paper that introduces a way to numerically validate an algorithm using not only simple cases with well-known parameters andwithout fundi, but also for highly convoluted surfaces with fundi.
Compared to surface-based approaches, PBT does not need extensive surface deformation or self-intersection tests, which are necessary
for both Freesurfer (Dale et al., 1999) and CLASP (Kim et al., 2005). In
contrast to FreeSurfer, PBT is able to use tissue segmentation images
produced using any segmentation approach, allowing a separate development of the segmentation algorithms and thus making this process
more transparent. Furthermore, this allows the use of segmentation images for other imaging modalities such as T2, PD (Ashburner andFriston, 2005; Zhang et al., 2001), DTI (Liu et al., 2007), and othermethods that take account of special contrast properties in diseasestates such as multiple sclerosis (Khayati et al., 2008; Wu et al., 2006),
white matter hyper-intensities (Admiraal-Behloul et al., 2005; Gibson
et al.), or tumors (Kaus et al., 2001; Prastawa et al., 2004), or for otherspecies (Andersen et al., 2002). As a result, the input of PBT and othersegment-based methods depends strongly on the results of thesegmentation. The tests with the spherical and brain phantoms were independent from the segmentation process, because the segmentation

N
C
O

637
638

U

636

E

C

T

Fig. 9. Collins phantom: Diagram (a) shows the mean number of defects per hemisphere for PBT (blue), Laplacian (red), and FreeSurfer (green). Subﬁgures (b–d) show the groundtruth surface of the Collins phantom noise test for PBT (left surface), Laplacian (middle surface), and FreeSurfer (right surface), in which the color map codes the cortical thickness ofthe ground truth surface (b), the mean thickness RMS error of all noise levels compared to the ground truth surface (c), and the mean distance RMS error of all noise levels to theground truth surface (d). A supplementary ﬁgure, including medial and lateral views of both hemispheres, is available online. (For interpretation of the references to color in thisﬁgure legend, the reader is referred to the web version of this article.)

images were directly simulated, whereas the Collins phantom and the 665
real dataset include a segmentation step. Evidence of the strong inﬂu- 666
ence of the segmentation algorithm on results may be seen with the 667
Collins phantom. Since the tissue boundaries are simulated, these 668
phantoms included artiﬁcially precise tissue classiﬁcation and resulted 669
in much more similar thickness measurements for all methods than 670
for the real data set (especially in the motor cortex).
671
Furthermore, PBT allows a direct voxel-based analysis, potentially 672
in combination with other voxel-based data (Hutton et al., 2009), and 673
it may also be used to measure the thickness of the WM and CSF 674[HBM2010]. The voxel-based thickness estimation of PBT and other 675
methods allows the easy creation of the central surface, which has 676
better properties than the WM or pial surface. Previous approaches 677
generally reconstruct a surface at a tissue boundary, which is either 678
the WM surface or the pial surface. In one sense, such a reconstruc- 679
tion makes sense, since the intensity gradient in these regions can 680
be used to estimate the position of the surface. However, due to the 681
PVE, the boundaries often contain voxels with more than one tissue 682
class and thus render it impossible to determine the precise location 683
of the surface within that voxel. In the approach suggested here, the 684
effect of PVE is somewhat reduced, since the central surface is 685
reconstructed simply at the 50% distance boundary between the 686
GM/WM and GM/CSF boundaries. This effect is responsible for the 687
constant RMS position error below the sampling resolution, whereas 688
thickness errors grow much stronger, because the PVE and neighbor 689
information can only code the exact position of one boundary. For in- 690
stance, for the 1D case of a voxel v and its left and right neighbors vl 691
and vr, where v = 2.25, vl = 2, and vr = 3, the WM boundary is exactly 692
described between v and vr, but if vl = 1, then there are two 693

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

D

P

R
O

O

F

12

703
704
705
706
707
708
709
710
711
712
713
714
715
716

Q8 717
718
719
720
721
722
723
724

E

R

R

701
702

O

699
700

C

698

N

696
697

boundaries within v and it is unclear how much GM is within v. It is
possible that there is only GM and WM in v, the WM boundary is at
the same position, and the CSF boundary is exactly between v and vl.
But it is also possible that there is some CSF in v, and v contains three tissue classes and both boundaries. As a result, thickness RMS errors growstrongly for thickness levels below the sampling resolution.
Independent of the chosen reconstruction method, the general
structure of the CS lends advantages that do not exist in the other surface
reconstructions. First, the CS has a lower “frequency” content, or fewer
ﬁnely detailed regions, since it is the average of the WM structure with
its strong gyri and the pial surface with its deep sulci. Due to this characteristic of the CS, it may have fewer topological defects and it tends to loseless anatomical detail when smoothed. Since brain surfaces usually mustbe smoothed to remove stair artifacts and noise, the CS provides a distinct
advantage over the other reconstructions. Secondly, another advantageof the CS is that it may be directly reconstructed from the data and thus
leads to a more uniform distribution of vertices across the surface,which may be perturbed in a method that uses a deformation process
to reconstruct a surface at a tissue boundary.
Before performing intersubject comparisons, the brain surface
meshes must usually be free of topological defects, and there are several
approaches available to retrospectively correct topological errors either
in volume space or directly on the surface (Kriegeskorte and Goebel,
2001; Segonne et al., 2007; Shattuck and Leahy, 2001; Yotter et al.,2009, in press). Despite the availability of these correction methods, it
is desirable to minimize both the size and number of topological defects,
since non-idealities in the correction step can often introduce errors. In
this respect, the PBT approach is the best choice, since it produced thelowest number of defects, and the defects were also relatively small. Despite using the same segmentation images, the Laplacian approach
resulted in a large number of defects, mostly due to overestimation of

U

694
695

C

T

E

Fig. 10. Real data: Diagram (a) show the mean number of defects per hemisphere for PBT (blue), Laplacian (red), and FreeSurfer (green). Shown in (b–d) are PBT (left column),
Laplacian (middle column), and FreeSurfer (right column) surfaces with (b) cortical thickness calculated for the average surface, (c) mean thickness RMS error of all scans compared to the thickness of the average surface, and (d) mean distance RMS error of all scans to the average surface. Strong differences are visible in the thickness measurements
for the motor cortex, which depended mostly on the segmentation algorithm (VBM8) and thus was similar for PBT and Laplacian, whereas FreeSurfer used internal routines. A supplementary ﬁgure, including medial and lateral views of both hemispheres, is available online. (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)thickness in sulcal regions and thus the formation of bridges. A detailed 725
discussion of the corrections of topology defects via spherical har- 726
monics can be found in (Yotter et al., in press).
727 Q9Necessity of a full phantom test suite

728

Comparing different software packages is never easy, because there 729
will always remain some differences in processing the data, i.e. the re- 730
striction of FreeSurfer to 1.0 mm resolution for all volumes whereas 731
PBT and the Laplacian approach can also use higher resolutions (here 732
0.5 mm). Especially, the different segmentation routines limited the 733
comparison between FreeSurfer and both other approaches. Further- 734
more, all methods based on different thickness deﬁnitions can also 735
lead to slightly different results (Lerch and Evans, 2005; MacDonald et 736
al., 2000).
737
Because visual inspection of surfaces gives only subjective, badly 738
reproducible, and often limited impressions of the reconstruction 739
quality (Kabani et al., 2001; Xu et al., 1999), we developed a complete 740
test suite containing several parameters that could be varied to fully 741
characterize both surface reconstruction and thickness measurement 742
approaches. Although previous approaches tested a small number of 743
phantom objects (Acosta et al., 2009; Das et al., 2009; Miller et al., 744
2000), it is apparent from our results that it is necessary to test several 745
parameters to gain information about an algorithm's performance, 746
especially for special cases such as sulcal blurring. It could be further 747
argued that simple geometrical objects provide only limited infor- 748
mation about performance that cannot be extrapolated to cortical 749
surfaces, thus it is appropriate to include pseudo-cortical surfaces 750
with constant thickness over the whole cortex in the test suite. Unlike 751
the previous methods (Liu et al., 2008), our cortical ribbon has an 752
equal thickness and a more realistic structure. This constant thickness 753

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

754
755

C

E

R

R

Conclusion

796

811

In this paper, we have presented a new method that allows for thesimultaneous reconstruction of the CS and measurement of cortical
thickness. Our PBT method is based on (probability) maps of a standard
CSF–GM–WM tissue segmentation and has several advantages over theprevious methods, such as direct estimation of the CS, comparable or
lower errors, and fewer topological defects. We introduce a framework
for thoroughly validating methods developed for surface reconstruction
and thickness estimation, which quantiﬁes the performance of themethods over a wide range of thickness levels and other parameterssuch as sampling resolution, noise, curvature, and PVE. The test framework explores performance both for the simple case of a sphere andalso for nearly normal folded cortices with uniform thickness. Finally,
we used real MR images from several scans of the same subject to compare both methods to FreeSurfer. The results indicate that the quality of
our CS reconstructions and thickness estimations is comparable, andmay be superior in certain respects, to other methods.

812

Acknowledgments

813

This work was supported by the following grants: BMBF 01EV0709
and BMBF 01GW0740.

799
800
801
802
803
804
805
806
807
808
809
810

814

U

797
798

N
C
O

795

815

References

D

P

R
O

O

F

Acosta, O., Bourgeat, P., Fripp, J., Bonner, E., Ourselin, S., Salvado, O., 2008. Automatic 816
delineation of sulci and improved partial volume classiﬁcation for accurate 3D 817
voxel-based cortical thickness estimation from MR. Med. Image Comput. Comput. 818
Assist. Interv. 11, 253–261.
819
Acosta, O., Bourgeat, P., Zuluaga, M.A., Fripp, J., Salvado, O., Ourselin, S., 2009. Automated 820
voxel-based 3D cortical thickness measurement in a combined Lagrangian–Eulerian 821
PDE approach using partial volume maps. Med. Image Anal. 13, 730–743.
822
Admiraal-Behloul, F., van den Heuvel, D.M., Olofsen, H., van Osch, M.J., van der Grond, J., 823
van Buchem, M.A., Reiber, J.H., 2005. Fully automatic segmentation of white matter 824
hyperintensities in MR images of the elderly. NeuroImage 28, 607–617.
825
Andersen, A.H., Zhang, Z., Avison, M.J., Gash, D.M., 2002. Automated segmentation of 826
multispectral brain MR images. J. Neurosci. Methods 122, 13–23.
827
Ashburner, J., Friston, K.J., 2005. Uniﬁed segmentation. NeuroImage 26, 839–851.
828
Aubert-Broche, B., Evans, A.C., Collins, L., 2006. A new improved version of the realistic 829
digital brain phantom. NeuroImage 32, 138–145.
830
Bouix, S.S., Kaleem, 2000. Divergence-based medial surfaces: ECCV, 1, pp. 603–618.
831
Cardoso, M.J., Clarkson, M.J., Ridgway, G.R., Modat, M., Fox, N.C., Ourselin, S., LoAd: A locally 832
adaptive cortical segmentation algorithm. Neuroimage 56, 1386–1397.
833 Q11
Collins, D.L., Zijdenbos, A.P., Kollokian, V., Sled, J.G., Kabani, N.J., Holmes, C.J., Evans, A.C., 834
1998. Design and construction of a realistic digital brain phantom. IEEE Trans. Med. 835
Imaging 17, 463–468.
836
Coupe, P., Yger, P., Prima, S., Hellier, P., Kervrann, C., Barillot, C., 2008. An optimized 837
blockwise nonlocal means denoising ﬁlter for 3-D magnetic resonance images. IEEE 838
Trans. Med. Imaging 27, 425–441.
839
Cuadra, M.B., Cammoun, L., Butz, T., Cuisenaire, O., Thiran, J.P., 2005. Comparison and val- 840
idation of tissue modelization and statistical classiﬁcation methods in T1-weighted 841
MR brain images. IEEE Trans. Med. Imaging 24, 1548–1565.
842
Dale, A.M., Fischl, B., Sereno, M.I., 1999. Cortical surface-based analysis. I. Segmentation 843
and surface reconstruction. NeuroImage 9, 179–194.
844
Das, S.R., Avants, B.B., Grossman, M., Gee, J.C., 2009. Registration based cortical thickness 845
measurement. NeuroImage 45, 867–879.
846
Eskildsen, S.F., Ostergaard, L.R., 2006. Active surface approach for extraction of the human 847
cerebral cortex from MRI. Med. Image Comput. Comput. Assist. Interv. 9, 823–830. 848
Eskildsen, S.F., Ostergaard, L.R., 2007. Quantitative comparison of two cortical surface 849
extraction methods using MRI phantoms. Med. Image Comput. Comput. Assist. 850
Interv. 10, 409–416.
851
Fischl, B., Dale, A.M., 2000. Measuring the thickness of the human cerebral cortex from 852
magnetic resonance images. Proc. Natl. Acad. Sci. U. S. A. 97, 11050–11055.
853
Fischl, B., Sereno, M.I., Dale, A.M., 1999. Cortical surface-based analysis. II: inﬂation, 854
ﬂattening, and a surface-based coordinate system. NeuroImage 9, 195–207.
855
Fjell, A.M., Walhovd, K.B., Reinvang, I., Lundervold, A., Salat, D., Quinn, B.T., Fischl, B., 856
Dale, A.M., 2006. Selective increase of cortical thickness in high-performing elderly– 857
structural indices of optimal cognitive aging. NeuroImage 29, 984–994.
858
Gibson, E., Gao, F., Black, S.E., Lobaugh, N.J., Automatic segmentation of white matter 859
hyperintensities in the elderly using FLAIR images at 3T. J. Magn. Reson. Imaging 860
31, 1311–1322.
861 Q12Haidar, H., Soul, J.S., 2006. Measurement of cortical thickness in 3D brain MRI data: 862
validation of the Laplacian method. J. Neuroimaging 16, 146–153.
863
Han, X., Pham, D.L., Tosun, D., Rettmann, M.E., Xu, C., Prince, J.L., 2004. CRUISE: cortical 864
reconstruction using implicit surface evolution. NeuroImage 23, 997–1012.
865
Hofman, M.A., 1989. On the evolution and geometry of the brain in mammals. Prog. 866
Neurobiol. 32, 137–158.
867
Hutton, C., De Vita, E., Ashburner, J., Deichmann, R., Turner, R., 2008. Voxel-based cortical 868
thickness measurements in MRI. NeuroImage 40, 1701–1710.
869
Hutton, C., Draganski, B., Ashburner, J., Weiskopf, N., 2009. A comparison between 870
voxel-based cortical thickness and voxel-based morphometry in normal aging. 871
NeuroImage 48, 371–380.
872
Jones, S.E., Buchbinder, B.R., Aharon, I., 2000. Three-dimensional mapping of cortical 873
thickness using Laplace's equation. Hum. Brain Mapp. 11, 12–32.
874
Kabani, N., Le Goualher, G., MacDonald, D., Evans, A.C., 2001. Measurement of cortical 875
thickness using an automated 3-D algorithm: a validation study. NeuroImage 13, 876
375–380.
877
Kaus, M.R., Warﬁeld, S.K., Nabavi, A., Black, P.M., Jolesz, F.A., Kikinis, R., 2001. Automated 878
segmentation of MR images of brain tumors. Radiology 218, 586–591.
879
Khayati, R., Vafadust, M., Towhidkhah, F., Nabavi, M., 2008. Fully automatic segmenta- 880
tion of multiple sclerosis lesions in brain MR FLAIR images using adaptive mixtures 881
method and Markov random ﬁeld model. Comput. Biol. Med. 38, 379–390.
882
Kim, J.S., Singh, V., Lee, J.K., Lerch, J., Ad-Dab'bagh, Y., MacDonald, D., Lee, J.M., Kim, S.I., 883
Evans, A.C., 2005. Automated 3-D extraction and evaluation of the inner and outer 884
cortical surfaces using a Laplacian map and partial volume effect classiﬁcation. 885
NeuroImage 27, 210–221.
886
Kriegeskorte, N., Goebel, R., 2001. An efﬁcient algorithm for topologically correct segmen- 887
tation of the cortical sheet in anatomical mr volumes. NeuroImage 14, 329–346.
888
Kuperberg, G.R., Broome, M.R., McGuire, P.K., David, A.S., Eddy, M., Ozawa, F., Goff, D., 889
West, W.C., Williams, S.C., van der Kouwe, A.J., Salat, D.H., Dale, A.M., Fischl, B., 890
2003. Regionally localized thinning of the cerebral cortex in schizophrenia. Arch. 891
Gen. Psychiatry 60, 878–888.
892
Lee, J., Lee, J.M., Kim, J.H., Kim, I.Y., Evans, A.C., Kim, S.I., 2006a. A novel quantitative val- 893
idation of the cortical surface reconstruction algorithm using MRI phantom: issues 894
on local geometric accuracy and cortical thickness. Med. Image Comput. Comput. 895
Assist. Interv. Int. Conf. Med. Image Comput. Comput. Assist. Interv. 9, 183–190.
896
Lee, J.K., Lee, J.M., Kim, J.S., Kim, I.Y., Evans, A.C., Kim, S.I., 2006b. A novel quantitative 897
cross-validation of different cortical surface reconstruction algorithms using MRI 898
phantom. NeuroImage 31, 572–584.
899

E

T

theoretically allows a direct comparison between different thickness
measurement algorithms.
756
Using phantoms with equal thickness has the fundamental advantage
757
that an equal ribbon allows theoretically similar thickness measurements,758
independent of the deﬁnition of the thickness measure. An illustration
Q10 759 may clarify this point. Let t be the simulated thickness of a convoluted
760
brain-like ribbon with equal thickness. First, for nearest-neighbor-based
761
methods (i.e. Tnear (MacDonald et al., 2000) for surface-based methods
762
or nearest voxel for voxel-based methods), it is obvious that the nearest
763
connection between both sides is given by the deﬁned thickness t. Second,764
the Tnormal (MacDonald et al., 2000) metric that measures the distance be765
tween both sides of the ribbon via the surface normal will measure the766same thickness t, because of the well-deﬁned structure of this ribbon,
767
i.e., both boundaries have the same curvature by deﬁnition. Third, the768
streamline of the Laplacian approach will be equal to the surface normal,
769
because they depend on the vector ﬁeld given by the Laplace ﬁlter, which
770
in turn depends on the curvature of both boundaries that are equal by771
deﬁnition. Fourth, the Tlink (MacDonald et al., 2000) metric is deﬁned772
for surfaces with equal numbers of vertices. Here, one surface is the re773
sult of a deformation of the other surface. The deformation is mostly
774
based upon a ﬁeld given by the intensity and/or by the surface normal
775
or another Laplace vector ﬁeld (Kim et al., 2005). Because the intensity776
is equal within the ribbon, only the surface normal or the vector ﬁeld
777
can be used for the deformation. As a result, the deformation is similar778
to the streamlines of the Laplacian approach that are similar to the sur779
face normal.
780
The PVE approximation of the phantom generation based on distance781
maps leads to errors that depend on the resolution, the intensity (given782
by the distance), and the angles of the voxel to the coordinate system.
783
The highest possible error for a resolution of 1×1×1 mm3 happens for784
a diagonal voxel within the middle slice, and is, with a volume error
785
below 0.05 mm 3, comparable to other approximation methods
786(Acosta et al., 2009) in which the object is rendered ﬁrst to787
0.1 × 0.1 × 0.1 mm 3 and then down-sampled back to 1 × 1 × 1 mm 3.
788
The advantage of using distance maps is the much lower memory de789
mand and faster computation.
790
In the approach used here, segmentation images were directly
791
simulated to avoid inﬂuences from the segmentation algorithm.
792
However, it is possible to simulate a T1 image based on the tissue793
maps (Aubert-Broche et al., 2006), which would be useful for testing794
other methodological approaches using this test suite.

13

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

D

P

R
O

O

F

Shattuck, D.W., Leahy, R.M., 2001. Automated graph-based analysis and correction of 943
cortical volume topology. IEEE Trans. Med. Imaging 20, 1167–1177.
944
Sowell, E.R., Thompson, P.M., Leonard, C.M., Welcome, S.E., Kan, E., Toga, A.W., 2004. 945
Longitudinal mapping of cortical thickness and brain growth in normal children. 946
J. Neurosci. 24, 8223–8231.
947
Sowell, E.R., Peterson, B.S., Kan, E., Woods, R.P., Yoshii, J., Bansal, R., Xu, D., Zhu, H., 948
Thompson, P.M., Toga, A.W., 2007. Sex differences in cortical thickness mapped in 949
176 healthy individuals between 7 and 87 years of age. Cereb. Cortex 17, 1550–1560. 950
Thompson, P.M., Hayashi, K.M., Sowell, E.R., Gogtay, N., Giedd, J.N., Rapoport, J.L., de 951
Zubicaray, G.I., Janke, A.L., Rose, S.E., Semple, J., Doddrell, D.M., Wang, Y., van Erp, 952
T.G., Cannon, T.D., Toga, A.W., 2004. Mapping cortical change in Alzheimer's disease, 953
brain development, and schizophrenia. NeuroImage 23 (Suppl. 1), S2–S18.
954
Tohka, J., Zijdenbos, A., Evans, A., 2004. Fast and robust parameter estimation for statis- 955
tical partial volume models in brain MRI. NeuroImage 23, 84–97.
956
Toro, R., Perron, M., Pike, B., Richer, L., Veillette, S., Pausova, Z., Paus, T., 2008. Brain size 957 Q13
and folding of the human cerebral cortex. Cereb. Cortex.
958
Tosun, D., Rettmann, M.E., Han, X., Tao, X., Xu, C., Resnick, S.M., Pham, D.L., Prince, J.L., 2004. 959
Cortical surface segmentation and mapping. NeuroImage 23 (Suppl. 1), S108–S118.
960
Van Essen, D.C., Drury, H.A., 1997. Structural and functional analyses of human cerebral 961
cortex using a surface-based atlas. J. Neurosci. 17, 7079–7102.
962
Van Essen, D.C., Drury, H.A., Dickson, J., Harwell, J., Hanlon, D., Anderson, C.H., 2001. An 963
integrated software suite for surface-based analyses of cerebral cortex. J. Am. Med. 964
Inform. Assoc. 8, 443–459.
965
Wu, Y., Warﬁeld, S.K., Tan, I.L., Wells III, W.M., Meier, D.S., van Schijndel, R.A., Barkhof, 966
F., Guttmann, C.R., 2006. Automated segmentation of multiple sclerosis lesion 967
subtypes with multichannel MRI. NeuroImage 32, 1205–1215.
968
Xu, C., Pham, D.L., Rettmann, M.E., Yu, D.N., Prince, J.L., 1999. Reconstruction of the 969
human cerebral cortex from magnetic resonance images. IEEE Trans. Med. Imaging 970
18, 467–480.
971
Yezzi Jr., A.J., Prince, J.L., 2003. An Eulerian PDE approach for computing tissue thickness. 972
IEEE Trans. Med. Imaging 22, 1332–1339.
973
Yotter, R.A., Dahnke, R., Gaser, C., 2009. Topological correction of brain surface meshes 974
using spherical harmonics. Med. Image Comput. Comput. Assist. Interv. 12, 125–132. 975
Yotter, R.A., Dahnke, R., Gaser, C., in press. Topological Correction of Brain Surface 976
Meshes Using Spherical Harmonics. Hum. Brain Mapp.
977 Q14
Zeng, X., Staib, L.H., Schultz, R.T., Duncan, J.S., 1999. Segmentation and measurement of 978
the cortex from 3-D MR images using coupled-surfaces propagation. IEEE Trans. 979
Med. Imaging 18, 927–937.
980
Zhang, K., Sejnowski, T.J., 2000. A universal scaling law between gray matter and white 981
matter of cerebral cortex. Proc. Natl. Acad. Sci. U. S. A. 97, 5621–5626.
982
Zhang, Y., Brady, M., Smith, S., 2001. Segmentation of brain MR images through a hidden 983
Markov random ﬁeld model and the expectation-maximization algorithm. IEEE 984
Trans. Med. Imaging 20, 45–57.
985

T

Lerch, J.P., Evans, A.C., 2005. Cortical thickness analysis examined through power analysis
and a population simulation. NeuroImage 24, 163–173.
Liu, T., Li, H., Wong, K., Tarokh, A., Guo, L., Wong, S.T., 2007. Brain tissue segmentation
based on DTI data. NeuroImage 38, 114–123.
Liu, T., Nie, J., Tarokh, A., Guo, L., Wong, S.T., 2008. Reconstruction of central cortical surface
from brain MRI images: method and application. NeuroImage 40, 991–1002.
Lohmann, G., Preul, C., Hund-Georgiadis, M., 2003. Morphology-based cortical thickness
estimation. Inf. Process. Med. Imaging 18, 89–100.
Luders, E., Thompson, P.M., Narr, K.L., Toga, A.W., Jancke, L., Gaser, C., 2006. A
curvature-based approach to estimate local gyriﬁcation on the cortical surface.
NeuroImage 29, 1224–1230.
MacDonald, D., Kabani, N., Avis, D., Evans, A.C., 2000. Automated 3-D extraction of
inner and outer surfaces of cerebral cortex from MRI. NeuroImage 12, 340–356.
Mietchen, D., Gaser, C., 2009. Computational morphometry for detecting changes inbrain structure due to development, aging, learning, disease and evolution. Front.
Neuroinform. 3, 25.
Miller, M.I., Massie, A.B., Ratnanather, J.T., Botteron, K.N., Csernansky, J.G., 2000. Bayesianconstruction of geometrically based cortical thickness metrics. NeuroImage 12,676–687.
Prastawa, M., Bullitt, E., Ho, S., Gerig, G., 2004. A brain tumor segmentation framework
based on outlier detection. Med. Image Anal. 8, 275–283.
Rocha, K.R., Yezzi Jr., A.J., Prince, J.L., 2007. A hybrid Eulerian–Lagrangian approach for
thickness, correspondence, and gridding of annular tissues. IEEE Trans. Image Process.
16, 636–648.
Rodriguez-Carranza, C.E., Mukherjee, P., Vigneron, D., Barkovich, J., Studholme, C.,
2008. A framework for in vivo quantiﬁcation of regional brain folding in premature
neonates. NeuroImage 41, 462–478.
Rosas, H.D., Salat, D.H., Lee, S.Y., Zaleta, A.K., Pappu, V., Fischl, B., Greve, D., Hevelone, N.,
Hersch, S.M., 2008. Cerebral cortex and the clinical expression of Huntington's
disease: complexity and heterogeneity. Brain 131, 1057–1068.
Rosenfeld, A.P., John, L., 1966. Sequential operations in digital picture processing. J. Assoc.
Comput. Maschiery 13 (4), 471–494.
Sailer, M., Fischl, B., Salat, D., Tempelmann, C., Schonfeld, M.A., Busa, E., Bodammer, N.,
Heinze, H.J., Dale, A., 2003. Focal thinning of the cerebral cortex in multiple sclerosis.
Brain 126, 1734–1744.
Schaer, M., Cuadra, M.B., Tamarit, L., Lazeyras, F., Eliez, S., Thiran, J.P., 2008. A surface-based
approach to quantify local cortical gyriﬁcation. IEEE Trans. Med. Imaging 27, 161–170.
Scott, M.L., Thacker, N.A., 2005. Robust tissue boundary detection for cerebral cortical
thickness estimation. Med. Image Comput. Comput. Assist. Interv. Int. Conf. Med.
Image Comput. Comput. Assist. Interv. 8, 878–885.
Segonne, F., Pacheco, J., Fischl, B., 2007. Geometrically accurate topology-correctionof cortical surfaces using nonseparating loops. IEEE Trans. Med. Imaging 26,
518–529.

C

900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942

R. Dahnke et al. / NeuroImage xxx (2012) xxx–xxx

E

14

U

N

C

O

R

R

E

987

Please cite this article as: Dahnke, R., et al., Cortical thickness and central surface estimation, NeuroImage (2012), http://dx.doi.org/10.1016/
j.neuroimage.2012.09.050

986