Next Article in Journal
A Scheme to Smooth Aggregated Traffic from Sensors with Periodic Reports
Next Article in Special Issue
On the Prediction of Flickr Image Popularity by Analyzing Heterogeneous Social Sensory Data
Previous Article in Journal
SVM-Based Spectral Analysis for Heart Rate from Multi-Channel WPPG Sensor Signals
Previous Article in Special Issue
Towards Building a Computer Aided Education System for Special Students Using Wearable Sensor Technologies
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Reconstruction of Undersampled Big Dynamic MRI Data Using Non-Convex Low-Rank and Sparsity Constraints

1
Department of Imaging and Interventional Radiology, The Chinese University of Hong Kong, Shatin 999077, N.T., Hong Kong, China
2
Department of Medicine and Therapeutics, The Chinese University of Hong Kong, Shatin 999077, N.T., Hong Kong, China
3
Chow Yuk Ho Technology Center for Innovative Medicine, The Chinese University of Hong Kong, Shatin 999077, N.T., Hong Kong, China
4
Department of Mathematics and Computer Science, Northeastern State University, Tahlequah, OK 74464, USA
5
Research Center for Medical Image Computing, The Chinese University of Hong Kong, Shatin 999077, N.T., Hong Kong, China
6
Shenzhen Research Institute, The Chinese University of Hong Kong, Shenzhen 518057, China
*
Authors to whom correspondence should be addressed.
Sensors 2017, 17(3), 509; https://doi.org/10.3390/s17030509
Submission received: 10 November 2016 / Revised: 16 February 2017 / Accepted: 20 February 2017 / Published: 3 March 2017
(This article belongs to the Special Issue Multisensory Big Data Analytics for Enhanced Living Environments)

Abstract

:
Dynamic magnetic resonance imaging (MRI) has been extensively utilized for enhancing medical living environment visualization, however, in clinical practice it often suffers from long data acquisition times. Dynamic imaging essentially reconstructs the visual image from raw (k,t)-space measurements, commonly referred to as big data. The purpose of this work is to accelerate big medical data acquisition in dynamic MRI by developing a non-convex minimization framework. In particular, to overcome the inherent speed limitation, both non-convex low-rank and sparsity constraints were combined to accelerate the dynamic imaging. However, the non-convex constraints make the dynamic reconstruction problem difficult to directly solve through the commonly-used numerical methods. To guarantee solution efficiency and stability, a numerical algorithm based on Alternating Direction Method of Multipliers (ADMM) is proposed to solve the resulting non-convex optimization problem. ADMM decomposes the original complex optimization problem into several simple sub-problems. Each sub-problem has a closed-form solution or could be efficiently solved using existing numerical methods. It has been proven that the quality of images reconstructed from fewer measurements can be significantly improved using non-convex minimization. Numerous experiments have been conducted on two in vivo cardiac datasets to compare the proposed method with several state-of-the-art imaging methods. Experimental results illustrated that the proposed method could guarantee the superior imaging performance in terms of quantitative and visual image quality assessments.

1. Introduction

Magnetic resonance imaging (MRI) has been established as an advanced non-invasive diagnostic imaging technique for visualizing the structure and functionality of the human body [1]. Many efforts have been made to dramatically improve the imaging speed and quality. In particular, the imaging of dynamic medical information plays an important role in enhancing medical living environment visualization [2,3,4,5,6,7]. Dynamic imaging has been becoming more and more important for several clinical applications, such as perfusion [8], cardiac [9] and functional imaging [10]. A large number of advanced image processing techniques (e.g., segmentation [11,12,13], classification [14,15,16], denoising [17,18] and deconvolution [19,20], etc.) have been developed based on MRI data to assist medical diagnosis.
MRI often suffers from long data acquisition times in clinical practice. Essentially, dynamic MRI, which reconstructs the visual image from raw (k,t)-space measurements, is commonly referred to as big data computing. The inherently slow acquisition speed often makes it difficult for dynamic MRI to capture high spatio-temporal resolutions, high signal-to-noise ratio, and proper volume coverage. The low spatio-temporal resolutions easily lead to erroneous diagnostic information and inaccurate estimation of tracer kinetic parameters, especially for dynamic susceptibility contrast MRI (DSC-MRI) [21]. The current undersampled reconstruction methods are of importance for the acceleration of imaging times in dynamic MRI. The spatio-temporal resolutions could be improved accordingly to capture the rapid tracer kinetics and enhance the diagnostic accuracy [22]. In particular, the accurate estimation of tracer kinetics plays a crucial role in quantifying the cerebral blood flow (CBF), cerebral blood volume (CBV) and mean transit time (MTT) to assess brain perfusion [23]. The high spatio-temporal resolutions could provide more medical information which is able to assist physicians in making an accurate diagnosis. Thus, there is a great potential to develop advanced imaging techniques to accelerate dynamic MRI in clinical practice.
Compressed sensing (CS) [24,25,26,27,28,29], which exploits the fact that an image can be sparsely represented in a certain transform domain, has been successfully used in dynamic MRI, such as k-t SPARSE [30], k-t FOCUSS [31,32] and k-t SPARSE-SENSE [33]. Note that the dynamic MRI sequence is both spatially and temporally correlated. The reconstructed results yielded by above CS-based methods may suffer from low spatial or temporal resolution. The purpose of this study is to accelerate dynamic MRI with high spatial and temporal resolutions. Recently, many efforts have been made to exploit the low-rank property of dynamic MRI sequence, instead of its only sparsity. Inherently, low-rank matrix completion can be considered as an effective extension of CS, which recoveries the missing or corrupted entries of a matrix under low-rank and incoherent conditions [34,35,36,37]. For instance, Lingala et al. [8,38] and Zhao et al. [39] have investigated the low-rank and sparse properties of the Casorati matrix (i.e., spatio-temporal MRI signal) and combined them together to further accelerate the dynamic MRI. The low-rank plus sparse decomposition model, referred to as robust principal component analysis (RPCA) [40,41], has also attracted increasing attention recently. Many results in the literature showed that it is possible to recover both low-rank and sparse components from a few incoherent observations under some assumptions [40,42,43].
Motivated by the theory of RPCA, Trémoulhéac et al. [44] proposed to reconstruct the dynamic MRI as the sum of low-rank plus sparse components. In this case, the dynamic MRI reconstruction problem was formulated as a least-squares optimization problem which was regularized by convex low-rank and sparsity constraints. Experimental results have demonstrated the effectiveness under different imaging conditions. The combination of low-rank with sparsity constraints can significantly improve the dynamic imaging result compared with the reconstruction models regularized by low-rank or sparsity constraint alone. It is generally thought that the superior performance of the combination version benefits from the full use of MRI data coherence and redundancy. Current studies [38,45,46,47] have shown that non-convex minimization could generate sparser solutions and provide better reconstruction in practice. It inherently means that the potential solutions can be sparsely represented in certain transform domains. Non-convex minimization has gained increasing attention in the fields of numerical optimization, image processing and computer version [48,49,50,51]. To make non-convex minimization more practical, a large number of numerical methods have been presented to solve the corresponding optimization problem. Inspired by success of non-convex optimization, there is a great potential to use the non-convex Schatten-p norm ( 0 < p < 1 ) as the low-rank constraint. In addition, to further enhance the sparse component recovery, non-convex Lq quasi-norm ( 0 < q < 1 ) can be effectively adopted to yield better result compared with commonly-used convex L1-norm [52,53,54,55,56]. In this paper, non-convex low-rank and sparsity constraints will be incorporated into our dynamic MRI reconstruction framework to improve imaging speed and quality.
In this work, we propose to formulate the undersampled dynamic MRI reconstruction as a least-squares optimization problem regularized by non-convex low-rank and sparsity constraints. The quality of the reconstructed images can be significantly improved by taking advantages of these constraints. However, due to the non-convex and non-smooth natures of the constraints, it is difficult to directly solve the resulting optimization problem through commonly-used numerical methods. To guarantee solution stability and efficiency, a numerical optimization algorithm based on Alternating Direction Method of Multipliers (ADMM) is proposed to solve the resulting optimization problem. In particular, the original non-convex non-smooth optimization problem will be decomposed into several subproblems by introducing several intermediate variables. These subproblems have simple closed-form solutions or could be efficiently handled using existing numerical methods. Thus, the main contributions of this paper mainly rely on the non-convex dynamic MRI reconstruction model and its numerical optimization algorithm. Experimental results on two in vivo cardiac MRI datasets have verified the effectiveness of the proposed method in terms of both quantitative and qualitative image quality evaluations.
The remainder of this paper is organized into several sections. The following section briefly introduces the basic concepts of dynamic MRI and robust principal component analysis. In Section 3, we tend to develop a non-convex minimization framework for accelerating dynamic MRI data acquisition using both non-convex low-rank and sparsity constraints. The resulting optimization algorithm is efficiently solved using an ADMM-based numerical algorithm. Numerous experiments on two in vivo cardiac datasets are performed in Section 4. Finally, we conclude this work in Section 5 by summarizing our contributions.

2. Background

2.1. Dynamic MRI from Partial Measurements

We denote the dynamic MRI to be reconstructed as a spatio-temporal signal I ( x , t ) , where x is the spatial coordinate and t denotes time. The imaging equation in dynamic MRI is defined as follows:
S ( k , t i ) = I ( x , t i ) exp ( j 2 π ( k · x ) ) d x + N ( k , t i )
for i = 1 , 2 , , T , where S ( k , t ) denotes the measured ( k , t ) -space signal, and N ( k , t ) is assumed to be additive complex-valued Gaussian noise. For the sake of simplicity, we consider a discrete image model in this dissertation. Given N t MR images of dimension N x × N y . The spatio-temporal signal I ( x , t ) can be rearranged into a matrix form, i.e.,
F = [ I ( 1 , t 1 ) I ( 1 , t N t ) I ( N x N y , t 1 ) I ( N x N y , t N t ) ] N x N y × N t
The Casorati matrix F , whose first and second directions respectively represent the spatial and temporal dimensions, is approximately low-rank due to the strong correlation existed between dynamic MR images [38]. Thus, a finite-dimensional spatio-temporal MRI model equivalent to Equation (1) can be written as:
d = 𝓐 ( F ) + n
where d P denotes the vector of the stacked ( k , t ) -space measurements, 𝓐 :   N x N y × N t P is the undersampled Fourier transform operator with P N x N y × N t , n P is the Gaussian noise vector. Recovering the matrix F from a limited number of measurements d is a typical ill-conditional problem. To cope with the ill-conditional nature, there has been considerable interest in exploiting low-rank and sparisty of the dynamic images to enhance the reconstruction accuracy [35]. Thus, there is a great potential to combine low-rank with sparsity constraints to improve dynamic MRI from partial measurements.

2.2. Robust Principal Component Analysis

Recent work has shown that robust principal component analysis (RPCA) [40,41] is capable of decomposing a Casorati matrix X into a low-rank component L and a sparse component S by constraining the rank of L and the sparsity of S simultaneously. It is worth noting that RPCA can be regarded as a robust extended version of traditional PCA [57]. PCA is performed based on the basic assumption of additive Gaussian noise and uses the sum of squared differences as the loss function. Theoretically, it works well as long as the value of noise is small enough. However, if the data are perturbed by high-level noise, it will be impossible to generate satisfactory reconstruction results since the traditional PCA could be easily corrupted by these gross errors [58]. It is well known that raw (k,t)-space data are often obtained from the MRI machines under complex imaging conditions. The raw data may be significantly corrupted by the outliers and large noise in clinical practice. If we directly adopt PCA to reconstruct the MRI images, it will be difficult to guarantee high-quality reconstruction due to the inherent limitation of PCA. In contrast, RPCA can perform well when the observed data are corrupted by severe perturbations because of its robust properties. The recent theoretical results indicate that the low-rank matrix can be exactly recovered by considering the following constrained convex minimization problem:
min L , S   L * + ρ S 1 such   that   X = L + S
where ρ is a weight parameter, the nuclear norm L * is defined as L * = i = 1 r σ i with σ 1 ,   σ 2 , , σ r being the singular values of L and r being the rank of L. In the literature, the classical ADMM [59] can be exploited to efficiently solve the RPCA Problem (3) based on the following augmented Lagrangian function:
𝓛 A =   L * + ρ S 1 + α 2 X ( L + S ) + Z α 2 2
where Z denotes the Lagrange multiplier and α is a positive penalty parameter. ADMM makes full use of the separable structure of Equation (4), and decomposes Equation (4) over primal variables L and S. The iterative procedure is detailedly summarized in Algorithm 1.
Algorithm 1. RPCA
1 Input: Casorati matrix X , decomposition parameter ρ > 0
2 Initialize: S 0 = Z 0 = 0 , k = 0
3 while stopping criterion is not satisfied do
4   L k + 1 = T α 1 ( X S k + α 1 Z k )
5   S k + 1 = S ρ / α ( X L k + 1 + α 1 Z k )
6   Z k + 1 = Z k + α ( X L k + 1 S k + 1 )
7 end while
8 Output: ( L * , S * )
In particular, ADMM minimizes 𝓛 A in Equation (4) over L and S separately, and then updates the Lagrange multiplier Z . The corresponding L-subproblem and S-subproblem are given by:
L k + 1 = min   L L * + α 2 ( X S k + Z k α ) L 2 2
S k + 1 = min   S ρ S 1 + α 2 ( X L k + 1 + Z k α ) S 2 2
where k represents the iteration number. The closed-form solutions to Equations (5) and (6) are respectively obtained using the singular value thresholding operator T τ ( Y ) = U S τ ( Ξ ) V H (where Y = U Ξ V H denotes the singular value decomposition) and shrinkage operator S τ ( Y ) = s i g n ( Y ) max ( | Y | τ , 0 ) , i.e.,
L k + 1 = T α 1 ( X S k + α 1 Z k )
S k + 1 = S ρ / α ( X L k + 1 + α 1 Z k )
In Equations (7) and (8), the definition of sign function s i g n ( · ) is given by:
s i g n ( s ) = { 1 ,   s > 0 0 ,   s = 0 1 ,   s < 0
Given the fixed values of L k + 1 and S k + 1 , the Lagrange multiplier Z is accordingly updated as Z k + 1 = Z k + α ( X L k + 1 S k + 1 ) . The parameter ρ can be seen as a trade-off between low-rank and sparse components. The theoretically suggested value introduced in [40] is set to ρ = 1 / max ( N x N y , N t ) , where N x N y and N t denote the number of pixels in each frame and the number of time-frames in Casorati matrix X , respectively.

3. k-t NCRPCA: Formulation

3.1. Joint Non-Convex Low-Rank and Sparsity Constraints

In this study, it is assumed that the Casorati matrix F can be decomposed as a low-rank component L and a sparse component S . Figure 1 shows the L + S decomposition of an axial cardiac MRI dataset, where L captures the corrected background between time frames and S captures the temporal or dynamic information.
In can be observed that the singular values of L tend to zero quickly as their index increases, and S can be effectively represented using a sparsifying transform. Based on the spatio-temporal MRI Model (2) and RPCA (3), Trémoulhéac et al. [44] proposed to reconstruct the dynamic MRI F as the sum of low-rank L plus sparse S components. The dynamic MRI reconstruction can be formulated as a least squares optimization problem regularized by convex low-rank and sparsity constraints, i.e.,
min L , S 1 2 𝓐 ( L + S ) d 2 2 + μ 1 L * + μ 2 F t ( S ) 1
where both μ 1 and μ 2 are positive regularization parameters, and F t represents the Fourier transform operator along the temporal direction. The assumption behind the use of this sparsifying transform is that the dynamic MRI in time exhibits strong correlation or periodicity. The unclear norm * is equivalent to the L1-norm 1 because the singular values are all nonnegative. Current studies in the literature [38,45,46,47] illustrated that non-convex minimization could yield sparser solutions and provide consistently better performance over convex minimization. Thus, there is a great potential to use non-convex Schatten-p norm ( 0 < p < 1 ) as the low-rank constraint. To improve the sparse component recovery, non-convex Lq quasi-norm ( 0 < q < 1 ) can yield better recovery results over commonly-used convex L1-norm [52,53,54,55]. Motivated by the success of non-convex Lq quasi-norm, we propose to replace the convex L1-norm F t ( S ) 1 by its non-convex version F t ( S ) q q in Equation (9) to further improve the quality of reconstructed image.
With the above notations, we propose to formulate undersampled dynamic MRI reconstruction as a least-squares optimization problem regularized by non-convex low-rank and sparsity constraints, i.e.,
min L , S 1 2 𝓐 ( L + S ) d 2 2 + μ 1 L p p + μ 2 F t ( S ) q q
where L p p = i = 1 min { N x N y , N t } σ i p ( L ) for some p ( 0 , 1 ) denotes the non-convex low-rank constraint in which σ i ( L ) is the ith largest singular value of L and min { N x N y , N t } is the rank of the matrix L . The non-convex sparsity constraint, defined as F t ( S ) q q = i = 1 N x N y N t | ( F t ( S ) ) i | q for some q ( 0 , 1 ) , is used to promote sparsity in Fourier transform domain. The positive regularization parameters μ 1 and μ 2 are adjusted to balance the trade-off between the data-fidelity term and the joint non-convex low-rank and sparsity constraints. The sparsifying transform F t is incoherent with the Fourier sampling operator 𝓐 , thus the proposed regularized reconstruction is well posed in practice.

3.2. Numerical Optimization Algorithm

The resulting optimization Problem (10) includes both non-convex low-rank and sparsity constraints. The introduced non-convex constraints make dynamic reconstruction problem difficult to solve. Thus, it is intractable to solve Problem (10) directly using current minimization schemes. To guarantee solution stability and efficiency, an ADMM-based optimization algorithm is proposed in this work. Two auxiliary variables P = L and Q = F t ( S ) are first introduced and the unconstrained minimization Problem (10) is then transformed into the following constrained version:
min P , Q , L , S 1 2 𝓐 ( L + S ) d 2 2 + μ 1 P p p + μ 2 Q q q subject   to   P = L ,   Q = F t ( S )
whose associated augmented Lagrangian function is given by:
𝓛 A = 1 2 𝓐 ( L + S ) d 2 2 + μ 1 P p p + μ 2 Q q q + α 1 2 L P + Z 1 α 1 2 2 + α 2 2 F t ( S ) Q + Z 2 α 2 2 2
where Z 1 and Z 2 denote the Lagrange multipliers, α 1 and α 2 represent positive penalty parameters that control the weights of penalty terms. As α 1 , α 2 , solution of the above Problem (12) tends to that of Problem (11).
Since the variables P , Q , L and S are coupled together in the augmented Lagrangian Function (12), it is computationally intractable to solve them simultaneously. ADMM minimizes 𝓛 A over P , Q , L and S separately leading to several subproblems which have closed-form solutions or could be efficiently solved using existing numerical methods. This decouples the individual updates of P , Q , L and S , therefore the original optimization task can be simplified as follows:
P k + 1 = min   P μ 1 P p p + α 1 2 P ( L k + Z 1 k / α 1 ) 2 2
Q k + 1 = min   Q μ 2 Q q q + α 2 2 Q ( F t ( S k ) + Z 2 k / α 2 ) 2 2
L k + 1 = 1 2 𝓐 ( L + S k ) d 2 2 + α 1 2 L ( P k + 1 Z 1 k / α 1 ) 2 2
S k + 1 = 1 2 𝓐 ( L k + 1 + S ) d 2 2 + α 2 2 F t ( S ) ( Q k + 1 Z 2 k / α 2 ) 2 2
The first P -subproblem (13) can be considered as a similar form of the standard unclear norm minimization problem. To achieve an efficient solution, the iterative singular value thresholding scheme [60,61], originally proposed for unclear norm minimization, can be extended to solve (13) which has a non-convex low-rank constraint. Let P ˜ k = L k + Z 1 k / α 1 , the solution P k + 1 is given by:
P k + 1 = i = 1 min { N x N y , N t } max ( σ i μ 1 σ i p 1 α 1 , 0 ) u i v i T
where the superscript T denotes the transpose (conjugate transpose) operator for real (complex) matrices or vectors. In Equation (17), u i , v i and σ i respectively represent the singular vectors and values of P ˜ k . Note that, if p = 1 , the expression in Equation (17) can be regarded as equivalent to a standard shrinkage algorithm for nuclear norm optimization problem.
It is difficult to directly solve the second subproblem (14) owing to the non-convex sparsity constraint. Motivated by the successful applications of soft thresholding [62] and iterative shrinkage/thresholding [63], a generalized iterated shrinkage algorithm (GISA) algorithm [64] is introduced in this study to efficiently solve the non-convex minimization problem. In particular, the Q -subproblem (14) can be decoupled into N x N y × N t independent and unconstrained subproblems. The standard version of these subproblems is given by:
y * = min   y μ ˜ | y | q + 1 2 ( y c ) 2
where μ ˜ is a positive parameter, y denotes the scalar variable that needs to be estimated, and c is a known scalar constant. As discussed in [64], the proposed GISA algorithm can be adopted to efficiently solve the non-convex Problem (18) with high accuracy. Therefore, the optimal solution of (18) can be obtained as follows:
y * = { 0 ,   if   | c | τ μ ˜ sign ( c ) S T q ( | c | ,   μ ˜ ) ,   if   | c | > τ μ ˜
where:
τ μ ˜ = ( 2 μ ˜ ( 1 q ) ) 1 / ( 2 q ) + μ ˜ q ( 2 μ ˜ ( 1 q ) ) ( q 1 ) / ( 2 q )
and S T q ( | c | ,   μ ˜ ) can be obtained by iteratively performing the following equation:
S T q ( | c | ,   μ ˜ ) | c | + μ ˜ q ( S T q ( | c | ,   μ ˜ ) ) q 1 = 0
Thus, solution Q k + 1   of the subproblem (14) is given by:
Q k + 1 = { 0 ,   if   | F t ( S k ) + Z 2 k / α 2 | τ μ 2 / α 2 sign ( F t ( S k ) + Z 2 k / α 2 ) S T q ( | F t ( S k ) + Z 2 k / α 2 | ,   μ 2 α 2 ) ,   if   | F t ( S k ) + Z 2 k / α 2 | > τ μ 2 / α 2
Both the third and fourth subproblems (15) and (16) are essentially quadratic and thus can be solved analytically as:
L k + 1 = ( 𝓐 𝓐 + α 1 I ) 1 ( 𝓐 d + α 1 P k + 1 Z 1 k 𝓐 𝓐 S k )
S k + 1 = ( 𝓐 𝓐 + α 2 I ) 1 ( 𝓐 d + F t ( α 2 Q k + 1 Z 2 k ) 𝓐 𝓐 L k + 1 )
For fixed values of P k + 1 , Q k + 1 , L k + 1 and S k + 1 , the Lagrange multipliers Z 1 and Z 2 are updated as follows:
Z 1 k + 1 = Z 1 k α 1 ( P k + 1   L k + 1 )
Z 2 k + 1 = Z 2 k α 2 ( Q k + 1   F t ( S k + 1 ) )
Based on these analytic solutions, the proposed image reconstruction procedure is detailedly summarized in Algorithm 2. The original non-convex optimization problem is decomposed into four subproblems. Each of these subproblems has a closed-form solution or could be efficiently solved using existing numerical method. Due to the non-convex and non-smooth constraints, it is difficult to yield the resulting theoretical convergence result. However, the quality of reconstructed image could be improved using the joint non-convex low-rank and sparsity constraints. In this study, the proposed method is referred to as k-t NCRPCA for dynamic MRI reconstruction from undersampled ( k , t ) -space measurements.
Algorithm 2. k-t NCRPCA
1 Input: Fourier transform 𝓐 ,   ( k , t ) -space data d and parameters ( μ 1 , μ 2 , α 1 , α 2 , p , q , τ α 2 , μ 2 ) .
2 Initialize: F 0 = L 0 = 𝓐 d , S 0 = Z 1 0 = Z 2 0 = 0 .
3 while stopping criterion is not satisfied do
4   P k + 1 = i = 1 min { N x N y , N t } max ( σ i μ 1 σ i p 1 α 1 , 0 ) u i v i
5   Q k + 1 = { 0 ,   if   | F t ( S k ) + Z 2 k / α 2 | τ μ 2 / α 2 sign ( F t ( S k ) + Z 2 k / α 2 ) S T q ( | F t ( S k ) + Z 2 k / α 2 | , μ 2 α 2   ) ,   if   | F t ( S k ) + Z 2 k / α 2 | > τ μ 2 / α 2
6   L k + 1 = ( 𝓐 𝓐 + α 1 I ) 1 ( 𝓐 d + α 1 P k + 1 Z 1 k 𝓐 𝓐 S k )
7   S k + 1 = ( 𝓐 𝓐 + α 2 I ) 1 ( 𝓐 d + F t ( α 2 Q k + 1 Z 2 k ) 𝓐 𝓐 L k + 1 )
8   Z 1 k + 1 = Z 1 k α 1 ( P k + 1   L k + 1 )
9   Z 2 k + 1 = Z 2 k α 2 ( Q k + 1   F t ( S k + 1 ) )
10   F k + 1 = L k + 1 + S k + 1
11   α 1 = 1.2 α 1 and α 2 = 1.2 α 2
12 end while
13 Output: ( F * , L * , S * )

4. Experimental Results and Discussion

Experimental results on in vivo axial and coronal cardiac MRI datasets were conducted in this study to demonstrate the superior performance of our proposed method in terms of both quantitative and visual quality evaluations.

4.1. Acquired Datasets

To evaluate the performance of dynamic MRI reconstruction, experiments were conducted on both in vivo axial and coronal cardiac datasets. The 2D cardiac cine imaging was performed in a healthy adult volunteer, with the approval from The Joint Chinese University of Hong Kong—New Territories East Cluster Clinical Research Ethics Committee (The Joint CUHK-NTEC CREC). The datasets were acquired from a clinical 3T MRI scanner (Achieva, Philips Medical Systems, Best, The Netherlands) with an eight-channel receiver coil. The relevant imaging parameters were as follows: TR/TE = 3.8 / 1.9 ms, flip angle = 45°, image matrix size = 128 × 128 , and temporal frames = 60. During cardiac imaging, the volunteer was instructed to hold the breath for as long as possible. As shown in Figure 2, the undersampled ( k , t ) -space measurements were obtained from 8 , 12 , 16 , 24 and 32 radial projections. To obtain a full Nyquist-sampled dataset, π 2 n projections ( 201 for n = 128 in our case) should be acquired theoretically. Therefore, the above radial projections respectively correspond to undersampling factors of ~ 25 , 16 , 12 , 8 and 6 . In addition, the radial sampling pattern shown in Figure 2 has uniformly spaced radial rays per frame and subsequent random rotations across time frames for maintaining sampling incoherence. All experiments were implemented in MATLAB® (The MathWorks Inc., Natick, MA, USA) using a computer with 3.1 GHz Intel Core i5-2500 PU, 4 GB RAM. For comparison, the high-quality results reconstructed from fully sampled ( k , t ) -space measurements were used as ground truth reference in our numerical experiments. The proposed method will be compared with three state-of-the-art reconstruction methods, i.e., k-t FOCUSS [31,32], k-t SLR [38] and k-t RPCA [44]. k-t FOCUSS is able to take full advantage of the sparse properties of dynamic MRI to enhance reconstruction quality. The last two reconstruction methods were proposed by simultaneously using the low-rank and sparsity constraints. The resulting optimization methods were effectively solved through the alternating direction methods. Note that k-t RPCA can be considered as a special case of our proposed method (10). The satisfactory reconstructed results were generated by these three methods with the corresponding parameters manually optimized.

4.2. Parameters Settings

The accuracy of reconstruction results will be evaluated using the signal to error ratio (SER) [38] and the structural similarity index (SSIM). SSIM is more consistent with the human visual perception. The SER in mathematical form is defined as follows:
SER =   10 log 10 F rec F truth 2 2 F truth 2 2   ( dB )
where F rec and F truth denote the recovered matrix and the ground truth fully-sampled noiseless matrix, respectively. Let f rec and f truth represent the time-frame magnitude images extracted from F rec and F truth , respectively. The definition of SSIM between f rec and f truth is given by:
SSIM = ( 2 μ f rec μ f truth + c 1 ) ( 2 σ f rec f truth + c 2 ) ( μ f rec 2 + μ f truth 2 + c 1 ) ( σ f rec 2 + σ f truth 2 + c 2 )
where μ f rec and μ f truth represent the local mean values, σ f rec and σ f truth denote the standard deviations, σ f rec f truth is the covariance value, and c 1 and c 2 are two predefined constants to avoid instability with weak denominator. Both SER and SSIM are able to provide quantitative indexes of reconstruction results, but they cannot be well correlated with perceptual quality. Thus, the specific reconstructed frames and the temporal profiles will also be shown to enhance visual comparisons.
The proposed method in Algorithm 2 mainly involves six parameters, i.e., μ 1 , μ 2 , α 1 , α 2 , p and q . Inspired by previous work [40], a suggested μ 2 is set to μ 2 = ρ μ 1 with ρ = 1 / max ( N x N y , N t ) , where N x N y and N t respectively denote the number of pixels in each frame and the number of time-frames in Casorati matrix F . The regularization parameter μ 1 is manually set as 2 × 10 2 for k-t NCRPCA. The convergence rate of k-t NCRPCA is highly dependent on the selection of α 1 and α 2 . It has been proven that lower values of α 1 and α 2 result in a much faster convergence rate. However, the proposed method would suffer from image quality degradation since constraints in Equation (4) are not satisfied. Higher values of α 1 and α 2 can guarantee that all constraints can be satisfied, but generally lead to slower convergence. To balance the trade-off between convergence rate and reconstruction accuracy, the proposed method in Algorithm 2 employs a continuous strategy where both parameters α 1 and α 2 are initialized to small values, and are gradually increased until the stopping criteria is satisfied.
The parameters p and q are respectively related to the non-convex low-rank and sparsity constraints. They play an important role in the improvement of reconstructed image quality. Exhaustive experiments in this section were performed to manually determine the optimal choices of p and q . Take the in vivo axial cardiac dataset as an example, the undersampled ( k , t ) -space measurements were generated using 12 radial projections. A series of searches within a predefined range of parameters was performed to select the optimal values. For the sake of simplicity, the predefined range [ 0.1 ,   0.2 , ,   1.0 ] was set for both p and q . As shown in Figure 3, the proposed k-t NCRPCA is more sensitive to p than q . It means that non-convex low-rank constraint proposed in method (10) can significantly affect the reconstructed image quality. According to the quantitative results, the resulting optimal parameters were set as p = 0.9 and q = 0.8 for axial cardiac dataset. These manually-selected parameters were also used for coronal cardiac imaging experiments. The reconstruction results under current parameter settings were consistently promising in our numerical experiments. Further study on automatically calculating the above-discussed parameters for k-t NCRPCA is our future work.
The standard for the stopping criteria in Algorithm 2 is that the relative change of F is sufficiently small, i.e.,
F k + 1 F k 2 / F k 2 ϵ   or   k > k max
where ϵ and k max , respectively, denote the predefined tolerance parameter and maximum number of iterations.
For all reconstruction experiments, ϵ = 10 4 and k max = 300 were respectively set in Equation (27). It is fair to compare the experimental results because the results of other competing methods are all generated using the authors’ codes with the parameters manually optimized.

4.3. Comparisons on In Vivo Axial Cardiac Dataset

In order to evaluate the reconstruction performance, the proposed k-t NCRPCA will be compared to other dynamic reconstruction techniques, i.e., zero-filled inverse Fourier transform (ZF-IDFT), k-t FOCUSS [31,32], k-t SLR [38] and k-t RPCA [44]. The good performance of k-t FOCUSS benefits from the sparse properties of dynamic MRI. k-t SLR took into account both sparsity and spectral priors to accelerate dynamic imaging. k-t RPCA formulated dynamic reconstruction as a least-squares optimization problem which was regularized by a low-rank plus sparse prior.
Figure 4 visually illustrate the time-frame magnitude images, local magnification views, temporal profiles and 1D profiles for the in vivo axial cardiac dataset with 8 radial projections (correspond to undersampling factor of ~ 25 ). It can be observed that ZF-IDFT generates the worst perceptual results. The loss of information details severely degrades the image quality. k-t FOCUSS and k-t SLR could reconstruct the main geometrical structures but seemed to overcome some fine details. As shown by the white ellipse in y-t temporal profiles, both k-t RPCA and k-t NCRPCA perform well in preserving fine image details to guarantee high image quality. For the sake of comparison, the reconstruction performance is further confirmed by the 1D profiles. It is easy to see that the intensity values of k-t NCRPCA are more structurally similar to the fully sampled reference image. In contrast, k-t FOCUSS, k-t SLR and k-t RPCA lead to accuracy loss in determining the temporal profiles due to reconstruction biases. The good performance of k-t NCRPCA over other methods can possibly be attributed to the fact that non-convex minimization benefits for improving image reconstruction from fewer measurements.
To further evaluate the reconstruction result, more 1D profiles for different number of radial projections are visually illustrated in Figure 5. Similar to our findings in Figure 4, we can find significant performance improvement of k-t NCRPCA in comparison to other competing methods. The quantitative results are given in Table 1. It can be concluded that our proposed k-t NCRPCA achieves the best reconstruction results in all cases.

4.4. Comparisons on In Vivo Coronal Cardiac Dataset

In Figure 6, we compared our k-t NCRPCA with ZF-IDFT, k-t FOCUSS, k-t SLR and k-t RPCA on an in vivo coronal cardiac dataset for 8 radial projections. It is impossible to satisfactorily reconstruct the main geometical structures using ZF-IDFT. As shown by the arrows in x-t and y-t temporal profiles, k-t FOCUSS, k-t SLR and k-t RPCA failed to achieve high quality temporal profiles due to over smoothing. The loss of useful visual information could lead to visual quality degradation. In contrast, k-t NCRPCA performs well in feature-preserving image reconstruction even for highly undersampled ( k , t ) -space measurements.
For 1D profiles, k-t NCRPCA is able to preserve more structural information. For the sake of visual comparison, we generated the error images resulting from the absolute differences between the fully sampled reference and reconstructed images shown in Figure 7. To achieve high quality reconstruction, the error images should include as little geometrical information as possible. As can be observed, the geometrical structures are noticeable in error images generated by ZF-IDFT. k-t FOCUSS, k-t SLR and k-t RPCA tend to oversmooth some image details. In contrast, our k-t NCRPCA generates error images with more random and smaller absolute differences, which result in benefit to image quality improvement. The advantage of k-t NCRPCA is further confirmed by the SER comparison in Table 2.
The main benefit of our proposed method is that it takes full advantage of non-convex low-rank and sparsity constraints. Current studies have illustrated that non-convex minimization could provide consistently better performance over convex minimization [32,39,40]. Thus, our method is able to effectively preserve fine details while suppressing undesirable artifacts for dynamic MRI reconstruction from undersampled ( k , t ) -space measurements.

4.5. Algorithm Convergence and Robustness

The convergence of the proposed k-t NCRPCA was investigated on in vivo axial and coronal cardiac datasets. The undersampled ( k , t ) -space measurements for these two datasets were generated from 8 , 16 and 32 radial projections. It is difficult to establish the theoretical convergence result because of the non-convex and non-smooth objective Function (10). For the sake of simplicity, the convergence property of the proposed algorithm was analyzed empirically. As shown in Figure 8, the objective quality metric SER is visually displayed as a function of iteration number. It is observed that with the increase of iteration number, SER values increase quickly at the first few iterations and then become stable. These observations illustrate that the convergence property of the proposed method, summarized in Algorithm 2, can be guaranteed for undersampled dynamic MRI reconstruction. Furthermore, it is obvious that the proposed method is convergent even for robust reconstruction from highly undersampled ( k , t ) -space measurements (i.e., 8 radial projection in experiments). The experiments on in vivo axial and coronal cardiac datasets fully illustrate the robustness of the proposed method.

5. Conclusions

In this study, a new constrained imaging method (termed k-t NCRPCA) was proposed to accelerate dynamic MRI. It effectively integrated both non-convex low-rank and sparsity constraints into a unified mathematical framework. The resulting non-convex and non-smooth optimization problem was effectively solved using an ADMM-based optimization algorithm. Numerous experiments have been conducted on two in vivo cardiac datasets to compare k-t NCRPCA with other state-of-the-art reconstruction methods. The experimental results have demonstrated the superior performance of k-t NCRPCA in terms of both quantitative and qualitative image quality evaluations. Therefore, there is a strong incentive to extend k-t NCRPCA to effectively accelerate dynamic MRI in clinical practice.

Acknowledgments

This work was partially supported by grants from the Research Grants Council of HKSAR (Project Nos. CUHK 416712 and 14113214), a grant from the Science, Technology and Innovation Commission of Shenzhen Municipality (Project No. CXZZ20140606164105361), grants from the Innovation and Technology Commission (Project Nos. GHP/028/14SZ and ITS/293/14FP), and grants from CUHK Technology and Business Development Fund (Project Nos. TBF16MED002 and TBF16MED004).

Author Contributions

The work presented in this paper corresponds to a collaborative development by all authors. R.W.L., L.S. and D.W. conceived and designed the experiments; R.W.L., L.S., S.C.H.Y., N.X. and D.W. performed the numerical experiments; R.W.L., S.C.H.Y., N.X. and D.W. analyzed the experimental results; R.W.L., N.X. and D.W. wrote the manuscript and improved this manuscript's English language and style.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Edelman, R.R.; Warach, S. Magnetic resonance imaging. N. Engl. J. Med. 1993, 328, 785–791. [Google Scholar] [CrossRef] [PubMed]
  2. Xiong, N.; Vasilakos, A.V.; Yang, L.T.; Song, L.; Pan, Y.; Kannan, R.; Li, Y. Comparative analysis of quality of service and memory usage for adaptive failure detectors in healthcare systems. IEEE J. Sel. Areas Commun. 2009, 27, 495–509. [Google Scholar] [CrossRef]
  3. O'connor, J.P.; Jackson, A.; Parker, G.J.; Roberts, C.; Jayson, G.C. Dynamic contrast-enhanced MRI in clinical trials of antivascular therapies. Nat. Rev. Clin. Oncol. 2012, 9, 167–177. [Google Scholar] [CrossRef] [PubMed]
  4. Xiong, N.; Jia, X.; Yang, L.T.; Vasilakos, A.V.; Li, Y.; Pan, Y. A distributed efficient flow control scheme for multirate multicast networks. IEEE Trans. Parallel Distrib. Syst. 2010, 21, 1254–1266. [Google Scholar] [CrossRef]
  5. Harada, T.; Tsuji, Y.; Mikami, Y.; Hatta, Y.; Sakamoto, A.; Ikeda, T.; Kubo, T. The clinical usefulness of preoperative dynamic MRI to select decompression levels for cervical spondylotic myelopathy. Magn. Reson. Imaging 2010, 28, 820–825. [Google Scholar] [CrossRef] [PubMed]
  6. Leach, M.O.; Morgan, B.; Tofts, P.S.; Buckley, D.L.; Huang, W.; Horsfield, M.A.; Whitcher, B. Imaging vascular function for early stage clinical trials using dynamic contrast-enhanced magnetic resonance imaging. Eur. Radiol. 2012, 22, 1451–1464. [Google Scholar] [CrossRef] [PubMed]
  7. Xia, Z.; Xiong, N.N.; Vasilakos, A.V.; Sun, X. EPCBIR: An efficient and privacy-preserving content-based image retrieval scheme in cloud computing. Inf. Sci. 2016, 387, 195–204. [Google Scholar] [CrossRef]
  8. Lingala, S.G.; DiBella, E.; Adluru, G.; McGann, C.; Jacob, M. Accelerating free breathing myocardial perfusion MRI using multi coil radial k-t SLR. Phys. Med. Biol. 2013, 58, 7309–7327. [Google Scholar] [CrossRef] [PubMed]
  9. Denney, T.R.; Prince, J.L. Reconstruction of 3-D left ventricular motion from planar tagged cardiac MR images: An estimation theoretic approach. IEEE Trans. Med. Imaging 1995, 14, 625–635. [Google Scholar] [CrossRef] [PubMed]
  10. Thesen, S.; Heid, O.; Mueller, E.; Schad, L.R. Prospective acquisition correction for head motion with image-based tracking for real-time fMRI. Magn. Reson. Med. 2000, 44, 457–465. [Google Scholar] [CrossRef]
  11. Mahapatra, D.; Buhmann, J.M. Prostate MRI segmentation using learned semantic knowledge and graph cuts. IEEE Trans. Biomed. Eng. 2014, 61, 756–764. [Google Scholar] [CrossRef] [PubMed]
  12. Xia, Z.; Wang, X.; Sun, X.; Liu, Q.; Xiong, N. Steganalysis of LSB matching using differences between nonadjacent pixels. Multimed. Tools Appl. 2016, 75, 1947–1962. [Google Scholar] [CrossRef]
  13. Chen, Y.; Zhang, H.; Zheng, Y.; Jeon, B.; Wu, Q.J. An improved anisotropic hierarchical fuzzy c-means method based on multivariate student t-distribution for brain MRI segmentation. Pattern Recognit. 2016, 60, 778–792. [Google Scholar] [CrossRef]
  14. Zhang, Y.; Dong, Z.; Wu, L.; Wang, S. A hybrid method for MRI brain image classification. Expert Syst. Appl. 2011, 38, 10049–10053. [Google Scholar] [CrossRef]
  15. Wen, X.; Shao, L.; Xue, Y.; Fang, W. A rapid learning algorithm for vehicle classification. Inf. Sci. 2015, 295, 395–406. [Google Scholar] [CrossRef]
  16. Li, Z.; Lai, Z.; Xu, Y.; Yang, J.; Zhang, D. A locality-constrained and label embedding dictionary learning algorithm for image classification. IEEE Trans. Neural Netw. Learn. Syst. 2017, 28, 278–293. [Google Scholar] [CrossRef] [PubMed]
  17. Liu, R.W.; Shi, L.; Yu, S.C.; Wang, D. A two-step optimization approach for nonlocal total variation-based Rician noise reduction in magnetic resonance images. Med. Phys. 2015, 42, 5167–5187. [Google Scholar] [CrossRef] [PubMed]
  18. Baselice, F.; Ferraioli, G.; Pascazio, V.; Sorriso, A. Bayesian MRI denoising in complex domain. Magn. Reson. Imaging 2017, 38, 112–122. [Google Scholar] [CrossRef] [PubMed]
  19. Pizzolato, M.; Boutelier, T.; Deriche, R. Perfusion deconvolution in DSC-MRI with dispersion-compliant bases. Med. Image Anal. 2017, 36, 197–215. [Google Scholar] [CrossRef] [PubMed]
  20. Koh, T.S.; Cheong, D.L.H.; Hou, Z. Issues of discontinuity in the impulse residue function for deconvolution analysis of dynamic contrast-enhanced MRI data. Magn. Reson. Med. 2011, 66, 886–892. [Google Scholar] [CrossRef] [PubMed]
  21. Han, S.; Paulsen, J.L.; Zhu, G.; Song, Y.; Chun, S.; Cho, G.; Ackerstaff, E.; Koutcher, J.A.; Cho, H. Temporal/spatial resolution improvement of in vivo DCE-MRI with compressed sensing-optimized FLASH. Magn. Reson. Imaging 2012, 30, 741–752. [Google Scholar] [CrossRef] [PubMed]
  22. Pain, F.; Besret, L.; Vaufrey, F.; Grégoire, M.C.; Pinot, L.; Gervais, P.; Ploux, L.; Bloch, G.; Mastrippolito, R.; Lanièce, P.; et al. In vivo quantification of localized neuronal activation and inhibition in the rat brain using a dedicated high temporal-resolution β+-sensitive microprobe. Proc. Natl. Acad. Sci. USA 2002, 99, 10807–10812. [Google Scholar] [CrossRef] [PubMed]
  23. Boschetto, D.; Di Prima, P.; Castellaro, M.; Bertoldo, A.; Grisan, E. Baseline constrained reconstruction of DSC-MRI tracer kinetics from sparse fourier data. In Proceedings of the IEEE International Symposium on Biomedical Imaging, Beijing, China, 29 April–2 May 2014; pp. 321–324.
  24. Candes, E.J.; Romberg, J.; Tao, T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans. Inf. Theory 2006, 52, 489–509. [Google Scholar] [CrossRef]
  25. Donoho, D.L. Compressed sensing. IEEE Trans. Inf. Theory 2006, 52, 1289–1306. [Google Scholar] [CrossRef]
  26. Lustig, M.; Donoho, D.; Pauly, J.M. Sparse MRI: The application of compressed sensing for rapid MR imaging. Magn. Reson. Med. 2007, 58, 1182–1195. [Google Scholar] [CrossRef] [PubMed]
  27. Lustig, M.; Donoho, D.L.; Santos, J.M.; Pauly, J.M. Compressed sensing MRI. IEEE Signal Process. Mag. 2008, 25, 72–82. [Google Scholar] [CrossRef]
  28. Zhang, Y.; Peterson, B.S.; Ji, G.; Dong, Z. Energy preserved sampling for compressed sensing MRI. Comput. Math. Method Med. 2014, 2014, 546814. [Google Scholar] [CrossRef] [PubMed]
  29. Zhang, Y.; Dong, Z.; Phillips, P.; Wang, S.; Ji, G.; Yang, J. Exponential wavelet iterative shrinkage thresholding algorithm for compressed sensing magnetic resonance imaging. Inf. Sci. 2015, 322, 115–132. [Google Scholar] [CrossRef]
  30. Cai, N.; Wang, S.; Zhu, S.; Liang, D. Accelerating dynamic cardiac MR imaging using structured sparse representation. Comput. Math. Methods Med. 2013, 2013, 160139. [Google Scholar] [CrossRef] [PubMed]
  31. Jung, H.; Sung, K.; Nayak, K.S.; Kim, E.Y.; Ye, J.C. k-t FOCUSS: A general compressed sensing framework for high resolution dynamic MRI. Magn. Reson. Med. 2009, 61, 103–116. [Google Scholar] [CrossRef] [PubMed]
  32. Jung, H.; Ye, J.C.; Kim, E.Y. Improved k-t BLAST and k-t SENSE using FOCUSS. Phys. Med. Biol. 2007, 52, 3201–3226. [Google Scholar] [CrossRef] [PubMed]
  33. Feng, L.; Srichai, M.B.; Lim, R.P.; Harrison, A.; King, W.; Adluru, G.; Dibella, E.V.; Sodickson, D.K.; Otazo, R.; Kim, D. Highly accelerated real-time cardiac cine MRI using k-t SPARSE-SENSE. Magn. Reson. Med. 2013, 70, 64–74. [Google Scholar] [CrossRef] [PubMed]
  34. Candes, E.J.; Recht, B. Exact matrix completion via convex optimization. Found. Comput. Math. 2009, 9, 717–772. [Google Scholar] [CrossRef]
  35. Otazo, R.; Candes, E.; Sodickson, D.K. Low-rank plus sparse matrix decomposition for accelerated dynamic MRI with separation of background and dynamic components. Magn. Reson. Med. 2015, 73, 1125–1136. [Google Scholar] [CrossRef] [PubMed]
  36. Liu, Q.; Lai, Z.; Zhou, Z.; Kuang, F.; Jin, Z. A truncated nuclear norm regularization method based on weighted residual error for matrix completion. IEEE Trans. Image Process. 2016, 25, 316–330. [Google Scholar] [CrossRef] [PubMed]
  37. Fang, X.; Xu, Y.; Li, X.; Lai, Z.; Wong, W.K. Robust semi-supervised subspace clustering via non-negative low-rank representation. IEEE Trans. Cybern. 2016, 46, 1828–1838. [Google Scholar] [CrossRef] [PubMed]
  38. Lingala, S.G.; Hu, Y.; DiBella, E.; Jacob, M. Accelerated dynamic MRI exploiting sparsity and low-rank structure: k-t SLR. IEEE Trans. Med. Imaging 2011, 30, 1042–1054. [Google Scholar] [CrossRef] [PubMed]
  39. Zhao, B.; Haldar, J.P.; Christodoulou, A.G.; Liang, Z.P. Image reconstruction from highly undersampled (k,t)-space data with joint partial separability and sparsity constraints. IEEE Trans. Med. Imaging 2012, 31, 1809–1820. [Google Scholar] [CrossRef] [PubMed]
  40. Candes, E.J.; Li, X.D.; Ma, Y.; Wright, J. Robust principal component analysis? J. ACM 2011, 58, 1–37. [Google Scholar] [CrossRef]
  41. Yuan, X.; Yang, J. Sparse and low-rank matrix decomposition via alternating direction methods. Pac. J. Optim. 2013, 9, 167–180. [Google Scholar]
  42. Tao, M.; Yuan, X.M. Recovering low-rank and sparse components of matrices from incomplete and noisy observations. SIAM J. Optim. 2011, 21, 57–81. [Google Scholar] [CrossRef]
  43. Wright, J.; Ganesh, A.; Min, K.R.; Ma, Y. Compressive principal component pursuit. Inf. Inference 2013, 2, 32–68. [Google Scholar] [CrossRef]
  44. Trémoulhéac, B.; Dikaios, N.; Atkinson, D.; Arridge, S.R. Dynamic MR image reconstruction-separation from undersampled (k,t)-space via low-rank plus sparse prior. IEEE Trans. Med. Imaging 2014, 33, 1689–1701. [Google Scholar] [CrossRef] [PubMed]
  45. Majumdar, A.; Ward, R.K. An algorithm for sparse MRI reconstruction by Schatten p-norm minimization. Magn. Reson. Imaging 2011, 29, 408–417. [Google Scholar] [CrossRef] [PubMed]
  46. Trzasko, J.; Manduca, A. Highly undersampled magnetic resonance image reconstruction via homotopic l0-minimization. IEEE Trans. Med. Imaging 2009, 28, 106–121. [Google Scholar] [CrossRef] [PubMed]
  47. Majumdar, A. Improved dynamic MRI reconstruction by exploiting sparsity and rank-deficiency. Magn. Reson. Imaging 2013, 31, 789–795. [Google Scholar] [CrossRef] [PubMed]
  48. Zhang, Y.; Sun, X.; Wang, B. Efficient algorithm for k-barrier coverage based on integer linear programming. China Commun. 2016, 13, 16–23. [Google Scholar] [CrossRef]
  49. Mehranian, A.; Reader, A. Non-convex joint-sparsity regularization for synergistic PET and SENSE MRI reconstruction. J. Nucl. Med. 2016, 57, 639. [Google Scholar]
  50. Ding, Y.; Selesnick, I.W. Artifact-free wavelet denoising: Non-convex sparse regularization, convex optimization. IEEE Signal Process. Lett. 2015, 22, 1364–1368. [Google Scholar] [CrossRef]
  51. Zhang, J.; Xiong, R.; Zhao, C.; Zhang, Y.; Ma, S.; Gao, W. CONCOLOR: Constrained non-convex low-rank model for image deblocking. IEEE Trans. Image Process. 2016, 25, 1246–1259. [Google Scholar] [CrossRef] [PubMed]
  52. Chartrand, R. Exact reconstruction of sparse signals via nonconvex minimization. IEEE Signal Process. Lett. 2007, 14, 707–710. [Google Scholar] [CrossRef]
  53. Chartrand, R.; Staneva, V. Restricted isometry properties and nonconvex compressive sensing. Inverse. Probl. 2008, 24, 035020:1–035020:14. [Google Scholar] [CrossRef]
  54. Yuan, C.; Sun, X.; Lv, R. Fingerprint liveness detection based on multi-scale LPQ and PCA. China Commun. 2016, 13, 60–65. [Google Scholar] [CrossRef]
  55. Gasso, G.; Rakotomamonjy, A.; Canu, S. Recovering sparse signals with a certain family of nonconvex penalties and DC programming. IEEE Trans. Signal Process. 2009, 57, 4686–4698. [Google Scholar] [CrossRef]
  56. Liu, R.W.; Shi, L.; Huang, W.; Xu, J.; Yu, S.C.H.; Wang, D. Generalized total variation-based MRI Rician denoising model with spatially adaptive regularization parameters. Magn. Reson. Imaging 2014, 32, 702–720. [Google Scholar] [CrossRef] [PubMed]
  57. Lai, Z.; Xu, Y.; Chen, Q.; Yang, J.; Zhang, D. Multilinear sparse principal component analysis. IEEE Trans. Neural Netw. Learn. Syst. 2014, 25, 1942–1950. [Google Scholar] [CrossRef] [PubMed]
  58. Zhou, X.; Yang, C.; Zhao, H.; Yu, W. Low-rank modeling and its applications in image analysis. ACM Comput. Surv. 2015, 47, 36. [Google Scholar] [CrossRef]
  59. Boyd, S.; Parikh, N.; Chu, E.; Peleato, B.; Eckstein, J. Distributed optimization and statistical learning via the alternating direction method of multipliers. Found. Trends Mach. Learn. 2011, 3, 1–122. [Google Scholar] [CrossRef]
  60. Cai, J.F.; Candes, E.J.; Shen, Z.W. A singular value thresholding algorithm for matrix completion. SIAM J. Optim. 2010, 20, 1956–1982. [Google Scholar] [CrossRef]
  61. Lu, Y.; Lai, Z.; Xu, Y.; Li, X.; Zhang, D.; Yuan, C. Low-rank preserving projections. IEEE Trans. Cybern. 2016, 46, 1900–1913. [Google Scholar] [CrossRef] [PubMed]
  62. Donoho, D.L. De-noising by soft-thresholding. IEEE Trans. Inf. Theory 1995, 41, 613–627. [Google Scholar] [CrossRef]
  63. Daubechies, I.; Defrise, M.; De Mol, C. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Commun. Pure Appl. Math. 2004, 57, 1413–1457. [Google Scholar] [CrossRef]
  64. Zuo, W.M.; Meng, D.Y.; Zhang, L.; Feng, X.C.; Zhang, D. A generalized iterated shrinkage algorithm for non-convex sparse coding. In Proceedings of the IEEE International Conference on Computer Vision, Sydney, Australia, 1–8 December 2013; pp. 217–224.
Figure 1. The low-rank plus sparse (L + S) decomposition of a fully sampled 2D axial cardiac MRI dataset. RPCA decomposes the x-t (or y-t) temporal profiles (i.e., Casorati matrix) F into a low-rank component L and a sparse component S. The singular values of both x-t and y-t L tend to zero quickly as their index increases. S can be effectively represented using a sparsifying transform.
Figure 1. The low-rank plus sparse (L + S) decomposition of a fully sampled 2D axial cardiac MRI dataset. RPCA decomposes the x-t (or y-t) temporal profiles (i.e., Casorati matrix) F into a low-rank component L and a sparse component S. The singular values of both x-t and y-t L tend to zero quickly as their index increases. S can be effectively represented using a sparsifying transform.
Sensors 17 00509 g001
Figure 2. Radial sampling trajectories with N number of uniformly spaced rays in each frame with subsequent random rotations across time frames. From left to right: different number of radial rays ranging from N = 8 , 12 , 16 , 24 and 32 was considered to simulate the undersampled ( k , t ) -space measurements.
Figure 2. Radial sampling trajectories with N number of uniformly spaced rays in each frame with subsequent random rotations across time frames. From left to right: different number of radial rays ranging from N = 8 , 12 , 16 , 24 and 32 was considered to simulate the undersampled ( k , t ) -space measurements.
Sensors 17 00509 g002
Figure 3. The influences of different parameters ( p , q ) on image reconstruction for axial cardiac dataset. The undersampled ( k , t ) -space measurements were generated using 12 radial projections. The optimal values of ( p , q ) were determined using the quality metric SER.
Figure 3. The influences of different parameters ( p , q ) on image reconstruction for axial cardiac dataset. The undersampled ( k , t ) -space measurements were generated using 12 radial projections. The optimal values of ( p , q ) were determined using the quality metric SER.
Sensors 17 00509 g003
Figure 4. Qualitative results for in vivo axial cardiac dataset with 8 radial projections. From top to bottom: time-frame magnitude images, local magnification views (corresponding to the white square on the reference image), x-t temporal profiles, y-t temporal profiles, and 1D profiles (according to the point A on the reference image). Left colormap refers to local magnification views. Right colormap refers to magnitude images and temporal profiles. (The images are best viewed in full-screen mode.)
Figure 4. Qualitative results for in vivo axial cardiac dataset with 8 radial projections. From top to bottom: time-frame magnitude images, local magnification views (corresponding to the white square on the reference image), x-t temporal profiles, y-t temporal profiles, and 1D profiles (according to the point A on the reference image). Left colormap refers to local magnification views. Right colormap refers to magnitude images and temporal profiles. (The images are best viewed in full-screen mode.)
Sensors 17 00509 g004
Figure 5. 1D profiles respectively generated by ZF-IDFT, k-t FOCUSS, k-t SLR, k-t RPCA and k-t NCRPCA for different number of radial projections (i.e., 12, 16, 24 and 32). (The images are best viewed in full-screen mode.)
Figure 5. 1D profiles respectively generated by ZF-IDFT, k-t FOCUSS, k-t SLR, k-t RPCA and k-t NCRPCA for different number of radial projections (i.e., 12, 16, 24 and 32). (The images are best viewed in full-screen mode.)
Sensors 17 00509 g005
Figure 6. Qualitative results for in vivo coronal cardiac dataset with eight radial projections. From top to bottom: time-frame magnitude images, x-t temporal profiles, y-t temporal profiles, and 1D profiles (according to the point A on the reference image). (The images are best viewed in full-screen mode.)
Figure 6. Qualitative results for in vivo coronal cardiac dataset with eight radial projections. From top to bottom: time-frame magnitude images, x-t temporal profiles, y-t temporal profiles, and 1D profiles (according to the point A on the reference image). (The images are best viewed in full-screen mode.)
Sensors 17 00509 g006
Figure 7. Visual appearance of error images respectively generated by ZF-IDFT, k-t FOCUSS, k-t SLR, k-t RPCA and k-t NCRPCA for different number of radial projections (i.e., 8, 12, 16, 24 and 32). (The images are best viewed in full-screen mode.)
Figure 7. Visual appearance of error images respectively generated by ZF-IDFT, k-t FOCUSS, k-t SLR, k-t RPCA and k-t NCRPCA for different number of radial projections (i.e., 8, 12, 16, 24 and 32). (The images are best viewed in full-screen mode.)
Sensors 17 00509 g007
Figure 8. Convergence property of the proposed k-t NCRPCA. Progression of the objective quality metric SER for in vivo axial (left) and coronal (right) cardiac datasets with respect to iteration number. The undersampled ( k , t ) -space measurements were generated from 8 , 16 and 32 radial projections.
Figure 8. Convergence property of the proposed k-t NCRPCA. Progression of the objective quality metric SER for in vivo axial (left) and coronal (right) cardiac datasets with respect to iteration number. The undersampled ( k , t ) -space measurements were generated from 8 , 16 and 32 radial projections.
Sensors 17 00509 g008
Table 1. SSIM/SER comparison of various reconstruction methods on in vivo axial cardiac MRI dataset for different number of radial projections.
Table 1. SSIM/SER comparison of various reconstruction methods on in vivo axial cardiac MRI dataset for different number of radial projections.
ProjectionsZF-IDFTk-t FOCUSSk-t SLRk-t RPCAk-t NCRPCA
80.2555/3.46030.7376/11.20020.8802/12.66340.9476/19.60960.9503/20.1538
120.3088/4.46260.7969/12.99510.9257/15.64410.9760/22.03930.9784/22.9995
160.3608/5.30260.8518/15.27400.9518/17.74280.9769/22.90160.9789/23.7973
240.4408/6.79890.9004/17.03830.9671/20.38970.9828/24.08020.9942/24.9276
320.4981/8.19530.9257/20.47910.9831/22.83580.9870/25.12640.9956/25.8280
Table 2. SSIM/SER comparison of different reconstruction methods on in vivo coronal cardiac MRI dataset for different number of radial projections.
Table 2. SSIM/SER comparison of different reconstruction methods on in vivo coronal cardiac MRI dataset for different number of radial projections.
ProjectionsZF-IDFTk-t FOCUSSk-t SLRk-t RPCAk-t NCRPCA
80.3281/7.31120.7041/16.18930.8328/16.96400.7948/17.57190.8268/17.8084
120.4306/8.75480.8099/17.91730.8894/18.85290.8682/20.21870.8937/20.6879
160.4974/9.94460.8618/19.65150.9300/20.88790.9134/21.92710.9384/22.6999
240.5929/11.99030.9064/21.96810.9546/23.18060.9383/23.53530.9768/25.1972
320.6753/13.94920.9489/24.29310.9757/25.76130.9674/25.47100.9887/27.2034

Share and Cite

MDPI and ACS Style

Liu, R.W.; Shi, L.; Yu, S.C.H.; Xiong, N.; Wang, D. Reconstruction of Undersampled Big Dynamic MRI Data Using Non-Convex Low-Rank and Sparsity Constraints. Sensors 2017, 17, 509. https://doi.org/10.3390/s17030509

AMA Style

Liu RW, Shi L, Yu SCH, Xiong N, Wang D. Reconstruction of Undersampled Big Dynamic MRI Data Using Non-Convex Low-Rank and Sparsity Constraints. Sensors. 2017; 17(3):509. https://doi.org/10.3390/s17030509

Chicago/Turabian Style

Liu, Ryan Wen, Lin Shi, Simon Chun Ho Yu, Naixue Xiong, and Defeng Wang. 2017. "Reconstruction of Undersampled Big Dynamic MRI Data Using Non-Convex Low-Rank and Sparsity Constraints" Sensors 17, no. 3: 509. https://doi.org/10.3390/s17030509

APA Style

Liu, R. W., Shi, L., Yu, S. C. H., Xiong, N., & Wang, D. (2017). Reconstruction of Undersampled Big Dynamic MRI Data Using Non-Convex Low-Rank and Sparsity Constraints. Sensors, 17(3), 509. https://doi.org/10.3390/s17030509

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop