Adjugate Diffusion Tensors for Geodesic Tractography in White Matter

One of the approaches in diffusion tensor imaging is to consider a Riemannian metric given by the inverse diffusion tensor. Such a metric is used for geodesic tractography and connectivity analysis in white matter. We propose a metric tensor given by the adjugate rather than the previously proposed inverse diffusion tensor. The adjugate metric can also be employed in the sharpening framework. Tractography experiments on synthetic and real brain diffusion data show improvement for high-curvature tracts and in the vicinity of isotropic diffusion regions relative to most results for inverse (sharpened) diffusion tensors, and especially on real data. In addition, adjugate tensors are shown to be more robust to noise.


Introduction
Geodesic tractography is one of the many existing approaches to perform tractography from diffusion images.The current state of the art tracking methods are based on high angular resolution diffusion imaging (HARDI), which can for example be described by multi-compartment [32] or higher order diffusion tensor models [37].HARDI tractography algorithms generally perform better than those arising from DTI, particularly in regions of complex fibre architecture such as crossings.Nevertheless, DTI is still widely used in a clinical research context since HARDI scanning protocols are by no means always available, and scanning and data processing times are considerably larger.It is therefore clinically useful to further improve existing methods and algorithms for DTI processing.
In the Riemannian framework for diffusion tensor imaging (DTI) [5], white matter is represented as a Riemannian manifold and neural fibres are conjectured to coincide with certain geodesic curves 1 (locally shortest paths in a non-Euclidean sense).In this way, the problem of tractography becomes one of finding geodesics.This is attractive from a practical point of view, as it obviates the need for ad hoc stopping and bending criteria necessary in traditional fibre-tracking algorithms.Another advantage with respect to other types of tracking algorithms is that geodesic tractography tends to be more robust to noise.Finally, it has the conceptual advantage that Riemannian geometry is a well understood and powerful theoretical machinery, facilitating mathematical modelling and algorithmics [2][3][4]6,12,[15][16][17][24][25][26]30,31,33,35].
However, there are problematic aspects to the existing formulation of the Riemannian paradigm, in which the metric is postulated to coincide with the inverse diffusion tensor [24,31].This idea is based on the transformation of anisotropic diffusion in Euclidean space to isotropic Brownian diffusion in a curved Riemannian space.As we show in this work this is however not achieved with such metric definition, despite claims in the diffusion MRI literature.
Another drawback of geodesic tractography based on the inverse diffusion tensor is the fact that geodesic curves tend to take shortcuts in the case of high-curvature tracts.A related problem is that the standard metric does not favour tracts through anisotropic diffusion regions over tracts through isotropic ones, making masking a necessary preprocessing step.
In this paper, we reconsider the relation between the DTI tensor and the Riemannian metric tensor, extending our preliminary work [17].We stipulate a novel Riemannian metric, given by the adjugate diffusion tensor (with and without sharpening), that does yield Brownian motion in the corresponding curved space.We investigate the practical implications of the proposed metric on geodesic tractography by performing experiments on clean and noisy synthetic data, and on real brain diffusion data.We compare our results with geodesic curves obtained from the inverse (sharpened) diffusion tensor, and with constrained spherical deconvolution tractography results.We evaluate our metric in relation to the problematic issues mentioned above, and show that these shortcomings are largely removed in our approach.
In recent work, both the inverse diffusion tensor and our metric have been extensively evaluated (using 40 subjects from the Human Connectome Project database) in combination with probabilistic shortest path tractography [36].It has been shown that our metric produces paths which agree most often with experts.

Related Work
The question of how to choose an appropriate metric in the context of DTI has already been posed in [20,25].Two types of modifications of the standard metric, given by the inverse diffusion tensor, have been proposed in the literature so far.In the first approach, "sharpening", diffusion tensors are raised to a certain power in order to increase their anisotropy where D and D sharp are the diffusion and sharpened tensor, and n > 1 is a constant integer.It has been shown that geodesics related to sharpened versions of the diffusion tensor follow the principal eigenvector directions more closely [19,38].A different sharpening strategy has been introduced by Descoteaux et al. [11], which relies on the transformation of diffusion tensors into so-called "fiber" tensors.They showed that geodesic tractography results improve by employing such deconvolution sharpening.Sharpening approaches seem to result in better tractography results.However, they decrease the robustness to noise.
Another downside of sharpening is that it makes use of parameters which have to be chosen in an ad hoc way: the power n in Eq. ( 1), and a factor controlling the fibre tensor sharpening in the case of deconvolution sharpening.
A second type of modification has been presented by Hao et al. [18,19].It is based on a conformal rescaling of the standard Riemannian metric given by the inverse diffusion tensor: Here the function α ≡ α(x) is chosen by requiring geodesic curves to follow more closely the diffusion tensor principal eigenvectors.White matter segmentation based on geodesic tractography with sharpened diffusion tensors and the modified metric (2) show similar improvement with respect to segmentation based on the standard metric.The metric modification that we propose in this work is somewhat similar to that by Hao et al., since both approaches rely on a conformal tensor rescaling.We address the relation between the two distinct metrics explicitly in Sect.3.4.

Preliminaries
First of all, our conventions and notation are summarized in Table 1.A (nondegenerate) diffusion generator L is a second order elliptic differential operator which in local coordinates {x i } on a manifold M has the form where a i j (x) and b i (x) are smooth functions and a i j (x) is a symmetric positive definite tensor.This is the general form of  a diffusion generator, and it captures many possible diffusionrelated processes.Several partial differential equations can be associated with L , for example, L u = 0 or ∂ t u − L u = 0, where u is the concentration of diffusing particles (or the distribution of heat in a given region over time).For example, ∂ t u −L u = 0 amounts to the diffusion-convection equation where a i j (x) represents a local diffusion tensor and b i (x) is the velocity of the medium.Note that the first term on the right-hand side corresponds to pure diffusion, while the second term relates to a convection process.We introduce the following Riemannian metric on M based on Eq. ( 3) where a i j denotes the inverse matrix of a i j .Then we can rewrite L as where B = B(a i j (x), b i (x)) is a smooth vector field and g is the Laplace-Beltrami operator w.r.t. the metric (5).This operator is given by where g = det g i j .It may be split as follows: The operator L is said to be an intrinsic Laplacian if B = 0.By definition, an intrinsic Laplacian generates Brownian motion on (M, g).For technical details, we refer to Liao [27] and Cohen de Lara [10].

Discrepancy
The standard anisotropic diffusion generator is given by where D i j is the diffusion tensor.Note that this is a special case of Eq. (3) with a i j = D i j and b i = ∂ j D i j .Therefore, a Riemannian metric g i j = D i j can be introduced, where D i j is the inverse diffusion tensor.The Laplace-Beltrami operator in this case reads with d = det D i j .Eq. ( 9) can thus be rewritten as: Comparing with Eq. ( 6) it is clear that The only scenario in which B = 0 is when the diffusion tensor determinant d is constant, which generally is not the case.We conclude that L 1 is not an intrinsic Laplacian and therefore the associated diffusion process is not a Brownian motion on (M, g).

Riemannian Framework Revisited
We propose to modify the Riemannian framework for DTI in such a way that the diffusion process associated with the diffusion generator is a Brownian motion on (M, g) for a certain Riemannian metric g, i.e.
The motivation behind this choice is to fully encode the anisotropic diffusion properties of the medium into the metric tensor.This is analogous to general relativity, where the properties of spacetime are completely described by the associated pseudo-Riemannian metric.Note that the standard metric given by the inverse diffusion tensor leads to a convection term in Eq. (11); the anisotropy information provided by this term will not be taken into account by this choice of metric.
Let us consider metrics which are conformally equivalent to g = D −1 , i.e.
where f ≡ f (x) is a positive scalar function.The corresponding Laplace-Beltrami operator, Eq. ( 7), is Here we have used the relation g = det gi j = f 3 d −1 .This expression is similar to the anisotropic generator given by Eq. ( 9), except for an overall scaling factor of 1/ f and the last term.The last term vanishes uniquely if f ∝ d, and so without loss of generality we set f = d so that gi j = d D i j , and By construction the generator L 2 is an intrinsic Laplacian.The diffusion process associated with L 2 is thus a Brownian motion on (M, g).Therefore, we postulate the following Riemannian metric in the context of DTI: Recall that, for a regular square matrix A, with adj(A) the adjugate matrix.Thus the proposed metric is the adjugate of the diffusion tensor, rather than the inverse (with det(A) = d, A i j = D i j in our case).Since the diffusion tensor is symmetric its adjugate equals the cofactor matrix.
Remark Diffusion generators in Eqs. ( 9), ( 16) are related by The following proposition concerning diffusion generators has been proven [27].Given a generator L which is an intrinsic Laplacian, and another generator L such that L = b L for some function b > 0, then L is an intrinsic Laplacian if and only if b is constant (in dimensions n = 2).We have shown that the generator L 2 is an intrinsic Laplacian by definition.From relation (19) and the above proposition it is clear that generator L 1 is an intrinsic Laplacian if and only if d is constant, a nongeneric case.

Relation to Previous Work
Note that the metric proposed in Hao et al. [18,19] is also of the form given by 2 Eq. ( 14).In their case, the function f (x) is determined by the equation where (grad g f ) i = ∂ j f D i j , V is the principal eigenvector field of the diffusion tensor and ∇ V V is the covariant derivative of V along itself: Here Γ i kl are the so-called Christoffel symbols: 2 Namely, (g Hao ) i j = f D i j = e α D i j .
In our case we have, by construction: Comparing Eqs. ( 20) and ( 23) it is clear that the conformal factors in Hao's and our metric satisfy rather different equations and the two metrics are therefore distinct.This is not surprising since, although both metrics are local rescalings of the inverse diffusion tensor, they arise from different considerations.In Hao et al. the local factor is chosen so that geodesic curves more closely follow the diffusion tensor principal eigenvectors.Our metric, on the other hand, relates anisotropic diffusion in Euclidean space to isotropic diffusion in the corresponding Riemannian space.
It is also important to note that our metric can be obtained from diffusion data in a straightforward way, since the (inverse) fitted diffusion tensor field needs just to be locally rescaled with the determinant.Hao's metric, however, is more difficult to deal with because its defining equation ( 20) is considerably more cumbersome.In addition, our metric has a simple and elegant interpretation as the adjugate diffusion tensor.

Geodesic Curves
Let us show how geodesic curves corresponding to the usual metric definition and our proposed metric relate to each other.Geodesic curves x(t) satisfy the geodesic equations ẍi + Γ i kl ẋk ẋl = 0, (24) where ẋi = dx i /dt and we consider unit speed parametrization g i j ẋi ẋ j = 1.It is straightforward to derive the relation between the Christoffel symbols corresponding to the standard metric g i j = D i j , and the adjugate metric gi j = d D i j : After some calculations the modified geodesic equations read: The term on the right-hand side represents the contribution of the conformal factor.Its second term merely affects the speed parametrization of g-geodesics and it can be parameterized away.The first term induces an effective force field which causes g-curves to bend in a different way compared to g-geodesics.It is clear that the modified geodesic equations reduce to the ones related to the standard metric g i j if the determinant d of the diffusion tensor is constant.This is consistent with the fact that a constant conformal factor does not modify geodesics.Equation (26) shows that modifying the metric in the way we propose does, in general, leads to geodesic curves which are intrinsically different than those obtained for the usual metric identification.

Method
We obtain geodesic curves from the inverse, inverse sharpened, and our newly proposed adjugate diffusion tensor.We also perform experiments for the adjugate sharpened tensor (see Appendix for the derivation).In particular, we use the normalized sharpened diffusion tensor in [13] where n > 1 is a constant.By using this normalization, the determinant of the sharpened tensor is still given by the determinant of the original diffusion tensor, i.e. tensor volume is preserved.We consider n = 2, 4, representative values employed in the literature [23,41].In order to find the shortest geodesic connecting a given target point to the seeding region, we use the optimization strategy proposed in [28], based on the fast sweeping algorithm introduced by Kao et al. [22].Our source code 3 implements the multi-threaded fast sweeping suggested by Zhao [43] and it is programmed in C++/ITK.Fast sweeping is based on dynamic programming, and it guarantees the convergence of iterative local computations to the globally shortest geodesic.Each iteration of the algorithm comprises an outer loop that considers all possible "sweeping directions", and an inner loop that traverses all voxels according to the given sweeping direction.The algorithm assigns to each visited voxel the minimum cost of reaching it from a set of neighbours following predefined spatial orientations, assuming that seeding points have zero cost.As local cost, we use the infinitesimal curve length function where again ẋi = dx i /dt, and g i j is the inverse, the inverse sharpened or the adjugate of the diffusion tensor.The minimum cost and the spatial direction chosen are stored at each voxel.The set of preferred orientations contains a vector field which is "back-traced" (integrated) from the targets to the seeding points, using an order two Runge-Kutta method, to retrieve the desired geodesics.
The number of iterations fast sweeping takes to converge depends on several factors, such as the number of neighbours 1.2 [21] considered in the inner loop, the total number of voxels to process and the curvature of the resulting geodesics, typically ranging from several tens to few hundred iterations.Recall that only a subset of geodesic curves corresponds to actual fibres; therefore we refer to geodesics either as "candidate fibres", or simply as tracts.

Results on Synthetic Data
We first demonstrate the method on a noiseless synthetic DTI data set.The fibres consist of rotated tensors with eigenvalues (λ 1 , λ 2 , λ 3 ) = (1.5, 0.5, 0.5) × 10 −3 , where rotation matrices are used to orient the tensor such that the principal eigenvector is parallel to the closest part of the centerline.Each voxel from which the distance to the centerline is smaller than 1.5 voxels is considered to be part of the fibre.The centerline is constructed by joining a half circle of radius 5 voxels, a horizontal straight line of length 5, a quarter circle of radius 8 and finally a straight vertical line of length 5.The surrounding tissue is comprised of isotropic tensors with eigenvalues (λ, λ, λ), where λ is taken to be 3λ 1 .
We regard realistic values for λ and λ 1 , based on a collection of experimental DTI measures in the literature.This case is representative of the interface between white matter in the corticospinal tract (CST) and the cerebrospinal fluid (CSF) in the ventricles.We consider axial diffusivity 4 values λ within the CST as principal eigenvalue λ 1 of the anisotropic fibre glyphs.On the other hand, we take the mean diffusivity in CSF as λ, since MD = (Tr D)/3 and in an isotropic region D = diag(λ, λ, λ).From Table 2 it can be seen that the relation λ = 3λ 1 is a good estimate in this scenario.Noiseless experiments are shown in Fig. 1.
Next we perform experiments on noisy synthetic data, obtained by adding Rician noise.We consider two different noise levels, σ = 0.15 and σ = 0.3, with σ the standard deviation of the underlying Gaussian distribution.We compare geodesic curves from the inverse and adjugate diffusion tensor, in three different cases: no sharpening, sharpening factor 4 By definition λ = λ 1 .
123  , 2, 3).We consider separately the U-fibre and the longer, less curved upward tract.In the first case, we seed from the lower middle voxel to the upper middle one; in the second case we use the latter as seed point.We visualize the obtained tracts using vIST/e [42].Hereby, we depict geodesics from inverse diffusion tensors in black, while we use magenta for those from adjugate diffusion tensors.Noisy data experiments are shown in Figs. 2 and 3.
For both the U-and the longer fibres, we see that geodesics from the inverse (unsharpened) diffusion tensor take a shortcut through the isotropic background, completely failing to describe the fibres in all situations (Figs.1a, 3a).When a sharpening factor n = 2 is used this is still the case for the inverse sharpened tensor, except for the longer fibre and σ = 0.15 (Fig. 2b).The latter geodesic, however, degrades for σ = 0.3, taking again a shortcut through the isotropic background (Fig. 3b).When the higher sharpening factor n = 4 is used geodesics from the inverse sharpened tensor nicely follow the synthetic tracts, although a slight degradation is observed for σ = 0.3 (Fig. 3c).This effect for sharpened metrics had been shown in [38], for a sharpening power n = 2 (and a slightly different normalization factor in Eq. ( 27)), and in [19] for n = 3.In the latter, however, the isotropic background has been masked out.Note that we obtain good sharpening results in the case n = 4 but not for n = 2, in contrast to [38].
On the other hand, geodesics from adjugate tensors, with or without sharpening, follow the synthetic fibres rather well in all scenarios and without taking shortcuts through the isotropic background.As in the n = 4 inverse tensor case, a slight degradation is observed for the n = 4 adjugate sharpened tensor σ = 0.3 (Fig. 3c).Comparing geodesics from the adjugate unsharpened tensor (Figs.1a, 3a) and the n = 4 inverse or adjugate sharpened tensor (Figs.1c, 3c), we observe that the sharpened ones follow the fibres more closely in the noiseless case.However, in the σ = 0.3 case these degrade by taking a shortcut of about one voxel, while the adjugate unsharpened ones remain almost unchanged.Sharpening thus appears to decrease the robustness to noise.This shortcoming of sharpened tensors had already been pointed out in [19].

Results on Real Data
We consider a diffusion MRI data set with 64 gradient directions and a b-value of 3000 s/mm 2 ; the dimensions are 128 × 128 × 60 and the voxel size is 1.75 × 1.75 × 2 mm 3 , corresponding to a patient with a tumour located next to the ventricles.We have segmented the cerebrospinal fluid inside the ventricles, together with the tumour.We seed from the cerebral peduncles to a number of target points in the motor cortex, and in the cingulum region.We visualize the obtained tracts using 3D Slicer [34].
We also present higher order model tractography experiments on this data set.In particular, we use the MRtrix package [40] to perform constrained spherical deconvolutionbased [39] deterministic and probabilistic fibre tracking.The algorithm either follows the peaks of the FOD, the fibre orientation distribution (deterministic) or uses orientations sampled from the FOD at each step (probabilistic).Optimal parameter values are selected by visual inspection and are given in Table 3.The cerebral peduncles are again used as a seed region and only fibres reaching the motor cortex are selected.Obtained results are visualized with vIST/e [42].
In Figs. 4, 5 and 6, we show candidate fibres reaching the trunk and foot motor area of the cortex (upward bundle) and the lip area (bundle bending to the left), which ought to correspond to the corticospinal and corticobulbar tracts.Results above the ventricles are clearly consistent with the left and right cingulum.In Fig. 4, we show tractography results for metrics given by the inverse and adjugate diffusion tensor, and the outcome for inverse sharpened diffusion tensors is given in Fig. 5. Results obtained with our approach, Fig. 4b, seem to better resemble the anatomy of the stipulated white matter bundles.Additionally, the curvature of the candidate fibres is smoother and the bundles are more coherent.A particularly interesting result is the fact that our candidate fibres circumvent the ventricles, known to be void of fibres, while most of the ones obtained with other approaches go through them.Note that for inverse sharpened tensors, Fig. 5, less bundles cross the CSF than in the original diffusion tensor case, Fig. 4a.Still, the problem is not completely overcome, as is the case in our approach, Fig. 4b.These results are consistent with our synthetic data experiments in the case of sharpening power 2, Fig. 5a, but worst than expected for sharpening power 4, Fig. 5b.This is likely to be explained by the presence of noise in real data, since sharpening decreases the robustness to noise as we have seen in the synthetic experiments.
In Fig. 4, we also see that our tracts do not go through the tumour.This is consistent with our findings concerning the CSF since diffusion in tumours is usually also isotropic.Our results may reflect real fibres being pushed aside by a tumour, or white matter integrity inside the tumour having been destroyed.In contrast to the ventricles case, however, fibres might be found within a tumour and therefore we cannot draw any decisive conclusions about the validity of our results in this sense.In Fig. 6, we show results for a metric given by the adjugate sharpened tensor.Note that the outcome improves drastically compared to that from inverse sharpening in Fig. 5, to the extent that none of the tracts cross isotropic diffusion regions in this case.Results appear to be very similar to those from the adjugate diffusion tensor in Fig. 4b.
In Fig. 7, we show CSD-based deterministic and probabilistic tractography results for the corticospinal tract, together with those obtained from geodesic tractography based on the adjugate diffusion tensor.In all cases, fibres circumvent the ventricles.However, a noticeable difference is that CSD fibres do not reach the anterior part of the (right) motor cortex.This could be achieved by tuning the parameters (step size 0.5 mm, minimum radius of curvature ≥5 mm, FOD amplitude cutoff ≤0.01 mm), but results in unrealistically straight fibres crossing the ventricles and jumping from one hemisphere to the other.

Intuition Behind Results
The rather different behaviour of the metrics given by the inverse and adjugate diffusion tensor can be intuitively explained by the following argument, where for simplicity we regard the case with no sharpening.Consider two neighbouring voxels with a typical diffusion tensor D = diag(λ, λ, λ) in an isotropic region and D = diag(λ 1 , λ 2 , λ 3 ), with λ 1 > λ 2 = λ 3 , in a vertically oriented fibre bundle (see Fig. 8).
To fix ideas, we take as example the case λ = λ 1 .Using the classical metric, the Riemannian cost (Eq.28) of (travelling along) an infinitesimal vertical line element scales with 1/λ 1 .In this case, the classical metric clearly assigns the same cost to the isotropic and anisotropic line elements, since directional diffusivities (lengths of the vertical lines) are equal in both cases.Using the newly proposed metric, however, the Riemannian cost of the same line elements scales with λ 2 λ 3 , the area of the orthogonal cross section indicated by the shaded equatorial planes.This is clearly smaller in the anisotropic case, leading to a smaller cost.Therefore, the adjugate tensor metric favours the anisotropic tract over the isotropic one.
This argument clearly holds as well when isotropic diffusivities are larger than the anisotropic ones (λ > λ 1 ) as in the presented synthetic experiments and (often) in real diffusion data, since the isotropic cost becomes even larger.In this case, the isotropic region is preferred by the classical metric since its Riemannian cost, related to 1/λ, is smaller than the anisotropic one, 1/λ 1 .In fact, the classical metric is only able to avoid isotropic regions when λ < λ 1 .Our metric favours anisotropic regions up to the limit λ ≤ λ 2 , λ 3 .In this scenario, the area of the orthogonal cross section in the isotropic case becomes equal to or smaller than the anisotropic one, and so does the Riemannian cost.However, such scenarios seemingly take place in real data only in the case of complex architecture, where DTI fails to describe the underlying diffusion profiles in any case.

Conclusion
We have proposed a new Riemannian metric in the context of diffusion tensor imaging, namely, the adjugate of the diffusion tensor.In the sharpening framework, this translates into a metric given by the adjugate of the (normalized) sharpened tensor.This is derived in a rigorous way from the relation between anisotropic diffusion in Euclidean space and isotropic diffusion in the corresponding curved space.Our metric represents solely diffusion, the process which is encoded in the diffusion MRI signal, in contrast to the standard DTI metric which leads to additional convection in the curved space.
We show results of geodesic tracking on synthetic and real brain diffusion data, based on our new metric and other established ways to extract the metric from the diffusion tensor, and also compare our results to constrained spherical deconvolution (CSD).Note that isotropic diffusion regions are not masked out in a preprocessing step.Moreover, we regard a realistic synthetic scenario based on experimental DTI measures in the literature.
Results on synthetic data show that geodesics from the inverse diffusion tensor fail to describe the fibres for both clean and noisy data when either a low sharpening factor (n = 2) or no sharpening is used.They do succeed for a higher sharpening factor (n = 4), although a slight degradation is observed in the presence of noise.Geodesics from adjugate tensors, with and without sharpening, follow the synthetic fibres rather well in all scenarios and without taking shortcuts through the isotropic background.Again, a slight degradation with noise is observed for a high sharpening factor (n = 4).Comparing geodesics from the adjugate unsharpened tensor and the n = 4 (inverse or adjugate) sharpened tensor, we observe that the sharpened ones follow the synthetic fibres more closely in the noiseless case.However, we see that sharpening decreases the robustness to noise, as has been pointed out in the literature.We observe that adjugate tensors are less sensitive to noise than inverse tensors, in particular for no sharpening or low sharpening factor.
In real brain data, tracts obtained with our adjugate metric, with and without sharpening, avoid isotropic diffusion regions such as ventricles.Experiments show that this is definitely not the case for the standard DTI metric, and only sometimes for metrics given by the inverse sharpened diffusion tensor.These results are consistent with the synthetic experiments outcome.The presence of noise in real data seemingly negatively affects inverse sharpened tensors, while this does not appear to be the case for adjugate tensors.We therefore conclude that the adjugate framework leads to better results, also in the case of sharpening.The positive performance of our adjugate approach on real diffusion data agrees with the recent literature [36].
Finally, we obtain comparable results for the corticospinal tract from the adjugate tensor method and (deterministic and probabilistic) CSD tractography.The only noticeable difference being that CSD tracts do not reach the anterior part of the motor cortex but one cannot draw strong conclusions from this since ground truth is not available.
In terms of practicalities, in our approach, there are no free parameters such as the sharpening power or those related to the fibre orientation distribution in CSD.These parameters have to be chosen in an ad hoc way and a globally satisfactory setting might not exist, which can be a disadvantage.On the other hand, such free parameters do offer some flexibility to model the diffusion data.
In future work, we will evaluate our adjugate method for geodesic tractography of subcortical U-fibres.Based on the performed experiments, we would expect to recover such fibres relatively well, which is not the case for classical diffusion tensor tractography methods [9].In addition, it has been shown that DTI geodesic tractography results improve by using a multivalued geodesic algorithm [38].This aspect could also be evaluated in our case by employing such an algorithm instead of fast sweeping.It would also be interesting to compare our method to the deconvolution sharpening in Descoteaux et al. [11], and to the different Riemannian approach in Hao et al. [19].
The method we propose is based on DTI, which is well known to suffer shortcomings in regions of complex fibre architecture.However, higher order diffusion models may benefit from our approach as well, provided one can define a proper metrical distance.For example, the framework proposed in [14] stipulates a Finsler metric for geodesic tractography in HARDI, which can in principle be adapted in a similar way to our modification of the Riemannian metric in DTI.
Here (D sharp ) i j denotes the inverse sharpened tensor.Note that ĝ = det ĝi j = d 2 .The corresponding Laplace-Beltrami operator, Eq. ( 7), is given by which may be split as The diffusion generator is, by construction, an intrinsic Laplacian.The diffusion process associated with L 3 is thus a Brownian motion on (M, ĝ).
Remark The adjugate sharpened diffusion tensor relates properly to Brownian motion if and only if the sharpened tensor is normalized such that its determinant coincides with that of the original diffusion tensor, as in Eq. ( 29).Adjugate versions of sharpened tensors such as those proposed in [23,38,41] will not lead to Brownian motion. 123

Fig. 1
Fig. 1 Results on synthetic noiseless data, for metrics given by the inverse diffusion tensor (black) and the adjugate diffusion tensor (magenta) in the cases a unsharpened diffusion tensor, b sharpening factor n = 2 and c sharpening factor n = 4. Background and fibres are RGB colour coded based on the direction of the diffusion tensor main eigenvector.Inverse tensor geodesics fail to describe fibres except for n = 4. Adjugate tensor geodesics follow the fibres well, and (higher) sharpening improves results further (Color figure online)

Fig. 2 Fig. 3 Table 3
Fig. 2 Results on synthetic data with Rician noise of σ = 0.15 (standard deviation of the underlying Gaussian distribution), for metrics given by the inverse diffusion tensor (black) and the adjugate diffusion tensor (magenta) in the cases a unsharpened diffusion tensor, b sharpening factor n = 2 and c sharpening factor n = 4. Colour coding as in Fig. 1.Again, inverse tensor geodesics fail to describe fibres except for n = 4, although the longer fibre n = 2 tracking does improve w.r.t. the noiseless case.Adjugate tensor geodesics follow the fibres well, and (higher) sharpening improves results further.The Ufibre adjugate n = 4 tracking degrades slightly w.r.t. the noiseless case (Color figure online)

Fig. 4
Fig. 4 Candidate fibres possibly corresponding to corticobulbar and corticospinal tracts (brown and blue, respectively), and cingulum (red), in an anterior view.No candidate fibres shown in-between since we do not consider target points in that part of the cortex.A tumour is located next to the ventricles on the left-hand side.Results for metric given by a inverse diffusion tensor and b adjugate diffusion tensor.Candidate fibres going through the ventricles or the tumour are indicated by yellow and white arrows, respectively.Bundles obtained with our approach, in b, avoid both the CSF in the ventricles and the tumour (Color figure online)

Fig. 5
Fig.5As in Fig.4, but now showing results for metric given by a inverse sharpened diffusion tensor d 1/3 D −2 and b inverse sharpened diffusion tensor d D −4 .Note that results from sharpened tensors improve compared to those without sharpening in Fig.4a(i.e. less tracts cross isotropic diffusion regions), but the problem is not completely overcome as in our approach (Color figure online)

Fig. 6
Fig. 6 As in Fig. 5, but now showing results for metric given by a adjugate sharpened diffusion tensor d 4/3 D −2 and b adjugate sharpened diffusion tensor d 2 D −4 .Note that results from adjugate sharpened tensors improve drastically compared to those from inverse sharpening in Fig. 5 (i.e.none of the tracts cross isotropic diffusion regions).The outcome is very similar to that from the adjugate diffusion tensor, Fig. 4b (Color figure online)

Fig. 7
Fig. 7 CSD-based tractography results for the corticospinal tract (yellow), together with those obtained by geodesic tractography from the adjugate diffusion tensor (blue).a Deterministic CSD and b Probabilistic CSD.In all cases, fibres circumvent the ventricles.CSD fibres do not reach the anterior part of the (right) motor cortex (Color figure online)

Fig. 8
Fig. 8 Graphical sketch of the quadratic forms corresponding to a typical diffusion tensor D in an isotropic region (left) and in a vertically oriented fibre bundle (right); the vertical axis corresponds to λ = λ 1 (Color figure online)

Table 1
Conventions and notation used through the paper.Tensor index notation is employed, with latin indices running from 1 to 3

Table 2
[29] and axial diffusivity values in CSF and white matter (WM) in the corticospinal tract, expressed in units of 10 −3 mm 2 /s.Literature references are indicated in the table.WM-CST diffusivity value in[29]corresponds to the posterior limb of the internal capsule (PLIC), which contains a.o.corticospinal fibres