EAP-Sparsity for Compressed Sensing of MS-HARDI in $$({\mathbf {k}},{\mathbf {q}})$$ -Space

-Space

Jiaqi Sun¹, Elham Sakhaee¹, Alireza Entezari¹ and Baba C. Vemuri¹

(1)

CISE Department, University of Florida, Gainesville, FL 32611, USA

Baba C. Vemuri

Abstract

Compressed Sensing (CS) for the acceleration of MR scans has been widely investigated in the past decade. Lately, considerable progress has been made in achieving similar speed ups in acquiring multi-shell high angular resolution diffusion imaging (MS-HARDI) scans. Existing approaches in this context were primarily concerned with sparse reconstruction of the diffusion MR signal $S({\mathbf {q}})$ in the ${\mathbf {q}}$ -space. More recently, methods have been developed to apply the compressed sensing framework to the 6-dimensional joint $({\mathbf {k}},{\mathbf {q}})$ -space, thereby exploiting the redundancy in this 6D space. To guarantee accurate reconstruction from partial MS-HARDI data, the key ingredients of compressed sensing that need to be brought together are: (1) the function to be reconstructed needs to have a sparse representation, and (2) the data for reconstruction ought to be acquired in the dual domain (i.e., incoherent sensing) and (3) the reconstruction process involves a (convex) optimization.

In this paper, we present a novel approach that uses partial Fourier sensing in the 6D space of $({\mathbf {k}},{\mathbf {q}})$ for the reconstruction of $P({\mathbf {x}},{\mathbf {r}})$ . The distinct feature of our approach is a sparsity model that leverages surfacelets in conjunction with total variation for the joint sparse representation of $P({\mathbf {x}}, {\mathbf {r}})$ . Thus, our method stands to benefit from the practical guarantees for accurate reconstruction from partial $({\mathbf {k}},{\mathbf {q}})$ -space data. Further, we demonstrate significant savings in acquisition time over diffusion spectral imaging (DSI) which is commonly used as the benchmark for comparisons in reported literature. To demonstrate the benefits of this approach, we present several synthetic and real data examples.

This research was funded in part by the AFOSR FA9550-12-1-0304 and NSF CCF-1018149 grants to Alireza Entezari and the NIH grant NS066340 to Baba C. Vemuri.

1 Introduction

Diffusion weighted MRI is a non-invasive way to probe the axonal fiber connectivity in the body by making the MR signal sensitive to water diffusion through tissue. In diffusion weighted MRI, the water diffusion is fully characterized by the diffusion Probability Density Function (PDF) called the ensemble average propagator (EAP) [1]. Under the narrow pulse assumption, the EAP denoted by $P(\mathbf {r})$ and the diffusion signal attenuation $E({\mathbf {q}})$ are related through the Fourier transform [1]:

$\begin{aligned} P({\mathbf {r}})=\int E({\mathbf {q}})\exp (-2\pi j{\mathbf {q}}\cdot {\mathbf {r}}) d{\mathbf {q}} \end{aligned}$

(1)

where, $E({\mathbf {q}})=S({\mathbf {q}})/S_0$ , $$S_0$$

is the diffusion signal with zero diffusion gradient, ${\mathbf {q}}$ is the vector along which the diffusion gradient is applied and ${\mathbf {r}}$ is the radial vector in the dual space defined through the Fourier relationship above. $P({\mathbf {r}})$ at each voxel, captures all the information needed to perform tractography since it is well known that the peaks of this distribution correspond to the local fiber orientations.

In order to estimate the $P({\mathbf {r}})$ , one normally acquires the diffusion-weighted MR data by sampling $E({\mathbf {q}})$ in the ${\mathbf {q}}$ -space along different diffusion sensitizing gradient directions, ${\mathbf {q}}_k$ (with $1\le k \le N$ ), spanning a unit hemisphere either over a single shell or multiple shells [2]. For every gradient direction ${\mathbf {q}}_k$ , a full 3-D acquisition in the ${\mathbf {k}}$ -space follows. In order to reconstruct $P({\mathbf {r}})$ with a reasonable angular accuracy, a substantial number of sensitizing gradient directions, on multiple shells, are necessary (e.g., $$N=180$$

). The time incurred in this extensive data acquisition is the key problem making high angular resolution diffusion imaging impractical for clinical use. Very recently however, novel techniques such as multi-band imaging have been implemented, in connection with the well known connectome project, to speed up the acquisition of MS-HARDI [3]. However, these techniques do not exploit the redundancy present in the $({\mathbf {k}},{\mathbf {q}})$ -space which is the main theme of our work in this paper. Thus, the methods presented in this paper maybe applied in addition to the multi-banding techniques to achieve further gains in acquisition time.

Compressed sensing has been applied to magnetic resonance image (MRI) acquisition quite successfully by under sampling in the ${\mathbf {k}}$ -space (frequency space) and still achieving accurate signal reconstruction from this sparse sampling [4]. In the context of diffusion MRI acquisition, there have been some attempts at applying compressed sensing concepts to diffusion spectral imaging (DSI) [5–7]. These techniques reported to use approximately 200 gradient directions to achieve accurate diffusion MR signal reconstruction and this amounts to over forty minutes of scan time which is not practical in many situations such as for movement disorder and Autism patients. As an alternative, there has been some ground breaking work reported in literature on reducing the number of directions along which the magnetic field gradients that are applied to acquire the data in order to achieve sparse reconstruction of the signal and the EAP [8–10]. They however did not apply the compressed sensing jointly to $({\mathbf {k}},{\mathbf {q}})$ -space.

More recently, Mani et al. [11] proposed compressed sensing in $({\mathbf {k}},{\mathbf {q}})$ -space by jointly under-sampling ${\mathbf {k}}$ and ${\mathbf {q}}$ spaces. This was achieved by under-sampling the ${\mathbf {k}}$ -space randomly for each direction q. Another recent development in the same vein was reported in [12], where joint $({\mathbf {k}},{\mathbf {q}})$ -space compressed sensing is proposed while the sparsity is enforced in the ${\mathbf {q}}$ -space. Naturally their reconstruction is again geared to recovering the $S({\mathbf {k}},{\mathbf {q}})$ signal, first. To achieve the EAP reconstruction their method employs the typical Fourier transform relationship between $S({\mathbf {k}},{\mathbf {q}})$ and $P({\mathbf {x}},{\mathbf {r}})$ post reconstruction of $S({\mathbf {k}},{\mathbf {q}})$ (using the dual spherical polar Fourier basis) and thus fails to exploit the incoherence between $P({\mathbf {x}},{\mathbf {r}})$ and $S({\mathbf {k}},{\mathbf {q}})$ .

In this paper, we present a novel technique based on advances in sampling theory to alleviate this time and cost expensive acquisition process that will make MS-HARDI a more viable imaging technique in the clinic. We pose the diffusion-weighted imaging problem as a six-dimensional sampling problem in the 6-dimensional $({\mathbf {k}},{\mathbf {q}})$ -space (i.e., $$(k_x, k_y, k_z)$$

and

). The diffusion sensitized MR signal and the EAP are related through the 6-dimensional Fourier transform given by,

$\begin{aligned} S({\mathbf {k}},{\mathbf {q}}) = \int _{\mathbb {R}^3}\int _{\mathbb {R}^3} P({\mathbf {x}},{\mathbf {r}})\exp ({-2\pi j({\mathbf {x}}^t{\mathbf {k}}+{\mathbf {q}}^t{\mathbf {r}})})\, d{\mathbf {r}}\, d{\mathbf {x}} \end{aligned}$

(2)

For simplicity, we omitted the scaling factor $S({\mathbf {x}},{\mathbf {0}})$ from the Fourier transform in the equation above.

In order to utilize the compressed sensing principles to achieve accurate reconstruction from partial data, the sparsity constraint is often enforced in the space domain while the sensing occurs in the (dual) frequency space. The notion of incoherent sensing formalizes the idea that sensing basis (e.g., Fourier) and representational basis (e.g., Dirac) are dual to each other; thus yielding full incoherence. Since $({\mathbf {k}},{\mathbf {q}})$ and $({\mathbf {x}},{\mathbf {r}})$ spaces are Fourier duals of each other and the acquisition occurs in $({\mathbf {k}},{\mathbf {q}})$ -space, we seek to reconstruct with sparsity constraints in $({\mathbf {x}},{\mathbf {r}})$ -space. The key distinction between our approach and existing approaches is that, enforcing sparsity in $P({\mathbf {x}}, {\mathbf {r}})$ entitles us to leverage incoherent sensing, not only in ${\mathbf {k}}$ , but also in the ${\mathbf {q}}$ -space simultaneously. Therefore, our approach presented here stands to benefit from practical guarantees for accurate reconstruction from partial $({\mathbf {k}},{\mathbf {q}})$ data. We then combine the $({\mathbf {k}},{\mathbf {q}})$ sampling with sparse reconstruction to exploit the principle of compressed sensing for reconstruction of $P({\mathbf {x}},{\mathbf {r}})$ . The key ingredient enabling sparse representation for $P({\mathbf {x}},{\mathbf {r}})$ is accomplished using surfacelet basis. The most attractive feature of surfacelet basis is the inherent directional selectivity that leads to a sparse representation in the ${\mathbf {r}}$ -space. For further details, see Sect. 2.2.

The rest of the paper is organized as follows, in Sect. 2, we present the theoretical formulation of the sampling and reconstruction problem. Section 3 contains several synthetic and real data experiments demonstrating the performance of our method. Finally, we wrap up in Sect. 4 with conclusions.

2 Formulation

In this section, we present the theoretical formulation for our full 6D Compressed Sensing (CS) and sparse reconstruction of the field of EAPs, $P({\mathbf {x}},{\mathbf {r}})$ .

2.1 Compressed Sensing

The significant achievement of the CS theory [13, 14] is the ability to reconstruct a function from partial data given the function has a sparse representation. The three ingredients of the CS framework necessary to guarantee accurate reconstruction are:

Sparsity: The function to be reconstructed needs to be sparsely representable, possibly in some transform domain.
Incoherent Sensing: The data for reconstruction must be acquired in a domain incoherent (e.g., dual) to the domain in which the function is sparsely representable.
Nonlinear Reconstruction: The reconstruction problem involves an (convex) optimization process.

In the case of diffusion MR imaging, with the presence of the Fourier dual relationship between $({\mathbf {k}},{\mathbf {q}})$ and $({\mathbf {x}},{\mathbf {r}})$ space, as illustrated in Eq. (2), the above conditions will be met when a proper sparsifying transform is applied to $P({\mathbf {x}},{\mathbf {r}})$ . In this work, we propose to use surfacelets as a choice of sparsifying basis for representation of EAPs.

2.2 The Surfacelet Transform

Measuring the diffusion of water molecules along several directions in HARDI acquisitions is an attempt to capture diffusion anisotropy. It is well known that EAPs capture this local information quite adequately. The key question then is, how best to represent the EAPs that are to be reconstructed from patially sensed data in $({\mathbf {k}},{\mathbf {q}})$ -space? Thus, the primary goal here (in accordance with the principles of CS described above) is to find a basis in which EAPs with their inherent directional information are sparsely representable.

Fig. 1.

Frequency partitioning of surfacelet transform

Wavelets, as a common choice of sparsifying transforms, lack directional sensitivity and exhibit inadequacy in efficiently capturing orientational features. As geometric generalizations to wavelets, directional decomposition methods, e.g., ridgelets [15], have been proposed to detect the orientational structures in a signal. Although these transforms can be generalized to higher dimensions, they are only optimal for 2D signals such as images. The three dimensional curvelet (3D-Curvelet) was suggested in [16] for detecting/representing directional information and geometry of the object; however, its high redundancy factor (i.e., ratio of the number of transformed coefficients to the number of signal elements) makes the problem size excessively large, hence, limits its application in diffusion MR image analysis.

Surfacelets [17], on the other hand, are real three-dimensional transforms and were shown to be particularly efficient for sparse approximation of volumetric data [18]. They have a low redundancy factor ( $\sim$ 4) and are able to capture directional information which is predominant in the ${\mathbf {q}}$ -space diffusion sensitized MR signal.

The surfacelet transform is implemented as a combination of a multi-scale pyramid with 3D directional filter banks (3D-DFB) [17]. The basis functions are a spatial domain representation of symmetric pyramids partitioning the frequency space. Figure 1 depicts the support of one surfacelet basis in the frequency domain.

Only gold members can continue reading. Log In or Register to continue