# A full bi-tensor neural tractography algorithm using the unscented Kalman filter

## Abstract

We describe a technique that uses tractography to visualize neural pathways in human brains by extending an existing framework that uses overlapping Gaussian tensors to model the signal. At each point on the fiber, an unscented Kalman filter is used to find the most consistent direction as a mixture of previous estimates and of the local model. In our previous framework, the diffusion ellipsoid had a cylindrical shape, i.e., the diffusion tensor's second and third eigenvalues were identical. In this paper, we extend the tensor representation so that the diffusion tensor is represented by an arbitrary ellipsoid. Experiments on synthetic data show a reduction in the angular error at fiber crossings and branchings. Tests on in vivo data demonstrate the ability to trace fibers in areas containing crossings or branchings, and the tests also confirm the superiority of using a full tensor representation over the simplified model.

## 1 Introduction

Diffusion-weighted magnetic resonance imaging has provided the opportunity for non-invasive investigation of neural architecture of the brain. Neuroscientists use this imaging technique to find out how neurons originating from one region in the brain connect to other regions and how well-defined those connections are. The quality of the results of such studies relies heavily on the chosen fiber representation and the reconstruction method, to trace neural pathways.

For studying the microstructure of fibers, we need a model to interpret the diffusion-weighted signal. There are two main categories for such models: parametric and non-parametric models. The simplest parametric model is the diffusion tensor describing a Gaussian estimate of the diffusion orientation and its strength in each voxel. Despite its robustness, this model is inadequate in cases where several fiber populations cross, join or split in one voxel [1, 2]. Several parametric models have been introduced for handling more complex diffusion patters such as mixtures of tensors , higher-order tensors , diffusion orientation transforms  and directional functions . When fitting the data, one must make certain assumptions about the model. For example, the number of components present in a particular voxel must be determined . As seen in this paper, incorporating information from neighboring voxels helps in this process .

Usually non-parametric models contain more information about the diffusion pattern. Instead of estimating a discrete number of fibers as in parametric models, non-parametric techniques estimate the orientation distribution function (ODF) describing arbitrary fiber configurations. Q-ball imaging  was invented to compute the ODF estimation numerically using the Funk-Radon transform, and later on, the spherical harmonics would simplify computations by providing an analytic form . An algorithm for online direct estimation of single-tensor and harmonic coefficients using a linear Kalman filter has recently been introduced by . Assuming a model for the signal response of a single fiber and using spherical deconvolution provides another approach for obtaining an ODF [2022, 11, 23]. Good reviews of both parametric and non-parametric models are found in [24, 25].

Different techniques try to reconstruct neural pathways using the mentioned models. For example, deterministic tractography methods directly follow the diffusion pathways. In the case of a single-tensor model, one continuously follows the principal diffusion direction . Most multi-tensor techniques, on the other hand, try to find the number of fibers present in a voxel and to detect branching pathways [27, 6, 28]. Particle and Kalman filters have been used for path regularization in single-tensor streamline tractography . A further approach for regularizing single-tensor tractography uses a moving least squares estimate which is weighted with the previous tensor . Probabilistic tractography is an alternative to deterministic path tracing methods. Typically these methods form probabilistic notions of the connectivity based on sampling individual paths; it is an approach that quantifies the uncertainty of connections [33, 34, 13, 35]. Based on these local models, a third class of algorithms has started to develop which attempts to optimize the fiber path instead of just the local model .

While parametric methods directly describe the principal diffusion directions, interpreting the ODFs from model independent representations usually involves finding the number and orientation of principal diffusion directions present . For example, one approach is to deconvolve with a sharpening kernel before extracting maxima , while another approach decomposes a high-order tensor into several rank-1 tensors .

Almost all of the listed approaches try to fit the model at each voxel independent of other voxels, but tractography is a causal process, the next position on the fiber is always based on the diffusion found at the previous position. Based on the unscented Kalman filter, Malcolm et al. presented a filtering strategy that treats model estimation and tractography as a causal process [43, 44]. As the signal is examined at every new position, the filter recursively updates the underlying local model parameters. Thereby the filter provides the variance of that estimate and the most consistent direction in which the tractography is to be continued. Using causal estimation in this way yields inherent path regularization and accurate fiber resolution at crossing angles which is not found in independent optimization approaches. The work proposed in this paper uses this filtering strategy as foundation.

### 1.1 Our contributions

We take the approach of Malcolm et al. but change the underlying fiber model. The framework based on the unscented Kalman filter is generic and can be applied to arbitrary fiber models with a finite dimensional parameter space. The unscented Kalman filter is a good choice for the filter since the signal reconstruction is non-linear. Tractography works with the filter first finding the estimates of the model parameters and then propagating in the most consistent direction.

The existing method of  models the diffusion signal with a subset of Gaussian tensors where the second and third eigenvalues are identical, i.e., the diffusion ellipsoids are cylindrical. In this work, we extend this framework such that diffusion can be modeled by an arbitrary ellipsoid (no assumption on the eigenvalues).

By transforming these parameters, it is possible to describe the diffusion with three Euler angles and three eigenvalues. In our experiments, we compare both methods; in particular, we test the behavior in regions where crossing fiber populations are present.

Our implementation can model the diffusion signal with one, two or three Gaussian tensors. The method also enables tractography directly from the raw signal data without separate preprocessing or regularization steps.

## 2 Approach

Our approach traces neural fibers by using estimations from previous positions to guide the estimation at the current position. An unscented Kalman Filter is used in a loop that estimates the model at the current position, moves one step in the most consistent direction, and then starts estimating again. Recursive estimation improves the accuracy of resolving individual orientations and yields inherently smooth tracts despite the presence of noise and uncertainty. Each iteration begins with a near-optimal solution provided by the previous position's estimation, and therefore, the convergence of model fitting is improved and many local minima are naturally avoided.

First we explain in Section 2.1 how we model the signal with a mixture of tensors and how our fiber model works. Then we show in Section 2.3 how this model can be estimated using the unscented Kalman filter framework introduced in [43, 44].

### 2.1 Modeling local fiber orientations

In diffusion-weighted imaging, the contrast is connected to the strength of water diffusion. Our goal is to relate these signals to the underlying fiber model.

The water diffusion is measured in each voxel along a set of gradients, ${u}_{1},\dots ,{u}_{n}\in {S}^{2}$ (on the unit sphere). For each gradient, we record the corresponding signal, s = [s1, ..., s n ]T n. For voxels containing several fiber populations, the diffusion pattern can generically be described by a mixture of several weighted Gaussian tensors. The signal values can then be written as:

${s}_{i}={s}_{0}\sum _{j}{w}_{j}{\mathsf{\text{e}}}^{-b{u}_{i}^{\mathsf{\text{T}}}{D}_{j}{u}_{i}}$
(1)

where s0 is the baseline signal which is obtained without applying any gradients, b is an acquisition-specific constant known as b-value (which is a measurement of the strength of the diffusion weighting), w j are convex weights, and D j is a tensor matrix characterizing the diffusion pattern.

Given the generic mixture model (1), we choose to work with two tensors. Our implementation also supports one or three tensors but several studies have shown that a two-fiber model performs best in our environment with b = 1000 s mm-2[4, 6, 28, 39, 7, 13]. Second, we decided to weigh both tensors equally as suggested in the study of . At a first glance, this might seem to limit flexibility, but the unscented Kalman filter adjusts the eigenvalues to better fit the signal. This has almost the same effect as scaling the tensors .

With those assumptions we end up with the following formulation of the fiber model:

${s}_{i}={s}_{0}\left(\frac{1}{2}{\mathsf{\text{e}}}^{-b{u}_{i}^{\mathsf{\text{T}}}{D}_{1}{u}_{i}}+\frac{1}{2}{\mathsf{\text{e}}}^{-b{u}_{i}^{\mathsf{\text{T}}}{D}_{2}{u}_{i}}\right)$
(2)

where D1, D2 can each be written as D = λ1mmT + λ2ppT + λ3qqT, with $m,p,q\in {S}^{2}$ forming an orthonormal basis aligned with the axes of the diffusion ellipsoid. Our implementation restricts each λ to be positive, and without loss of generality, we assume the eigenvalues in decreasing order, i.e., λ1λ2λ3. m is the principal diffusion direction since it corresponds to largest eigenvalue λ1. The model parameters are m 1, p 1, q 1, λ11, λ21, λ31, m2, p2, q2, λ12, λ22 and λ32. Malcolm et al. use a simplified model in [43, 44] where the second and third eigenvalues are equal, which yields the following formulation: D = λ1mmT + λ2 (ppT + qqT). Our model overcomes this simplification and allows for a better fit to the signal. At first, it appears that the large number of free parameters in our new model might be a drawback (the simpler model only has the parameters m1, λ11, λ21, m2, λ12 and λ22). But as we are about to show in Section 2.2 it is possible to represent the orientation of the diffusion tensor with Euler angles. Then we end up with six parameters per tensor, three Euler angles plus three eigenvalues (the existing simpler model stores the principal diffusion direction and two eigenvalues per tensor).

The formulation for the two-tensor model can directly be extended to a three-tensor version:

${s}_{i}={s}_{0}\frac{1}{3}\sum _{j=1}^{3}{\mathsf{\text{e}}}^{-b{u}_{i}^{\mathsf{\text{T}}}{D}_{j}{u}_{i}}$
(3)

which uses the additional parameters m3, p3, q3, λ13, λ23, λ33.

### 2.2 Tensor representation with Euler angles

We use singular value decomposition to rewrite the diffusion tensor matrix D (in this special case it is identical to using the eigenvalue decomposition since D is real, symmetric and positive-definite) which yields:

$D=Q\Lambda {Q}^{\mathsf{\text{T}}}$
(4)

with Q being a rotation matrix whose columns are the orthonormal eigenvectors of D and Λ being a diagonal matrix containing the eigenvalues:

$\Lambda =\left(\begin{array}{ccc}\hfill {\lambda }_{1}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\lambda }_{2}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\lambda }_{3}\hfill \end{array}\right)$
(5)

As in , we use the ZYZ Euler angle convention that allows splitting Q up into three individual rotations around the z-axis, the y-axis and around the z-axis again. The amount of the rotations is given by the Euler angles ϕ, θ, ψ:

$Q=\left(\begin{array}{ccc}\hfill {Q}_{11}\hfill & \hfill {Q}_{12}\hfill & \hfill {Q}_{13}\hfill \\ \hfill {Q}_{21}\hfill & \hfill {Q}_{22}\hfill & \hfill {Q}_{23}\hfill \\ \hfill {Q}_{31}\hfill & \hfill {Q}_{32}\hfill & \hfill {Q}_{33}\hfill \end{array}\right)={R}_{z}\left(\varphi \right){R}_{y}\left(\theta \right){R}_{z}\left(\psi \right)$
(6)

where R y , R z are the rotation matrices around the y- and the z-axis, respectively:

${R}_{y}\left(\Omega \right)=\left(\begin{array}{ccc}\hfill cos\left(\Omega \right)\hfill & \hfill 0\hfill & \hfill sin\left(\Omega \right)\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill -sin\left(\Omega \right)\hfill & \hfill 0\hfill & \hfill cos\left(\Omega \right)\hfill \end{array}\right)\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{R}_{z}\left(\Omega \right)=\left(\begin{array}{ccc}\hfill cos\left(\Omega \right)\hfill & \hfill -sin\left(\Omega \right)\hfill & \hfill 0\hfill \\ \hfill sin\left(\Omega \right)\hfill & \hfill cos\left(\Omega \right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right)$
(7)

Fully written out, Q looks as follows:

$\left(\begin{array}{c}\hfill {Q}_{11}\hfill \\ \hfill {Q}_{21}\hfill \\ \hfill {Q}_{31}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill cos\left(\varphi \right)cos\left(\theta \right)cos\left(\psi \right)-sin\left(\varphi \right)sin\left(\psi \right)\hfill \\ \hfill sin\left(\varphi \right)cos\left(\theta \right)cos\left(\psi \right)+cos\left(\varphi \right)sin\left(\psi \right)\hfill \\ \hfill -sin\left(\theta \right)cos\left(\psi \right)\hfill \end{array}\right)$
(8)
$\left(\begin{array}{c}\hfill {Q}_{12}\hfill \\ \hfill {Q}_{22}\hfill \\ \hfill {Q}_{33}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill -cos\left(\varphi \right)cos\left(\theta \right)sin\left(\psi \right)-sin\left(\varphi \right)cos\left(\psi \right)\hfill \\ \hfill sin\left(\varphi \right)cos\left(\theta \right)sin\left(\psi \right)+cos\left(\varphi \right)cos\left(\psi \right)\hfill \\ \hfill sin\left(\theta \right)cos\left(\psi \right)\hfill \end{array}\right)$
(9)
$\left(\begin{array}{c}\hfill {Q}_{13}\hfill \\ \hfill {Q}_{23}\hfill \\ \hfill {Q}_{33}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill cos\left(\varphi \right)sin\left(\theta \right)\hfill \\ \hfill sin\left(\varphi \right)sin\left(\theta \right)\hfill \\ \hfill cos\left(\theta \right)\hfill \end{array}\right)$
(10)

Extracting the Euler angles from Q can be carried out as explained in . θ = cos-1(Q33) is trivial, and then, if θ ≠ 0 the other two angles can be obtained through ϕ = atan2(Q23, Q13) and ψ = atan2(Q32, Q31), where atan2 is a function commonly defined in programming languages. Otherwise if θ = 0, ϕ and ψ are not uniquely determined. One solution is ψ = 0, and ϕ = atan2(-Q12, Q22) .

### 2.3 Estimating the fiber model

Given the scanned signal at a particular voxel, we want to estimate the underlying model parameters that best explain this signal. Every fiber is treated as a trajectory of a particle. At each step, we examine the measured signal, estimate the underlying model parameters, and propagate forward in the most consistent direction, which is the principal diffusion direction most aligned with the incoming vector.

Following we define the components used by the unscented Kalman filter:

1. 1.

The system state x: the model parameters

2. 2.

The state transition function f[·]: how the model changes as we trace the fiber

3. 3.

The observation function h[·]: how the signal appears given a particular state

4. 4.

The measurement y: the actual signal obtained from the scanner

The state vector is given by:

$x=\left[{\varphi }_{1}\phantom{\rule{0.5em}{0ex}}{\theta }_{1}\phantom{\rule{0.5em}{0ex}}{\psi }_{1}\phantom{\rule{0.5em}{0ex}}{\lambda }_{1}1\phantom{\rule{0.5em}{0ex}}{\lambda }_{2}1\phantom{\rule{0.5em}{0ex}}{\lambda }_{3}1\phantom{\rule{0.5em}{0ex}}{\varphi }_{2}\phantom{\rule{0.5em}{0ex}}{\theta }_{2}\phantom{\rule{0.5em}{0ex}}{\psi }_{2}\phantom{\rule{0.5em}{0ex}}{\lambda }_{1}2\phantom{\rule{0.5em}{0ex}}{\lambda }_{2}2\phantom{\rule{0.5em}{0ex}}{\lambda }_{3}2{\right]}^{\text{T}}$
(11)

where ϕ, θ, ψ, and λ +. Since the local fiber configuration does not drastically change from one position to the next, we assume identity dynamics for the state transition function f [·]. Our observation h[·] is the signal reconstruction, s = [s1, ..., s n ]T which is given as the left hand side of (2). Using (6) and (4) allows to calculate the diffusion tensor matrix D from the Euler angles, which is then plugged into (2). The measurement y is the actual signal from the diffusion-weighted images. At subvoxel positions, we interpolate directly on the signal.

The signal reconstruction is a non-linear process, and therefore, we use an unscented Kalman filter. A Kalman filter tries to reconcile the predicted state of the system with the measured state. The prediction and measurement processes can be non-linear in the unscented version, whereas they are only linear in the classic Kalman filter. The extended Kalman filter would have been an alternative. However, the extended Kalman filter is only a first-order approximation whereas the unscented version is accurate up to the second-order moment of the state distribution. See [46, 47] for more details about the unscented Kalman filter, particularly  shows the superiority of the unscented Kalman filter over the extended Kalman filter.

Particle filters would be another approach for non-linear estimation, but their number of particles used is exponential to the state dimension n. In contrast, the unscented Kalman filter requires a linear number of particles (2n + 1 sigma points) and is therefore computationally less complex [43, 44] (especially in the three-tensor case). Further, the distribution of the state is not likely to be highly complex (which would warrant the use of a particle filter), and therefore, it can be sufficiently captured by the unscented Kalman filter.

The system of interest is at time t, and we have a Gaussian estimate of its current state with mean, x t n, and covariance, P t n × n. Prediction begins with the formation of a set sample states called sigma points, X t = {χ i } n of 2n + 1. Each of these states has an associated convex weight, w i . The covariance, P t , is used to deterministically distribute the sigma points around the current state:

$\begin{array}{c}{\chi }_{0}={x}_{t}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{w}_{0}=\kappa ∕\left(n+\kappa \right)\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}{w}_{i}={w}_{i+n}=\frac{1}{2\left(n+\kappa \right)}\\ {\chi }_{i}={x}_{t}+{\left[\sqrt{\left(n+\kappa \right){P}_{t}}\right]}_{i}\phantom{\rule{1em}{0ex}}{\chi }_{i+n}={x}_{t}-{\left[\sqrt{\left(n+\kappa \right){P}_{t}}\right]}_{i}\end{array}$
(12)

with [A] i denoting the ith column of matrix A and κ is an adjustable scaling parameter (we use κ = 0.01 in all our experiments). In a next step the sigma points are propagated through the state transition function, $\stackrel{^}{\chi }=f\left[\chi \right]\in {ℝ}^{n}$, and a new set of predicted sigma points is obtained: ${X}_{t+1|t}=\left\{f\left[{\chi }_{i}\right]\right\}=\left\{{\stackrel{^}{\chi }}_{i}\right\}$.

As mentioned, we assume that the fiber configuration does not change abruptly from one voxel to the next. This is modeled with an identity transition function: xt+1|t= f [x t ] = x t . Next, the predicted system mean state and covariance are predicted:

$\begin{array}{lll}\hfill {\stackrel{̄}{x}}_{t+1|t}& =\sum _{i}{w}_{i}{\stackrel{^}{\chi }}_{i}\phantom{\rule{2em}{0ex}}& \hfill \text{(1)}\\ \hfill {P}_{xx}& =\sum _{i}{w}_{i}\left({\stackrel{^}{\chi }}_{i}-{\stackrel{̄}{x}}_{t+1|t}\right){\left({\stackrel{^}{\chi }}_{i}-{\stackrel{̄}{x}}_{t+1|t}\right)}^{\mathsf{\text{T}}}+Q\phantom{\rule{2em}{0ex}}& \hfill \text{(2)}\\ \hfill \text{(3)}\end{array}$
(13)

where Q is the Kalman filter's process noise which ensures a non-null spread of sigma points and a positive-definite covariance. The explained sampling technique is known as the unscented transform, and it is used to estimate the behavior of a non-linear function. The sigma points are spread based on the current uncertainty, propagated and then their spread is measured.

The predicted observation is obtained by applying the unscented transform again, this time using the predicted states, Xt+1|t, to estimate what we expect to observe from the hypothetical measurement of each state: $\gamma =h\left[\stackrel{^}{\chi }\right]\in {ℝ}^{m}$. Since the observation is the reconstructed signal, this step estimates the diffusion-weighted signal. Then, we calculate the predicted set of observations, ${Y}_{t+1|t}=\left\{h\left[{\stackrel{^}{\chi }}_{i}\right]\right\}=\left\{{\gamma }_{i}\right\}$, and we obtain their mean and covariance,

$\begin{array}{lll}\hfill {\stackrel{̄}{y}}_{t+1|t}& =\sum _{i}{w}_{i}{\stackrel{^}{\gamma }}_{i}\phantom{\rule{2em}{0ex}}& \hfill \text{(1)}\\ \hfill {P}_{yy}& =\sum _{i}{w}_{i}\left({\stackrel{^}{\gamma }}_{i}-{\stackrel{̄}{y}}_{t+1|t}\right){\left({\stackrel{^}{\gamma }}_{i}-{\stackrel{̄}{y}}_{t+1|t}\right)}^{\mathsf{\text{T}}}+R\phantom{\rule{2em}{0ex}}& \hfill \text{(2)}\\ \hfill \text{(3)}\end{array}$
(14)

where R is the observation noise of the Kalman filter. The cross-correlation between the estimated state and observation are calculated as:

${P}_{xy}=\sum _{i}{w}_{i}\left({\stackrel{^}{\chi }}_{i}-{\stackrel{̄}{x}}_{t+1|t}\right){\left({\stackrel{^}{\gamma }}_{i}-{\stackrel{̄}{y}}_{t+1|t}\right)}^{\mathsf{\text{T}}}$
(15)

As in the classic linear Kalman filter, the final step uses the Kalman gain, $K={P}_{xy}{P}_{yy}^{-1}$, to correct the prediction and provide the final estimated system mean and covariance,

${x}_{t+1}={\stackrel{̄}{x}}_{t+1|t}+K\left({y}_{t}-{\stackrel{̄}{y}}_{t+1|t}\right)$
(16)
${P}_{t+1}={P}_{xx}-K{P}_{yy}{K}^{\mathsf{\text{T}}}$
(17)

where y t m is the actual signal measurement taken at this time.

In a few cases erroneous estimation of the reconstruction signal may result in fibers which do not exist in reality. A solution for detecting and removing such false positives was introduced in .

### 2.4 The algorithm

In summary we use the unscented Kalman filter to estimate the local model parameters as we trace each fiber. For each fiber, we store the position at which we are currently tracing it and the current estimate of its model parameters (mean and covariance). Each iteration of the algorithm predicts the new state (which is identity in our case): xt+1|t= x t . The observation y t in (16) is the diffusion-weighted signal coming from the scanner, s. With these, we use the above equations to find the new estimated model parameters, xt+1. Last, we use first-order forward Euler integration to move a small step in the most consistent principal diffusion direction, m, which is given as the first column of the rotation matrix Q, see (6). Afterward, the same is repeated for the new position. Fractional anisotropy (FA)  is a measurment of a diffusion tensor's anisotropy, and it is used as break condition for the tractography. A fiber is followed until the FA falls below a certain level (in our experiments meaning full FA thresholds lay in [0.1, 0.15]). The full procedure is presented in Algorithm (1).

Algorithm 1 Main tractography loop repeated for each fiber

1. 1.

Initialize x0 and P0

2. 2.

repeat

3. 3.

Form the sigma points X t around x t

4. 4.

Predict the new sigma points Xt+1|tand observations Yt+1|t

5. 5.

Compute weighted means and covariances, e.g. ${\stackrel{̄}{x}}_{t+1|t},{P}_{xy}$

6. 6.

Update estimate (xt+1, Pt+1) using scanner measurement y t

7. 7.

Obtain the tensors' principal diffusion directions from the state by putting the Euler angles into (8)

8. 8.

Use the principal diffusion direction m j most aligned to the incoming vector to proceed

9. 9.

until estimated model appears isotropic

For initializing x0 before the main loop of the algorithm is started, we use a one-tensor model which is easily derived from (1):

${s}_{i}={s}_{0}{\mathsf{\text{e}}}^{-b{u}_{i}^{\mathsf{\text{T}}}D{u}_{i}}$
(18)

This is an overdetermined system (assuming s0 is known and that more than 6 gradients are used, what commonly is the case), and the diffusion tensor matrix D can be obtained by using a least squares method. No matter how many tensors are used, the initial state of each tensor is the same. P0 is initialized as a diagonal matrix with small values (0.01I in our experiments, with I being the identity matrix).

## 3 Experiments

First, we use synthetically generated diffusion-weighted images to validate our technique against ground truth. We perform tractography through crossing fiber fields of different angles and examine the resulting orientations and branchings. Our new full two-tensor algorithm is compared with the existing simpler model of Malcolm et al. [43, 44] on the same synthetic crossings. For the comparison we measure several statistics to show that the full model performs better. Lastly in Section 3.3, we examine a real data set to demonstrate that our full model works on in vivo data too and that it finds fibers that do not appear with the simpler model (but are known to exist anatomically).

Note the importance of the matrices that inject process noise Q and observation noise R (see (13) and (14)). Their magnitude can be determined by the data itself by running our synthetic experiments with different values and choosing those that yield the best results. Usually the noise matrices are diagonal matrices. The process noise determines how much variance is allowed in the model: high values allow for more variation but if it gets too high the estimation could become inaccurate. We found that the best results are obtained with ,θ, ψ, [0.001, 0.002] and with qλ 1,2,3≈ 100 (those values are placed on Q's diagonal in the same order as the variables appear in the state x). These values allow an appropriate amount of angular and diffusive flexibility. The injected observation noise governs how much variance is expected in the measurement: higher values mean we expect more variance and hence trust our measurement less. For our experiments, we found r s [0.01, 0.03] to work quite well (R is also a diagonal matrix and all the values on the diagonal are r s ) but experimentation might be necessary since this value depends on the amount of physical noise present, which again varies depending on the scanner, protocol or pre-processing.

### 3.1 Synthetic fiber crossings

We generated synthetic MR signals according to (2) using {1700, 500, 300}μ m2 msec-1 (FA = 0.73) as the tensor's eigenvalues to form an anisotropic tensor at b = 1000 s mm-2. The used eigenvalues are typical for brain white matter. We use 81 gradient directions uniformly spread on the hemisphere. We assume s0 = 1, and we perform each experiment with two different levels of Rician noise (SNR ≈ 5 dB for low noise and SNR ≈ 20dB for high). We also ran all the experiments with b = 3000 s mm -2.

Since the algorithm depends on previous estimations, it is not enough to inspect individual voxels for testing. We construct a 2D field with crossings which the fiber has to navigate through. Figure 1a depicts a schematic of one such field with a 60° crossing. The fibers start on the left in a region with only one true fiber population present, and they try to find their way to the right side. In the middle the fibers encounter a region with crossing fibers at a fixed angle. Figure 1b shows the resulting fibers. In blue we present the simpler model and in red our new full model. The simpler model drifts off shortly after entering the crossing region while our full model maintains a straight path for longer. Outside the crossing region, where only one fiber exists, the second component is aligned with the first. The filter begins estimating a single-tensor until it hits the crossing region where the two components start to point in different directions. After the crossing both components will realign again. We generated a set of similar fields for crossings angles of [0°, 5°, ..., 90°].

In Figure 1c we take a closer look at several points along the fibers as they enter the crossing region (the close-up area is marked with a yellow frame in Figure 1b). At every fiber point the principal diffusion directions of both tensors are shown. Again the simpler model is in blue and the full one in red. The principal diffusion directions adapt stepwise, until they are aligned 60° angle in the crossing region. Our full model adapts quicker than the simple model.

### 3.2 Measurements

Having verified the underlying behavior, we then began a more comprehensive evaluation. We measured the following three different statistics: the average angular error in the non-crossing regions, the average angular error in the crossing region, and the average absolute FA error. With angular error we mean the difference in angle between the principal diffusion directions of the tensors. The estimated values of each fiber point in a certain region (crossing or non-crossing) are compared to the ground truth, the difference is calculated, and the average is taken. The FA error is calculated over both regions at once.

Figure 2 shows the angular error in the non-crossing regions, Figure 3 depicts the angular error in the crossing region and Figure 4 displays the absolute FA error. Each graph plots the crossing angle, from 0° to 90°, versus the error. As mentioned we ran every experiment for two different noise levels and for b = 1000 s mm-2 and b = 3000 s mm-2. We seeded from 18 different voxels in each test scenario, and also, our crossing region is quite large as you can see in Figure 1. This guarantees a fair number of used tensors for averaging. In each graph the trendlines indicate the mean error while the bars indicate the standard deviation.

The graphs for the angular errors show that our new full model generally performs better. In the non-crossing region our model is more stable, especially for the cases with higher crossing angles of 60°-90°. For these angles the simpler model's error raises quickly while our full model is almost constant. In the crossing region the full model performs better for lower angles 0°-45°. Working with different b-values or with a different noise levels does not significantly change the errors. The absolute FA error is lower and yields a lower variance with our new full tensor algorithm (see Figure 4). This is the expected result of having independent second and third eigenvalues. In the simpler model the second and third eigenvalues have to be identical, which makes it more difficult to match the signal.

### 3.3 In vivo tractography

We tested our approach also on a real human brain scan which was weighted with 51 gradients directions. The diffusion-weighted images have a voxel size of 1.66 × 1.66 × 1.7 mm2 and b = 900 s mm-2. Figure 5 shows the results from a tractography started in the right Thalamus of the human brain. The blue fibers were obtained with the simpler fiber model whereas the red fibers were traced with our new full model. At a first glance both methods seem similar but the separated close-ups in Figure 5b, c show that our approach finds fibers in certain areas where the simpler model stops.

## 4 Conclusions

We used the tractography framework of Malcolm et al. [43, 44] which is built on an unscented Kalman filter and extended the used fiber model. By changing the fiber model and allowing it to have three different, independent eigenvalues, the tractography procedure was improved and the angular and FA errors were minimized. By representing the orientation of the diffusion ellipsoid with Euler angles, the Kalman filter's state only contains one variable more per tensor than the simpler model. The small increase in calculation time due to this additional variable is negligible. We confirm that using a causal filter for tractography performs much better than independent alternatives. We believe that exploring more alternative fiber models and plugging them into the filtering framework will provide new insights into neural pathways and that they, ultimately, will enhance non-invasive diagnosis of human brains. Furthermore, exploring filtering techniques other than the unscented Kalman filter might yield promising results.

## References

1. Alexander DC, Barker G, Arridge S: Detection and modeling of non-Gaussian apparent diffusion coefficient profiles in human brain data. Magn Reson Med 2002, 48: 331-340. 10.1002/mrm.10209

2. Frank L: Characterization of anisotropy in high angular resolution diffusion-weighted MRI. Magn Reson Med 2002, 47: 1083-1099. 10.1002/mrm.10156

3. Alexander A, Hasan K, Tsuruda J, Parker D: Analysis of partial volume effects in diffusion-tensor MRI. Magn Reson Med 2001, 45: 770-780. 10.1002/mrm.1105

4. Tuch D, Reese T, Wiegell M, Makris N, Belliveau J, Wedeen V: High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity. Magn Reson Med 2002, 48: 577-582. 10.1002/mrm.10268

5. Parker G, Alexander DC: Probabilistic anatomical connectivity derived from the microscopic persistent angular structure of cerebral tissue. Phil Trans R Soc B 2005, 360: 893-902. 10.1098/rstb.2005.1639

6. Kreher B, Schneider J, Mader I, Martin E, Hennig J, Il'yasov K: Multi-tensor approach for analysis and tracking of complex fiber configurations. Magn Reson Med 2005, 54: 1216-1225. 10.1002/mrm.20670

7. Peled S, Friman O, Jolesz F, Westin CF: Geometrically constrained two-tensor model for crossing tracts in DWI. Magn Reson Med 2006,24(9):1263-1270.

8. Hlawitschka M, Scheuermann G: HOT-lines: tracking lines in higher order tensor fields. Visualization 2005, 27-34.

9. Özarslan E, Shepherd T, Vemuri B, Blackband S, Mareci T: Resolution of complex tissue microarchitecture using the diffusion orientation transform. NeuroImage 2006,31(3):1086-1103. 10.1016/j.neuroimage.2006.01.024

10. McGraw T, Vemuri B, Yezierski B, Mareci T: Von Mises-Fisher mixture model of the diffusion ODF. International Symposium on Biomedical Imaging 2006, 65-68.

11. Kaden E, Knøsche T, Anwander A: Parametric spherical deconvolution: inferring anatomical connectivity using diffusion MR imaging. NeuroImage 2007, 37: 474-488. 10.1016/j.neuroimage.2007.05.012

12. Rathi Y, Michailovich O, Shenton ME, Bouix S: Directional functions for orientation distribution estimation. Med Image Anal 2009, 13: 432-444. 10.1016/j.media.2009.01.004

13. 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. 10.1016/j.neuroimage.2006.09.018

14. King MD, Gadian DG, Clark CA: A random effects modelling approach to the crossing-fibre problem in tractography. NeuroImage 2009, 44: 753-768. 10.1016/j.neuroimage.2008.09.058

15. Tuch D: Q-ball imaging. Magn Reson Med 2004, 52: 1358-1372. 10.1002/mrm.20279

16. Anderson A: Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magn Reson Med 2005,54(5):1194-1206. 10.1002/mrm.20667

17. Hess C, Mukherjee P, Han E, Xu D, Vigneron D: Q-ball reconstruction of multimodal fiber orientations using the spherical harmonic basis. Magn Reson Med 2006, 56: 104-117. 10.1002/mrm.20931

18. Descoteaux M, Angelino E, Fitzgibbons S, Deriche R: Regularized, fast, and robust analytical Q-ball imaging. Magn Reson Med 2007, 58: 497-510. 10.1002/mrm.21277

19. Poupon C, Roche A, Dubois J, Mangin JF, Poupon F: Real-time MR diffusion tensor and Q-ball imaging using Kalman filtering. Med Image Anal 2008,12(5):527-534. 10.1016/j.media.2008.06.004

20. Jian B, Vemuri B: A unified computational framework for deconvolution to reconstruct multiple fibers from diffusion weighted MRI. Trans Med Imag 2007,26(11):1464-1471.

21. Jansons K, Alexander DC: Persistent angular structure: new insights from diffusion MRI data. Inverse Probl 2003, 19: 1031-1046. 10.1088/0266-5611/19/5/303

22. Tournier JD, Calamante F, Gadian D: A Connelly, Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution. NeuroImage 2004, 23: 1176-1185. 10.1016/j.neuroimage.2004.07.037

23. Kumar R, Barmpoutis A, Vemuri BC, Carney PR, Mareci TH: Multi-fiber reconstruction from DW-MRI using a continuous mixture of von Mises-Fisher distributions. Mathematical Methods in Biomedical Image Analysis (MMBIA) 2008, 1-8.

24. Alexander DC: Multiple-fiber reconstruction algorithms for diffusion MRI. Annal NY Acad Sci 2005,1046(1):113-133.

25. Descoteaux M, Deriche R, Knoesche T, Anwander A: Deterministic and probabilistic tractography based on complex fiber orientation distributions. Trans Med Imag 2009,28(2):269-286.

26. Basser PJ, Pajevic S, Pierpaoli C, Duda J: A Aldroubi: in vivo fiber tractography using DT-MRI data. Magn Reson Med 2000, 44: 625-632. 10.1002/1522-2594(200010)44:4<625::AID-MRM17>3.0.CO;2-O

27. Hagmann P, Reese T, Tseng WY, Meuli R, Thiran JP, Wedeen VJ: Diffusion spectrum imaging tractography in complex cerebral white matter: an investigation of the centrum semiovale. International Symposium on Magnetic Resonance in Medicine (ISMRM) 2004, 623.

28. Guo W, Zeng Q, Chen Y, Liu Y: Using multiple tensor deflection to reconstruct white matter fiber traces with branching. International Symposium on Biomedical Imaging 2006, 69-72.

29. Gössl C, Fahrmeir L, utz BP, Auer L, Auer D: Fiber tracking from DTI using linear state space models: detectability of the pyramidal tract. NeuroImage 2002, 16: 378-388. 10.1006/nimg.2002.1055

30. Björnemo M, Brun A, Kikinis R, Westin CF: Regularized stochastic white matter tractography using diffusion tensor MRI. Medical Image Computing and Computer Assisted Intervention (MICCAI) 2002, 435-442.

31. Zhang F, Hancock E, Goodlett C, Gerig G: Probabilistic white matter fiber tracking using particle filtering and von Mises-Fisher sampling. Med Image Anal 2009, 13: 5-18. 10.1016/j.media.2008.05.001

32. Zhukov L, Barr A: Oriented tensor reconstruction: tracing neural pathways from diffusion tensor MRI. Visualization 2002, 387-394.

33. Parker G, Alexander DC: Probabilistic Monte Carlo based mapping of cerebral connections utilizing whole-brain crossing fiber information. Information Processing in Medical Imaging (IPMI) 2003, 684-696.

34. Hosey T, Ansorge R: Inference of multiple fiber orientations in high angular resolution diffusion imaging. Magn Reson Med 2005, 54: 1480-1489. 10.1002/mrm.20723

35. Iturria-Medina Y, Canales-Rodríguez EJ, Melie-García L, Valdés-Hernández PA, Martńez-Montes E, Alemán-Gómez Y, Sánchez-Bornot JM: Characterizing brain anatomical connections using diffusion weighted MRI and graph theory. NeuroImage 2007, 36: 645-660. 10.1016/j.neuroimage.2007.02.012

36. Fillard P, Poupon C, Mangin JF: A novel global tractography algorithm based on an adaptive spin glass model. Medical Image Computing and Computer Assisted Intervention (MICCAI) 2009, 927-934.

37. Jbabdi S, Woolrich M, Andersson J, Behrens T: A bayesian framework for global tractography. NeuroImage 2007, 37: 116-129. 10.1016/j.neuroimage.2007.04.039

38. Kreher B, Madeer I, Kiselev V: Gibbs tracking: a novel approach for the reconstruction of neuronal pathways. Magn Reson Med 2008, 60: 953-963. 10.1002/mrm.21749

39. Zhan W, Yang Y: How accurately can the diffusion profiles indicate multiple fiber orientations? A study on general fiber crossings in diffusion MRI. J Magn Reson 2006, 183: 193-202. 10.1016/j.jmr.2006.08.005

40. Seunarine K, Cook P, Hall M, Embleton K, Parker G, Alexander DC: Exploiting peak anisotropy for tracking through complex structures. Mathematical Methods in Biomedical Image Analysis (MMBIA) 2007, 1-8.

41. Jian B, Vemuri B, Özarslan E, Carney PR, Mareci TH: A novel tensor distribution model for the diffusion -weighted MR signal. NeuroImage 2007,37(1):164-176. 10.1016/j.neuroimage.2007.03.074

42. Schultz T, Seidel H: Estimating crossing fibers: a tensor decomposition approach. Trans Vis Comput Graph 2008,14(6):1635-1642.

43. Malcolm JG, Shenton ME, Rathi Y: Neural tractography using an unscented kalman filter. Inf Process Med Imaging (IPMI) 2009, 21: 126-138.

44. Malcolm JG, Shenton ME, Rathi Y: Filtered multi-tensor tractography. IEEE Trans Med Imaging 2010, 29: 1664-1675.

45. Koay CG, Chang LC, Pierpaoli C, Basser PJ: Error propagation framework for diffusion tensor imaging via diffusion tensor representations. IEEE Trans Med Imaging 2007,26(8):1017-1034.

46. Julier S, Uhlmann J: Unscented filtering and nonlinear estimation. IEEE 2004,92(3):401-422. 10.1109/JPROC.2003.823141

47. van der Merwe R, Wan E: Sigma-point Kalman filters for probabilistic inference in dynamic state-space models. Workshop on Advances in Machine Learning 2003.

48. Rathi Y, Malcolm JG, Bouix S, Westin CF, Shenton ME: False positive detection using filtered tractography. International Symposium on Magnetic Resonance in Medicine (ISMRM) 2010.

49. Basser PJ, Pierpaoli C: Toward a quantitative assessment of diffusion anisotropy. Magn Reson Med 1996,36(6):893-906. 10.1002/mrm.1910360612

## Author information

Authors

### Corresponding author

Correspondence to Stefan Lienhard.

### 5 Competing interests

The authors declare that they have no competing interests.

## Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

## Rights and permissions

Reprints and Permissions

Lienhard, S., Malcolm, J.G., Westin, CF. et al. A full bi-tensor neural tractography algorithm using the unscented Kalman filter. EURASIP J. Adv. Signal Process. 2011, 77 (2011). https://doi.org/10.1186/1687-6180-2011-77

• Accepted:

• Published:

• DOI: https://doi.org/10.1186/1687-6180-2011-77

### Keywords

• Fractional Anisotropy
• Kalman Filter
• Euler Angle
• Extended Kalman Filter
• Orientation Distribution Function
{"ft_lang_label":"__label__en","ft_lang_prob":0.83850574,"math_prob":0.96459866,"size":36840,"snap":"2023-14-2023-23","text_gpt3_token_len":8650,"char_repetition_ratio":0.14214356,"word_repetition_ratio":0.01888112,"special_character_ratio":0.23797503,"punctuation_ratio":0.14451522,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.99372804,"pos_list":[0,1,2],"im_url_duplicate_count":[null,3,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2023-04-01T23:32:40Z\",\"WARC-Record-ID\":\"<urn:uuid:abbdd8db-0281-4cd0-88fb-20e32f960aee>\",\"Content-Length\":\"448521\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:d71a94a6-95f7-4494-a5bc-a679caa65f0e>\",\"WARC-Concurrent-To\":\"<urn:uuid:abcf76f7-5b93-40a3-af0a-7fbbbb4cdc69>\",\"WARC-IP-Address\":\"146.75.32.95\",\"WARC-Target-URI\":\"https://asp-eurasipjournals.springeropen.com/articles/10.1186/1687-6180-2011-77\",\"WARC-Payload-Digest\":\"sha1:HFEGE6LC44ASLOFYTVUIJKWHEZFAVJVL\",\"WARC-Block-Digest\":\"sha1:PVNMIBNTYSPKLDVHZ7ZCRVT5GO3GOJ4Q\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2023/CC-MAIN-2023-14/CC-MAIN-2023-14_segments_1679296950363.89_warc_CC-MAIN-20230401221921-20230402011921-00492.warc.gz\"}"}