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.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.⁷ 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 projection^4,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.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.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.