Spatially regularized compressed sensing of diffusion MRI data

The present paper introduces a method for substantial reduction of the number of diffusion encoding gradients required for reliable reconstruction of HARDI signals. The method exploits the theory of compressed sensing (CS), which establishes conditions on which a signal of interest can be recovered from its under-sampled measurements, provided that the signal admits a sparse representation in the domain of a linear transform. In the case at hand, the latter is defined to be spherical ridgelet transformation, which excels in sparsifying HARDI signals. What makes the resulting reconstruction procedure even more accurate is a combination of the sparsity constraints in the diffusion domain with additional constraints imposed on the estimated diffusion field in the spatial domain. Accordingly, the present paper describes a novel way to combine the diffusion- and spatial-domain constraints to achieve a maximal reduction in the number of diffusion measurements, while sacrificing little in terms of reconstruction accuracy. Finally, details are provided on a particularly efficient numerical scheme which can be used to solve the aforementioned reconstruction problem by means of standard and readily available estimation tools. The paper is concluded with experimental results which support the practical value of the proposed reconstruction methodology.

signals. What makes the resulting reconstruction procedure even more accurate is a combination of the sparsity constraints in the diffusion domain with additional constraints imposed on the estimated diffusion field in the spatial domain. Accordingly, the present paper describes a novel way to combine the diffusion-and spatial-domain constraints to achieve a maximal reduction in the number of diffusion measurements, while sacrificing little in terms of reconstruction accuracy. Finally, details are provided on a particularly efficient numerical scheme which can be used to solve the aforementioned reconstruction problem by means of standard and readily available estimation tools. The paper is concluded with experimental results which support the practical value of the proposed reconstruction methodology.

Introduction
The human brain consists of about 10 11 nerve cells that can be subdivided into about 1000 different cell types, a complexity that far exceeds that of other organs of the body. A further complexity is evident in the way in which the component cells of the brain interconnect and function [1]. In contrast to other types of the cells, each neuron communicates with many target cells by means of its distinctive protoplasmatic protrusion, called an axon. Axons with similar destinations, in turn, tend to form bundles -known as neural fibre tractswhich play a pivotal role in the determination of brain connectivity. Through reconstructing the pattern of connectivity of the neural tracts in both healthy and diseased subjects, it is therefore possible to obtain an abundance of valuable diagnostic information that could be used for early diagnostics of brain-related disorders, for assessing the damage caused to the brain by stroke, tumours or injuries, as well as for planning and monitoring of neurosurgeries [2].
Central to MRI is the notion of contrast, which is typically defined by the biochemical composition of interrogated tissue as well as by the morphology of its associated parenchyma. Prevalent in MRI practice are the contrasts determined by the T 1 /T 2 relaxation times and proton density (PD). Despite their exceptional importance to clinical diagnosis, none of the above contrast mechanisms has demonstrated effectiveness in delineating the morphological structure of the white matter. It is only with the advent of diffusion MRI (dMRI) that scientists have been able to perform quantitative measurements of the diffusivity of white matter, based on which its structural delineation has become possible [2][3][4][5][6][7][8]. Although over the last two decades dMRI has developed into an established technique with a great impact on health care and neurosciences, like any other MRI technique it remains subject to artifacts and pitfalls [3]. While many of such artifacts can be overcome by using advanced hardware designs and/or more sophisticated imaging protocols [2,Ch.2], [9], one particularly critical limitation of dMRI stems from the physics of the acquisition of diffusion MR images, and therefore is impossible to resolve by operational means. Specifically, since collecting the dMRI data requires a repetitive acquisition of MR responses from the same volume for a number of diffusion-encoding gradi-ents, it is the relatively long acquisition times that greatly impair the practical value of this important imaging modality. Particularly, longer acquisition times entail a higher probability for the patient to exercise involuntary motion (typically caused by fatigue and/or stress related tremors, swallowing, uncontrollable sighing or coughing), which severely affects the quality of dMRI data. Moreover, since during the whole duration of the scan the patients are required to remain motionless, it is currently deemed ineffective to apply dMRI-based diagnosis in paediatrics as well as to patients with dementia or post-traumatic injuries, where non-compliance is typical. The problem of long acquisition times also hampers the application of dMRI for intra-operational imaging, where it could be an irreplaceable tool to use for neurosurgical planning and decision-making support [10][11][12]. Lastly, relatively long scanning times required by dMRI aggravate the problem of accessibility to MRI equipment. All the above arguments suggest that the practical value of dMRI could be improved by shortening the scanning times required for acquisition of dMRI data. A particular method to achieve such an improvement is detailed in this paper.
In this work, we adopt a general diffusion model in which each voxel within a region of interest (ROI) is allowed to support more than one fibre tract. In this case, under some general assumptions (see, e.g., [13,Sect. 3.1] for more details), the diffusion signal s(u | r) originating from a voxel with spatial coordinate r ∈ R 3 containing M (r) fibres can be modelled as [13,14] s(u | r) = s 0 (r) where u denotes the spherical coordinate, i.e.
u ∈ S 2 := v ∈ R 3 | v 2 = 1 , and α i (r) > 0 are positive weights obeying M (r) i=1 α i (r) = 1. In (1), s 0 denotes the diffusion signal obtained in the absence of diffusion encoding (i.e. the socalled "b0 image"), b is defined as a function of the shape and amplitude of diffusion-encoding gradients [15,Eq. 3.18], and {D i (r)} M (r) i=1 are 3 × 3 diffusion tensors associated with the M (r) neural fibres passing through the coordinate r. In practical settings, the spherical coordinate u is sampled at K distinct points {u k } K k=1 over the unit sphere. In this case, for each u k , MR measurements are acquired in the form of a diffusion-encoded image s k (r) := s(u k | r). As a result, a typical dMRI data set consists of a collection of such diffusion-encoded images {s k (r)} K k=1 , whose size K determines the accuracy with which the directions of local diffusion flows can be estimated.
Provided unlimited scanning time, one could measure the diffusion in thousands of orientations, making it possible to identify the directions of dominant diffusion with very high precision. For the reasons explained earlier, however, scanning times are always limited, which necessitates restriction of K to a reasonably small value. This brings us to the central question addressed in the present paper: what is a sufficient number K of diffusion-encoding directions to use? It turns out that, for some realizations of dMRI, the above question can be answered in a rigorous manner. In particular, in diffusion tensor imaging (DTI) [3-5, 7, 8, 16, 17], the reconstruction is carried out under the assumption that each voxel can support only one diffusion flow as most, which corresponds to setting M (r) = 1 for all r in (1). Accordingly, a minimum of K = 7 diffusionencoded images are theoretically sufficient to measure S 0 (r) and recover the six non-repetitive components of the symmetric tensors D(r) by means of leastsquare fitting. (In practice, however, a larger number of gradient directions is employed to reduce the estimation variance, with a typical K being between 25 and 30 [18].) Unfortunately, the accuracy of DTI is known to deteriorate dramatically at the sites where the neural fibres (or bundles thereof) cross, touch upon each other, or diverge [13,[19][20][21][22][23][24][25][26][27].
The fibre crossing problem in DTI has prompted efforts to develop dMRI methodologies which are capable of detecting multiple diffusion flows (or, equivalently, neural fibre tracts) within a given voxel. One of such techniques is High Angular Resolution Diffusion Imaging (HARDI) [14,20,[22][23][24][26][27][28], which is capable of capturing multi-modal diffusion patterns by sampling the spherical shell at a much greater number of points (usually between 60 and 100) as compared to the case of DTI. Increasing K makes it possible to describe the diffusion measurements using much more accurate models. Among these are parametric models [29][30][31][32][33][34] which allow HARDI signals to be expressed in terms of a relatively small number of prototype signal forms. Unfortunately, fitting a parametric model often entails minimization of non-convex functionals, which is a noise-sensitive and computationally intensive task, prone to the problem of local minima. The need to predetermine the optimal number of fitting terms is known to be another disadvantage of using the models of the above type.
The problems associated with parametric modelling of HARDI signals can be overcome by using non-parametric models, in which case the signals are recovered by projecting the observed data onto properly defined functional subspaces. In particular, the applicability of spherical Fourier analysis [35] to dMRI has been demonstrated in [25][26][27][28], where HARDI signals are approximated by truncated series of spherical harmonics (SH). Despite its stability and computational efficiency, however, the SH-parameterization involves a relatively large number of SHs, which suggests that the SHs cannot be an adequate basis for sparse representation of HARDI signals. The main reason for this is rooted in the fact that the energies of elementary signals (1) are concentrated alongside the great circles of S 2 , whereas the energy of SHs is spread all over S 2 , and, as a result, a relatively large number of SHs are needed to effectively "encode" each of d i . The inability of the basis of SHs to sparsely represent diffusion signals has led to the proposal of spherical ridgelets in [36,37], where it was shown that it only takes 6 to 8 spherical ridgelets on average to represent the HARDI signals with a precision exceeding the precision of their representation with 45 SHs.
The present work takes the ideas of [36,37] one step further and shows that the availability of a sparsifying basis for HARDI signals can be used to reduce the number of diffusion gradients required for data acquisition. In particular, we suggest to use the theory of compressed sensing (CS) [38][39][40][41][42][43][44] to recover the HARDI signals using the number of spherical samples K in a range of values typical for DTI (i.e. K ∈ [16,24]), thus allowing a multi-fibre analysis of dMRI data to be performed at the "acquisition cost" of a standard DTI.
It is worthwhile noting that the ideas of CS have already paved their way into the field of diffusion imaging [45][46][47][48][49]. In this regard, conceptually the closest to the proposed approach is the method reported in [46]. In spite of this similarity, however, there are two principal distinctions which make the present method a more powerful alternative. In particular, the basis functions used in [46] are limited to represent an average diffusivity and anisotropy of the white matter, thereby neglecting both intra-and inter-voxel variability of tensors D i (r) in (1). The ridgelet representation, on the other hand, is a multiresolution technique, which possesses an intrinsic ability to deal with a continuum of different diffusion scales. Second, the approach in [46] is applied in a "voxel-by-voxel" manner and it therefore does not take into consideration the spatial regularity of diffusion field. The present paper, on the other hand, proposes a novel formulation of the problem of CS-based reconstruction of diffusion signals, in which the sparsity constraints enforced in the diffusion domain are augmented by regularity constraints enforced in the spatial domain. The resulting reconstruction problem has the format of a convex minimization problem, which is solved using a specially adapted version of the split Bregman algorithm [50][51][52]. As will be shown below, the proposed algorithm results in a particularly advantageous computational structure which allows the solution to be computed via a sequence of simple and easily parallelizable steps.
The rest of the paper is organized as follows. Section II provides additional comments on the input-output structure of the proposed algorithm. The construction of spherical ridgelets is briefly outlined in Section III, whereas Section IV gives a formal description of the proposed reconstruction methodology. Some principal details on the numerical implementations of the proposed algorithm are summarized in Section V, with the results of our experimental studies reported in Section VI. Section VII finalizes the paper with a discussion and conclusions.

Problem Statement
In the centre of our considerations is the diffusion signal s(u | r) which, when normalized by its related b0-image s 0 (r), quantifies the attenuation of MR readout caused by the diffusion of water molecules in the direction u ∈ S 2 through the spatial position r ∈ R 3 . In practical settings, both u and r are discretized. Specifically, restricting u to a discrete set of orientations {u k } K k=1 prescribes the acquisition of diffusion data in the form of K diffusion-encoded images {s k (r)} K k=1 , with each s k (r) corresponding to a given u k . In this case, for a fixed r 0 , the vector [s 1 (r 0 ), s 2 (r 0 ), . . . , s K (r 0 )] T ∈ R K represents a discretization of s(u | r 0 ). Note that such a discretization follows a linear measurement model, since each sample s k (r 0 ) can be represented as an inner product of s(u | r 0 ) with a sampling function. In particular, let {ϕ k (u)} K k=1 be a Dirac basis of sampling functions defined by where δ denotes the Dirac delta function and the dot stands for the Euclidean dot product. Then, formally, with dη being the standard rotation invariant measure on S 2 .
Next, given a collection of M spherical ridgelets {ψ m (u)} M m=1 (defined below), the signal s(u | r) is assumed to be expandable as with c(r) ∈ R M being a vector of spherical ridgelet coefficients which depend on the spatial coordinate r. It is important to note that the set of spherical ridgelets is allowed to be overcomplete, implying dim Span{ψ m (u)} M m=1 < M . A practical consequence of this fact is that the definition of coefficients c(r) in (5) is, in general, not unique. This non-uniqueness is further aggravated by the fact that c(r) will have to be recovered from an under-sampled set of diffusion measurements, in which case K M . Overcoming such a severe underdetermination in the problem of estimating the ridgelet coefficients c(r) will be possible based on the fundamental premise of the theory of CS, which states that an accurate estimation of c(r) is possible if the latter is sufficiently sparse and if the sampling and representation bases are sufficiently decorrelated. While the sparsity of c(r) is rooted in the very design of spherical ridgelets [37], the incoherency between the Dirac sampling functions (3) and spherical ridgelets stems from the fact that the former have an infinitely small support, whereas the latter are "smeared" all over the unit sphere. The above properties of the basis of spherical ridgelets yield conditions for an effective application of CS, in which case one can obtain a faithful reconstruction of diffusion signals using as few as K = 16 diffusion-encoding gradients.
The proposed algorithm produces an estimate of the ridgelet representation coefficients c(r) in (5). Once available, the coefficients provide an access to the analytical definition of diffusion signals by virtue of (5). This can be used in a number of ways. One possibility could be to use the ridgelet coefficients to compute their associated orientation distribution functions [14], based on which a multi-fibre tractography analysis can be done [53,54]. Alternatively, (5) can be used to evaluate the diffusion signals over an arbitrarily fine grid of orientations. Subsequently, such refined "measurements" could be fitted using a different representation model, whose application to the original data would not have been possible without causing severe underestimation errors. Deconvolving the refined data to estimate the underlying fibre orientation functions [55][56][57][58] would be another important option to follow. In this paper, we refrain from questioning which of the above possibilities is more advantageous over the others. Our sole objective here is to specify a signal processing algorithm which can be used to recover HARDI signals, while using the number of diffusion-encoded images typical for a standard DTI.
Finally, it should be noted that the primarily purpose of the proposed methodology is to improve the value of HARDI in terms of its time efficiency. Since the improvement is achieved through merely decreasing the number of diffusion-encoding gradients, the proposed method by no means abrogates the use of fast imaging protocols [59][60][61] to further accelerate the data acquisition. Furthermore, an additional speed-up can be achieved via applying CS to reconstruct the diffusion encoded images s k (r) from their subcritical samples in the spectral domain [62][63][64][65][66]. Generally speaking, we believe this is a combination of such software and hardware technologies which will eventually lead to substantantial improvements in the practical value of HARDI-based diagnosing. In this paper, however, we confine our contribution to showing one particular way of attaining this important objective.

Spherical ridgelets
It is the property of spherical ridgelets to provide sparse representation of diffusion signals described by (1) which makes them an unparalleled tool for CSreconstruction of HARDI data. To avoid repetitions, in what follows, we present only the most principal points of ridgelets design, while their detailed description can be found in [37].
Spherical ridgelets are constructed using the fundamental principles of wavelet theory [67,68]. Specifically, let x ∈ R + and ρ ∈ (0, 1) be a positive scaling parameter. Further, let κ(x) = exp{−ρ x (x + 1)} be a Gaussian function, which we subject to a series of dyadic scalings to result in where j ∈ N := {0, 1, 2, . . .}. Subsequently, the Gaussian-Weierstrass scaling function χ j,v : S 2 → R at resolution j ∈ N and orientation v ∈ S 2 can be defined as given by [69,70] where P n denotes the Legendre polynomial of order n. It is worth noting that the L 2 energy of χ j,v is concentrated around the spherical point v, with this concentration becoming more and more localized when j approaches infinity. The spherical ridgelets are designed with the help of the Funk-Radon transform which, for an arbitrary continuous function f : S 2 → R, is defined as with σ(v) denoting the great circle perpendicular to direction v, i.e. σ(v) := u ∈ S 2 | u · v = 0 . Subsequently, following [37], the semi-discrete frame U of spherical ridgelets can be defined as where the spherical ridgelet functions ψ j,v are obtained from χ j,v according to with χ −1,v (u) ≡ 0. Using (7), the ridgelets (10) can be redefined in a closed form as (see [35] for details) where κ −1 (n) = 0, ∀n and The set U in (9) is infinite-dimensional, and hence is not suitable for practical computations. To define a discrete counterpart of U, one has first to restrict the values of the resolution index j to a finite set {−1, 0, 1, . . . , J}, where J defines the highest level of "detectable" signal details. Additionally, the set of all possible orientations v ∈ S 2 of spherical ridgelets needs to be discretized as well. To find a proper discretization scheme, we first note that the construction in (11) suggests that the bandwidth of the spherical ridgelets (and therefore the dimensionality of the functional space they belong to) increases proportionally to 2 j . Since the space of spherical harmonics of degree n has a dimension of (n + 1) 2 , it seems to be reasonable to define the number of ridgelet orientations at resolution j to be equal to M j = (2 j+1 m 0 + 1) 2 , with m 0 being the smallest spherical order resulting in κ 0 (m 0 ) ≤ for some predefined 0 < 1 (e.g. = 10 −6 ). Consequently, for each j, a total of M j orientations {v i j } Mj i=1 are chosen so that a discrete counterpart of U can now be defined as where the subscript d stands for "discrete". It should be noted that, although the set U d is composed of continuously defined functions, its dimension is finite, since U d consists of a total of M = J j=−1 (2 j+1 m 0 + 1) 2 spherical ridgelets. To slightly simplify our notation, in what follows, the spherical ridgelets in U d will be indexed as ψ m (u), with m = 1, 2, . . . , M being a combined index accounting for both different resolutions and orientations.
Given a sampling set of K diffusion-encoding orientations {u k } K k=1 , one can use (11) to compute the values of the spherical ridgelets in U d over the sampling set 1 . The resulting values can be stored into a K × M matrix A defined as Then, for a given vector s(r) : T ∈ R K of the measured values of a diffusion signal s(u | r) at the spatial location r, the model (5) asserts the existence of representation coefficients c(r) ∈ R M such that where e(r) accounts for both model and measurement noises. Clearly, the nonnegligibility of e(r) along with the fact that K M makes the problem of recovering the representation coefficients c(r) from s(r) a very challenging inverse problem, our solution to which is presented next.

Proposed reconstruction framework
Let Ω represent the volume within which diffusion measurements are acquired. Also known as a region of interest, Ω is assumed to be a bounded rectangular subdomain of R 3 , i.e. Ω : Let Ω d be a discrete subset of Ω, which represents the spatial locations at which the diffusion signal is measured. Specifically, Ω d is assumed to be a uniform lattice which can be formally defined as where 0 ≤ i 1 < N x , 0 ≤ i 2 < N y , and 0 ≤ i 3 < N z are sampling indices in the direction of x, y and z coordinates, respectively. Let K be the number of diffusion-encoding gradients used for HARDI data acquisition, and let the corresponding gradient orientations be denoted by {u k } K k=1 , where u k ∈ S 2 . For each of these values of u k , MRI measurements result in its associated diffusion-encoded image s k , which can be formally viewed as a mapping from Ω d to R. For the sake of convenience, each s k can be stored and manipulated as an N x × N y × N z array of real numbers, namely s k ∈ R Nx×Ny×Nz . Alternatively, at a given coordinate r ∈ Ω d , one can combine the values s 1 (r), s 2 (r), . . . , s K (r) into a column vector s(r) := [s 1 (r), s 2 (r), . . . , s K (r)] T ∈ R K (as it was already done in (15)). This vector can then be regarded as a vector of discrete measurements of an associated HARDI signal s(u | r) measured at orientations {u k } K k=1 . It is worth noting that, according to the above notations, the value s k (r) admits a twofold interpretation, viz. either as the k th coordinate of vector s(r) or as the value of image s k at spatial position r.
When combined together, the continuum of vectors s(r) can be regarded as a discrete vector field s : Ω → R K , in which case s(r) has a natural interpretation of the value of s corresponding to position r. The vector space V of such vector fields can be endowed with the standard inner product being the standard inner product between the scalar-valued images s 1 k and s 2 k . Accordingly, congruent to the definition in (17), where · 2 and · F denote the Euclidean vector and Frobenius matrix norms, respectively.
Another norm on V that we shall make use of is the total variation (TV) semi-norm which is defined as follows. First, let us define the total variation of the k th component s k of the field s in a standard manner as where Consequently, the TV norm of the discrete vector field s can be defined in terms of the TV-norms of its K components as Thus, for example, α = 2 was used in the TV-denoising method reported in [71].
In this paper, we use α = 1. Now, let {ψ m } M m=1 be a set of spherical ridgelets defined by (13), which is assumed to be rich enough so that each HARDI signal can be expressed according to (5). Analogously to the discrete measurements s(r), the representation coefficients c(r) corresponding to different voxels r can be aggregated into a vector field c ∈ U, where U : Ω d → R M (with c(r) being the value of c observed at r). Although it is possible to endow the vector space U with both the 2and TV-norms by analogy with (18) and (20), it will be particularly useful to consider the 1 -norm of c which can be defined as Using the definitions of the vector fields V and U as well as the definition of A in (14), a connection between V and U is established by means of a linear map A : U → V that is defined as given by Consequently, using A, one can define the HARDI data formation model as where e ∈ V is supposed to account for both measurement noise and modelling errors, and it is assumed to have a relatively small 2 -norm e V,2 ≤ . The model (23) suggests a reduction of the problem of estimation of HARDI signals to the problem of estimation of their corresponding representation coefficients c from the discrete and noisy measurements s. Furthermore, as our main intension has been to recover the coefficients c using as few diffusion-encoding gradients as possible (implying K M ), there is an infinite number of solutions which would fit the constraint A{c} − s V,2 ≤ . However, if it is known a priori that, for each r ∈ Ω d , the vector of representation coefficients c(r) is sparse, then a useful estimate of c can be obtained as a solution to the following convex optimization problem [38][39][40][41][42][43][44] It should be noted that the optimization problem (24) is equivalent to solving and therefore, under the assumption of spatially homogeneous noise e, the problem (24) is separable in the spatial coordinate r. This means that an optimal field c can be recovered by solving for its components min c(r) independently at each r ∈ Ω d . While computationally attractive, the above solution is suboptimal, since it completely disregards the dependencies which are likely to exist between spatially adjacent HARDI signals. A possible way to take such dependencies into consideration is to require the noise-free version of the measured signal s to possess a minimal TV norm among all possible candidate solutions [72]. This requirement can be translated into the following minimization problem where the role of γ > 0 is to balance the relative influence of the sparse and TV terms in the above cost function. The optimization problem (30) can be rewritten in its equivalent Lagrangian form for some optimal values of λ > 0 and µ > 0 [73]. Below, we are going to specify a particular, computationally efficient method for solving (31). In this connection, it is instructive to outline the following two instances of (31).

Sparse-only reconstruction
When µ = 0, the functional in (31) becomes separable in the spatial variable r in the sense that, in such a case, an optimal c can be recovered by solving min c(r) for each c(r) independently. Note that (32) can be considered to be a Lagrangian form of the optimization problem (28). There exist a broad spectrum of methods which could be used for solution of (32). Some particularly attractive algorithms seem to be those exploiting the principle of iterative shrinkage (aka iterated thresholding) [74][75][76][77]. While the non-differentiability of 1 -norm in (32) rules out the applicability of gradient-based optimization tools, iterative shrinkage methods are capable of finding a solution of (32) through iterative application of a first-order, fixed-point update rule. Specifically, many algorithms of this kind find an optimal solution as a stationary point of the sequence of estimates produced by where S τ {t} = sign(t)(|t| − τ ) + denotes the operator of soft thresholding and ν is chosen to obey ν > AA T . In the present paper, a modification of the iterative update in (33), known as the fast iterative shrinkage thresholding algorithm (FISTA) [78], was employed due to its considerably faster convergence as compared with many alternative "accelerated" methods. It should be emphasized that, while being suboptimal from the viewpoint of spatial-domain regularity, the solution of (32) through iterative shrinkage is advantageous in two important practical ways. First, it suggests considerable storage reduction, since the thresholding operator in (33) sets to zero the representation coefficients with amplitudes less or equal to λ/ν in absolute value. It makes it possible to use sparse data formats to store and manipulate the representation coefficients. Second, the fact that the estimation of c(r) is performed at each voxel independently suggests a natural way to speed up the overall estimation process though parallel computing on a multicore system.

TV-only reconstruction
When λ = 0, solving the optimization problem (31) is equivalent to simultaneously solving K optimization problems of the form where k = 0, 1, . . . , K − 1 and [A{c}] k denotes the k-th component of the vector field A{c} ∈ V. Let [A{c}] k be denoted by u k , i.e. u k := [A{c}] k . Then, reformulated with respect to u k , the problem (34) can be rewritten as in which case it can be recognized as the problem of TV-denoising of the diffusion-encoded image s k [72]. It is important to note that, in contrast to (34), the problem (35) can be solved for each k independently, in which case we say that the estimation becomes separable in the diffusion direction.
The current arsenal of methods which can be used for solving (35) is broad. Originated from the work of Rudin et al [72], TV-denoising methods currently include gradient-based optimization methods [79], first-order methods [80], iterative shrinkage methods [81,82], and Bregman-type iterative algorithms [50,51,83]. The first-order methods appear to be a particularly attractive option when relatively large data arrays need to be processed, which is the case relevant to dMRI. Consequently, in the present paper, we employ the algorithm of [84] for the simplicity and elegancy of its implementation as well as for its outstanding convergence properties.

Solution using Split Bregman Algorithm
Directly solving the original problem (31) is difficult because of the compound nature of the regularization it involves. The split Bregman approach [50,51] allows reducing (31) to a simpler form through introduction on an auxiliary variable u ∈ V, which can be viewed as a noise-free version of the data field s. Particularly, using u (31) can be redefined as min c,u s.t. A{c} − u 2 V,2 = 0.
Then, starting from an arbitrary b 0 ∈ V, the Bregman algorithm [85] finds optimal c and u through the following iterations u t+1 , c t+1 = arg min c,u for some γ > 0 2 . The functional in (37) is supposed to be minimized over two variables, i.e. u and c. However, due to the way the 1 and TV components of this functional have been split, the minimization can now be performed efficiently by iteratively minimizing with respect to u and c separately. The resulting iteration steps are Step 1: c t+1 = arg min Step 2: u t+1 = arg min Note that the functional at Step 2 contains two quadratic terms which can be combined together to result in Step 2: u t+1 = arg min To cause a substantial reduction in the value of the cost functional in (37), Step 1 and Step 2 should be applied recursively for a predefined number of times before the Bregman parameter b t is updated according to (37). It was argued in [51], however, that the extra precision gained through such a repetitive application of Step 1 and Step 2 is likely to be "wasted" when b t is updated. Consequently, it was suggested in [51] to perform these steps only once per iteration cycle. It is interesting to note that, in this case, the split Bregman algorithm transforms into the alternating directions method of multipliers (ADMM) [52], whose convergence is guaranteed by the Eckstein-Bertsekas theorem [86] (see also Theorem 3.1 in [52]).
The final algorithm is summarized below. Lines 3-4 of Algorithm 1 correspond to Step 1 in (38), while lines 5-6 correspond to Step 2. An even more important fact to notice is that the optimization problem in line 4 is separable in the spatial coordinate r. This optimization, therefore, can be performed at each voxel independently as discussed in Section 4.1. Moreover, the optimization problem in line 6 is separable in the diffusion coordinate k, and hence it amounts to applying TV-denoising to each of the K components of u independently as discussed in Section 4.2.

Technical details of the experimental study
Both the choice of diffusion-encoding directions {u k } K k=1 and of the orientations of spherical ridgelets require the use of a sampling scheme. In dMRI, one of the standard methods used to distribute a given number of spherical points over S 2 in a quasi-uniform manner is by means of the electrostatic repulsion algorithm (also known as the Thomson problem). However, since the diffusion signals are symmetric (implying s(u | r) = s(−u | r)), it is a unit hemisphere, not the entire S 2 , which actually needs to be discretized. In view of the absence of a formulation of the Thomson problem for hemispherical domains, a common practice is to run the standard procedure for twice as many points as needed, followed by keeping only a half of the resulting configuration. However, as the retained points are not explicitly constrained to lie on a hemisphere, they may include nearly antipodal pairs which are likely to introduce undesirable dependencies between the diffusion measurements as well as between the basis functions. This limitation can be overcome through adapting a different sampling strategy. Particularly, in this paper, both the diffusion-encoding directions and the orientations of spherical ridgelets have been defined by using the method of generalized spiral points [87,88], in which the sampling points are arranged along a spiral in such a way that the distance between the points along the spiral is approximately equal to the distance between its coils. This method is easily adaptable for sampling of the "northern" hemisphere (i.e. u ∈ S 2 | u · [0, 0, 1] T ≥ 0 ), providing a nearly uniform, unique and analytically computable coverage which is in no respect inferior to the one produced by solving the Thomson problem.
To assess the performance of the proposed algorithm under controllable conditions, experiments with simulated data sets have been performed. In this case, the HARDI signals were generated according to model (1) with different values of M (r), D i (r), and s 0 (r) = 1, ∀r. The resulting signals were contaminated by variable levels of Rician noise, giving rise to a set of different SNRs. In this work, we adapt the standard definition of the SNR as SNR = 20 log 10 where s ands denote an original signal and its noise-contaminated version, respectively, and the norms are computed as defined by (18). It is worthwhile noting that the optimal values of regularization parameters λ and µ in (31) are normally a function of the noise level. In the present paper, however, no attempts have been extended to optimize these values for different SNRs. Instead, it was found that λ = 0.03 and µ = 0.05 provided acceptable estimation results in all the simulation scenarios, and hence these values have been used throughout the whole study. Following [37], the scaling parameter ρ in (6) was set to 0.5 and the resolution parameter J in (13) was set to be equal to 1, corresponding to a total of three resolution levels. The number of spherical ridgelet orientations were predefined with m 0 = 4, resulting in M −1 = 16, M 0 = 49 and M 1 = 169 ridgelets spanning the resolution levels j = −1, j = 0 and j = 1, respectively. Thus, the total number of spherical ridgelets used in the reconstruction was equal to 234.
To quantitatively compare the reconstruction results produced by the proposed and references methods for different numbers of sampling directions K and various SNRs, three performance measures were used. The first of the three was the normalized mean-squared error (NMSE) defined as with s(r) being a reference HARDI signal corresponding to location r andŝ(r) being its estimate. Depending on the nature of a specific experiment, the reference signal can be either a simulated signal discretized at 642 spherical points obtained by the 3rd order tessellation of the icosahedron or a signal reconstructed using a maximum possible number of diffusion-encoding orientations. One of the most valuable outcomes of HARDI is in providing an access to computation of orientation distribution functions (ODFs) -the functions whose modes are likely to coincide with the direction of local diffusion flows [14]. Both the SH-based [28] and ridgelet-based [37] methods of reconstruction of HARDI signals come with analytical expressions which relate the HARDI signals to their corresponding ODFs. The latter can in turn be used to recover the directions of local diffusion flows (or, equivalently, the orientations of their related fibre tracts) using, e.g., the steepest ascent procedure detailed in [37]. Suppose u 0 is the true direction of a diffusion flow andũ is its estimate. Then, the angular orientation error δ can be defined (in degrees) as δ = 180 π arccos(u 0 ·ũ).
In this paper, as a performance measure, we use an average angular orientation error which is obtained by averaging the values of δ computed for all "fibres" within a specified Ω d . The last performance measure used in this work is the probability P d of false fibre detection. To define P d , let M (r) be the true number of fibre tracts passing through voxel r (as defined by model (1)). Also, letM (r) be an estimated number of fibres, which is equal to the number of modes (maxima) of the ODF recovered at position r. Then, one can define In addition to the quantitative comparison, the reconstruction results will be evaluated through visual comparison as well. In this paper, we choose to visualize spherical functions by means of 3-D surface plots. Such a plot tends to project away from the origin of R 3 in the directions along which a spherical function is maximized, while passing near the origin in the directions where the function approaches zero.
Finally, our choice of reference methods was motivated by the scope of the main statements made in this paper. First, since we argue that the frame of spherical ridgelets is optimally suited for CS-based reconstruction of HARDI signals, its performance has to be compared with that of alternative representation systems. In particular, the basis of spherical harmonics up to the order 8 inclusive has been used for a different definition of the sensing matrix A in (14). (Note that, in the case of a real and symmetric analysis, this SH-basis consists of 45 functions.) Additionally, the representation system proposed in [46,47] has been exploited in the comparative study as well. This system is formed by applying a set of rotations to a Gaussian kernel of the form d(u) = exp{−b (u T D 0 u)}, with b defined as in (1) and D 0 equal to a (diagonal) diffusion tensor having a mean diffusivity of 766 mm 2 /s and a fractional anisotropy of 0.8. Following [47], the number of rotations (and hence the number of Gaussian basis functions) was set to be equal to 253. For the convenience of referencing, the CS-based reconstruction methods using the spherical ridgelets, the 8th-order spherical harmonics, and the rotated Gaussian kernels will be referred below to as the RDG, SH8 and GSS algorithms, respectively.
To assess the significance of the proposed spatial regularization, all the above algorithms have been applied with two different values of µ in (31), viz. µ = 0 and µ = 0.05. Note that, in the first of these cases, the spatial regularity is ignored, which leads to the sparse-only reconstruction discussed in Section 4.1.
In the second case, on the other hand, the spatial regularity is taken into account and the reconstruction is performed by means of the split Bregman algorithm of Section 5.

In silico experiments
To assess the performance of the proposed and reference methods under controllable conditions, two simulated data sets were used. The first set (referred to below as Phantom #1) had a spatial dimension of 12 × 12 pixels, and consisted of two "fibres" crossing each other at the right angle as it is shown in the upper row of subplots of Fig. 1. In addition, each pixel in the set was assigned an extra diffusion flow in the direction perpendicular to the image plane. As a result, the number of diffusion components M (r) in Phantom #1 varied between 1 and 3. Subsequently, model (1) was used to generate corresponding diffusion-encode images {s k } K k=1 for a range of different values of K. Two different values of b, namely b = 1000 s/mm 2 and b = 3000 s/mm 2 were used for data generation. The diffusion tensors D i (r) in (1) were obtained by applying rotations to a tensor of the form D 0 = diag ([α, β, β]), where α and β were equal to 1700 · 10 −6 and 300 · 10 −6 , respectively. Note that the mean diffusivity and fractional anisotropy of D 0 are equal to 766 mm 2 /s and 0.8, respectively. Thus, the same diffusion tensors were used for data synthesis and for the construction of basis functions in the GSS algorithm, thereby allowing the latter to perform under the best possible conditions.  The lower row of subplots in Fig. 1 depict four examples of the diffusionencoding images obtained for Phantom #1 before their contamination by Rician noise. One can see that the images are piecewise constant functions, which appears to be in a good agreement with the bounded-variation model suggested by (31). However, real images may be more complicated than that. Accordingly, to test the robustness of the proposed regularization scheme, a different in silico phantom was designed. Phantom #2 had a spatial dimension of 16 × 16 pixels and it was obtained through supplementing the configuration of Phantom #1 by an additional circular "fibre" as shown in the upper row of subplots in Fig. 2. The lower row of subplots of the figure show a subset of the resulting diffusionencoded images, which can be seen to no longer exhibit a piecewise constant behaviour characteristic for Phantom #1.
The simulated diffusion-encoded images were contaminated by three different levels of Rician noise, giving rise to SNR of 24, 18 and 12 dB. Some typical examples of the resulting images are demonstrated in Fig. 3, where the upper row of subplots depict a noise-free version of one of the diffusion-encoded images pertaining to Phantom #1 along with its noise-contaminated counterparts. The lower row of subplots in Fig. 3 depict an analogous set of examples for Phantom #2. Observing the figure, one can see that the SNR values have been chosen so as to cover a range of possible noise scenarios, which could be characterized as moderate-to-severe contamination.
As it was mentioned earlier, in our in silico study we compared the performance of three different representation bases, i.e. spherical harmonics (SH8), Gaussian kernels (GSS) and spherical ridgelets (RDG). All the resulting algorithms have been further subdivided into two different types, depending on whether or not the spatial regularization was engaged. Thus, in the absence of the spatial regularization (corresponding to µ = 0), the reconstruction has been performed on a voxel-by-voxel basis, as detailed in Section 4.1. For the convenience of referencing, the corresponding algorithms will be referred to below as SH8-CS, GSS-CS, and RGD-CS. In the case of µ > 0, the estimation has been carried out using the split Bregman method of Section 5. The corresponding algorithms will be referred below as SH8-TV, GSS-TV, and RGD-TV.
The upper subplot of Fig. 4 shows the original field of ODFs of Phantom #1 (corresponding to b = 3000 s/mm 2 ), which have been computed based on Tuch's approximation [14] (i.e. by applying the Funk-Radon transform to the diffusion signals). At the same time, the middle row of subplots of Fig. 4 show the ODFs recovered by (from left to right) SH8-CS, GSS-CS and RDG-CS with K = 16 and SNR = 18 dB. One can see that the inability of the SH basis to sparsely represent HARDI signals results in a poor performance of SH8-CS. A better result is obtained with GSS-CS, which uses a basis of rotated Gaussian kernels, and therefore has a potential to represent the HARDI signals in a sparse manner. Unfortunately, the excessive correlation between the Gaussian basis functions adversely affects the ability of this method to withstand the effect of noise. Consequently, the reconstruction obtained using GSS-CV suffers from sizeable errors. The RDG-CS method, on the other hand, provides an estimation result of a much higher quality, albeit some inaccuracies are still noticeable in the   central part of the phantom. The reconstruction accuracy improves dramatically when the spatial regularization is "switched on", as it is demonstrated by the bottom row of subplots in Fig. 4. Specifically, while SH8-TV is still unable to provide a valuable reconstruction, the estimates obtained using GSS-TV and RDG-TV represent correctly the "flow structure" of Phantom #1. Moreover, among the latter two methods, RDG-TV is clearly the best performer, resulting in a close-to-ideal recovery of the original ODFs. The superiority of RDG-TV over the alternative methods is further evident in the results presented by Fig. 5, which depicts the reconstructions obtained for Phantom #2 (with the same values of b, K and SNR as above).
In general, the reconstruction results obtained using SH8-CS and SH8-TV have been observed to be of a lower quality in comparison to the other methods under consideration. For this reason, in what follows, only the GSS and RDG methods are compared. Thus, Fig. 6 contrasts the performances of GSS-CS, GSS-TV, RDG-CS and RDG-TV in terms of the NMSE criterion. One can see that the best performance here is attained by the RDG-TV algorithm, which results in the smallest values of NMSE for both phantoms and for all the tested values of b, SNR and K. It is also interesting to note that the incorporation of spatial regularization allows GSS-TV to outperform RDG-CS, with the effect of the regularization becoming more pronounced at lower SNRs. On the whole, all the NMSE curves demonstrate an expected behaviour, with the error values increasing proportionally with a decrease in SNR, while going down with an increase in the number of diffusion-encoding gradients K. However, as opposed to the others, the NMSE curves obtained with RDG-TV are characterized by a relatively low rate of convergence, which indicates a reduced sensitivity of RDG-TV to the value of K.
The above algorithms have been also compared in terms of the angular error (42). The results of this comparison are summarized in Fig. 7, which again indicates that the most accurate reconstruction is obtained using the RDG-TV method. In general, the angular error tends to grow as SNR decreases and to converge to a minimum as K increases. As opposed to the case of NMSE, however, there is an additional dependency of the angular error on the type of a phantom in use as well as on the b-value. In particular, the errors obtained for Phantom #2 are (on average) greater than those obtained for Phantom #1. This discrepancy is rooted in the fact that Phantom #2 has a more complex "fibre structure" as compared to Phantom #1. In particular, while the "fibers" of Phantom #1 are designed to cross each other at the right angle, the "fibres" of Phantom #2 are allowed to decussate at much smaller angles, which makes them much harder to resolve. Moreover, this effect becomes more noticeable with a decrease in the b-value, which reduces the resolution of q-ball imaging. Finally, we notice that, on average, GSS-TV performs better than RDG-CS (though still worse than RDG-TV), which justifies the value of spatial regularization.
The comparison in terms of the rate of false fibre detection P d (43) was last in the line of our in silico performance tests; its results are shown in Fig. 8. One can see that, in the case of Phantom #1, RDG-TV yields a virtually zero false detection rate for both values of b, whereas the other methods result in considerably higher values of P d (mainly due to the detection of spurious local maxima in the estimated ODFs). The situation is different for Phantom #2, where all the compared methods yield sizeable errors (especially for b = 1000 s/mm 2 ). However, in comparative terms, the most accurate reconstruction is still obtained by means of the proposed RDG-TV algorithm.

In vivo results
As the next validation step, experiments with real HARDI data were carried out. The proposed algorithm was tested on human brain scans acquired on a 3-Tesla GE system using an echo planar imaging (EPI) diffusion-weighted image sequence. A double echo option was used to suppress eddy-current related distortions. To improve the spatial resolution of EPI, an eight channel coil was used to perform parallel imaging by means of the ASSET technique with a speed-up factor of 2. The data were acquired using 51 gradient directions (quasi-uniformly distributed over the northern hemosphere) with b = 1000 s/mm 2 . In addition, eight baseline (b0) scans were acquired, averaged and used for normalization. The following scanning parameters were used: TR = 17000 ms, TE = 78 ms, FOV = 24 cm, 144 × 144 encoding steps, and 1.7 mm slice thickness. All scans had 85 axial slices parallel to the AC-PC line covering the whole brain.
The main question addressed through the in vivo experiments has been whether or not it is possible to supersede the spatial regularization by prefiltering of HARDI signals. To this end, the RDG-CS algorithm was applied first to the HARDI data containing the full set of K = 51 diffusion gradients. (Note that such dense reconstruction is analogous to the one reported in [37], where the latter is shown to outperform the SH-based estimation [28].) The resulting ODFs have been used as a fiducial against which different reconstruction results were compared.
As the next step, three different subsets of 16, 24 and 32 spherical points were composed out of the original set of 51 diffusion gradients. Within each of these subsets, their corresponding points were chosen so as to result in a quasiuniform coverage of the northern hemisphere. Accordingly, the HARDI data were rearranged into three data sets of size 144 × 144 × 85 × 16, 144 × 144 × 85 × 24 and 144 × 144 × 85 × 32 to emulate compressed sensing data acquisition. The above sets were used to assess the performance of different reconstruction methods. Unfortunately, we have not succeeded to find conditions under which the SH8 and GSS algorithms would provide stable reconstruction results (either with or without pre-filtering). For this reason, only the RDG-CS and RDG-TV algorithms are compared below.
The upper row of subplots in Fig. 9 show the generalized anisotropy (GA) [14] image of a coronal cross-section of the brain along with the reference field of ODFs corresponding to the region indicated by the yellow rectangular. Anatomically, this region is expected to contain the fibre bundles of corona radiata as well as those of superior longitudinal and arcuate fasciculi. The middle row of subplots in the same figure depict the ODFs reconstructed by RDG-CS using    One can see that the quality of reconstruction progressively improves as K increases. It is important to note that, before applying the RDG-CS algorithm, the diffusion-encoded images had been pre-processed by a TV filter to reduce the effect of measurement noises on the estimation result. However, this pre-processing appears to be not nearly as effective as the spatial regularization of the RDG-TV algorithm, whose reconstruction results are shown in the bottom row of subplots in Fig. 9. The above conclusion is further supported by an additional example of Fig. 10, which shows the reconstructions pertaining to the indicated area within an axial cross-section of the brain. (The relevant fibre bundles here are those of cingulum and corpus callosum). As in the previous example, one can see that the most accurate reconstruction is attained by means of the proposed RDG-TV method. The superiority of RDG-TV is also confirmed by the quantitative figures of Table 1, which summarizes the NSME obtained by the compared algorithms for different values of K.

Discussion and Conclusions
When considered as a whole, the HARDI signals which pertain to a given volume of interest can be described as multi-valued (or, more generally, measure-valued) functions from a subset of R 3 to the space of square-integrable spherical functions L 2 (S 2 ). Such functions can be thought of as if they had two "modes of variation" -one in the spatial and another in the diffusion domain. Although applying various inverse problems (a particular instance of which is addressed by the theory of CS) along the spatial and diffusion coordinates independently is not new to the community of medical imaging scientists, formulating a CS reconstruction problem in both domains simultaneously has not been proposed before. Accordingly, the present paper introduced the RDG-TV algorithm which exploits the above idea and can be used for reliable reconstruction of HARDI signals from as few as K = 16 diffusion-encoded scans (as compared to 60-80 scans required by existing reconstruction tools). The algorithm exploits the fact that HARDI signals can be sparsely represented by spherical ridgelets in the diffusion domain, while their associated diffusion-encoded images have bounded variation in the spatial domain. Moreover, it has been shown experimentally that either using different representation bases or excluding the spatial regularization would result in considerably less accurate reconstruction results. At the practical level, the reconstruction is implemented based on the split Bregman approach (with 20 being the maximum number of Bregman iterations used in this study). The resulting algorithm alternates between two estimation stages (38): first, a sequence of basis pursuit de-noising problems are solved independently on a voxel-by-voxel basis, followed by applying a TV filter to a total of K discrete images. Such computations are straightforward to accelerate using a multicore processing, which is another advantage of the proposed reconstruction method. We believe that the algorithm presented in this paper can be improved in a number of ways. First, the square metric used to assess the model fidelity could be replaced by a different metric, which would be more specific to the nature of Rician noise. Second, the fact that diffusion signals are positive-valued could be explicitly incorporated into the reconstruction process in the form of additional constraints. Lastly, the bounded variation model could be substituted by an alternative model, which could (possibly) provide a better account for the spatial regularity of HARDI signals. Exploring the above options constitutes essential part of our ongoing research.
Finally, as the experimental study reported in this paper was comparative in its nature, it was not really important what method to use for approximation of ODFs. Specifically, the present results have been obtained using Tuch's approximation [14]. However, more accurate computation of ODFs is possible based on the solid angle formulation as detailed in [89,90] 3 . It should also be noted that the technique proposed in [89] can be applied to multi-shell HARDI data (i.e. the data acquired for a range of b-values). Until recently, collecting such data has been deemed impractical due to extremely long acquisition time required. We believe, however, that the proposed method for CS-based reconstruction of HARDI data has a potential to help multi-shell HARDI develop into a clinically relevant tool of diagnostic imaging.