Abstract
We propose a structured low rank matrix completion algorithm to recover a time series of images consisting of linear combination of exponential parameters at every pixel, from undersampled Fourier measurements. The spatial smoothness of these parameters is exploited along with the exponential structure of the time series at every pixel, to derive an annihilation relation in the k − t domain. This annihilation relation translates into a structured low rank matrix formed from the k−t samples. We demonstrate the algorithm in the parameter mapping setting and show significant improvement over state of the art methods.
Index Terms: structured low rank, Toeplitz, smoothness penalty, parameter mapping
1. INTRODUCTION
The recovery of linear combination of damped exponentials from few of their measurements is a classical problem in signal processing, starting with the seminal work of Prony [1]. Early work, which focused on uniform sampling, relied on the existence of a filter that annihilates the uniform measurements. Recently, several researchers have extended this idea to the non-uniform sampling setting, which is a more efficient measurement strategy. These methods pose the estimation as the completion of a structured matrix (e.g Hankel or Toeplitz), whose entries are the signal samples, from few of its measurements. Specifically, the annihilation filters are the null space vectors of the Toeplitz matrix, or equivalently the matrix is low-rank.
The above exponential estimation problem is of high significance in several medical imaging applications, including MR spectroscopic imaging and MR parameter mapping. These methods perform a pixel-by-pixel estimation of the exponential parameters (e.g. frequencies, decay rate, amplitudes) from an image time series. The parameters provide valuable clues about the abnormal metabolic activity and tissue micro-structural changes, which are early markers for neurological disorders. The parameters often vary smoothly in space since they depend on the underlying tissue microstructure. Unfortunately, the acquisition of the image time series at high spatial resolution results in prohibitively long scan time. A common approach to accelerate the acquisition is to acquire under sampled data, followed by reconstruction using low rank and sparsity priors.
We introduce a novel formulation that directly exploits the exponential structure of the time series and the spatial smoothness of the exponential parameters. Specifically, we couple the time series estimation problems at all the pixels into a single annihilation relation, involving a spatially smooth annihilation filter. This approach exploits the spatial smoothness naturally; we do not require additional spatial regularization priors. The annihilation relations translate to a low-rank penalty on a Toeplitz matrix, whose entries are the k − t space samples of the time series.
A challenge with the above formulation is the large size of the Toeplitz matrix. The direct implementation of the structured low-rank problem will be prohibitive for high resolution applications. We generalize our generic iteratively reweighted annihilating filter (GIRAF) framework, originally introduced for piecewise smooth images [2], to accelerate the computations. In particular, we use an iterative reweighted formulation, where we exploit the Toeplitz structure of the matrix to avoid its direct computation and storage. The GIRAF scheme approximates the linear convolutions resulting from the Toeplitz multiplications with circular convolutions, which allowed their evaluation using Fast Fourier Transform (FFT). These approximations were enabled by the fast decay of the image Fourier coefficients towards the edges of the observed region. Unfortunately, this approach breaks down in our setting since we rely on annihilation relations in k − t space, where the signal has considerable intensity at the first image. Hence, we modify the GIRAF steps using a hybrid circular-linear convolution strategy. These generalizations provide a fast algorithm that can be readily applied to large-scale exponential estimation problems. Our validations using MRI data clearly demonstrate the benefits offered by the proposed algorithm, which can simultaneously exploit the exponential structure and the spatial smoothness of the parameters.
This work has similarities to recent works on structured low-rank priors for the estimation of piecewise smooth signals and exponential signals. For example, [3] exploited the smoothness of the kx − t signal using a low-rank prior on a wavelet transform weighted Hankel matrix. This work does not exploit the exponential structure of the signal. In addition, the recovery of each ky − t slice is decoupled into a separate problem to keep the computational complexity reasonable [3]. This decoupled approach is less constrained and cannot handle efficient 2-D sampling patterns; Cartesian patterns with fully sampled kx lines were considered in [3]. In [4], the linear predictability of the time series at each pixel is exploited individually. Since this decoupled strategy is less constrained, the authors rely on additional low-rank and wavelet sparsity regularization penalties; the optimization of several regularization parameters is also challenging. While the piecewise smoothness and sparsity has also been exploited by several researchers using the structured low-rank setting [5] [6], [7] they do not exploit the exponential signal structure.
2. PROBLEM FORMULATION
We focus on the recovery of a 3-D volume ρ, consisting of a linear combination of damped exponentials at every pixel, from noisy and undersampled measurements denoted by b. We will now introduce the annihilation relations, and the resulting structured low-rank priors.
2.1. Annihilation of spatially smooth exponentials
We model the temporal signal at the spatial location as :
(1) |
Note that the temporal signal ρ[r, n] is a linear combination of L exponentials, are the amplitudes and are the exponential parameters that are dependent on the physiological process. For example, in T2 mapping applications, the exponential parameters , where ΔT is the time between two image frames and T2,i is the relaxation parameter of the ith tissue component (e.g. gray matter, CSF, white matter) at the voxel indexed by r.
It is well known that the exponential signal in (1), corresponding to a specific spatial location r, can be annihilated by a filter h[r, n] [1]:
(2) |
where (2) represents a 1-D convolution and θ is a valid index set. The coefficients of h[r, n] at r are specified by . Computing the 2-D Fourier transform (along the spatial dimension, denoted by r) of (2), we obtain
(3) |
Here, are the spatial (2-D) Fourier coefficients of the measurements ρ[r, n] and ⊗ denotes 3D convolution. Similarly, c[k, n] denotes the spatial (2-D) Fourier coefficients of h[r, n]. We assume that the parameters βi(r) describing the physiological process are bandlimited functions of the spatial variable r; this implies that c[k, n] is a 3D FIR filter. In particular, the extent of the filter c[k, n] in the spatial frequency (k) dimension controls the spatial smoothness of the filter h[r, n], while its extent in the temporal dimension (n) depends on the number of exponentials in the model (1). We express (3) in a matrix form as
(4) |
where is a linear operator that maps a 3D volume x to a lifted matrix . Here the lifted matrix has a multi-fold Toeplitz structure, such that Qc corresponds to the 3-D convolution between and the FIR filter c[k, n]. The number of columns of the matrix is equal to the assumed support of the filter. Likewise, if the measurements are restricted to a cube shaped volume in (k, n), the rows correspond to the valid convolutions.
In practice, the support of the filter c[k, n] is unknown. Let us assume that (3) is satisfied by a minimal filter c[k, n] of support . In this case, we observe that
(5) |
where d = c ⊗ e and e[k, n] is any FIR filter. Note that the support of d, denoted by Θ is larger than Δ; i.e, Δ ⊂ Θ. Hence, if we over-estimate the support of the filter c, there will be many linearly independent annihilating filters in the nullspace of , or equivalently is low rank.
To avoid oversmoothing of the parameter maps, we choose the spatial support of the 3D filter c to be large. This implies that the lifted matrix Q is rectangular (more columns than rows) in the parameter mapping setting.
2.2. Recovery using structured low-rank matrix priors
In this paper, we focus on the recovery of T2 weighted MR images from their undersampled multichannel encoded measurements, denoted by
(6) |
where is a linear operator representing Fourier under sampling and multiplication of coil sensitivities with . We pose the recovery of the signal (1) as the structured low rank matrix recovery problem:
(7) |
where λ is a regularization parameter and ∥X∥p is the Schatten p norm, defined as are the singular values of X. Here, denotes the structured multifold Toeplitz matrix, whose entries are the samples of ; corresponds to the 2-D Fourier coefficients of ρ.
3. OPTIMIZATION ALGORITHM
We rely on the iterative re-weighted least squares (IRLS) based algorithm to solve (7). The basic idea is to use the identity , where . We use an alternating minimization algorithm, which alternates between the following two sub-problems, to solve (7):
(8) |
(9) |
where ∊n → 0 is added to stabilize the inverse. We will now focus on an efficient implementation of the subproblems.
3.1. Least squares-Update
Let the rows of be denoted by . Substituting for in (9), we obtain
(10) |
Note that the term is the linear convolution between the 3-D sequences h(i) and . In the GIRAF algorithm [2], we relied on the approximation of the linear convolution by circular convolutions to accelerate its computations using FFT. We exploited the fast decay of the Fourier coefficients in GIRAF, which made the approximations valid. Here we rely on the annihilation relations in the k − t domain, where the signal does not decay rapidly in the temporal dimension. Hence we use a hybrid approach that consists of performing linear convolutions along time and circular convolutions along the spatial dimension.
Let h(i) and consist of Nt and T frames respectively and let k denote T − Nt + 1. Let each frame of the filter h(i) be of dimension N1 × N2. Now can be simplified as,
(11) |
In the above equation, represents a Toeplitz matrix formed from ; the product corresponds to the 2-D convolution between the ith frame of the filter h(i) and jth frame of . We safely approximate each of the 2-D convolutions by circular convolutions as in GIRAF and compute them efficiently using Fast Fourier transforms.
3.2. Weight-Update
We consider the Gram matrix R in (8) that has the following structure:
(12) |
with Nt column and row partitions and Ri,j is a matrix block of dimension N1N2 × N1N2. We obtain the general expression for the matrix block Rp,q as.
(13) |
We observe that the Toeplitz matrix can be expressed as , where is a matrix that restricts the convolution onto a valid index set represented by Γ, is a circulant matrix formed from and is a matrix representing zero-padding outside the filter support Λ; here P, Q are the spatial dimensions of . Now can be simplified as
(14) |
where the entries of Circ(g) are generated from the array g given by , conj denotes the conjugate operation and ◦ denotes point-wise multiplication. We evaluate the weight matrix H as
where UΛU* is the eigen decomposition of R. Hence, one choice of the matrix square root is .
4. EXPERIMENTS
A fully sampled axial 2-D dataset was acquired on a Siemens 3T Trio scanner using 12 coils and a turbo spin echo sequence. The following scan parameters were used in the acquisition: Matrix size:128×128, FOV: 22×22 cm2, TR = 2500 ms and slice thickness = 5 mm. The T2 weighted images were obtained for 12 equispaced echo (TE) times ranging from 10 to 120 ms. Post reconstruction, the T2 maps were obtained by fitting a mono-exponential model to each voxel.
We first study the effect of the filter size (dimensions of the block Toeplitz matrix) on the SNR of the images recovered from 8-fold undersampled multichannel Fourier data in Table. 1. We observe that the spatial dimensions of the filter have the largest influence on the SNR, as seen from Table. 1.(b). Specifically, we observe that filters with smaller support provides improved results, which demonstrates the benefit of exploiting spatial smoothness; a filter with size 128×128×10 fails to exploit smoothness. We also observe from the Table. 1.(a) the benefit of using temporal annihilation relations. Specifically, we obtain a 0.5 dB improvement over the filter with size 122×122×1, which only exploits joint sparsity, by increasing the length of the filter along time.
Table 1:
(a) Varying temporal dimension | (b) Varying spatial dimensions | ||
---|---|---|---|
filter size | SNR (dB) | filter size | SNR (dB) |
122x122x11 | 31.80 | 128x128x10 | 29.9 |
122x122x10 | 32.1 | 122x122x10 | 32.1 |
122x122x7 | 32 | 118x118x10 | 32.44 |
122x122x2 | 31.84 | 116x116x10 | 32.49 |
122x122x1 | 31.56 | 114x114x10 | 32.54 |
In Fig. 2, we compare the proposed approach with k − t low rank algorithm on the recovery of single coil (coil compressed) T2 weighted data from 30% uniform random measurements. We chose a filter of size 122 × 122 × 2 and a Schatten p = 0.6. We observe that the proposed scheme provides lower errors (see caption for details).
In Fig. 3, we compare the two methods on the recovery of multi-channel T2 weighted data from 12-fold (3-fold variable density + 4-fold 2×2 Cartesian) undersampled data. We chose a filter of size 114 × 114 × 10 and a Schatten p = 0.7 for the proposed scheme (see caption for details).
5. CONCLUSION
We introduced a novel structured low rank algorithm to recover a 3-D volume consisting of a linear combination of exponentials, from undersampled Fourier measurements. A convolution relation was obtained between the k − t Fourier samples and a 3-D FIR filter by exploiting the exponential structure of the time series at every pixel and the smoothness of the parameters. To speed up the computations, a hybrid approach was employed which consisted of performing circular and linear convolutions along the spatial and temporal dimensions respectively. The algorithm was applied in the context of MR parameter mapping and the reconstructed images and maps were sharper and had fewer errors than those obtained from state of the art methods.
Acknowledgments
This work is supported by NIH 1R01EB019961-01A1
6. REFERENCES
- [1].Stoica P and Moses RL, Introduction to spectral analysis. Prentice hall Upper Saddle River, 1997, vol. 1. [Google Scholar]
- [2].Ongie G and Jacob M, “GIRAF: A fast algorithm for structured low-rank matrix recovery.” [Online]. Available: http://arxiv.org/abs/1609.07429 [DOI] [PMC free article] [PubMed]
- [3].Lee D, Jin KH, Kim EY, Park S-H, and Ye JC, “Acceleration of mr parameter mapping using annihilating filter-based low rank hankel matrix,” MRM, 2016. [DOI] [PubMed] [Google Scholar]
- [4].Peng X and et al. , “Accelerated exponential parameterization of T2 relaxation with model-driven low rank and sparsity priors (MORASA),” MRM, 2016. [DOI] [PubMed] [Google Scholar]
- [5].Haldar J, “Low-rank modeling of local k -space neighborhoods (loraks) for constrained mri,” IEEE Trans. Med. Imaging, no. 3, pp. 668–681, March 2014. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [6].Jin KH, Lee D, and Ye JC, “A general framework for compressed sensing and parallel MRI using annihilating filter based low-rank hankel matrix,” 2015. [Online]. Available: http://arxiv.org/abs/1504.00532
- [7].Ongie G and Jacob M, “Off-the-Grid Recovery of Piecewise Constant Images from Few Fourier Samples,” SIAM J Imaging Sciences, in press, 2015. [DOI] [PMC free article] [PubMed] [Google Scholar]