|Home | About | Journals | Submit | Contact Us | Français|
Cortical thickness is an important biomarker for image-based studies of the brain. A Diffeomorphic registration based cortical thickness (DiReCT) measure is introduced where a continuous one-to-one correspondence between the gray matter-white matter interface and the estimated gray matter-cerebrospinal fluid interface is given by a Diffeomorphic mapping in the image space. Thickness is then defined in terms of a distance measure between the interfaces of this sheet like structure. This technique also provides a natural way to compute continuous estimates of thickness within buried sulci by preventing opposing gray matter banks from intersecting. In addition, the proposed method incorporates neuroanatomical constraints on thickness values as part of the mapping process. Evaluation of this method is presented on synthetic images. As an application to brain images, a longitudinal study of thickness change in frontotemporal dementia (FTD) spectrum disorder is reported.
A number of methodologies for cortical thickness measurements have appeared in neuroimaging literature over the past decade (Zeng et al., 1999; Miller et al., 2000; Jones et al., 2000; Fischl and Dale, 2000; Lohmann et al., 2003; Srivastava et al., 2003; Yezzi and Prince, 2003; Thompson et al., 2004; Lerch et al., 2005; Barta et al., 2005; Scott and Thacker, 2005; Hutton et al., 2008). Nonetheless, no gold standard has emerged that one can evaluate a measurement against. One reason is the difficulty of manually measuring thickness values in 3D images, unlike tissue segmentation where voxel-by-voxel manual marking can serve as a gold standard. Thickness measured directly using postmortem brains is also not considered an absolute metric because of possible tissue shrinkage. What complicates matters more is the lack of a consistent definition of cortical thickness. Some methods require explicit point associations between the white matter (WM)/gray matter (GM) surface and the GM/cerebrospinal fluid (CSF) surface. Some of these methods also require explicit construction of surface meshes (Fischl and Dale, 2000; MacDonald et al., 2000), typically using the extracted WM surface as a model that is fit to the GM surface by deformation, thus establishing node-to-node associations. Alternatively, definitions of thickness based on nearest point (Miller et al., 2000; Fischl and Dale, 2000) and distance along the surface normal (Das et al., 2007) do not require explicit point associations. However, regardless of whether a given method requires a priori correspondence maps or not, eventual computation of thickness is always based on some measure of distance between two points. This range of work points to a common definition of thickness that is based on two components: first, a principled point correspondence and, second, a distance measure between the points. Theoretical distinctions between algorithms occur in the definition of correspondence and/or the distance measurement between corresponding points. Practical differences usually occur in the underlying data representation.
Some researchers have argued that methods that require explicit surface extraction may suffer from inaccuracies due to surface generation (Srivastava et al., 2003; Thompson et al., 2005). One popular class of methods that establishes a priori point-to-point correspondence, without an explicit surface representation, is the gamut of PDE based methods (Jones et al., 2000; Yezzi and Prince, 2003; Rocha et al., 2007). Jones et al. (2000) models the cortical mantle as a dielectric and solves Laplace’s equation to compute electric field lines through the cortex in order to establish correspondence. There are still other image-based methods that rely on image operations such as morphological filters (Lohmann et al., 2003), geodesic distance transform (Srivastava et al., 2003), edge detection (Scott and Thacker, 2005) or explicit geometric modeling of the cortical sheet (Barta et al., 2005). The method introduced here is also image-based, but uses Diffeomorphic image registration based point-to-point correspondence.
One important feature of computational methods in neuroanatomy is the ability to include prior knowledge about the anatomy of interest. Such constraints increase reliability and performance in the presence of noise and when image resolution is sub-optimal. von Economo (1929) reported cortical thickness to be between 1.2–4.5 millimeters (mm) from ex vivo measurements, and ≈5 mm is generally reported to be the maximum observed value from in vivo measurements (Fischl and Dale, 2000; Kabani et al., 2001; Hutton et al., 2002). Some surface reconstruction based methods employ proximity constraints on the distance between two surfaces to ensure thickness lies within an anatomic range (Zeng et al., 1999; MacDonald et al., 2000). However, many methods, including Jones et al. (2000), provide no inherent anatomically motivated constraints on the computation of thickness. Note that the method based on solving Laplace’s equation within the cortical sheet provides thickness metrics based on curved paths. Such thickness values cannot be directly related to existing post-mortem knowledge, as provided by von Economo and others. Furthermore, straight-line/Euclidean distances are more easily interpreted in terms of a sheet-like model, as relevant for the cortex. Thus, the majority of thickness definitions use Euclidean distances that are consistent with what might be measured post-mortem. For these clinically relevant reasons – and for the purpose of using clinically motivated priors – we also focus on Euclidean distance as a thickness metric, though our technique applies as easily to curved thickness.
A related anatomical problem in thickness computation is the presence of buried cortex, which refers to the parts of the cortical surface hidden within the deep sulci. Within these sulci, thin strands of CSF are often mislabeled as GM due to partial volume effects, limited resolution and noisy tissue intensity levels. This mislabeling may lead to overestimation of thickness. Note that we use the terms buried cortex and buried sulci to refer to these regions of unresolved CSF, as has been used by other studies of cortical thickness (Kim et al., 2005; Hutton et al., 2008), as opposed to referring to the cortex in deep sulcal regions in general (with CSF separating the sulcal banks resolved or not) that is termed buried rather than cortex visible on the outer surface, as described in Van Essen et al. (1998); Rettmann et al. (2006). Several existing approaches explicitly address this problem: by image enhancement (Jones et al., 2000) to recover CSF, using morphological processing (Lohmann et al., 2003), stochastic modeling (Barta et al., 2005) and explicit image-based labeling of sulcal points (Hutton et al., 2002). The ACE method in Han et al. (2004) uses a level set based framework that corrects for buried sulci using fuzzy CSF membership information. Our own prior work introduced a novel topology preserving segmentation method (Das et al., 2007) that is able to recover some deep sulci. A similar method that does not preserve topology uses the Laplacian to measure thickness on segmentations with digitally recovered sulci (Hutton et al., 2008).
This paper introduces a new method, which we call Diffeomorphic Registration based Cortical Thickness (DiReCT). DiReCT reduces the problems of buried cortex and unconstrained thickness measurement by using a prior-constrained estimate of the distance between the gray/white interface and the gray/cerebrospinal fluid interface. This strategy, in essence, assigns a distance between opposing faces of a thin, sheet-like structure, here, the cortical mantle. Our approach begins with a hypothetical, infinitesimally thin, cortical mantle with ε thickness that lies along the white matter/gray matter interface. The initial model thus has “open” sulci. This thin cortical mantle is then allowed to expand under a one-to-one, differentiable and invertible map (a diffeomorphism) toward the edges of the data-derived gray matter probability. This mapping gives a correspondence field that allows an estimate of the gray/csf tissue interface, and thus thickness. An overview of this process is depicted in Figure 1.
The correspondence problem that we solve is made well-posed by the diffeomorphic constraint, which seeks a minimally deforming solution (Dupuis et al., 1998). The topology preservation provided by diffeomorphisms, along with explicit thickness priors, provide shape guidance to this mapping that prevent sulci from unfolding and neighboring banks of gray matter from intersecting. Thus, a maximum a posteriori estimate of the location of buried sulci is gained. Additionally, our thickness measures are guaranteed to stay within a range that is set by user-defined, spatially varying prior constraints. This shape-constrained mapping is able to represent buried sulci as infinitesimally close, opposing banks of gray matter, in the continuous domain, thus yielding sub-voxel resolution. Further, each gray matter bank is connected Diffeomorphically to a unique white matter bank.
In summary, the advantage of our approach, relative to prior work, is that we directly use standard probabilistic segmentation maps – with no further processing – to make a continuous, sub-voxel estimate of cortical thickness, constrained by thickness priors and accurate within buried sulci, within a single algorithmic framework.
The rest of the paper is organized as follows. Section 2.1 introduces various terminologies used. Section 2.2 describes the method for defining correspondences from Diffeomorphic maps and computing cortical thickness. Section 2.3 describes the construction of synthetic images used for evaluation. Methods used for conducting clinical study are discussed in Section 2.4. Section 3.1 presents results of evaluation of our method on synthetic data. Evaluation on a 2D brain slice as well as results of a longitudinal study of atrophy using these methods in a diseased population is presented in Section 3.2. An application of our thickness measurements is presented in the form of a longitudinal study of atrophy in a patient population diagnosed with FTD spectrum disorder. The small number of studies of gray matter atrophy over time in FTD (Avants et al., 2007; Brambati et al., 2007) — compared to a substantial number of similar studies in case of Alzheimer’s Disease (AD) — have either used volumetric measures of atrophy or voxel based measures like the Jacobian. Various issues pertaining to the proposed cortical thickness measure are discussed in Section 4. Section 5 concludes with advantages of the proposed method and future work.
Methods for generating a DiReCT map from an MRI volume begin with a set of probabilistic segmentation images, labeling the local mixture of tissue probabilities for each of three tissues, the white matter (WM), gray matter (GM) and the cerebrospinal fluid (CSF), in the brain. A variety of algorithms may be used to derive such images, but we rely on FAST as given in Zhang et al. (2001). Each voxel of the probability images, which are of the same size and spacing, is constrained to sum to 1, such that Pg(x) + Pw(x) + Pc(x) = 1 throughout all x in Ω where Ω is the domain of the image, where Pg(x), Pw(x) and Pc(x) are the GM, WM and CSF probability images respectively. We denote an image containing the sum of any two of these probability images as, in the case of WM and GM, Pwg = Pw + Pg (see Figure 2). Finally, the gray/white interface (GWI) is defined as voxels within Pg where Pg(x) >= 0.5 and at least one neighboring voxel within the 26-connected neighborhood of x has Pw >= 0.5. Note that no explicit surface mesh is extracted from the WM segmentation, thus avoiding complexities of surface generation. Instead, voxels representing the GWI are identified in the image domain.
Our goal is, for every point along the GWI, to find a corresponding GM point as given by Pwg that closely estimates the likely edge of the gray matter cortex. As argued in MacDonald et al. (2000), the cortical GM surface should have the same topology and a coarsely similar shape as the GWI. Note that this GM surface edge is not necessarily bounded, within the segmentation probability images, by neighboring CSF values with high probability, due to partial voluming and the problem of buried cortex. We denote the correspondence map as (x, t): Ω × [0, 1] → Ω and specify that (x, t) is a diffeomorphism, or a differentiable map with a differentiable inverse. Note that, while we are interested in mapping the GWI to a likely location of the cortical surface, we embed the mapping of the GWI within a volumetric diffeomorphism. We choose the Diffeomorphic transformation space because curves are guaranteed to remain curves and surfaces are guaranteed to remain surfaces under a diffeomorphic mapping. Thus, diffeomorphisms provide a natural shape constraint on the deformation of the GWI towards the cortical surface, which we augment with a thickness constraint. We detail, below, our approach for computing this mapping.
Define the velocity field, v(x, t) : Ω × t [0, 1] → d, where t is the “time” parameter, d is dimension, and v is a regularized vector field. We restrict the velocity field v to exist in the space of continuous and square-integrable functions with bounded derivatives by choosing a linear operator, L, which measures the energy of v, integrated over time. In practice, we choose L = 2 + (where is the identity operator) which has a Green’s kernel that is closely approximated by Gaussian smoothing.
The fundamental theorem of ODEs guarantees that this diffeomorphism is unique and varies continuously with initial conditions, (x, 0). Here, (x, 0) = x for all computations. Standard numerical integration techniques, such as Runge-Kutta methods, may be used to perform the integration. A diffeomorphic map may be applied to deform an image, I, through I((x)).
Given this transformation space, our goal is to generate a time one mapping, (x, 1), that estimates correspondence between the GWI and the outer edge of the cortical sheet, with given thickness constraints. Let M(x) be the moving image that represents the initial model of the gray matter sheet, derived from Pw(x), that is to be deformed under this mapping. Thus, we minimize the following variational energy, which captures both the notion that the velocity field should be smooth and the mapping should be a diffeomorphism, with constrained Euclidean displacement T(x, t),
The model image M(x) can be taken as the probabilistic WM segmentation image Pw(x), or a binary white matter segmentation image derived by thresholding Pw(x) at 0.5. In the latter case, one would have non-zero moving image gradient only at the GWI. Our experiments, using FAST (Zhang et al., 2001) as the segmentation method, have shown that the resultant mapping of the GWI, and hence the thickness measures, are nearly identical in either case. In all experiments reported here, we have used M(x) = Pw(x). M ((x, 0)) represents the GWI as described in section 2.1. Then, M((x, 1)) represents the other side of the cortical sheet, the estimated gray/CSF interface, in the same way. The function τ: Ω → + provides a prior on the particle x’s Euclidean distance, from its origin, as it transports along the diffeomorphic path. Note that, given a Diffeomorphic correspondence, thickness T(x, t) at a GWI point x is defined as the Euclidean length of the displacement that (x, t) imposes on x: T(x, t) = |(x, t) − x|. The inverse mapping −1(x, 1) is easily computed in the Diffeomorphic space by integrating the negative velocity field, as in Avants et al. (2006). Note that because −1(((x, 1), 1) = (identity map), we have an invertible correspondence map between the GM and GWI surface points.
Most critically, since we solve for a Diffeomorphic (differentiable, one-to-one and onto) map, the topology of the moving image (the WM segmentation in this case) cannot change. This is consistent with our treatment of the WM segmentation as the space in which we transport our GWI model toward the GM segmentation. This is one way our model preserves the shape of the resolved sulcus in the WM segmentation, thus recovering extremely narrow sulci that might have been mislabeled as GM. In addition to preserving topology, our optimization scheme seeks a minimally deforming solution (first term in equation 2) for a given level of image similarity between the deformed image and the target image (the data term in registration, second term in equation 2). Consider possible solutions to the registration problem for a closed sulcus such as that shown in Figure 1(a). A solution that maps the deepest points in the sulcus in the moving image (white/gray interface points) towards the “gyral” region (Figure 3(a), we’ll call it solution A) provides a comparable level of image similarity as one that maps these points towards the center of the sulcus (Figure 3(b), we’ll call it solution B). This is because, solution B maps the gray/white interface points from opposing banks of sulci to points that are infinitesimally close to each other in the continuous domain – thus the deformed image matches the target image as closely as it does for solution A. In other words, the deformed image “engulfs” the entire closed sulcus equally well in both cases in the continuous domain. In the discrete domain, however, we can see the single-pixel thick resolved sulcus for solution B. In fact, a range of solutions in between these two solutions exist, that resolve a progressively shallower sulcus; all of these solutions have similar image matching quality. Solution B, however, is the one that imposes the least amount of total deformation, and is hence the preferred solution, providing an estimate of a likely gray/CSF boundary. Finally, an explicit thickness prior in our optimization limits the maximum amount of deformation a white/gray interface point can undergo – which may also prevent solutions such as solution A which entails larger deformations. Thus, our shape-constrained mapping prefers to generate solutions in the class of panel (b) as opposed to panel (a) in terms of the classic buried sulcus recovery problem. These solutions are still far from perfect, particularly when CSF is not visible. At the same time, they represent a principled approach in terms of estimating thickness in the absence of perfect segmentation.
In the absence of a prior constraint on the thickness, one may optimize equation 2 by gradient descent on the Euler-Lagrange equations given by Beg et al. (2005) or, relatedly, a greedy approach as described in Avants et al. (2006). Optimization of the energy with the prior constraint can be performed with only a slight modification to use a constrained integration of the ODE that respects the thickness constraint τ. We detail this optimization, assuming a standard integration technique, below, and n steps in the temporal discretization of the velocity field. We take τ(x) and the probability images as inputs to the method and only T(x, 1) is the desired output.
Thus, we perform a greedy minimization of E within the space of diffeomorphisms while maintaining T(x) < τ(x) explicitly.
The map (x, t) defines the deformation path a GWI point follows through the GM mantle as a function of time t, reaching its corresponding GM surface point at t = 1. Let Tsurf (x) be the deformation norm image which has the value of measured thickness at every GWI point, and zero elsewhere. In order to propagate the thickness values Tsurf (x) volumetrically, we sample (x, t) at time intervals of t in [0, 1] and successively apply the transformations (x, t) such that, . This corresponds to transporting the thin thickness valued surface through the same velocity field that generates (x, t) and labeling each point in space with the passing value. In theory, this mapping would provide a dense field of thickness values throughout the domain traversed by Tsurf ((x, t)), as each point would be passed only once. However, due to discretization, care must be taken in how this integration is performed. We detail the discretized algorithm as,
This method iteratively accumulates the propagating level sets over time in the same image space. The thickness values on the WM surface thus propagate along the deformation path specified by (x, t) into the GM volume, generating the dense volumetric thickness map image Tvol, which we refer to as the DiReCT map.
We tested our registration based thickness measurements on a set of 3D phantom images at different image resolutions and noise levels. The gray matter mantle in these images is defined either as a spherical shell or one with cortex-like undulations. Similar phantoms have been used by Srivastava et al. (2003) in 2D. Mathematically, the phantom is given by a family of surfaces S : R2 R3 parameterized by θ and , the azimuth and elevation angles respectively – centered at c :
where R(γ, α) represents the radial distance from the center c of the volume enclosed by the surface S. R determines the shape of the surface. We consider different forms of the function f(γ, α, θ, ) for different shapes of phantoms that were studied.
In the simplest case, f(·) = 0, R = r and S reduces to a sphere of radius r centered at c. If SW (θ, ; R = rw) and SG(θ, ; R = rg) are nominally defined as the WM and GM surfaces respectively (rg > rw), the cortical gray matter would be a spherical shell of radius rg−rw, which also represents the theoretical cortical thickness value at every WM surface point on SW.
Next we consider a surface with ribbon-like undulations mimicking the cortex. This is realized as
where γ controls the magnitude and α controls the frequency of the undulations along the surface. Again, we construct the cortical mantle by defining the WM and GM surfaces as SW (θ, ; R = rw + fr(γ, α, )) and SG(θ, ; R =rg + fr(γ, α, ))(rg > rw) respectively. An example of such a phantom is shown in Figure 4(a).
The undulations in the phantom above are only along the (elevation) axis as can be seen from Figure 4(a). Any cross-section of the surface parallel to the z axis is a circle, since fr is not a function of θ, the azimuth angle. Now we introduce a third phantom shape with undulations along both θ and where
In this case, the phantom has bigger variations in curvature along the surface, resulting in deep sulcus-like structures in some areas. Cortical mantle is constructed the same way as described in section 2.3.2 with SW (θ, ; R =rw + fs(γ, α, θ, )) and SG(θ, ; R = rg + fs(γ, α, θ, ∀))(rg > rw). An example of such a phantom is shown in Figure 4(b).
For each phantom shape described above, we created phantoms with the same values of rw = 7 mm, rg = 10 mm, γ = 2, α = 5, but with different spatial resolution and noise levels. Binary WM and GM masks are generated by spanning the parameter spaces as = −π(Δ )π, θ = 0(Δ θ)π and letting r = 0(Δr)rw and r = rw(Δr)rg respectively and discretizing at the appropriate resolution. Voxels thus labeled as WM and GM are then set to intensity values obtained from mean WM and GM intensity values in sample brain images. The phantom image I obtained in this manner is then corrupted by noise as Inoisy = Gμ (I + νσ) where νσ is an additive Gaussian noise with zero mean and variance σ and Gμ is a Gaussian smoothing kernel with spread μ. Table 1 summarizes the parameters used. Three noise levels are nominally taken as none, low and high with increasing values of both μ and σ. To measure thickness, the WM and GM probability maps are generated by applying the FAST (Zhang et al., 2001) segmentation algorithm to the noisy phantoms.
We applied the methodology to a longitudinal study of gray matter thickness changes in a patient cohort diagnosed with frontotemporal dementia. Previous neuroimaging studies on FTD patient populations have identified locally focused atrophy (Chan et al., 2001) and these observations have been validated through correlation of annualized gray matter loss with cognitive decline (Avants et al., 2005) and in patients with autopsy-confirmed disease (Grossman et al., 2007; Murray et al., 2007).
The dataset contains twenty elderly individuals in total, diagnosed with FTD spectrum disorders, imaged at two time points each. The patients were right-handed, high school-educated (education = 16.0 ± 2.7 years) native English-speakers. The age at the first assessment was 61.8 ± 6.8 years with a mean MMSE (mini-mental state examination) score of 21.4±7.3. The second assessment was conducted after an average of 12.4±3.3 months, ranging from 7 to 21 months, with a mean MMSE score of 18.1 ±9.4. Patients were identified clinically in the Department of Neurology at the University of Pennsylvania School of Medicine. Initial clinical diagnosis was established by an experienced neurologist (author MG) using a modification of published criteria (Price et al., 2001). Subsequently, at least two trained reviewers of a consensus committee confirmed the presence of specific diagnostic criteria based on an independent review of the semi-structured history, mental status examination and neurological examination. Rare disagreements were resolved through discussion.
All images were acquired from a Siemens Trio 3.0 T MRI scanner. Each study began with a rapid sagittal T1-weighted scan to determine patient position. A T1 structural scan was then acquired using the MP-RAGE sequence with TR (repetition time) = 1620 ms, TE (echo time) = 3 s, slice thickness = 1 mm, in-plane resolution = 0.9766 × 0.9766 mm and field of view (FOV) = 256 × 256 × 192. The same protocol was used at baseline (younger) and followup (older) time points.
Groupwise morphometric studies are conventionally performed in a normalized coordinate system. We derive a template for this purpose from an expanded dataset (n=117) of FTD patients combined with healthy age-matched controls using symmetric Diffeomorphic normalization (SyN) methods (Avants and Gee, 2004). The template is one that is specific for the scanner, age and the MR protocol. This represents the average shape of the aged brain more accurately than a generic template and aids normalization accuracy. For each subject, we first find the mapping between the baseline (earlier time) image and the followup (later time) image. We then find the mapping between the baseline image and the template. We then compose these two maps to obtain the mapping between the followup image and template (Avants et al., 2005). In each case, the inverse maps were also generated. Once we derive thickness maps for each subject at both time points, we can warp them into the template coordinate system using these maps and perform statistical analysis in the template space. Generating maps between baseline image and template as well as followup image and template this way eliminates inconsistencies that may occur in the conventional approach which directly registers both the baseline and followup images separately to the template.
The semiautomated ITK-SNAP segmentation method (Yushkevich et al., 2006) is used to label the brain in the whole head template. Individual brain parenchyma are then extracted by mapping the template brain mask to each brain volume using the inverse maps obtained by SyN. Fuzzy segmentations for gray matter (GM), white matter (WM) and cerebrospinal fluid (CSF) are obtained for each patient scan using the FAST algorithm (Zhang et al., 2001). However, spatially varying three-tissue priors using tissue probability maps from the template brain are used to initialize the segmentation procedure.
We use the SyN map between the respective subject space at earlier (or later) time and the template space to map the individual DiReCT maps to the template space. The thickness value for the later time is then annualized by obtaining a linearly interpolated estimate of thickness at every voxel one year after baseline scan from the two values at earlier and later time. The assumption about linear progression of atrophy over time may not be valid in general as progression is likely to vary with disease stage. However, with only two time points available and in the absence of any other known temporal model, it is the only well-constrained assumption (Avants et al., 2005) that may be made. Once we have transported the DiReCT maps in the template space (Figure 9(d)–(e)), it is possible to perform groupwise longitudinal analysis of thickness changes. In order to study mean thickness change in Brodmann regions, we registered a template brain with Brodmann area labels to our local template.
Results of application of our thickness measurement algorithm in both phantom images and real brain images are presented below. For the phantom images derived from equation 4, we evaluate the correspondences produced by our method and discuss how the thickness measurements depend on image resolution and varying noise levels. We examine various properties of the method on a high resolution 2D slice of the brain. Finally, we describe results of a longitudinal study of thickness changes in a clinical population.
Since we measure thickness as the norm of the deformation field that maps a WM surface point to its corresponding GM surface point, a natural way to evaluate our method is to compare the correspondence map produced with some theoretically prescribed correspondence map. In the simplest case of the 3D phantoms used, the sphere, the corresponding point for any point on the WM surface would be the closest point on the GM surface, which is also the point on the GM surface along the direction of the WM surface normal, as well as the point on the GM surface in the radial direction. The thickness value is expected to be the same everywhere, and is given by T = rg − rw = 3 mm. Figure 5 shows the correspondence map produced by our algorithm.
Notice that the estimated thickness varies roughly between 2.75 – 3.25 mm because of partial volume effects at the low resolution. As the resolution increases, the effects of partial volume diminish and the estimated thickness approaches the theoretical value of 3 mm everywhere. This can be seen from the color of the arrows representing deformation field as well as from the gradual sharpening of the distribution of thickness values as resolution increases.
In case of the Ribbon image, as shown in Figure 6, there is no obvious choice for what kind of correspondence map one should use for measuring thickness. Correspondence based on surface normal direction may lead to inconsistencies at valleys with high curvature (Yezzi and Prince, 2003). Correspondence based on the closest GM surface point (Figure 6(a)) is another option. This can be easily defined on arbitrary surfaces once the surfaces have been extracted, but it doesn’t guarantee a smoothly varying deformation field for arbitrary shape, nor does it ensure one-to-one mapping. In case of geometric objects like the ribbon, one could use the radial direction to define correspondence (Figure 6(b)). However, this definition will not work for surfaces with arbitrary shapes. Registration based correspondence (Figure 6(c)), on the other hand, establishes one-to-one mapping and is applicable to arbitrarily shaped objects. It also generates smoothly varying deformation fields. In addition, as we shall see next, it can recover some deep sulci lost due to segmentation errors.
Figure 7 shows a slice through the Ribbon phantom at different noise levels and how it affects the correspondences found. As one would expect, with increasing noise, the histogram of the distribution of thickness widens. Yet, the correspondence map remains fairly robust across different noise levels. Measured thickness, as given by the color of correspondence arrows, varies in few locations. Note that the differences across noise levels may not necessarily reflect the level of reliability of the thickness measurement per se, but may mostly arise from the variability of the segmentation.
Importantly, even though the 3D phantoms described in this section are constructed using an analytical equation (Equation 4), only for certain simple choices of the thickness definition one could obtain an analytical ground truth for the thickness value, such as a definition based on nearest point. In fact, ground truth thickness values for a phantom of a given shape could vary substantially depending on how thickness is defined, as demonstrated in Figure 6. For Diffeomorphic correspondence based definition of thickness, no analytical ground truth is available for a direct evaluation. Nonetheless, it is informative to examine how a given correspondence model (and hence the spatial distribution of measured thickness) captures an intuitive notion of thickness for a thin ribbon-like object (Figure 6). It is also useful to evaluate how a method performs under limited image resolution (Figure 5), as is often the case in practice, and in the presence of noise (Figure 7). Hence, these are the aspects of evaluation we considered using synthetic images, rather than direct comparison with mathematical ground truth. Recovery of buried sulci is another important aspect of evaluation, and is considered in detail in the following subsection.
The behavior of the proposed method under different noise conditions and spatial resolution has been described with phantom images in section 3.1. Next, different properties of the method are illustrated with the help of a 2D high resolution image of a brain slice in Figure 8. Segmentations in (a), (c), (e) and (i) are drawn such that progressively thicker cortex and narrower sulci are visible. This way, we can examine the solutions obtained by DiReCT under slightly different GM segmentations. We begin with the segmentation in panel (a), where all the sulci are fully resolved (we also separate the gyral crowns by CSF even though they touch in the original image). Panel (b) shows the DiReCT correspondence map of (a). The smoothly varying deformation field, and consequently the smoothly varying thickness values are able to capture the variation of thickness visible in the segmentation of the cortical sheet.
In Figure 8(c), GM in a sulcus gets thicker and no CSF is visible except a few voxels in the fundus. The DiReCT map in (d) is able to resolve this sulcus accurately, and the yellowish color in the region marked A reflects the thicker gray matter there. In (e), the sulcus is completely closed. The correspondence map in (f) is still able to resolve a narrow sulcus. However, since there is no information available about the depth of the sulcus, as was provided by the few CSF voxels in (c), the method produces a thicker fundus, marked B.
Panels (e)–(h) in Figure 8 demonstrate the effect of the thickness prior τ. Sulcal CSF can remain unresolved in the segmentation, as in (e), for various reasons. For a sulcus narrower than a voxel, partial volume effect is often responsible for the CSF therein being mislabeled as GM. In this scenario, setting the prior to a high value compared to the thickness of the sulcal GM such that it essentially doesn’t constrain the measurements there would lead the diffeomorphic mapping to (correctly) resolve an infinitesimally thin sulcus. This is the case in (f), where the prior is set to a high value. Both sulcal banks have similar thickness in this case, which is the most unbiased solution. However, other factors, including segmentation errors, may label even a wider sulcus completely as GM. In this case, an unconstrained measurement would, erroneously, give rise to banks of GM on both sides that are too thick. Also, the two sulcal banks may have asymmetric thickness. If such prior information is known, the proposed method provides the flexibility to use it to take us one step closer to what is neuroanatomically plausible. This is demonstrated in (g) and (h). In (g), a slightly lower thickness prior than the thickest gyral regions (marked C) is used as a prior, leading to a solution with a thinner gyrus, with shape guidance provided by Diffeomorphic mapping. In (h), on the other hand, a spatially varying prior is used, with a prior of about one-fourth of the total thickness of two sulcal banks used for the sulcal bank marked D. This leads to a smooth solution with asymmetric thickness on both sides of the sulcus. In summary, thickness prior may be rendered ineffective by setting it to a very high value, yielding unconstrained thickness measurements, or it can be used in conjunction with Diffeomorphic mapping to get better estimates – if the user so desires and if reliable prior information is available. Using an arbitrary threshold as prior in thickness measurement on populations may bias group analysis. Therefore, we advocate using known, possibly spatially varying, neuroanatomical prior information only to the extent that this information is available. Indeed, when segmentation is perfect, the prior is unnecessary, and can be ignored.
Panels (i)–(l) in Figure 8 demonstrate the effect of the width σ of the smoothing kernel K in step 4 of the optimization method in section 2.2.1. In (i), in addition to the closed sulcus in (e), two gyral crowns touch each other. This leads to a solution in (j) that is different than (f), most notably in gyral region marked E, yet the two gyri are separated correctly. As the size of the smoothing kernel is doubled in (k) and quadrupled in (l), we get a smoother and smoother solution, but the GM/CSF boundaries do not match well, even though thickness prior is set to the same high value for (j)–(l). σ is set to 1 voxel width for (b)–(j). As expected, a larger kernel focuses more on larger scale features and is less sensitive to finer scale variations in cortical shape. In general, setting it of the order of one voxel width ensures a sufficiently smooth solution as well as good image similarity.
Next, we describe results of our clinical study.
We construct DiReCT maps for a longitudinal dataset of elderly individuals diagnosed with FTD spectrum disorders at two different time points. τ is set to 5 mm for results reported here. Figure 9 shows examples of DiReCT maps. Note the recovery of deep sulci as demonstrated in (f) and (g). In Figure 9(f) top panel, the original segmentation is able to recover only some CSF in the sulcus. The rest of the sulcus is resolved by our algorithm as shown in Figure 9(f) bottom panel. In Figure 9(g), the sulcus is completely closed. In this case, an infinitesimally narrow sulcus is recovered. These examples show the ability of our shape-constrained Diffeomorphic mapping to recover buried sulci in 3D brain images.
Mean thickness change in Brodmann areas is calculated for every subject. Pairwise t-tests for longitudinal thickness reduction are performed across all subjects in the template space. Figure 10 shows areas of statistically significant change at a false discovery rate (FDR) of 5% to correct for multiple comparisons overlaid on the rendered SyN template. Table 3 lists most significant Brodmann areas. Mean thickness changes are also shown in Figure 10.
Areas of significant atrophy include large parts of the temporal lobe, as well as some areas of parietal and frontal lobe. Jacobian-based voxelwise measures of atrophy on the same dataset revealed many of the same atrophied regions (Avants et al., 2007). The agreement of these independently gained results help support the robustness of these disease effects. These same areas have been shown to be atrophied compared to healthy age-matched controls in cross-sectional studies (Grossman et al., 2004). In addition, in Table 2, we compare the mean cortical thickness in different lobes of the brain in the baseline scans of our clinical population with published values (Du et al., 2007) in an FTD population of similar age group (61.7 ± 7.5 years as opposed to 61.8 ± 6.8 years for our dataset). The most significant atrophied regions (Table 3) show annual mean thickness reduction of more than 0.10 mm. In contrast, average reduction in the occipital lobe is 0.033 mm (not significant). Cortical thinning due to normal aging is usually reported to be far less (Salat et al., 2004; Preul et al., 2006) – less than 0.10 mm per decade – than the rate of decrease in thickness we found, suggesting atrophy related to neurodegeneration. Even though such corroborative evidence does not conclusively prove that the thickness changes reported here are entirely due to progression of disease, it is also very unlikely that such high rates of cortical thinning in areas known to be affected by the disease could be explained by other factors like normal aging. A voxelwise map of thickness reduction is shown in Figure 11.
These areas showing significant change over time are directly relevant to the patients participating in this study. Correlations of imaging studies with independent measures of clinical functioning and pathology validate the importance of these areas in FTD spectrum disorders. The atrophied areas in the left frontal and temporal regions are implicated in primary progressive aphasia, one presentation of these conditions; and the atrophied areas in the right hemisphere are limited to the disorder of social comportment and personality that is another presentation of these diseases. Parietal lobe atrophy may be associated with the worsening apraxia and short-term memory disorder of corticobasal degeneration (Murray et al., 2007).
In this section, we discuss some pertinent issues in the context of cortical thickness measurements as they relate to the proposed methodology.
We alluded to the lack of a gold standard for cortical thickness measurements in section 1. As we have demonstrated in section 3.1, common thickness measurements may be thought of as a distance measure between a pair of points; depending on how the pair of points is defined, measured thickness can vary greatly. Thus, it is not obvious how to compare clinical studies where thickness measurements are based upon different definitions of thickness. On the other hand, it is challenging to prove the validity of any single thickness definition. Labeling three dimensional brain images – a point operation – is an extremely challenging task for a human. Labeling points on surfaces that are likely to occur in different slices and different orientations – a line operation – is even more challenging, thus making it difficult to generate reliable ground truth thickness measures. This is one reason why we, and other researchers, have not made a quantitative comparison of our evaluation results to a ground truth measurement. Rather, we focused on discussing the variability of our measurements as a function of relevant imaging parameters such as image resolution and noise levels, as well as qualitative visual assessment of the derived correspondence maps. Also, we have demonstrated that our measurements are sensitive enough to be used in a clinical study.
Furthermore, any thickness measure is only as good as the quality of tissue segmentation it is based upon. This makes it difficult to assess the reliability of thickness measures in the presence of noise. While some studies have examined the reliability as a function of image acquisition parameters (Haidar and Soul, 2006; Sato et al., 2003; Han et al., 2006), over repeated scans (Sowell et al., 2004) and as a function of the particular distance measure used to define thickness (Lerch and Evans, 2005; Haidar and Soul, 2006) – no systematic study of the variability of thickness measurements as a function of the tissue segmentation algorithm used has been reported. The issue of recovering deep sulci is also a pertinent one in this context – since it is often segmentation errors coupled with partial volume effects that lead to buried sulci. We have taken a model-based approach to minimize the effects of segmentation errors. As long as the WM segmentation is locally neuroanatomically correct, our constrained Diffeomorphic mapping ensures that the derived GM correspondence preserves the initial anatomical model, in effect often correcting for erroneous GM over-segmentation. It can be noted here that methods that define thickness as distance between corresponding nodes on WM and GM surface meshes (MacDonald et al., 2000), although still dependent on the quality of segmentation, can be evaluated directly in the presence of noise by taking an initial surface reconstruction as the ground truth, corrupting it, and recovering it again (Kabani et al., 2001). The reliability of thickness measurements as a function of segmentation method used deserves to be studied further.
The problem of detecting buried sulci to improve thickness measurements has been addressed explicitly by many researchers (Hutton et al., 2002; Jones et al., 2000; Scott and Thacker, 2005; Lohmann et al., 2003; Barta et al., 2005; Han et al., 2004), often in the form of a pre or post-processing step. On the other hand, labeling cortical sulci directly from MRI volumes (Tu et al., 2007), or from surface models (Sandor and Leahy, 1997) is a well-studied problem in itself. It is hard to say whether any particular approach will produce the result that is closest to “reality” in the context of a given thickness measurement method, including ours. Nonetheless, the proposed solution in this paper is a natural consequence of our constrained Diffeomorphic mapping, and thus directly built into the definition and optimization of our point correspondence. We deal with thickness constraints directly in our optimization, as opposed to having to deal with the problem of thickness over-estimation separately. Once we find a topology preserving, thickness constrained fit of the GWI to the GM surface, the buried sulci are automatically recovered.
Many areas in the cortex are known to have asymmetric thickness on opposing sulcal banks. If spatially varying thickness prior information about such asymmetry is available, the proposed method will recover asymmetric thickness even if no CSF is visible. Also, note that whenever even some amount of CSF (reflected in lowered gray matter probability) is visible in the segmentation, this will guide the registration process (For example, in Figure 9 (f)). However, if no thickness prior information is available and no CSF is visible, our method will still be able to recover an infinitesimally narrow sulcus (For example, in Figure 9 (g)), but produce a symmetric thickness map on either sulcal bank, which is the most unbiased solution under those circumstances. We note that the estimated gray/CSF interface between opposing sulcal banks will be influenced by the specific cortical geometry, in particular the curvature of the fundus. In a “tighter” sulcus (which would usually imply greater curvature at the deepest part of the sulcus), the opposing banks will meet more quickly. Thus – in the absence of a prior – the estimated thickness at the depth of the sulcus will be lower than if the sulcus was wider. This implies the following relationship for sulci with zero resolved CSF: the narrower the sulcus, the lower the estimated thickness at the fundus will be. We believe that this is a natural relationship. However, it is only a “best estimate” as determined by our methodology.
We do not claim that ours is the best approach to recover deep sulci per se, but it is a principled one consistent with how we define point associations for measuring thickness. At the same time, our method naturally builds in neuroanatomic constraints on thickness, with the intention of providing a more plausible and accurate estimate that agrees with what is known about the in vivo brain and the folded, ribbon-like structure of gray matter hidden in buried sulci.
For clinical applications, a statistical framework for carrying out groupwise morphometric studies on measured thickness values is required. When measurements are made on explicitly constructed cortical surface meshes (Lerch et al., 2005; Chung et al., 2005; Fjell et al., 2006), or a hybrid approach where measurements are made in image volumes but then transferred to surface meshes by high dimensional mapping (Thompson et al., 2005) – surface based statistics can be used. Surface based smoothing is often used to increase the signal-to-noise ratio. Similar smoothing schemes can be used in the image domain as well (Avants and Gee, 2003), as was done in Das et al. (2007). In the longitudinal population study presented in section 3, because of the inherent smoothness of the deformation fields, and hence in the DiReCT maps, explicit smoothing was not necessary. Although image registration-based normalization as we have used here for group study poses an inherent limitation to accuracy, since we look at mean region-based quantities – statistics based in the template volume should be adequate. Moreover, for a longitudinal comparison, intra-subject registration error is minimal – even voxelwise statistics would be much more accurate than if it were a cross-sectional study. We note that even though volumetric thickness maps were used in this study, it is obviously possible to use surface based statistics in the SyN template – either on the WM surface in the image domain or by extracting the WM surface mesh on each individual. In that case, one would not need to propagate the thickness measurements to the GM volume. In the clinical study presented here, we use propagation along the dense and Diffeomorphic deformation path to ensure that all GM voxels are assigned a thickness. This is also useful in visualization.
In this paper, we introduced a principled approach to finding correspondence between WM and GM surfaces using shape-constrained Diffeomorphic mapping. We evaluated this methodology with 3D phantom data at different spatial resolutions and noise levels. We constructed examples to show the merit of buried sulcus recovery, the use of thickness priors and the effect of the size of smoothing kernel in a 2D brain slice data. Finally, as an application of the methodology to real brain images, we present a groupwise longitudinal study of thickness changes in a patient population.
The advantages of DiReCT in comparison to other approaches include:
In summary, the proposed method integrates many advantages of various existing approaches. DiReCT incorporates the model-based, deformation-controlled and regularized approach of methods that use surface deformation, yet works in the volumetric image space. DiReCT also provides an anatomically motivated and prior-constrained, subvoxel measurement of thickness along with an estimate of buried sulci, within the same algorithmic framework. Finally, we showed that DiReCT produces clinically meaningful measurements in a study of cortical thickness within Brodmann areas defined in an FTD cohort.
Our future work will continue to investigate the effect of our new thickness constraints on clinical studies, will incorporate the estimate of GM probability in the thickness estimation step and will contrast surface-based population studies with the volumetric approach used here. We intend to extend the proposed method to use spherical topology representations which are required for inflation and flattening, as in Das et al. (2007). We will also investigate the feasibility of labeling protocols for generating gold standard thickness measurements.
This work was partially supported by NIH grants EB006266, NS045839, DA022807, DA14129 and HD046159.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.