Abstract
White matter fiber clustering aims to get insight about anatomical structures in order to generate atlases, perform clear visualizations, and compute statistics across subjects, all important and current neuroimaging problems. In this work, we present a diffusion maps clustering method applied to diffusion MRI in order to segment complex white matter fiber bundles. It is well known that diffusion tensor imaging (DTI) is restricted in complex fiber regions with crossings and this is why recent high-angular resolution diffusion imaging (HARDI) such as Q-Ball imaging (QBI) has been introduced to overcome these limitations. QBI reconstructs the diffusion orientation distribution function (ODF), a spherical function that has its maxima agreeing with the underlying fiber populations. In this paper, we use a spherical harmonic ODF representation as input to the diffusion maps clustering method. We first show the advantage of using diffusion maps clustering over classical methods such as N-Cuts and Laplacian eigenmaps. In particular, our ODF diffusion maps requires a smaller number of hypothesis from the input data, reduces the number of artifacts in the segmentation, and automatically exhibits the number of clusters segmenting the Q-Ball image by using an adaptive scale-space parameter. We also show that our ODF diffusion maps clustering can reproduce published results using the diffusion tensor (DT) clustering with N-Cuts on simple synthetic images without crossings. On more complex data with crossings, we show that our ODF-based method succeeds to separate fiber bundles and crossing regions whereas the DT-based methods generate artifacts and exhibit wrong number of clusters. Finally, we show results on a real-brain dataset where we segment well-known fiber bundles.
1. Introduction
Recent work shows that diffusion magnetic resonance imaging (dMRI) can help recovering complex white matter fiber bundles. However this is still an open problem due to the structural complexity of the fiber bundles, which can have crossing configurations. Diffusion tensor imaging (DTI) [1] is restricted in these conditions due to the hypothesis that the diffusion within a voxel follows a Gaussian distribution, a model that cannot model intravoxel crossings. Q-ball imaging (QBI) [2], a recent high-angular resolution diffusion imaging (HARDI) technique, overcomes this limitation by reconstructing a diffusion orientation distribution function (ODF), a spherical function that has its maxima agreeing with the underlying fiber populations. The ODF reconstruction from QBI is attractive because it is model-free and has been recently shown possible with a regularized and analytical solution [3], which produces a robust and very fast ODF reconstruction. In fact, the ODF estimation is, in practice, as fast as a standard least-square diffusion tensor (DT) estimation.
Efficient segmentation of fiber tracts in dMRI images is an important problem in neuroimaging problem because it has many potential applications. For example, it could potentially provide important information on diseases that affect fiber tracts. Alteration of the fiber tracts may provide new biomarkers in white matter pathologies and segmentation of these tracts can also improve our understanding of the functional role these tracts have and the cognitive consequences of their disruption.
The goal of this work is to provide a segmentation method that can separate the main white matter fiber bundles in the brain. We propose a new method that can segment fiber bundles and deal with fiber crossings while also requiring a minimum number of hypothesis from the data and a small number of algorithmic parameters. Spectral embedding and clustering methods have recently proved to be effective in image segmentation [4, 5]. However, classical approaches require restrictive hypotheses that are difficult to meet in real applications. For instance, N-Cuts [4] and Laplacian eigenmaps [6] require data within each cluster to be uniformly sampled, which produces artifacts when this hypothesis is not met. Moreover, classical approaches for image segmentation also assume that the scale within each cluster is the same using a single-scale parameter for the whole dataset. In order to overcome these limitations, we propose to use diffusion maps [7] as spectral embedding method. This method looses the dependence on the sampling of the elements to cluster. Moreover, we propose to use an adaptive scale-space parameter in order to deal with space-scale differences across different clusters. Finally, we propose two approaches to automatically determine the number of clusters by analyzing the spectra of the image embedding.
Another contribution of this paper is to show that the Q-ball ODF clustering using diffusion maps can reproduce the DT clustering using N-Cuts on simple synthetic images without crossings. On more complex data with crossings, we show that our method succeeds to separate fiber bundles and crossing regions on synthetic data, where the DT-based methods generate artifacts and exhibit wrong number of clusters. Finally, we successfully segment the fiber bundles in a real-human brain dataset in different regions with fibers crossing.
2. Methods
The main goal of this work is to produce a segmentation algorithm able to segment white matter fiber bundles from dMRI data. In order to represent intravoxel crossings with the ODF, we need at least 15 real coefficients when a spherical harmonic basis is used [3, 8, 9]. This leads to 3D images with a high dimensional element at each voxel. This high dimensionality makes previous diffusion imaging segmentation approaches based on level set methods such as [10–12] computationally expensive. Moreover, these methods require an initialization step. In order to perform the segmentation in an initialization-free manner and with a lower-dimensionality image, we use spectral clustering methods [4, 5] which perform dimensionality reduction before performing the segmentation and do not need initialization. The segmentation is then performed on the statistics within each cluster and the fiber crossings can be identified.
In this section, we present the three main parts of our algorithm. First, the estimation of the Q-ball diffusion ODF and its compact representation using spherical harmonics. Second, the metric used to measure distances between Q-ball ODFs. Last, the diffusion aps spectral clustering technique used to segment the ODF image into the background and the different fiber bundles.
2.1. ODF Estimation from QBI
QBI [2] reconstructs the diffusion ODF directly from the HARDI measurements on a single sphere by the Funk-Radon transform (FRT). In practice, the FRT value at a given spherical point is the great circle integral of the signal on the sphere defined by the plane through the origin with normal vector. The FRT is qualitatively illustrated in Figure 1. The ODF is intuitive because it has its maximum(a) aligned with the underlying population of fiber(s). However, computing statistics on a large number of discrete ODF values on the sphere is computationally heavy and infeasible to integrate into a segmentation algorithm of the whole brain. A more compact representation of the ODF is thus needed. In [3, 8, 9, 13] a simple analytic spherical harmonic (SH) reconstruction of the ODF is proposed. For completeness of the article, we now review and develop the main parts of our regularized analytical ODF reconstruction solution. The idea is to first estimate HARDI signal on the sphere with a regularized spherical harmonics approximation and then do a simple linear transformation of the harmonics to obtain the desired regularized ODF.
(a)
(b)
Spherical Harmonic (SH) Estimation of the HARDI Signal
The SH,
normally indicated by ( denotes the order and the phase
factor), are a basis for complex functions on the unit sphere. Explicitly, they
are given as follows:
where obey physics
convention () and is an
associated Legendre polynomial. For and , we define the new index and define our
modified basis with elements such that
where Re() and Img() represent the real and imaginary parts of , respectively.
The basis is designed to be symmetric, real, and orthonormal. Symmetry is
ensured by choosing only even order SH and the ratios in front of each term
also ensure that the modified basis is real and orthonormal with respect to the
inner product , where denotes
integration over the unit sphere and is the complex
conjugate of for and complex
functions on the sphere. We thus approximate the signal at each of the gradient
directions as
where is the number
of terms in the modified SH basis of order . Letting be the × 1 vector
representing the input signal for every encoding gradient direction, the × 1 vector of
SH coefficients , and is the × matrix constructed
with the discrete modified SH basis
We can write
the set of equations as an overdetermined linear system . We want to solve for the SH series coefficients , where . At this point, instead of simply evaluating the
integrals directly as done in [14] or performing a straightforward least-squared
minimization as in [15, 16],
we add local regularization directly into our
fitting procedure. This is to be able to
use a high-order estimation without
overmodeling the small perturbations due
to noise in the input diffusion MRI
signal. We thus define a measure, E, of the deviation from smoothness of a
function defined on the
unit sphere as , where is the
Laplace-Beltrami operator. Using the orthonormality of the modified SH basis,
where we have , the above functional can be
rewritten straightforwardly [3, 13] as
where is simply the matrix with
entries along the
diagonal ( is the order associated with the th coefficient,
that is, for ). We thus
obtain a closed-form expression for the regularization term. Therefore, the
quantity we wish to minimize can be expressed in matrix form as
where is the weight
on the regularization term. The coefficient vector minimizing this expression
can then be determined just as in the standard least-squares fit () from which we obtain the generalized expression
for the desired spherical harmonic series coefficient vector
From this SH
coefficient vector we can recover the signal on the Q-ball for any as Intuitively,
this approach penalizes an approximation function for having higher-order terms
in its modified SH series. This eliminates most of the higher-order terms due to
noise while leaving those that are necessary to describe the underlying
function. However, obtaining this balance depends on choosing a good value for
the parameter . We use the L-curve numerical method [17] and
experimental simulations to determine a good smoothing parameter [3, 13, 18].
Here, is used as in [3, 13, 18].
Analytical ODF Estimation
The true
diffusion orientation distribution function (ODF) in a unit direction , , is given by the radial projection of the probability
distribution function (PDF) of the diffusing water molecule. Tuch [2] showed
that this diffusion ODF could be estimated directly from the raw HARDI signal on a single
sphere of Q-space by the Funk-Radon transform (FRT) (Figure 1). In [3, 13], we showed how this FRT can be evaluated
analytically with an elegant corollary to the Funk-Hecke theorem [19]. The
final ODF reconstruction on the sphere then becomes a simple linear
transformation of the SH coefficients describing the
input HARDI signal ,
where are the SH
coefficients describing the ODF and because is always even
in our modified SH basis. We see that the SHs are eigenfunctions of the
Funk-Radon transform with eigenvalues depending only on the order of the SH
series.
Hence, by using an SH estimation of the HARDI signal,
we have showed that the QBI can be solved analytically. This was also showed in
[8, 9]. An important contribution in favor of our approach is that this solution
can be obtained while imposing a well-defined regularization criterion. The
accuracy of the modified SH series approximation with the Laplace-Beltrami
smoothing was established in [18] and our regularized ODF solution was also
shown to have better fiber detection properties and shown to be more robust
to noise than similar solutions [8, 9].
2.2. Distances between ODFs
Once the ODF are computed, we want to capture similarities and dissimilarities between two ODFs, that is, two spherical functions that can be represented by real-SH vectors of length , and , as shown in (8) in the previous section. Since the ODFs come from real physical diffusion measurements they are bounded and form an open subset of the space of real-valued spherical functions with an inner product defined as Again, because of the orthonormality of the spherical harmonic basis, the cross-terms cancel and the expression is simply Therefore, the induced norm giving us the distance metric between two ODFs is
The Euclidean distance was also used successfully for ODF segmentation in [12] and for DTI segmentation in [11] even though more appropriate metrics exist such as the J-divergence [11, 20] and Riemannian geodesic distances [11]. Similarly, one can think of choosing another metric to compare ODFs. For instance, since the ODF can be viewed as a probability distribution function (PDF) of fiber orientations, one can use the Kullback-Leibler distance between two PDFs, as done in [2]. However, in that case the problem quickly blows up computationally because one needs to use all discrete HARDI data on the sphere instead of the SH coefficients ().
2.3. Diffusion Maps-Based Clustering
We now want to segment white matter fiber bundles in a Q-ball image. One of the open questions in Q-ball image analysis and clustering is that which metric should be used to compare Q-ball ODFs. Here, we describe a clustering algorithm that infers an embedding and a metric to compare ODF images. We derive an affinity measure incorporating the Euclidean distance and the spatial location distance between ODFs. This affinity measure then used in a spectral embedding framework. As mentioned in [7], the Euclidean distance within this embedding actually represents an intrinsic metric of the data, which can be used to perform statistics in the embedded space and can thus be used to segment Q-ball ODF images into white matter fiber bundles.
Spectral Embedding and Clustering
In recent
years, spectral manifold learning and clustering techniques [4, 6, 21–23]
have become one of the most popular modern clustering
family of methods. They
are simple to implement,
they can be solved efficiently by standard linear
algebra software, and they very often outperform
traditional manifold learning
and clustering algorithms such as the classical principal component
analysis
(PCA) [24] and -means [25] algorithms. Moreover, due to the dimensionality
reduction properties, they are especially well suited to work with
high-dimensional data. These techniques have been recently used to cluster
various types of images [4, 5] and white matter fiber tracts [26]. In our case,
we perform the spectral clustering for two different types of elements: the DT
and the ODF. In the DT case, the element is represented by a -dimensional
vector corresponding to the upper (or lower)
triangular part of the DT
symmetric matrix. In the ODF case, the element is represented by the -dimensional
vector corresponding to the 4th-order spherical
harmonic ODF estimation.
Spectral clustering reduces the clustering problem to a graph partitioning problem. Each element to be clustered is represented as a node in a graph and the edges joining the vertex are a measure of affinity between the elements. This affinity measure lies between and , being the less affine case. A spectral decomposition of this graph is taken by calculating the eigenvalue decomposition (EVD) of the graph Laplacian [27]. Then a low-dimensional Euclidean manifold embedding is inferred from this decomposition. Finally, the clustering is performed in the inferred Euclidean manifold.
All the above techniques rely on three hypotheses.
(1)Preservation of the distance relationship: after a distance is defined between elements, the learned manifold should preserve the distance relation.(2)Uniform sampling of the elements: the density of the extracted elements changes if and only if these elements belong to anatomically different bundles.(3)Convexity of the elements: if two elements are in the dataset, almost all of the intermediate tracts obtained by the interpolation that can be inferred from the metric used to build the affinity matrix are in the dataset.
It is not easy to guarantee that the data to be embedded and clustered will adhere to these hypotheses. Donoho and Grimes, in [13], analyze when a spectral embedding algorithm is able to recover the true parameterization of a set of images. As medical images represent the discretization of a continuous space, hypotheses 1 and 3 are plausible. However, there is no indication that within a fiber bundle the distribution of the elements (DT or ODF) are uniformly sampled. Moreover, in [29], it is shown that different sampling frequencies within one cluster leads the N-Cuts and Laplacian eigenmaps methods to subdivide the cluster in several parts. In order to overcome this limitation and to be resilient to sampling frequency differences within a cluster, we use the diffusion maps [7] spectral embedding technique. We now describe the three steps involved in the diffusion maps algorithm in turn.
Step1 (Computing the Affinity Matrix). Letting represent the set of all ODF elements to cluster, the main idea is to look for a representation between the elements of that is more representative than (recall that ODFs are ) and reduces the dimensionality of the problem. With keeping this in mind, a fairly good way of representing any set of elements with an affinity function , is a weighted graph, , where the weight of the edge between two vertices represents the affinity of the elements connected by this edge. More formally, for an edge,1, the weight of the edge is . Hence, each element of the adjacency matrix of or conversely the affinity matrix of is Taking this in account, the weighted graph can be also noted as .
Usually, a distance function instead of an affinity function is given. The distances can be easily converted into affinities by applying a kernel to the distance function where is an adaptive scale-space parameter that may depend on the elements and . In this work, the adaptive scale-space parameter is taken following [30]. A “neighbor-number” is given as parameter to the algorithm and then , where is the th closest neighbor according to the distance function of element . This choice of a scaling parameter for each point allows self-tuning of the point-to-point distances according to the local statistics of the neighborhoods surrounding points and .
As in image segmentation, the spatial position of each element is important, the spatial dependency should be incorporated into the affinity matrix. Following [5, 31], we use Markovian relaxation to incorporate this information. In order to represent the affinity of all the elements that can be reached within one spatial step, the affinity matrix is modified in the following way: where are spatial coordinates of element in the image where is a diagonal matrix with , usually called the row-sum or degree matrix of .
Then, obtaining the affinities of elements that can be reached within spatial steps is enough to elevate to the power of , as stated in [31]. Moreover, can be chosen to be the smallest positive integer which results in nonzero elements in the whole matrix in order to represent the weakest connected induced graph. The diagonal adjustment forces the inherent random walk to a uniform steady state, hence every part of the Markov field will be explored at the same speed. For the sake of clarity, will be referred to as affinity matrix in the rest of the paper.
Step 2 (Performing the Embedding). The algorithm
must embed the elements of into an -dimensional
Euclidean space . This is done by applying eigenvalue decomposition to
the Laplacian of the affinity matrix. This embedding must be compliant with
hypothesis 1. As in [6, 7, 27], this is done by performing the
spectral decomposition of the graph Laplacian of the graph induced by ,
where is number of
elements to be clustered.
In order to overcome the necessity of hypothesis 2,
we prenormalize the affinity matrix as done in
[7].
This is done by normalizing the weight of each edge of the graph,
, by the probability density of both elements relating
through the edge,
where , the probability density function of the elements in , is not known but can be approximated up to a
multiplication factor by
Due to the necessity of having a uniform behavior of
the clustering algorithm without minding the scale of the affinity measure
taken, a doubly stochastic matrix normalization is performed:
As is a double
stochastic symmetric matrix, the eigenvalue decomposition
of (16) can
be calculated by taking the singular value decomposition (SVD)
Finally, the Euclidean coordinates of an element in the -dimensional
embedding manifold are
where
is the
eigenvector column matrix and the corresponding eigenvalues are, . The first eigenvector is not taken into account as a component in the
embedding because it is constant and hence meaningless as shown in [6, 7, 27].
Step 3 (Clustering). Once the embedding has been performed, several techniques have been proposed for the clustering step [4, 6, 32].
The first step in this process is to determine the number of clusters, this can be done in two ways. The first, as in [33], is choosing the number of clusters according to the “elbow”. This is present in the eigenvalue plot. For instance, if the slope of the eigenvalue plot changes noticeably at eigenvector , the number of clusters should be . The second way is reordering the affinity matrix rows and columns following the second eigenvector as proved in [34], which shows the block structure of the matrix as squared blocks along the matrix diagonal. Then, the number of clusters is the number of blocks. Their commended number of dimensions for the embedding is the same as the number of clusters. Finally, the clustering is performed by running a -means clustering algorithm on the space spanned by . A formal justification for this approach can be found in [6, 32].
2.3. Q-Ball Data Generation and Acquisitions
Synthetic Data
We generate
synthetic HARDI data using the multitensor model which is simple and leads to
an analytical expression of the ODF [2, 18]. For a given -factor and
noise level, we generate the diffusion-weighted signal
where is the th gradient direction on the
sphere, is the number
of fibers, and is the volume
fraction of each fiber. In practice, we use from a 3rd-order tessellation of the
icosahedron, ?s/mm2, and or . is the
diffusion tensor with standard eigenvalues ×10-2?mm2/s oriented in direction , which agree with reported physiological values [35]. Finally, we add complex Gaussian noise with standard
deviation of , producing a signal with signal-to-noise ratio of 35.
We generate three synthetic data example, two simple examples: one with a ring of sinusoidal-shaped fibers, one with fibers with different sizes and scales, and the other with complex crossing areas simulating the “U”-fibers (corticocortical fibers) that can occur in the brain. These synthetic datasets help understand the behavior of the different spectral clustering methods when confronted with simple and complex fiber geometries.
Human Brain Data
Diffusion-weighted data and high-resolution T1-weighted images
were acquired on a
whole-body 3 Tesla Magnetom Trio scanner
(Siemens, Erlangen) equipped with an
8-channel head array coil
[36].
The spin-echo echo-planar-imaging sequence, ?ms, ?s, 128 × 128 image matrix,
FOV = 220 × 220?mm2,
consists of 60 diffusion
encoding gradients [37]
with a -value of 1000?s/mm2. Seven images without any
diffusion weightings are placed at the beginning
of the sequence and after each
block of 10 diffusion-weighted images as anatomical
reference for offline
motion correction. The measurement of 72 slices
with 1.7?mm2 thickness (no gap)
covered the whole brain. Random variations
in the data were reduced by
averaging 3 acquisitions,
resulting in an acquisition time of about 45 minutes.
No cardiac gating was employed to limit
the acquisition time. The issue of
cardiac gating is discussed in
[38].
Additionally, fat saturation was employed
and we used 6/8 partial Fourier imaging,
a Hanning window filtering, and
parallel acquisition (generalized autocalibrating
partially parallel
acquisitions, reduction factor = 2) in the axial plane.
The brain is peeled from the T1-anatomy, which was aligned with the Talairach stereotactical coordinate system [39]. The 21 images without diffusion weightings distributed within the whole sequence were used to estimate motion correction parameters using rigid-body transformations [40], implemented in [41]. The motion correction for the 180 diffusion-weighted images was combined with a global registration to the T1 anatomy computed with the same method. The gradient direction for each volume was corrected using the rotation parameters. The registered images were interpolated to the new reference frame with an isotropic voxel resolution of 1.72?mm2 and the 3 corresponding acquisitions and gradient directions were averaged.
Distance Functions between Elements to Cluster
In order to
implement the diffusion maps
spectral clustering method a distance function for
each data type is chosen. This distance
function is used to calculate the
affinity matrix as expressed by
(13).
In the DT case, following [42],
we use the Riemannian tensor distance. In the ODF case we use the distance
shown in (11).
3. Results
3.1. Synthetic Data Experiments
Diffusion Maps versus N-Cuts
The first
experiment shows the difference in performance
between the diffusion maps and
N-Cuts approach. The N-Cut algorithm does not
perform the sampling-based
normalization described by
(17) and is thus
sensitive to sampling frequency differences
within the clusters. In order to
show this sampling hypothesis problem,
we used a ring fiber bundle with
different sampling frequencies. Within the
ring, the fibers have a sinusoidal
shape and the frequency of the modulating sine
function is times bigger in
the lower half of the ring. More formally,
the fibers follow the angular
function , where for the upper
half of the ring and for the lower
half. Two clusters are expected, the ring and the background. The results of
both clustering techniques are shown in Figure 2, where the background
has been masked out. Figure 2(a)
shows the plot of the first 10
eigenvalues for the N-Cuts method, shown in
Figures 2(b) and 2(c), and
the slope between the line joining and and the line
joining and changes
drastically. This elbow at indicates that
there are 2 clusters. Figure 2(d)
shows the plot of the first 10
eigenvalues for the diffusion maps method whose clustering results are shown in
Figure 2(e). The N-Cuts exhibits frequency-dependent clustering
artifacts while the diffusion maps method clearly shows two clusters. In the
diffusion maps, the clustering has correctly
segmented the background and the
ring.
(a) N-Cuts eigenvalue plot
(b) N-Cuts, 2 clusters (blue and black)
(c) N-Cuts, 3 clusters (blue, orange, black)
(d) Diffusion maps eigenvalue plot
(e) Diffusion maps, 2 clusters (orange and black)
ODF versus DT Images
In Figure 3, a single
fiber scenario with no fiber crossing is
shown. The DT-based and ODF-based image clustering
produce the same results.
Hence, ODF clustering reproduces DT-based
results on a simple fiber population
example.
Finally, Figure 4
shows a fiber crossing
scenario with two overlapping fiber bundles
that have different geometries.
Segmentation was performed over the DT
and the ODF image shown in
Figure 5.
Note that the cluster number is correctly estimated only in
the ODF image. Moreover, the ODF N-Cuts
segmentation exhibits artifacts not
present in the ODF diffusion maps segmentation.
The ODF diffusion maps
effectively identify the two different
fiber bundles as well as the fiber
crossing areas.
(a) DTI, 3 clusters (purple, green, black)
(b) ODF, 3 clusters (orange, blue, black)
(a) DTI ellipsoids
(b) ODF spherical functions
(a) DTI, 4 clusters (purple, green, blue, black) 2
(b) ODF, 4 clusters (black, purple, blue, green)
3.2. Real Data
The real-data experiment presented in this section shows the segmentation and labeling of a cropped axial and coronal slice. The cropped slices were chosen by an expert in regions of known fiber crossings where the DT model is normally limited. The ROIs show intersection of several fiber bundles. Hence, our segmentation algorithm is confronted with elements that have different orientation and different diffusion characteristics.
In order to show that ODF data segments the white matter fiber bundles better than the DT data in real cases, we analyze the evolution of the affinity matrix as the scale-space parameter changes in the axial cropped slice shown in Figure 6. Affinity matrices were computed with varying scale-space parameter between , , , and of the quantity of elements () to cluster, respectively. In order to show the block structure of the affinity matrices, they were reordered using the second (Fiedler) biggest eigenvector [34]. It can be seen in Figure 7 that as the scale diminishes, the DT data shows a high correlation between all the elements of the slice. This makes clustering very difficult because the blocks are small and highly correlated. On the other hand, the ODF data shows a very clear block structure across all scales. This block structure shows a high correlation of the elements within each block and a low interblock correlation, giving a much better input to the clustering algorithm than the DT data.
(a) Axial GFA slice with axial slice marked
(b) Coronal GFA slice with cropped region marked
(a) DTI reordered affinity matrices as the scale-space parameter decreases 0
(b) ODF reordered affinity matrices as the scale-space parameter decreases
In Figure 6, the location of the cropped axial slice is shown in the axial slice, Figure 6(a), and coronal slice, Figure 6(b). As it can be seen in the segmented and labeled axial slice, Figure 8, the segmentation also allows to identify and label some of the main white matter structures, Corpus Callosum (CC), Anterior Corona Radiata (ACR), Forceps Major (fmajor) and Forceps Minor (fminor).
In Figure 9, the location of the cropped coronal slice is shown in the axial slice, Figure 9(a), and coronal slice, Figure 9(b). As it can be seen in the segmented and labeled coronal slice, Figure 9(c), the segmentation allows to identify and label main white matter structures: Corpus Callosum (CC), Cingulum (CG), Corona Radiata (CR), Superior Longitudinal Fasciculus (SLF). Note that the segmentation is resilient to crossing areas such as seen at the interface between CR and CC.
(a) Axial GFA slice with coronal slice marked
(b) Coronal GFA slice with cropped region marked
(c) Labeled ODF visualization of the cropped region
4. Discussion
We have presented an algorithm to perform Q-ball imaging segmentation of white matter fiber bundles. The proposed method combines state-of-the-art HARDI reconstruction and state-of-the-art spectral clustering techniques. Our algorithm is initialization-free and has only two parameters. A scale-space parameter and the number of regions (clusters) are to be found. Regarding this number of clusters parameter, we have proposed to estimate it automatically. We have introduced a spectral embedding technique that does not require uniform sampling of the elements. To do so, the affinity measure used incorporates an Euclidean distance measure between the spherical harmonic coefficients describing the Q-ball ODFs and also incorporates the spatial location distance between ODFs. The affinity measure and the metric induced in the embedded space is then used to cluster Q-ball ODF images into multilabel segmentation representing the fiber bundles. Spectral embedding has already been applied to dMRI (e.g., [5]). However, to our knowledge, this is the first work using the diffusion maps that avoids the high dependence on element sampling. It is also the first work attempting Q-ball ODFs.
We have illustrated that the ODFs are the desirable elements to use for clustering in the white matter because the classical DT model is limited in regions of fiber crossings. The ODF is even more attractive because of the recent analytical spherical harmonic solution to the ODF reconstruction [3, 8, 9, 13]. The analytical solution is in fact as fast as a standard DT least-square estimation. In this work, we believe that we have used the state-of-the-art ODF reconstruction method [13], which is regularized, robust and very simple to implement.
The spectral embedding performed by the diffusion maps technique is at the heart of our segmentation algorithm. Whereas other spectral embedding techniques have a tendency to produce artifacts in the presence of different sampling characteristics within a cluster, the technique used in this work greatly reduces this tendency by performing the simple linear algebra calculation shown in (17).
Spectral embedding techniques produce a representation of the embedded data based on element-to-element affinities. This leads to the fundamental issue: how to choose the affinity measure? It is a challenge to find a measure that incorporates similarities between elements as well as the spatial location difference between elements. For similarities between elements, we chose the Euclidean distance between spherical harmonic coefficients describing the ODFs. This approach is simple and very efficient because it allows to process the ODFs directly on the SH coefficients. The Euclidean distance has also been used successfully in a level set segmentation framework [12] and it would be interesting to compare our spectral clustering approach against it. For spatial location difference, we chose Markovian relaxation in order to be consistent with the graph theoretical representation of the diffusion maps technique. Although this way of representing the distance involves an artificial elimination of all the nonneighboring relations of the ODF elements in the affinity matrix and an adjustment of the diagonal elements, we believe that the resulting affinity relations represent the affinity better. The affinity of two neighboring elements at the beginning of the Markovian relaxation algorithm is represented by a function of the Euclidean distance between them. This affinity can be interpreted as the probability that a random walker has of going from the first element to the second. The affinity of two elements at the end of the relaxation is the probability of a random walker starting from one element and reaching the second in a certain number of steps.
The final step of our algorithm is -means clustering. We believe that there is room for improvement in this last part of the algorithm. In the first place, the -means algorithm needs an explicit number of clusters to find. This can be heuristically determined by analyzing the eigenvalue plot or the reordered affinity matrix structure, as shown in this work. However, an automatic method that could find the number of clusters would considerably improve the algorithm. In the second place, the -means algorithm and its variants, for instance, k-medians, k-medioids, search for isotropic clusters in the embedding space [25]. These methods are able to perform clustering on convex structures. This could also improve the last clustering phase of our algorithm
Finally, in order to analyze the importance of the difference between our diffusion maps algorithm and the widely used N-Cuts, we used synthetic simulations. In these simulations, we generated a synthetic image with a single cluster within which the sampling of the elements changed. We showed that when this sampling changes, the N-cuts algorithm produces artifacts while our diffusion maps method does not. As uniform sampling within a cluster is a difficult property to guarantee in the white matter fiber bundles, our diffusion maps method is better suited for this task. We thus believe that diffusion maps are the right clustering method to be used on dMRI processing problems.
5. Conclusions
In this work, we have presented two contributions. First, we have shown that in order to perform spectral clustering on complex dMRI with crossing fiber bundles, an HARDI technique such as Q-ball imaging is better than the classical DTI technique. This is because the ODF reconstructed from QBI is able to recover multiple crossing fiber populations. Second, a diffusion maps-based technique for image segmentation was introduced to reduce artifacts arising from the widely used N-Cuts image segmentation. We have illustrated the advantages of the ODF diffusion maps segmentation algorithm, and showed on a real dataset that our algorithm is able to identify important and complex white matter fiber bundles.
Finally, the diffusion maps technique has been shown to be more robust to sampling frequency variations within each object to be segmented. In order to cluster the elements, we have used an adaptive scale-space parameter and we have used Markovian relaxation in order to incorporate spatial dependencies. Overall, the approach is theoretically sound with the graph-based representation which lies at the heart of spectral clustering methods.
Therefore, we have an algorithm to perform fiber bundle clustering for a single brain. It is now important to study the behavior over several subjects in order to assess the reproducibility of the algorithm. In time, this will enable to perform multisubject statistics within bundles in the embedded space. This will help characterize the white matter fiber bundles of several subjects and study if the alteration of these segmented tracts can provide new biomarkers for white matter diseases.
Acknowledgments
The authors would like to thank the Max Planck Institute for Human Cognitive and Brain Sciences (Leipzig, Germany) and, in particular, to Timm Wetzel for providing the diffusion-weighted MR datasets and Alexandre Gramfort and Marc Niethammer for useful comments. Part of this work was supported by PAI Procope and the Diffusion MRI ARC. Demian Wassermann would like to acknowledge CONICET (Argentina).