Geodesic Uncertainty in Diffusion MRI

Abstract

We study theoretical and operational issues of geodesic tractography, a geometric methodology for retrieving biologically plausible neural fibers in the brain from diffusion weighted magnetic resonance imaging. The premise is that true positives are geodesics in a suitably constructed metric space, but unlike traditional first order methods these are not a priori constrained to connect nongeneric points on subdimensional manifolds, such as the characteristics in traditional streamline methods. By virtue of the Hopf-Rinow theorem geodesic tractography furnishes a huge amount of redundancy, ensuring the a priori existence of at least one tentative fiber between any two points and permitting additional tractometric and data-extrinsic constraints for (fuzzy or crisp) classification of true and false positives. In our feasibility study we consider a hybrid paradigm that unifies existing ideas on tractography, combining deterministic and probabilistic elements in a way naturally supported by metric geometry. Particular attention is paid to an analytical prediction of geodesic deviation on numerically computed geodesics, a ‘tidal’ effect induced by small perturbations resulting from data noise. Taking these effects into account clarifies the inherent uncertainty of geodesics, while simultaneosuly offering a dimensionality reduction of the tractography problem.

1 Introduction

A hydrogen nucleus (proton) behaves like a tiny magnet when interacting with an external magnetic field. In this situation the Zeeman splitting effect discloses two quantum eigenstates of its nuclear magnetic moment, distinguishable by a relative energy gap. These states are known in the trade as “spin up”, or “parallel”, and “spin down”, or “anti-parallel”. The attributes refer to the relative alignment with the external field, but should be taken with a grain of salt for quantum systems of individual nuclei. In a typical macroscopic system, however, the Zeeman effect induces a relatively small yet measurable “classical” magnetization [Bloch (1946); Cowan (1997)], for which these attributes can be taken literally.

Magnetic resonance imaging (MRI) allows one to manipulate the magnetization of a hydrogen spin system in a way that enables non-invasive in-vivo acquisition of informative localized signals, i.e. “images”. Although tuned in on water-bound hydrogen, because of its abundance in the human body, the recorded signals are affected by the physicochemical environment and therefore carry (indirect) anatomical and/or functional information depending on the measurement protocol. MRI is a versatile modality with endless options for image formation, cf. (Brown et al. (2014)).

The empirical feature of interest for our purpose is the anisotropic diffusivity of mobile (water-bound) hydrogen spins in the human brain¹. The MRI protocol that allows one to infer this from measurable nuclear magnetic resonance and relaxation effects is known as diffusion weighted imaging (DWI) (Torrey, 1956; Moritani et al., 2009). Insight in local water diffusivity patterns is particularly interesting for studying fibrous tissues, since it is stipulated that anisotropic media impart non-random barriers for diffusion, inducing relatively long mean free paths along “preferred orientations” (candidate fibers). In particular, DWI is the only non-invasive in-vivo technique for probing white matter pathways (bundles of myelinated axons or nerve fibers, aka “tracts”) in the brain. As such it plays a pivotal role in connectomics, the 21st century’s grand challenge that aims to disclose structure and function of the human brain.

DWI produces massive, non-visual data, calling for the right questions to be posed in order to render such data insightful in congruity with our visual inclination. An important question pertains to tractography, the inverse problem that aims to infer white matter pathways from the observed local diffusivity patterns. Such tracts furnish the wiring between various grey matter areas containing nerve cell bodies, the brain’s computational units so to speak. For this reason tractography also plays an essential role in clinical applications, such as neurosurgery, in which spatial maps of anatomical tracts and functional networks are relevant for effective patient management.

However, tractography has inherent limitations that should be taken into account as well, e.g. inter- and intra-rater variablity makes comparisons of tractograms a challenging endavour (Maier-Hein et al., 2017; Schilling et al., 2019). In a recent paper by Schilling et al. (Schilling et al., 2020) the authors argue that tractography is anatomically plausible if furnished by appropriate constraints, notably end point conditions and no-goes. This observation argues for a tractography paradigm that is compatible with end point conditions imposed by the user and that does not exclude the possibility of multiple connections, two properties not shared by classical tracing (streamline) methods based on diffusion tensor imaging (DTI). Moreover, no-go conditions should be at the discretion of an expert, and not a result from a priori methodological limitations. An example is the fractional anisotropy (FA) stopping criterion of classical streamline methods, which results in incomplete trajectories due to destructive interference of anisotropic diffusivities at fiber crossings or severely reduced effective anisotropies in grey matter.

The conceptual issues alluded to above are gracefully handled by geodesic tractography, the mathematical details of which are explained in section 2. It is important to realize that this is not a one-on-one replacement of streamline tractography, but a genuine generalization that (intentionally) raises a new problem: Any two points in the brain are connected by at least one geodesic (“geodesic completeness”). Thus geodesics should not be confused with (streamlines intended to reflect) biological fibers, but should instead be seen as elements of a highly redundant data representation. The premise is, of course, that these can be pruned using data intrinsic and external knowledge to (hopefully) produce the desired true positives without incurring false negatives. That is, the assumption is that biological tracts are geodesics, but not vice versa. Interestingly, one can show that classical streamlines emerge as geodesics when the two minor eigenvalues of the DTI tensor are set to zero. This Gedanken experiment connects the two tractography paradigms conceptually, but also shows (by an argument of continuity) that there exist “streamline-like” geodesics in brain regions characterized by high FA, with “streamline-unlike” infillings in isotropic regions.

For the reasons above we focus on geodesic tractography, and present some feasibility studies in Section 3. More specifically, we aim to establish a conceptual and operational framework geared towards clinical needs in a neurosurgical workflow for brain tumor patients, in collaboration with the Department of Neurosurgery of Elisabeth Tweesteden Hospital (ETZ), a highly specialised referral center for brain tumor surgery in Netherlands (Meesters et al., 2019; Meesters et al., 2021; Rutten et al., 2014). Our approach is motivated by conceptual as well as practical demands in this context, viz. 1) to overcome the rigidity and shortcomings of classical tracing methods in return for a more versatile one, and 2) to connect to clinically viable MRI protocols, notably DTI, while keeping tabs on the development of more sophisticated DWI models for future refinement.

We sketch the geodesic tractography rationale in all generality, which in principle requires so-called Finsler geometry, but work out details for its application to DTI only. This case is mathematically considerably less mind-boggling since it allows us to stay within the realm of well-established Riemannian geometry. Another practical issue is the effect of noise and errors in (automatic or manual) seeding, which can be naturally accounted for by considering the concept of geodesic deviation, furnishing geodesic tracks with volumetric counterparts in the form of tubular uncertainty neighbourhoods (“geodesic tubes”) (Sengers et al., 2021).

2 Theory

Taking the empirical adequacy of DWI for tractography for granted, something must be optimal along biologically meaningful curves. This suggests a variational approach, in which tracts arise as extrema of some functional (Arnold, 1989; Lovelock and Rund, 1988; Rund, 1973; Sagan, 1992; Weinstock, 1974).

In geodesic tractography, the functional of interest is constructed on the basis of a norm or metric. The relevant geometry is known as Riemann-Finsler geometry, or Finsler geometry for short (Bao et al., 2000; Cartan, 1934; Shen and Shen, 2016; Szilasi et al., 2014). Inner product induced norms are quite special, and often used out of habit rather than necessity. Such norms fall within the much narrower realm of Riemannian geometry, an area much better understood and (therefore) exploited than the general framework (do Carmo, 1976; do Carmo, 1993; Cartan, 1963; Jost, 2011; Koenderink, 1990; Misner et al., 1973; Spivak, 1975). The essential constraint is that inner product norms are square roots of quadratic forms, characterized by a small number of degrees of freedom, viz. six per point in three spatial dimensions for specifying the independent parameters of the defining positive-definite symmetric Gram matrix. Such norms are especially suited for applications in which the observable of interest is compatible with the quadratic restriction, notably diffusion tensor imaging (DTI). General DWI models, however, do not have this limitation. Extension of the variational approach based on metric geometry beyond DTI would therefore render a genuine Finslerian approach natural and compulsory.

Although Riemannian geodesic tractography does not cover traditional “streamline tractography”, the latter may be viewed as a singular limit that can be inferred from a sequence of Riemannian metrics. In turn, Finsler geodesic tractography may be related to an orientation-parametrized family of Riemannian metrics (Astola, 2010; Astola and Florack, 2011; Astola et al., 2011; Florack and Fuster, 2014; Florack et al., 2015; Dela Haije, 2017; Dela Haije et al., 2019). In this sense Riemannian geodesic tractography plays an intermediate role, being a generalisation of the original streamline case and a special instance of the Finslerian case. Its limitations are offset by the inspiration that can be drawn from the vast body of literature on Riemannian geometry, as well as by the particular clinical appeal of DTI.

Positive definite matrix fields play a central role in DWI, and in DTI in particular, and so do, a fortiori, positivity preserving operators. In a differential geometric approach, differential operators are of prime interest. Unfortunately, positivity is not manifest in “standard” calculus. A more suitable framework for differentiation and integration is provided by multiplicative calculus, which, in the non-commutative case of positive symmetric matrices, is a non-trivial extension in which positivity is manifest (Bashirov et al., 2008; Florack and van Assen, 2012; Gill and Johansen, 1990).

To make this article self-contained the most relevant mathematical concepts are loosely explained, de-emphasizing proofs, in subsections 2.1 (calculus of variations), 2.2 (multiplicative calculus) and 2.3 (Riemann-Finsler geometry). Literature pointers should help the reader interested in more detail.

2.1 Calculus of Variations

In the Lagrangian formalism of the calculus of variations the point of departure is an action functional,

The integrand is called the Lagrangian. In many cases Lagrangians only depend on the m-dimensional independent variable x, the n-codimensional dependent variable ϕ, and its first order x-derivatives, ∇ϕ:

The variational principle relates the physical manifestation of a Lagrangian system to critical points ϕ_crit of the action functional, formally defined in terms of a vanishing functional derivative:

Although the left hand side can be rigorously defined, it suffices to appreciate its intuitive meaning as the “rate of change” of the action functional S(ϕ) with respect to any small perturbation δϕ of the dependent variable ϕ relative to a fiducial instance ϕ_crit. Eq. 3 yields a system of n differential equations known as the Euler-Lagrange system.

Our case of interest will be m = 1, n = 3. In this situation the domain of definition is a one-parameter interval, and one usually writes t (“time”) or s (“length”) instead of x, and L instead of . Since in this case we have an ordinary rather than partial derivative operator, the ∇-symbol now represents d/dt or d/ds. Our codomain will be a 3-dimensional tangent space, representing the linear space of instantaneous “velocities’” within (any tangent space of) 3-dimensional space. For this reason one typically uses the spatial coordinate symbol x instead of ϕ. Its first order derivatives are then conveniently written either as (Newton’s dot-notation) or as x′ = dx/ds (Lagrange’s prime notation). In tractography the independent parameter serves as an arbitrary label for disambiguating points along a tract; the connotation with “time” or “length” is metaphorical.

Thus the action functional that concerns us here takes the formin which x = x(t) represents the coordinate triple along a parametrized tract over an interval . Stationarity entailswhich yields an Euler-Lagrange system in the form of three ordinary differential equations, one for each component , i = 1, 2, 3, of the tangent vector relative to a fiducial local frame of basis vectors {e₁, e₂, e₃}:

The precise form of the Lagrangian is constructed on the basis of geometric considerations and will be specified later. It is important to realize that the spatial variability of the local frame (more precisely, its parallel transport²) affects the form of the Euler-Lagrange equations.

For a thorough introduction to the calculus of variations, cf. Lovelock and Rund (1989). Rund elaborates on the profound implications of one-homogeneous Lagrangians Rund (1973). The book by Neuenschwander (2011) is particularly interesting in case the Lagrangian exhibits certain symmetries, explaining Noether’s seminal work on this issue.

2.2 Multiplicative Calculus

Multiplicative calculus (Bashirov et al. (2008)) provides a natural framework in problems in which positive images or positive symmetric matrix fields and positivity preserving operators are of interest. In the commutative case it is a convenient, albeit arguably redundant framework, which can be restated in terms of standard (additive) calculus via a log-exp commuting diagram (Arsigny et al., 2006; Fillard et al., 2007; Florack and Astola, 2008; Pennec et al., 2006). Figure 1 illustrates this for the case involving a single independent variable , with the multiplicative derivative of a positive function f defined as

FIGURE 1

The multiplicative antiderivative is then naturally defined as

The left hand side notation is due to Gill (Gill and Johansen, 1990) in analogy with Leibiz “continuous-sum” symbol.

For a product antiderivative or indefinite product integral we havemimicking its familiar standard counterpart

Definite product integrals are introduced via a spatial partitioning and limiting procedure,akin to the standard Riemann sum approximation in the additive case,in which △x_i = x_i − x_i−1, ξ_i ∈ (x_i−1, x_i) and x₀ = a, x_N = b. The relationship between Eq. 9, and Eq. 11 is formalized by the fundamental theorem of multiplicative calculus:recall the well-known classical counterpart relating Eq. 10 and Eq. 12:

Numerical approximations of Eq. 11 are obtained by partitioning the integration domain and suppressing the limiting procedure.

It requires minor efforts to generalize foregoing results to the multivariate case. The non-commutative case, on the other hand, is highly nontrivial and not without ambiguity. The complication is apparent from the ambiguous ordering of non-commuting factors on the right hand side of Eq. 11. There are multiple options to establish an unambiguous definition, cf. Gantmacher (2001) and Slavík (2007)) We will briefly discuss some and subsequently pick the one suited for our purpose.

The log-exp framework may seem the most straightforward option, defining the multiplicative derivative of a positive symmetric matrix function according to Eq. 7, in which log and exp are replaced by their well-defined matrix counterparts. The chain rule for matrix functions, however, complicates matters. Declining from the assumption of commutativity we have, e.g.for an arbitrary (square) matrix field g. For a positive symmetric matrix field f we furthermore havein which I is the identity matrix. Both identities can be readily proven with the help of the defining Taylor series expansions of the exp and ln functions, cf. Higham (2008). In the commutative case (only) Eqs 15 and 16 reproduce the familiar chain rule formulasandwell-known from classical calculus. In view of this one may avoid the cumbersome expression of Eq. 16 resulting from the log-exp trick of Eq. 7 by defining the multiplicative derivative in a different way, viz. by replacing the expression (ln f(x))′ on the right hand side in Eq. 7 formally by one of the alternative expressions in Eq. 18. As these two options (or any alternative construct obtained by sandwiching factorizations of f(x)⁻¹ = f(x)^α−1f(x)^−α around f′(x)) are clearly different, the application context should determine one’s bias. We will adhere to the following definition:

It can be shown that, with this choice, Eq. 11 is applicable, provided we use a right-to-left ordering of factors on the right hand side:

The reason for choosing Eqs 19 and 20 is that these provide the correct derivative/antiderivative for recasting a linear matrix differential equation f′ = A f, with positivity constraint on f—which will be of interest later, notably in Eq. 48 on page 13—, into its natural multiplicative counterpart with built-in positivity preservation, viz. f* = exp A. The solution then takes the trivial form of an antiderivative, , according to³Eq. 20.

Another application of multiplicative calculus relevant for our purpose is the so-called log-Euclidean paradigm for symmetric positive-definite matrix functions⁴ (Arsigny et al., 2006; Fillard et al., 2007; Pennec et al., 2006). We consider it in the context of multiscale representations (Florack and Astola, 2008). Figure 2 shows how to construct a consistent multiscale representation f_t of a raw symmetric postive-definite matrix function f ≐ f₀, with t > 0 denoting (quadratic) scale, according to the following multiplicative convolution recipe:in which ϕ_t denotes the normalized isotropic Gaussian of (quadratic) scale t and ∗ acts component-wise:

FIGURE 2

This scheme ensures closure and commutativity of blurring and inversion: If and , with g = f⁻¹, then at all scales t > 0.

2.3 Riemann-Finsler Geometry

Riemann-Finsler geometry is a generalization of Riemannian geometry. Both are metric geometries, concerned with spaces endowed with a distance concept. The generalization entails the replacement of an inner product induced norm (square root of a quadratic form reflecting the Pythagorean rule) by a general norm. For this reason Riemann-Finsler geometry is sometimes referred to as ‘Riemannian geometry without the quadratic restriction’. First studied by Finsler (1918), it has matured towards its present form mainly due to Cartan (1934). The mathematics is rather mind-boggling. For modern in-depth accounts cf. Bao et al. (2000), Shen and Shen (2016), and Antonelli et al. (1993); (Antonelli and Zastawniak, 1999) in relation to some applications in physics and biology. Riemann-Finsler geometry is also attracting considerable attention in the context of (quantum) gravity, see for example (Fuster and Pabst, 2016; Gibbons et al., 2007; Girelli et al., 2007; Pfeifer and Wohlfarth, 2012).

For our purpose a norm is needed to introduce an abstract notion of “length” of a curve, one for which units of length are proportional to local mean free path lengths of diffusing water molecules in the underlying fibrous tissue, as observed via DWI. This implies that if local mean free paths are indeed longest along biological tracts, then such tracts would define “locally shortest” paths, or geodesics, along a Riemann-Finsler manifold. The attribute “local” means that any sufficiently small perturbation of a fiducial geodesic with fixed endpoints would incur an increase of length. This leaves room for multiple geodesics connecting the same pair of endpoints. The Hopf-Rinow theorem (Jost, 2011) ensures geodesic completeness, i.e. the existence of a “minimal geodesic”, the length of which corresponds to (or, if you want, defines) the distance between its endpoints. In the case of multiple local minima between two endpoints, the global one is not necessarily the biologically most plausible one.

The Riemann-Finsler manifold is thus constructed in order to turn the tractography problem into a geodesic problem, which falls within the scope of the calculus of variations, recall section 2.1. The inherently positive nature of distance naturally calls for multiplicative calculus, recall section 2.2.

The most general functional expressing the Riemann-Finsler length of a tentative tract, arbitrarily parametrized as x = x(t), has the following form, recall Eq. 4:in which the general Lagrangian has been replaced by the so-called Finsler norm , expressing the Riemann-Finsler length of the tangent vector at point x in a parameter invariant way. The latter not only means that it does not explicitly depend on the curve parameter t ∈ I, but also that it must be absolutely homogeneous of degree one, i.e.⁵for any tangent vector rescaling by .

One can show that, under mild conditions besides Eq. 24, there exists a positive-definite symmetric Riemann-Finsler metric tensor , such that⁶

The relation between and is, remarkably, one-to-one, with

The Riemann-Finsler metric tensor is homogeneous of degree zero,

Thus a Riemann-Finsler metric has an orientation dependency, but magnitude and “polarity” of the tangent vector are immaterial.

The quadratic restriction of Riemannian geometry entails a severe restriction on the Finsler norm, viz. independent of orientation. Equation 26 may be interpreted as a family of Riemannian metrics, one for each orientation. All members of this family clearly coincide in the Riemannian case, so that Eq. 25 reduces to its Riemannian counterpart, the Pythagorean form

The Euler-Lagrange equations for Eq. 25 are known as the geodesic equations, and take the form

The right hand side can be seen as a pseudoforce acting along the geodesic and as such merely reflects the arbitrariness of “time”-parametrization, contributing to an artificial acceleration along the trajectory. It has no effect on the relevant (parameter-invariant) geometry of the curve and can be effectively removed via (nonlinear) reparametrization. Any new parameter thus obtained is then coined an affine parameter, yielding “constant speed” geodesics characterized by . The geodesic spray coefficients can be viewed as forces inducing an acceleration and bending of the trajectory that is completely determined by the inhomogeneous nature of the metric tensor, cf. Bao (Bao et al., 2000) or Shen and Shen Shen and Shen (2016) for technical details. In the Riemannian limit they assume a quadratic form in the velocities .

Figure 3 illustrates the rationale of constructing the action functional Eq. 23 as a data-induced curve length in the context of the variational principle for the Riemannian case applicable to DTI.

FIGURE 3

Apart from the geodesic equations, Eq. 29, we will be interested in “tidal effects’, i.e. the relative separation acceleration⁷ of neighbouring geodesics, aka geodesic deviation. This is relevant for understanding how stable a geodesic is with respect to small perturbations. By quantifying this one can analytically predict (to linear approximation) the behaviour of sufficiently narrow bundles of geodesics around any fiducial geodesic. Understanding geodesic deviation in a general metric space requires some fairly sophisticated mathematical machinery. It suffices here to convey the gist of it.

The geodesic deviation equations form a linear system of ordinary differential equations for a deviation vector ξ = ξ(t) defined along a given geodesic curve x = x(t). This vector can be seen as connecting points on two ‘infinitesimally’ neighbouring geodesics at the same parameter value t (one must, to begin with, first agree on an unambiguous affine parametrization). We have

The δ-derivative is a so-called covariant modification of the d-derivative from standard calculus, involving geometric ‘correction terms’ that depend on partial derivatives of the metric tensor with respect to both variables. The tensor field is one of the curvature tensors that can be defined on a Finsler manifold. In the Riemannian case Eq. 30 simplifies considerably, with reducing to the so-called Riemann curvature tensor (a function of x only). In the particular case of interest, with n = 3 spatial dimensions, it may be reduced further to involve only the symmetric Ricci tensor field R_ij (x), with six local degrees of freedom

The interested reader is referred to the literature for technical details (Bao et al., 2000; Shen and Shen, 2016; Rund, 1959).

2.4 Riemann-DTI Paradigm: Geodesic Tractography

DTI tensors are represented by symmetric positive definite matrices obtained by fitting their six degrees of freedom at a fiducial point x to the DWI local signal attenuation factor according to the anisotropic Stejskal-Tanner equation (Stejskal and Tanner (1965); Stejskal (1965)),in which q denotes the diffusion sensitizing magnetic gradient direction⁸ and b a dimensionful constant related to diffusion time and gradient magnitude used in the scanning protocol. A Riemannian metric may then be constructed proportional to the inverse of the DTI tensor (O’Donnell et al., 2002; Lenglet et al., 2004; Fuster et al., 2016)

The metric tensor is now only dependent on the position x and not on the orientation , so that the Riemannian length of a parameterized curve x = x(t) is given byrecall Eqs 23, 28. Recall that geodesic tractography entails finding curves x = x(t) for which S(x) is locally minimal. This endeavour may be tackled in two different ways, either by directly minimizing Eq. 34, or by solving its Euler-Lagrange equations, Eq. 29. In the Riemannian case at hand, and using an affine (“constant speed”) parametrization, these equations reduce toin which the so-called Christoffel symbols of the second kind (Spivak, 1975) are given by linear combinations of derivatives of g_ij:

Equation 35 may be reformulated in purely geometrical terms by defining a so-called covariant differential operator D (Spivak, 1975), given in terms of the standard differential operator d and Γ-correction terms (a special case of the δ-operator used in Eq. 30), accounting for the non-constant nature of g_ij (x). For a vector function v = v (t) defined along the curve x = x (t) we haveso that Eq. 35 is equivalent to

For disambiguation of its solution, Eq. 35 requires a pair of initial or boundary conditions

The existence of a geodesic constrained by either Eq. 39a or 39b is guaranteed by the Hopf-Rinow theorem (Jost, 2011). For this reason the x-manifold is called geodesically complete. Uniqueness is not provided for the boundary problem, as there may be multiple curves (locally) minimizing Eq. 34 for sufficiently distant endpoints x₀ and x_T allowing for multiple potential fiber tracts connecting these endpoints.

2.5 Riemann-DTI Paradigm: Geodesic Deviation

Let us consider a single geodesic x = x (t) obtained either by minimization of Eq. 34 or by integration of Eq. 35. In order to understand the stability of this geodesic with respect to small perturbations of its initial or boundary conditions, Eq. 39, we may employ perturbation theory in the form of the linear geodesic deviation equations, recall Eq. 30, restricted to the Riemannian case of interest. This allows us to predict the behaviour of neighbouring geodesics without explicitly solving their nonlinear geodesic equations, as long as these are sufficiently close to the fiducial geodesic x = x (t), as described in Sengers et al. (2021), in which families of neighbouring geodesics are defined via continuous perturbations of the side conditions, Eq. 39. Here we extend our uncertainty quantification by including the effect of DTI data noise, i.e. perturbations of the metric field g_ij (x) along (a full neighbourhood of) a fiducial geodesic, which requires a generalization of Eq. 30 by accounting for an inhomogeneous “force” term. We consider general perturbations of the metric field, as to capture the intrinsic data-induced uncertainty, arising when using the Riemann-DTI tractography paradigm. In practice, these perturbations include perturbations in the DTI tensor and the underlying DWI data.

Thus we consider two types of perturbations, those affecting the intitial or boundary conditions in Eq. 39 and those affecting the metric field, viz.respectively

The formal parameter 0 ≤ ɛ≪ 1 is dimensionless. We are interested in effects; terms are considered irrelevant and will be suppressed⁹.

For an arbitrary geodesic x = x (t) satisfying Eq. 35, consider the parameterized family of neighbouring tracks,induced by any of the perturbations in Eqs 40a, 40b, 41. The perturbative, vector-valued function y = y (t), defined along the geodesic x (t), must be such that the perturbed path represents itself a geodesic with respect to the metric for any sufficiently small value of ɛ. Thus the perturbed path, Eq. 42, satisfies Eq. 35 up to , with the replacement . The vectors y(t) are referred to as deviation vectors, measuring the extent to which deviates from . Substitution of (42) into Eq. 35, with , yields (cf. Eq. 30)where D/dt denotes the covariant derivative operator along x (t), recall Eq. 37, and is the Riemann curvature tensor given, in terms of the Christoffel symbols of Eq. 36, by

Writing the first order expansion of the Christoffel symbols asakin to Eq. 41, the inhomogeneous term in Eq. 43 is given by

In order to appreciate the Eqs 43–46, consider the hypothetical case of a homogeneous, noise free DTI tensor image, inducing a constant metric tensor field. In the absence of noise we have no metric perturbations, so that fⁱ = 0, while the x-independent nature of the metric also induces a vanishing Riemann tensor, , as well as a trivial covariant derivative, D/dt = d/dt. Equation 43 then reduces to , implying a linear evolution of the deviation vector along the geodesic xⁱ = xⁱ(t). Both the fiducial geodesic as well as all geodesics induced by this type of deviation are then straight paths by virtue of Eqs 35–38 and Eq. 42. This is what one would expect of a deviation vector in a flat space. In a curved space, with a nonzero Riemann tensor in Eq. 43, neighbouring geodesics may be locally attracted towards or repelled away from the fiducial one, depending on the nature of curvature. Random data noise, producing a nonzero force term fⁱ in Eq. 43, adds further jitter to the geodesic paths.

Defining the block matrix and the column vectors by

Equation 43 can be written as an inhomogeneous linear first order system of matrix differential equations

This equation may be solved with non-commutative multiplicative calculus (Bashirov et al., 2008; Florack and van Assen, 2012; Gill and Johansen, 1990) using the product integral defined in Eq. 20:

The block structure in M, F and Y induces a similar structure in the product integrals:with the help of which the solution y (t) to the initial value problem may be expressed compactly aswith

If one considers instead the boundary value problem, the initial tangent vector must be treated as an unknown and upon setting t = T can be solved for in terms of y₀, y_T and ∫Θ(T, σ) f(σ) dσ. The solution of the boundary value problem presents itself as

3 Experimental Results

So far the theory has been formulated in a continuous setting. In clinical applications of tractography we have a discretized DTI and ditto metric field. This raises the question of sub-voxel tensor field interpolation and discrete differentiation for evaluating Christoffel symbols, recall Eq. 36. These issues are simultaneously handled by using a multiscale representation of the field, Eq. 21, combined with Log-Euclidean interpolation (Arsigny et al. (2006)). For technical details, cf. Florack et al. (2021). Supplied with initial conditions (39a), we may solve for a geodesic by time-integrating Eq. 35 over a specified interval t ∈ (0, T), e.g. using the Runge-Kutta fourth order integration scheme. If instead boundary conditions (39b) are supplied, both direct minimization of Eq. 34 and a non-linear solver for Eq. 35 are viable solution strategies. Finally, we may discretize the explicit analytical solutions to obtain geodesic deviation vectors from Eq. 53.

To illustrate our perturbative approach, we perform experiments on a DWI dataset from the Human Connectome Project. In our experiments the metric tensor of our Riemann-DTI paradigm is the adjugate of the DTI matrix D, , cf. Fuster et al. (2016). We restrict ourselves to perturbations of the metric given by Eq. 41, while keeping boundary conditions fixed, Eq. 39b i.e. we only consider the influence of perturbations on the DTI tensor field. The effect of boundary variations of type Eq. 40a have recently been considered elsewhere (Sengers et al. (2021)). These two types of perturbations are invariably present in any tractography algorithm due to DTI domain and codomain uncertainties.

We simulate DTI field perturbations from D_ij by a bootstrapping procedure (Whitcher et al., 2008). Corresponding metric perturbations are then determined as according to Eq. 41 with ɛ = 1, recall footnote 10. Subsequently, the deviation y(t) along an arbitrary track x (t) can be determined by Eq. 53, with y (0) = y(T) = 0, resulting in a perturbed track . Recall that, apart from discretization, Eq. 53 is an analytical result, rendering this perturbation procedure quite efficient. This data noise simulation procedure is repeated so as to obtain a large collection of perturbed tracks.

For the sake of comparison we explicitly compute, for each simulated metric perturbation , the geodesic subject to and . We thus obtain two collections of tracks, viz. a set of (computationally intensive) explicitly computed geodesics , and a set of (analytical, thus computationally efficient) approximations to geodesics , where k = 1, … , N pertains to the track associated with the kth perturbed DTI field. In our experiments we take N = 100.

Figure 4 illustrates the behaviour of analytical perturbations vs. simulated geodesics. In the first case the analytically computed tracks almost completely coincide with the simulation results, while in the second case a small degree of misalignment is clearly observable. Even though the DTI perturbations are of similar order everywhere, the linear approximation of Eq. 42 apparently does not suffice in all cases, from which we may infer that the mathematical notion of a “sufficiently small” perturbation is case dependent.

FIGURE 4

In Figure 5 we extend the experiment of Figure 4 to a bundle of geodesics. For each of 4,237 (unperturbed) geodesics starting in the brain stem and ending in the precentral gyrus, their (simulated) perturbed geodesics and (analytical) approximations are determined, again with k = 1, … , N = 100, producing two sets of 423,700 tracks each. For each of these sets we construct isosurfaces of the densities obtained by counting the number of tracks passing through each voxel. Subsequently, we threshold the density maps and construct two additional ones, so called difference maps. The difference is set to ±1 for voxels in which fibers are present in only one of the two density maps, and to 0 otherwise. A significant overlap between the analytical and simulated isosurfaces can be observed, confirming the feasibility of our approach.

FIGURE 5

4 Discussion

We have presented an extension to the Riemann-DTI paradigm in the form of an uncertainty quantification due to DTI data noise. To appreciate the Riemmanian framework an overview of some of the neccesary tools upon which it is built is given as well. Variational calculus lies at the basis of geodesic tractography as the mathematical machinery allowing to deduce the geodesic equations from the length functional in the form of Euler-Lagrange equations. The solutions of the geodesic equations are the locally “shortest” paths between two fiducial points in the brain, when “distance” is measured in units of mean free path length by the (inverse) of the DTI tensor.

Linear perturbation theory is applied to the geodesic equations, taking into account perturbations in the metric field induced by DTI noise, and leading to an inhomogeneous geodesic deviation equation. Together with multiplicative calculus, notably product integrals, this allow us to express the first order behaviour of the perturbations along the tracks in closed-form. As long as the perturbations in the DTI tensor are sufficiently small (which is case dependent), the analytically computed perturbed tracks do not significantly differ from the explicitly computed tracks, which are geodesics with respect to the perturbed DTI metric.

This is confirmed by our experiments, both in the case of a single (unperturbed) geodesic as well as for bundles of geodesics. In the experiments we compared both the analytically computed perturbations as well as the simulated geodesics with respect to the perturbed metric. In practice, however, one may prefer the analytical approach, as it offers a strong computational advantage. By doing brute force simulations we need to solve non-linear differential equations to find the perturbed tracks, while in our analytical approach we employ a linear differential equation with a closed-form solution that can be efficiently computed.

We recall that geodesic completeness ensures (at least) one connecting track for any pair of sufficiently distant endpoints. In case of multiple connections external knowledge may be necessary to select the most plausible one from the biological point of view. For example, tractometric features can be used to classify and subsequently prune connections according to their diffusion-related characteristics (Colby et al. (2012); De Santis et al. (2014)). We will address this important point in future work.

Our perturbative analysis may be extended past the Riemannian paradigm by using Finsler geometry, which appears to be a natural extension for general DWI models beyond DTI. This extension removes the quadratic restriction in the metric and allows for a more accurate modelling of the DWI signal capturing additional features such as intra-voxel heterogeneity, e.g. due to partial volume effects or crossing fibers. The Finslerian approach yields stronger descriptive power, at the cost of a mathematically similar but more cumbersome representation. The Finslerian counterpart of the Riemannian geodesic deviation equation may likewise inform us on the effect of perturbations on general DWI data.

References

• 1

  AntonelliP. L.IngardenR. S.MatsumotoM. (1993). The Theory of Sprays and Finsler Spaces with Applications in Physics and Biology, Vol. 58 of Fundamental Theories of Physics. Dordrecht, Netherlands: Kluwer Academic Publishers.

  • Google Scholar

• 2

  AntonelliP. L.ZastawniakT. J. (1999). Fundamentals of Finslerian Diffusion with Applications, Vol. 101 of Fundamental Theories of Physics. Dordrecht, Netherlands: Kluwer Academic Publishers. 10.1007/978-94-011-4824-5

  • CrossRef
  • Google Scholar

• 3

  ArnoldV. I. (1989). Mathematical Methods of Classical Mechanics, Vol. 60 of Graduate Texts in Mathematics. second edn. New York: Springer-Verlag.

  • Google Scholar

• 4

  ArsignyV.FillardP.PennecX.AyacheN. (2006). Log-Euclidean Metrics for Fast and Simple Calculus on Diffusion Tensors. Magn. Reson. Med.56, 411–421. 10.1002/mrm.20965

  • CrossRef
  • Google Scholar

• 5

  AstolaL.FlorackL. (2011). Finsler Geometry on Higher Order Tensor fields and Applications to High Angular Resolution Diffusion Imaging. Int. J. Comput. Vis.92, 325–336. 10.1007/s11263-010-0377-z

  • CrossRef
  • Google Scholar

• 6

  AstolaL. J. (2010). “Multi-Scale Riemann-Finsler Geometry: Applications to Diffusion Tensor Imaging and High Angular Resolution Diffusion Imaging,” (Eindhoven, Netherlands: Eindhoven University of Technology, Department of Mathematics and Computer Science). Ph.D. thesis.

  • Google Scholar

• 7

  AstolaL.JalbaA.BalmashnovaE.FlorackL. (2011). Finsler Streamline Tracking with Single Tensor Orientation Distribution Function for High Angular Resolution Diffusion Imaging. J. Math. Imaging Vis.41, 170–181. 10.1007/s10851-011-0264-4

  • CrossRef
  • Google Scholar

• 8

  BaoD.ChernS.-S.ShenZ. (2000). An Introduction to Riemann-Finsler Geometry, Vol. 2000 of Graduate Texts in Mathematics. New York: Springer-Verlag.

  • Google Scholar

• 9

  BashirovA. E.KurpınarE. M.ÖzyapıcıA. (2008). Multiplicative Calculus and its Applications. J. Math. Anal. Appl.337, 36–48. 10.1016/j.jmaa.2007.03.081

  • CrossRef
  • Google Scholar

• 10

  BlochF. (1946). Nuclear Induction. Phys. Rev.70, 460–474. 10.1103/physrev.70.460

  • CrossRef
  • Google Scholar

• 11

  BrownR. W.HaackeE. M.ChengY.-C. N.ThompsonM. R.VenkatesanR. (2014). Magnetic Resonance Imaging: Physical Principles and Sequence Design. New York: John Wiley & Sons.

  • Google Scholar

• 12

  BurgethB.DidasS.FlorackL.WeickertJ. (2007). A Generic Approach to Diffusion Filtering of Matrix-fields. Computing81, 179–197. 10.1007/s00607-007-0248-9

  • CrossRef
  • Google Scholar

• 13

  CartanE. (1934). Les Espaces de Finsler. Paris: Hermann.

  • Google Scholar

• 14

  CartanE. (1963). Leçons sur la Géométrie des Espaces de Riemann. second edn. Paris: Gauthiers-Villars.

  • Google Scholar

• 15

  ColbyJ. B.SoderbergL.LebelC.DinovI. D.ThompsonP. M.SowellE. R. (2012). Along-tract Statistics Allow for Enhanced Tractography Analysis. NeuroImage59, 3227–3242. 10.1016/j.neuroimage.2011.11.004

  • CrossRef
  • Google Scholar

• 16

  CowanB. (1997). Nuclear Magnetic Resonance and Relaxation. Cambridge: Cambridge University Press.

  • Google Scholar

• 17

  De SantisS.DrakesmithM.BellsS.AssafY.JonesD. K. (2014). Why Diffusion Tensor MRI Does Well Only Some of the Time: Variance and Covariance of white Matter Tissue Microstructure Attributes in the Living Human Brain. Neuroimage89, 35–44. 10.1016/j.neuroimage.2013.12.003

  • CrossRef
  • Google Scholar

• 18

  Dela HaijeT. C. J. (2017). “Geometry in Diffusion Weighted MRI,” (Eindhoven, Netherlands: Eindhoven University of Technology). Ph.D. thesis.

  • Google Scholar

• 19

  Dela HaijeT.SavadjievP.FusterA.SchultzR. T.VermaR.FlorackL.et al (2019). Structural Connectivity Analysis Using Finsler Geometry. SIAM J. Imaging Sci.12, 551–575. 10.1137/18M1209428

  • CrossRef
  • Google Scholar

• 20

  LovelockD.RundH. (Editors) (1988). Tensors, Differential Forms, and Variational Principles (Mineola, N.Y.: Dover Publications, Inc.).

  • Google Scholar

• 21

  LovelockD.RundH. (Editors) (1989). Tensors, Differential Forms, and Variational Principles (Mineola, N.Y.: Dover Publications, Inc.).

  • Google Scholar

• 22

  do CarmoM. P. (1976). Differential Geometry of Curves and Surfaces. Mathematics: Theory & Applications. Englewood Cliffs, New Jersey: Prentice-Hall.

  • Google Scholar

• 23

  do CarmoM. P. (1993). Riemannian Geometry. Mathematics: Theory & Applications. second edn. Boston: Birkhäuser.

  • Google Scholar

• 24

  FillardP.PennecX.ArsignyV.AyacheN. (2007). Clinical DT-MRI Estimation, Smoothing, and Fiber Tracking with Log-Euclidean Metrics. IEEE Trans. Med. Imaging26, 1472–1482. 10.1109/tmi.2007.899173

  • CrossRef
  • Google Scholar

• 25

  FinslerP. (1918). “Ueber Kurven und Flächen in allgemeinen Räumen,” (Göttingen, Germany: University of Göttingen). Ph.D. thesis.

  • Google Scholar

• 26

  FlorackL.Dela HaijeT.FusterA. (2015). “Direction-controlled DTI Interpolation,” in Visualization and Processing of Tensors and Higher Order Descriptors for Multi-Valued Data. Mathematics and Visualization. Editors HotzI.SchultzT. (Berlin, Germany: Springer-Verlag), 149–162. 10.1007/978-3-319-15090-1_8

  • CrossRef
  • Google Scholar

• 27

  FlorackL.FusterA. (2014). “Riemann-Finsler Geometry for Diffusion Weighted Magnetic Resonance Imaging,” in Visualization and Processing of Tensors and Higher Order Descriptors for Multi-Valued Data. Mathematics and Visualization. Editors WestinC.-F.VilanovaA.BurgethB. (Berlin, Germany: Springer-Verlag), 189–208. 10.1007/978-3-642-54301-2_8

  • CrossRef
  • Google Scholar

• 28

  FlorackL. M. J.AstolaL. J. (2008). “A Multi-Resolution Framework for Diffusion Tensor Images,” in CVPR Workshop on Tensors in Image Processing and Computer Vision, Anchorage, Alaska, USA, June 24–26, 2008. Editors Aja FernándezS.de Luis GarciaR. (IEEE). Digital proceedings.

  • Google Scholar

• 29

  FlorackL.SengersR.MeestersS.SmoldersL.FusterA. (2021). “Riemann-DTI Geodesic Tractography Revisited,” in Anisotropy across Fields and Scales. Mathematics and Visualization. Editors ÖzarslanE.SchultzT.ZhangE.FusterA. (Berlin, Germany: Springer-Verlag), 225–243. 10.1007/978-3-030-56215-1_11

  • CrossRef
  • Google Scholar

• 30

  FlorackL.van AssenH. (2012). Multiplicative Calculus in Biomedical Image Analysis. J. Math. Imaging Vis.42, 64–75. 10.1007/s10851-011-0275-1

  • CrossRef
  • Google Scholar

• 31

  FusterA.Dela HaijeT.Tristán-VegaA.PlantingaB.WestinC.-F.FlorackL. (2016). Adjugate Diffusion Tensors for Geodesic Tractography in white Matter. J. Math. Imaging Vis.54, 1–14. 10.1007/s10851-015-0586-8

  • CrossRef
  • Google Scholar

• 32

  FusterA.PabstC. (2016). Finslerpp-waves. Phys. Rev. D94, 104072–1–104072–5. 10.1103/PhysRevD.94.104072

  • CrossRef
  • Google Scholar

• 33

  GantmacherF. R. (2001). The Theory of Matrices. Providence, Rhode Island: American Mathematical Society.

  • Google Scholar

• 34

  GibbonsG. W.GomisJ.PopeC. N. (2007). General Very Special Relativity Is Finsler Geometry. Phys. Rev. D76, 081701. 10.1103/PhysRevD.76.081701

  • CrossRef
  • Google Scholar

• 35

  GillR. D.JohansenS. (1990). A Survey of Product-Integration with a View toward Application in Survival Analysis. Ann. Stat.18, 1501–1555. 10.1214/aos/1176347865

  • CrossRef
  • Google Scholar

• 36

  GirelliF.LiberatiS.SindoniL. (2007). Planck-scale Modified Dispersion Relations and Finsler Geometry. Phys. Rev. D75. 10.1103/physrevd.75.064015

  • CrossRef
  • Google Scholar

• 37

  SaganH. (Editor) (1992). Calculus of Variations (New York: Dover Publications, Inc.).

  • Google Scholar

• 38

  HighamN. J. (2008). Functions of Matrices: Theory and Computation. Philadelphia, PA: SIAM.

  • Google Scholar

• 39

  JostJ. (2011). Riemannian Geometry and Geometric Analysis. Universitext. sixth edn. Berlin: Springer-Verlag.

  • Google Scholar

• 40

  KoenderinkJ. J. (1990). Solid Shape. Cambridge: MIT Press.

  • Google Scholar

• 41

  LengletC.DericheR.FaugerasO. (2004). “Inferring white Matter Geometry from Diffusion Tensor MRI: Application to Connectivity Mapping,” in Proceedings of the Eighth European Conference on Computer Vision (Prague, Czech Republic, May 2004). Editors PajdlaT.MatasJ. (Berlin: Springer-Verlag), 127–140. vol. 3021–3024 of Lecture Notes in Computer Science. 10.1007/978-3-540-24673-2_11

  • CrossRef
  • Google Scholar

• 42

  Maier-HeinK. H.NeherP. F.HoudeJ. C.CôtéM. A.GaryfallidisE.ZhongJ.et al (2017). The challenge of Mapping the Human Connectome Based on Diffusion Tractography. Nat. Commun.8, 1349. 10.1038/s41467-017-01285-x

  • CrossRef
  • Google Scholar

• 43

  MeestersS.LandersM.RuttenG.-J.FlorackL. (2021). “Automatic Tractography for Brain Tumor Surgery: Towards Clinical Application,” in Proceedings of the ISMRM & SMRT Annual Meeting & Exhibition (online, May 15–20, 2021).

  • Google Scholar

• 44

  MeestersS.RuttenG.-J.FusterA.FlorackL. (2019). “Automated Tractography of Four white Matter Fascicles in Support of Brain Tumor Surgery,” in 2019 OHBM Annual Meeting, Rome, Italy, June 6–13 2019 (Rome: Organization for Human Brain Mapping). Abstract nr. Th768.

  • Google Scholar

• 45

  MisnerC. W.ThorneK. S.WheelerJ. A. (1973). Gravitation. San Francisco: Freeman.

  • Google Scholar

• 46

  MoritaniT.EkholmS.WestessonP.-L. (2009). Diffusion-Weighted MR Imaging of the Brain. second edn. Berlin, Germany: Springer-Verlag. 10.1007/978-3-3540-78785-3

  • CrossRef
  • Google Scholar

• 47

  NeuenschwanderD. E. (2011). Emmy Noether’s Wonderful Theorem. Baltimore: The Johns Hopkins University Press.

  • Google Scholar

• 48

  O’DonnellL.HakerS.WestinC.-F. (2002). “New Approaches to Estimation of white Matter Connectivity in Diffusion Tensor MRI: Elliptic PDEs and Geodesics in a Tensor-Warped Space,” in Proceedings of the 5th International Conference on Medical Image Computing and Computer-Assisted Intervention—MICCAI 2002, Tokyo, Japan, September 25–28 2002. Editors DohiT.KikinisR. (Berlin, Germany: Springer-Verlag), 459–466. vol. 2488–2489 of Lecture Notes in Computer Science.

  • Google Scholar

• 49

  PennecX.FillardP.AyacheN. (2006). A Riemannian Framework for Tensor Computing. Int. J. Comput. Vis.66, 41–66. 10.1007/s11263-005-3222-z

  • CrossRef
  • Google Scholar

• 50

  PfeiferC.WohlfarthM. N. R. (2012). Finsler Geometric Extension of Einstein Gravity. Phys. Rev. D85. 10.1103/PhysRevD.85.064009

  • CrossRef
  • Google Scholar

• 51

  WeinstockR. (Editor) (1974). Calculus of Variations with Applications to Physics & Engineering (New York: Dover Publications, Inc.).

  • Google Scholar

• 52

  RundH. (1959). The Differential Geometry of Finsler Spaces. Berlin, Germany: Springer-Verlag.

  • Google Scholar

• 53

  RundH. (1973). The Hamilton-Jacobi Theory in the Calculus of Variations. Huntington, N.Y.: Robert E. Krieger Publishing Company.

  • Google Scholar

• 54

  RuttenG. J. M.KristoG.PigmansW.PelusoJ.VerheulH. B. (2014). Het gebruik van MR-tractografie in de dagelijkse neurochirurgische praktijk (the use of MR-tractography in daily neurosurgical practice). Tijdschrift voor Neurologie & Neurochirurgie115, 204–211. With English abstract.

  • Google Scholar

• 55

  SchillingK. G.NathV.HansenC.ParvathaneniP.BlaberJ.GaoY.et al (2019). Limits to Anatomical Accuracy of Diffusion Tractography Using Modern Approaches. NeuroImage185, 1–11. 10.1016/j.neuroimage.2018.10.029

  • CrossRef
  • Google Scholar

• 56

  SchillingK. G.PetitL.RheaultF.RemediosS.PierpaoliC.AndersonA. W.et al (2020). Brain Connections Derived from Diffusion MRI Tractography Can Be Highly Anatomically Accurate-If We Know where white Matter Pathways Start, where They End, and where They Do Not Go. Brain Struct. Funct.225, 2387–2402. 10.1007/s00429-020-02129-z

  • CrossRef
  • Google Scholar

• 57

  SengersR.FlorackL.FusterA. (2021). “Geodesic Tubes for Uncertainty Quantification in Diffusion MRI,”in Proceedings of the Twenty-Seventh International Conference on Information Processing in Medical Imaging–IPMI 2021, Bornholm, Denmark. Editors SommerS.FeragenA.SchnabelJ.NielsenM. (Berlin: Springer-Verlag), 279–290. vol. 12729 of Lecture Notes in Computer Science. 10.1007/978-3-030-78191-0_22

  • CrossRef
  • Google Scholar

• 58

  ShenY.-B.ShenZ. (2016). Introduction to Modern Finsler Geometry. Singapore: World Scientific.

  • Google Scholar

• 59

  SlavíkA. (2007). Product Integration, its History and Applications. Prague: Matfyzpress.

  • Google Scholar

• 60

  SpivakM. (1975). Differential Geometry, Vol. 1–5. Berkeley: Publish or Perish.

  • Google Scholar

• 61

  StejskalE. O.TannerJ. E. (1965). Spin Diffusion Measurements: Spin Echoes in the Presence of a Time‐Dependent Field Gradient. J. Chem. Phys.42, 288–292. 10.1063/1.1695690

  • CrossRef
  • Google Scholar

• 62

  StejskalE. O. (1965). Use of Spin Echoes in a Pulsed Magnetic‐Field Gradient to Study Anisotropic, Restricted Diffusion and Flow. J. Chem. Phys.43, 3597–3603. 10.1063/1.1696526

  • CrossRef
  • Google Scholar

• 63

  SzilasiJ.LovasR. L.KertészD. C. (2014). Connections, Sprays and Finsler Structures. Singapore: World Scientific.

  • Google Scholar

• 64

  TorreyH. C. (1956). Bloch Equations with Diffusion Terms. Phys. Rev.104, 563–565. 10.1103/physrev.104.563

  • CrossRef
  • Google Scholar

• 65

  WhitcherB.TuchD. S.WiscoJ. J.SorensenA. G.WangL. (2008). Using the Wild Bootstrap to Quantify Uncertainty in Diffusion Tensor Imaging. Hum. Brain Mapp.29, 346–362. 10.1002/hbm.20395

  • CrossRef
  • Google Scholar