Using A Hyperspherical Harmonic Basis for Sparse Multi-Voxel Modeling of Diffusion MRI
Evan Schwab1,2, Hasan Ertan Cetingul2, Rene Vidal3, and Mariappan Nadar2

1Electrical and Computer Engineering, Johns Hopkins University, Baltimore, MD, United States, 2Medical Imaging Technologies, Siemens Healthcare, Princeton, NJ, United States, 3Biomedical Engineering, Johns Hopkins University, Baltimore, MD, United States

### Synopsis

For diffusion magnetic resonance imaging (dMRI), 3D signals are acquired at each voxel to estimate neuronal fiber orientation in the brain. Traditionally, dMRI signals are reconstructed using the same basis for each voxel with added spatial regularization and sparsity constraints. By repeating the same basis for each voxel, there exist millions of redundant parameters to represent an entire brain volume. In this work, we reconstruct dMRI signals jointly across multiple voxels to reduce the number of parameters needed to represent a brain volume by 65%. Sparse dMRI representation is important for storage, information extraction, and reduction of clinical scan times.

### Introduction

Traditionally dMRI signals are reconstructed per-voxel.1,2,3 The signal measured in each voxel is reconstructed with the same basis, (e.g., spherical harmonics, SHORE, BFOR, spherical ridgelets) and spatial regularity and sparsity constraints may be added to the optimization using compressed sensing (CS) to reduce the number of coefficients in each voxel.4,5,6 However, by repeating the same basis for each voxel, there are inherent redundancies, and by reconstructing multiple signals simultaneously we can more greatly reduce the total number of signal measurements needed. Our goal is to develop a multi-voxel basis representation of dMRI data that will drastically reduce the total number of parameters needed to reconstruct dMRI signals.

### Methods

Our work is motivated by Hosseinbor et al.4,5,6 who model multiple disjointed 3D anatomical shapes with a single 4D basis, $Z_{nl}^m(\beta,\theta,\phi)$, the hyperspherical harmonic (HSH) basis on $\mathbb{S}^3$, instead of modeling each 3D shape separately using spherical harmonics.7 This representation captures the spatial relationship between the disjoint objects whereas modeling them separately ignores potential redundancies which can be exploited for sparsity. In this vein, we can view 6D dMRI signals $S(x,q) \in \mathbb{R}^3 \times \mathbb{S}^2$ as multiple disjointed sets of 3D points centered at each voxel of the MRI volume which can be represented jointly with the HSH basis to model redundancies between neighboring voxels. Figure 1 illustrates this parameterization for a patch of 4 voxels of single-shell HARDI signals with the origin $O$ and the coordinates of $S(x,q) = S(g) \in \mathbb{R}^3$ which depend on the voxel location $x$, the gradient vector $q$, and the physical voxel dimensions. Admittedly, this particular parameterization and basis representation are more useful from a purely signal processing perspective (whereby once the signals are measured one may disregard their biophysical meaning in q-space). Though this may be seen as a model weakness, it does provide a generality for which it can be applied to any dMRI acquisition scheme such as multi-shell or a cartesian grid for DSI since the proposed representation can be applied to any set of points in 3D space. This means that unlike prior work which need to choose a per-voxel basis that depends on the biophysics of the acquisition, the proposed framework is universal. Our 3D signals $S(g)$ are projected onto the 4D sphere using inverse stereographic projection4,5,6 and we can write $S(\beta,\theta,\phi) = \sum_{nlm} c_{nl}^m Z_{nl}^m (\beta,\theta,\phi)$ a linear combination of HSH basis functions (see Figure 2 for an example of signal basis atoms after converting to ODFs). In discrete form, for an arbitrary patch of signals, we wish to solve: $$\mathbf{c} = arg\min_\mathbf{c} ||\mathbf{S} - \mathbf{Z}\mathbf{c}||_2^2 \ \ s.t. \ \ ||\mathbf{c}||_0 \le k.$$ In this work we use Orthogonal Matching Pursuit (OMP) to solve the $L_0$ problem. We also compare against least squares $L_2$ solution.

### Results

We evaluated our multi-voxel method on synthetic, phantom, and real HARDI brain data compared to single voxel reconstruction with SH basis. Figure 3 (top) shows qualitative results on a synthetic dataset (SNR=5) consisting of quadrants of signals corresponding to 0-, 1-, and 2-fiber orientation distribution functions (ODFs). Figure 3 (bottom) shows the corresponding HSH coefficients for each quadrant. It is observed that the HSH representation produces ODFs that are more similar to the ground truth while using only 368 coefficients, as opposed to 4,500 coefficients in the case of the per-voxel SH representation. In particular, the patch of isotropic ODFs are properly represented with a single isotropic basis atom, while patches for 1- and 2-fibers require 121 and 125 atoms respectively, an average of ~5 atoms per voxel. For the phantom dataset we used the ISBI 2013 HARDI Reconstruction Challenge with SNR=10. The real human brain HARDI data was acquired with 127 gradient directions and b = 1000s/mm$^2$. Table 1 summarizes the reconstruction results for each dataset using single-voxel SH reconstruction and non-overlapping patch based HSH reconstruction in terms of the minimum number of coefficients needed to achieve a certain level of accuracy calculated by mean square error (MSE). We used the sparsity and HSH parameters to empirically find a suitable tradeoff between complexity and accuracy, and we observe a drastic decrease in the number of coefficients while maintaining error.

### Conclusion

We have developed a novel framework of HARDI reconstruction that removes the long standing tradition of voxel-wise computation and drastically reduces the total number of coefficients needed to encode a HARDI dataset. This framework opens the doors for new patch-based methods for HARDI processing and analysis such as multi-voxel fiber orientation estimation, feature extraction, segmentation and CS-based signal measurement reduction to speed up HARDI acquisition.

### Acknowledgements

Evan Schwab performed this work while at Siemens Healthcare.

### References

1. Gramfort A, Poupon C, Descoteaux M. Denoising and fast diffusion imaging with physically constrained sparse dictionary learning. Med Image Anal. 2014; 18(1):36-49.

2. Michailovich O, Rathi Y. Fast and accurate reconstruction of HARDI data using compressed sensing. In: MICCAI 2010; 607-14.

3. Merlet S, Caruyer E, Ghosh A, et al. A computational diffusion MRI and parametric dictionary learning framework for modeling the diffusion signal and its features. Med Image Anal. 2013; 17(7):830-43.

4. Hosseinbor A, Chung M, Schaefer S, et al. 4D hyperspherical harmonic (HyperSPHARM) representation of multiple disconnected brain subcortical structures. In: MICCAI 2013; 598-605.

5. Hosseinbor A, Kim W, Adluru N, et al. The 4D hyperspherical diffusion wavelet: a new method for the detection of localized anatomical variation. In: MICCAI 2014; 65-72.

6. Hosseinbor A, Chung M, Wu Y, et al. A 4D hyperspherical interpretation of q-space. Med Image Anal. 2015; 21(1):15-28.

7. Styner M, Oguz I, Xu S, et al. Framework for the statistical shape analysis of brain structures using SPHARM-PDM. The Insight Journal. 2006; 1071:242.

### Figures

Figure 1: Illustration of multi-voxel parameterization for a patch of 4 single-shell HARDI signals (top view) with center $p_0$, and q-space vector $q$. The global coordinate $g = p_0 + q$ depends on voxel size and gradient direction $q$. Any array of 3D points can be reconstructed with this parameterization.

Figure 2: Examples of HSH signal basis elements (after conversion to ODFs for visualization) for a 1×2 patch of voxels. Notice that all different types of ODF configurations are presented for a given HSH atom and that neighboring ODFs, regardless of similarity, can be jointly modeled with very few atoms.

Figure 3: (Top) $L_2$ reconstruction of a 10×10×3 synthetic HARDI signals with SNR=5. SH order L=4 requires 15×(10×10×3)=4,500 coefficients. For HSH order N=16, we used four 5×5×3 non-overlapping patches requiring 1 (0-fiber quad) + 2×121 (1-fiber quads) + 125 (2-fiber quad) = 368 coefficients. (Bottom) HSH coefficients for each quadrant.

Table 1: Results of experiments on synthetic, phantom and real HARDI datasets comparing the total number of coefficients and the average MSE of the entire datasets using single-voxel SH and multi-voxel HSH $L_2$ and $L_0$ reconstruction. We can see a drastic decrease in the number of coefficients while maintaining error.

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