Fast magnetic resonance imaging simulation with sparsely encoded wavelet domain data in a compressive sensing framework

Abstract. Randomly encoded compressive sensing (CS) has potential in fast acquisition of magnetic resonance imaging (MRI) data in most naturally compressible images. However, there is no guaranteed good performance for general applications by any of the traditional CS-MRI theoretical schemes developed so far. On the other hand, recent research demonstrates that adaptive sampling exploiting the tree structure of nonzero wavelet coefficients of signals allows more control over the sensing procedure in the form of feedback and improve the CS performance significantly. Following recent research strategies in CS-MRI, well-known adaptive sampling strategies in the wavelet domain, as used in image compression, to encode the MRI data yielding good reconstruction quality are introduced. Based on this underlying characteristic, adaptive k-space trajectories are designed with tailored spatially selective RF excitation pulses created by Battle-Lemarie wavelet functions. The input vectors formed from these significant samples of multilevel wavelet decomposed images are used in a CS framework for reconstruction of MR images. This MR image reconstruction uses a CS algorithm based on the minimization of total-variation regularized signal to provide stable results. The simulated results show that this approach can reduce almost 70% of MR image acquisition time and achieve good reconstructed image quality.


Introduction
In conventional magnetic resonance imaging (MRI), scanners acquire samples of the encoded image in a spatial frequency (Fourier) domain, which is called k-space. The MR image of the object in the spatial domain can be reconstructed through inverse Fourier transform of the k-space data. 1 A well-known problem in MRI is its long scan time that depends on the required number of acquired samples determined through Shannon-Nyquist sampling theory. Potential reconstruction of MR images from a reduced number of acquired samples without degrading the image quality would offer an effective method to reduce scan time while improving the resolution of current MR imagers. However, insufficient sampling of k-space data violates the Shannon-Nyquist criterion and produces aliasing artifacts and noise in reconstructed MR images. Several common approaches to alleviate these artifacts and noise are to exploit redundancies in k-space and to use temporal filtering, such as partial-Fourier, variable-density sampling, unaliasing by Fourier-encoding the overlaps using the temporal dimension (UNFOLD), etc. [2][3][4][5] However, these methods still do not overcome low sampling rate artifacts. Investigations on the feasibility of wavelet encoding in MRI demonstrated better flexibility and adaptability in MR data acquisition than Fourier encoding. [6][7][8] Unfortunately, the practicality of wavelet encoding in faster acquisition of MRI data has not been demonstrated in the above efforts.
Through recent developments in the theory of compressive sensing (CS), MR images with a sparse representation in some transform domain, under the constraint of incoherent measurement, can be recovered from "randomly" under sampled k-space. 9,10 CS theory here exploits such sparsity of MRI data and reconstructs an MR image from very few incoherent measurements through a nonlinear procedure. Several methods have applied CS to Fourier-encoded MRI to reduce scan time. [10][11][12][13][14][15][16][17] These methods are based on two characteristics of conventional MRI: MR images are naturally compressible in certain transform domains, and the sampling pattern is incoherent with respect to some sparsifying transformations. The natural fit of CS to MRI is discussed 10 by reviewing the constraints required for successful reconstruction. Lustig et al. 11 have applied CS to recover MR images from a subset of Fourier-encoded k-space acquired by pseudo-random variable-density under sampling of phase-encoding to achieve fast imaging. The discrete cosine transforms (DCT), wavelet transform and finite-difference transform were used to exploit the sparsity of brain and angiogram MR images. The incoherence between these sparsifying transforms and Fourier operator were measured by the point spread function (PSF), and transform PSF (TPSF), whereas variable density random under sampling was used to improve the degree of incoherence. As their results showed, this CS approach achieved high reduction of imaging time in three-dimensional (3-D) Fourier-encoded imaging with good reconstruction performance. However, in two-dimensional (2-D) Fourier-encoded MRI, this approach does not achieve such benefit because only one-dimensioanl (1-D) sparsity was exploited.
In mathematical CS theory, Fourier basis is maximally incoherent with the canonical basis. 9,18 In traditional MRI, because the k-space data are encoded by Fourier basis, the incoherent measurement constraint is only satisfied when the MR signal is sparse in the spatial domain. Consequently, it can be seen 11 that CS-MRI is more suitable for angiogram imaging than for brain imaging because angiograms are sparse in their pixel representation. Therefore, Fourierencoded MRI is not a good encoding scheme for stable and accurate CS reconstruction in most scenarios where MR signals are not sparse in the spatial domain. In addition, some existing crucial limitations remain unresolved in Fourierencoded MRI, such as optimization of sampling trajectories and the computation time of reconstruction algorithms.
Several studies have also investigated non-Fourier MRI and CS schemes to reduce the image acquisition time while improving the reconstruction quality. [19][20][21] Random encoding of k-space is achieved by tailoring spatially selective RF pulses to satisfy CS constraints in non-Fourier domains, 19 achieving impressive reconstruction results. The CS constraints are satisfied better by making the energy of kspace spread out in non-Fourier domain and pseudo 2-D random sampling. 20,21 Despite the potential of the CS approaches in non-Fourier MRI to outperform traditional Fourier MRI, no general CS approach has yet been developed with strict theoretical justification.
Recent research on application of CS to fast MRI suggests that there is no theoretical performance guarantee for CS-MRI with or without random encoding, and new CS based methods should be developed for specific applications. 19 In recently published CS papers, adaptive sampling in the wavelet domain has been proposed to improve the signal recovery performance. 22,23 Adaptive sampling schemes exploiting the tree structure of nonzero wavelet signal coefficients have been used to replace the 'universal' acquisition of random or pseudo-random sampling in traditional CS methodology. 22,23 These schemes use nonrandom sampling and allow more control over the sensing procedure in the form of feedback to improve the CS performance significantly with fewer measurements and higher reconstruction quality. Haupt et al., have demonstrated that the adaptive sampling in CS, called compressive distilled sensing, has the advantage of significantly improving the error bounds compared to traditional CS schemes without adaptive sampling. 24 These adaptive strategies have the potential to achieve more accurate and robust signal recovery in some practical applications. 25 Therefore, if adaptive sampling is used in wavelet-encoded CS-MRI, it is possible to improve the reconstruction quality of images. We introduced well-known adaptive sampling in the wavelet domain, as used in image compression, [26][27][28] to encode the MRI data yielding good reconstruction quality.
In this paper, an efficient implementation of a specific CS approach in a simulator that allows good MRI reconstruction from a sparse wavelet-encoded k-space following the wellknown embedded zero-tree structure for image coding is presented. [26][27][28] With this approach, sparsity of image data in the spatial domain is not required. In wavelet-encoded MRI, the k-space data are encoded in multiple levels in the wavelet domain. 6 The multilevel decomposed image is then transformed into a vector of sparse wavelet coefficients, and these few significant coefficients are encoded to reconstruct the image with high fidelity. 26 Therefore, it is possible to represent the k-space data with only a few significant samples and still reconstruct an MR image with a reduced scan time. Additionally, the reconstruction quality of such under sampled k-space may be guaranteed by using the input vector of sparse wavelet coefficients in a specific CS framework. The following sections will discuss the proposed CS-MRI approach in the wavelet domain to achieve fast imaging and assess the reconstruction performance at different sampling rates, noise and sparsity levels. Section 2 introduces the procedure of sparse acquisition in wavelet-encoded MRI through the description of RF pulse definition, pulse sequence implementation and adaptive k-space trajectories in a simulator. Section 3 presents an efficient CS algorithm based on the minimization of total-variation regularization and a least squares measurement of signal to reconstruct MR image from these sparsely encoded k-space samples in the wavelet domain. Section 4 discusses the simulation results and reconstruction performance for phantom and brain images. Section 5 describes feasibility of future applications of our work.

Implementation of Sparse Data Acquisition in
Wavelet-Encoded k-Space

Sparsity of Image/Signal in Wavelet Domain
The wavelet transform of natural images produces a large number of coefficients of values with zero or near-zero magnitudes and a small number of significant values (e.g., the ones with larger magnitude than a determined threshold T). If we localize and encode these significant values, it is possible to reconstruct images with high fidelity. In addition, wavelet transform provides a compact multiresolution representation of the image, and the significant wavelet coefficients are well-organized in hierarchical trees. Thus, the high-resolution detail tends to be significant only if significant details can be found at all resolutions from the lowest level. This property is also commonly used to predict the positions of significant information across scales in wavelet-based compression algorithms. [26][27][28] Therefore, for many natural or medical images, it is possible to represent the wavelet-encoded k-space using only a few tree structured significant samples. These samples capture most of the signal information and can be used to reconstruct the MR image with insignificant distortion. This sparse acquisition of wavelet-encoded k-space data is achieved by designing spatially selective RF excitation pulses for generating adaptive k-space trajectories with reduced scan time.
In particular, assume that x-axis is the wavelet encoding direction and y-axis is the frequency encoding direction. In each subspace, the resolution of the y-axis is constant, and the resolution of the x-axis is varied with the decomposition level j. In the j-level subspace, the resolution along the x-axis can be calculated by n j ¼ N · 2 −j , where N is the highest resolution when j ¼ 0. In order to simplify the discussion, we useρðxÞ to represent one line of fully sampled frequency encoding at the location x, e.g.,ρðxÞ ¼ P y ρðx; yÞ expð−i2πk y · yÞ. Therefore, the wavelet-encoded MRI can be considered as wavelet transform of one dimension signalρðxÞ. For level J, each subspace is constructed by wavelet coefficients as follows: Here, φ −J;l ðxÞ ¼ 2 −J∕2 φ 0;0 ð2 −J x − lÞ, and ψ −j;k ¼ 2 −j∕2 ψ 0;0 ð2 −j x − kÞ, j ¼ 1; 2; : : : ; J. The coefficient setĀ −J in V −J in terms of the scaling function φ −J;l ðxÞ represents a coarse approximation of the signal. The coefficient sets {D −j , j ¼ 1; 2; : : : ; J} in terms of the wavelet function {ψ −j;k ðxÞ, j ¼ 1; 2; : : : ; J} represent the detail or difference between the coarse approximation and a finer resolution approximation. From Eqs. (2) and (3), we observe that the sparse acquisition of a signal can be achieved by tailoring spatially selective RF excitation pulses.
A hierarchical binary tree structure exists in the wavelet coefficients set, where all wavelet coefficients can be arranged by a "parent-child" relationship. 28,29 The coefficient inD −J is called "parent," and all coefficients corresponding to the same spatial location in {D −j , j ¼ 1; 2; : : : ; J − 1} are called "child." If the "parent" is insignificant, all its "children" are insignificant. All "parents" are tested first for significance, and then the positions of all significant "children" are contained in their trees. It is proposed that the wavelet-encoded MRI data be acquired starting with the lowest resolution k-space W −J and then the next finer scale k-space {W −j , j ¼ 1; 2; : : : ; J − 1} be undersampled based on the tested results in W −J to acquire only significant samples. The support of coefficient setD −J is defined to be the positions where its coefficientsd −J;k are larger in magnitude than an application dependent threshold T. The support ofD −J is thus represented by Here, T ≥ 0 determines the number of samples in k-space.

2.3
Simulation of Sparse MR Data Acquisition in the Wavelet Domain Most MRI simulators are developed based on the solutions of the Bloch equation, and the generated MR signal is encoded in Fourier domain. 30,31 Wavelet domain simulators were developed in earlier years but did not demonstrate any advantage in MRI scanning time. [6][7][8] In this paper, a simulator generating MR signal in wavelet domain using adaptive sampling in a CS framework has been developed by designing RF excitation pulses and pulse sequences as described in detail in the following sections. This simulator generating MR signals in the wavelet domain is developed for sparse data acquisition as shown in Fig. 1

Simulator overview
An overview of the proposed MRI simulator to acquire sparse wavelet-encoded k-space data is shown in Fig. 1. In the proposed data acquisition scheme, we do not need the information from fully sampled data. We actually acquire the wavelet-encoded MRI data starting with the lowest resolution k-space W −J by designing spatially localized wavelet shaped RF excitation pulses for generating kspace trajectory by replacing the traditional phase encoding trajectory. Only the lowest resolution trajectory is fully sampled (requiring a small number of RF excitation pulses and the least data acquisition time) to identify for the significant parent wavelet coefficients and to allow full sampling across the levels if desired. The subsequent finer scale k-space {W −j , j ¼ 1; 2; : : : ; J − 1} data are generated by a multiple (e.g., two to eight for four level wavelet encoding based on the parent-child tree structure of wavelet coefficients) of RF excitation pulses as shown in the encoding overview in Fig. 1 (see Sec. 4.2 for details of implementation). Using a virtual object, scanning parameters, such as B 0 , G x , G y (magnetic gradient along x and y axis), flip angle, TE (time echo), TR (time repetition), are initialized first. This simulator generates k-space data in multiple levels. For J-level imaging, the level j of k-space acquisition starts from J to 1. When j equals to J, two RF excitation pulses, which are designed, respectively by the J'th-level scaling and wavelet functions, are used to generate spin echo signals. After wavelet-encoding and frequency-encoding, the subspaces V −J and W −J are acquired, and then significant samples are tested in W −J to construct {suppðD −j Þ, j ¼ 1; 2; : : : ; J − 1}. When level j is from J − 1 to 1, RF excitation pulses are designed by j'th-level wavelet function and indexed by supp (D −j ) to acquire under sampled k-space W -j . Finally, the sparse-encoded k-space S is constructed by these subspaces V −J , and W −J to W −1 to contain all significant samples. To simulate realistic images, noise is added to the k-space.

Design of RF excitation pulses
Excitation profiles in this wavelet-encoded MRI simulator are varied and shaped by functions of wavelet transform. At small flip angles (less than 30 deg), RF pulses are shaped by Fourier transform of excitation profiles. Here, we utilize the Battle-Lemarie wavelet basis functions as shown in Fig. 2 to design RF pulses. Figure 2(a) shows the wavelet scaling and basic functions φ 0;0 ðxÞ and ψ 0;0 ðxÞ in the spatial domain. The major reason to choose this specific basis function set is that their Fourier transforms are smooth and rapidly decay to zero as shown in Fig. 2(b). This basis set provides short RF pulses with relatively precise spatial excitation profiles. 6 For J-level wavelet encoding, assume that the gradient strength along x-axis is G x , the FOV along x-axis is X centimeters, and the highest resolution is N. According to Eqs.      Fig. 3. After reconstruction of MR images in the spatial domain, the imaging resolution of 1 mm can be achieved. Any deviation of G x (e.g., B 1 inhomogeneity along xaxis) would result in a distortion of excitation profiles, especially occurring in high static field (B 0 ≥ 3T) and with tailored RF pulses. In practice, this distortion can be corrected using a prescan with the same sequence. A fully sampled image acquired for each spatially selective pulse is used as a standard to compute compensation parameters for each excitation profiles. 19 The inhomogeneity maps generated from 3-D images are acquired by prescans and the carrier frequency and amplitude of pulses are adjusted to excite a uniform slice. 32 In this paper, if the inhomogeneity maps (e.g., G 0 x ) is acquired using prescans, the distortion of excitation profiles can be calibrated by adjusting the carrier frequency of each pulses. In the simplified version of this simulator, we neglect the inhomogeneity in G x and there is no distortion of excitation profiles.

Pulse sequences
During MR imaging, the timing of pulse sequences determines the RF signal acquisition and k-space trajectories. Pulse sequences contain RF pulses and magnetic field gradients. In wavelet-encoded MRI simulator, the timing of pulse sequences is defined with spatially selective RF excitation and adaptive k-space trajectories to yield a sparse encoding scheme.
As shown in Fig. 1, at the beginning, RF pulses {Φ −J;k , k ¼ 1; 2; : : : ; N · 2 −J } and {Ψ −J;k , k ¼ 1; 2; : : : ; N · 2 −J } with the carrier frequency offset by {k · 2 J Δω x } are used to generate J-level wavelet-shaped excitation profiles along x axis. In the following step, a precession with application of magnetic gradient along y axis is specified by its duration and the gradient magnitude to fill the k'th line of frequency-encoding in k-space. After N · 2 −J RF pulses and lines of frequency-encoding, the k-space V −J and W −J are fully sampled. In W −J , the positions of significant samples, suppðD −J Þ, are tested by Eq. (4) and then {suppðD −j Þ, j ¼ 1; 2; : : : ; J − 1} are constructed. When j is from J − 1 to 1, RF pulses {Ψ −j;k , k ∈ suppðD −j Þ} with the carrier frequency offset by {k · 2 j Δω x } and magnetic gradient are applied to achieve the trajectory consisting only of significant samples in W −j . Finally, the sparsely encoded k-space S constructed by these subspaces V −J , W −J to W −1 contains the significant samples.

Noise
In real MRI scanners, the acquired raw k-space data are considered to be corrupted by complex white Gaussian noise with the same variance in the real and imaginary parts. After reconstruction of MR images from the corrupted kspace data, the noise embedded in the image intensity has Rayleigh statistics in the background and Rician statistics in the signal regions. 33,34 Therefore, this noise model is used in this simulator with different signal to noise ratio (SNR), with SNR is defined as where A signal is the average root mean square (RMS) amplitude of the signal, and A noise is the average RMS amplitude of the noise.

CS Scheme in Wavelet-Encoded MRI
In the proposed wavelet-encoded MRI simulator, if k-space data are fully sampled based on Shannon-Nyquist sampling theory, an MR image could be reconstructed by inverse wavelet transform along the wavelet-encoding direction and inverse Fourier transform along the frequency-encoding direction. If k-space data are undersampled, this method cannot reconstruct MR images accurately. CS theory offers a potential to reconstruct a compressible signal from far fewer samples than that required by Shannon-Nyquist sampling theory and provides a new approach to overcome conventional limitations in signal sampling. 9,18 Therefore, CS may be applied to reconstruct MR image from under sampled k-space. The reconstruction quality of such under sampled k-space is guaranteed by using the input vector of sparse wavelet coefficients in a specific CS framework. Because the frequency encoding is fully sampled, ifρ is recovered, MR image ρ can be reconstructed by inverse Fourier transform ofρ. The main task of CS reconstruction is to reconstructρ from the sparse wavelet-encoded k-space replacing the traditional phase encoding direction. According to Eq. (2) and (3),Ā −J andD −j can be represented asĀ where Φ −J is the N2 −J × N matrix whose rows are φ −J;l ðxÞ for l ¼ 1; 2; : : : ; N2 −J , and Ψ −j is the N2 −j × N matrix whose rows are ψ −j;k ðxÞ for k ¼ 1; 2; : : : N2 −j . We denotē A ¼Ā −J ,D ¼ ½D −J ;D −Jþ1 ; : : : ;D −1 , Φ ¼ Φ −J and Ψ ¼ ½Ψ −J ; Ψ −Jþ1 ; : : : ; Ψ −1 . Therefore, the measurement y of wavelet-encoded k-space can be represented as where E is the sensing matrix, I is the identity matrix, F ∈ R M×ðN−N2 −J Þ is a rectangular matrix that selects M elements of any (N − N2 −J )-dimensional vector (M < N − N2 −J ) and is determined by {supp suppðD −j Þ, j ¼ 1; 2; : : : ; J}. We consider the global sensing matrix A and the sparse wavelet coefficients set α in Eq. (8) to be Because the sensing matrix E involves partial wavelet transform, it is not suitable to exploit the sparsity ofρ by wavelet coefficient set α due to the "incoherent" constraint in the CS framework. Most real MR signals are approximately piecewise smooth and have a small total variation (TV), 11,19 which is defined as TVðxÞ ¼ . The use of TV regularization could promote sparsity of finite differences in MR images and reduce ringing artifacts near edges that are caused by tailoring of the wavelet coefficients. 8 Therefore, the reconstruction ofρ, denoted asρ, can be computed by solving the following problem: 35 ρ ¼ arg min TVðρÞ s:t: kAα − yk l 2 ≤ ε 2 ; where l 2 -norm is defined as kxk l 2 ¼ ð P i jx i j 2 Þ 1∕2 and parameter ε 2 controls the fidelity ofρ to the measurements y and is determined by estimating the noise variance. Nesterov's algorithm can be used to solve the Lagrangian form of (10) to achieve the reconstruction by 36 minimize λTVðρÞ þ Aα − y l 2 ; where λ > 0 is used to balance TV regularization, and the measurement consistency is experimentally set. Large λ leads to staircase effect. To solve the CS problems in Eq. (11), the CVX software package by Grant, Boyd, and Ye was used (http:// www.stanford.edu/~boyd/cvx). Although the CVX package is a slow code, we can acquire the raw k-space data first and then implement reconstruction off-line. This procedure still reduces the scanning time and provides the potential to achieve fast imaging.

Reconstruction Stability
In order to guarantee the performance of reconstruction by Eq. (10), the global sensing matrix A in Eq. (9) must satisfy the restricted isometry property (RIP). 37,38 If A satisfies the RIP of order k, for all k-sparse vectors α, where δ k is the isometry constant, which should be a small number to satisfy RIP constraint. Because the degree of coherence may be considered as isometry constant of order 2, the RIP-based guarantees are generally stronger than incoherence-based guarantees for stable and accurate reconstruction. 39 The reconstruction performance is improved as δ k is made smaller. For example, from recent CS papers, [38][39][40] if δ 2k < ffiffi ffi 2 p − 1 in a noiseless environment, the k-sparse vector α can be reconstructed exactly by Eq. (11) with ε 2 ¼ 0 with fewer than k nonzero entries. However, given any matrix A, it is often infeasible to compute practically useful RIP-based guarantees for all sparse vectors α due to computational problems. Therefore, we investigate the RIP guarantees for A, which is designed adaptively based on a sparse vector α.
In Eq. (8), the measurement y is divided into two parts y 0 and y 1 , where y 0 ¼Ā is the fully sampled scaling coefficient set and y 1 ¼ FD is the under sampled wavelet coefficient set D. Therefore, the measurement consistence constraint kAα − yk l 2 ≤ ε 2 equals to kFD − y 1 k l 2 ≤ ε 2 , and the RIP condition in Eq. (12) can be rewritten as Larger values of kFDk l 2 result in smaller δ k . For example, if kFDk l 2 ¼ kDk l 2 , δ k ¼ 0. Thus, all nonzero entries inD are contained in the under sampled k-space y, andρ can be recovered exactly. Therefore, sparse acquisition of waveletencoded k-space data can result in a smaller value of δ k than other encoding methods by including as many significant samples as possible in such under sampled k-space. In addition, for each MR image, specific global sensing matrix A is constructed adaptively by investigation of sparsity in wavelet-encoded k-space.

CS Reconstruction Experiment
The numerical experiments for CS reconstruction of 1D piecewise smooth signal {x½n, n ¼ 1; 2; : : : ; 256} are performed by 4-level sparse wavelet encoding. The proposed CS reconstruction result from significant wavelet coefficients only, denoted as x Ã , is compared with the best sparse approximation result, denoted as x 0 , which includes the zero-filled insignificant coefficients as well. Figure 4 shows the reconstruction results x Ã and x 0 with different measurement number M ¼ 60, 88, 116 and 144. The measurement y at each sampling number is acquired by sparse encoding in the wavelet domain. CS reconstruction performs better than the best sparse approximation reconstruction does with x Ã approaching the original signal x closely under all conditions. Table 1 shows RIP constants and reconstruction SNRs at different measurement number M. Larger values of M result in smaller values of δ k and higher reconstruction SNRs by CS with better sparse approximation results.

Experimental Results
The Shepp-Logan phantom and a 3-D digital brain phantom from the McConnell Brain Imaging Center, Montreal Neurological Institute, McGill University are used as virtual objects in the MRI simulator to investigate the properties of the proposed wavelet-encoded CS-MRI scheme (WCS-MRI). 41,42 The digital brain phantom has a spatial resolution of 1 mm 3 and contains most of the relevant tissue types. Based on these phantoms, we first simulate the fully sampled wavelet-encoded MRI, and then we simulate the sparse acquisition of k-space data at different sampling rates and reconstruct MR images by CS. These simulation results are compared with the results of the commonly used Fourier-encoded CS-MRI scheme (FCS-MRI). In Fourierencoded k-space, most of the energy is concentrated close to the center and rapidly decays toward the periphery. Therefore, a variable-density sampling scheme is implemented in FCS-MRI, where the phase encoding locations are randomly spaced but cover the low-frequency portion near the center of k-pace, matching the energy distribution in k-space. The CS reconstruction of MR image is performed from the undersampled k-space data.

Implementation of Wavelet-Encoded MRI
Simulation Data are collected for MR image reconstruction on a 256 × 256 voxel grid using the Shepp-Logan phantom and two different brain slices (slice 1 and slice 2) from the brain phantom. Figures 5 and 6 show fully sampled Fourierencoded and 4-level wavelet-encoded MRI results, respectively. These results are generated with zero noise and homogeneous magnetic field (B 0 ¼ 3T). The imaging parameters are TE∕TR ¼ 25∕500 ms, read out time ¼ 10 ms and flip angle ¼ 30 deg. For each image, there are 256 RF excitation pulses required, and its spatial resolution is 256 × 256. It is seen visually that wavelet-encoded MRI acquires images with the same quality and resolution as Fourier-encoded MRI does. The quantitative assessment of the quality of the reconstructed images has been computed by peak signal to noise ratio (PSNR) and structural similarity (SSIM) as shown in Sec. 4.2. Figure 7 shows the energy of k-spaces acquired by Fourier-and wavelet-encoded MRI for brain image 1, respectively. In Fig. 7(a), most of the energy in Fourierencoded k-space is concentrated close to the center of  k-space, at low spatial frequencies. However, in Fig. 7(b), the wavelet-encoded k-space is represented in multiple levels with a high degree predictability between levels. The coarse scale subspace V −4 has its spectrum localized in the low frequency band, while these fine scale subspaces W −1 to W −4 have their energy widely dispersed in the high frequency band and well-organized in trees along wavelet encoding direction.

Simulation of CS-MRI
The proposed WCS-MRI scheme has been simulated with the phantoms shown in Fig. 5 with different sampling rates. The implementations of sampling patterns reported [22][23][24] are not the same as our approach. Our data acquisition time is reduced by sparse acquisition of only significant samples across multiple levels by designing appropriate RF excitation pulses at each level. Such wavelet-encoded MRI provides the flexibility and adaptivity to represent the k-space data, especially in multiple levels.
According to the wavelet decomposition theory, it is possible to achieve the sparse acquisition of k-space data by exploiting the hierarchical binary tree structure in the wavelet coefficients set. [6][7][8]28 Our implementation is based on the structural characteristics of wavelet coefficients and uses a desired percentage of the full samples in our proposed CS-MRI making the acquisition time directly proportional to the sampling rate. Our data acquisition time is reduced by sparse acquisition of only significant samples across multiple levels by designing appropriate RF excitation pulses at each level. For 4-level wavelet-encoded MRI, there are 16  , the corresponding number of RF pulses in W −3 , W −2 and W −1 are 2m, 4m and 8m, respectively. Therefore, the total sampling rate (SR) can be calculated by ð32 þ 14mÞ∕256, including 32 samples from W −4 and V −4 . In this simulation, if m ¼ f1; 2; : : : ; 8g, the corresponding sampling rate SR ¼ f18%; 23%; 29%; 34%; 40%; 45%; 51%; 56%g. The scan time is determined by SR· T total , where T total is the imaging time in fully encoded MRI with the same resolution. Finally, MR images are reconstructed from these under sampled k-spaces through CS. In order to compare with WCS-MRI, the same number of RF pulses is used to acquire k-space data in FCS-MRI. To make the simulation close to real MRI acquisition, Gaussian noise is added. Figures 8 and 9 show the reconstruction of test images by WCS-and FCS-MRI schemes without noise; sampling rates are 23%, 40%, and 56% (e.g., 60, 102 and 144 RF excitation pulses used), respectively. Large-scale ringing artifacts are apparent at low sampling rate (e.g., SR ¼ 23%) in the FCS-MRI results. However, artifacts in WCS-MRI are much less than those in FCS-MRI, and more detail and high frequency components are kept in WCS-MRI than in FCS-MRI. When the sampling rate is 40%, no artifacts are apparent in WCS-MRI, while artifacts are still visible in FCS-MRI. When the sampling rate is 56%, both reconstruction schemes perform well without apparent artifacts. Because the Shepp-Logan phantom image is much sparser in pixels than brain images, CS reconstruction of the Shepp-Logan phantom image can achieve a much better performance.

Simulations under noisy conditions
Simulations are also performed to illustrate reconstruction performance with noise. Figures 10 and 11 show low-SNR simulation results (e.g., SNR ¼ 10 dB) of the Shepp-Logan phantom reconstruction and high-SNR simulation results (e.g., SNR ¼ 30 dB) of brain image 2 reconstruction by WCS-and FCS-MRI schemes with sampling rates 23%, 40%, and 56%. In Fig. 10, FCS-MRI scheme is unable to reconstruct MR images exactly at low sampling rate. However, the performance of MR images by WCS-MRI is better than that by FCS-MRI. In Fig. 11, WCS-MRI can reconstruct MR images without apparent artifacts and loss of high-resolution components, while in FCS-MRI, there are much apparent artifacts and loss of high-resolution components in reconstruction results. As expected, improved SNR leads to improved reconstruction quality for these two schemes. These figures indicate that WCS-MRI performs better than FCS-MRI regardless of noise level.

Quantitative evaluation of image reconstruction performance
The fully sampled 256 × 256 Fourier-encoded MR images in Fig. 4 are used as the gold standards to evaluate the image reconstruction performance. In order to measure the similarity between the reconstructed imageρ and the gold-standard image ρ ref , the PSNR and SSIM index are used. 43,44 PSNR is the most widely used image quality assessment metric and has clear physical meaning The SSIM metric is used to measure visual reconstruction quality by capturing the similarity between the original image and the reconstructed image. SSIM models any distortion as a combination of three different factors: loss of correlation, luminance distortion, and contrast distortion. The dynamic range of SSIM is ½−1; 1. SSIM represents the degree of similarity between ρ ref andρ. The higher the SSIM is, the more similar ρ ref andρ are. The maximum value of 1 is achieved only if ρ ¼ ρ ref . Figures 12 and 13 show PSNR and SSIM as functions of data sampling rate in the ideal condition, respectively. As shown in Figs. 12 and 13, with SR ¼ 29%, WCS-MRI can achieve very good imaging quality with PSNR > 30dB and SSIM > 0.975 in all tested images. In low sampling rates (i.e., SR < 29%), WCS-MRI outperform FCS-MRI with higher PSNRs and SSIMs. Even with SR ¼ 18%, PSNR in WCS-MRI can reach more than 25 dB and SSIM more than 0.9. With all tested sampling rates and images, PSNRs and SSIMs are higher in WCS-MRI than in FCS-MRI. Figures 14 and 15

Improvement in CS reconstruction
The use of sparse-encoding in wavelet-encoded k-space in this work was motivated by the desire to improve RIP constants to achieve good CS reconstruction. It is difficult to compute RIP constants for all signal construction, but we can calculate them in the FCS-MRI scheme by Eq. (12) and in WCS-MRI scheme by Eq. (13) for each tested MR image. Table 2 shows the value of RIP constant δ K (K is determined by SR) for each test image. When the value of δ K is smaller, the reconstruction performance is better in the corresponding image. The value of δ K in WCS-MRI is less than that in FCS-MRI for most conditions. This implies that WCS-MRI provides better reconstruction stability and accuracy than FCS-MRI.      the performance of MR image reconstruction in ideal and noisy conditions. In the case of low-SNR and low sampling rate, WCS-MRI achieved much higher values of PSNR and SSIM than FCS-MRI. In the case of high-SNR and high sampling rate, FCS-MRI and WCS-MRI both achieved similar values. That means WCS-MRI has advantages in reduction of scan time. When the sampling rate is more than 30%, WCS-MRI reconstruction achieves SNR ≥ 30 dB. The acquisition time of MR images can be reduced significantly (by almost 70%) without much distortion when the proposed WCS-MRI method is applied.

Conclusions
This work describes the development of a simulator for MRI with a wavelet-encoding scheme for stable and accurate MR image reconstruction by acquiring subsampled data below the Nyquist sampling rate following the embedded hierarchical tree structure of significant wavelet coefficients within a CS framework. According to our simulation results, the proposed scheme improves the existing Fourier-encoded CS-MRI method by introducing this specific wavelet encoding instead of commonly used random encoding in the k-space data acquisition. The hierarchical wavelet tree structure is applied to select sparse encodings in order to improve the precision of the CS reconstruction by satisfying the requirement of RIP. The acquired under sampled k-space contains many significant samples and the reconstruction quality of such under sampled k-space is guaranteed by using the input vector of sparse wavelet coefficients in a specific CS framework, which is based on the minimization of totalvariation (TV) regularization signal and a least squares measurement. Based on the flexibility and adaptivity in waveletencoded MR data acquisition, it is feasible to implement sparse encoding of k-space data with tailored spatially selective RF excitation pulses without much modification of MRI machinery. Therefore, this proposed encoding scheme may provide a practical method for shortening the patient scan time by reducing the required number of RF excitation pulses without decreasing the resolution of reconstructed MR images. In addition, this work may be applied in future to functional MRI (fMRI), where the three-and four-dimensional data would be even more compressible than the two dimensional data used in the wavelet transform domain. For example in echo planar imaging (EPI), using wavelet-shaped slice-selection excitation pulse with sparse encoding of slice data may accelerate fMRI with improved CS reconstruction.