Abstract
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 PDEbased method is presented for the automated cortical reconstruction that is computationally efficient and scales well with grid resolution, and thus is particularly suitable for highresolution 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 topologypreserving 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 crosscomparison 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.
1. 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 stateoftheart methods; please refer 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 secondorder smoothing term [6] and by a mesh selfintersection 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 T1weighted 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 topologycorrected 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 timeconsuming check of mesh selfintersections, which is arguably less critical for finding the central surface, compared to the pial surface. Kim et al. [2] presented a different deformable meshbased 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 selfproximity are included to regularize the deforming mesh and prevent from mesh selfintersection 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 selfproximity term may become prohibitive for highresolution 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 topologypreserving 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 highresolution 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 selfintersection and selfproximity; (2) similarly to CRUISE, it uses an efficient narrowband algorithm for the level set evolution; (3) in contrast to CRUISE that requires solving a system of three secondorder PDEs in GGVF, our method solves just one secondorder 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 crosscomparison of our method’s cortical reconstruction results with those obtained using FreeSurfer [3, 4] on standard resolution data for 20 healthy young subjects (testretest repeated scans) from the OASIS database [15].
2. 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 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 topologycorrected 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 T1weighted volumetric MR image is (optionally) aligned with the stereotaxic coordinate system, interpolated to isotropic voxel size, and is preprocessed with a brainpeeling 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 maximumprobability rule), and brain stem and cerebellum are (optionally) removed from the WM segmentation. (4) A topologycorrected 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 3Tesla invivo images, and we applied algorithms described in Kruggel et al. [18] for the analysis of exvivo highresolution images. In step 4, for exvivo images, we used manual editing for filling ventricles and correcting large topological defects, and we applied a topological regiongrowing algorithm similar to the one in Kriegeskorte and Goebel [19] to obtain a genus zero WM binary object. In crossvalidation 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 topologicallycorrected WM segmentation); therefore, the crossmethod 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 crossvalidation study on the OASIS data sets.
2.2. 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 topologypreserving geometric deformable model (similar to [12]), which is described in detail in Section 2.6. For smoothing, we typically run 23 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 .
2.3. 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 and (), where is a 3D point. Let denote a potential field, a scalar function defined over . Let denote another scalar function, called permittivity and computed from class probabilities as follows: where is the maximum permittivity of the insulating layer ( should be in order to emphasize the inhomogeneity of the dielectric layer; was used, and was tested with similar results). Thus, permittivity is close to 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 , which is computed by thresholding of the WM chamfer distance transform as , where the distance threshold can be set at the lower bound on cortical thickness (≈1 mm), just enough to ensure a “highpermittivity” transition via 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 and , 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 can be obtained as a steadystate solution of a corresponding nonstationary equation:
Equation (3) can also be viewed as describing the diffusion in inhomogeneous medium, where is a spatially varying but stationary diffusion coefficient and is the concentration of the diffusing substance. This allows for a different physical interpretation of the model: we seek a steadystate spatial distribution of “particles” diffusing from 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].
2.4. 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 denote the length of a line segment originating at some point in WM boundary and ending in point . If, for any point , there is one and only one streamline passing through it, then 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 must satisfy the following PDE: with the boundary condition . 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 and a “source” point in the WM boundary . These correspondence functions can be found as solutions of the following PDE (see [24]): with boundary conditions , where .
The firstorder 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 onetoone correspondence property of the flow by clamping several streamline paths into one point on a grid, so the solutions and 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.
2.5. 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 are defined as a set of points, where spatial derivatives of the field are not wellbehaved, that is, the gradient 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 , 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 , as described by Han et al. [1]. The observation is that a centered finite difference numerical scheme will produce values of that are significantly lower than 1 on the shock points and are close to unity elsewhere. Then, the skeleton can be detected as , where is a minimum distance parameter set at the lower bound on cortical thickness and is a specified threshold value (; values mm and 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 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 , there are correspondences to source points that are “distant” between themselves. More formally, the skeleton can be determined as , where . We used voxels and 6 adjacent points in our computations.
2.6. Geometric Deformable Model
The geometric deformable model uses an implicit representation of a surface, embedding it into a level set function . The evolving interface is represented by the zerolevel set (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 is the advection velocity vector field, is the mean curvature, and are weights of the respective terms . The mean curvature of the interface embedded in the level set function is [16] where the subscripts denote partial derivatives. The advection velocity vector field is derived from the gradient of the potential or distance field : where is a stopping/directionreversal factor computed from the GM/WM class probabilities. For example, this factor can have a form of a logistic function: where is the constant controlling the steepness of the slope of the sigmoid curve and is the GM class probability threshold value that determines the “setpoint” of the deformable model. Figure 3 illustrates how the factor depends on GM and WM probability . In our experiments, a moderately steep sigmoid curve with and the threshold were used. For spatial regularization, the combined GM and WM class probability can be calculated as a weighted sum over the (closed) neighborhood of the point : where are the neighborhood weights (e.g., , where or 26, and for the central point ), and the skeleton of the sulcal medial surfaces 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 distanceconstraining factor: where the sign function is an “OR” combination of two signs: and the distanceconstraining factor can also have a form of a logistic function:
(a)
(b)
In (13), 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 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 narrowband algorithm [12, 16, 25]. The initial level set function is computed as a signeddistance 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 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 narrowband 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 builtin 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 (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 until convergence (i.e., until the relative amount of change in the SDF per iteration becomes small, for example, lower than ) or for a specified number of time steps and then smoothed by a few iterations with the curvature term, similarly to the inner surface.
3. Experiments and Results
Our algorithm was implemented in C++ in the Linux environment and ran on a PC with 2.5 GHz AMD64 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 T1weighted MR images of human brains, and on highresolution (submm) MR images of extracted brain hemispheres.
3.1. 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 (, 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.
(a)
(b)
(c)
(d)
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 crosssection 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. Sidebyside 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 cutoff 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 Gaussiandistributed noise and intensity inhomogeneity (see [13, 14]).
3.2. HighResolution MR Images
Our method’s performance is illustrated by results for highresolution 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 highresolution (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 AMD64 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 crosssections of a highresolution MR image are shown in Figure 6.
Lateral views of pial surfaces of three brain samples (3D rendering of thickness maps) are shown in Figure 7, left column. Measured thickness values (mean 2.2 mm; stdev 0.7 mm) are in good agreement with the literature. Inflated maps (Figure 7 middle and right column) are intended for better visualization of the surface inside sulci; they were produced with 20 iterations of Laplacian smoothing of the mesh. Maps in the right column are colorcoded with convexity values that were computed as vertex travel distances during smoothing/inflation, similarly to FreeSurfer [4]. On convexity maps, gyral crowns appear in blue color and sulcal fundi appear in yelloworange. 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.
3.3. CrossValidation with FreeSurfer: TestRetest Precision
Our method was validated by crosscomparison of cortical reconstruction results with those obtained using FreeSurfer. Standard resolution images for 20 righthanded healthy young subjects (age 19–34, average 23; 8 males/12 females) were selected from the crosssectional 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 T1weighted standard resolution images acquired per session. This relatively short delay between two consecutive scan sessions makes data sets suitable for the assessment of testretest 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 halfvoxel spacing. Implicit level set surfaces were tessellated using connectivityconsistent marching cubes algorithm [12], and triangular meshes were simplified down to 300,000 faces by a topologypreserving variant of the mesh simplification method [30]. DELFMAP processing took approximately 30 min per brain hemisphere (at halfvoxel 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.
(a)
(b)
The geometric precision or testretest 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 surfacewise mean and standard deviation were computed, as well as the groupwise statistics. In addition, we evaluated the testretest 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 pointcorrespondences 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 groupwise 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 and 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 values. For FreeSurfer WM surfaces, reproducibility is characterized by mean absolute error 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 mm (). For DELFMAP pial surfaces, reproducibility is characterized by a mean absolute error mm (L/R) that is similar in FreeSurfer (L: , R: , see details in Table 1). The standard deviation of the absolute distance is much lower in DELFMAP than in FreeSurfer (L: , R: ) which can be interpreted as a “tighter” reconstruction of pial surfaces in DELFMAP. Table 1 summarizes the statistics of the testretest analysis. The mean absolute difference of the cortical thickness is similar in both methods (L: , R: ), but the corresponding standard deviation is again much smaller in DELFMAP than in FreeSurfer (L: , R: ). To summarize, testretest 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.
3.4. CrossValidation with FreeSurfer: Intermethod Accuracy
The geometric accuracy of our method was evaluated by crosscomparison 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.2) was computed on pial surface. Next, maps of intermethod geometric differences (, , , ) were computed as signed distances (A.1) between white or pial surfaces reconstructed with FreeSurfer and DELFMAP (Wd, Gd). On these geometricdifference maps (40 sets, four maps per set), surfacewise statistics D_{mean}, D_{stdev}, , and (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 ). The 40 individual maps were resampled to a common template and averaged into groupwise maps of mean difference and standard deviation. The groupwise maps of intermethod cortical thickness measurement differences allow to assess and visualize any 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 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 ( 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: mm, L/R; positive sign here means that DELFMAPmeasured thickness is larger w.r.t. FreeSurfer). The intermethod accuracy, characterized by the mean absolute difference ( 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 sidebyside, 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, 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 groupaverage maps of geometric (Figure 11) and cortical thickness differences (Figure 12), where the above outlined three trends are also noticeable.
4. Discussion
We presented a novel PDEbased 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 highresolution 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 selfintersection (e.g., FreeSurfer) or mesh selfproximity (CLASP) term, which does not scale linearly with increasing mesh resolution. Although the computational cost of the straightforward mesh selfproximity term [2], which is quadratic on the number of faces , 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 selfintersection detection, which is based on intersection of bounding boxes, is [35]. In contrast to this, the computational complexity of the narrowband level set algorithm used in our method (and in CRUISE) is linear with respect to the size of the interface ( 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 secondorder PDE in DELFMAP versus a system of three secondorder 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 distanceconstraining 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 (a “setpoint”) and the maximal distance (which has strong influence only if set below the upper bound on cortical thickness).
The results from three highresolution data sets demonstrate that the method is capable of reconstructing the outer 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 crossvalidation 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. Groupaverage 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 over or 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.
Appendix
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 and was measured as: where is the closest orthogonal projection operator projecting a vertex from the first mesh onto one of the triangles 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 with the first surface outward normal at vertex ; 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 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 ; 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:
Acknowledgments
The authors gratefully acknowledge the Athinoula A. Martinos Center for Biomedical Imaging for providing FreeSurfer software and the Open Access Series of Imaging Studies (OASIS) project for providing MRI data sets, which were used for the validation of our new method. The OASIS project is made available by Dr. Randy Buckner at the Howard Hughes Medical Institute (HHMI) at Harvard University, the Neuroinformatics Research Group (NRG) at Washington University School of Medicine, and the Biomedical Informatics Research Network (BIRN). The authors would like to thank Professors T. Arendt and M. K. Brückner (PaulFlechsigInstitut für Hirnforschung, Leipzig) for providing specimens for exvivo imaging, and Professor C. J. Wiggins (MPI für Kognitions und Neurowissenschaften, Leipzig) for highresolution MRI acquisition, under project C15 supported by a grant from the Interdisziplinäres Zentrum für Klinische Forschung (IZKF), University of Leipzig.