Thinking Outside the Voxel: A Joint Spatial-Angular Basis for Sparse Whole Brain HARDI Reconstruction

Evan Schwab^{1}, Rene Vidal^{2}, and Nicolas Charon^{3}

A HARDI signal $$$S(v,\vec{g})\in \Omega \times \mathbb{S}^2$$$ is a function of two variables: $$$v$$$, and $$$\vec{g}$$$ the gradient direction, or angular location, on the unit sphere $$$\mathbb{S}^2$$$. Let $$$\Gamma$$$ be a basis on $$$\mathbb{S}^2$$$ then we can write $$$S(v,\vec{g}) = \sum_i a_i(v) \Gamma_i(\vec{g})$$$. We notice that each $$$a_i$$$ is a 3D volume so they can each be represented by a spatial basis for $$$\Omega$$$, say $$$\Psi$$$, as $$$a_i(v) = \sum_j c_{i,j} \Psi_j(v)$$$. Combining equations we have $$$S(v,\vec{g}) = \sum_i \sum_j c_{i,j} \Psi_j(v)\Gamma_i(\vec{g})$$$. Therefore we arrive at a joint spatial-angular basis representation of the HARDI signal. In matrix notation, $$$s = \Phi c$$$ where $$$s$$$ is the vectorized HARDI tensor $$$S$$$, $$$c$$$ is the vectorized $$$C = [c_{i,j}]$$$ and $$$\Phi = \Psi \times \Gamma$$$, the Kronecker product. To find a sparse global parameterization $$$c$$$ we aim to solve $$c = arg\min_c || s - \Phi c || \ \ s.t \ \ ||c||_0 < k.$$ We use an Orthogonal Matching Pursuit (OMP) algorithm designed to exploit the structure of $$$\Phi$$$. In particular, we choose $$$\Psi$$$ to be an orthonormal 3D wavelet basis since MRIs are known to be sparse in this domain. With $$$\Psi$$$ orthonormal we need not compute explicitly $$$P = \Phi^\top \Phi$$$ for OMP but can pre-compute the smaller $$$G = \Gamma^\top \Gamma$$$ and sparsely populate $$$P(i_{\Gamma},j_{\Psi}) = G(i_{\Gamma},j_{\Gamma})$$$ if the next chosen index $$$i’_{\Psi} = i_{\Psi}$$$ and zero otherwise.
For this work we explore two angular bases, $$$\Gamma$$$, the spherical harmonic (SH) basis and the over-complete spherical ridgelet (SR) basis^{1} known to provide a sparse ODF representation in the spherical wavelet (SW) domain.

[1] Tristán-Vega, Antonio, and Carl-Fredrik Westin. "Probabilistic ODF estimation from reduced HARDI data with sparse regularization." Medical Image Computing and Computer-Assisted Intervention (MICCAI) 2011, pp. 182-190.

[2] Merlet, S. L. and Caruyer, E. and Ghosh, A. and Deriche, R. "A Computational Diffusion MRI and Parametric Dictionary Learning Framework for Modeling the Diffusion Signal and its Features." Medical Image Analysis. 2013, 17(7), pp. 830-843.

[3] Gramfort, A. and Poupon, C. and Descoteaux, M. "Denoising and Fast Diffusion Imaging with Physically Constrained Sparse Dictionary Learning." Medical Image Analysis. 2014, 18(1) 36-49.

[4] Ye, W. and Vemuri, B. C. and Entezari, A. "An over-complete dictionary based regularized reconstruction of a field of ensemble average propagators." International Symposium on Biomedical Imaging (ISBI) 2012, pp. 940-943.

[5] Gramfort, A. and Poupon, C, and Descoteaux, M. "Sparse DSI: Learning DSI structure for denoising and fast imaging." Medical Image Computing and Computer-Assisted Intervention (MICCAI) 2012. 288-296.

[6] Wu, Y. and Feng, Y. and Li, F. and Westin, C. F. "Global consistency spatial model for fiber orientation distribution estimation.” International Symposium on Biomedical Imaging (ISBI) 2015, pp. 1180-1183.

[7] Cheng, J. and Shen, D. and Basser, P. J. and Yap, P. T. "Joint 6D kq Space Compressed Sensing for Accelerated High Angular Resolution Diffusion MRI." Information Processing in Medical Imaging (IPMI) 2015, pp. 782-793

[8] Ourselin, S. and Alexander, D. C. and Westin, C.-F., and E Cardoso, M. J. "Leveraging EAP-Sparsity for Compressed Sensing of MS-HARDI in (k, q)-Space." Information Processing in Medical Imaging (IPMI) 2015, pp. 375-386.

[9] Ouyang, Y., Chen, Y., Wu, Y., & Zhou, H. M. "Total Variation and Wavelet Regularization of Orientation Distribution Functions in Diffusion MRI." Inverse Problems and Imaging, 2013, 7(2): 565-583.

[10] Ouyang, Y. and Chen, Y. and Wu, Y. "Vectorial total variation regularisation of orientation distribution functions in diffusion weighted MRI."International Journal of Bioinformatics Research and Applications, 2014, 10(1): 110-127.

Figure 1: Comparison between the proposed method (Angular: SH basis, Spatial: 2D Haar wavelets using k=250,500,2500 atoms (Col. 1-3)) and the GT (SH per-voxel (Col. 4)). The second row offers a closeup within the bounding box. We notice a slower progression of detail focus than using SR basis (Figure 2).

Figure 2: Comparison between the proposed method (Angular: SR basis, Spatial: 2D Haar wavelets using k=250,500,2500 atoms (Col. 1-3)) and the GT (SH per-voxel (Col. 4)). The second row offers a closeup within the bounding box. We notice

Figure 3: Normalized mean square error (NMSE) between the proposed method using an angular SH or SR basis and the GT per-voxel SH method for increasing sparsity level. We end at 2500 atoms (93.33% sparse), to show what performance we achieve with an average of only 1 atom per voxel.

Figure 4: Visualization of the first 3 atoms chosen by the OMP algorithm for the proposed SR and 2D Haar wavelet product basis, illustrating the spatial redundancies captured by a single spatial-angular atom. Top: Angular atom (ODF). Bottom: Spatial-angular atom.

Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)

2041