Diffusion tensor magnetic resonance imaging (DT-MRI, shortened as DTI) produces, from a set of diffusion-weighted magnetic resonance images, tensor-valued images where each voxel is assigned a 3×3 symmetric, positive-definite matrix. This tensor is simply the covariance matrix of a local Gaussian process with zero mean, modelling the average motion of water molecules. We propose a three-dimensional geometric flow-based model to segment the main core of cerebral white matter fibre tracts from DTI. The segmentation is carried out with a front propagation algorithm. The front is a three-dimensional surface that evolves along its normal direction with speed that is proportional to a linear combination of two measures: a similarity measure and a consistency measure. The similarity measure computes the similarity of the diffusion tensors at a voxel and its neighbouring voxels along the normal to the front; the consistency measure is able to speed up the propagation at locations where the evolving front is more consistent with the diffusion tensor field, to remove noise effect to some extent, and thus to improve results. We validate the proposed model and compare it with some other methods using synthetic and human brain DTI data; experimental results indicate that the proposed model improves the accuracy and efficiency in segmentation.
Diffusion tensor magnetic resonance imaging (DTI; Basser et al. 1994, 2000) is a modality that permits non-invasive quantification of the water diffusion in living tissues. It adds to conventional magnetic resonance imaging (MRI) the capability of measuring the random motion of water molecules, referred to as diffusion. DTI is extremely useful in identifying the neural connectivity patterns of the human brain and allows one to distinguish the anatomical structures of the cerebral white matter such as the corpus callosum, the superior longitudinal fasciculus or the corticospinal tract. Recent development of the Q-Ball imaging (e.g. Tuch 2004; Descoteaux & Deriche 2007) technique enables us to reveal more complicated tissue structures by detecting multi-fibres passing a voxel. The segmentation method proposed in this work is for DTI data only, but the idea can be generalized to higher angular resolution diffusion-weighted MRI data.
There have been many models developed for DTI segmentation in the literature, which might be outlined as three classes. One of them is based on a clustering technique. Wiegell et al. (2003) applied a K-mean algorithm to segment thalamic nuclei from DTI data. Ziyan et al. (2006) used spectral clustering to segment corticothalamic/thalamocortical striations within each nucleus.
The second class is region-based active contour methods. The approaches in this class can be viewed as the extension of region-based active contours for scalar images to matrix-valued images or to the probability density functions (p.d.f.s) of water diffusion, which are Gaussian distributions with zero mean and the corresponding diffusion tensors as their covariance matrices. Work in this streamline includes Lengthlet et al. (2004, 2006), Rousson et al. (2004), Wang & Vemuri (2004) and Awate et al. (2007).
The third class of models is based on a front propagation approach, which uses tensor information near the front (edges) rather than the regional tensor statistics for segmentation. The models in this class can be considered as an extension of the geometric active contour or diffusion equations for scalar image segmentation to DTI segmentation. For instance, Zhukov et al. (2003) proposed a level set method to segment the images of the eigenvalues of the diffusion tensor. A drawback with Zhukov's method is that the directional information of the water diffusion is ignored. This can lead to low discrimination capability, potentially yielding the segmentation of mixed structures. Feddern et al. (2003) extended geodesic active contours to DTI segmentation by using gradient magnitude of the matrix-valued image as the stopping criteria. Pichon et al. (2003) proposed a diffusion equation that replaces the stopping criteria in geodesic active contour by an alignment penalty of the curve tangent to the dominant eigenvector. Jonasson et al. (2005) presented a geometric flow, where the propagation force depends on the similarity of the tensor near the front.
The proposed model belongs to the third class. We propose a geometric flow-based model, where the propagation force is determined by both similarity of the diffusion tensors near the front and the consistency or coherence between the evolving front and the tensor field. More precisely, we extend the work of Jonasson et al. (2005) by introducing a propagation force that is a linear combination of two measures: a similarity measure of any type and a consistency measure. To our knowledge, this is the first work to incorporate the relationship between the evolving front and the tensor field into a geometric flow, although this kind of idea has been applied to fibre tracking. We aim at using the consistency measure to obtain DTI segmentation results more accurately and efficiently.
The organization of this paper is the following. In §2, we first briefly present the concept of DTI and review certain similarity measures of interest, then provide basic theories on geometric flows and their level set implementation. Section 3 focuses on the description of the proposed model. Section 4 illustrates experimental results for both synthetic and real human brain DTI data. Comparison between models with and without the consistency measure is presented to show the advantage of the proposed model. Discussion is presented in §5.
(a) Diffusion tensor imaging
At a given location, the diffusion of water molecules in tissues over a time interval with length t can be described by a p.d.f. (pt(r)) on the displacement r. The standard methodology employed in most DTI experiments is the Stejskal–Tanner pulsed-gradient spin-echo method (Stejskal & Tanner 1965). Two magnetic field gradient pulses of strength G and duration δ with a temporal separation of t between the onset of the pulses are applied to the simple spin-echo sequence. If the duration of the pulses δ is negligible compared with t, the attenuation of the MR signal s(q) with respect to the diffusion sensitizing gradient q measures the Fourier transformation of the average p.d.f. (pt(r)) on a spin displacement r over diffusion time t (Basser et al. 1994). That is(2.1)where q=(2π)−1γδG, γ is the gyromagnetic ratio of protons in water, and s0 is MRI signal in the absence of any gradient. In 1994, Basser et al. proposed to model the local p.d.f. of the three-dimensional molecular motion by a Gaussian distribution with zero mean and covariance D, i.e.(2.2)By substituting this density function into (2.1), we obtain the following Stejskal–Tanner equation:(2.3)where u=q/|q| and b=4π2|q|2(t−δ/3) is the diffusion weighting factor.
To estimate the diffusion tensor D from the measurements s(q) and s0, we require the acquisition of diffusion weighted images in at least six different sampling directions (6 q′s). There have been many algorithms on the estimation and regularization of this tensor field: Westin et al. (2002), Mangin et al. (2003), Tschumperlé & Deriche (2003), Wang et al. (2003), Chefd'hotel et al. (2004) and Coulon et al. (2004).
At each voxel, the diffusion tensor D is a 3×3 symmetric, positive definite matrix; it can be characterized by its three positive eigenvalues λ1, λ2 and λ3 (assuming λ1≥λ2≥λ3) and corresponding eigenvectors e1,e2 and e3:(2.4)
Water diffusion is highly anisotropic and oriented in areas of compact nerve fibre organization, but isotropic anywhere else. The eigenvector e1, corresponding to the largest eigenvalue λ1, is usually called the principal eigenvector (PE) and provides the preferred direction of diffusion. The eigenvalues can be used to compute a scalar measure of the degree of anisotropy. A commonly used anisotropy index is fractional anisotropy (FA) that ranges from 0 to 1(2.5)where is the mean diffusion (λ1+λ2+λ3)/3. The FA index is quantitative and dimensionless. For an isotropic medium, all the three eigenvalues are the same and FA will be 0. When , water diffuses mainly along e1 and FA will be close to 1.
(b) Certain similarity/dissimilarity measures
In this section, we will recall a few similarity or dissimilarity measures between the two tensors D1 and D2 used in DTI segmentation.
One of the most intuitive similarity measures would be the following Euclidean distance:(2.6)
Several other similarity measures for tensors have been designed in Alexander et al. (1999) to perform elastic matching of diffusion tensor images. These measures make use of both magnitudes and the direction information of the diffusion tensor. One example is the tensor scalar product (TSP), which is defined as(2.7)It measures the overlap between two tensors. TSP is large when the two tensors have both high eigenvalues and similar eigenvectors. To remove the effect of relative size (in the sense of the magnitude of order of entries) of the two tensors, a normalized TSP, shortened as NTSP, is used more often to emphasize the shape and orientation of the tensor:(2.8)
Another type of dissimilarity measures of two tensors is based on the distance of their corresponding p.d.f.s (pt(r)). Note that the diffusion tensor is the covariance matrix of a local Gaussian distribution. A natural measure of dissimilarity between the two tensors D1 and D2 would be the distance between two Gaussian distributions with zero mean and the covariance matrices D1 and D2. Therefore, by using Kullback–Leibler (KL) divergence defined for the two general p.d.f.s f and g as(2.9)the KL divergence of two Gaussian p.d.f.s mentioned above can be computed as(2.10)KL(f, g) is not symmetric. The J divergence, defined as mean of KL(f, g) and KL(g, f), is a symmetrized KL divergence and is sometimes used.
For more general p.d.f.s, f and g, the Bhattacharyya distance is defined as(2.11)and is used sometimes because it is easier to calculate than KL or J divergence.
The Bhattacharyya distance for the two Gaussian p.d.f.s mentioned above is(2.12)
Kailath (1967) compared the properties of J divergence and the Bhattacharyya distance. It was found that the latter often works as well as the former.
Lengthlet et al. (2006) considered the Riemannian manifold of three-dimensional normal distributions with zero mean parametrized by the six components of their covariance matrix D (the diffusion tensor), where a Riemannian metric was introduced in terms of the Fisher information matrix. By this view, the geodesic distance Dg on was defined as follows:(2.13)where ηi denotes the ith eigenvalue of matrix . The geodesic distance was compared with the Euclidean distance and J divergence in their work; it was claimed that the geodesic distance outperforms the other two in segmentation tasks.
(c) Geometric flows and level set method
A geometric flow-based method is to evolve a curve or surface (called a front) with a velocity depending on the external properties determined by the image features or other factors. One example is curvature shortening flow that evolves a curve or surface at each point along its normal with a velocity depending on the mean curvature at that point. This process results in the smoothing of the front to eliminate the effects of noise. It becomes an important regularization tool in computer vision.
A general form of a geometric flow for a closed curve/surface can be described as(2.14)where F is an image-based speed function, κ is the curvature/mean curvature of the curve/surface S and α is a balancing factor. N is the normal to the curve/surface and t is time.
To solve the time-dependent partial differential equation (2.14), we use the celebrated level set method introduced by Osher & Sethian (1988), where the evolving curve/surface is considered as the zero level set of a function ϕ that is of one dimension higher. By doing this, a numerically stable algorithm that easily handles topology changes of the evolving curve/surface is obtained. The level set formulation of (2.14) reads as(2.15)
For DTI segmentation tasks, most of the work in the literature has defined F as a dissimilarity measure of tensors near the front. For example, Jonasson et al. (2005) defined F as a mean of two NTSP:(2.16)where Di is the diffusion tensor at the current voxel i and Di−p is the diffusion tensor at voxel i−p that is p voxels backwards along the normal of the front (figure 1a). NTSP(Di,Di−1) is the normalized TSP defined in (2.8). The fundamental assumption of the segmentation technique in Jonasson et al. (2005) is that adjacent voxels in a tract have similar diffusion property. Our proposed model is still under this assumption; however, we will also take the consistency between the propagating front and the tensor field into account to obtain a better segmentation result in a shorter time.
3. Proposed model
We first explain the motivation of this work through a two-dimensional tensor field as shown in figure 1b. At each location, a black arrow represents the PE of the diffusion tensor. To segment the arch-shaped tensor field from the homogeneous background, an initial curve (shown as a black ellipse in figure 1b) is set inside the object and then evolves along its normal direction. PEs are almost horizontal for points on the ellipse. The vertical and horizontal arrows show normal directions at corresponding locations on the ellipse. Intuitively, propagation speed along horizontal arrows should be larger than that along vertical arrows; the reason is that the evolving contour is farther away from the boundary at locations with horizontal arrows, while closer to the boundary at locations with vertical arrows. Further observation gives that normal directions at points with horizontal arrows are more consistent with (parallel to) the PE field than those at locations with vertical arrows. The normal directions at the locations with vertical arrows are almost perpendicular to the nearby PE field. Thus, propagation should be faster at locations where the normal directions are more consistent with the orientations of the PEs of the diffusion tensors. We introduce a consistency measure CONS to quantify how consistent the evolving front is with the orientation of the PE of the diffusion tensor at a voxel.
An intuitive consistency at location i is defined to be(3.1)with PE the principal eigenvector of the current diffusion tensor Di and N the unit normal direction of the front at i. Absolute value is taken to remove the sign ambiguity of the PE. If the normal direction of the front is parallel to the principal direction of the tensor, the evolution direction is consistent with that of the PE, the evolution speed is higher along that direction. A drawback of CONS1 is that it totally relies on the PE, which is sensitive to noise.
The above Di*N is a vector that will turn towards the PE of the tensor Di. When the tensor has lower than full rank, it will act as a projection operator. In particular, when the tensor is spherical or when N is one of the eigenvectors of Di, N will not turn, i.e. , which implies obtains maximum value 1 in all these situations. To distinguish between a highly anisotropic tensor and an isotropic tensor that both lead to the high value of , FA is multiplied in the front of it to form CONS2. Lazar et al. (2003) claimed that CONS2 is less sensitive to noise. We will apply CONS1 to synthetic data that are not very noisy and CONS2 to the noisy human brain data.
(b) A new evolution force
We propose a new force F that is a linear combination of a similarity measure (SIM) and a consistency measure (CONS) with a balancing factor β, i.e.(3.3)SIM can be any of the measures defined in §2b, CONS could be chosen as CONS1 or CONS2. The systematic strategy for choosing β is that we select the one such that the mean value of the βCONS over the mean value of the SIM is approximately 0.5. This is to make sure the βCONS term does not overweigh the SIM term because the similarity term is more important in determining the main structure of the to-be-segmented object. Based on this strategy, β is set to be 0.8 and 0.5 in the synthetic and human brain data, respectively.
One specific example of propagation force at location i, for instance, is(3.4)In comparison with the work of Jonasson et al. (2005), where the evolution force is defined as(3.5)we not only add the consistence measure in the force, in the situation when the tensor field is noisy, the calculated F is noisy also, the segmentation based on which will be noisy. A simple filter, for example, a convolution with a Gaussian kernel, can be applied to F to obtain smoother results. The added term CONS, which does not cost much time to calculate, speeds up the propagation of the evolving surface. Therefore, the proposed model takes much less iterations and thus much less time.
(c) Termination and initials
Termination criteria are carefully chosen based on two threshold values for F and similarity SIM, which are called ThF and ThS, respectively: stop propagation of the front at locations where F<ThF and SIM<ThS in two successive steps. The whole algorithm stops when each point on the front terminates. In comparison, Jonasson et al. (2005) only used one threshold ThS for similarity to terminate. The advantage of the former over the latter will be presented in §4.
Initial front is obtained by a quick and thin delineation inside the to-be-segmented object. This makes segmentation more efficient than taking a circle/sphere as the initial front.
4. Experimental results
This section aims to compare the segmentation quality and convergence speed to show the advantages of the proposed model that uses the consistency measure over some models that do not use the consistency term.
(a) Synthetic results
The synthetic data shown in figure 2a are created by replacing the middle part of a randomly generated two-dimensional tensor (2×2 positive definite matrix) field of size 55×55 with an arch-shaped tensor field. At each point (x, y), the diffusion tensor equals V′..V. If (x, y) is inside the arch-shaped object, D is a diagonal matrix with diagonal entries 0.2×10−3 and 0.8×10−3 and V is defined as in (4.1) with ; if (x, y) is on the background, D is a diagonal matrix with diagonal entries 0.7×10−3 and 0.4×10−3, V is defined as in (4.1) with θ a random angle between 0 and π/2:(4.1)
To demonstrate the advantage of the consistency term in the proposed model, we create a tensor field based on the above field. First, inside the dashed square in figure 2a, as zoomed in and shown in figure 2b, there are two perturbed tensors such that similarities between them and their neighbouring tensors drop below the threshold ThS that is appropriate for Jonasson et al.'s (2005) termination. The front will stop at these two locations by Jonasson's model and end up with insufficient evolution of the front or a hole. But the normal of the front is consistent with the tensor, the consistency term CONS1 is still high, so it will bring up the sum F; therefore when an appropriate threshold ThF is chosen such that F>ThF, then even though similarity is less than ThS, the proposed model will keep propagating. To this end, under the situation that the consistency term is high at a noisy region, if appropriate threshold values are chosen, the proposed model is able to remove the noise effect at that location.
Jonasson et al.'s (2005) model withand the proposed model withare applied to these synthetic data. The same initial contour (solid white curve in figure 2a) is used in both models, no convolution is applied because the tensor field inside the arch-shaped object is smooth except at two locations. ThS is set to the same value of 0.64 in both models, ThF equals 1.5 and time step size dt is 0.4, α=1.0. Implementation is worked on Matlab, using a PC with a 2.9 GHz processor and 3 GB RAM. It takes 117 and 81 steps in 3.1 and 2.0 s, respectively, for Jonasson et al.'s and the proposed model to stop propagation. It is observed from figure 2c,d that a hole exists in the result of Jonasson's model (figure 2c), while no hole shows up in the result of the proposed model (figure 2d). Results in other regions are similarly good. The overall accuracy rate of the proposed model compared with the ground truth is 98.25%, which is higher than the 94.18% accuracy rate of Jonasson et al. (2005). In words, the proposed model is able to obtain a more accurate result in a shorter time.
(b) Human brain results
Human brain diffusion-weighted magnetic resonance images were acquired using a 1.5 T GE Signa Neuro Vascular Interactive (NV/i) system at Nanjing Brain Hospital, China. Diffusion gradients were applied in 15 non-collinear directions and b=1000 s mm−2. The thickness of each slice was 3 mm without gaps, and 39 axial slices were acquired in parallel to the anterior commissure–posterior commissure (AC–PC) plane. Parameters are TE=104.4 ms; TR=800 ms; FOV is 24×24 cm2; and the size of the acquisition matrix is 128×128. We use the idea of Westin et al. (2002) to estimate the tensor field that is noisy in general.
Geometric flow models withandare applied to these human brain data to segment the corpus callosum. Convolution with a Gaussian kernel of zero mean and 0.45 variance is applied to both models. dt=0.1, α=0.05, ThS=0.4 and ThF=0.95. The same initial (figure 4a) is applied to both models. Figure 3 shows an intermediate and the final result from these two models. Figure 3a,b represents the results at steps 60 (intermediate) and 550 (final), respectively, of the model without consistency term, while figure 3c,d shows the results of the proposed model at step 60 (intermediate) and 120 (final), respectively. At step 60, there are two big gaps on the corpus callosum using the force without the consistency term, while there is only a small hole in the one with the consistency term. By using a consistency term, a good result has been obtained at step 120, while there is still a small hole in the model without consistency after 550 iterations. This shows that the model without the consistency term takes more than four times longer to obtain a result, which is still not as good as the one from the model with consistency measure. Figure 4 shows the initial and the final results of the above model with consistency measure from another view.
In figure 5, cross-sectional segmentation contours are superimposed on the corresponding FA images for slice 7 (figure 5a(i)–(iii)) and slice 10 (figure 5b(i)–(iii)). Figure 5a(i)b(i), a(ii)b(ii), a(iii)b(iii) present those of the models with β=0, 0.2 and 0.5, respectively. Bright regions correspond to high FA values. Diffusion in the corpus callosum region is highly anisotropic, and thus FA is higher there. A gap and a leak are observed on the result of the model without consistency term (see the bottom right and bottom left portions of slice 7 and 10, respectively, in figure 5a,b); when β is increased to 0.2, the results are better, but not as good as those when β is 0.5. This demonstrates that the model with consistency CONS2 and β=0.5 captures more accurate corpus callosum details than the model without consistency term.
Lastly, we apply the proposed model withto the same human brain dataset; the same initial as mentioned previously is applied to both models, results of which are shown in figure 6. Convolution is applied to both models. dt=0.1, α=0.05, ThS=0.45 and ThF=1.0. It takes more than three times for the model without any consistency term to obtain a result, which is still not as good as that obtained by using a consistency term (see bottom left and bottom right regions).
A new geometric flow-based method is proposed to segment diffusion tensor images. Consistency between the propagating front and the tensor field is incorporated into the evolution force to speed up the propagation, to compress noise effects and thus to obtain more accurate segmentation results in a shorter time. The proposed model using the new evolution force together with the new proposed termination criteria is able to handle images with abrupt change of tensor field caused by noise.
Experimental results for synthetic and human brain data confirmed the advantages of the model.
The idea of the proposed work can be extended to HARD MRI data. One approach is to use the HARD MRI data to approximate DTI data, based on which the proposed model can be applied directly. Another approach is to design a new consistency measure using the HARD MRI data directly.
The authors would like to thank the research centre for learning science of Southeast University at China for providing data. Yunmei Chen was partially supported by NIH R01NS052831-01 A1. Thanks also go to the reviewers and the editors for their valuable comments.
One contribution of 13 to a Theme Issue ‘Mathematical and statistical methods for diagnoses and therapies’.
- © 2008 The Royal Society