Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2893231

Formats

Article sections

Authors

Related links

Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2010 June 28.

Published in final edited form as:

Med Image Comput Comput Assist Interv. 2009; 12(Pt 1): 894–902.

PMCID: PMC2893231

NIHMSID: NIHMS208694

See other articles in PMC that cite the published article.

We describe a technique to simultaneously estimate a weighted, positive-definite multi-tensor fiber model and perform tractography. Existing techniques estimate the local fiber orientation at each voxel independently so there is no running knowledge of confidence in the estimated fiber model. We formulate fiber tracking as recursive estimation: at each step of tracing the fiber, the current estimate is guided by the previous. To do this we model the signal as a weighted mixture of Gaussian tensors and perform tractography within a filter framework. Starting from a seed point, each fiber is traced to its termination using an unscented Kalman filter to simultaneously fit the local model and propagate in the most consistent direction. Further, we modify the Kalman filter to enforce model constraints, *i.e.* positive eigenvalues and convex weights. Despite the presence of noise and uncertainty, this provides a causal estimate of the local structure at each point along the fiber. Synthetic experiments demonstrate that this approach significantly improves the angular resolution at crossings and branchings while consistently estimating the mixture weights. *In vivo* experiments confirm the ability to trace out fibers in areas known to contain such crossing and branching while providing inherent path regularization.

The advent of diffusion weighted magnetic resonance imaging has provided the opportunity for non-invasive investigation of neural architecture. Using this imaging technique, neuroscientists can investigate how neurons originating from one region connect to other regions, or how well-defined these connections may be. For such studies, the quality of the results relies heavily on the chosen fiber representation and the method of reconstructing pathways.

To begin studying the microstructure of fibers, we need a model to interpret the diffusion weighted signal. Such models fall broadly into two categories: parametric and nonparametric. One of the simplest parametric models is the diffusion tensor which describes a Gaussian estimate of the diffusion orientation and strength at each voxel. While robust, this model can be inadequate in cases of mixed fiber presence or more complex orientations, and so to handle more complex diffusion patterns, various alternatives have been introduced: weighted mixtures [1,2,3,4], higher order tensors [5], and directional functions [6]. In contrast, nonparametric techniques estimate an orientation distribution function (ODF) describing an arbitrary configuration of fibers. For this estimation, several techniques have been proposed, among them Q-ball imaging [2], spherical harmonics [7,8], spherical deconvolution [9,10,11,6], and diffusion orientation transforms [12].

Based on these models, several techniques can be used to reconstruct pathways. Deterministic tractography using the single tensor model simply follows the principal diffusion direction, while multi-fiber models use various techniques for determining the number of fibers present or when pathways branch [3,13]. While parametric methods directly describe the principal diffusion directions, interpreting the ODFs from model independent representations typically involves a separate algorithm to determine the number and orientation of diffusion patterns present [14,9,8,15]. Several filtering approaches have been proposed. For example, Kalman and particle filters [16,17,18], as well as a moving least squares approach [19], have been used with single tensor streamline tractography, but these have been used for path regularization and not to estimate the underlying fiber model. One approach has used a linear Kalman filter, although this method was applied to estimate each voxel independently during acquisition [20].

Of the approaches listed above, nearly all fit the model at each voxel independent of other voxels; however, tractography is a causal process: we arrive at each new position along the fiber based upon the diffusion found at the previous position. In this paper, we treat model estimation and tractography as such by placing this process within a causal filter. As we examine the signal at each new position, the filter recursively updates the underlying local model parameters, provides the variance of that estimate, and indicates the direction in which to propagate tractography.

To begin estimating within a finite dimensional filter, we model the diffusion signal using a weighted mixture of two tensors. This enables estimation directly from the raw signal without separate preprocessing or regularization. Because the signal reconstruction is nonlinear, we use the unscented Kalman filter to perform local model estimation and then propagate in the most consistent direction (Fig. 1). Further, we use a constrained version of the unscented Kalman filter to ensure the tensor eigenvalues are positive and the mixture weights are non-negative and convex. Using causal estimation in this way yields inherent path regularization, consistent partial volume estimation, and accurate fiber resolution at crossing angles not found with independent optimization.

Section 2.1 provides the necessary background on modeling the measurement signal using tensors and defines the specific weighted two-fiber model employed in this study. Then, Section 2.2 describes how this model can be estimated using an unscented Kalman filter and further how the constraints are enforced.

In diffusion weighted imaging, image contrast is related to the strength of water diffusion, and our goal is to accurately relate these signals to an underlying model of putative fibers. At each image voxel, diffusion is measured along a set of distinct gradients, **u**_{1},…, **u**_{m}^{2} (on the unit sphere), producing the corresponding signal, **s** = [*s*_{1}, …, *s _{m}*]

From that general mixture model, we choose a restricted form with only two weighted components. This choice is guided by several previous studies which found two-component models to be sufficient at *b* = 1000 [2,3,13,14,4,21]. Also, we assume the shape of each tensor to be ellipsoidal, *i.e.* there is one dominant principal diffusion direction **m** with eigenvalue *λ*_{1} and the remaining orthonormal directions have equal eigenvalues *λ*_{2} = *λ*_{3} (as in [4,6]). These assumptions leave us with the following model used in this study:

$${s}_{i}={s}_{0}{w}_{1}{e}^{-b{\mathbf{\text{u}}}_{i}^{T}{D}_{1}{\mathbf{\text{u}}}_{i}}+{s}_{0}{w}_{2}{e}^{-b{\mathbf{\text{u}}}_{i}^{T}{D}_{2}{\mathbf{\text{u}}}_{i}},$$

(1)

where *w*_{1}, *w*_{2} are convex weights and *D*_{1}, *D*_{2} are each expressible as *D* = *λ*_{1}**mm*** ^{T}* +

Given the measured signal at a particular voxel, we want to estimate the underlying model parameters that explain this signal. As in streamline tractography, we treat the fiber as the trajectory of a particle which we trace out. At each step, we examine the measured signal at that position, use that measurement to update our model parameters within the filter, and propagate forward in the most consistent direction. Fig. 1 illustrates this filtering process.

To use a state-space filter for estimating the model parameters, we need the application-specific definition of four filter components:

Algorithm 1. Unscented Kalman Filter |

1: Form weighted sigma points
${\mathbf{\text{X}}}_{t}={\{{w}_{i},{\mathbf{\text{x}}}_{i}\}}_{i=0}^{2n}$ around current mean x and covariance _{t}P with scaling factor _{t}ζ |

${\mathbf{\text{x}}}_{0}={\mathbf{\text{x}}}_{t}\phantom{\rule{2em}{0ex}}{\mathbf{\text{x}}}_{i}={\mathbf{\text{x}}}_{t}+{[\sqrt{\zeta {P}_{t}}]}_{i}\phantom{\rule{2em}{0ex}}{\mathbf{\text{x}}}_{i+n}={\mathbf{\text{x}}}_{t}-{[\sqrt{\zeta {P}_{t}}]}_{i}$ |

2: Predict the new sigma points and observations |

${\mathbf{\text{X}}}_{t+1|t}=f[{\mathbf{\text{X}}}_{t}]\phantom{\rule{2em}{0ex}}{\mathbf{\text{Y}}}_{t+1|t}=h[{\mathbf{\text{X}}}_{t+1|t}]$ |

3: Compute weighted means and covariances, e.g. |

${\overline{\mathbf{\text{x}}}}_{t+1|t}=\text{\u2211}_{i}{w}_{i}{\mathbf{\text{x}}}_{i}\phantom{\rule{2em}{0ex}}{P}_{x\mathit{\text{y}}}=\text{\u2211}_{i}{w}_{i}({\mathbf{\text{x}}}_{i}-{\overline{\mathbf{\text{x}}}}_{t+1|t}){({\mathbf{\text{y}}}_{i}-{\overline{\mathbf{\text{y}}}}_{t+1|t})}^{T}$ |

4: Update estimate using Kalman gain K and scanner measurement y_{t} |

${\mathbf{\text{x}}}_{t+1}={\overline{\mathbf{\text{x}}}}_{t+1|t}+K({\mathbf{\text{y}}}_{t}-{\overline{\mathbf{\text{y}}}}_{t+1|t})\phantom{\rule{2em}{0ex}}{P}_{t+1}={P}_{xx}-K{P}_{yy}{K}^{T}\phantom{\rule{2em}{0ex}}K={P}_{xy}{P}_{yy}^{-1}$ |

- The system state (
**x**): the model parameters - The state transition function (
*f*): how the model changes as we trace the fiber - The observation function (
*h*): how the signal appears given a particular model state - The measurement (
**y**): the actual signal obtained from the scanner

For our state, we directly use the parameters for the two-tensor model in Eq. 1:

$$\mathbf{\text{x}}={\left[{\mathbf{\text{m}}}_{1}\phantom{\rule{0.3em}{0ex}}{{\lambda}_{1}}_{1}\phantom{\rule{0.3em}{0ex}}{{\lambda}_{2}}_{1}\phantom{\rule{0.3em}{0ex}}{w}_{1}\phantom{\rule{0.3em}{0ex}}{\mathbf{\text{m}}}_{2}\phantom{\rule{0.3em}{0ex}}{{\lambda}_{1}}_{2}\phantom{\rule{0.3em}{0ex}}{{\lambda}_{2}}_{2}\phantom{\rule{0.3em}{0ex}}{w}_{2}\right]}^{T},\phantom{\rule{1em}{0ex}}\mathbf{\text{m}}\in {\mathbb{S}}^{2},\lambda \in {\mathbb{R}}^{+},w\in [0,1].$$

(2)

For the state transition we assume identity dynamics; the local fiber configuration does not undergo drastic change as it moves from one location to the next. Our observation is the signal reconstruction, **y** = *h*[**x**] = **s** = [*s*_{1},…, *s _{m}*]

Since our signal reconstruction in Eq. 1 is nonlinear, we employ an unscented Kalman filter to perform estimation. Similar to classical linear Kalman filtering, the unscented version seeks to reconcile the predicted state of the system with the measured state and addresses the fact that these two processes–prediction and measurement–may be nonlinear or unknown. In Algorithm 1 we present the standard version of this filter; for more thorough treatments, see [22,23]. It is important to note that while particle filters are a common approach to nonlinear estimation, we chose instead the unscented Kalman filter primarily for its low computational complexity. With respect to state dimension, particle filters require the number of particles to be exponential to properly explore the state space. In contrast, the unscented filter requires 2*n* + 1 particles (sigma points) for a Gaussian estimate of the *n*-dimensional state.

In this standard formulation, we have ignored the constraints on our model. This results in instabilities: the diffusion tensors may become degenerate with zero or negative eigenvalues, or the weights may become negative. To enforce appropriate constraints, one can directly project any unconstrained state **x** onto the constrained subspace [23]. In other words, we wish to find the state **x^** closest to the unconstrained state **x** which still obeys the constraints, *A***x^** ≤ **b**. Using *P _{t}* as a weighting matrix, this becomes a quadratic programming problem:

$$\underset{\widehat{\text{x}}}{\text{min}}\phantom{\rule{0.2em}{0ex}}{(\mathbf{\text{x}}-\widehat{\mathbf{\text{x}}})}^{T}{P}_{t}^{-1}(\mathbf{\text{x}}-\widehat{\mathbf{\text{x}}})\phantom{\rule{0.4em}{0ex}}\text{subject to}\phantom{\rule{0.4em}{0ex}}A\widehat{\mathbf{\text{x}}}\le \mathbf{\text{b}}.$$

(3)

This projection procedure is applied within unscented Kalman filter procedure to correct at every place where we move in the state-space: after spreading the sigma points **X*** _{t}*, after transforming the sigma points

In this study, for voxels that can be modeled with only one tensor, we found it preferable to have both the tensor components similarly oriented. Upon encountering a region of dispersion, the second component is poised and ready to begin branching instead of having zero weight and arbitrary orientation. To favor such solutions, we require the weights of each of the components to be not just non-negative but also greater than 0.2, and so, in our current implementation, *D* and **b** are constructed to encode the following state constraints:

$${{\lambda}_{1}}_{1},{{\lambda}_{2}}_{1},{{\lambda}_{1}}_{2},{{\lambda}_{2}}_{2}>0\phantom{\rule{1em}{0ex}}{w}_{1},{w}_{2}\ge 0.2\phantom{\rule{1em}{0ex}}{w}_{1}+{w}_{2}=1.$$

(4)

We first use experiments with synthetic data to validate our technique against ground truth. We confirm that our approach accurately recognizes crossing fibers over a broad range of angles and consistently estimates the partial volumes (Section 3.1). We then examine a real dataset to demonstrate how causal estimation is able to pick up fibers and branchings known to exist *in vivo* yet absent using other techniques (Section 3.2).

In these experiments, we compare against two alternative techniques. First, we use sharpened spherical harmonics with peak detection as described in [8] (order *l* = 8, regularization *L* = 0.006). This provides a comparison with an independently estimated nonparametric representation. Second, when performing tractography on real data, we also compare against single-tensor streamline tractography for a baseline.

Following the experimental method of generating multi-compartment synthetic data found in [2,8,15], we averaged the eigenvalues of the 300 voxels with highest fractional anisotropy (FA) in our real data set: {1200, 100, 100}*μ*m^{2}/msec. We used these eigenvalues to generate synthetic MR signals according to Eq. 1 at *b* = 1000 with 81 gradient directions on the hemisphere and introduced Rician noise (SNR ≈ 5 dB).

While the independent optimization techniques can be run on individually generated voxels, care must be taken in constructing reasonable scenarios to test the causal filter. For this purpose, we constructed a set of two-dimensional fields through which to navigate. In the middle is one long pathway where the filter starts at one end estimating a single tensor but then runs into voxels with two crossed fibers at a fixed angle and weighting. In this crossing region we calculated error statistics to compare against sharpened spherical harmonics.

From these synthetic sets, we examined detection rate, angular resolution, and estimated volume fractions and we plot the results in Fig. 2. Each column looks at a different primary-secondary weighting combination, and each row looks at a different metric. In the top row, we count how many times each technique distinguishes two separate fibers. The filtered approach *(black)* is able to detect two distinct fibers at crossing angles far below that using spherical harmonics *(red)*. Further, the filtered approach maintains such relatively high detection rates even at 80/20 partial voluming *(far right column)*. In the middle row, we look at where each technique reported two fibers and we record the error in estimated angles. From this, we see that spherical harmonics result in an angular error of roughly 15° at best and fails to detect a second component at angles below 60°. In contrast, the filtered approach has an error between 5-10° and is able to accurately estimate down to crossing angles of 30°. In the bottom row, we look at the primary fiber weight estimated by the filter. As expected, this estimate is most accurate closer to 90° *(blue line indicates true weight)*.

This study focuses on fibers originating in the corpus callosum. Specifically, we sought to trace out the lateral transcallosal fibers that run through the corpus callosum out to the lateral gyri. It is known that single-tensor streamline tractography only traces out the dominant pathways forming the U-shaped callosal radiation while spherical harmonics only capture some of these pathways [8,15].

We begin by seeding each algorithm up to thirty times in voxels at the intersection of the mid-sagital plane and the corpus callosum. To explore branchings found using the proposed technique, we considered a component to be branching if it was separated from the primary component by less than 40° with FA≥0.15 and weight above 0.3. Similarly, with sharpened spherical harmonics, we considered it a branch if we found additional maxima over the same range. We terminated fibers when either the generalized fractional anisotropy [2] of the estimated signal fell below 0.1 or the primary component FA fell below 0.15 or weight below 0.3.

We tested our approach on a human brain scan using a 3-Tesla magnet to collect 51 diffusion weighted images on the hemisphere at *b* = 900 s/mm^{2}, a scan sequence comparable those of [8,15]. Fig. 3 shows tracts originating from within a few voxels intersecting a chosen coronal slice. Confirming the results in [8,15], sharpened spherical harmonics only pick up a few fibers intersecting the U-shaped callosal radiata. In contrast, our proposed algorithm traces out many pathways consistent with the apparent anatomy. Fig. 4 shows a view of the whole corpus callosum from above. The filtered approach is able to pick up many transcallosal fibers throughout the corpus callosum as well as infiltrating the frontal gyri to a greater degree than either alternate technique. To emphasize transcallosal tracts, we color as blue those fibers exiting a corridor of ±22 mm around the mid-sagittal plane.

Filtered tractography picks up many fiber paths consistent with the underlying structures. Both single-tensor streamline and sharpened spherical harmonics are unable to find the majority of these pathways. Seed region indicated in yellow.

In this work, we demonstrated that using the unscented Kalman filter provides robust estimates of the fiber model with much higher accuracy than independent estimation techniques. Specifically, the proposed approach gives significantly lower angular error (5-10°) in regions with fiber crossings than using sharpened spherical harmonics (15-20°), and it reliably estimates the partial volume fractions.

1. Alexander A, Hasan K, Tsuruda J, Parker D. Analysis of partial volume effects in diffusion-tensor MRI. Magnetic Resonance in Medicine. 2001;45:770–780. [PubMed]

2. Tuch D, Reese T, Wiegell M, Makris N, Belliveau J, Wedeen V. High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity. Magnetic Resonance in Medicine. 2002;48:577–582. [PubMed]

3. Kreher B, Schneider J, Mader I, Martin E, Hennig J, Il'yasov K. Multitensor approach for analysis and tracking of complex fiber configurations. Magnetic Resonance in Medicine. 2005;54:1216–1225. [PubMed]

4. Peled S, Friman O, Jolesz F, Westin CF. Geometrically constrained two-tensor model for crossing tracts in DWI. Magnetic Resonance in Medicine. 2006;24(9):1263–1270. [PMC free article] [PubMed]

5. Basser P, Pajevic S. Spectral decomposition of a 4^{th}-order covariance tensor: Applications to diffusion tensor MRI. Signal Processing. 2007;87:220–236.

6. Kaden E, Knøsche T, Anwander A. Parametric spherical deconvolution: Inferring anatomical connectivity using diffusion MR imaging. NeuroImage. 2007;37:474–488. [PubMed]

7. Anderson A. Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magnetic Resonance in Medicine. 2005;54(5):1194–1206. [PubMed]

8. Descoteaux M, Deriche R, Anwander A. Technical Report 6273. INRIA; 2007. Deterministic and probabilistic Q-ball tractography: from diffusion to sharp fiber distributions.

9. Jian B, Vemuri B. A unified computational framework for deconvolution to reconstruct multiple fibers from diffusion weighted MRI. Trans on Medical Imaging. 2007;26(11):1464–1471. [PMC free article] [PubMed]

10. Jansons K, Alexander D. Persistent angular structure: New insights from diffusion MRI data. Inverse Problems. 2003;19:1031–1046. [PubMed]

11. Tournier JD, Calamante F, Gadian D, Connelly A. Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution. NeuroImage. 2004;23:1176–1185. [PubMed]

12. Özarslan E, Shepherd T, Vemuri B, Blackband S, Mareci T. Resolution of complex tissue microarchitecture using the diffusion orientation transform. NeuroImage. 2006;31(3) [PubMed]

13. Guo W, Zeng Q, Chen Y, Liu Y. Using multiple tensor deflection to reconstruct white matter fiber traces with branching. Int Symp on Biomedical Imaging. 2006:69–72.

14. Zhan W, Yang Y. How accurately can the diffusion profiles indicate multiple fiber orientations? A study on general fiber crossings in diffusion MRI. J of Magnetic Resonance. 2006;183:193–202. [PubMed]

15. Schultz T, Seidel H. Estimating crossing fibers: A tensor decomposition approach. Trans on Visualization and Computer Graphics. 2008;14(6):1635–1642. [PubMed]

16. Gössl C, Fahrmeir L, Pütz B, Auer L, Auer D. Fiber tracking from DTI using linear state space models: Detectability of the pyramidal tract. NeuroImage. 2002;16:378–388. [PubMed]

17. Björnemo M, Brun A, Kikinis R, Westin CF. Regularized stochastic white matter tractography using diffusion tensor MRI. In: Dohi T, Kikinis R, editors. MICCAI 2002. Vol. 2488. Springer; Heidelberg: 2002. pp. 435–442. LNCS.

18. Zhang F, Goodlett C, Hancock E, Gerig G. Probabilistic fiber tracking using particle filtering. Medical Image Computing and Computer Assisted Intervention (MICCAI) 2007:144–152. [PubMed]

19. Zhukov L, Barr A. Oriented tensor reconstruction: Tracing neural pathways from diffusion tensor MRI. Visualization. 2002:387–394.

20. Poupon C, Roche A, Dubois J, Mangin JF, Poupon F. Real-time MR diffusion tensor and Q-ball imaging using Kalman filtering. Medical Image Analysis. 2008;12(5):527–534. [PubMed]

21. Behrens T, Johansen-Berg H, Jbabdi S, Rushworth M, Woolrich M. Probabilistic diffusion tractography with multiple fibre orientations: What can we gain? NeuroImage. 2007;34:144–155. [PubMed]

22. Julier S, Uhlmann J. Unscented filtering and nonlinear estimation. IEEE. 2004;92(3):401–422.

23. Simon D, Simon D. Kalman filtering with inequality constraints for turbofan engine health estimation. IEE Proc–Control Theory and Appl. 2006;153(3):371–378.