Introduction¶
- This study aims to improve the acquisition time and/or quality of reconstruction of tractographies based on diffusion weighted MRI data.
- Current approaches suffer from:
- Low resolution sampling but high resolution desired outcome (2-3mm scan resolution, good for major white matter pathways, but not for micrometer axons)
- Long acquisition times
- Reconstruction inaccuracies due to uncertainty in gradient directions, dropped directions, low SNR, model unknowns (# directions etc.)
- This project:
- Single-shell HARDI imaging: at each voxel, presented with values sampled at a constant b-value but in varying directions
- Q-space is sampled sparsely: saves acquisition time
- Alternative approach: diffusion spectrum imaging [Eleftherios] or QBI
- In which gradient directions should Q-space be sampled? Determined by the gradient directions, and controlled by us. Typical grids:
- Minimum energy grid
- Random grid
- Optimal quadrature points (introduced by this study)
- Random grid: advantage of arbitrary # of points, points can be dropped, variable scan duration
- Optimal quadrature points: "optimal" representation, as far as integration over the sphere goes
- Integration over the sphere is important when computing coefficients of spherical harmonics via inner product
- Higher b-value, lower SNR, higher directional selectivity
- ^ Also, error in control of b-directions due to D2A conversion
- Typical b-values of around 2000
- SNR loss with increased b: b proportional to G^2 t^3. We'd like short t, have to generate stronger G. b of 2000 (relatively strong), t_E around 110ms...signal already fairly decayed.
- Prefer to jump around in gradient directions. Heat induction in coils, try not to hit same parts each time to allow optimal cooling.
- Eddy currents induced around activated coil. Oxford group, Jesper.
Mathematical formulation¶
- Under certain moderate assumptions, Stejskal and Tanner [?] shows that signal attenuation is related to the diffusion PDF by the 3-D Fourier Transform:
$\frac{S(\mathbf{q}, \tau)}{S_0} = \int_\mathcal{R^3} P(\mathbf{r} | \mathbf{r_0},\tau)\exp (-2 \pi i \mathbf{q}^\mathrm{T}\mathbf{R})\mathrm{d}\mathbf{r} = \mathcal{F}[ P(\mathbf{r} | \mathbf{r_0}, \tau) ]$
with
$\mathbf{q}=\gamma \delta \mathbf{G} / 2 \pi$, $\delta$ the gradient pulse duration, $\gamma$ the nuclear gyromagnetic ratio of water protons, $\mathbf{G}$ the applied gradient vector, and $S_0$ the baseline image acquired without any gradient ($b=0$). $P$ is the diffusion PDF.
Discuss PDF structure $P(\mathbf{r} | \mathbf{r_0},\tau)$
- An important special case occurs when the PDF is assumed to be Gaussian:
$S(\mathbf{q}, \tau) = S_0 e^{-\tau \mathbf{q}^\mathrm{T} D \mathbf{q}}$
or, with normalized $\mathbf{q}$ and $b = \tau |\mathbf{q}|^2$,
$S(b, \mathbf{g}) = S_0 e^{-b \mathbf{g}^\mathrm{T} D \mathbf{g}}.$
We use this relationship in construction artificial data-sets for our simulations (diffusion tensor model).
We are interested in the Orientation Distribution Function, related to the diffusion PDF as follows
$\mathrm{ODF}(\hat{u}) d\Omega = \int_{r=0}^{r=\infty}P(r\hat{u}) dv = \int_{r=0}^{r=\infty}P(r\hat{u})r^2 dr d\Omega$
i.e. PDF over directions (surface of unit sphere) of diffusion, s.t.
$\mathrm{ODF}(\hat{u}) = \int_0^\infty P(r\hat{u})r^2 dr$
The above $r^2$ term is important to include; was neglected in the widely cited paper by Tuch et al.
- Based on the 3-D FFT relationship, we will proceed to model the reconstruction as a linear algebra problem, based on discretization of spherical harmonics.
Signal model¶
Aganj shows that the ODF, with the above $r^2$ factor taken into account, can be expressed in terms of $P$'s 3D Fourier transform $E$, as
$\mathrm{ODF}(\hat{u}) = \frac{1}{4\pi} - \frac{1}{8\pi^2} \int_0^{2\pi} \int_0^\infty \frac{1}{q}\nabla_b^2 E(\mathbf{q}) dq d\phi$
where $E(q_0 \hat{u}) = S(\hat{u}) / S_0$ and $|\mathbf{q}| = q_0$.
Since we only have values for $E$ on a sphere in Q-space, we need to model values elsewhere. If we assume the model
$E(\mathbf{q}) \approx E(q_0 \hat{u})^{q^2 / q_0^2} = \tilde{E}(\hat{u})^{q^2 / q_0^2}$
(useful because correctly has $E(0) = 1$)
we can show that
$\mathrm{ODF}(\hat{u}) = \frac{1}{4 \pi} + \frac{1}{16 \pi^2} \mathrm{FRT} \\{ \Delta_b^2 \ln(-\ln \tilde{E}(\hat{u})) }.$
Here, $\mathrm{FRT}$ is the Funk-Radon Transform (integral of function on great circle perpendicular to evaluated vector). While it can be tricky to evaluate directly, once the function is expressed as a combination of spherical harmonics it can be computed analytically.
From the above, note that there are two important operators we need to be able to compute: FRT and spherical Laplacian.
Spherical harmonics representation of ODF¶
Consider the subspace of spherical harmonics of degree N:
$\mathcal{H_N} = \mathrm{span} \\{ Y_l^m : |m| \leq l, 0 \leq l \leq N \\}$
A function in this subpace can be represented as a linear combination of the spherical harmonics (associated Legendre polynomials). It can be shown (see notes by Cory) that, for a function $f \in \mathcal{H_N}$, we have a reproducing kernel
$K(\mu) = \sum_{n=0}^{N} \frac{2n + 1}{4 \pi} P_n(\mu)$
such that
$f(\Omega) = \int_\mathbb{S^2} K(\Omega \cdot \Omega') f(\Omega') d\Omega'.$
If an exact quadrature exists on the sphere, with nodes and weights $\\{\Omega_i, w_i \\}_{i=1}^M$, we have
$f(\Omega) = \sum_{i=1}^{M} w_i K(\Omega \cdot \Omega_i) f(\Omega_i).$
In other words, the function can be computed at any position by adding up weighted kernels at the measurement positions. Or: given measurements at sparse measurement positions, we have knowledge of the entire function.
Note that we can also restrict the kernel above to always be symmetric (applicable in the case of the ODF), which yields:
$K_e(\mu) = \sum_{n=0}^{\lfloor N/2 \rfloor} \frac{ 2(2n) + 1}{4\pi} P_{2n}(\mu).$
Operators on the function¶
We need to investigate the effect of the inverse radon transform and inverse spherical Laplacian operators on functions represented by SPH; but first, let's see examine the effect of these two operators on the associated Legendre polynomials in general.
Spherical Laplacian¶
The spherical harmonics are eigenfunctions of the spherical Laplacian, s.t.
$\Delta Y_n^m(\Omega) = -n(n+1)Y_n^m(\Omega).$
Therefore, computing the Laplacian of the kernel becomes:
$\Delta K_e(\Omega \cdot \Omega') = - \sum_{n=0}^{\lfloor N/2 \rfloor} (2n)(2n+1)\left[ \frac{2(2n)+1}{4 \pi} \right] P_{2n}(\Omega \cdot \Omega').$
¶
Funk-Radon Transform¶
$\mathcal{G}(\Delta K_e(\Omega \cdot \Omega_i))(\Omega) = \int_\mathbb{S^2} \delta(\Omega \cdot \Omega') \Delta K_e (\Omega' \cdot \Omega_i)d\Omega'$
After some manipulation (refer to Cory's notes):
$= - \sum_{n=0}^{\lfloor N/2 \rfloor} (2n)(2n+1)[2\pi P_{2n}(0)]\left[ \frac{2(2n) + 1}{4\pi}\right]P_{2n}(\Omega \cdot \Omega_i)$
Inverse transforms¶
Note that, since the operators simply multiply the expression by additional factors, we can easily move in forward or reverse directions. I.e., given the FRT of the Laplacian of the function, we can reconstruct the function, or vice versa. This means that we can freely travel between ODF and Signal Space.
We can define a kernel based on the inverse spherical Laplacian and inverse Funk-Radon transform:
$H(\mu) = - \sum_{n=1}^{\lfloor N/2 \rfloor} \frac{2(2n)+1}{8 \pi^2 P_{2n}(0)(2n)(2n+1)}P_{2n}(\mu)$,
i.e.
$\mathcal{G}\Delta H(\Omega \cdot \Omega_i) = K_e(\Omega \cdot \Omega_i)$.
This gives us two kernels: one associated with signal space and another with ODF space. These kernels allow us to move between signal and ODF space:
$\mathrm{ODF}(\Omega) = \frac{1}{4\pi} + \frac{1}{16\pi^2}\mathcal{G}\left\{ \Delta \ln \left[ -\ln(\tilde{E}) \right] \right\}$
Now assume that
$\ln \left[-\ln \tilde{E}(\Omega) \right] \approx \sum_{i=1}^M \phi_i H(\Omega \cdot \Omega_i)$
then we also have
$\mathrm{ODF}(\Omega) \approx \frac{1}{4\pi} + \frac{1}{16\pi^2} \sum_{i=1}^M \phi_i K_e(\Omega \cdot \Omega_i).$
Extending to linear algebra¶
Note that our signal approximation above gives a set of linear equations:
$\ln \left[ -\ln(\tilde{E}(\Omega_j)) \right] \approx \sum_{i=1}^M \phi_i H(\Omega_j \cdot \Omega_i),\quad j=1,2,3,\ldots,M'$
or
$\hat{E} \approx A \mathbf{x}.$
Where each column of the matrix A contains the kernel $i$, evaluated at the different sampling points $\Omega_j$. Note that the evaluation points $\Omega_i$ are chosen arbitrarily.
Also: this allows us to evaluate either the signal or the ODF at arbitrary points by simply performing a matrix multiplication of weights times relevant kernel.
Sparse fitting theory¶
We note that, while the signal space is densely populated, the ODF can be represented fairly sparsely. E.g., it should be possible to accurately represent peaks in four directions using non-zero coefficients on only two even kernels. Of course, that would require for the kernels to be appropriately positioned, but the general idea holds.
So, how is it possible to have the same set of coefficients represent both dense and sparse functions? Examine the kernel shapes (ODF vs Q-space):
from sphdif.kernel import kernel_plot, even_kernel, inv_funk_radon_even_kernel
from sphdif.plot import get_mlab
mlab = get_mlab()
kernel_plot(even_kernel, N=14)
mlab.show()
This is the idea behind compressed sensing. Think, for example, of the sum of three cosines. In the time-domain, this signal has positive magnitude values for most values of time, but in the frequency domain, there are only 6 non-zero coefficients to estimate.
As such, we solve the linear system above with an $\ell_1$-penalization to enforce sparsity:
$\hat{E} = A\mathbf{x} + \lambda |\mathbf{x}|_1, \quad x_i >= 0.$
It can be solved using various optimization techniques, including the lasso from scikit-lean.
Solving for directions¶
Given the signal, what are the fiber directions and how many are there? In order to simulate the problem, we model a signal as the sum of single-fiber diffusion:
$S = S_0 \sum_i e^{-b \mathbf{g_i}^\mathrm{T} D \mathbf{g_i}}.$
Future work¶
- Automated determination of number of angles (K-means, BIC?)
- Investigate sensitivity to dropped directions, true vs. real gradients on reconstruction
- Include random grid in comparison (useful: no problem with dropped data, acquisition as long as patient is comfortable)
- Comparison of real-world data to model fit--but how? Visual? Prediction?
References¶
- I. Aganj et al., Reconstruction of the Orientation Distribution Function in Single- and Multiple-Shell q-Ball Imaging Within Constant Solid Angle, Magnetic Resonance in Medicine 64:554-566 (2010).
- M. Descoteaux, High Angular Resolution Diffusion MRI: from Local Estimation to Segmentation and Tractography, PhD thesis, University of Nice-Sophia Antipolis, 2008.
- E. Garyfallidis, Towards an Accurate Brain Tractography, PhD thesis, Cambridge University, 2012.
- R. Heidemann et al., k-space and q-space: Combining ultra-high spatial and angular resolution in diffusion imaging using ZOOPPA at 7T, preprint: Neuroimage (2011).
- E. Stejskal and J. Tanner, Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient, Journal of Chemical Physics, 42:288-292 (1965).