MRI Simulation-based evaluation of an efficient under-sampling approach

1.   Introduction

In the medical imaging field, MRI has revolutionized the diagnosis of diseases through images, based on the magnetic resonance phenomenon of hydrogen nuclei in tissues of imaged objects. In principle, objects are stimulated with radio pulses (RF) and resonator signals are received using RF coils. Rapid imaging in MRI is an important issue to improve the quality and resolution of images, to avoid the physiological effects on patients or to meet the time requirements when the imaged structures are dynamic.

Lauterbur proposed to use the concept of k-space since 1973 for interpreting MR acquisition as Fourier encoding in 2D or 3D spaces ^[1]. In particular, information gathered by an MRI scanner is samples of the spatial Fourier change of an image. Subsequently, so as to get an image without aliasing artifacts, k-space samples need to fulfill the Nyquist sampling criterion. The development of parallel MRI (pMRI) enabled to speed up the image acquisition of MRI by utilizing an array of RF coils for simultaneously receiving MR signals ^[2,3]. Many studies tried to improve the pMRI with novel reconstruction algorithms such as the Generalized Autocalibrating Partially Parallel Acquisition (GRAPPA) ^[4], the multiphasic contrast-enhanced volume-interpolated sequence ^[5], a central elliptical cylinder in k-space repeating n times (keyhole) with a random acquisition (CENTRA) ^[6], a combination of contrast enhanced time robust angiography, keyhole, and viewsharing techniques ^[7], and a combination of pseudorandom sampling and temporal viewsharing ^[8].

Besides the pMRI, in the field of signal processing and information theory, there is a breakthrough that is compressed sensing (CS)-MRI ^[1,9], which indicates that sparse and compressive signals can be recovered from a tiny amount of random measurements. CS has been successfully applied for MRI fast acquisition in many studies ^{[10,11,12,13,14]}. The common strategy of these approaches is to exploit the redundancy of image data for reducing the sampling rate. Thus, the CS-MRI enables fast image acquisition by reducing the scanning time, so that it offers many advantages such as reducing patient burden, motion artifacts, and contrast washout ^[15]. Since the first introduction of CS by Lustig et al. ^[14], many other studies have been developed using various sparsifying transforms and optimization algorithms. In the traditional CS-MRI methods, random under-sampling approaches with fixed compression ratios are simply performed for obtaining the sampled horizontal lines in the binary mask, which are completely based on the power law ^[12]. Interest signals can be reconstructed using nonlinear approaches such as basic pursuit (l₁-BP) ^[1], orthogonal matching pursuit (OMP) ^[16], and non-linear gradient conjugate (NGC) ^[14,17]. Difference from the simple conventional CS methods, recent the state-of-the-art CS methods have been reported with advanced reconstruction algorithms such as the structured Hankel matrix for solving the image quality degradation issues ^[18], the deep de-aliasing generative adversarial networks (DAGAN) can reduce aliasing artifacts ^[19,20] and the stochastic deep CS for improving the reconstruction of diffusion tensor cardiac MRI ^[21]. Although these advanced CS approaches offer undeniable advantages, the highly computational steps of the reconstruction algorithms are still required. The developments of new CS techniques utilizing the well-developed reconstruction algorithms without increasing the computational complexity are still highly recommended, which will ease the clinical workflow in real MRI systems.

More recently, other advanced techniques based on the artificial intelligence (AI) for improving MRI quality have been also raising attention from science societies. From the fast development of Generative Adversarial Network (GAN), many studies based deep learning techniques have been performed to obtain super-resolution quality for real-time MRI reconstruction ^{[22,23,24,25,26]}. However, in comparison with the CS, these advanced methods are considered as completely different technology, which are mainly based on highly computational complexity of AI algorithms. Thus CS-MRI techniques have still their own merits if more effort has been investigated.

In this study, we propose a new hybrid random under-sampling approach for improving the CS-MRI, which combines the conventional random under-sampling with the definite sampling. With a same compression ratio, the total number of sampled horizontal lines is divided into two parts: (ⅰ) the large first part is still based on the power-law; (ⅱ) the small second part is enhanced with the remaining lines that are near the origin of k-space. Because the amount of encoded information is concentrated at the origin of the k-space, the proposed method suggests that the amount of useful information will be collected more and therefore the MRI image reconstruction will be more accurate. Our results based simulation studies have confirmed the significant improvement of image reconstruction fidelity from the original under-sampling methods. The utilization of well-developed reconstruction algorithms in this proposed under-sampling technique also avoid the highly computational complexity as required in many state-of-the-art CS methods.

2.   Methodology

2.1.   MRI image acquisition

It is supposed that m(x, y) is a 2D image of a slice in the interest object. The analogue obtained signal by a collected coil in k-space is expressed as follows:

where k_x and k_y are locations of the encoded information in axes x and y of the interesting slice. Using 2D-Fourier transform, m(x, y) can be obtained. The discrete form of Eq (1) is shown as follows:

where N_x and N_y are the total numbers of pixels along axes x and y of the slice, respectively. In this paper, the Cartesian trajectory is used for 2D imaging and the encoded information density of the k-space is followed by the power-law (Figure 1). Figure 1a presents the k-space of an MR image of a brain slice. It can be seen that higher amplitudes of signals mainly concentrate around the origin of the k-space as shown by the bright domain at the center, which indicates more information is distributed at the origin. In practice, the most encoded information is concentrated at the origin due to the low frequencies of phase coding steps and the density of the k-space follows a power law ^[12,27]. A full Cartesian sampling of the k-space followed the Nyquist criterion is indicated by red dots (samples) as in the Figure 1a. By applying the 2-D FFT to the signal in k-space, the MR image can be reconstructed as shown in Figure 1b.

2.2.   Fundamental compressive sampling

Supposed that $x \in {R^N}$ is the interest signal and it has a sparse linear representation in some domain with x = Φs, where $s \in {R^N}$ and $\Phi \in {R^{N \times N}}$ are respectively a L-sparse vector (the exact number of nonzero values in s is L) and the sparsifying matrix. In addition, x is assumed to be sensed by using a linear system $\Psi \in {R^{M \times N}}$ , and then, the acquired measurements, $y \in {R^M}$ , are defined by y = Ψx ^[28]. To recover x from y, which also corresponds to the recovery of s from y, because the obtained measurements can be expressed as y = Θs, where Θ = ΦΨ as shown in Figure 2 illustrating the basic principle of compressive sampling technique.

In compressed sensing, measurement matrix, Ψ or Θ, is commonly underdetermined. One of the important conditions of CS is the restricted isometry property (RIP), which allows the robust recovery of certain input signals. The RIP condition is satisfied if the number of measurements M $\ge$ c*K*log(N/K) in which c is a constant ^[29]. To exactly recover x, the RIP condition is ensured. This is equivalent to the problem that Φ is incoherent with Ψ ^[30]. When this condition is satisfied, s can be also dependably reconstructed from y, using sparse approximation approaches, such as l₁-BP ^[1] or OMP ^[16]. In other work, compressed sensing for MRI imaging has been successfully applied by designing a random measurement matrix for data acquisition ^[31,32].

2.3.   Compressed sensing in MRI

For MR images, a high degree of the sparsity is needed because it suggests that data substance can be represented by a little information. One of these cases is MR angiography, where foundation tissue is ignored but the vessels are shown. It means that in the image domain, this kind of MR image is sparse. For other kinds of MR images, they are not sparse in the image domain but the transform domains ^[33]. There are several transform techniques such as the discrete wavelet transform (DWT), the discrete cosine transform (DCT), the fast Fourier transform (FFT), and finite difference operations, which can be used for representing the sparsity of these images. In this paper, we focus on MR images that are sparse in the wavelet and frequency domains. If the MR signal in the k-space is acquired with a small number of samples, the acquisition time will be reduced. Figure 3 illustrates realizable sampling designs that have been used in the previous studies. Both Figure 3a, b are classified as Cartesian sampled k-space but in Figure 3b, the number of samples is reduced twice (and the under-sampling ratio r is 0.5), which are different from the non-Cartesian method of radial under-sampling technique as shown in Figure 3c.

In this study, we only focus on the 2D Cartesian sampling. It has been found that the regular Cartesian under-sampling will cause the artifacts to manifest as coherent copies of the image structure as appeared in Figure 4a. The incoherent criterion is needed for the successful image reconstruction using compressive sensing. Thus, an outstanding candidate is random Cartesian k-space under-sampling, which satisfies the incoherent criterion. In this scheme, some portions of the phase-encoding steps are arbitrarily skipped, coming about in incoherent artifacts (Figure 4b).

The low-frequency components of Fourier basis functions, which locate at the origin of k-space, are highly correlated with the basic functions of most sparsifying transformations. Consequently, by collecting encoded information that locates around the origin of k-space, we are able to improve the performance of the MR image reconstruction. Because under-sampling the frequency-encoding (k_x) dimension does not influence to the acquisition time, hence, only the phase-encoding (k_y) dimension can be exploited in 2D Cartesian imaging.

However, if a small value of under-sampling ratio r is applied, the quality of the reconstructed image will be low (see section 4) because the random Cartesian k-space under-sampling scheme will be lost the important portion of the phase-encoding steps around the center of k-space even that the center of the scheme is at the center of k-space.

3.   Proposed scheme for MRI acquisition

Random under-sampling is used to break the regularity of Cartesian k-space. As shown in Figure 1a, the signal intensity reduces from the focal point of k-space to the periphery. Consequently, it is supposed that random schemes with their sampling distributions should be dense at the center and sparse at the periphery. In this paper, we concerned to two under-sampling approaches: (a) the traditional random under-sampling approach: It acquired Gaussian random measurements along the phase-encoding (k_y) in the k-space (Algorithm 1); and (b) Our proposed hybrid under-sampling approach (Algorithm 2). The proposed hybrid under-sampling approach is implemented as follows: For a certain value of the under-sampling ratio r (0 < r < 1), we divide the number of measurements in the (k_y) dimension into two portions: 70% of the measurements is for random under-sampling and 30% is for definite under-sampling taken near the origin of the k-space (see Algorithm 2). Because the amount of encoded information is concentrated at the origin of the k-space, the proposed method suggests that the amount of useful information will be collected more and therefore the MRI image recovery will be more accurate. Thus, instead of fully random under-sampling as described in the Algorithm 1, the proposed under-sampling described in Algorithm 2 takes the major part of measurements with the random pattern while the minor part of measurements is taken with a definite pattern without coinciding with any taken data from the random pattern (a conditional loop is performed by steps 4.1 to 4.4 in Algorithms 2). The advantage of our proposed algorithm is that encoded information near the k-space is always taken by the definite pattern.

The traditional random under-sampling approach and the proposed hybrid one are outlined in Algorithm 1 and Algorithm 2, respectively.

Image recovery of MRI can be implemented by using some iterative algorithms such as gradient descent, conjugate gradient, etc. ^[34,35]. These algorithms begin with an aliased MR image and continuously update to remove this artifact. In these algorithms, the regularization parameter λ is used to trade-off between information consistency (ℓ₂-norm) and the advancement of sparsity (ℓ₁-norm). In this paper, we choose a Nonlinear Conjugate Gradient (NCG) ^[14,17]. Assuming that m is the interest object, the recovered object $\hat m$ is achieved by solving this issue:

where F_u is the Fourier operator, y is the obtained measurements, and Ψ is the sparsifying transform operator.

To evaluate the performance of proposed methods, the normalized error is used to compare the error between the recovered object and the initial object. Assuming that m is an N × M initial object and $\hat m$ is the recovered object. The normalized error $\epsilon$ can be expressed as follows:

Another performance index, the universal image quality index (Q), is also used for evaluating the proposed under-sampling method as introduced by Wang and Bovik ^[36]. This index represents the distortion based on three different components: Loss of correlation, luminance distortion, and contrast distortion. The Q index is defined as:

where $\bar x$ and $\bar y$ are the mean of the original image and the reconstructed one, respectively; $\sigma _x^2$ and $\sigma _y^2$ are the variances of $x$ and $y$ ; and ${\sigma _{xy}}$ is the covariance between $x$ and $y$ . The Q index varies between −1 and 1 (Q index reaches to 1 if two images are identical).

4.   Numerical simulation results and discussions

4.1.   Numerical simulation results

To illustrate the advantage of the proposed method, the compression ratio r of 0.15 is first selected for evaluating the ε from reconstructed images of both two methods. The data source used in the numerical simulation is original brain MR slice with a 128 × 128 image size as shown in Figure 5.

Figures 6 and 7 present the binary masks that illustrate the random under-sampling implemented in the k-space based on the power-law using the traditional and proposed approaches. In the traditional approach, with the compression ratio r of 15%, the total number of sampled horizontal lines in the binary mask is rounded to 21 (horizontal bright lines as seen in Figure 6). Differently, in the proposed approach, the compression ratio r of 15% is divided into two parts: 11.5% for the traditional compression ratio r₁ and 4.5% for the definite compression ratio r₂. Thus, the numbers of traditional and definite sampled horizontal lines in the binary mask are rounded to 16 and 5, respectively (Figure 7). It means that, in the proposed approach, there are four sampled horizontal lines which are always taken around the origin of k-space for this compression ratio.

The histograms of the traditional random under-sampling and the hybrid approaches are shown in Figure 8. They both have the Gaussian shape that agrees with the literature. It can be seen that our proposed approach is more focused on the central area of the k-space than the previous one.

Figure 9 presented the reconstructed brain MR slice images using the traditional and the proposed approaches with under-sampling ratios of 0.15, 0.25, and 0.35. It can be seen that the reconstruction quality of the hybrid under-sampling approach is noticeably better than the previous one with the under-sampling ratios of 0.15 and 0.25. However, at the under-sampling ratios of 0.35, the reconstruction quality between two methods is more or less the same.

Figure 10 presents the averaged values of the ε calculated by 100 simulations of the traditional method and the newly proposed method for different compression ratios (r) from 0.05 to 0.5. It can be seen that the ε calculated by the new method is significantly lower than the value calculated by the traditional method, especially at the compression ratio of 0.15. However, there is almost no difference in the ε value between the proposed and traditional methods for r from 0.35 to 0.5. Figure 11 shows the results of the Q index comparison. We can see that, for compression ratios that are less than 0.35, the image reconstructed by the proposed approach offers a better value of Q than that reconstructed from traditional one.

To evaluate the time complexity of the proposed under-sampling method, the processing times of both under-sampling methods are estimated for comparison. The averaged processing time of each under-sampling approach is shown in Table 1, which indicates that the processing time of the proposed method is only marginally higher than the processing time of the conventional random method. Thus, our proposed under-sampling method can be considered as not increasing the time complexity of the algorithmic processes.

4.2.   Discussions

The compressed sensing reconstruction implements sparsity of the solution to suppress the incoherent aliasing artifacts and maximizes data consistency between the solution and the available under-sampling data. As shown in Figure 10, the average error of the proposed method decreased by 2.19 times (about 55%) compared to the traditional method at the compression ratio of 0.15. The Q-indexes of both methods described in Figure 11 also confirm the improvement of the proposed under-sampling approach with a maximum Q-index value of about 77% enhancement. At the compression ratio of 0.25, the proposed method decreased by 1.66 times to the traditional method. It can be explained that the traditional method samples the k-space followed by the power law, which has a similar energy distribution as the energy distribution of the k-space. The energy distribution of the k-space is followed by the Gauss distribution, so more high-energy points are located around the center of the k-space. When r is smaller than 0.35, only a few lines of k-space are sampled. The proposed method will outweigh the traditional method because more lines that are near the center of k-space are taken. However, when r is higher than 0.35, the number of sampling lines of k-space for both the traditional and proposed methods are more or less the same, but the sampled lines of the traditional method will match better with the energy distribution of the k-space. To prove the efficiency of the proposed under-sampling algorithm, both under-sampling methods are repeatedly simulated with an original knee MR image of the same size. The simulation results (not shown in the main text) have indicated consistent findings. The proportion of the random under-sampling over the define sampling one for the proposed under-sampling method is also investigated to find out the optimal proportion of 70 and 30% for the random and definite sampling patterns respectively. Moreover, the estimation of algorithmic processing time of the proposed under-sampling method confirms that this proposed approach does not increase the computational complexity for processing as compared to the traditional random method (Table 1).

5.   Conclusions

For the first time, a simple and efficient under-sampling approach based simulation is suggested in this paper for quality improvements of MRI image reconstruction using compressive sensing. Based on the fact that more information is concentrated at the origin of the k-space, instead of using the traditional random under-sampling approach that samples the phase-encoding (k_y) in the k-space followed by a power-law, we propose to use a large portion of the phase-encoding (k_y) in the k-space in the traditional way of sampling and strengthen the rest in sampling around the origin of k-space. Therefore, the more amount of useful information is permanently collected for reconstructing and imaging. The numerical simulation consequences have exhibited the efficiency of the appointed scheme. The obtained result of this study also indicates that instead of requiring the highly computational complexity of reconstruction algorithms, our proposed under-sampling method still has its own merits, which are able to be applied for reducing the image acquisition time of MRI systems.

Conflict of interest

The authors declares that they has no conflict of interest.