We present a novel open-source module, included in the DIPY (Diffusion Imaging in Python) software library¹, implementing the contextual PDE framework of Portegies et al.² for processing HARDI data. Here we propose an efficient implementation based on multithreaded processing and pre-computed lookup tables that enable large scale experiments. We demonstrate its potential in crossing-preserving enhancement of ODF/FOD fields where the aim is to enhance the alignment of elongated structures such that crossings are maintained.

The contextual enhancement method is based on a linear (hypo-elliptic) PDE in the domain of coupled positions and orientations³ $$$\mathbb{R}^3 \rtimes S^2$$$. This domain carries a non-flat geometrical differential structure, that allows including a notion of alignment between neighboring points². Let $$$({\bf y},{\bf n}) \in \mathbb{R}^3\rtimes S^2$$$ where $$${\bf y}$$$ and $$${\bf n}$$$ denote the spatial and angular part, respectively. Let $$$W:\mathbb{R}^3 \rtimes S^2 \times \mathbb{R}^+ \to \mathbb{R}$$$ be the function representing the evolution of FOD/ODF field. Then, the contextual PDE with evolution time $$$t\geq 0$$$ is given by:

$$\begin{cases} \frac{\partial}{\partial t} W({\bf y},{\bf n},t) &= ((D^{33}({\bf n} \cdot \nabla)^2 + D^{44} \Delta_{S^2})W)({\bf y},{\bf n},t) \\ W({\bf y},{\bf n},0) &= U({\bf y},{\bf n}), \end{cases}$$

where $$$D^{33}>0$$$ is the coefficient for the spatial smoothing which occurs strictly in the direction of $$$n$$$; $$$D^{44}>0$$$ is the coefficient for the angular smoothing ($$$\Delta_{S^2}$$$ denotes the Laplace-Beltrami operator on the sphere $$$S^2$$$); $$$U:\mathbb{R}^3\rtimes S^2 \to \mathbb{R}$$$ is the initial condition. The equation is solved via a shift-twist convolution (denoted by $$$\ast_{\mathbb{R}^3\rtimes S^2}$$$) with its corresponding kernel $$$P_t:\mathbb{R}^3\rtimes S^2 \to \mathbb{R}^+$$$:

$$W({\bf y},{\bf n},t) = (P_t \ast_{\mathbb{R}^3 \rtimes S^2} U)({\bf y},{\bf n}) = \int_{\mathbb{R}^3} \int_{S^2} P_t (R^T_{{\bf n}^\prime}({\bf y}-{\bf y}^\prime), R^T_{{\bf n}^\prime} {\bf n} ) U({\bf y}^\prime, {\bf n}^\prime) \hspace{0.2em} d\sigma ({\bf n}^\prime) d{\bf y}^\prime.$$

Here, $$$R_{\bf n}$$$ is a 3D rotation that maps the vector $$$(0,0,1)$$$ onto $$${\bf n}$$$ and $$$d\sigma({\bf n^\prime})$$$ is the usual surface measure on the sphere $$$S^2$$$. The kernel $$$P_t$$$ has a stochastic interpretation, see Fig. 1. It can be seen as the limiting distribution obtained by accumulating random walks of particles in the position/orientation domain, where in each step the particles can (randomly) move forward/backward along their current orientation, and (randomly) change their orientation. In practice, as the exact analytical formulas for the kernel $$$P_t$$$ are yet unknown, we use the recent approximation given by Portegies et al.⁴ The shift-twist convolutions are implemented using the Cython language and multithread processing via the OpenMP library. Following ideas in Rodriguez et al.⁵, further speedup is achieved by computing lookup-tables containing rotated versions of the kernel $$$P_t$$$ sampled over a discrete set of orientations. In order to ensure rotationally invariant processing, the discrete orientations are required to be equally distributed over a sphere, obtained by an electrostatic repulsion model⁶.

An example illustrating the method is performed on the Stanford HARDI dataset⁷ (150 orientations, b=2000s/mm^2) where Rician noise is added. CSD⁸ is first applied to create the FOD field in Fig. 2b. The result of the enhancement for $$$D_{33}=1$$$, $$$D_{44}=0.02$$$ and $$$t=1$$$ is depicted in Fig. 2c. For those parameters, the kernel is sampled on a 7x7x7 spatial grid and 100 orientations. Optionally, in order to recover sharp angular distributions, one can apply the Spherical Deconvolution Transform (SDT)⁹ to obtain the results in Fig. 2d. For comparison we applied CSD to the original data (the data without the added noise), see Fig 2a. The better alignment of glyphs shows that a significant enhancement of the coherent structures is achieved. This is also supported quantitatively by computation of the average angular error¹⁰ w.r.t to the original data, which is 23.5° for the enhanced data and 33.8° for the noisy data. Regarding computation time, after a one-time computation of the lookup-table (taking several minutes) for the same experiment a speed up of 209% is achieved w.r.t. the finite difference implementation proposed by Creusen et al¹¹ (explicit scheme, $$$\Delta t = 0.05$$$). Although contextual PDEs had proven to be useful in different clinical applications^2,12, their use had been limited due to the lack of flexible implementations or their inherent high computational demands. Our new open source module available in DIPY makes the framework presented here widely available to the neuroimage community and flexible to be included within the existing processing pipeline of DIPY. The novel implementation of the shift-twist convolution via multithreading and lookup-tables allow for large scale experiments. A similar framework can be applied to other contextual PDE applications, e.g. morphological operations in ODF fields¹³. Fiber-to-bundle coherence measures² can be obtained using the same routines and will be included in DIPY.