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 algorithm^1,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 basis^{4,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.

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

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