Compressed Sensing for Thoracic MRI with Partial Random Circulant Matrices

The use of circulant matrix as the sensing matrix in compressed sensing (CS) scheme has recently been proposed to overcome the limitation of random or partial Fourier matrices. Aside from reducing computational complexity, the use of circulant matrix for magnetic resonance (MR) image offers the feasibility in hardware implementations. This paper presents the simulation of compressed sensing for thoracic MR imaging with circulant matrix as the sensing matrix. The comparisons of reconstruction of three different type MR images using circulant matrix are investigated in term of number of samples, number of iteration and signal to noise ratio (SNR). The simulation results showed that circulant matrix works efficiently for encoding the MR image of respiratory organ, especially for smooth and sparse image in spatial domain.


Introduction
Magnetic resonance imaging (MRI) is one of the alternatives for intra-cardiac visualization.Compared with other imaging modalities such computed tomography (CT), MRI is relatively safe since it does not using ionizing radiation (x-rays).It is also able to obtain sequential 2D slices or 3D volumes with high spatio-temporal resolution.Despite all the advantages of MRI, some major problems are inevitable.Some issues and problems of brain MRI such as segmentation, localization, correction, and classification are discussed in [1].Long scan time in data acquisition is also found to be a major problem in MRI, especially to capture thoracic or respiratory organ.Thoracic MR image provides representation of heart, valves and major vessels anatomy and function.It is also useful to diagnose any cardiovascular problems.Ideally, physicians want to be able to observe real-time motion of respiratory organ.Unfortunately, it is impossible to get real-time motion using current MR scanner.Moreover, acquiring MR image of respiratory organ can lead to quality image degradation due to the motion during inhale and exhale of the patient.Breath holding during acquisition is one of the methods to reduce image artifact.An overview of reducing motion artifacts principles such as retrospective gating or projection reconstruction were discussed in [4].Hardware based approach to reduce imaging time using parallel data acquisition (P-MRI) method is proposed in ISSN: 1693-6930 TELKOMNIKA Vol. 10, No. 1, March 2012 : 147 -154 148 [6].Clearly, the needs to reduce imaging time which yields high quality image is the holy grail in MR imaging technology.
Compressed Sensing (CS) is a relatively new method introduced by Donoho and Candes [2][3] and expected to answer one of the issue in MR image acquisition, fast imaging time acquisition without losing its quality.It is built upon the fact that a signal (image) has a sparse representation in a priori known basis (or compressible).A signal ū with length n can be transformed using an orthogonal basis Ψ (e.g.Fourier or Wavelet basis) such that Ψ. ‫ݑ‬ ത= ‫̅ݔ‬ .
Here ‫̅ݔ‬ can be seen as a sparse representation of signal ‫ݑ‬ ത, and it has k numbers of nonzeros (ksparse), where k < n.Once known that vector ‫̅ݔ‬ is sparse or compressible, it can be compressed by projecting it onto a measurement matrix (Φ).Matrix Φ is an m x n matrix where the number of its row is much smaller than the size of the signal (m << n).The projection is simply the inner product of Φ and vector ‫ݑ‬ ത (y ത=Φ‫ݑ‬ ത).Since m << n, the equation y ത = Φ‫ݑ‬ ത leads to an ill-posed problem (underdetermined system) and it is impossible to recover ‫ݑ‬ ത (sparse representation of original signal ‫̅ݔ‬ ) from y ത.However, the sparsity of ‫ݑ‬ ത changes the impossible recovery into a perfect recovery of ‫ݑ‬ ത.Another condition needs to be satisfied to reconstruct ‫ݑ‬ ത from y ത is restricted isometry property (RIP) of matrix Φ which is obeyed by many types of matrices (e.g Random Gaussian, Bernoulli, Partial Fourier Matrices) [7].The reconstruction process can be done by greedy method (fast, but requires more measurements) or solving convex optimization problem (slow but requires less measurement data).
Lustig et al in [5] showed that MR image is suitable for CS implementation since it has sparse representation in certain transform domain.The DCT and Wavelet transform have good performance of recovering brain and angiogram MR images with 5-10% coefficients.Moreover, wavelet transform also proven to be effective for real-time data transmission for real image (JPEG2000) [13] showing that both medical and real image are actually sparse in when transformed into certain domain.But the success of reconstruction does not only merely depend on the sparse representation.Another factor is to satisfy RIP, which is a tool for analyzing the performance how efficient a matrix measures sparse signal [8].Recent widely measurement matrices that satisfy RIP with high probability is independent and identically distributed (i.i.d) random Gaussian matrix.Although Gaussian matrix provides incoherence with any sparse signal and the number of measurements for signal recovery is minimal, it requires huge memory storage and high complexity computation.Using Gaussian matrix as sensing matrix to acquire a 256x256 image with 50% sampling will require nearly gigabyte space and giga-flop operation, leads to impractical application of CS scheme for MR device.Structurally random matrices, such as Toeplitz or Circulant matrix have been shown to be fast computable and satisfy RIP property as well for almost all orthonormal matrices [9].But, the efficiency using Circulant matrix to reconstruct MR images have not been known.
In this paper, we investigated the efficiency of using circulant matrix as sensing matrix for MR Image of respiratory organ compared to partial Fourier matrix.The reconstruction is done via conventional l1-norm minimization.
The contents of this paper are organized as follows: in Section II we review the basic theory of compressed sensing and MR image of respiratory organ motion.In Section III, we discuss the methodology we use in the experiment and also the simulation results.We conclude this paper in Section IV.

Basic Theory
There are three stages in CS scheme: encoding, sensing and decoding (reconstruction).Encoding is a stage of sparsifying input signal (or image) using certain basis representation (such as Fourier or wavelet).Sensing phase is a process to measure the sparse signal representation and reducing its dimension using sensing matrix and at last, decoding is a process to reconstruct the sensing signal.
This section is divided into two parts.First part discusses about basic CS theory and its reconstruction using Total Variation minimization.Second part reviews partial random circulant matrix as sensing matrix.

CS with Total Variation Minimization
The encoding phase starts with an input of 1-D discrete signal ‫ݑ‬ ത with length n that can be represented in term of an orthonormal basis {ψ} ୧ୀଵ ୬ as: Vector ‫̅ݔ‬ is a sparse representation of vector ‫ݑ‬ ത and has k non-zeros coefficients.We say that vector ‫̅ݔ‬ is highly compressible if k << n, or in other word, non-zeros or significant coefficients much smaller than length of vector ‫ݑ‬ ത.
The sensing phase of CS needs measurement matrix for projecting vector ‫̅ݔ‬ into smaller dimension.Let is a measurements matrix with m x n dimension (where k < m << n), the linear projection of sparse vector ‫̅ݔ‬ is given by: The measurement matrix has to be independent of ‫̅ݔ‬ (non-adaptive/fixed). From (2) we have vector of measurements y ത with length m.The decoding stage is a stage to reconstruct ‫̅ݔ‬ from y ത and matrix .Since there are m equation and n unknowns to be solved, the reconstruction becomes an ill-posed problem underdetermined system) with infinitely many solutions.CS theory [1,2] states that vector ‫̅ݔ‬ can be recovered by solving convex problem.As long as ‫ݕ‬ ത is highly sparse, finding the sparsest solution will be the "best guest" to recover vector ‫̅ݔ‬ .
Unfortunately, ( 3) is a NP-hard problem and computationally impractical.An effective way to recover ‫̅ݔ‬ is to solve ℓ ଵ -norm instead of ℓ -norm.This method is also known as Basis-pursuit.
By solving ℓ ଵ -norm optimization problem, vector ‫̅ݔ‬ can be faithfully recovered.Another assumption in decoding phase is the number of measurements, ݉ ≥ ܿ݇ ‫݈݃‬ ݊, for some small c where the measurements matrix is Gaussian or ݉ ≥ ݇ ‫݈݃‬ ݊ where the measurements matrix is random partial Fourier.
In practical application, MR image acquisition collects data in the frequency domain (kspace).Lustig et al in [4] proposed a reconstruction model for MR Image.Let be the undersampled Fourier transform, the reconstruction is given by: where represents as threshold parameter for noise level.Total Variation (TV) as the sum of the absolute variations in the image can be also employed.The reconstruction model in equation ( 5) can be written as: Here, is a positive parameter.Better reconstruction in MR image from a small number of Fourier coefficients is achieved by additional TV.

Partial Random Circulant Matrix as Sensing Matrix
In equation ( 2), A is called sensing matrix.Sensing matrix plays an important role in decoding (reconstruction) stage.Recent results show that stable reconstruction from k-sparse and compressible signal can be achieved when the sensing matrix satisfies restricted isometry property (RIP) [7].Current widely used sensing matrices are random matrices and partial ISSN: 1693-6930 Fourier.Random matrices (such as Gaussian or Bernoulli random matrices) has optimal performance but requires big memory storage and complexity computation while partial Fourier enables a fast computation of Fast Fourier Transform (FFT) but lack of universality (only incoherence with sparse signal in time domain but not for smooth signal).
Circulant matrices are orthonormal matrices that are nearly incoherence with other orthonormal matrices [10], which provide universality with many types of signal.It is formed by starting with a vector with N components.The subsequent rows are acquired by shifting the previous row to the right.Given a vectorc = ሺc , c ଵ , c ଶ , … , c ିଵ ሻ ∈ ℝ , and ϕ = ϕሺcሻ ∈ ℝ ୶ , the general form of circulant matrices are given by: Consider an index set θ ⊂ ሼ1, … , nሽ as a subset of randomly locations between (1...n) with cardinality m.The rows of ϕሺcሻ can be restricted by the element of and formed ϕሺcሻ as partial circulant matrix which contains m number of rows.
A new paradigm of random filter was proposed by Tropp in [11] to captures signal s by convolution operation of random-tap FIR filter h, followed by subsampling.Figure 1 shows a block diagram of the measurement process through random filtering.

Figure 1. Block diagram of measurement process through random filtering
Using this paradigm, equation ( 2) can be written as: where D ↓ is subsample of convolution operation ሺc * u തሻ .Using the FFT to implement the convolution, equation ( 8) can be expressed as: where ℱ and ℱ ିଵ is Fourier matrix and its adjoint. is diagonal matrix of c ∈ ℝ Fourier transformation.The theoretical analysis of restricted isometries for partial random circulant matrices is given by Rauhut et al in [10].

Results and Analysis
To demonstrate the effectiveness of partial random circulant matrix in MR image of respiratory organ, we used three models of MR images of respiratory organ from different patients.The MR images were acquired using a 1.5T Achieva Nova-Dual (Philips Medical Systems, Best, NL) whole-body scanner with a 16ch SENSE TORSO XL Coil.The size of each image is 256x256.The visual result of the reconstruction using partial random circulant matrix for image in Figure 2a, 2b and 2c are shown in Figure 3, 4 and 5 respectively.As shown in Figure 3, at the sample ratio of 10% and 25%, the reconstruction images still suffer from aliasing artifact due to the sharpness original image.For reconstruction of Figure 2b and 2c (as shown in Figure 4, 5), the aliasing artifact are less at sample ratio of 25% and has better quality at sample ratio of 50% compared to Figure 3.The reconstruction works better for sparse and smooth image as in Figure 2b and 2c.
Table 1 shows the comparison of each image based on number of samples, number of iteration and reconstruction quality (SNR).The best reconstruction is gained by 50% sampling of left side MR image of respiratory organ, as the sparsest in spatial among the images.As shown in [3], exact reconstruction of MR Images using CS scheme is achievable by employing random partial Fourier matrix.To validate our results, we compared the quality of MR Images reconstruction using circulant matrix and partial Fourier matrix.The number of sampling we picked for the comparison is 35%, 50% and 75%.This is due to the unstable reconstruction using random partial Fourier matrix when the number of sampling is below 35%.Table 2 shows the comparison of image reconstruction using circulant matrix (CM) and partial Fourier matrix (PFM) in term of SNR (dB).
As shown in Table 2, MR images reconstruction using random partial circulant matrix as sensing matrix has better SNR compared with reconstruction using partial Fourier matrix.Table 3 showed that the time required to reconstruct the image using Circulant Matrix with lower sampling has better speed compared with random partial Fourier matrix (35% and 50% sampling).In the other side, when the number of sampling is higher (75%), shorter

Conclusion
We applied CS using partial random circulant matrix as sensing matrix for MR image of respiratory organ.The comparison of three different MR image of respiratory organ using partial random circulant matrix as sensing matrix has been performed.It has been found that partial random circulant matrix can efficiently reconstruct MR image of respiratory organ especially for sparse image in spatial domain.As partial random circulant matrix has better image quality and faster reconstruction time when using lower sampling compared with partial Fourier matrix, its implementation for MR scanner is very promising.Further study can be carried out by combining the CS scheme using partial random circulant matrix with particular medical image techniques such as segmentation [14].

Compressed
Sensing for Thoracic MRI with Partial Random Circulant …. (Windra Swastika) 151 Experiments were performed under Windows 7 Professional running on Dell PC Desktop with an Intel® Core ™ 2 Quad CPU, 2.66GHz and 4GB RAM.The reconstruction model we used in the experiments is fast algorithm proposed in [12].This algorithm is implemented in a MatLab program, RecPC [15].Three original images are shown in Figure 2a, 2b, and 2c respectively.Each image represents different variation of respiratory organ.Figure 2a represents front side of respiratory organ, Figure 2b represents left side and Figure 2c represents back side.

Figure 2 .
Figure 2. Original MR images of respiratory organ.(a) Front side; (b) Left side; (c) Back side

TELKOMNIKA
ISSN: 1693-6930 Compressed Sensing for Thoracic MRI with Partial Random Circulant …. (Windra Swastika) 153 reconstruction time was achieved by employing partial Fourier matrix, yet the reconstruction quality (SNR) is far below the circulant matrix.

Table 1 .
Comparison of three different image reconstructions in term of number of samples, number of iteration and image quality