A Random Riemannian Metric for Probabilistic
Shortest-Path Tractography
Søren Hauberg1 , Michael Schober2 , Matthew Liptrot1,3 , Philipp Hennig2 , Aasa
Feragen3 ,1

Cognitive Systems, Technical University of Denmark, Denmark
Max Planck Institute for Intelligent Systems, Tübingen, Germany
3
Department of Computer Science, University of Copenhagen, Denmark2

Abstract. Shortest-path tractography (SPT) algorithms solve global
optimization problems defined from local distance functions. As diffusion
MRI data is inherently noisy, so are the voxelwise tensors from which local distances are derived. We extend Riemannian SPT by modeling thestochasticity of the diffusion tensor as a “random Riemannian metric”,
where a geodesic is a distribution over tracts. We approximate this distribution with a Gaussian process and present a probabilistic numerics
algorithm for computing the geodesic distribution. We demonstrate SPT
improvements on data from the Human Connectome Project.

1

Introduction

Diffusion weighted imaging enables inference of structural brain connectivity,
assuming that water diffuses along brain fibers, not across. At any brain location,a local brain fiber orientation distribution function (ODF) is estimated frommeasured diffusion along a fixed set of directions. Integration of the field of local
ODFs yields a heuristic estimate of the most likely path for a brain fiber.
We study the stochasticity of the ODF estimate and its effect on the distribution of estimated tracts through Riemannian shortest-path tractography (SPT)in diffusion tensor imaging (DTI), where the ODF is a second order tensor describing a Gaussian diffusion process at any given location. In Riemannian SPT,
tracts are estimated as geodesics on a Riemannian manifold, where the Riemannian metric is constructed from the diffusion tensor [9, 14, 15]. As the diffusion
tensor is estimated from data, it is subject to noise. The tensor, and the induced
Riemannian metric, should therefore be treated as stochastic variables, not exact
representatives of the true underlying diffusion process. This is not possible incurrent Riemannian tractography models, which assume a deterministic metric.
We present a solution to probabilistic SPT for DTI which treats diffusion
tensors as stochastic variables, and returns a distribution over the tracts connecting two given points in the brain, approximated by a Gaussian process (GP).
Algorithmically, we extend recent probabilistic solutions for ordinary differential equations (ODEs) [12, 18] by allowing “noisy” evaluations of the ODE due
to the noisy metric. The resulting ODE solver is extendable to other problems.
The combined model yields both quantitative and qualitative improvements overstandard SPT algorithms.
1

2

Shortest-Path Tractography and Random Geometries

In shortest-path tractography (SPT) the ODFs induce a local Riemannian metric [9, 11, 14, 15], where the most probable fiber connecting two points x and yis re-interpreted as the geodesic, or shortest path, connecting x and y.
Why SPT? Walk-based streamline methods searching for the most likely paths
from a seed point to anywhere in the brain [3,16] suffer from path-length dependency, i.e. long paths are harder to estimate than short ones. When two endpoints
are known, SPT avoids this by looking for the shortest connecting path. This
ensures that long and short tracts are estimated with the same sample size. As
endpoints are not always known, these are complementary methods.
Solving Riemannian SPT. In Riemannian SPT, the diffusion tensor Dp at
p is commonly interpreted as the inverse of a Riemannian metric Mp [14, 15].
Geodesics returned by this ad hoc interpretation are typically “too straight” orare attracted to high diffusivity cortical spinal fluid (CSF) regions. Hao et al. [11]suggest a per-voxel scaling Mp = αp D−1
p that avoids this issue, while Fuster et
al. [9] show that the adjugate metric Mp = det(Dp )D−1
p gives a metric underwhich free Brownian motion matches the diffusion implied by Dp . This reduces
attraction to CSF regions [9, 19], so we use the adjugate metric.
Given a Riemannian metric Mp , the Riemannian geodesic c : [0, 1] → R3
connecting two points x and y is found by solving the geodesic ODE
c̈d (t) = fd (t, c, ċ) = −Γ|d · (ċ(t) ⊗ ċ(t)) , d = 1, . . . , D = 3,
D2
∂ vec Mp
1 X h −1 i
Mc(t)∈ RD ×1 ,
Γd =
2
∂pk
d,k
p=c(t)

(1)
(2)

k=1

subject to the boundary conditions c(0) = x and c(1) = y. Here ⊗ is the
Kronecker product, c(t) = [c1 (t); c2 (t); c3 (t)] is a point on the shortest path c,and ċ and c̈ denote the first and second derivative along the path, respectively.
Estimating a Random Geometry. The Riemannian adjugate metric Mp isestimated from finite noisy data and should be considered a stochastic variable.
This, however, complicates the geometric interpretation as the resulting object
is now a random Riemannian metric [13], and not a deterministic metric. Given
two endpoints x and y, our interest is in finding a connecting geodesic. Since
the metric is stochastic, there is a distribution of geodesics connecting x and y.
We approximately solve the geodesic ODE (1) using a probabilistic description of the uncertainty over f (t, c, ċ), the stochastic variable arising from Eq. 1
with a stochastic metric. We approximate f with a Gaussian distribution (using
the shorthand notation ft,c,ċ ≡ f (t, c, ċ))p(ft,c,ċ ) = N (ft,c,ċ ; mt,c,ċ , Ct,c,ċ ) .

(3)

The mean and covariance that form this approximation are computed from
md = −E(Γd )| (ċ ⊗ ċ)|

Cd2 = (ċ ⊗ ċ) cov(Γd ) (ċ ⊗ ċ) ,

and
2

(4)

Fig. 1. Left: A GP regressor(blue) and its derivative (orange). Right: A numerical solver
seeks the smooth function (blue)that best fits observed derivatives (black) while meeting theboundary constraints (red).

Predicted curve derivative

Predicted curvewhere mt,c,ċ = [m1 ; m2 ; m3 ] and Ct,c,ċ = diag C12 , C22 , C32 .
We estimate the mean and covariance of Γd empirically, subsampling the gradient directions used to generate the local diffusion tensors into S = 20 batches
each containing 80% of the directions. For each subsample we fit a diffusion
tensor field and compute Γd (2). Finally,
we estimate
  the moments|E(Γd ) ≈PSPS  (s)
(s)
(s)
1
1
Γd − E(Γd ) .
s=1 Γd and cov(Γd ) ≈ S−1
s=1 Γd − E(Γd )S

3Regressing an ODE

In the Riemannian setting, geodesics are found as the solution to the geodesic
ODE (1). This is a smooth curve c(t) : [0, 1] → R3 , which must be estimatednumerically. We use probabilistic numerics that solves the ODE using Gaussian process (GP) regression [12, 18, 20] We extend previous work [19] to handle
uncertainty in the ODE due to a noisy metric.
GP Regression. A GP c(t) ∼ GP(µ(t), k(t, u)) [17] is a probability measureover real-valued functions c : R → R such that any finite restriction to function
values {c(tn )}N
n=1 has a Gaussian distribution. GPs are parameterized by a meanfunction µ : R → R and a covariance function k : R × R → R that determines
the regularity of the paths. GPs are closed under linear transformationsp(c) = GP(c; µ, k)⇒p(Ac) = GP(Ac; Aµ, AkA| ),(5)

where A| denotes application of operator A to the left. Given observations{t, Y} = {(t1 , y1 ), . . . , (tN , yN )} of likelihood p(yi | ti ) = N (yi ; Ac(ti ), σ 2 I), theposterior over c is a Gaussian process GP(c̃; µ̃, k̃) withµ̃(t) = µ(t) + k(t, t)A| (Aktt A| + σ 2 I)−1 (Y − Aµ(t))
k̃(t, u) = k(t, u) − k(t, t)A| (Aktt A| + σ 2 I)−1 Ak(t, u),(6)where (ktt )ij = k(ti , tj ) is the N × N covariance of input locations [17, §2.2],
and similarly for k(t, t). Differentiation is a linear operation so (by Eq. 5) a GP
belief over c implies a GP belief over ∂c = ċ as well (see left panel of Fig. 1).
Beliefs over multi-output functions c(t) = [c1 (t); c2 (t); c3 (t)] can be constructed through vectorization. If the covariance structure is assumed to factorize
between inputs and outputs, cov(ci (t), cj (u)) = [V]ij k(t, u), for covariance V,
then the belief over c can be written as p(c) = GP(c; µc , V ⊗ k).
3

GP ODE solvers. Numerical ODE solvers evalute the ODE at finitely many
points and estimate a smooth curve that fits these observations along with either
initial or boundary values (see right panel of Fig. 1). This is statistical regression.
This view gives rise to probabilistic ODE solvers [5,12,18–20] implemented using
GP regression. In these solvers, uncertainty represents the approximation of notevaluating the ODE at the true solution, but at the current approximation.
At runtime, the solver repeatedly uses the current posterior mean estimatec̃(ti ) for the true solution c(ti ) to construct approximate noisy “observations”
yi = f (ti , c̃(ti ), ∂c̃(ti )) = c̈(ti ) + ηi .

(7)This describes that c̃(ti ) is only an approximation to the true solution c(ti ) as theobservation yi is corrupted by an error ηi . Here, ηi is assumed to be Gaussian,and the current uncertainty over c is propagated through the algorithm:
At step i, assume a posterior p(c) = GP(c; µ̃i , k̃ i ). We construct the estimatesic̃ and ∂c̃i as the current “best guess”, the mean [c̃i ; ∂c̃i ] = [µ̃i (ti ); ∂ µ̃i (ti )]. An
estimate for the error of this approximation is provided by the local variance∂ k̃i (ti ,t) ik̃ i (ti , ti )∂tc̃ (ti )t=ti 
=: Σ i .
(8)cov=  ∂ k̃i (t,ti )∂ 2 k̃i (t,u)∂c̃i (ti )
∂t

t=ti

∂t∂u

t,u=ti

Assuming we have upper bounds U > ∂f/∂c and U̇ > ∂f/∂ ċ on the gradients off , Σ i can be used to estimate the error on yi as (U, U̇)| Σ i (U, U̇) =: Λi , i.e.
ηi ∼ N (0, Λi ). This gives an observation likelihood functionp(yi | c̈(ti )) = N (yi ; c̈(ti ), Λi ).

(9)

The belief is updated with Eq. 6 to obtain µ̃i+1 , k̃ i+1 , and the process repeats.
The repeated extrapolation to construct c̃(ti ) in the GP solver is similar to
explicit Runge-Kutta methods as it defines a Butcher tableau. For some GP
priors, the posterior mean even coincide with results from such methods [18].
Regressing a Noisy ODE. When the metric is uncertain, the geodesic ODE
can only be evaluated probabilistically as p(f ). To handle this situation, we extend the probabilistic ODE solver to cope with Gaussian uncertain observations.
Due to noise, the curvature yi = f (ti , c̃i , ∂c̃(ti )) can only be estimated up to
Gaussian noise p(yi | mt,c̃,∂c̃ , Ct,c̃,∂c̃ ) = N (yi ; mt,c̃,∂c̃ , Ct,c̃,∂c̃ ). As the normal
distribution is symmetric around the mean, this is a likelihood for yi ,p(mt,c̃,∂c̃ | yi , Ct,c̃,∂c̃ ) = N (mt,c̃,∂c̃ ; yi , Ct,c̃,∂c̃ ).

(10)Since both likelihoods (9, 10) are Gaussian, the latent yi can be marginalized
analytically, giving the complete observation likelihood for mt,c̃,∂c̃
p(mt,c̃,∂c̃ | c̈(ti )) = N (mt,c̃,∂c̃ ; c̈i , Cti ,c̃i ,∂c̃i + Λi ) .

(11)

Solutions to noisy ODEs are then inferred by replacing Eq. 9 with Eq. 11, usingmt,c̃,∂c̃ = E[ft,c,ċ ] in place of the (inaccessible) function evaluations yi = ft,c̃,∂c̃ ,such that the approximation error is modeled by the additive uncertainty Ct,c̃,∂c̃ .
4The ODE solver rely on two noise terms: ηi in Eq. 7 captures the numerical error, while C in Eq. 10 describes the uncertainty arising from the data.
While the terms are structurally similar they capture different error sources; e.g.
changes in the number of grid points imply a change in η, while C is unaffected.
Both noise terms are approximated as Gaussian to ensure efficient inference.
Adaptation to Tractography. To compute geodesics in DTI we subsample
the diffusion gradients, estimate the mean and covariance of Γd , and finallysolve the noisy geodesic ODE (1) numerically.
For the prior covariance we use
a Gaussian kernel k(ti , tj ) = exp −(ti − tj )2 /(2λ2 ) with length scale λ. To
estimate V and λ we maximize their likelihood with gradient descent on thepositive definite cone. We initialize the prior mean with the Dijkstra path under
the deterministic metric, which we pre-process with a GP smoother.

4

Experimental Results

Tractography was performed on pre-processed diffusion data of 20 subjects fromthe Q3 release of the Human Connectome Project (HCP) [7, 8, 10, 21]. The
pre-processed HCP diffusion data contains 270 diffusion directions distributed
equally over 3 shells with b-values = 1000, 2000 and 3000 s/mm2 [21]. Diffusion
directions were uniformly subsampled into 20 batches each containing 80% of
the directions, where DTI tensors were computed with dtifit [2]. Segmentation was performed with FAST [23]. The cortico-spinal tract (CST) and inferior
longitudinal fasiculus (ILF) used for experiments were obtained from the probabilistic expert-annotated Catani tract atlas [4]. ROI atlases were constructed in
MNI152 “template space” by overlapping the tract atlas with regions from theHarvard-Oxford atlas [6]. The CST ROIs are the overlap with the brainstem,
the hippocampus and the amygdala for one region, and the overlap with thesuperior frontal gyrus, the precentral gyrus and the postcentral gyrus for thesecond region. The ILF ROIs are the overlap with the temporal pole for one
region, and the overlap with the superior occipital cortex, the inferior occipitalcortex and the occipital pole for the second region. The constructed ROIs were
warped from “template space” to “subject space” using warps provided by HCP.
As a first illustration, Fig. 2 shows the density of a single geodesic within the
CST projected onto a slice. The center column shows that the geodesic density
roughly consists of two certain vertical line segments with an uncertain connection between them (green box). The bottom row shows the standard deviation of
Γd . The geodesic uncertainty appears related to data noise. This is not attainedwith other GP ODE solvers [19] as they model constant observation noise.
Next we sample 250 endpoint pairs and compute the corresponding geodesics.
Fig. 3 shows the resulting heatmaps. The GP solution provides a more coherent
picture of the tract compared to the picture generated with Dijkstra’s algorithm.
To compare solution qualities we compute the set of voxels each geodesic
passes through. Taking the Catani atlas as a reference we measure the percentage
of voxels which are classified as part of the tract by at least one expert. Figure 4
shows the results for 20 subjects. In ILF the median accuracy is comparable
5GP density
std(Γ1 , Γ2 , Γ3 )Fig. 2. Top row: The density of a single GP geodesic under the random metric. The
density heatmap is projected into axis-aligned slices; the background image is the
expected metric trace; and the outline is where at least one expert annotated the CST.
Bottom row: Standard deviation of Γd at the same slices as the top row.
Dijkstra in CST

GP in CST

Dijkstra in ILF

GP in ILF

Fig. 3. Example shortest paths in the CST and ILF using both Dijkstra’s algorithm
on a deterministic metric, and a GP solver with a random Riemannian metric.
ILF

CST

Dijkstra
GP

0

0.5

1

0.2

0.4

0.6

0.8

1

Fig. 4. Agreement with at least one expert in the Catani atlas as estimated by Dijkstras
algorithm and the GP solver.

6

for Dijkstra and GP paths, but the GP error bars are significantly smaller. For
the CST, we observe similar quality results for Dijkstra and GP paths. This is
expected, as the ILF is generally considered the harder tract to estimate.

5

Discussion and Conclusion

SPT methods are advantageous for tracts that connect two brain regions, e.g. forstructural brain networks. However, a long-standing problem is that SPT only
finds a single connection with no insight to its uncertainty. This is in contrast to
walk-based methods, which are often probabilistic in nature. In this paper, we
provide a first fully probabilistic SPT algorithm that models stochastic diffusion
tensors and returns a distribution over the shortest path. While uncertainty
propagates in probabilistic walk-based methods, the uncertainty of probabilistic
SPT only depends on local data uncertainty, not the seed point distance.
Through experiments we visualize the estimated geodesic densities between
brain regions, and validate that the estimated geodesic densities are less certain
in areas where the estimated diffusion tensor is uncertain. We see that the visualized geodesic densities from the probabilistic SPT yield smoother and more
coherent tracts than the corresponding Dijkstra solutions.
While Riemannian SPT only applies to second order tensor models, a Finsler
geometry framework has emerged enabling continuous SPT with higher ordertensors for HARDI data [1]. Our proposed numerical tools also extend to thesemodels, and can provide a probabilistic interpretation of Finsler models. Further
future work includes shape analysis on the resulting estimated geodesic densities,
which are suitable for GP-based tract shape analysis [22].
Acknowledgements. Data were provided [in part] by the Human Connectome
Project, WU-Minn Consortium (Principal Investigators: David Van Essen and
Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centersthat support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. S.H. is funded inpart by the Danish Council for Independent Research (DFF), Natural Sciences.
A.F. is funded in part by the DFF, Technology and Production Sciences.

References
1. Astola, L., Jalba, A., Balmashnova, E., Florack, L.: Finsler streamline tracking withsingle tensor orientation distribution function for high angular resolution diffusion
imaging. Journal of Mathematical Imaging and Vision 41(3), 170–181 (2011)
2. Basser, P.J., Mattiello, J., LeBihan, D.: Estimation of the effective self-diffusion
tensor from the NMR spin echo. J Magn Reson B 103(3), 247–254 (1994)
3. Basser, P.J., Pajevic, S., Pierpaoli, C., Duda, J., Aldroubi, A.: In vivo fiber tractography using DT-MRI data. Magnetic resonance in medicine 44(4), 625–632 (2000)
4. Catani, M., de Schotten, M.T.: A diffusion tensor imaging tractography atlas forvirtual in vivo dissections. Cortex 44(8), 1105–1132 (2008)
5. Chkrebtii, O., Campbell, D., Girolami, M., Calderhead, B.: Bayesian uncertainty
quantification for differential equations. arXiv stat.ME 1306.2365 (Jun 2013)7

6. Desikan, R.S., et al.: An automated labeling system for subdividing the humancerebral cortex on MRI scans into gyral based regions of interest. NeuroImage
31(3), 968–980 (2006)
7. van Essen, D., et al.: The Human Connectome Project: a data acquisition perspective. NeuroImage 62(4), 2222–2231 (2012)
8. van Essen, D., et al.: The WU-Minn Human Connectome Project: an overview.
NeuroImage 80, 62–79 (2013)
9. Fuster, A., Tristan-Vega, A., Haije, T., Westin, C.F., Florack, L.: A novel Riemannian metric for geodesic tractography in DTI. In: Computational Diffusion MRI
and Brain Connectivity, pp. 97–104. Math. and Visualization, Springer (2014)
10. Glasser, M., et al.: The minimal preprocessing pipelines for the human connectome
project. NeuroImage 80, 105–124 (2013)
11. Hao, X., Zygmunt, K., Whitaker, R.T., Fletcher, P.T.: Improved segmentation ofwhite matter tracts with adaptive Riemannian metrics. Medical Image Analysis
18(1), 161–175 (2014)12. Hennig, P., Hauberg, S.: Probabilistic solutions to differential equations and theirapplication to Riemannian statistics. In: AISTATS. vol. 33. JMLR, W&CP (2014)
13. LaGatta, T., Wehr, J.: Geodesics of random riemannian metrics. Communicationsin Mathematical Physics 327(1), 181–241 (2014)
14. Lenglet, C., Deriche, R., Faugeras, O.: Inferring white matter geometry from diffusion tensor MRI: Application to connectivity mapping. In: Pajdla, T., Matas, J.
(eds.) Computer Vision – ECCV. Lecture Notes in Computer Science, vol. 3024,
pp. 127–140. Springer (2004)
15. O’Donnell, L., Haker, S., Westin, C.F.: New approaches to estimation of whitematter connectivity in diffusion tensor MRI: Elliptic PDEs and geodesics in atensor-warped space. In: Dohi, T., Kikinis, R. (eds.) Medical Image Computing
and Computer-Assisted Intervention – MICCAI. pp. 459–466. Lecture Notes in
Computer Science, Springer (2002)
16. Parker, G., Haroon, H., Wheeler-Kingshott, C.: A framework for a streamlinebased probabilistic index of connectivity (PICo) using a structural interpretationof MRI diffusion measurements. J. Magn. Reson. Imaging 18(2), 242–254 (2003)
17. Rasmussen, C., Williams, C.: GPs for Machine Learning. MIT Press (2006)
18. Schober, M., Duvenaud, D., Hennig, P.: Probabilistic ODE solvers with RungeKutta means. In: NIPS (2014)19. Schober, M., Kasenburg, N., Feragen, A., Hennig, P., Hauberg, S.: Probabilistic shortest path tractography in DTI using Gaussian Process ODE solvers. In:
Golland, P., Hata, N., Barillot, C., Hornegger, J., Howe, R. (eds.) Medical Image Computing and Computer-Assisted Intervention – MICCAI. Lecture Notes in
Computer Science, vol. 8675, pp. 265–272. Springer (2014)20. Skilling, J.: Bayesian solution of ordinary differential equations. Maximum Entropy
and Bayesian Methods, Seattle (1991)
21. Sotiropoulos, S.N., et al.: Effects of image reconstruction on fiber orientation mapping from multichannel diffusion MRI: reducing the noise floor using SENSE. Magn
Reson Med 70(6), 1682–89 (2013)
22. Wassermann, D., Rathi, Y., Bouix, S., Kubicki, M., Kikinis, R., Shenton, M.E.,
Westin, C.F.: White matter bundle registration and population analysis basedon Gaussian Processes. In: Proceedings of the 22nd International Conference on
Information Processing in Medical Imaging. pp. 320–332. Springer (2011)
23. Zhang, Y., Brady, M., Smith, S.: Segmentation of brain MR images through ahidden Markov random field model and the expectation-maximization algorithm.
IEEE Trans Med Imaging 20(1), 45–57 (2001)

8