Improving Fiber Alignment in HARDI by Combining Contextual PDE Flow with Constrained Spherical Deconvolution

We propose two strategies to improve the quality of tractography results computed from diffusion weighted magnetic resonance imaging (DW-MRI) data. Both methods are based on the same PDE framework, defined in the coupled space of positions and orientations, associated with a stochastic process describing the enhancement of elongated structures while preserving crossing structures. In the first method we use the enhancement PDE for contextual regularization of a fiber orientation distribution (FOD) that is obtained on individual voxels from high angular resolution diffusion imaging (HARDI) data via constrained spherical deconvolution (CSD). Thereby we improve the FOD as input for subsequent tractography. Secondly, we introduce the fiber to bundle coherence (FBC), a measure for quantification of fiber alignment. The FBC is computed from a tractography result using the same PDE framework and provides a criterion for removing the spurious fibers. We validate the proposed combination of CSD and enhancement on phantom data and on human data, acquired with different scanning protocols. On the phantom data we find that PDE enhancements improve both local metrics and global metrics of tractography results, compared to CSD without enhancements. On the human data we show that the enhancements allow for a better reconstruction of crossing fiber bundles and they reduce the variability of the tractography output with respect to the acquisition parameters. Finally, we show that both the enhancement of the FODs and the use of the FBC measure on the tractography improve the stability with respect to different stochastic realizations of probabilistic tractography. This is shown in a clinical application: the reconstruction of the optic radiation for epilepsy surgery planning.


Introduction
Diffusion weighted magnetic resonance imaging (DW-MRI) is a non-invasive technique for the characterization of biological tissue microstructure [1]. In brain white matter, water molecules diffuse predominantly along axonal fibers. This results in an observable macroscopic orientation dependence in the DW signal, that is measured by scanning the tissue in multiple orientations and gradient strengths. To model the angular anistropy of the diffusion profile, diffusion tensor imaging (DTI) [2] is widely used, but this has the limitation that only a single fiber direction can be estimated per voxel [3]. It is estimated in [4] that more complex fiber configurations occur in approximately 90% of the white matter voxels. To overcome this, high angular resolution diffusion imaging (HARDI) techniques are used, that can describe more complex (crossing) fiber configurations. An overview of HARDI techniques can be found in [5]. Here we use the method of constrained spherical deconvolution (CSD) [6], that from the initial diffusion data constructs a fiber orientation distribution (FOD), which models the distribution of fibers along different directions.
Tractography methods are often used in the DW-MRI pipeline to provide insight in the structural connectivity of the white matter bundles. Independently of the model used for interpreting the DW-MRI data, noise originating from the scanner, acquisition artifacts and partial volume effects [7] are likely to result in spurious (aberrant) fibers in the tractography output. To improve the data on which the tractography is performed, different regularization methods can be used. Methods exist that apply filtering for the reduction of noise directly on the DW-MRI data [8][9][10], other methods aim to regularize the DTI tensor fields [11][12][13][14][15]. On HARDI data the regularization can be performed on individual voxels [16][17][18] or in combination with the local spatial information [19][20][21][22][23].
We introduce two new strategies based on the same underlying principle to improve fiber alignment in tractography results, in order to have more reliable information on the structural connectivity of brain. First we perform contextual regularization to the FOD obtained with CSD, see Fig 1A, and secondly we introduce a fiber to bundle coherence (FBC) measure that can be applied to any fiber bundle to classify and remove spurious fibers, see Fig 1B. Both approaches are based on a partial differential equation (PDE) framework introduced in [24][25][26][27], where the Fokker-Planck equation of a stochastic process for enhancement of elongated structures is considered. These type of PDE-based enhancement methods have been widely used for the processing of 2D-images. In this framework, images are represented in the extended space of positions and orientations via a stable invertible orientation score [28], that associates to every location an orientation distribution of the local image features (lines and contours). Then, the stochastic processes for contour completion [29][30][31][32][33][34] and contour enhancement [35][36][37] (see Fig 2A) on this extended space R 2 ⋊ S 1 induce crossing preserving completion and enhancement of lines [28].
The DW-MRI data that we use is naturally defined on the coupled space R 3 ⋊ S 2 of 3D positions and orientations. As in the 2D case, crossing preserving enhancement of line structures is required, for which we use the 3D extension of the 2D stochastic process for contour enhancement, introduced in [25]. The linear PDE corresponding to this stochastic process can be solved by convolution of the initial condition with the kernel of the PDE. This kernel is also a function on the position-orientation space and can be seen as a transition distribution from the origin (in position and orientation) to neighboring elements. From the stochastic point of view, the kernels arise as limits of the accumulation of infinitely many sample paths drawn from the stochastic process, illustrated in Fig 2A. For mathematical details of the underlying stochastic processes of the PDEs, see [27 §10.1]. The general idea needed for this article is sketched in Fig 2. In Fig 2B and 2C we show the contour enhancement kernel using glyph visualization on a grid, each glyph being a polar (red, 2D) or spherical (blue, 3D) graph plot where in every orientation the (spherical) radius is proportional to the value of the kernel. This type of visualization is used throughout the paper for functions defined on the space of positions and orientations. The proposed pipeline of the paper. CSD is used to estimate an FOD from DW-MRI data. The FOD is enhanced (A) with PDE techniques. Then a deterministic or probabilistic tractography is applied to the (enhanced) FOD (probabilistic shown here, with coloring indicating the fiber direction). In the lower right figure, we applied our coherence quantification method (B), based on the same PDE framework, which shows that blue fibers are well aligned (high Fiber to Bundle Coherence (FBC)) and yellow fibers are spurious (low FBC). The green arrows indicate the steps in which the contextual PDEs are used.
Recently, many authors [23,25,27,34,[38][39][40][41][42] demonstrated the advantages of contextual processing of DW-MRI data. The general rationale behind contextual processing is to include alignment of local orientations and their surroundings (i.e. the context) on the coupled space of positions and orientations. For this alignment of local orientations, roto-translations are needed, which imposes a non-Euclidean structure in the PDE-based processing as we explain in Section 2.2. More details on the embedding of R 3 ⋊ S 2 in the roto-translation group SE(3) can be found in [25]. This demonstrates how either the completion or enhancement PDEs can be used to extrapolate DTI information to increase the angular resolution and resolve some fiber crossings. This idea was shown to be promising in clinical experiments [38,39], but in some cases extreme parameters had to be set to obtain clear maxima at crossings (where DTI data is inadequate). Therefore in this paper we introduce and test the combination of CSD with contextual enhancements. The method proposed in [34] uses an advection-diffusion equation (that we called contour completion above) to improve HARDI data to obtain connectivity measures. In our work we rely on a purely diffusive process, contour enhancement, which in contrast to contour completion does not suffer from singularities [27] and is less sensitive to small perturbations of the initial conditions. This property makes the enhancement process more suited to be combined with the sharp angular distributions produced by CSD. As the methods mentioned above still result in broad angular distributions, they need to be combined with some sharpening method. To this end, a geometric morphological sharpening based on erosions was presented in [23,27,42]. Another related method presented in [40,41] is the so-called fiber continuity model in which purely spatial regularization is considered in combination with spherical deconvolution as alternative to the non-negativity constraint in the classical CSD [43]. In Section 2.2 we demonstrate the importance of including also an angular regularization term.

Contributions
The first contribution of this article is to study the combination of the widely used CSD method with a regularization induced by the enhancement PDE acting on the FOD. Since the FOD obtained with CSD consists of sharp angular profiles, it is well-suited as an initial condition for the enhancement PDE, that typically has a smoothing effect on the orientation distributions. The contextual regularization method reduces non-aligned crossings in the FOD, allowing for a better alignment of fibers when tracking is applied on the enhanced FOD. We show that this method is therefore useful to reduce the number of false positive fibers, but mainly to find more true positives in the tractography output. Although in this paper we compare to the classical CSD method, the PDE enhancements can also be applied to extensions of this method [44][45][46][47][48].
The second contribution of this article is to introduce the fiber to bundle coherence (FBC) measure. The motivation for this measure is that, especially probabilistic, tracking methods typically produce spurious fibers that should be removed from the tractography. In contrast to the first approach, this method serves as a post-processing tool. For the computation of the FBC we regard the fiber bundle as a set of oriented points, by considering for every fiber point also the local tangent to the fiber. We construct a density using the enhancement PDE with an initial condition that is a sum of superposed δ-distributions at every oriented point in the bundle. The construction of such a density from tracks relates to track density imaging [49] and track orientation density imaging [50], though here the use of the contour enhancement kernels, Fig 2, allows to use a sparse set of fiber tracks. The FBC, a measure for spuriousness of fibers, is computed by efficient integration of this fiber-based density. Fibers that are most spurious according to the FBC can be removed from the tractography, resulting in a better aligned fiber bundle. Complementary to the first method, this FBC measure has the purpose to remove false positives in a tractography.

Structure of the Article
Section 2 covers theory of the individual parts of the pipeline as outlined in Fig 1, consisting of CSD, PDE enhancements, tractography and coherence quantification in Sections 2.1-2.4, respectively. In Section 3 we provide extensive validation of the combination of CSD and PDE enhancements and the FBC, using three experiments: 1. First we use the Tractometer evaluation system [51,52] on the ISBI 2013 HARDI reconstruction challenge dataset [53], a digital phantom with known ground truth, to demonstrate how contour enhancement improves both the local FOD reconstruction and the global connectivity of fiber bundles compared to CSD, see Section 3.1.
2. In Section 3.2 we show on a human DW-MRI dataset, containing different crossing bundles, that CSD combined with enhancements yields an FOD that is more robust with respect to the b-value and the number of gradient directions used in the acquisition. Furthermore, we make a comparison with earlier work involving erosions and nonlinear diffusion of FODs directly applied to a DTI-model [23,27], that was based on the same data. We show that with our method the glyphs are sharper at the locations where bundles cross.
3. Finally in Section 3.3, we show an experiment with clinical data in which we reconstruct the optic radiation (OR) to determine the position of the tip of the Meyer's loop, that is of interest in epilepsy surgery planning [23,[54][55][56][57]. Accurate estimation of this position is difficult due to the presence of spurious fibers in the reconstruction of the OR. We show that both the FOD enhancement and the FBC measure, see Fig 1, and in particular the combination of the two allow for a more stable determination of the tip of the Meyer's loop. Here 'more stable' means less variation with respect to stochastic realizations in the probabilistic tractography results.
Conclusions and a discussion can be found in Section 4.

Methods
In this paper it is assumed that we have HARDI data as input, from which we derive an FOD U that models the orientation of fibers in each voxel, i.e. U : R 3 × S 2 ! R + . For this we use CSD [6], concisely described in Section 2.1, as it gives sharp angular profiles and is able to distinguish multiple fiber directions within a voxel. Then we use the enhancement PDE for diffusion of the FOD U, coupling spatial and angular information. The combination of CSD and such enhancement is a powerful method to obtain an enhanced FOD in which the coherence inherent in the data is included, providing a more coherent input for the tractography. The enhancement technique is explained in Section 2.2.
We use the MRtrix algorithm [6] for both deterministic and probabilistic tractography to estimate the structural connectivity in the brain. In the deterministic tractography, fiber tracks are obtained by integrating a directional field, given an initial position and direction. The directional field is given by the locally maximal orientations in the glyphs. In contrast to deterministic tractography, the probabilistic tractography method of MRtrix samples the orientations from the entire FOD and does not use just the maxima. More difficult paths can be reconstructed than with deterministic tracking, but typically also many spurious fibers are produced due to the probabilistic sampling. Both the deterministic and the probabilistic method are explained in more detail in Section 2.3.
In Section 2.4 we introduce our new technique to quantify the coherence of fibers with respect to all the fibers in a bundle, based on the same PDE theory as employed for the contextual enhancement in Section 2.2. We explain how the kernel of the enhancement PDE is used to construct a tractography-based density, how the FBC is computed and how this measure is able to classify spurious fibers in a tractography.

A Brief Review of CSD
In CSD it is assumed that at each voxel position y the measured signal S y : S 2 ! R can be represented by a spherical convolution of the FOD f y : S 2 ! R with a response function K : S 2 ! R, that is estimated from the data [58]. Since the spherical deconvolution to determine the FOD is ill-posed, a non-negativity constraint is included as in [43,59]. Then, given the signal S y (n) for a sample of orientations n 2 S 2 , the solution of CSD is found by iteratively solving the minimization problem: for i = 1,. . .,i max , with i max the maximum number of iterations. Here K 2 L 2 (S 2 ) is aligned with and symmetric around the z-axis, the convolution Ã S 2 is the usual S 2 spherical convolution [60], dσ(n) is the Jacobian of the surface measure in orientation n and λ is a parameter to influence the trade-off between the data driven term and regularization term. The linear operator L h : L 2 (S 2 ) ! L 2 (S 2 ) in the regularization term gives the non-negativity constraint and is defined by: where H is the Heaviside function and τ h is a threshold equal to a fixed factor τ times the mean of h. The initial function f 0 y for the iteration is computed by taking only the data driven term of Eq (1). The iteration stops when successive iterations yield the same result, typically after 5 to 10 iterations [59]. Throughout the paper, we call U the FOD obtained by In practice CSD is performed using spherical harmonics with a maximal spherical harmonic order of 8 (l max = 8) as discussed in [61].
Improvements to the original CSD exist to modify and improve the response function, either by recursive calibration or auto-calibration [44,45], by using multiple acquisition shells [46] or by including anatomical data [47,48]. The latter two methods aim to reduce the partial volume effects, where CSD is likely to produce spurious fiber orientations. These partial volume effects can occur when in a voxel multiple tissues or multiple bundles with different orientation are present. Here we use the classical CSD as it is the basic technique available in several neuroimaging packages. However, we stress that our method is not restricted to this type of CSD. In any case, our method aims to reduce non-aligned crossings in the FOD, also the ones induced by partial volume effects, as we will show in several experiments in this paper. Further improvement of the methodology can be expected when including recently extended and more elaborate CSD techniques [45,46,48], but this is left for future work.

Contour Enhancement (Step A)
To improve alignment of neighboring glyphs of the FOD U, recall the glyph field visualization in Figs 1 and 2C, we apply contextual enhancements. Before we specify the PDE we consider for this enhancement, we first need to express the notion of alignment in mathematical terms. To this end, let us consider Fig 3, where it is shown that the notion of alignment cannot be supported by a decoupled, flat Cartesian product R 3 × S 2 with the combined Euclidean distance. It is clear that the green bar at (y 1 ,n 1 ) is better aligned with the gray bar at (y 0 ,n 0 ) than the orange bar at (y 2 ,n 2 ), even though the distances in the space R 3 × S 2 are equal, i.e. ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d 2 . This means that in order to appropriately describe the concept of alignment, we must consider more than just the amount of spatial displacement and the amount of change in orientation. Coupling these two types of motion (via rigid body motions) is a solution to this problem [25]. The coupling follows very naturally by expressing the motion of an oriented particle (y,n) in terms of a moving frame of reference determined by its orientation. That is, spatial movement along the orientation n should be much cheaper than spatial movement in the plane orthogonal to n. This creates a natural anisotropy for spatial movement. For angular motion we need isotropy. This extra structure can be obtained by embedding the space of positions and orientations in the rigid body motion group. This means that an element (y,n) 2 R 3 × S 2 is identified with the rigid body motion (y, R n ), where R n is any rotation matrix such that R n e z = n, with e z 2 S 2 pointing to the north pole. We denote this space of coupled positions and orientations by R 3 ⋊ S 2 , so we have The group R 3 ⋊ S 2 is equipped with the following (non-commutative) group product: This product moves oriented elements in a shift-twist fashion, rather than by a rotation followed by an independent translation. Due to this shift-twist group product in Eq (5), we automatically express motion of oriented particles in terms of a moving frame in R 3 ⋊ S 2 , which makes this space well-suited for the application of our contextual enhancements. Nevertheless, in the remainder of this article this space can still be regarded as the Cartesian 5D space R 3 × S 2 , where we secured the coupling of positions and orientations via our specific choice of differential operators and diffusions that are applied.
To improve alignment of FOD glyphs, we use a particular diffusion process called contour enhancement that uses both spatial and angular diffusion in the extended space of positions and orientations [25]. Given a structure (think of a fiber bundle) in this space, see Fig 4, we apply spatial diffusion only in the direction of the structure, not in the spatial plane perpendicular to it. Angular diffusion is applied in the plane tangent to S 2 at the point n. This diffusion process enhances elongated structures, while preserving crossing structures, and is given by a Fokker-Planck type of system, a linear diffusion equation on R 3 ⋊ S 2 . For t ! 0, y 2 R 3 , n 2 S 2 this system can be expressed as: @ t Wðy; n; tÞ ¼ ðD 33 ðn Á r y Þ 2 þ D 44 D S 2 ÞWðy; n; tÞ; Wðy; n; 0Þ ¼ Uðy; nÞ: [62]. The symbol r y denotes the gradient with respect to the spatial variables and Δ S 2 is the Laplace-Beltrami operator on the sphere. Parameters D 33 > 0 and D 44 > 0 are related to the amount of spatial and angular diffusion, respectivly. Parameter t ! 0 is the diffusion time of the contour enhancement process. It can be seen as a Brownian motion process, recall Fig 2A, where particles are allowed to spatially move back and forth in the direction they are heading, or change their direction, but are not allowed to step aside (comparable to the movement of a car). We refer to the solution of Eq (6) as the enhanced FOD. It can be obtained via a finite difference scheme [63], or via a convolution with a kernel p t : R 3 ⋊ S 2 ! R + : A basic approximation to the exact Green's function of the contour enhancement PDE is known [25] and can be written as the product of Green's functions p R 2 ⋊S 1 t in the following way: The R 2 ⋊ S 1 kernels are given by To avoid numerical errors, we use the estimate y=2 tan ðy=2Þ % cos ðy=2Þ 1Àðy 2 =24Þ for jyj < p 10 . This approximation is easy to use and allows for efficient implementation [64].
From the approximation kernel in Eq (8) it can be seen that problems could occur when either D 33 = 0 or D 44 = 0. To this end, a necessary and sufficient condition for the existence of a smooth solution kernel for the evolution process in Eq (6) is given by the Hörmander requirement [65]. This condition applies to more general situations than the one here, see e.g. [25], but for the specific case of contour enhancement the requirement is satisfied iff D 33 ,D 44 > 0. Setting D 44 = 0 would result in a singular non-smooth kernel, which has numerical disadvantages. More importantly, apart from this theoretical issue the need for both spatial and angular diffusion can also be argued from a practical point of view, as is illustrated in  Improving Fiber Alignment in HARDI via Contextual PDE Flow with CSD diffusion the peak is redirected and the glyphs lie better aligned with the fiber bundle. Hence D 44 > 0 is needed to ensure the crucial interaction between different orientations. Finally we recall the relation between Tikhonov regularization and diffusion, see e.g. [25], which allows us to connect diffusion with D 44 = 0 with the fiber continuity model in [40,41]. This model does not suffer from the inconvenience of considering only spatial regularization, as they represent the FOD in a truncated spherical harmonic basis. When the enhancements are used in combination with probabilistic tractography, we first apply a standard sharpening deconvolution transform to the FOD as described in [18], to maintain the sharpness of the FOD.

Tractography
As the next step in the pipeline we use the MRtrix tractography algorithm [6], as implemented in http://www.brain.org.au/software/index.html#mrtrix, version 0.2.12. It allows us to perform deterministic and probabilistic fiber tracking on spherical harmonic representations of the (enhanced) FOD. To have a fair comparison between trackings on the FOD and the enhanced FOD, we use the parameter settings as explained next.
• In the deterministic tracking of MRtrix, seed points are randomly selected from a seed region. The initial direction is sampled randomly and every next step follows the direction of the most aligned FOD maximum. If this maximum is below a threshold value, the fiber terminates. This threshold (cutoff) is set to 10% of the maximal angular response of the FOD.
There is no constraint on the maximal curvature of the fibers. To prevent that fibers have an initial direction that is not aligned with the fiber bundle, we force the initial direction to be approximately in the direction of the maximal FOD peak, by setting the initial cutoff to 0.9. The step size is set to 1/10th of the voxel size as is suggested in [6]. Tracks proceed in both directions from the seed point and terminate either when they hit the boundary of the volume or mask (if applicable), or due to the threshold stopping criterion.
• In the probabilistic case, starting from the seed region, every next step follows a direction randomly sampled from the FOD. Here we set the minimal radius of curvature to 1 mm, the default value in the MRtrix algorithm. Optionally, a target region of interest is used to select only those fibers that cross this region.
The methods proposed in this paper are not tied to the type of tractography used. In the Experiment and Results section we validate our methods combined with both deterministic and probabilistic tractography, on synthetic and real data. In the phantom experiment described in Section 3.1, our preference for deterministic tractography is due to the fact that deterministic tractography was reported in [52] to perform better in the considered dataset than probabilistic tractography. Since only a seed region is specified and no additional information is used to filter out spurious fibers, it is difficult to assess the quality of the probabilistic tractography results with respect to the used measures. The same holds for the experiment on real data in Section 3.2. On the other hand, on the optic radiation application of Section 3.3 where both a seed and end region are specified, a probabilistic method is required for a more thorough delineation of the fiber bundle. Actually, in many clinical applications in which the fiber configuration can be even more complex than the phantom data we have considered, it might be preferable to use probabilistic tractography, because more information can be provided to the neurologists. Then they can decide, in combination with all of the other clinical information at their disposal, which aspects of the tractography result are to be trusted.
Streamlines from a probabilistic tractography result that are anatomically implausible can be removed with scoring methods [23,56] or by imposing anatomical constraints. Even when using these methods, the filtered tractography output can still contain fibers that deviate from the fiber bundle and are likely to be spurious. In the next section, we propose a coherence measure for fibers in a fiber bundle in order to classify these spurious fibers.

Coherence Quantification of Fiber Bundles (Step B)
In this section we introduce our second contribution of the paper, a fiber to bundle coherence (FBC) measure to quantify the coherence of each fiber with respect to all other fibers in the bundle, recall Fig 1B. A spurious fiber, as schematically shown in Fig 6, is isolated from or poorly aligned with the bulk of the tracks and is therefore unlikely to represent the underlying brain structure. Fibers with low coherence, i.e. a low FBC, can then be classified as spurious.
To classify a fiber as spurious, we first construct a density by regarding each fiber as a superposition of δ-distributions in R 3 ⋊ S 2 and convolving this distribution with the kernel in Eq (8). This density is independent of the underlying data and is based purely on the collection of fibers Γ. Integration of this density along a part of length α of a fiber gives a local measure for the coherence of that part.
Next we explain the mathematical techniques that support the idea in Fig 6. We denote the fibers from a tractography output by y i (s) 2 R 3 , 1 i N, 0 s l i , with s the arc length parameter, l i the total length of fiber i and N the number of fibers. Now let n i ðsÞ ¼ _ y i ðsÞ be the tangent of the fiber, so that γ i (s) = (y i (s),n i (s)) forms a curve (fiber) in R 3 × S 2 . By construction, n i (s) points in the forward direction of the fiber. Since in DW-MRI data antipodal orientations are identified, we also consider g i ðsÞ ¼ ðy i ðsÞ; Àn i ðsÞÞ. The complete fiber bundle is defined as A discrete formulation of a fiber i with N i points is given by: This way there are N tot ¼ 2 P N i¼1 N i elements in Γ. Now we regard every point g j i as a δ-distribution in R 3 ⋊ S 2 centered around ðy j i ; n j i Þ. A density for the entire bundle is then constructed as follows: with index j running over points within a fiber, i running over all fibers and σ taking care of including forward and backward orientations. We use the same evolution process as in Eq (6) in which F = F Γ now serves as initial condition, to create a diffused density (y,n) 7 ! W F (y,n,t): @ t W F ðy; n; tÞ ¼ ðD 33 ðn Á r y Þ 2 þ D 44 D S 2 Þ W F ðy; n; tÞ; W F ðy; n; 0Þ ¼ Fðy; nÞ: with the shift-twist convolution as given in Eq (7). This is illustrated in Fig 6 in the 2D case.
We can now define the FBC for fiber γ i with respect to the bundle Γ as the integral of this density along the fiber: This results in a global property of the fiber, but spurious fibers often only locally deviate from the bundle as in Fig 6. To this end, we compute for each fiber the minimum of such integrals along the fiber over intervals of length α: The parameter α defines the scale over which spuriousness of fibers can be detected and is much smaller than the average fiber length. Our primary interest is not the FBC α value itself, but rather how it compares to the average coherence of fibers in the bundle, so finally we define the relative fiber to bundle coherence (RFBC) as: Here AFBC(Γ) is the average fiber to bundle coherence indicating the overall coherence of the N fibers in the bundle Γ, defined as To summarize, the RFBC(γ i ,Γ) of a fiber γ i in a bundle Γ is a measure for how well aligned the least aligned part of γ i is, compared to the average coherence of the total bundle. In practice, we evaluate the convolution in Eq (14) only in the fiber points. We compute the LFBCðg k i ; GÞ, the diffused density in the oriented point g k i ¼ ðy k i ; n k i Þ, recall the notation in Eq (11), as follows: where R n l j is any rotation matrix such that R n l j e z ¼ n l j , index q sums the contributions along a fiber, index j runs over all the fibers and σ as before. The FBC α can then be computed as where a,α 2 N in this discrete case, so the LFBC is summed along short intervals of the fiber. Likewise, the AFBC can be computed as We apply this method in Section 3.3 for quantifying the coherence of tractography results of the optic radiation and classifying the spurious fibers.

Experiments and Results
In this section we extensively test the performance of our CSD enhancement method • We use the HARDI Reconstruction Challenge dataset [65], which is artificial data with known ground truth, to quantitatively evaluate the CSD enhancement method (A) on deterministic tractography in Section 3.1.
• In Section 3.2 we show on DW-MRI human brain data that the enhancement (A) have a positive effect on deterministic tractography, for different acquisition protocols of the data. Furthermore, on this DW-MRI dataset and on the phantom dataset we compare our method to previous work [27], where a DTI-based FOD is used in combination with nonlinear PDE flow.
• In the third and last experiment, we reconstruct the optic radiation in human clinical data, see Section 3.3. We include an extensive evaluation of our methods, the enhancement of the FOD (A) and the use of the FBC to classify and remove spurious fibers (B), and the combination of both methods. We show that the reproducibility of the probabilistic tractography has increased, resulting in a more stable localization of the tip of the Meyer's loop.
For all datasets Mathematica [67] was used to perform the contour enhancement algorithm and the CSD, which in practice produces the same results as the MRtrix CSD implementation when the same deconvolution kernel is used. MRtrix software [6] was used to perform fiber tractography. The coherence quantification was implemented in C++. In Section 3.1 we make use of the Tractometer [51,52]

HARDI Reconstruction Challenge
The following experiment is performed on a digital phantom dataset that was designed for the ISBI 2013 Reconstruction Challenge [66,69]. It is used in combination with the Tractometer [51,52], as a benchmark to compare different reconstruction and tracking methods. The phantom is inspired by the Numerical Fiber Generator [70] and the code to reproduce it is freely available as part of the Python package Phantomas (http://www.emmanuelcaruyer.com/ phantomas.php). This synthetic dataset is of size 50 × 50 × 50 voxels with a resolution of 1 × 1 × 1 mm 3 . It consists of 27 simulated white matter bundles, designed to resemble challenging branching, kissing and crossing structures at angles between 30 and 90 degrees, with various curvature and bundle diameters ranging from 2 mm to 6 mm. An image indicating the ground truth fiber configuration is shown in the centre of Fig 7. The idea behind the signal simulation is that every voxel is subdivided into multiple subvoxels, each one with its own attenuation profile. The final signal arrives from integrating the contribution of all the sub-voxels. Then, it is possible to combine multiple compartment types in every voxel with added Rician noise. This allows for modelling complex configurations as well as taking into account partial volume effects. While the Numerical Fiber Generator uses a tensor-like model to simulate the signal in the sub-voxels, Phantomas uses a CHARMEDbased model [71]. The CHARMED model based on the Söderman-Jönsson cylinder model [72] captures well the non-Gaussian behaviour of the diffusion signal for large b-values. The main reason why we selected the ISBI phantom is that it is linked with the Tractometer that allows for performing quantitative evaluations of the tractography results, using global metrics as demonstrated in the subsequent experiments.
For the experiments presented in this section we used 64 uniformly distributed gradient directions using a b-value of 3000 s/mm 2 with different signal to noise ratios (SNRs). We use spherical harmonics in CSD with maximal order 8, resulting in 45 estimated coefficients on each position. We then enhance the resulting FOD functions using our contour enhancement algorithm with varying parameters. From the evolutions described in Eq (6) we see by a basic rescaling argument that it is sufficient to vary t and the ratio D 33 /D 44 . The larger this ratio, the more preference the spatial diffusion gets over the angular diffusion, resulting in elongated kernels (visualized by thin glyphs). A smaller ratio D 33 /D 44 is better suited in regions where the curvature of bundles is higher (visualized by thicker glyphs). The higher the diffusion time t, the more context is taken into account. When t is too large, fiber bundles with high curvature can be damaged or false positives could be created. Taking this into consideration, we choose our parameters as follows: we fix spatial diffusivity parameter D 33 = 1.0, we take the angular diffusivity parameter D 44  Tractography results for the entire dataset are shown in Fig 7. We can recognize the positive effect of the enhancements on deterministic tractography: we see less dropouts, better aligned fibers and better continuation of fibers at crossings. An extensive quantification of the performance of our method is done at the voxel level using the FODs and at the macroscopic level using tractography in Sections 3.1.1 and 3.1.2, respectively. Both sections support the results summarized in Fig 8. 3.3.1 Local Metrics. We compare reconstructed FODs locally with the ground truth using only the orientation of the peaks. Let M be the set of voxels v in the white matter mask, then we denote the ground truth number of peaks in a voxel v by N v and the orientations corresponding to the peaks by n v i;true , i = 1,. . .,N v . Maxima of the constructed FOD are found by evaluating the FODs on a 60 th order icosahedron tessellation with 18606 antipodally symmetric points, giving an angular resolution of less than 1 degree. Maxima are taken into account only if it exceeds a threshold of 0.1, the same value we use as threshold in the tractography. Let O v est be the set of peak orientations in voxel v estimated from the FOD. The average angular error in degrees can then be computed by: In the top row of Fig 8 we show the effects contour enhancement for different ratios of D 33 and D 44 upon variation of the diffusion time. The results are given for substantially low SNR levels 10, 6 and 4 and 2. These SNRs are computed w.r.t. the non-DW image. Specifically, if the b = 0 intensity is 1 then the standard deviation of the Rician noise distribution is 1/SNR. In all cases a clear improvement is found compared to CSD without enhancements and the more noise, the more the angular error is decreased. Higher diffusion times give better results and around t = 5 the angular error is almost stable. It can also be seen that the combination of CSD with enhancements at lower SNRs gives lower angular errors than just CSD for the higher SNRs.   Improving Fiber Alignment in HARDI via Contextual PDE Flow with CSD MRtrix tractography is used as described in Section 2.3, with seeds randomly selected in the white matter mask. The tracks have a minimum length of 10 mm and new seed points are chosen until 10000 streamlines are selected. For every FOD, the tractography is repeated five times with the same settings, to average out the variability in the tracking algorithm output. We then use the Tractometer [52] to perform a fiber tracking analysis based on the ground truth and the five results are averaged. The Tractometer outputs values for various metrics, from which we use the Valid Connections (VC), Invalid Connections (IC) and No Connections (NC). They indicate the percentage of tracks that correctly connect, incorrectly connect or do not connect gray matter areas in the dataset, respectively. We also use the Average Bundle Coverage (ABC), the percentage of voxels in a bundle that is crossed by a valid streamline, averaged over all bundles. We combine the (VC), (IC) and (NC) in two metrics introduced in [73]: • Connection to Seed Ratio (CSR), which represents the probability that a generated fiber actually connects two gray matter areas, computed as 100%−NC.
• Valid Connection to Connection Ratio (VCCR), the probability that a connecting fiber is correct, computed as VC/(VC+IC).
The results for the ABC, CSR and VCCR with the same enhancement parameters and SNRs as for the local metric are given in Fig 8. Similar remarks hold for the global metrics as for the angular error. For all three metrics and all SNRs the enhancements lead to an improvement compared to CSD, the only exception being the ABC for SNR = 10 and D 44 = 0.02. Furthermore, as the SNR decreases, the larger diffusion times are beneficial and the more significant the improvement is. The best results are obtained for D 44 = 0.005. We expect that truncation of the spherical harmonics already introduces some angular smoothing of the FODs on this artificial dataset, explaining the small effect of D 44 in the experiments. Furthermore, we see that the diffusion time t truly acts as a regularization parameter, resulting in a robustness for the metrics with respect to the SNRs: the higher the diffusion time, the smaller the differences in the metrics between the different SNRs.
Seeding from the white matter voxels can lead to an over-representation of the number of fibers in longer fiber bundles with respect to the shorter bundles [74]. The longer bundles thereby have a larger contribution to the global metrics than the shorter bundles, which could lead to an overestimation of the fiber bundles. As proposed in [75], we compared the global metrics when seeding from the gray/white matter interface for CSD and one specific set of enhancement parameters. The global metrics for that seeding strategy were slightly lower for CSD and comparable when including enhancements. For the sake of comparing our enhancement method with CSD, we therefore believe it is fair to use seeding from the white matter mask.
The convincing improvement in the global metrics is supported by Fig 9, that shows a selection of the fiber bundles in the dataset. It can be seen that after enhancements, there are more valid connections in the green bundle and less wrong exits in the red bundle, leading to a higher (VCCR) and a better bundle coverage. The glyphs in the top row show that the enhancements improve alignment of glyphs, especially at the boundary of the fiber bundles, where the original CSD result tends to suffer from partial volume effects.

Evaluation and Comparison on DW-MRI Data
In this experiment we consider a DW-MRI dataset of a part of a human brain, previously used in [27]. The study was approved by the local ethical commitee of Maastricht University, and informed written consent was obtained from the subject. Although the dataset consists of only 10 axial slices, the corpus callosum, corona radiata and superior longitudinal fasciculus are (partly) present in the data. We show that the combination of CSD and enhancement is well-suited for different combinations of the b-value and the number of gradient directions used in the acquisition. Furthermore, we make a qualitative comparison with the DTI-based method of [27] on this dataset and conclude with a brief quantitative comparison with this method on the dataset of 3.1.
3.2.1 Robustness with Respect to the Acquisition Parameters. The acquisition was performed on a 3T Siemens Allegra scanner, with FOV 208x208mm and voxel size 2x2x2mm. During the data acquisition, a brain region consisting of 10 axial slices was scanned with the following combinations of b-values and N o , the number of orientations: b = 1000 s/mm 2 with N o = 49, b = 1000 s/mm 2 with N o = 121 and b = 4000 s/mm 2 with N o = 49. The SNR in the non-DW image was estimated to be approximately 3 using the approach of [76] as implemented in Dipy [77]. We use again CSD with spherical harmonics up to order 8. The higher Improving Fiber Alignment in HARDI via Contextual PDE Flow with CSD b-value is obtained by using a stronger gradient pulse, making the acquisition more sensitive to detail in the tissue structure, but also inducing a lower SNR. Increasing the number of gradient directions gives a better angular resolution. We use deterministic tractography, with three seed regions manually selected in the middle of the corpus callosum, corona radiata and superior longitudinal fasciculus.
In the right column of Fig 10 we show that after enhancements, the FOD allows for a more coherent reconstruction of the three bundles. Especially in the region where the three bundles come together, it can be seen that the fibers have a better propagation through the crossings. Moreover, the FODs after enhancements are very similar to each other, visible in the glyph visualization, leading to three tractography results supporting similar fiber bundles. This is an improvement with respect to CSD without enhancement, shown in the left column of Fig 10. There we find more noisy FODs with more variation between the different protocols. This is also reflected in the tractography results, that contain more spurious fibers than after the enhancements.
We conclude, just like in the first experiment on the phantom data, that applying enhancements induces more robust tractography also on real DW-MRI data, in this case in the sense that it is less sensitive to the acquisition parameters b and N o .

Comparison with a DTI-based FOD.
In the next experiment we compare the performance of our combination of CSD with enhancements with the method in [27] which proposed to combine DTI with non-linear PDE-based enhancement obtained from successively applying erosions and diffusions. Let us briefly describe this method, for details we refer to [27], and an implementation of the PDE enhancements can be found in the HARDI package for Mathematica available at (http://bmia.bmt.tue.nl/people/RDuits/HARDIAlgorithms.zip). First an FOD on positions and orientations that we call U DTI was constructed via a transformation of the tensor field D fitted to the data [2], according to the following definition [27,78]: This FOD is then sharpened with PDE erosions, a type of morphological enhancement adapted from [15], on R 3 ⋊ S 2 and regularized with nonlinear diffusions to find crossing structures from DTI. Previously in [27], the same dataset as in Fig 10 for acquisition parameters b = 1000 s/mm 2 and N o = 49 was processed. Here we compare the FOD obtained with CSD, that we call U CSD here, with U DTI in the top and bottom figures, respectively, of Fig 11. Unlike DTI, which is limited by the Gaussian assumption of the diffusion profile, CSD can estimate multiple fiber orientations within a voxel. Furthermore, we see that the large glyphs in the Centrum Semiovale in the bottom figure are not apparent in U CSD . Applying (linear) enhancements, as explained in Section 2.2, to U CSD gives the second figure, and the approach in [27] using erosions/(nonlinear) enhancements applied to U DTI gives the second figure from below. It can be seen that also the enhanced DTI glyphs supports multiple fiber directions within voxels via extrapolation [27,38], but at the cost of high regularization. Another noticeable difference is the fact that the glyphs in the CSD case are slimmer and crossings are more clearly defined. Whether two separate maxima are visible at a crossing is less dependent on the diffusion parameters in the PDE diffusion.
Besides the visual comparison of the FOD glyphs, we provide deterministic tractography results for both procedures in the middle of Fig 11. It can be observed that both methods produce reasonable results, although the one obtained from the enhanced DTI dataset seems oversmoothed and outliers (indicated with the yellow arrow) can occur. This is due to the extreme diffusion parameters needed to perform the FOD extrapolation. We find that visually the  Improving Fiber Alignment in HARDI via Contextual PDE Flow with CSD combination of CSD and linear enhancements yields better tractography than DTI combined with erosions and nonlinear enhancements.
To provide a more quantitative and complete comparison of DTI, DTI and nonlinear enhancements, CSD and CSD with linear enhancements, we also include results of the experiment in Section 3.1 for the DTI methods, see Table 1. We heuristically determined good parameter settings for the nonlinear enhancement of DTI: erosions [27] with D 11 = 0.5, D 44 = 0.2, t = 2.0 and diffusion [27 Eq. (55)] with D 11 = 0.2, D 33 = 1.0, D 44 = 0.02 and t = 3. In Table 1 is shown that applying enhancements for contextual regularization of the FOD is beneficial for both DTI and CSD. The lower the SNR, the more evident the improvements become. Furthermore, we see that in terms of the local metric, the angular error θ of the peak orientations, the DTI methods can compete with the CSD based methods. However, the global metrics are significantly higher for CSD based methods. The quantitative results on the phantom data in Table 1 are in line with the qualitative comparison on real data in Fig 11.

Improved Reconstruction of the Optic Radiation
The optic radiation (OR) is a white matter fiber bundle connecting the primary visual cortex and the lateral geniculate nucleus (LGN), see Fig 12. The most anterior part of the OR is called the Meyer's loop (ML), of which the exact location is of interest for treatment of temporal lobe epilepsy [23,54,56,57]. During neurosurgery, a part of the temporal lobe is resected. To ensure that the OR remains intact to prevent visual field defect, it is crucial to know the distance from the tip of the Meyer's loop to the Temporal Pole (ML-TP) [55], which shows large interpatient variability [79].
We use DW-MRI scans of four subjects, performed on a 3.0T Philips Achieva MR scanner, with b = 1000 s/mm 2 , N o = 32 and a spatial resolution of 2x2x2 mm. All subjects gave written informed consent; the study was approved by the Medical Ethics Committee of Maastricht University Medical Center (N 43386.068). The data is acquired from healthy volunteers, and ground-truth ML-TP distance is not known. Therefore accuracy of this measure of our methods cannot be checked, instead we focus on consistency and reproducibility. We apply CSD to the data to construct the FOD, with spherical harmonics up to order 6 requiring the estimation of 28 coefficients (as 32 directions are insufficient to estimate the 45 coefficients when a spherical harmonic order 8 is used, when not using super-resolution as in [59]). We seed from the LGN and include all fibers that reach the primary visual cortex. Both regions of interest are Table 1. For two SNR values, the results are shown for the DTI method described in Section 3.2, with or without nonlinear enhancements. We compare with CSD and a specific instance of enhanced CSD with parameters D 33 = 1, D 44 = 0.01,t = 2. For local metric θ lower is better, for the other metrics higher is better. In boldface are the best results for the DTI and CSD methods. selected manually on a T1-weighted image. We use probabilistic fiber tracking as described in Section 2.3. We demonstrate the effect of the enhancement of CSD and the use of the FBC measure in Sections 3.3.1 and 3.3.2, respectively, in this relevant clinical setting. A quantitative comparison of the four methods CSD (O), CSD + enhancement (A), CSD + FBC (B) and CSD + enhancement + FBC (A+B) is provided in Section 3.3.3. We show that the enhancement and/or the removal of spurious fibers, but in particular the combination of both methods, allows for a more stable computation of the ML-TP distance than the original tractography result.
3.3.1 Effect of the Enhancement of CSD on Tractography of the OR. In this section, we apply the PDE enhancement (step A) to the CSD FOD as before, with parameter settings D 33 = 1, D 44 = 0.01 and t = 2. After the enhancement we apply the sharpening deconvolution transform [18] and probabilistic tractography with 10000 streamlines. We compare the results of the tractography on the subjects both before and after the enhancement in Fig 13. We see that the tracking on enhanced data generally shows less spurious fibers, and has a better pronounced tip of the Meyer's loop. However, the optic radiation is a highly curved structure, where the advantage of the enhancement of elongated structures cannot be fully exploited. To further reduce the spurious fibers, we explore our other approach in the next section.
3.3.2 Effect of the FBC measure on Tractography of the OR. In this section, we apply probabilistic tractography on subject 1, with 20000 streamlines and including state of the art data scoring as in [23] (only relying on the data term, i.e. λ = 0 in [23]), see Fig 14. The kernel parameters for the coherence quantification (step B) are set to D 33 = 1, D 44 = 0.04 and t = 1.4 for the convolution [27]. Let Γ be the set of the 1000 most anterior fibers in a tractography of the OR, that roughly form the Meyer's loop. We compute the LFBC and subsequently the RFBC for all the fibers in Γ. Fig 12. A reconstruction of the optic radiation and its positioning in the brain. The left figure shows how the OR is positioned in the brain, the close-up on the right shows how the OR wraps around the ventricular system. The probabilistic tractography outputs many spurious fibers. The tip of the Meyer's loop, indicated by the orange sphere, is localized on a spurious fiber and is therefore very dependent on the realization of the tractography. As a result, the distance from the Meyer's loop to the Temporal pole (ML-TP) that is used in temporal lobe resection surgery, shows a high variation among different tractography outcomes. RFBCðg; GÞ, the RFBC corresponding to the "central" fiber, in the sense that it is most coherent with the fiber bundle. We define the filtered set Γ as G ≔fg 2 Gj RFBCðg; GÞ ! g; This means the parameter acts as a threshold parameter and can be set such that fibers with a high spuriousness are removed. The fiber point in Γ that is closest to the temporal pole defines   spheres) is estimated at different locations. When we set the threshold = 0.1 max , removing in these cases between 6% and 8% of the most spurious fibers, we obtain the results as shown in the bottom row of Fig 14. It can be seen that the resulting fiber bundles are very similar to each other, demonstrating less variation in the localization of the tip.

Quantitative Comparisons on Four Subjects.
To support our claims of the two previous sections, we test the effect of our methods on the stability of the ML-TP distance under different stochastic realizations. Here we perform probabilistic tractography with 10000 fibers ten times with the same settings, for each of the four subjects and each of the four methods (CSD, CSD + enh, CSD + FBC and CSD + enh + FBC). The FBC measure is computed from the 1000 most anterior fibers as in the previous experiment and the threshold is set to = 0.05 max . We compare the mean ML-TP distance and sample standard deviation determined from the tracking results of each of the methods. The results are summarized in the boxplots in Fig 15. The figure strongly supports the application of the enhancements methods. For subjects 1-3 the ML-TP distance shows much less variation when including the FBC. For all subjects also (CSD + enh) gives more stable results than just CSD. Moreover, in all cases the combination (CSD + enh + FBC) outperforms CSD and for all but subject 1 the combined method (CSD + enh + FBC) also gives better results than the enhancement or FBC individually. It should be remarked that higher up the graph indicates a larger resection if used for pre-surgical evaluation, which is not necessarily positive. However, we prefer to have a stable and reproducible method that can be used with a safety margin, then a method that is more conservative, but shows large variations.

Conclusions and Discussion
We have proposed two new tools to improve alignment of fibers in tractography results: (A) the combination of CSD with contextual PDE enhancements and (B) a fiber to bundle coherence measure to classify spurious fibers. Both approaches rely on the same contextual processing via PDEs on the space of coupled positions and orientations. We validate our methodology with a variety of experiments on synthetic and human data.
In the first experiment we consider a digital phantom [66] that simulates DW-MRI data of a challenging configuration of multiple neural-like fiber bundles for different noise levels, see Fig  7. The combination of CSD with enhancements and subsequent deterministic tracking was extensively tested for varying enhancement parameters, see Fig 8. The enhanced FOD peaks were compared with the ground truth fiber orientations, showing for all SNRs that the maxima of the enhanced FOD coincide better with the ground truth peaks than without application of enhancement. Also, this improvement is particularly high for very low SNR values. To quantitatively evaluate the impact of the enhancement on the tractographies we used the Tractometer evaluation system [52]. The results, shown in Fig 8 confirm the benefit, for all the metrics considered, of including the enhancement. Also an improved stability of the metrics with respect to different enhancement parameters is observed. Furthermore, we found that data with a lower SNR requires more regularization, obtained by choosing a higher diffusion time t in the enhancement. These quantitative evaluations of local and global metrics are supported by the qualitative results in Figs 7 and 9, where we saw that after enhancement fibers are better aligned and propagate better through crossings.
The second experiment is performed on human data of a representative area of the brain with crossing fiber bundles. We evaluate our combination of CSD and enhancement for three different (single-shell) acquisition protocols, corresponding to different b-values and number of gradient directions. We observed, see Fig 10, that whereas tractography on CSD without enhancement showed notable differences between the three acquisition protocols, tractography after our enhancement lead to a qualitatively similar reconstruction in all cases. This implies that the application of enhancement in the processing pipeline makes the tractography results less dependent on the scanning protocol used.
We use the same dataset and the phantom dataset to compare our method qualitatively and quantitatively with previous work [23,27] in which sharpening methods and nonlinear enhancement PDEs are applied to DTI. We observed qualitatively on real data in Fig 11 and quantitatively in Table 1 the advantage of CSD, that allows to use linear enhancements with less extreme regularization parameters than with the DTI based method, resulting in a more reliable tractography.
For our second approach to improve fiber alignment, we introduced a fiber to bundle coherence measure that can be used for detecting and filtering spurious fibers. The fiber to bundle coherence (FBC) is computed from a tractography based density that we constructed using the same PDE foundation as in the first method. As an application we considered the reconstruction of the optic radiation, a fiber bundle of which the position of the anterior extent (the Meyer's loop) is of interest for temporal lobe resection surgery. Accurate and stable localization of the tip of the Meyer's loop is difficult due to the presence of spurious fibers, as shown in Fig  12. We demonstrated in Figs 13, 14 and 15 that either by enhancement of the CSD FOD, or by removing the most spurious fibers using the FBC measure leads to a robust probabilistic tractography. In particular, the combination of both methods in one pipeline allows for a more stable localization of the tip of the Meyer's loop and a more stable determination of the Meyer's loop to Temporal Pole distance.
Our experiments show that our PDE enhancement methods for contextual processing are an effective and widely applicable tool to both enhance CSD data and to remove spurious fibers from tractographies. While we used CSD to construct an FOD, the PDE enhancement can be applied to an FOD obtained with any other method. We have seen that both our methods improve fiber alignment in tractography results and hence provide information on structural connectivity of the brain white matter more robustly. In the future, we aim to improve this framework by using data-adaptive smoothing, for example using local gauge frames [79].