Probabilistic white matter fiber tracking using particle filtering and von Mises–Fisher sampling

Abstract

Standard particle filtering technique have previously been applied to the problem of fiber tracking by Brun et al. [Brun, A., Bjornemo, M., Kikinis, R., Westin, C.F., 2002. White matter tractography using sequential importance sampling. In: Proceedings of the ISMRM Annual Meeting, p. 1131] and Bjornemo et al. [Bjornemo, M., Brun, A., Kikinis, R., Westin, C.F., 2002. Regularized stochastic white matter tractography using diffusion tensor MRI, In: Proc. MICCAI, pp. 435–442]. However, these previous attempts have not utilised the full power of the technique, and as a result the fiber paths were tracked in a goal directed way. In this paper, we provide an advanced technique by presenting a fast and novel probabilistic method for white matter fiber tracking in diffusion weighted MRI (DWI), which takes advantage of the weighting and resampling mechanism of particle filtering. We formulate fiber tracking using a non-linear state space model which captures both smoothness regularity of the fibers and the uncertainties in the local fiber orientations due to noise and partial volume effects. Global fiber tracking is then posed as a problem of particle filtering. To model the posterior distribution, we classify voxels of the white matter as either prolate or oblate tensors. We then construct the orientation distributions for prolate and oblate tensors separately. Finally, the importance density function for particle filtering is modeled using the von Mises–Fisher distribution on a unit sphere. Fast and efficient sampling is achieved using Ulrich–Wood’s simulation algorithm. Given a seed point, the method is able to rapidly locate the globally optimal fiber and also provides a probability map for potential connections. The proposed method is validated and compared to alternative methods both on synthetic data and real-world brain MRI datasets.

Introduction

Diffusion tensor MRI (DTI) has become a popular tool for non-invasive exploration of the anatomical structure of the white matter in vivo (Basser et al., 1994). It endows each voxel with a $3\times 3$ symmetric positive-definite matrix, which characterises the local water diffusion process. It is based on a local Gaussianity assumption concerning the probability of water molecule motion in a defined time period. White matter fiber tracking or “tractography” estimates likely fiber paths by tracing the local tensor orientations (Mori and van Zijl, 2002, Parker, 2004). In this paper, we present a new and fast probabilistic fiber tracking algorithm which utilises the particle filtering technique and von Mises–Fisher sampling.

Broadly, the fiber tracking methods described in the literature can be classified as belonging to two groups. The first group of methods is based on local line propagation techniques (or streamline techniques) (Basser et al., 2000, Lazar et al., 2003, Mori et al., 1999). Step-by-step, they integrate a fiber pathway from a predefined seed point along the principal diffusion directions, which correspond to the principal eigenvectors of the diffusion tensors. The main difference among the methods in this group is the way in which local information is incorporated to locate smooth fiber paths. For instance, Lazar et al. (2003) use the entire diffusion tensor to deflect the estimated fiber trajectory in the desired directions. The main drawback of line propagation methods is that errors accumulate as the propagation takes place over a long distance.

The second group of methods are based on global optimisation techniques (Gossl et al., 2002, Parker et al., 2002, Pichon et al., 2005, Staempfli et al., 2006). Starting from a seed point, they attempt to locate an improved estimate of the true fiber pathway using energy minimisation techniques. For instance, Parker et al. (2002) apply the fast marching technique to propagate connection paths determined by the principal eigenvectors of the tensors. Gossl et al. (2002) apply Kalman filtering to track globally optimal paths, according to a fiber smoothness criterion. Pichon et al. (2005) determine the optimal path between two voxels by solving the Hamilton–Jacobi–Bellman equation using dynamic programming. Prados et al. (2006) use Riemannian geometry and control theory to trace the neural fiber bundles by computing the geodesic distances between seed and end point locations. More recently, Fletcher et al. (2007) develop a volumetric approach for quantitatively studying region-to-region white matter connectivity from diffusion tensor MRI. They use the Hamilton–Jacobi equation to formulate the minimal path problem between two regions.

One common feature of the above methods is that local fiber orientations are determined in a purely deterministic way. However, due to both noise (Macovski, 1996) and ambiguities for voxels where multiple fibers cross or branch (partial volume effects (Alexander et al., 2001)), the local fiber orientations measured by DTI are not completely reliable. For instance, in the case of partial volume effects, the diffusion process within voxels is no longer Gaussian. As a result, the diffusion tensor is an incomplete model of the diffusion signal. A natural way to deal with this problem is to measure uncertainty of the local fiber orientation measurement probabilistically at each voxel. Jones (2003) quantifies the uncertainties in the fiber orientations using the bootstrap method and then visualises the orientational field using uncertainty cones.

As a result, probabilistic fiber tracking methods have received considerable interest recently as a means of incorporating orientational uncertainties (Brun et al., 2002, Bjornemo et al., 2002, Behrens et al., 2003, Behrens et al., 2007, Friman et al., 2006, Lazar and Alexander, 2005, Parker et al., 2003, Parker and Alexander, 2003, Parker and Alexander, 2005). Instead of reconstructing fiber pathways deterministically, they aim to measure the probability of connectivity between brain regions.

These methods have two stages. In the first stage, they model the uncertainty in DTI measurements at each voxel using a probability density function (PDF) for the fiber orientations (Behrens et al., 2003, Friman et al., 2006). Behrens et al. (2003) were the first to formalise the PDF for local fiber orientations using a Bayesian framework. They present three models for describing the local diffusion process of water molecules, and their models use different levels of complexity. The parameters of the models are estimated using a Gibbs sampler. They applied a Markov Chain Monte Carlo (MCMC) procedure to sample from a single fiber orientation distribution. They then used the set of samples to model the uncertainty in orientation. Friman et al. (2006) proposed an alternative Bayesian method based on a simplified and more tractable diffusion tensor model. They replace the continuous PDF of fiber orientation by a set of uniform samples on a unit sphere. However, their PDF does not consider the uncertainties due to partial volume effects. Broadly speaking, the diffusion tensor allows the reliable definition of the principle diffusion direction, which corresponds to the local fiber orientation. However, in voxels with more complex configurations such as fiber crossings, the orientation conveyed by a diffusion tensor is ambiguous. To overcome this difficulty, a multi-tensor (or multi-compartment) extension of the single diffusion tensor model has been proposed (Tuch et al., 2002, Frank, 2001, Behrens et al., 2007). The basic idea here is to approximate the diffusion process using a mixture of Gaussian densities (Frank, 2001). Behrens et al. (2007) determine the complexity of the fiber structure at each voxel using automatic relevance determination. They then extend their previous work (Behrens et al., 2003) to deal with multi-orientations within a voxel by using a multi-compartment description of the diffusion process. More sophisticated methods for modeling PDFs of multi-fiber orientations at a voxel have also been developed for characterising the observed diffusion. For instance, Tuch (2004) describes the observed diffusion for high angular resolution diffusion imaging (HARDI) using the q-space framework, which is able to resolve fiber crossings. Tournier et al. (2004) express the diffusion weighted signal as a spherical convolution of the response function. The approximate PDF for fiber orientations is then recovered using spherical deconvolution. McGraw et al. (2006) model the PDF for fiber orientations using a mixture of von Mises–Fisher distributions. They use this model for segmenting high angular resolution diffusion MRI. Bhalerao and Westin (2007) build the PDFs from HARDI using a hyperspherical von Mises–Fisher mixture model. They map the fiber orientation samples to a 5D representation which can avoid the ambiguities associated with the sign flips of directions in 3D. They then fit their model using Maximum Likelihood estimation. Although these sophisticated models can better represent diffusion weighted signals, new imaging methods require much finer angular resolution and more scanning time than DTI (often 50 or more gradient directions are required). Therefore, in this paper, we focus on estimating the PDF of fiber orientations using the widely adopted diffusion tensor imaging model.

In the second stage, probabilistic tracking algorithms simply repeat a streamline propagation process (typically 1000–10000 times) with propagation directions randomly sampled from the PDF for fiber orientation. The fraction of the streamlines that pass through a voxel provide an index of the strength of connectivity between that voxel and the seed point. By contrast with the first stage, few methods have been developed to efficiently sample fiber paths from orientation PDF. The main difference between existing methods for probabilistic fiber tracking is found in the way that the PDF for fiber orientation is modeled. Most methods estimate the connectivity map by sampling directly from the PDF for fiber orientation. For instance, Parker et al. (2003) modeled the uncertainty of the fiber orientations using the normal distribution, and the parameters of the distribution is controlled in a heuristic way. They then sequentially sample fiber paths from the normal distribution. Parker and Alexander (2003) have further developed their previous work to deal with multi-fiber crossings. They distinguish between prolate and oblate tensors using the spherical harmonic parameters of the diffusion weighted signals. They fit the DWIs using both single and multi-tensor models. The normal distribution is used to control the uncertainties and sample fiber paths in a efficient way. Cook et al. (2004) further improved the method by capturing the orientation errors using the Watson distribution. The advantage of these methods is that they are able to efficiently sample fiber paths in a probabilistic way. Although the normal distributions are easily sample from, the choice is motivated by computational expediency and it is not clear that it accurately models the sources of uncertainty present. The main differences between these sequential sampling methods of Parker et al., 2003, Parker and Alexander, 2003, Cook et al., 2004 and our method are the weighting and resampling stages of particle filtering. On the other hand, the uncertainties of the fiber orientations have been modeled in a theoretically principled Bayesian framework by Behrens et al., 2003, Behrens et al., 2007, Friman et al., 2006. However, the fiber orientation distributions are difficult to simulate directly. As a result, the sampling process is not easily tractable. Thus it is necessary to resort to Markov Chain Monte Carlo methods (Behrens et al., 2003) or to evaluate the PDF discretely with low angular resolution. One drawback of these tracking methods is their computational complexity (often more than several hours on a high-end PC (Behrens et al., 2003, Friman et al., 2006)), and this is unacceptable in practice.

Particle filtering was originally investigated for probabilistic fiber tracking by Brun et al., 2002, Bjornemo et al., 2002. In a short abstract paper, Brun et al. (2002) were the first to sketch out the idea of applying the sequential importance sampling and resampling mechanism (also known as particle filtering) for tractography. Bjornemo et al. (2002) further detailed and developed the ideas in Brun et al. (2002). Specifically, they discussed in detail the conceptual model of importance sampling and resampling. They then construct the uncertainties of fiber orientations using a Gaussian noise model. Stochastic fiber paths are generated from the noisy distribution in a goal directed way. Finally, The connections between the proposed tracking model and the importance sampling and resampling framework are discussed. Since the weights of the sampled paths are constant throughout the propagation process, the weighting and resampling technique are not incorporated into their tracking model.

Developing on the work of Brun et al., 2002, Bjornemo et al., 2002, in this paper we present a fast and novel probabilistic method for white matter fiber tracking based on particle filtering and von Mises–Fisher sampling. We first formulate fiber tracking using a non-linear state space model and recursively compute the posterior distribution using particle filtering (Doucet et al., 2000, Doucet et al., 2001, Gordon et al., 1993). This technique has been successfully used in computer vision for visual tracking (Isard and Blake, 1998) and contour extraction (Perez et al., 2001). It provides a sound statistical framework for propagating a sample-based approximation of the posterior distribution. There are almost no restrictions on the type of model that can be used. As a result, we can model the fiber orientation distributions for prolate tensors and oblate tensors separately. The proposed tracking model can capture both smoothness regularity of the fibers and the uncertainties of the local fiber orientations due to both noise and partial volume effects. Since samples from the posterior path distribution are maintained at each propagation step, different decision criteria can be used to identify the optimal fiber. This procedure is similar to the active testing and tree pruning method for maximum a posteriori (MAP) road tracking developed by Geman and Jedynak (1996). Given a seed point, our method is able to rapidly locate the global optimal fiber path and also provide a probabilistic connectivity map between the seed point and all other voxel locations.

Compared to the methods proposed by Brun et al., 2002, Bjornemo et al., 2002, we make the following two novel contributions.

First, our proposed method generates reliable fiber paths by fully utilising the weighting and resampling technique of particle filtering. We model the probability density for fiber orientations in a theoretically justified way. In order to overcome the difficulty of sampling the theoretical fiber orientation distribution, we simulate paths using a simpler approximating distribution (namely the importance density function) for fiber orientations which is efficient to sample from. We then sequentially evaluate and adjust the sampled paths according to the true orientation distribution. In this way, the method achieves both efficiency and accuracy. Second, to implement an effective particle filtering algorithm for fiber tracking, we first apply the von Mises–Fisher distribution to model the prior and importance density of our tracking process. Fast sampling is realised on the unit sphere and fiber paths are efficiently generated using the simulation algorithm for the von Mises–Fisher distribution developed by Ulrich, 1984, Wood, 1994.

A preliminary version of the algorithm reported here was first described in a conference paper (Zhang et al., 2008). However, here we consolidate the work and expand the description of the method. Moreover, we provide additional qualitative and quantitative results to validate the method. The outline of the paper is as follows. In Section 2, we first formally develop the global fiber tracking model based on a non-linear state space model, and we also show how particle filtering can be used to recursively compute the posterior distribution, so as to compute optimal fiber paths and a probabilistic connectivity map. Section 3 provides the model ingredients necessary to implement the global tracking technique described in Section 2. Specifically, we describe how to construct the observation density, the prior density, the importance density function and the simulation for von Mises–Fisher distribution. In Section 4, we evaluate both qualitatively and quantitatively the performance of the algorithm on synthetic data and real-world diffusion MRI datasets. We also compare the results of our method with those obtained using alternative methods described elsewhere in the literature. Section 5 concludes the paper and discusses directions for future research.