Fast Fiber Orientation Estimation in Diffusion MRI from kq-Space Sampling and Anatomical Priors

High spatio-angular resolution diffusion MRI (dMRI) has been shown to provide accurate identification of complex neuronal fiber configurations, albeit, at the cost of long acquisition times. We propose a method to recover intra-voxel fiber configurations at high spatio-angular resolution relying on a 3D kq-space under-sampling scheme to enable accelerated acquisitions. The inverse problem for the reconstruction of the fiber orientation distribution (FOD) is regularized by a structured sparsity prior promoting simultaneously voxel-wise sparsity and spatial smoothness of fiber orientation. Prior knowledge of the spatial distribution of white matter, gray matter, and cerebrospinal fluid is also leveraged. A minimization problem is formulated and solved via a stochastic forward–backward algorithm. Simulations and real data analysis suggest that accurate FOD mapping can be achieved from severe kq-space under-sampling regimes potentially enabling high spatio-angular resolution dMRI in the clinical setting.


Introduction
Diffusion Magnetic Resonance Imaging (dMRI) is a unique non-invasive technique to infer the microscopic architecture of tissues in vivo. In recent years, dMRI has gained a lot of attention in neuroscience since it enables the mapping of the white matter fiber paths, revealing the existing connections between different brain areas [1][2][3]. In clinics, dMRI has shown to provide insights into many neurodegenerative diseases, such as Schizophrenia and Alzheimer's disease [4,5]. Diffusion MRI enables to map the restricted diffusion of the water molecules comprising the white matter tissue. The information captured by dMRI is then processed in order to infer the connectivity and the integrity of the white matter pathways.
A typical approach to trace the complex pathways of the white matter fiber bundles from dMRI signals relies on piecing together local fiber orientation information. Such local information is obtained by processing a multitude of MRI signals generated by applying various diffusion gradients during the acquisition. Each diffusion-weighted (DW) signal is sensitive to diffusion along a specific direction and at a specific intensity, identified by a so-called q-space point defined by the diffusion gradient applied [6].
Diffusion Tensor Imaging (DTI) [7] is the most widely-used technique in clinics for the intra-voxel fiber orientation estimation since it requires less than 7 diffusion volumes, i.e., 7 q-points. DTI relies on a Gaussian description of the water molecules diffusion phenomenon, which enables the presence of only one preferred diffusion direction. For this reason, the presence of multiple fiber populations cannot be distinguished within the same voxel by using DTI, making this technique inappropriate for complex fiber architecture estimation.
High Angular Resolution Diffusion Imaging (HARDI) methods [8] have shown to provide more accurate fiber representations for complex fiber arrangements, by probing the observed brain with a large number of diffusion gradients. For instance, Diffusion Spectrum Imaging [9] considers 515 diffusion gradients allocated on a Cartesian grid, while q-ball approaches [10] use between 30 and 100 diffusion gradients distributed over a single shell. However, acquisitions characterized by large numbers of diffusion gradients result in long scan times, which are unfeasible in clinical applications. This limitation becomes even more obvious when the brain microstructure is investigated, encompassing not only fiber orientation, but also other microstructure parameters (e.g., diameter, etc.).
To this purpose, more elaborate acquisition protocols are used, requiring more diffusion gradients, typically distributed over multiple shells in q-space [11,12]. Consequently, numerous HARDI approaches have been proposed in order to optimize the number and the spatial distribution of the q-space points to be considered for an accurate description on the white matter fiber arrangement [13,14]. More specifically, depending on the microstructure parameters of interest, different approaches have been developed to reduce the dMRI acquisition time while preserving the reconstruction quality.
We focus here on the problem of estimating white matter fiber orientations. A generalization to the study of additional microstructure parameters is beyond the scope of the present work. In this context, some of these methodologies are based on Spherical Deconvolution (SD) approaches, modeling the HARDI signal as a spherical convolution between the Fiber Orientation Distribution (FOD), representing the few active fiber orientations, and a kernel, representing the response signal of a single fiber [15][16][17]. Recently, SD has been utilized to recover the FOD from a reduced number of diffusion gradients (q-space under-sampling) [18][19][20].
The FOD recovery problem can be formulated as an inverse problem, which is illposed and highly sensitive to noise. In this context, an efficient approach consists in leveraging the convex optimization theory by defining the FODs as a solution to a convex regularized minimization problem, incorporating a priori information on the target FODs. In particular, the coefficients of interest are imposed to be non-negative and real valued, since FODs are functions on the unit sphere expressing the orientation and the volume fraction of the fiber populations contained in a single voxel [21].
When expressed as set of discrete orientations, FODs are characterized by coefficients summing up to one in each voxel. Finally, leveraging the fact that, at the imaging resolution available currently, each single voxel is assumed to be populated by only a few fiber bundles, sparsity priors can be incorporated in the model (see, e.g., [18][19][20][22][23][24][25][26]). In this context and inspired by the Compressed Sensing (CS) theory [27,28], many methods have obtained promising results recovering the FOD coefficients from a reduced number of diffusion gradients by leveraging sparsity-based priors.
Constraint Spherical Deconvolution (CSD) [21] represents the first attempt to regularize the FOD recovery problem. Based on the fact that the FOD coefficients correspond to the fiber volume fractions, they are assumed to be non-negative. CSD proposes to iteratively solve a sequence of weighted 2 norm problems whose weights, which depend on the FOD coefficients estimated at the previous iteration, are used to drive to zero negative and small coefficients. Inspired by the CS theory, refs. [18,23] suggested the use of the 1 norm for the recovery of the volume fractions coefficients. However, ref. [22] showed that the use of the 1 norm, meant as simple sum of the coefficients value, is inconsistent with the fact that the volume fractions coefficients sum up to 1 by definition.
Additionally, ref. [22] reformulated the FOD recovery in order to approach the 0 minimization. This is done by solving a sequence of weighted 1 norm problems whose weights at each iteration correspond to the inverse of the values of the solution of the previous problem [29]. At convergence, these weights lead the values of the coefficients to be independent from the magnitude of the non-zero values, thus, mimicking the results of the 0 pseudo-norm.
Ref. [19] proposed the use of spatial regularization to solve the SD ill-posed problem. In this work a piece-wise smoothness of the FOD coefficients is promoted while encouraging the coefficients to provide high contrast. Ref. [19] observed that the contribution of each FOD coefficient is highly correlated to the coefficients associated with its neighbor directions. They proposed to regularize the FOD recovery by penalizing the presence of coefficients that exhibit large variations along similar orientations. In addition, they provided a way to discard the noisy contributions by driving to zero the coefficients that are not sufficiently distant from the mean value of the coefficients of interest.
In addition to the angular resolution, the spatial resolution of the acquired DW volumes is also important to achieve an accurate identification of the white matter paths. In principle, by considering small voxel sizes, the complexity of the inner fiber arrangements can be reduced. Thus, considering both high angular and spatial resolutions would be ideal to ensure an accurate fiber recovery [30][31][32].
Currently, single-shot (SS) DW echo-planar imaging (EPI) is the most popular technique used to perform the dMRI acquisition in clinics. SS-DW-EPI enables extremely fast scan times, which make DW signals nearly immune to patient motion and mostly sensitive to the water molecule movements. The SS-DW-EPI scheme is typically applied in 2D, where each slice is sequentially excited, diffusion-encoded and then collected using a unique trajectory. Despite being very fast, 2D SS-DW-EPI is very inefficient due to the procedure required for the generation of intense diffusion gradients responsible for sensitizing the MRI signal produced by the examined tissues with respect to the diffusion phenomenon occurring along the gradient direction. Indeed, more than 50% of the acquisition time is dedicated to the preparation of the diffusion gradients to be applied, which, in the 2D setting, is required for each single slice.
A number of methods have been proposed in order to acquire DW signals more efficiently and devote the most of the acquisition time to achieve higher resolution and enhance the image quality [33][34][35]. Parallel imaging and 3D diffusion MRI acquisitions are some of them [36]. In particular, moving from the 2D to the 3D approach enables reducing the number of times the diffusion preparation process is applied to cover the same brain volume.
Ideally, a single signal preparation would be enough to encode an entire 3D DW volume. However, this is not the case in practice due to the fast decay DW signals are subject to. Consequently, 3D DW EPI acquisitions typically perform the collection of multiple sub-volumes, which still requires the signal preparation to be applied fewer times than in the 2D setting. In this context, 3D acquisitions can be used in dMRI to reduce the acquisition time required to scan each DW volume or, alternatively, achieve higher signals quality within the same scan time.
In order to accelerate the dMRI acquisition process, previously proposed q-space methods have mainly focused on reducing the number of diffusion gradients from fully sampled DW images. However, the recovery of white matter structures can further benefit from the incomplete acquisition of each k-space volume.
Traditional approaches aiming at exploring the white matter microstructure from under-sampled kq-space data relies on two separate steps: first, DW images are recovered from incomplete k-space data, followed by the recovery of the fiber architecture from a limited number of q-space points.
Recently, methods recovering the diffusion profile directly from under-sampled kqspace data have started to gain significant attention [26,[37][38][39]. The work of [37] developed the first approach to simultaneously recover the diffusion signal and the diffusion propagator from data sub-sampled in both the k and q 3D spaces. The regularization of the problem is performed both in the angular and in the spatial domain. More specifically, the proposed recovery problem promotes the sparsity of the DW images in the wavelet domain, the sparsity of the diffusion signal in the angular domain and the smoothness of the DW images.
In the work of [37], the dictionary used to sparsify the signal in the angular domain was generated by a Dictionary Learning technique that enabled the creation of an adaptive dictionary to fully characterize the signal in each single voxel. The surfacelet transform has been shown to efficiently represent the directional information of the diffusion propagator by using only few coefficients. Inspired by these observations, ref. [40] proposed to recover the full diffusion propagator from measurements that are under-sampled in kq-space by leveraging two sparsity priors. The first prior promotes sparsity of the diffusion propagator in the surfacelet domain through the 1 norm. The second prior promotes sparsity of the diffusion propagator in the gradient domain by means of the TV penalty.
In the work of [41], they used rotation-invariant dictionaries in which only few atoms are active by adapting their orientation to the diffusion MRI data. The recovery framework proposed by [41] consists of alternating between the estimation of the coefficients identifying the dictionary atoms, the rotational transformation matrix, the phase contamination and the DW images. The sparsity of the vectors containing the coefficients individuating the dictionary atoms is leveraged to regularize the reconstruction. In addition, the corresponding DW images are recovered by promoting their sparsity in the wavelet domain.
As the method presented in this paper, ref. [26] focuses on the recovery of the FOD coefficients. In this work, the sparsity of the FOD coefficients is promoted through the 1 norm. Furthermore, a total variation (TV) penalty, acting on the DW images, indirectly promotes spatial fiber regularity within neighbor voxels.
In the present work, we develop a method to recover intra-voxel fiber configurations (FODs) at high spatio-angular resolution relying on a 3D kq-space under-sampling scheme. For each of a reduced number of q-space points, the time available for k-space sampling is used to acquire a sub-sampled k-space of a 3D sub-volume, rather than the full k-space of a 2D slice, thus, providing a potentially significant acceleration over the state-of-the-art.
Anatomical constraints are leveraged in order to recover the FOD coefficients. First, prior knowledge of the spatial distribution of the brain tissues, which can be inferred from the image acquired in the absence of diffusion, is explicitly incorporated in the FOD recovery problem. Secondly, the FOD estimate is defined as a solution to a regularized minimization problem, using the structured sparsity prior proposed in [24]. This regularization promotes simultaneously voxel-wise sparsity and spatial smoothness of fiber orientation.
Following [22,24], the resulting non-convex minimization problem is solved via a re-weighting scheme [29] involving a sequence of convex minimization problems with weighted sparsity prior. Specific to our contribution, a stochastic Forward-Backward (FB) algorithm is used to solve each of these convex problems, with convergence guarantees toward the minimum of the corresponding convex objective [42]. One of the advantages of the stochastic approach is that it offers the possibility to handle efficiently multi-coil acquisitions involving large volumes of data, minimizing both the memory requirement per iteration and the reconstruction time.
Results from both simulated and real data experiments showed that the proposed approach outperformed, in terms of its potential of accelerating the acquisition while preserving the FOD reconstruction quality, both the existing method in kq-space [26] and the traditional reconstruction approaches considering the TV prior for the reconstruction of the DW signals before the recovery of the FOD coefficients. Importantly, our experiments reveal that the optimal reduction of samples was achieved, not by simply reducing the number of q-space points while fully sampling the k-space, but rather by a combined kq-space under-sampling.
The remainder of the paper is organized as follows. The proposed approach is described in Section 2. In particular, the FOD reconstruction inverse problem is described in Section 2.2, and the proposed minimization problem and algorithm are given in Section 2.3. In Section 3, we provide the description of the experimental setup, and, in Section 4, we present the results obtained from both synthetic and real data. Finally, we discuss the achieved reconstruction performances, and we conclude in Section 5.

k and q Spaces Overview
The dMRI signal is captured by a multitude of MRI images, each of which is sensitive to water molecule diffusion occurring at a specific angle. As for the standard MRI, the dMRI signal is spatially encoded by the application of three spatial encoding gradients spanning the 3D space. The k-space corresponds with good approximation to the 2D or 3D Fourier transform of the signal. In addition to the application of the spatial encoding gradients, dMRI signals must be encoded for the diffusion process.
To this aim, additional gradients are applied during the acquisition, that are the diffusion gradient. Similarly to the k-space, encoded by the spatial gradients, the q-space is the 3D space defined by the diffusion gradients applied during the dMRI acquisition process. Each dMRI volume is associated with diffusion along a specific directionq and, at a specific intensity |q|, characterizing a point q = |q|q. The commonly used b-value parameter [10] associated with each DW volume is defined as b = 4π 2 |q| 4 t, with t as the diffusion time.

FOD Recovery via SD Framework
In this section, we describe the SD framework for FOD recovery, from q-space under-sampled measurements. We refer to [24] for further details of the global FOD recovery problem.
Let X ∈ R + (n+2)×N be the unknown matrix containing the FOD coefficients of interest, where N is the number of imaged voxels and n is the number of dictionary atoms. Each column of X contains the n + 2 coefficients of the corresponding voxel. The objective is to find an estimate of the FOD field associated with all the voxels of the brain, from the degraded measurements Y ∈ R M×N , given by Y = ΦX + η, where Φ ∈ R M×(n+2) is the observation matrix and η ∈ R M×N models the acquisition noise. Note that, for SNR > 2 the acquisition noise consists in a realization of an additive random i.i.d. Gaussian noise [43], where the Signal-to-Noise Ratio on the image acquired in absence of diffusion gradients (i.e., s 0 image) is defined as SNR = s 0 σ , where · denotes the mean of its argument and σ is the standard deviation of the noise.
More precisely, Φ is a known dictionary whose columns (called dictionary atoms) correspond to the response signals on M q-points to n single fiber orientation [19]. In addition, two atoms representing isotropic diffusion (typically in gray matter or CSF) are considered in the dictionary. The corresponding atoms are invariant under rotation in q-space. The procedure adopted to generate Φ is provided in Section 3.3, based on the method developed in [22]. In particular, the first row of Φ is dedicated to atom values in the absence of diffusion gradient, equal to 1: The rows of the data matrix Y ∈ R M×N span the unfolded diffusion volumes, acquired with M gradients and normalized by the intensities of the volumes acquired in the absence of diffusion, denoted by s 0 . The s 0 volume is assumed to be obtained from a separate acquisition and is used here as prior information. The first row of Y is devoted to the normalized s 0 volume, i.e., all its coefficients are equal to 1. Thus, for each voxel v ∈ {1, . . . , N}, In summary, introducing unitary lines in both Y and Φ intrinsically forces the sum of the FOD coefficients to be equal to 1, injecting directly this prior information into the inverse problem.

Proposed Measurement Model
In this section, we describe the proposed under-sampled kq-space measurement model. We consider acquisitions obtained from C coil receivers, with M diffusion gradients, sampling K k-space points. The measurement matrixŶ ∈ C MC×K is expressed as follows: whereŶ q,c ∈ C 1×K is the under-sampled k-space of the diffusion volume acquired with gradient indexed by q and by the receiver coil indexed by c, and η q,c ∈ C 1×K is a realization of an independent and identically distributed (i.i.d) Gaussian noise. Indeed, while the noise contamination on the magnitude signals is characterized by a Rician distribution, the original noise in k-space is Gaussian [44]. The matrix X ∈ R (n+2)×N represents the unknown FOD field of interest, and the linear operator A q,c is given by: The various terms defining A q,c are described as follows. Φ q ∈ R 1×(n+2) is the qth row of the dictionary Φ that spans the response of a single fiber oriented along n different directions, to which two isotropic compartments are added, to represent the gray matter and the CSF. In order to implicitly force the FOD coefficients of each voxel to sum up to one, the Fourier coefficients of the s 0 volume are intentionally introduced in the first row ofŶ. Here again, the s 0 volume is assumed to be obtained from a separate acquisition and is used as prior information.
Measurements cannot be normalized as the diffusion volumes themselves are not accessible but only an incomplete k-space counterpart. Thus, a matrix S 0 ∈ R N×N , whose diagonal elements correspond to the s 0 volume coefficients, is explicitly introduced in the full measurement operator. The acquisition of the diffusion signal from multiple channels is taken into account through the diagonal matrix U (c) ∈ C N×N , which contains the sensitivity map of the corresponding receiver coil c. Moreover, motion and magnetic field inhomogeneities generate phase distortions that are taken into account in the diagonal matrix H (q,c) ∈ C N×N . We refer to Sections 3.5 and 3.4 for the description of the procedures used to estimate the coil sensitivity and the phase contamination, respectively. Finally, F ∈ C N×N represents the 3D Fourier matrix and M (q) ∈ R N×K q is a realization of a binary mask that undersamples each slice of the acquired volume. It is important to notice that different realizations of M (q) may be considered for each applied diffusion gradient. In addition, M (1) is chosen to be the identity matrix of R N , to fully acquire the Fourier coefficients of the s 0 signal (stored in Y 1,c ) for normalization, segmentation and calibration purposes.

Tissue Segmentation Constraints
The spatial distribution maps of the different brain tissues can be obtained from the segmentation of a structural image of the imaged brain. Such images are generally acquired by default in clinical practice; alternatively, the spatial distribution maps can be estimated from the s 0 signal, which is always fully acquired for normalization and calibration purposes. An illustration of a tissue distribution is given in the first row of Figure 1. In our approach, we propose to explicitly take advantage of the prior knowledge of tissue distributions over the space in order to further regularize the ill-posed FOD estimation problem. We assume that the number of voxels containing white matter, gray matter and CSF tissues, are known and indicated by N 1 , N 2 and N 3 , respectively.
In order to fully exploit the tissue distribution information, we propose a novel formulation for the FOD global problem in the kq-space. To this aim, we define three different variables S 1 ∈ R + n×N 1 , S 2 ∈ R + 1×N 2 and S 3 ∈ R + 1×N 3 , related to the white matter tissue, the gray matter tissue and the CSF, respectively, as shown in Figure 1.
More specifically, S 1 contains the effective FOD coefficients associated with the white matter fiber, S 2 models the isotropic behavior characterizing the gray matter tissue, and S 3 takes into account the isotropic behavior of the CSF. The object X is fully characterized by S = (S 1 , S 2 , S 3 ) through the linear mapping X = Z (S), where the operator Z concatenates the matrices resulting from the expansions of S 1 , S 2 and S 3 with zero-valued columns in the places of the voxels that are known to not contain the corresponding tissue.
A schematic representation of X = Z (S) is reported in Figure 1. In accordance with the white matter distribution maps (blue map on the top of Figure 1), the columns of S 1 (blue columns on the bottom) span the FOD coefficients of the voxels that are characterized by nonzero values in the map. In an analogous way, the coefficients of S 2 and S 3 (green and yellow columns on the bottom of Figure 1, respectively) represent the isotropic compartments associated with the voxels corresponding to the nonzero values in the gray matter and CSF maps (green and yellow map on the top of Figure 1, respectively). Using the proposed notation, problem (1) can be rewritten as: It is worth noticing that, in the particular case when the global FOD estimation formulation is applied only to white matter voxels, we have X = S = S 1 , and Z is the identity operator. Each column of S 1 spans the FOD coefficients of the voxels containing white matter tissue (blue columns). Coefficients modeling gray matter and CSF are stored in S 2 (green columns) and in S 3 (yellow columns), respectively. The FOD reconstruction is performed by estimating S 1 , S 2 and S 3 , whose active columns are provided by the tissue segmentation maps (top).

Constrained Weighted-1 Minimization Problem
Recently, convex and non-convex optimization methods have gained a great deal of attention to solve ill-posed or ill-conditioned inverse problems. In this context, the estimate of the sought object is defined to be a minimizer of a sum of two functions: the datafidelity term associated with the signal model, and the regularization term incorporating a priori information of the target solution. On the one hand, when having an additive i.i.d. Gaussian noise, the data-fidelity of choice is the least squares criterion. On the other hand, sparsity-aware models, especially leveraging a CS framework [45], are particularly suitable for solving ill-posed or ill-conditioned inverse problems.
The most suitable way to promote sparsity consists in using the 0 pseudo-norm, which counts for the number of non-zero coefficients of the signal of interest [46]. However, the 0 pseudo-norm, being neither smooth nor convex, can be difficult to handle efficiently and it is often replaced by its convex relaxation, namely the 1 norm [45]. In the context of FOD reconstruction, the signal X is sparse in the q-space as a consequence of the low number of fiber populations (i.e., non-zero FOD coefficients) that are expected to be contained within a voxel.
However, when SD problems are considered, it was demonstrated in [22] that the use of the 1 norm is inconsistent with the physical constraint that the volume fractions of each voxel sum up to one. To overcome this difficulty, [22] proposed to use a reweighting 1 approach [29] to approximate the 0 pseudo-norm in the context of the voxel-wise FOD estimation. Such an approach was subsequently generalized in [24] to solve the inverse problem described in Section 2.1.2. This consists of solving sequentially several constrained weighted 1 problems. The weights are chosen to simultaneously promote voxel-wise sparsity and spatial smoothness of fiber orientation. Formally, the problem is written as follows: min where B + 1,W (κ) denotes the intersection of the real positive orthant R (n+2)×N + with the weighted 1 ball of radius κ 0 centered in 0 with weighting matrix · 1,W is the weighted 1 norm given by In this work, we generalize this weighting scheme to solve the kq-space inverse problem (3). In particular, we aim to solve where A : R (n+2)×N → C MC×k consists of the concatenation of the operators A q,c , for all (q, c).

Algorithm for Constrained Weighted 1 Minimization
As proposed by the reweighting framework, problem (7) is solved several times, considering different weighting matrices W, in order to mimic the 0 pseudo-norm. Note that the reweighting framework is underpinned by recent theoretical results [47][48][49][50] showing that sequentially minimizing convex problems with weighted-1 priors corresponds to solving for a critical point of a non-convex problem with a log-sum prior, itself a close approximation to the target 0 prior. Considering one weighting cycle, a simple algorithm to solve problem (7) is the FB algorithm [51].
by performing a gradient step on the data-fidelity term, and three projection steps to handle the three constraints described in (7). However, computing the gradient step can be computationally unaffordable when processing high dimensional data, which is particularly the case with multi-coil kq-space raw data.
To overcome this issue, we rely on a more evolved algorithm, namely the stochastic FB algorithm [42], where only an approximation of the gradient of interest is computed at each iteration. Specifically, we propose to select randomly at each iteration a subset of the available coils, and approximate the gradient using only the data from the selected coils. The proposed method is described in Algorithm 1.
Steps 6-10 describe the computation of the approximated gradient. Below is an explanation of these steps. LetŶ c be the data measured by coil c and A c be its associated measurement operator. At iteration j ∈ N, the gradient of the data fidelity term A Z (S (j) ) −Ŷ 2 2 , denoted by G (j+1) , can be written as the sum Using the stochastic approach described in [42], the gradient G (j) can be approximated by only updating a subset S c . More specifically, we define the approximation G (j) of G (j) as follows: ; otherwise, note that (·) † denotes the adjoint operator of its argument.

Algorithm 1 Stochastic FB to solve
end for 10: 11: Steps 12-14 describe the projection steps. The projection of a variable x ∈ R N onto a nonempty, closed, convex set C ⊂ R N is given by P C (x) = min y∈C x − y 2 2 . This corresponds to the closest point of x belonging to C, using an Euclidean distance [51].
Step 12 is the projection onto B + 1,W (κ), related to the white matter. In accordance with the 0 approximation, κ is chosen to represent the S 1 sparsity level. Steps 13 and 14 are the projections related to the gray matter and the CSF, respectively, performed onto R 1×N 2 + and R 1×N 3 + . Let (S (j) ) j∈N be a sequence generated by Algorithm 1. According to [51], if 0 < γ < 2/ A 2 S , where · S denotes the spectral norm, then (S (j) ) j∈N converges to a solution S of problem (7). In our simulations, we consider that Algorithm 1 has converged when S (j+1) − S (j) 2 < ν S (j) 2 , where ν = 10 −3 is a tolerance, fixed by the user.

Weights Computation and Reweighting Procedure
We propose to use an 1 reweighting procedure to approximate the 0 pseudo-norm [29]. To this aim, the matrix of weights W in Step 12 of Algorithm 1, is chosen using a similar approach as described in [24]. The reweighting method is described in Algorithm 2.
2: Let W (0) = 1 n×N 1 . 3: Iterations: 4: For t = 0, 1, . . . , T − 1 5: : compute W (t+1) as per Equations (8) and (9)  7: end for At each cycle t ∈ {0, . . . , T − 1} of the reweighting scheme, problem (7) is solved using Algorithm 1 for a particular choice of the weighting matrix (see Step 5). The weighting matrix W (t+1) = (W (t+1) ) 1≤d≤n+2,1≤v≤N 1 is updated in Step 6 of Algorithm 2, as follows: where and with τ > 0. The parameter τ (t) is a stability parameter that avoids the weights going to infinity when B in such a way that large weights progressively force to zero spurious peaks, while small weights favor the presence of the FOD coefficients. Consequently the weighting matrix, associated with the 1 norm, promotes sparsity. In addition, the smoothness of the variation of the FOD coefficients across spatial neighborhoods is also promoted in (9) by averaging over neighbor voxels. Due to the discretization of the dictionary, fiber contributions are usually spread over a small angular support.
That is the reason why (9) contains a summation (rather than averaging) over neighbor directions. If a spatio-angular neighborhood is reasonably homogeneous, the weight in (8) will be reasonably small and the corresponding FOD coefficient will be preserved at the next iteration. On the contrary, spurious peaks isolated in their spatio-angular neighborhood will be associated with large weights. These large weights will tend to prevent such isolated peaks at the next iteration to not violate the weighted 1 -norm constraint.
Finally, the reweighting process stops when the number of allowed iterations T is reached or there is no substantial relative variation between successive estimates of S 1 , i.e., S

Data Post-Processing
Once the solution is found, a post-processing procedure is performed along the columns of S 1 in order to extract the directions of the fibers within each voxel. The identification of the highest peaks is performed among all the directions contained within a cone of 30 • for each different direction. No more than eight local peaks are assumed per voxel, and peaks smaller than 20% of the maxima are disregarded in order to suppress spurious contributions [22,24].

Experimental Setting
In this section, we provide the description of the different experimental settings. In particular, the specifications of the sampling schemes considered in the q-and in the kspaces are provided in Sections 3.1 and 3.2, respectively. The method proposed to generate the dictionary is described in Section 3.3. Sections 3.4 and 3.5 provide the procedures used to estimate the phase and the coil sensitivities, respectively. A description of the metric used to quantitatively estimate the FOD reconstruction is provided in Section 3.6.

q-Space Under-Sampling
In this section, we describe the scheme adopted for sampling the signal in q-space. In the case of the synthetic data, the diffusion signal is uniformly sampled over a single shell using 30 diffusion gradients with b = 1000 s/mm 2 . In order to assess the FOD recovery in the presence of various kq-space regimes, the data provided was retrospectively under-sampled by selecting 6, 10, 15 and 20 uniformly distributed q-points, as reported in Figure 2. In the case of the real data, the signal in q-space is sampled by 60 gradients over eight shells, considering b-values going from 550 to 4800 s/mm 2 . In particular, 3, 6, 4, 3, 12, 12, 6 and 13 q-points are considered on sphere 1 (b = 550 s/mm 2 ), sphere 2 (b = 1050 s/mm 2 ), sphere 3 (b = 1600 s/mm 2 ), sphere 4 (b = 2150 s/mm 2 ), sphere 5 (b = 2650 s/mm 2 ), sphere 6 (b = 3200 s/mm 2 ), sphere 7 (b = 4250 s/mm 2 ) and sphere 8 (b = 4800 s/mm 2 ), respectively.
Three under-sampled datasets were created by retrospectively selecting the signal associated with 15, 30 and 45 q-points in order to compare the FOD recovery in the presence of different q-space under-sampling regimes. Each dataset was created by retrospectively sampling q-points from different shells following a Gaussian distribution centered in the fifth shell, that is to say, q-points with b-values that ranges between 2000 s/mm 2 and 3000 s/mm 2 are sampled with higher probability [52]. Five different realizations were created for each q-space setting.
For all the synthetic and real datasets, the acquisition of the signal in the absence of diffusion is required in order to obtain s 0 . Examples of this image, for both the synthetic and real data, are provided in Figure 3A,D, respectively.

k-Space Under-Sampling
As already explained, the objective of this work is to provide a robust reconstruction method when considering both q-and k-space under-sampling. The description of the schemes used to sample the signal in k-space is provided in this section.
In both synthetic and real data, for each gradient q ∈ {1, . . . , M}, selection masks M (q) are built such that the signal in k-space is sampled within a Cartesian grid by a continuous trajectory. We call k-space under-sampling factor the ratio between the total amount of grid points in k-space, for a given field of view and resolution, and the number of samples considered for the reconstruction.
The potential of the proposed method is fully exploited in the 3D setting where the acquisition time can be significantly reduced by limiting both the number of excitations to cover each DW volume and the number of DW sub-volumes. As a proof of concept, the proposed FOD recovery method is evaluated for a specific class of 3D Fourier sampling patterns where the same 2D k-space under-sampling pattern is probed at all spatial frequencies in the third dimension of a selected sub-volume. By a simple inverse Fourier transform in this third dimension, the corresponding data can be cast, without loss of generality, as a slice-by-slice-identical 2D k-space under-sampling. The scheme chosen for the 2D k-space sampling consists in fully sampling the central lines of the k-space while regularly skipping lines at the periphery, with the central zone of the k-space is required to estimate the phase affecting each shot (i.e., each sub-volume).
We argue that a two-phase EPI acquisition protocol could be used to that effect. EPI commonly considers regularly sampled lines to enable a simple calibration of the eddy currents affecting each different k-space line. The proposed k-space under-sampling scheme could result from the combination of two uniform EPI schemes, one for the center and one for the periphery, with the eddy currents calibration critical only at the interface between the two regions. A schematic representation of the under-sampling scheme used with both synthetic and real data is reported in Figure 4. An experimental evaluation of a precise acquisition protocol is beyond the scope of the present work. Figure 4. Illustration of the acquisition of the 3D volume containing the brain. Each sub-volume is acquired considering a 3D two-phase EPI scheme (bottom line). The presented results were obtained considering a sequence of identical slice-by-slice 2D k-space under-sampling schemes (top line), which can be brought back to a specific class of 3D k-space under-sampling scheme to which a Fourier transform in the third dimension was applied. The proposed pattern results from the combination of two different regular schemes: the fully sampled central zone (yellow lines) and the regularly under-sampled periphery (white lines).
In the remainder of this work, all our experiments are designed with this k-space sampling scheme. Note, again without loss of generality, that the performance of the proposed approach is evaluated in a setting including a single DW (sub-)volume at each q-space point. Multiple sub-volumes would simply be acquired sequentially and reconstructed separately.
Note that the proposed method requires the complete k-space acquisition of the s 0 signal to implicitly force the FOD coefficients to sum up to one in each voxel. Furthermore, s 0 is used for normalization, segmentation and coil sensitivities calibration purposes as described in detail in Sections 3.4 and 3.5.

Dictionary Generation
In this section, we provide the description of the method used to create Φ for both the simulated and the real data. The elements of the dictionary Φ are generated by relying on the Gaussian Mixture Model of the q-space signal [8]. More specifically, the dictionary where q is the index that explores the diffusion gradients q (with b-value b and orientationq) considered for the acquisition, and d is the index associated with the orientations chosen for the discretization of half of the unit sphere. We consider n = 500 points for the discretization of the dictionary. Moreover, two additional atoms are considered to model the isotropic diffusion in the gray matter and CSF. The d-th atom of the dictionary is indicated by Φ q,d 1≤q≤M , and this corresponds to the response of a single fiber, oriented along the d-th direction, subject to M different diffusion gradients. The diffusion tensor D = diag(λ 1 , λ 2 , λ 3 ) characterizes each fiber population, where λ 1 is the longitudinal diffusivity, and λ 2 and λ 3 are the 2 transverse diffusivity coefficients. In particular, D (d) is the rotated version of D along the direction d.

Phase Estimation
In theory, diffusion images are real-valued. However, in practice, they are often contaminated by phase factors during the diffusion encoding process. This phase contamination is mostly due to magnetic field inhomogeneities and biological motion (e.g., physiological and involuntary patient motion).
Usually, methods dealing with signals directly in q-space overcome this difficulty by simply taking the modulus of the complex diffusion signal, in order to obtain real diffusionweighted images. However this method cannot be used in k-space and kq-space. Indeed, since the phase contamination breaks the Hermitian symmetry of the images, it cannot trivially be removed. In particular, in the linear operator expressed in (2), the diagonal matrix H (q,c) takes into account this phase factor, for a fixed gradient q and coil receiver c.
In the case of synthetic data, the phase contamination produced by the magnetic field inhomogeneities is provided within the contest data set. In addition, effects due to motion are mimicked by adding a different linear phase map to the signal associated with each sub-volume of the phantom and coil receiver. The linear phase maps are generated in such a way that the corresponding k-space shift is constrained within 10 phase-encoding lines. In each voxel v ∈ {1, . . . , N}, the intensity of the signal in q-space is multiplied by e iθ v , θ v as the sum of the phases due to both motion and magnetic field inhomogeneities associated with v. In Figure 3C, an example of estimated phase map due to motion and magnetic field inhomogeneities is provided.
For both synthetic and real data, phase-distortions are estimated from the central portion of the k-space data [53,54], which is always fully sampled, (see Section 3.2). By performing the inverse Fourier transform of a zero-padded version of the fully sampled central part of the k-space, we obtain a complex-valued signal whose phase provides the diagonal elements of H (q,c) .

Coil Sensitivity Estimation
In this section, we provide the method to estimate the coil sensitivity coefficients stored in the diagonal matrix U (c) , appearing in the linear measurements described in Equations (1) and (2) for both synthetic and real data. Images acquired in absence of diffusion from different coil receivers are combined through the sum of squares [55] method in order to obtain a baseline image. Subsequently, the images associated with each receiver coil are divided by the baseline image in order to determine the corresponding sensitivity map.
For synthetic data, the acquisition from four different coil receivers was simulated by considering the toolbox at http://bigwww.epfl.ch/algorithms/mri-reconstruction/ (acessed on 11 August 2021). Examples of the coil sensitivity, when C = 4, are shown in Figure 3B.
Examples of the estimated sensitivity maps obtained in the real data experiment are provided in Figure 3E.

Evaluation Criteria
The quality of the fiber reconstruction was evaluated by using the metrics proposed in [24]. For each single voxel, the fiber reconstruction evaluation takes into account the number of fibers correctly identified and the angular accuracy of the recovered direction.
First, we define the success rate (SR) index representing the proportion of voxels in which the number of fibers is correctly identified. More precisely, when all the estimated fibers fall within a tolerance cone of 30 • around the true fibers, we have a success. The SR index depends on the false positive and negative rates that represent the average over all voxels of the number of overestimated and underestimated fibers per voxel.
Secondly, we define the mean angular error (θ) as the average of the angular errors associated with each true fiber. The angular error is defined as θ = 180 π arccos(|d true · d estimated |), where d true and d estimated are the true fiber direction and the direction of the closest estimated fiber, respectively.

Results
In this section, we present the results obtained when using the proposed method with both synthetic and real data. Specifically, in Sections 4.1 and 4.2 we discuss the results obtained for synthetic and real data, respectively.

Synthetic Data
A first analysis of the FOD reconstruction using the proposed kq-space under-sampling method was performed relying on the fiber configuration of the numerical phantom proposed in the ISMRM Tractography Challenge 2015 [56]. The phantom consists of a volume of N = 90 × 90 × 108 voxels, acquired by using M = 32 diffusion gradients distributed over a single shell (see Section 3.1 for more details).
We consider the acquisition from four coil receivers, where the sensitivity maps are simulated and estimated using the process provided in Section 3.5. Phase contamination due to motion and magnetic field inhomogeneities was generated as described in Section 3.4 and incorporated in the data. The q-space signal was then converted to k-space through the Fourier transform and contaminated with Gaussian noise with zero mean and standard deviation σ = s 0 SNR . Note that SNR = 30 was chosen for the numerical simulations. Lastly, selection masks, built as described in Section 3.2, were applied to the k-space of the diffusion-weighted images to obtain the under-sampled data in kq-space. The dictionary Φ was generated using the procedure described in Section 3.3.
In the following section, the proposed recovery scheme is compared to three different reconstruction frameworks, which represent some of the state-of-the-art approaches for the FOD recovery in the presence of measurements under-sampled both in k and q-spaces separately and in kq-space jointly. Among the state-of-the-art methods investigated, we can distinguish between two types of approaches: the two-step and the one-step approaches.
The two-step approaches are the conventional approach to recover the FOD coefficients in the presence of under-sampled kq-space data, where the DW images are first reconstructed, and the FOD coefficients are subsequently estimated. As already mentioned in Section 3.4, two-step approaches avoid to model the phase contamination due to motion and magnetic field inhomogeneities, by recovering complex images, whose imaginary part is discarded by taking the magnitude. Among all the k-space methods, the TV-regularized approaches were considered for this comparison. The one-step approaches consists on the fiber orientation estimation from the under-sampled kq-space signal directly. Contrary to the two-step approaches, the one-step approaches require the phase contamination to be modeled. In particular, the proposed approach is a one-step method. Below is a summary of the considered methods.

1.
TV − L2 + is a two-step approach consisting in a first step where the DW volumes are recovered by using the TV + regularization and in a second step where FOD coefficients are recovered by relying on the non-negative least squares problem.

STEP 1: min
whereŶ ∈ C MC×N contains the signal in kq-space. Complex DW volumes are recovered in I ∈ C M×N from the TV regularized problem, where F (I) = (F q,c (I)), with F q,c (I) = IU (c) FM (q) for q ∈ {1, . . . , M} and c ∈ {1, . . . , C}. To exclude the remaining phase contamination, the magnitude of the complex images is considered for the FOD recovery (i.e.,Ĩ = |I| S 0 ).

STEP 2: min
where D q (X) = Φ q X for all q ∈ {1, . . . , M}. Positivity of the FOD coefficients is imposed in this step.

2.
TV − STR + is a two-step approach consisting in the recovery of both DW volumes and FOD coefficients from regularized problems. In the first step, the TV prior is considered while, in the second step, the structured sparsity prior proposed in [24] is exploited.
TV + is a one-step approach consisting in the recovery of the FOD coefficients from the kq-space signal applying the TV prior to the images in order to implicitly promote a smooth variation of the FOD coefficients within neighbor voxels.

4.
STR + − TISS is the proposed one-step approach consisting in the recovery of the FOD coefficients by solving the problem proposed in (7) where the structured sparsity prior and the spatial distribution of the different tissues is taken into account.
The Primal-Dual [57] and the FB algorithms were used to, respectively, solve the first and the second step of the two-step approaches. More specifically, in the second step of TV − STR + , the FB algorithm was used multiple times to solve 10 weighted-1 problems. For all the frameworks, the optimal parameters were found using a gridsearch approach. The evaluation of the image reconstruction in the case of the two-step approaches was performed relying on the SNR index. On the other hand, the SR index was chosen to evaluate the quality of the fiber orientation reconstruction and determine the best parameters for both the second step of the two-step approaches and the onestep approaches.
The parameters obtained in such a way are as follows: 1 × 10 +3 ≤ λ ≤ 2.1 × 10 +3 and 3.9N ≤ κ ≤ 4.2N. TV + is solved using the Primal-Dual algorithm, considering 4 × 10 +5 ≤ λ ≤ 9.8 × 10 +5 . On the other hand, STR + − TISS was solved using the Algorithm 2, with κ = 4N wm , with N wm as the number of white matter voxels. However, in the case of the synthetic data, where only four coil receivers are considered, the gradient of the data fidelity term is fully computed and not approximated. Hence, at each iteration j ∈ N, we choose S (j) c = {1, 2, 3, 4} in Algorithm 1, which reduces to the FB algorithm. We emphasize that additional tests were performed for TV − L2 + and TV − STR + with λ = 0, corresponding to remove the TV regularization for the DW images in the least-squares minimization problem. These tests have shown that the use of the TV prior provides better results.
In Figure 5, the quantitative evaluation of the fiber reconstruction is provided in the case of various kq-space under-sampling regimes for the four different methods described above. Performances obtained considering different q-space under-sampling ratios are reported in different colors, while performances obtained considering different k-space under-sampling factors are reported along the x-axis. The quality of the reconstruction was evaluated considering noisy data with SNR = 30, and 5 different q-space undersampling factors (i.e., M ∈ {6, 10, 15, 20, 30}), and 6 different k-space under-sampling factors (i.e., N/K ∈ {1,2,4,6,8,10}).
The obtained SR index, mean angular error, and the rate of false positives and negatives are shown in Figure 5, from the first to the last row, respectively. The error bars correspond to the standard deviation of the evaluation indices obtained from k-space data corrupted by different noise realisations. Performances obtained with the two-step approaches are reported in the graphs in the first two columns of Figure 5. In the presence of k-space under-sampling regimes, both TV − L2 + and TV − STR + are characterized by SR rates lower that the proposed STR + − TISS method. For TV − L2 + , the low performances result from the high number of spurious peaks that is estimated in each voxel as can be observed in the diagram in the first column of Figure 5. When using the structured sparsity prior to regularize the FOD coefficients in TV − STR + the number of false positive drastically decreases, as shown in the corresponding graph in the second column of Figure 5. However, the number of false negative moderately rises up, preventing the SR index to increase. In general, the direct regularization of the FOD provides better performances than TV − L2 + .
The last two columns of Figure 5 give the graphs of the two one-step approaches: the TV + method (third column), and the proposed method STR + − TISS (fourth column). TV + does not reach high performances for any of the kq-space under-sampling regime. The SR index remains below 0.31, that is reached considering 30 q-points with no k-space under-sampling. These low SR values are due to the fact the TV regularization is incapable of excluding the numerous spurious peaks that affect the FOD reconstruction.
However, the performances obtained with the one-step approach TV + , are shown to be quite robust to the k-space under-sampling. Indeed, the SR decreases only of 0.15 (resp. 0.15) and θ increases of 4 • (resp. 5 • ) going from no k-space under-sampling to a k-space under-sampling factor of 10 with 30 q-points (resp. 10 q-points). As discussed in Section 2.3.1, and shown in [22], since fiber coefficients are implicitly forced to sum up to one in each voxel, imposing 1 norm in this context is ineffective. Additional simulations confirmed this behavior showing that using 1 norm to regularize the FOD coefficients in addition to TV did not lead to a significant difference in the results.
Concerning the proposed STR + − TISS method, STR + − TISS outperforms both the two-step and the one-step approaches here above presented. When compared to TV + , STR + − TISS shows to substantially benefit from the direct regularization of the FOD coefficients, through the structured sparsity prior. Considering no k-space under-sampling with 6-points (resp. 30 q-points), SR = 0.31 for TV + (resp. 0.48), against SR = 0.84 for STR + − TISS (resp. 0.86). For a k-space under-sampling factor of 10 with 6-points (resp. 30 q-points), SR = 0.06 for TV + (resp. 0.34), against SR = 0.62 for STR + − TISS (resp. 0.75). Thus STR + − TISS exhibits much higher performances when compared to TV + . However, both methods are characterized by a similar amount of false negative.
For the sake of brevity, we do not provide the comparison considering or ignoring the tissue spatial distribution. Additional simulations have shown that, when the regularization parameter κ is known, FOD configurations recovered by the two methods are comparable. However, in practice κ needs to be tuned, and using prior knowledge of the tissue distribution improves the reconstruction. This fact can be understood as follows. In the presence of under-estimated κ, coefficients associated with isotropic tissues compartments take precedence over the fiber compartments (since isotropic behavior can more easily fit the data). When the tissue distribution is taken into account, the coefficients associated with the fiber compartments are processed separately thus preventing them from being annihilated by the dominant presence of the isotropic compartment coefficients.
In Figure 6, we provide an illustrative example, considering 15 diffusion gradients with a k-space under-sampling factor of 4, of the performance obtained from the four considered methods. In particular, we show the mean angular error (first row), the SR (second row), the false positive (third row) and the false negative (fourth row) maps. These results are confirming the quantitative results described in Figure 5. In particular, the maps evaluating the fiber configurations recovered from TV − L2 + and TV − STR + show the worst performances. The SR maps in Figure 6(1B,2B) mostly show red pixels, indicating the presence of voxels with an uncorrected number of estimated fibers. For the TV − L2 + method, the majority of such voxels is characterized by over-estimated fiber populations, as it is evidenced by the map in Figure 6(1C).
By comparing the maps in Figure 6(1C,2C), we can observe that the use of the structured sparsity prior significantly decreases the amount of false positive. Nevertheless, the map in Figure 6(2D) highlights the presence of higher false negatives rates. For the TV + method, the SR remains low, as it can be observed from the large areas with red pixels in Figure 6(3B). However, the angular accuracy of the estimated fiber configurations is significantly improved when compared to the two-step approaches (see maps on the first row of Figure 6). In Figure 6(4B), the SR map obtained with the STR − TISS + method shows the higher rate of yellow pixels indicating the presence of a large number of voxels with the correct number of estimated fibers. In addition, the angular accuracy achieved with STR − TISS + remains low, as displayed by the map in Figure 6(4A).

Real Data
The kq-space under-sampling framework was tested for the recovery of the fiber orientation from real data. The real data set consists of a volume of 100 × 100 × 60 voxels acquired on a 3T Magnetom Trio system (Siemens, Germany) with a pixel size of 2.2 × 2.2 × 2.2 mm 3 . The dataset was acquired considering 17 different coil receivers, and 60 diffusion gradients distributed over multiple shells, as described in Section 3.1. The complete signal was retrospectively under-sampled both in q-and k-space (as described in Sections 3.1 and 3.2, respectively) in order to compare the reconstruction performances considering different under-sampling regimes.
In the case of real data, Φ was generated by considering n = 500 for the discretization of the sphere using the procedure described in Section 3.3. The tissue segmentation maps were estimated from the segmentation of the s 0 signal and H (q,c) was estimated from the low resolution image as described in Section 3.4. The diagonal of the matrix S 0 is filled with the s 0 values, which were recovered from the multi-coil signals while correcting for the estimated phase through a least squares approach. Ultimately, U (c) was estimated from the s 0 signal as described in Section 3.5. Note that the real data was corrected for the EPI Nyquist ghosting prior to the reconstruction.
Results were obtained by solving problem (7) using Algorithm 2, where the reweighing process is performed 10 times at most, and κ = 5 × N wm . The computation of the gradient in Algorithm 1 is approximated choosing 12 out of 17 coils at each iteration, i.e., |S c | = 12. In particular, four of these coils are fixed (chosen at each iteration), and eight are randomly selected. In addition, a Nesterov acceleration [58] is considered in Algorithm 1. Although no theoretical result ensures the convergence of the resulting method, additional simulations have shown that Nesterov acceleration leads to much faster convergence while the same minimum is achieved using both the methods (with or without Nesterov acceleration).
Considering the fiber configuration recovered from 60 diffusion gradients with full k-space as reference, quantitative evaluation maps are computed and provided in Figure 7, for the fiber configurations obtained from four different kq-space under-sampling settings. For each different setting, the FOD reconstruction was performed 5 times, each considering a different realisation of q-space under-sampling. The results presented here corresponds to the mean and the standard deviation of the evaluation indices obtained considering different q-space under-sampling realizations.
In the first column of Figure 7 the performances obtained considering 15 diffusion gradients with full k-space is reported. The performances achieved with 30 diffusion gradients considering a k-space under-sampling factor of 2 and 45 diffusion gradients considering a k-space under-sampling factor of 3 are shown in the second and third columns of Figure 7, respectively. Ultimately, the evaluation maps of the fiber geometries recovered from 60 diffusion gradients with k-space under-sampling factor of 4 are shown in the fourth column of Figure 7.
The comparison between the evaluation maps provided in Figure 7 shows that, considering the same overall kq-space under-sampling factor, the recovery from data that were under-sampled only in q-space results in lower quality reconstruction (SR = 0.18 ± 0.06, θ = 24.5 • ± 5.1), as highlighted by the larger amount of red pixels and extended areas with high angular error in Figure 7A,E, respectively. On the other hand, the fiber configuration recovered from 60 diffusion gradients with 4-fold k-space under-sampling exhibits the highest reconstruction quality (SR = 0.59 ± 0.01,θ = 7.5 • ± 0.2) as can be observed in the fourth column of Figure 7.
In Figure 8, we provide a zoomed view of the fiber configurations recovered from the real data in the presence of the same kq-space under-sampling settings described above. The fiber configuration obtained considering 60 diffusion gradients with full k-space is provided in Figure 8A as golden reference. Fiber geometries reported in Figure 8B-D appear to be very close to the golden reference. In contrast, the fiber configuration estimated in Figure 8E fails to represent the main fiber crossing revealed in the golden reference.
The analysis of the quantitative and qualitative results strongly suggests that strategies with combined kq-space sampling are advisable, when compared to q-space only strategies. Figure 7. Quantitative evaluation of the fiber orientation reconstruction performed on the real data by using the proposed kq-space under-sampling framework. Reconstructions considering 15 diffusion gradients with full k-space (first column), 30 diffusion gradients with k-space under-sampling factor of 2 (second column), 45 diffusion gradients with k-space under-sampling factor of 3 (third column) and 60 gradient with k-space under-sampling factor of 4 (fourth column) were compared with the fiber configuration obtained using 60 q-points and full k-space. The maps of the mean angular error (first row), success rate (second row), false positives (third row) and negatives (fourth row) are considered for the evaluation.
In Table 1, we compare the computational time required to process the real data with or without stochastic approach. Specifically, we investigate the computational time per iteration either when all the 17 coils are considered at each iteration (deterministic approach), or when 12 of the 17 coils are selected randomly at each iteration (stochastic approach). On the one hand, we can observe that the computational time per iteration when considering 12 coils per iteration corresponds (approximately) to 12/17 times the computational time per iteration when considering 17 coils. This suggests that the computation of the gradient step in Algorithm 1 is much heavier than the computation of the projection step.  Thus, the computational time per iteration is lower using the stochastic approach than the deterministic approach. On the other hand, it can be noticed that both the deterministic and the stochastic approaches necessitate approximately the same number of iterations to reach the stopping criterion. Consequently, the use of a stochastic approach to recover the fiber configurations considering 12 coils per iteration enables to reduce both the memory requirement and the total computational time.

Discussion
We developed a method to accelerate high angular and spatial resolution dMRI acquisition relying on a 3D kq-space under-sampling scheme. We provided a novel formulation to estimate the FOD coefficients from highly under-sampled sub-volumes. The proposed approach used two types of anatomical priors. On the one hand, a brain tissue segmentation constraint, which does not required additional acquisitions and can be inferred from the no-diffusion-weighted image, was explicitly imposed for the FOD recovery. On the other hand, structured sparsity, developed previously in [24], was leveraged to promote simultaneously voxel-wise sparsity and spatial smoothness of fiber orientation.
The resulting minimization problem was approached via a reweighting scheme to solve a sequence of convex minimization problems, using a stochastic FB algorithm structure. By considering the stochastic variation, multi-coil data can be processed while minimizing both the memory requirement and the reconstruction time. Through synthetic and real data experiments, we demonstrated that the proposed recovery framework outperformed the existing method in kq-space and the traditional two-step reconstruction approaches recovering sequentially the diffusion weighted images and the FODs. Furthermore, we observed that, for equal overall under-sampling ratios, the proposed kq-space approach performed better when the k-space under-sampling was exploited rather than heavily under-sampling in q-space only.
Some limitations of our work are to be highlighted. First, a two-phase uniform EPI scheme was proposed to simultaneously fully sample the center and under-sample the periphery of the k-space data in order to minimize the effects arising from the irregular sampling. In this context, the use of more advanced models, taking into account the correction of the geometrical distortions typically affecting the EPI data and the use of acquisition schemes less prone to artifacts, would certainly contribute to improving the performance of the proposed method.
Secondly, one of the critical steps of the proposed kq-space under-sampling method is the motion-induced phase estimation. We considered the phase estimated from the low resolution image associated with a fully sampled central k-space region. However, the use of more sophisticated methods to calibrate the motion-induced phase might further improve the performances of the kq-space approach and will be the focus of future investigations.
Finally, future work should address the validation of our approach when additional microstructure parameters are considered, where the necessity for accelerating the acquisition is even stronger. The proposed approach generalizes naturally to this context, with the main difference lying in the definition of a larger and more general dictionary Φ, accounting not only for fiber orientation, but also the diameter, etc.

Institutional Review Board Statement:
The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of the Lausanne University Hospital (CHUV) and Centre for Biomedical Imaging (CIBM) as a pilot study (protocol code 20160429, date of approval 29 April 2016).

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to size reasons.

Acknowledgments:
The authors would like to thank the Center for Biomedical Imaging (CIBM) of the Geneva-Lausanne Universities for all the support received for the real data acquisition.

Conflicts of Interest:
The authors declare no conflict of interest.