Cortical Surface Reconstruction from High-Resolution MR Brain Images

Reconstruction of the cerebral cortex from magnetic resonance (MR) images is an important step in quantitative analysis of the human brain structure, for example, in sulcal morphometry and in studies of cortical thickness. Existing cortical reconstruction approaches are typically optimized for standard resolution (~1 mm) data and are not directly applicable to higher resolution images. A new PDE-based method is presented for the automated cortical reconstruction that is computationally efficient and scales well with grid resolution, and thus is particularly suitable for high-resolution MR images with submillimeter voxel size. The method uses a mathematical model of a field in an inhomogeneous dielectric. This field mapping, similarly to a Laplacian mapping, has nice laminar properties in the cortical layer, and helps to identify the unresolved boundaries between cortical banks in narrow sulci. The pial cortical surface is reconstructed by advection along the field gradient as a geometric deformable model constrained by topology-preserving level set approach. The method's performance is illustrated on exvivo images with 0.25–0.35 mm isotropic voxels. The method is further evaluated by cross-comparison with results of the FreeSurfer software on standard resolution data sets from the OASIS database featuring pairs of repeated scans for 20 healthy young subjects.


Introduction
Cortical reconstruction, the derivation of a computerized representation of the cerebral cortical layer based on threedimensional (3D) magnetic resonance (MR) images of the brain, is an important step in quantitative analysis of the human brain structure, for example, in the analysis of cortical folding patterns, in brain morphometry, and in cortical thickness studies. Cortical surface models typically serve as a reference basis for all further analysis and therefore must be geometrically accurate and topologically correct in order to provide valid and accurate quantitative measures of brain structure [1].
The cerebral cortex, considered at the spatial scale of MR images, is a thin layer of neural tissue, called gray matter (GM), located on the outer side of the white matter (WM), and surrounded by the cerebrospinal fluid (CSF). The cortex has a complex geometry of a highly folded layer with spatially varying curvature and thickness (thickness range 1-5 mm, average ≈2.5 mm, see [1]). The cortical layer on a brain hemisphere can be represented as the inner space between two cortical surfaces (i.e., an inner surface at the WM/GM and an outer or pial surface at the GM/CSF interface, see Figure 1). It is a useful simplification to consider each surface as topologically equivalent to a 3D sphere. In practice, limited spatial resolution of MR images, noise, intensity inhomogeneities, and partial volume effects can all be the sources of geometrical inaccuracies and topological errors in the reconstructed cortical model. In particular, the opposite banks of gray matter in deep sulci are not always resolved as separate and can appear as fused together (Figure 1), leading to invalid models of the cortical layer and propagating errors further into quantitative measurements (e.g., cortical thickness). This may present a particular challenge for an automated reconstruction algorithm, requiring specific means for an automatic detection and correction of topologically and geometrically problematic cases.
Reconstruction of cortical surface models received considerable attention in neuroimaging research. Here, we only briefly overview some state-of-the-art methods; please refer Contours of the inner and pial surface are marked in red and green. Due to partial volume effects and limited resolution, adjacent banks of gray matter in some sulci may appear as fused together, creating either a "bridged" sulcus or an unresolved sulcal fundus (a "buried" sulcus). Note that a "bridged" sulcus creates a topological defect, a handle, which may be corrected by a topology-preserving model, whereas a "buried" sulcus does not change the topology.
to Han et al. [1] and Kim et al. [2] for additional discussion. A suite of algorithms for automated cortical reconstruction is implemented in the popular and freely available FreeSurfer software [3,4]. FreeSurfer includes an algorithm for finding and correcting the topological defects in the initial WM/GM surface [5] and a method to deform the mesh for reconstructing the inner and pial surfaces. The deformable model is constrained by a second-order smoothing term [6] and by a mesh self-intersection prevention routine [3], which both help to resolve the boundaries between adjacent banks in tight sulci. The FreeSurfer automated toolchain is optimized for standard resolution T1-weighted MR images and conforms input data to 1 mm isotropic voxel size, as a rule. This is consistent with the fact that mesh selfintersection detection and prevention is computationally expensive (see [1,6]) and does not scale well with increasing mesh resolution. Xu et al. [7] developed a deformable mesh model for reconstruction of the central cortical surface. The model deforms the topology-corrected initial WM/GM interface by forces derived from a smoothed gradient field [8] that was computed from a GM class membership function. The model does not perform a time-consuming check of mesh self-intersections, which is arguably less critical for finding the central surface, compared to the pial surface. Kim et al. [2] presented a different deformable mesh-based approach for reconstruction of a pial surface, which is called constrained Laplacian anatomic segmentation using proximity, or CLASP. The algorithm computes a Laplacian field mapping between the GM/WM interface and the skeleton of the partial volume classification of the CSF. The Laplacian map is then integrated into the deformable model's objective function, driving mesh vertices into locations with higher values of the Laplacian field. Terms for stretch and self-proximity are included to regularize the deforming mesh and prevent from mesh self-intersection inside sulci. The method by Kim et al. depends on accurate extraction of the CSF skeleton and therefore relies on an elaborate partial volume tissue classification algorithm. However, the accuracy of the Laplacian mapping may be compromised at locations, where the fused GM sulcal banks are not resolved. In addition, the computational cost of the self-proximity term may become prohibitive for high-resolution meshes. Zeng et al. [9] used implicit surfaces in a level set framework for simultaneous reconstruction of the inner and outer cortical surfaces coupled by the minimal and maximal distance constraint. However, this approach did not gain widespread use, because it does not preserve the topology of the evolving surfaces and, in some areas, the distance coupling term may suppress the data attachment term, resulting in geometrical inaccuracies [10]. Han et al. [1] described a method for automated reconstruction of cortical surfaces, called CRUISE, which is built around a geometric deformable model using level sets. To help resolve the cortical banks in sulci, a thin digital separating barrier is constructed using the anatomically consistent enhancement algorithm ACE [1,11], which finds a skeleton of the weighted distance function computed from the Eikonal equation with a speed function modulated by the CSF class membership. At the core of the CRUISE method is a topology-preserving geometric deformable surface model, TGDM [1,11,12], which models the evolution of a level set function under the influence of signed pressure forces computed from tissue class membership values and curvature forces defined by the surface geometry. The central surface of the cortex is reconstructed by a TGDM with GGVF advection forces similar to those in Xu et al. [7].
We present a method, henceforth, designated dielectric layer field mapping, or DELFMAP, for the automated reconstruction of the cortical compartment from MR images, which is based on several partial differential equation (PDE) modeling stages. Our method is inspired by the work of Han et al. and uses a similar level set framework, but introduces a different perspective, consolidating all algorithmic stages around the key mathematical model of a potential field in an inhomogeneous dielectric medium. Our method scales well with image resolution and has an advantage over other existing methods in reconstruction from high-resolution MR images with submillimeter voxel sizes, because (1) in contrast to deformable mesh models in FreeSurfer or CLASP, it avoids the computational cost of testing for mesh self-intersection and self-proximity; (2) similarly to CRUISE, it uses an efficient narrow-band algorithm for the level set evolution; (3) in contrast to CRUISE that requires solving a system of three second-order PDEs in GGVF, our method solves just one second-order PDE and does not need an intermediate step of reconstructing a central cortical surface.
Preliminary results of this work were presented in two conference publications [13,14]. This report expands on the methodology and experimental results and adds a validation study that performs cross-comparison of our method's cortical reconstruction results with those obtained using FreeSurfer [3,4] on standard resolution data for 20 healthy young subjects (test-retest repeated scans) from the OASIS database [15].

Methods
The DELFMAP method proceeds as follows. A potential field is computed using the mathematical model of an electric field in an inhomogeneous dielectric medium, where International Journal of Biomedical Imaging 3 the segmented WM poses as a charged conductive object and the classified GM poses as an inhomogeneous dielectric layer with permittivity proportional to GM class probability values. This electrostatic model serves the purpose of concentrating the flux of the mapping flow in a layer of voxels classified as GM and helps to identify the separating barriers between cortical banks in sulci, where the mapping flow collides. Correspondence trajectories following the lines of the potential field and geodesic distances from WM boundary are determined using PDEs, and a digital skeleton of the sulcal medial surface separating GM sulcal banks is derived by finding collisions in the correspondence trajectories and shocks in the distance field. The computed electric field retains the desired laminar properties of the Laplacian mapping in the bulk of the cortical layer and is used as the potential flow that maps the inner surface to the outer. The outer (pial) cortical surface is reconstructed using a geometric deformable model level set framework [16] with an advection along the gradient of the potential field, which is constrained by the identified skeleton of the sulcal medial surfaces and (optionally) by a maximal distance/proximity constraint.
2.1. Image Processing Chain. DELFMAP takes as input a set of volumetric images containing WM and GM tissue class probability/membership functions and a refined WM model, supplied either as a topology-corrected WM binary segmentation or as a WM/GM interface level set function. The overall chain of general image processing steps is outlined as follows ( Figure 2) (1) A T1-weighted volumetric MR image is (optionally) aligned with the stereotaxic coordinate system, interpolated to isotropic voxel size, and is preprocessed with a brain-peeling algorithm that derives a mask of voxels related to the cerebral tissues only. (2) The brain image is corrected for intensity inhomogeneities and is classified into WM,GM,CSF/background probability images. (3) A raw WM binary segmentation is derived from the class probability images (by thresholding or a maximum-probability rule), and brain stem and cerebellum are (optionally) removed from the WM segmentation. (4) A topology-corrected WM volume is obtained from the raw WM binary segmentation by an automated algorithm or by manual editing, or a combination of both. (5) DELFMAP uses the output of step 2 and step 4 to reconstruct the inner and outer cortical surfaces. We note that steps 1-4 are common to many brain MR image processing workflows, therefore DELFMAP can be easily integrated with a wide variety of toolchains. More specifically, we used processing steps described in Yang and Kruggel [17] in our experiments with 3-Tesla in-vivo images, and we applied algorithms described in Kruggel et al. [18] for the analysis of exvivo high-resolution images. In step 4, for exvivo images, we used manual editing for filling ventricles and correcting large topological defects, and we applied a topological region-growing algorithm similar to the one in Kriegeskorte and Goebel [19] to obtain a genus zero WM binary object. In cross-validation with FreeSurfer on the OASIS data sets, we used the FreeSurfer's processing toolchain for the initial steps that are common between the two methods (i.e., steps 1-4 that lead to a topologically-corrected WM segmentation); therefore, the cross-method comparison of cortical reconstructions is not confounded by differences in preprocessing approaches. Finally, we emphasize that, in all our experiments involving DELFMAP, the tissue classification was performed by a modified version (see [18]) of the adaptive fuzzy clustering algorithm [20] augmented with a spatial regularization term [1]; this also applies to GM and WM tissue classification that was used by DELFMAP in cross-validation study on the OASIS data sets.

Inner Cortical Surface.
The inner cortical surface is reconstructed by a deformable model (Figure 2, step 5.0) that smooths the initial WM/GM interface, which is determined by the corrected WM segmentation. For this purpose, we use a topology-preserving geometric deformable model (similar to [12]), which is described in detail in Section 2.6. For smoothing, we typically run 2-3 iterations of the deformable model with the mean curvature term only. We will denote the "inside" region of the level set function representing the inner cortical surface by Ω w .

Electric Field Model.
A potential field is found as a solution to the PDE modeling an electric field around a charged conductive object (WM) insulated by a dielectric layer (GM) having spatially inhomogeneous electric permittivity, which is set proportional to GM tissue class probability ( Figure 2, step 5.1). In such a model, the flux of the electric field is confined in regions of higher permittivity, that is, where GM class probability is higher; therefore, trajectories following the lines of the electric field trace through the GM layer before exiting into the background space. Thus, the flux of the mapping flow is concentrated in a layer of voxels classified as GM. Let Ω denote the 3D image domain with the boundary Γ(Ω). We will denote WM and GM tissue class probability images by P w ( r ) and P g ( r ) ( r ∈ Ω), where r = (x, y, z) is a 3D point. Let ϕ( r ) denote a potential field, a scalar function defined over Ω. Let ε( r ) denote another scalar function, called permittivity and computed from class probabilities as follows: where ε max is the maximum permittivity of the insulating layer (ε max should be 1 in order to emphasize the inhomogeneity of the dielectric layer; ε max = 100 was used, and ε max = 1000 was tested with similar results). Thus, permittivity is close to ε max when WM and/or GM class probabilities are high and is close to 1 when they are low. Note that the WM probability is included above only to ensure a proper transition of the field near the WM/GM interface, where some border voxels can be classified with low GM but high WM probability, for example, because a smoothed interface can slightly deviate from the initial WM segmentation. The inclusion of WM probability is therefore limited by the constraint field C d , which is computed by thresholding of the WM chamfer distance transform D cmf as C d = {1 if D cmf < d min , 0 otherwise}, where the distance threshold d min can be set at the lower bound on cortical thickness (≈1 mm), just enough to ensure a "high-permittivity" transition via (3) Preprocessing, generic steps 1-4 Cortical reconstruction, steps 5.0-5.3 boundary WM voxels to the layer of GM voxels. The potential field is found as a solution of Maxwell's equation for an electric field inside inhomogeneous dielectric medium in the absence of free charges: Equation (2) assumes that the dielectric medium has linear and isotropic properties; therefore, ε is a scalar, not a tensor. Boundary conditions are specified as ϕ( r ∈ Ω w ) = 1 and ϕ( r ∈ Γ(Ω)) = 0, that is, the potential is set to one in the WM core and is set to zero on the boundary of the image volume. The solution of the PDE ϕ( r ∈ Ω \ Ω w ) can be obtained as a steady-state solution (∂ϕ/∂t → 0) of a corresponding nonstationary equation: Equation (3) can also be viewed as describing the diffusion in inhomogeneous medium, where ε( r ) is a spatially varying but stationary diffusion coefficient and ϕ( r, t) is the concentration of the diffusing substance. This allows for a different physical interpretation of the model: we seek a steady-state spatial distribution of "particles" diffusing from International Journal of Biomedical Imaging 5 WM source into the medium with a diffusivity proportional to the GM class probability. Qualitatively, it is expected that "particles" would diffuse more freely in GM; therefore, the lines of the gradient field ∇ϕ would tend to concentrate in the GM compartment. Equation (3) can be discretized and solved iteratively as described by [21], for example, using the Jacobi method [22].

Distance Field and Correspondence Functions.
Lines of the potential field ϕ are defined as a family of curves that are at each point tangent to the gradient ∇ϕ. Let d( s, r ) denote the length of a line segment originating at some point in WM boundary s ∈ Γ(Ω w ) and ending in point r ∈ Ω \ Ω w . If, for any point r, there is one and only one streamline passing through it, then d( r ) defines a distance field. It is possible to compute the distance field by integrating trajectories explicitly in a Lagrangian framework. Alternatively, using the method described in Yezzi and Prince [23], the distance field can be found as a solution of a PDE in an Eulerian framework on a fixed grid. We note that ∇ϕ/ ∇ϕ is the unit tangent field of the potential field ϕ. Then, it can be shown that the distance field d must satisfy the following PDE: with the boundary condition d( r ∈ Γ(Ω w )) = 0. Correspondences along streamline trajectories can be computed in a similar way. More specifically, let denote a vector of correspondence functions, which establishes a correspondence between a point in the field domain r ∈ Ω \ Ω w and a "source" point in the WM boundary ψ ∈ Γ(Ω w ). These correspondence functions ψ i can be found as solutions of the following PDE (see [24]): with boundary conditions The first-order PDEs (4) and (5) can be solved using the numerical implementation described by Yezzi and Prince [23]. In principle, finite spatial discretization may violate the one-to-one correspondence property of the flow by clamping several streamline paths into one point on a grid, so the solutions d( r ) and ψ( r ) may experience numerical convergence problems in some grid locations. In practice, we found that such problematic points are very sparse and do not impede numerical convergence in the computational domain at large. These points are usually detected among other "shocks" in the distance field by a skeletonization method (Figure 2, step 5.2), which is described next.

Skeleton of the Sulcal Medial
Surface. Inside sulci, streamlines originating from opposite cortical banks collide (due to spatial discretization), which results into shocks in the distance field and into "discontinuities" in the correspondence functions. Shocks or singularities of a distance field d are defined as a set of points, where spatial derivatives of the field are not well-behaved, that is, the gradient ∇d is not well defined. Such shocks appear as discontinuities or sinks in the field. Note that even though the potential field in our model should be, in theory, free from the sinks (because there are no free charges), they may appear in the distance field due to spatial discretization. Let S ⊂ Ω \ Ω w , called a skeleton of the distance field, denote a set of points on a grid, where shocks are detected by a numerical procedure. Such numerical procedure can be based on finite difference approximations to ∇d, as described by Han et al. [1]. The observation is that a centered finite difference numerical scheme will produce values of ∇d that are significantly lower than 1 on the shock points and are close to unity elsewhere. Then, the skeleton can be detected as where d min is a minimum distance parameter set at the lower bound on cortical thickness and T is a specified threshold value (T < 1; values d min = 1 mm and T = 0.8 can be used, similarly to ACE in [1]). We found that the skeleton can be robustly detected by a novel algorithm based on the analysis of the correspondence function [14]. Recall that ψ( r 0 ) is a vector with coordinates of the streamline's source point at WM boundary. A streamline collision can be detected if, in the neighborhood of r 0 , there are correspondences to source points that are "distant" between themselves. More formally, the skeleton can be determined as We used D min = 4 voxels and 6 adjacent points N 6 ( r ) in our computations.

Geometric Deformable Model.
The geometric deformable model uses an implicit representation of a surface, embedding it into a level set function φ( r, t)( r ∈ Ω). The evolving interface is represented by the zero-level set Φ(t) = { r | φ( r, t) = 0} (see [16]), and it can be retrieved with subvoxel resolution by an isosurface algorithm (e.g., marching cubes). In our model, evolution of the level set function is described by the following PDE that has an advection and a mean curvature term: where V is the advection velocity vector field, κ is the mean curvature, and w κ are weights of the respective terms (w α , w κ ≥ 0). The mean curvature of the interface embedded in the level set function is [16] 6 International Journal of Biomedical Imaging The stopping/reversal factor The distance-constraining factor γ(d) where the subscripts x, y, z denote partial derivatives. The advection velocity vector field V ( r ) is derived from the gradient of the potential ϕ or distance field d: where β( r ) is a stopping/direction-reversal factor computed from the GM/WM class probabilities. For example, this factor can have a form of a logistic function: where K is the constant controlling the steepness of the slope of the sigmoid curve and P 0 is the GM class probability threshold value that determines the "set-point" of the deformable model. Figure 3 illustrates how the factor β depends on GM and WM probability P gw . In our experiments, a moderately steep sigmoid curve with K = 40 and the threshold P 0 = 0.8 were used. For spatial regularization, the combined GM and WM class probability P gw ( r 0 ) can be calculated as a weighted sum over the (closed) neighborhood of the point r 0 : where w i are the neighborhood weights (e.g., w i = 0.5/n, where n = 18 or 26, and for the central point w 0 = 0.5), and the skeleton of the sulcal medial surfaces S is used for masking of class probability values in separating barriers. As an option, the stopping factor β in (8) can be modified to include the distance-constraining factor: where the sign function is an "OR" combination of two signs: and the distance-constraining factor γ can also have a form of a logistic function: In (13), d max is a parameter constraining the maximum distance of advection along the streamlines of the gradient field (i.e., a proximity constraint that limits the thickness of the reconstructed cortical layer). We used d max = 6 mm (see Figure 3) in the reported cortical reconstructions, that is, the maximum distance constraint was set above the anatomically plausible upper bound on cortical thickness and therefore was affecting only the artefactual or noncortical gray matter areas.
Our numerical implementation for solving the level set (6) is based on the narrow-band algorithm [12,16,25]. The initial level set function is computed as a signed-distance function (SDF) of the initial interface in the corrected WM image using the fast marching method (FMM, [16,26]). By standard convention, "inside" points are represented by negative values of the SDF. During the evolution, the level set function φ( r, t) is maintained close to the SDF by periodic  reinitialization with the FMM. The advection term in (6) is discretized based on the upwind differencing scheme (for details, see [16]), and the curvature term is discretized along the lines of (7) using the central differencing scheme [22]. A pseudocode outlining the narrow-band algorithm is described elsewhere (e.g., in [12,25]). In Algorithm 1 pseudocode we focus on the core part that deals with the timestep update of the level set function. The update algorithm is novel in the way it uses the skeleton of the sulcal medial surface to create barriers for the evolving interface. In addition, the algorithm has a built-in rule preserving the digital topology of the deformed model [1,12] that is based on the concept of simple points [27] (function IsSimple() in Algorithm 1, see details in [13]), which guarantees that the deformed surface retains the same topology as the initial WM/GM surface.
As already mentioned, the inner cortical surface is reconstructed by a few iterations of the model with the curvature term only (w α = 0, w κ = 1) (Figure 2, step 5.0). In step 5.3 of Figure 2, the outer cortical surface is first reconstructed by a model using the advection term only (w α = 1, w κ = 0) until convergence (i.e., until the relative amount of change in the SDF per iteration becomes small, for example, lower than 10 −4 ) or for a specified number of time steps and then smoothed by a few iterations with the curvature term, similarly to the inner surface.

Experiments and Results
Our algorithm was implemented in C++ in the Linux environment and ran on a PC with 2.5 GHz AMD-64 CPU and 4 GB RAM, unless otherwise noted. The algorithm's performance was evaluated on simulated test cases with a simplified geometry of a sulcus, on simulated MRI datasets, on standard resolution T1-weighted MR images of human brains, and on high-resolution (sub-mm) MR images of extracted brain hemispheres.

Simulated Data.
The first test case is intended to illustrate the effect of the inhomogeneous dielectric model used in DELFMAP and shows the difference between the field produced with a nonuniform permittivity and the field computed with the uniform permittivity (ε = 1, the Laplacian field). Test images simulate a simplified 3D geometry of a sulcal fold and contain two WM stalks separated by the sulcal space (with a curvature radius of 10 mm); the WM is covered by a layer of GM having unequal thickness at the opposing banks and a smoothly varying thickness at the fundus (Figure 4(a)). Figure 4 shows the lines of the Laplacian field (Figure 4(b)) and the lines (Figure 4(d)) and isocontours (Figure 4(c)) of the field in the DELFMAP model. It can be seen that the "ridge" (where the field lines concentrate and the isocontours converge) of the DELFMAP field is close to the sulcal center line, whereas the "ridge" of the Laplacian field is at the geometric center.
The second test case demonstrates how the model resolves the barrier separating the two opposing cortical banks inside a sulcus. Test images simulate a fully resolved sulcus (with two banks fully separated by background), a sulcus with an unresolved fundus, and a sulcus with two banks bridged by unresolved voxels (the top row in Figure 5: left, middle, and right, resp.). The middle row in Figure 5 shows the cross-section of the sulcal medial surface (white lines) that was identified by the DELFMAP method. It can be seen that the method is capable of reconstructing the boundary surface separating the two cortical banks and finds a geometrically plausible solution in incompletely resolved cases. Side-by-side comparison of the results of our method and those of ACE (the bottom row in Figure 5) shows that skeletons produced by DELFMAP have a more regular structure, whereas ACE skeletons can have small extraneous branches and discontinuities. Our method does not produce spurious detections very close to WM and thus does not require a minimum distance cut-off parameter, which is needed in ACE. In addition, our method is more robust with respect to noise (see [14]): skeletons produced by DELFMAP show very little degradation even at the highest noise level, while ACE skeletons are significantly affected by strong levels of noise. Cortical reconstruction results for simulated brain phantom MR images [28] showed good reproducibility across various levels of simulated Gaussian-distributed noise and intensity inhomogeneity (see [13,14]).

High-Resolution MR Images.
Our method's performance is illustrated by results for high-resolution exvivo images, where, contrary to FreeSurfer, our method does not need to conform images to standard 1 mm isotropic voxel size.
The algorithm was evaluated on three high-resolution (0.25-0.35 mm isotropic voxel size) images of explanted brain left hemispheres. DELFMAP reconstruction at 0.35 mm resolution took 67 min on a PC with 2.5 GHz AMD-64 CPU and 4 GB RAM. We tried to process the same 0.35 mm data with the recently released CRUISE plugin for MIPAV [29] on a cluster node with four Opteron 285 2.6 GHz cores and 32 GB RAM. Reconstruction of the inner surface took 28 min using 4.9 GB RAM, computation of GGVF took 32 min using 3.5 GB RAM, while reconstruction of the central and pial surfaces took 49 and 52 min using 5.3 and 5.1 GB, respectively, but did not produce adequate results with the default settings. DELFMAP computations at 0.25 mm resolution required 4.7 GB RAM and were successfully completed after 3 h 20 min. Examples of the reconstructed cortical surfaces overlaid on orthogonal cross-sections of a high-resolution MR image are shown in Figure 6.  convexity maps, gyral crowns appear in blue color and sulcal fundi appear in yellow-orange. Thickness and convexity maps demonstrate noticeable correlation (Pearson's correlation coefficient computed over the entire surface mesh is 0.24, 0.22, and 0.28 for the three brain samples shown, that is, significantly different from zero at the 0.05 level), which is in good agreement with the known anatomical fact that cerebral cortex is generally thicker on gyral crowns and thinner in sulcal depths.

Cross-Validation with FreeSurfer: Test-Retest Precision.
Our method was validated by cross-comparison of cortical reconstruction results with those obtained using FreeSurfer. Standard resolution images for 20 right-handed healthy young subjects (age 19-34, average 23; 8 males/12 females) were selected from the cross-sectional OASIS database [15]. For each subject, data are available from two scan sessions (test and retest) separated by a short delay (1-89, average 21 days), with four T1-weighted standard resolution images acquired per session. This relatively short delay between two consecutive scan sessions makes data sets suitable for the assessment of test-retest reproducibility (i.e., precision) of the analysis by comparing measurements between scan sessions.
First, we analyzed data sets using the default automated pipeline in FreeSurfer and obtained 40 cortical reconstructions (two per subject), each including a pial and a white surface mesh. Next, we exported images of extracted brains (without any intensity normalization/correction) and corrected WM segmentations from FreeSurfer, ran our tissue classification algorithm on images of extracted brains, and used these results in the DELFMAP toolchain to obtain another set of 40 cortical reconstructions. For a subvoxel resolution of a digital skeleton, solutions of PDE in (3)-(6) were computed on a grid with half-voxel spacing. Implicit level set surfaces were tessellated using connectivity-consistent marching cubes algorithm [12], and triangular meshes were simplified down to 300,000 faces by a topology-preserving variant of the mesh simplification method [30]. DELFMAP processing took approximately 30 min per brain hemisphere (at half-voxel 0.5 mm res. grid) and was twice faster than FreeSurfer's deformable model step (mris make surfaces program, took ≈70 min at 1 mm res.). FreeSurfer computes cortical thickness at each vertex as the average of the closestpoint distance (Figure 8(a)) measured between the surfaces both ways using linked vertices [6]. Since vertices on pial and white surfaces are not linked in DELFMAP, which is not based on a deformable mesh model, for the crossmethod comparison, we recomputed cortical thickness using an orthogonal projection distance measure [31] (Figure 8(b) and the Appendix) that is robust and universally applicable to results from both methods. We verified that the two cortical thickness measures were in close agreement on all 40 reconstructions obtained with FreeSurfer.
International Journal of Biomedical Imaging The geometric precision or test-retest reproducibility of cortical reconstruction was evaluated independently for FreeSurfer and for DELFMAP as follows. For each subject, test and retest MR images (averages of 4 aligned scans from the first and the second session, resp.) were rigidly registered to each other using FSL FLIRT [32]. The obtained rigid transformation was applied to the first set of surface meshes, aligning the test surfaces to the retest ones. Next, signed and absolute distances (the Appendix, (A.1) and (A.2)) were measured between aligned test and retest white/pial surface meshes, and surface-wise mean and standard deviation were computed, as well as the group-wise statistics. In addition, we evaluated the test-retest precision of cortical thickness measured with FreeSurfer and with our method using the standard methodology described in the cortical thickness reproducibility study in Han et al. [33], which consists of the following four steps: (1) rigid registration of two repeated scans of each subject; (2) computation of a thickness difference map for each subject (on the first surface, using point-correspondences established according to closest Euclidean distance in registered space); (3) resampling the thickness difference map to a common template (e.g., any subject surface or the FreeSurfer's average template); (4) computing the group-wise mean and standard deviation of the differences at every vertex of the template mesh. Resampling to a common template relies on FreeSurfer's intersubject registration by nonlinear surface morphing [34].
Results of both methods, the absolute distance measure AD mean and AD stdev between test and retest cortical surfaces (the Appendix, (A.3)), per subject hemisphere, were compared statistically using a Wilcoxon signed rank test, and results are reported as P values. For FreeSurfer WM surfaces, reproducibility is characterized by mean absolute error 0.19(Δ0.06) mm (where the Δ value in parentheses indicates a statistical spread for the group, equal to two stdev). For DELFMAP WM surfaces, mean absolute error is 0.24(Δ0.06) mm (P = 9.5 × 10 −5 ). For DELFMAP pial surfaces, reproducibility is characterized by a mean absolute error 0.24/0.25(Δ0.04) mm (L/R) that is similar in FreeSurfer (L: P = 0.37, R: P = 0.16, see details in Table 1). The standard deviation of the absolute distance AD stdev is much lower in DELFMAP than in FreeSurfer (L: P = 8.2 × 10 −5 , R: P = 3.2 × 10 −4 ) which can be interpreted as a "tighter" reconstruction of pial surfaces in DELFMAP. Table 1 summarizes the statistics of the test-retest analysis. The mean absolute difference of the cortical thickness is similar in both methods (L: P = 0.10, R: P = 0.28), but the corresponding standard deviation is again much smaller in DELFMAP than in FreeSurfer (L: P = 1.9 × 10 −6 , R: P = 1.0 × 10 −4 ). To summarize, test-retest precision of cortical thickness measurement is similar in DELFMAP and FreeSurfer in terms of the mean error, which is close to a quarter of the voxel size, but is "tighter" in DELFMAP in terms of surfacewise variance in absolute differences.

Cross-Validation with FreeSurfer: Intermethod Accuracy.
The geometric accuracy of our method was evaluated by cross-comparison with FreeSurfer as follows. For each cortical reconstruction (two per subject), white (W) and pial (G) surfaces (Wf, Gf) were exported from FreeSurfer and a cortical thickness map A GfWf (A.2) was computed on pial surface. Next, maps of intermethod geometric differences (D WfWd , D WdWf , D GfGd , D GdGf ) were computed as signed distances (A.1) between white or pial surfaces reconstructed with FreeSurfer and DELFMAP (Wd, Gd). On these geometric-difference maps (40 sets, four maps per set), surface-wise statistics D mean , D stdev , AD mean , and AD stdev (A.3) were computed. In addition, maps of intermethod thickness differences were built using the cortical thickness reproducibility analysis steps 2-4 [33] as described in the previous section, except for using two pial surfaces from both methods in step 2 (we emphasize that for both FreeSurfer and DELFMAP, the compared thickness maps were measured by the same method, that is, as A GW ). The 40 individual maps were resampled to a common template and averaged into group-wise maps of mean difference and standard deviation. The group-wise maps of intermethod cortical thickness measurement differences allow to assess and visualize any Table 1: Precision analysis: summary of the group-average statistics for the signed distance (SD) and absolute distance (AD) measure (in mm) between test and retest surfaces (surface: DF-DELFMAP, FS-FreeSurfer; L/R: left/right hemisphere; mean: a group average of a surface-wise mean of the distance; stdev: a group average of a surface-wise stdev of the distance; "> X mm (%)": (group-average) percentage of surface points where AD was greater than X mm; values in parentheses indicate the statistical spread within the group, measured by the group-wise stdev).   regional patterns of agreement/disagreement between the two methods. The intermethod geometric accuracy analysis statistics is summarized in Table 2 (averaged over 40 image sets, two per subject). It can be seen from the mean signed distance SD mean that, on average, DELFMAP has a very small outward bias in pial surfaces (−0.08/−0.07(Δ0.08) mm, L/R; negative sign means FreeSurfer' surface is "inside" w.r.t. DELFMAP' surface). The intermethod accuracy can be characterized by the mean absolute distance AD mean (0.40/0.42(Δ0.04) mm, L/R), which is less than a half of the voxel size. The share of pial surface vertices where the AD was larger than 1 mm is less than 10%; less than 1% of pial vertices had an AD larger than 2 mm. The intermethod accuracy analysis of cortical thickness measurements, summarized in Table 3, is in good agreement with the above observations. On average, there is a small bias towards thicker values in DELFMAP (mean signed difference: 0.12/0.11(Δ0.16) mm, L/R; positive sign here means that DELFMAP-measured thickness is larger w.r.t. FreeSurfer). The intermethod accuracy, characterized by the mean absolute difference (0.35/0.34(Δ0.06) mm, L/R), is less than a half of the voxel size. The share of pial surface vertices where the absolute difference between thickness measurements was larger than 1 mm is less than 6%, and less than 1% of pial vertices had an absolute difference larger than 2 mm. An example comparing DELFMAP and FreeSurfer pial surface reconstructions side-by-side, for one subject, is shown in Figure 9 (colored with cortical thickness; see colorbar for color map and range of values). Overall, a good correspondence is visible, but some patterns of thickness difference are noticeable: (1) for FreeSurfer, thickness is larger (indicated as yellow) in the superior region of the frontal lobe and in some temporal regions (lateral view); (2) for DELFMAP, thickness is larger (indicated as orange) in the inferior occipitotemporal region (medial view, where the cerebellum is found); (3) for FreeSurfer, thickness is smaller (indicated as blue) in the medial orbitofrontal cortex (mOFC) region (medial view). These differences can be attributed and traced to the following segmentation trends in either of the two methods: (1) oversegmentation, by FreeSurfer, into meningeal space in superior frontal region and in temporal region (see Figure 10); (2) oversegmentation, by DELFMAP, into cerebellar gray matter in the inferior occipitotemporal region; (3) too conservative segmentation, 14 International Journal of Biomedical Imaging by FreeSurfer, in the mOFC region (too thin, less than 1.5 mm). Regional patterns of intermethod geometric differences in pial cortical reconstructions are visible on group-average maps of geometric ( Figure 11) and cortical thickness differences (Figure 12), where the above outlined three trends are also noticeable.

Discussion
We presented a novel PDE-based approach for reconstructing the cerebral cortex from MR images. We developed an accurate and scalable method that works on MR images with a high spatial resolution. Because high-resolution MRI begins to attract considerable attention in brain imaging research, a method that readily scales with imaging resolution is highly valuable. This scalability is achieved by using an implicit deformable surface model in a fast marching framework guided by a novel, computationally efficient model using potential field mapping. Our method requires much lower computational resources and has a much faster computation times than conventional methods. These demonstrated advantages come not only from an efficient practical implementation, but also from the design of our algorithms. For instance, other existing approaches that are based on deformable mesh models incur a significant computational cost associated with the mesh self-intersection (e.g., FreeSurfer) or mesh self-proximity (CLASP) term, which does not scale linearly with increasing mesh resolution. Although the computational cost of the straightforward mesh self-proximity term [2], which is quadratic O(N 2 /2) on the number of faces N, is significantly reduced in a mesh selfintersection prevention algorithm utilizing a spatial cache [3], it nevertheless remains supralinear. Similarly, the cost of another known efficient algorithm for mesh self-intersection detection, which is based on intersection of bounding boxes, is O(Nlog 3 2 N) [35]. In contrast to this, the computational complexity O(Nk) of the narrow-band level set algorithm used in our method (and in CRUISE) is linear with respect to the size of the interface N (k is the width of the narrow band). This difference between a linear and a quadratic or supralinear algorithmic complexity, which can be tolerated when dealing with standard resolution images and meshes, becomes quite large at high resolutions. As to the comparison with the available CRUISE MIPAV software, our method's dramatic gain in speed is most likely due to differences in implementation but, at least in part, can be attributed to a smaller algorithmic cost of our method (e.g., solving one second-order PDE in DELFMAP versus a system of three second-order PDEs in GGVF, and not using an intermediate step of reconstructing a central cortical surface).
Although some algorithmic building blocks of our method were previously known to the medical image processing community (e.g., [1,21,23]), the central aspect of our method, that is, the use of the model of the potential field in the inhomogeneous dielectric layer introduced here, is novel and has attractive advantages. The novelty of our method is also in the newly introduced skeletonization algorithm that is based on the analysis of correspondence trajectories and in several novel aspects of the geometric deformable model (e.g., the constraint of the evolution by the medial surfaces, the maximal distance constraint of the advection, and the novel form of the advection stopping/directionreversal factor β and the distance-constraining factor γ). We note that most of the design parameters introduced in Section 2 remain fixed, and the method is sensitive only to two settings, which can be easily tuned: the GM probability threshold P 0 (a "set-point") and the maximal distance d max (which has strong influence only if set below the upper bound on cortical thickness).
The results from three high-resolution data sets demonstrate that the method is capable of reconstructing the outer Figure 10: Contours of reconstructed cortical surfaces overlaid on the axial (left) and coronal (right) slice (red: inner surface; green: DELFMAP pial surface; yellow: FreeSurfer pial surface; note that the yellow contour appears jagged because it is displayed from FreeSurfer's volumetric signed distance function sampled at 1 mm grid, whereas red and green contours are from level set functions sampled at a finer resolution; left and right images are not to scale). On the left image at the cross-line cursor position (superior frontal region), the yellow contour of FreeSurfer's pial reconstruction oversegments into meningeal space, and a similar trend is noticeable next to cursor on the right image (temporal region). cortical boundary with good geometric precision and accuracy, while guaranteeing the preservation of the initial surface topology. The method's performance is illustrated on synthetic images and on standard resolution MR brain images, where it compares favorably to existing methods in both quality and speed.
The precision and accuracy of our method was assessed by cross-validation in standard resolution datasets with the widely accepted approach implemented in the available FreeSurfer software. Using a database of consecutive examinations in healthy subjects, the precision of both methods was evaluated using pointwise geometric distances of reconstructed surfaces and differences in cortical thickness. Both methods are similar in terms of the mean absolute error in position and mean absolute error in cortical thickness. However, DELFMAP has a much lower variance than FreeSurfer. In a second study, we evaluated the accuracy of our method by quantifying the intermethod reproducibility of reconstructed cortical surfaces, measured by pointwise geometric distances and differences in cortical thickness measurement between the two methods. Results demonstrate that the accuracy of our method, using FreeSurfer as a reference, is better than half of a mm in terms of both mean absolute error in geometric position and mean absolute error in measured cortical thickness. Group-average analysis of the spatial distribution of geometric and thickness differences between the two methods reveals some surface regions, where one of the two methods has a tendency to systematically overor undersegment the cortical ribbon, resulting in patterns of small (subvoxel) but measurable differences. Thus, crosscomparison of the two methods allows detection of existing regional patterns in intermethod differences, benefiting the study of accuracy of both approaches and highlighting some potentially problematic areas for further improvement of both methods.  Figure 11: Group-average maps of intermethod (DELFMAP-FreeSurfer) geometric differences in pial surface reconstructions, resampled to FreeSurfer's average template (left column: signed distance mean, colorbar range ±1 mm, negative/positive values mean FreeSurfer' surface is inside/outside of DELFMAP' surface, resp.; middle column: absolute distance mean, colorbar range 0-1 mm; right column: absolute distance stdev., colorbar range 0-1 mm; rows 1-4: lateral/medial surface of left/right hemisphere, resp.).  Figure 12: Group-average maps of intermethod (DELFMAP-Freesurfer) cortical thickness differences, resampled to FreeSurfer's average template (left column: signed difference mean, colorbar range ±1 mm, negative/positive values mean thickness measured with DELFMAP is smaller/larger than measured with FreeSurfer, resp.; middle column: absolute difference mean, colorbar range 0-1 mm; right column: absolute difference stdev., colorbar range 0-1 mm; rows 1-4: lateral/medial surface of left/right hemisphere, resp.).

Distance Measure between Two Surfaces
The orthogonal projection method [31] was adapted to define a measure of geometric distance between two meshes that was used throughout our validation study. We note that a similar approach was proposed in Tosun et al. [36] for accuracy and precision analysis of cortical surface reconstructions. The signed distance between two triangulated meshes M 1 = {V 1 , F 1 } N1 and M 2 = {V 2 , F 2 } N2 was measured as: where P F2 ( v 1,i ) is the closest orthogonal projection operator projecting a vertex v 1,i from the first mesh onto one of the triangles F 2 in the second mesh, along the normal to that triangle (Figure 8(b)). The sign of the distance measure is determined by the innerproduct of the projection difference vector p 12,i with the first surface outward normal n 1,i at vertex v 1,i ; thus, a positive/negative value signifies that the second surface is outside/inside of the first surface, respectively. For the signed distance, mean and standard deviation (stdev) are computed on d 12,i values. We define the absolute distance measure as: In this notation, the cortical thickness measure of Kruggel and von Cramon [31] is defined as the absolute distance from GM to WM mesh A GW ; this orthogonal projection measure should not be confused with the distance along surface normal [37], which was shown to be less reliable compared to other distance measures [38]. For the absolute distance, twoway mean and standard deviation are computed (see also [36]) as: