Thinking Outside the Voxel: A Joint Spatial-Angular Basis for Sparse Whole Brain HARDI Reconstruction
Evan Schwab1, Rene Vidal2, and Nicolas Charon3

1Electrical and Computer Engineering, Johns Hopkins University, Baltimore, MD, United States, 2Biomedical Engineering, Johns Hopkins University, Baltimore, MD, United States, 3Applied Mathematics and Statistics, Johns Hopkins University, Baltimore, MD, United States

### Synopsis

Sparse modeling of dMRI signals has become of major interest for advanced protocols such as HARDI which require a large number of q-space measurements. With few exceptions, prior work have developed bases to sparsely represent q-space signals per voxel with additional constraints of spatial regularity. In this work, we propose a single basis to represent an entire HARDI dataset by modeling spatial and angular domains jointly, achieving an unprecedented level of sparsity. With a globally compressed representation we can then redefine HARDI processing, diffusion estimation, feature extraction and segmentation, and drastically reduce acquisition time and data storage.

### Introduction

Diffusion magnetic resonance imaging (dMRI) is a 6D neuroimaging modality that measures 3D q-space signals at every voxel in a brain MRI volume to estimate the orientation of neuronal fiber bundles, in vivo. Most every dMRI reconstruction algorithm1,2,3 finds a basis to represent a q-space signal and adds terms to enforce sparsity, orientation distribution non-negativity or spatial regularity based on similar basis representations within neighboring voxels. Some recent methods have considered sparse optimization over multiple voxels or an entire volume jointly to reduce redundancies and invoke implicit regularization but all of these methods still only operate on a q-space basis4,5,6,7,8,9,10. In contrast, we propose a single basis representation that models q-space and the spatial domain jointly resulting in a drastically compressed model of an entire brain dMRI dataset. This sparse global representation has the potential to push the limits of dMRI acquisition time, data storage, diffusion estimation, fiber segmentation, and other processes. In this work we consider high angular resolution diffusion imaging (HARDI) which simplifies general positions in 3D q-space to 2D locations on the unit q-sphere, but our main idea can be used for any 6D dMRI protocol. We propose a single joint spatial-angular basis for sparse whole brain HARDI reconstruction and orientation distribution function (ODF) estimation.

### Methods

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) basis1 known to provide a sparse ODF representation in the spherical wavelet (SW) domain.

### Results

For our experiments we used the ISBI 2013 HARDI Challenge Phantom dataset (50x50x50 volume with 64 gradient directions, SNR=30). For simplicity we chose a single slice (z=25). As a full, non-sparse "ground truth" (GT) representation we used the traditional SH (order $L=4$) reconstruction per voxel by least squares, which exudes a total of 15x50x50 = 37,500 atoms. Figures 1 and 2 show a qualitative comparison between the GT and the proposed method using the SH and SR bases, respectively with a fraction of the total atoms. We notice an inherent denoising due to the spatial wavelets. In Figure 3 we present a quantitative comparison of the normalized mean squared error (NMSE) between the proposed SH and SR reconstructions and the GT as we increase the sparsity level. Then Figure 4 illustrates how spatial and angular redundancies are encoded by the global atoms over the entire dataset as the OMP algorithm progresses. As an important note, using 2500 atoms for this dataset achieves a comparably accurate reconstruction and amounts to an average of only 1 atom per voxel, which is unachievable with a sparse angular basis alone. By exploring new spatial-angular bases and dictionary learning techniques we believe we can achieve even more accurate reconstructions with much fewer atoms.

### Conclusion

We have thus presented a framework for encoding an entire HARDI dataset with a single joint spatial-angular basis which achieves an unprecedented level of sparsity. This framework is important for globally subsampling HARDI signals to speed up acquisition, enforcing ODF non-negativity globally, spatial-angular feature extraction, fiber segmentation and data storage.

### Acknowledgements

This work funded by JHU start-up funds.

### References

[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.

### Figures

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