Online dynamic cardiac imaging based on the elastic-net model

Abstract

Purpose: The purpose of this work was to develop an online dynamic cardiac MRI model to reconstruct image frames from partial acquisition of the Cartesian k-space data, which utilizes structural knowledge of consecutive image frames. Materials and methods: Using an elastic-net model, the proposed algorithm reconstructs dynamic images using both L₁ and L₂ norm operations. The L₁ norm enforces the sparsity of the frame difference, while the L₂ norm with motion-adaptive weights catches the internal structure of frame differences. Unlike other online methods such as the Kalman filter (KF) technique, the new model requires no assumption of Gaussian noise, and can faithfully reconstruct the dynamic images within a compressive sensing framework. Results: The proposed method was evaluated using simulated dynamic phantoms with 40 frames of images (128 × 128) and a cardiac MRI cine of 25 frames (256 × 256). Both results showed that the new model offered a better performance than the online KF method in depicting simulated phantom and cardiac dynamics. Conclusion: It is concluded that the proposed imaging model can be used to capture a large variety of objects in motion from highly under-sampled k-space data, and being particularly useful for improving temporal resolution of cardiac MRI.

Keywords:

• Dynamic cardiac imaging
• online reconstruction
• elastic-net model
• MRI
• inverse problem

Introduction

Dynamic MRI is an important imaging technique in clinical applications, such as interventional and functional imaging, that need to capture objects such as heart and intake agent in motion with high spatial and temporal resolutions. Ideally, this can be achieved by acquiring full k-space data of each frame and then reconstructing the image via inverse Fourier transform. In this straightforward method, the temporal resolution hinges on the speed of acquisition, which is limited by physical or physiological constraints. To improve temporal resolution or to speed-up scan procedure, one effective method is to acquire partial k-space data of each frame. In recent years, several strategies [1] have been proposed to reconstruct images from partially acquired k-space data via effective exploitation of spatial and/or temporal correlations of dynamic images. For example, dynamic models have been previously built to exploit temporal correlation between frames. Using the block match algorithm from video compression, k–t FOCUSS (FOCal Underdetermined System Solver) with ME/MC (Motion Estimation/Motion Compensation) [2] estimates the motion by relocating blocks of pixels from a reference frame into current frame and then encodes the sparse residual with k–t FOCUSS [3]; MASTeR (Motion-Adaptive Spatio-Temporal Regularization) [4] constructs a motion-adaptive linear dynamic system and recovers images based on the sparsity of forward and backward differences. However, they both work offline, that is, the acquisition and reconstruction of each frame cannot proceed alternately, because the reconstruction of current frame depends on the information of subsequent frames. As a result, those offline methods limit their usages in interventional applications, such as MRI-guided surgery.

To monitor the motion of objects during scan, several online reconstruction techniques have been proposed. Bilen [5] uses optical flow to estimate the object’s motion, which relocates each pixel from previous frame into current frame, and then recovers images by compressed sensing (CS). However, because both quality reference frame and current frame are not available in the process of imaging as opposed to video compression, the motion estimation by simply relocating pixels from reference frame may incur loss of quality for reconstructed images. To avoid this problem, the Kalman filter (KF) technique has been used to estimate unknown variables recursively from measurements observed over time. In dynamic MRI, the CS scheme is incorporated into the framework of KF to reconstruct the sparse residual [6]. Sümbület [7,8] modelled dynamic MRI as a discrete linear system, and recovered images out of the model frame by frame using KF. However, to make KF computationally tractable, the covariance matrix in KF is simplified by ignoring off-diagonal elements, which can only be possible for non-Cartesian sampling trajectory. For Cartesian trajectory, Feng [9] broke down an image into columns with much lower dimensions along the phase-encoding direction for each read-out position, and then reconstructed each column independently using KF. Hong [10] improved the models proposed in [9] by formulating KF as an optimization programme and then imposing additional sparsity constraints on it. However, in both models of [9] and [10], the system noise of KF model is actually the difference between two consecutive frames, which is usually not Gaussian and hence violates the condition of the KF.

In this work, with the consideration of the fact that the frame difference is caused by object motion and is usually sparse, an elastic-net [11] model is proposed that describes the object dynamics in an online manner. The method reconstructs images using a hybrid optimization procedure that combines the L₁ and L₂ norms, in which L₁ norm is used to promote sparsity and L₂ norm with motion-adaptive weights is employed to catch the internal structure of frame differences. The proposed imaging model does not require Gaussian assumption of the system and measurement noise, as does the KF. Simulated phantom and cardiac MRI cine were used to demonstrate improved performance of the proposed method.

Theory

Overview of the discrete linear system model for dynamic MRI

This section briefly reviews the discrete linear system for dynamic MRI. The temporal correlation of consecutive frames of the imaging objects can be described by a discrete linear system. When the k-space is sampled in a Cartesian trajectory, dynamic MRI can be described by(1a) $${\mathbf{x}}_{\text{t}}={\mathbf{x}}_{\text{t}-1}+{\mathbf{u}}_{\text{t}}$$(1a) (1b) $${\mathbf{z}}_{\text{t}}={\mathbf{F}}_{\text{u},\text{t}}{\mathbf{x}}_{\text{t}}+{\mathbf{v}}_{\text{t}}$$(1b)

where ${\mathbf{x}}_{\text{t}}$ is an image column vector assuming the row vector is along the phase-encoding direction, ${\mathbf{z}}_{\text{t}}$ is the corresponding column vector after applying an inverse Fourier transform to row vectors of acquired k-space data, ${\mathbf{F}}_{\text{u},\text{t}}$ is the under-sampled Fourier transform, ${\mathbf{u}}_{\text{t}}$ is the frame difference, ${\mathbf{v}}_{\text{t}}$ is the measurement noise of MRI scanners. The subscript $\text{t}$ denotes the time index of a frame. For an N by N MR image, ${\mathbf{x}}_{\text{t}}$ and ${\mathbf{u}}_{\text{t}}$ are vectors of size N and ${\mathbf{F}}_{u,t}$ is a matrix of K by N, where K is the number of sampled lines. The relationship of ${\mathbf{x}}_{\text{t}}$ and ${\mathbf{z}}_{\text{t}}$ is shown in Figure 1.

Figure 1. The reconstruction of a 3D cine MRI can be broken into the reconstruction of multiple 2D x–t images along the read-out direction. The symbol ⊙ denotes element-wise multiplication.

Given ${\mathbf{z}}_{\text{t}}$ and ${\mathbf{F}}_{\text{u},\text{t}}$, ${\mathbf{u}}_{\text{t}}$ can be reconstructed by a model formulized in terms of ${\mathbf{u}}_{\text{t}}$ and ${\mathbf{v}}_{\text{t}}$ with a prior structure. The estimate for ${\mathbf{x}}_{\text{t}}$ is then obtained by adding reconstructed ${\mathbf{u}}_{\text{t}}$ to the previous frame ${\mathbf{x}}_{\text{t}-1}$. With different assumptions of the structure of frame difference ${\mathbf{u}}_{\text{t}}$, various modalities had been proposed to recover images from Equation (1). For example, reference [12] modelled ${\mathbf{u}}_{\text{t}}$ as a rank-deficient matrix and images are recovered from the partially sampled k-space data via nuclear norm minimization; references [13,14] observed that ${\mathbf{u}}_{\text{t}}$ is sparse and thus the dynamic frames can be recovered with the help of CS technique; reference [9] assumed both ${\mathbf{u}}_{\text{t}}$ and ${\mathbf{v}}_{\text{t}}$ are ‘Gaussian’ signal, thus KF was applied to recover ${\mathbf{u}}_{\text{t}}$ for Cartesian trajectory. The KF is a recursive algorithm that estimates the state ${\mathbf{x}}_{\text{t}}$ of a process in a way that minimizes the variance of the estimation error when ${\mathbf{u}}_{\text{t}}$ and ${\mathbf{v}}_{\text{t}}$ represent independent white-noise sequences.[16,19] It incorporates a new measurement ${\mathbf{z}}_{\text{t}}$ into the a priori estimate $({\mathbf{x}}_{\text{t}-1}+{\mathbf{u}}_{\text{t}})$ to obtain an improved a posteriori estimate ${\mathbf{x}}_{\text{t}}$.[16] Compared with other filters such as Wiener filter [19], the KF is both memory and computation efficient and is well suited to online/real-time execution.

Specifically, given ${\mathbf{u}}_{\text{t}}\sim \text{N}\left(0,{\mathbf{Q}}_{\text{t}}\right)\text{and}{\mathbf{v}}_{\text{t}}\sim \text{N}\left(0,{\mathbf{R}}_{\text{t}}\right)$, the KF technique produces a recursive solution to the above discrete linear system, that can be described using the following equations [16](2a) $${\mathbf{K}}_{\text{t}}={\mathbf{P}}_{\text{t}}{\mathbf{F}}_{\text{u},\text{t}}^{\text{H}}{[{\mathbf{F}}_{\text{u},\text{t}}{\mathbf{P}}_{\text{t}}{\mathbf{F}}_{\text{u},\text{t}}^{\text{H}}+{\mathbf{R}}_{\text{t}}]}^{-1}$$(2a) (2b) $${\mathbf{x}}_{\text{t}}={\mathbf{x}}_{\text{t}-1}+{\mathbf{K}}_{\text{t}}\left[{\mathbf{z}}_{\text{t}}-{\mathbf{F}}_{\text{u},\text{t}}{\mathbf{x}}_{\text{t}-1}\right]$$(2b) (2c) $${\mathbf{P}}_{\text{t}+1}=\left[\mathbf{I}-{\mathbf{K}}_{\text{t}}{\mathbf{F}}_{\text{u},\text{t}}\right]{\mathbf{P}}_{\text{t}}+{\mathbf{Q}}_{\text{t}}$$(2c)

where ${\mathbf{P}}_{\text{t}}$ is defined to be the estimation error covariance matrix and ${\mathbf{K}}_{\text{t}}$ be the Kalman gain [16] at the moment $\text{t}$. Substituting Equation (2a) into Equation (2b), we have(3) $${\mathbf{x}}_{\text{t}}-{\mathbf{x}}_{\text{t}-1}={\mathbf{P}}_{\text{t}}{\mathbf{F}}_{\text{u},\text{t}}^{\text{H}}{[{\mathbf{F}}_{\text{u},\text{t}}{\mathbf{P}}_{\text{t}}{\mathbf{F}}_{\text{u},\text{t}}^{\text{H}}+{\mathbf{R}}_{\text{t}}]}^{-1}[{\mathbf{z}}_{\text{t}}-{\mathbf{F}}_{\text{u},\text{t}}{\mathbf{x}}_{\text{t}-1}]$$(3)

If the measurement noise covariance matrix ${\mathbf{R}}_{\text{t}}$ is diagonal and the time-invariant, that is ${\mathbf{R}}_{\text{t}}={\mathit{\sigma }}^2\mathbf{I}$, Equation (3) is actually equivalent to the following equation(4) $${\mathbf{x}}_{\text{t}}-{\mathbf{x}}_{\text{t}-1}={\text{argmin}}_{\mathbf{u}}||{\mathbf{z}}_{\text{t}}-{\mathbf{F}}_{\text{u},\text{t}}({\mathbf{x}}_{\text{t}-1}+\mathbf{u}){||}_2+2{\mathit{\sigma }}^2{\mathbf{u}}^{\text{H}}{\mathbf{P}}_{\text{t}}^{-1}\mathbf{u}$$(4)

Because ${\mathbf{P}}_{\text{t}}$ is symmetric and semi-definite, and it can be decomposed as ${\mathbf{P}}_{\text{t}}={(\sqrt{{\mathbf{P}}_{\text{t}}})}^{\text{H}}\sqrt{{\mathbf{P}}_{\text{t}}}$, so Equation (4) reveals that KF minimizes $||{\left(\sqrt{{\mathbf{P}}_{\text{t}}}\right)}^{-\text{H}}\mathbf{u}{||}_2^2$ with the fidelity constraint $||{\mathbf{z}}_{\text{t}}-{\mathbf{F}}_{\text{u},\text{t}}({\mathbf{x}}_{\text{t}-1}+\mathbf{u}){||}_2<\mathit{\epsilon }$.

Elastic-net imaging model

To further reveal the structure of frame difference, we depict one of frame differences with the most non-zero elements for a cardiac cine, as shown in Figure 2, indicating that the frame difference reflects the motion of the heart (Figure 2(a)), and that it is compressible and sparse (Figure 2(b)).

Figure 2. One frame difference of a cardiac cine (a) and its sparsity (b).

In dynamic MRI, large similarity between adjacent frames usually exists; in other words, the differences of them are generally small, as demonstrated by Figure 2(a). So, the following imaging model can be initially set up(5) $${\widehat{\mathbf{u}}}_{\text{t}}={\text{argmin}}_{\mathbf{u}}||{\mathbf{z}}_{\text{t}}-{\mathbf{F}}_{\text{u},\text{t}}({\mathbf{x}}_{\text{t}-1}+\mathbf{u}){||}_2+{\mathit{\alpha }||\mathbf{u}||}_1+{\mathit{\beta }||\mathbf{u}||}_2$$(5)

Equation (5) is a mixed-norm model and is known as the elastic-net [11] in the field of statistics. It can also be rewritten to(6) $${\widehat{\mathbf{u}}}_{\text{t}}={\text{argmin}}_{\mathbf{u}}||\left(\begin{array}{c}{\mathbf{z}}_{\text{t}}-{\mathbf{F}}_{\text{u},\text{t}}{\mathbf{x}}_{\text{t}-1}\\ 0\end{array}\right)-\left(\begin{array}{c}{\mathbf{F}}_{\text{u},\text{t}}\\ \mathit{\beta }\mathbf{I}\end{array}\right){\mathbf{u}||}_2+\mathit{\alpha }{|\left|\mathbf{u}\right||}_1$$(6)

In the original elastic-net model, $\mathit{\beta }$ is a constant scalar. That means all elements of $\mathbf{u}$ are penalized with the same weights towards zero. However, Figure 2(a) shows that the elements of $\mathbf{u}$ should be penalized adaptively to catch the regional variation of the motion of the imaging object. Let ${\mathbf{\beta }}_{\text{t}}$ denote the adaptive weights of $\mathbf{u}$, then Equation (5) can be re-expressed as(7) $${\widehat{\mathbf{u}}}_{\text{t}}={\text{argmin}}_{\mathbf{u}}||{\mathbf{z}}_{\text{t}}-{\mathbf{F}}_{\text{u},\text{t}}({\mathbf{x}}_{\text{t}-1}+\mathbf{u}){||}_2+{\mathit{\alpha }||\mathbf{u}||}_1+\mathit{\gamma }||\text{diag}({\mathit{\beta }}_{\text{t}})\mathbf{u}{||}_2$$(7)

where $\mathit{\gamma }$ is the regularization constant for L₂ norm term, and Equation (7) is equivalent to(8) $${\widehat{\mathbf{u}}}_{\text{t}}={\text{argmin}}_{\mathbf{u}}||\left(\begin{array}{c}{\mathbf{z}}_{\text{t}}-{\mathbf{F}}_{\text{u},\text{t}}{\mathbf{x}}_{\text{t}-1}\\ 0\end{array}\right)-\left(\begin{array}{c}{\mathbf{F}}_{\text{u},\text{t}}\\ \mathit{\gamma }diag({\mathit{\beta }}_{\text{t}})\end{array}\right){\mathbf{u}||}_2+\mathit{\alpha }{|\left|\mathbf{u}\right||}_1$$(8)

Equation (8) is the proposed imaging model and can be solved using various sparse optimization solvers.

Relationship to existing work

In [10], to enforce the sparsity of the KF solution $\mathbf{u}$, Equation (4) is solved with an L₁ norm constraint as(9) $${\widehat{\mathbf{u}}}_{\text{t}}={\text{argmin}}_{\mathbf{u}}||{\mathbf{z}}_{\text{t}}-{\mathbf{F}}_{\text{u},\text{t}}({\mathbf{x}}_{\text{t}-1}+\mathbf{u}){||}_2+\mathit{\alpha }{|\left|\mathbf{u}\right||}_1+\mathit{\beta }{\mathbf{u}}^{\text{H}}{\mathbf{P}}_{\text{t}}^{-1}\mathbf{u}$$(9)

Equation (9) was then solved using FOCUSS [15] and the results show that this sparsity-enforced KF technique statistically improves the original KF.[9]

Comparing Equation (9) with Equation (7), we show that the difference between the KF technique and the proposed method is that the KF weights $\mathbf{u}$ using the inverse of ${\mathbf{P}}_{\text{t}}$, while the proposed method uses an estimated vector of weights. Specifically, Equation (9) is different from Equation (7) in two aspects:

(1)

Equation (7) does not require $\mathbf{u}$ to be Gaussian, nor the measurement noise ${\mathbf{v}}_{\text{t}}$ to be white Gaussian with a constant variance. In dynamic MRI, $\mathbf{u}$ is actually the frame difference and is largely due to object motion,[13] which very unlikely is random Gaussian. This violates the condition of the KF.

(2)

The ${\mathbf{P}}_{\text{t}}$ in Equation (9) is derived by minimizing the a posteriori error covariance[16] and is updated independently from the state variables by Equation (2a) and Equation (2c). It is difficult to interpret ${\mathbf{P}}_{\text{t}}$ from Equation (2) and how the frame difference is weighted. While in Equation (7), ${\mathbf{\beta }}_{\text{t}}$ is known by design to penalize the $\mathbf{u}$ towards zero.

Reconstruction

To solve Equation (8), ${\mathbf{\beta }}_{\text{t}}$ should first be estimated. As stated earlier, ${\mathbf{\beta }}_{\text{t}}$ is a vector of weights to penalize the element of ${\mathbf{u}}_{\text{t}}$ towards zeros. Obviously, stationary voxels should be penalized more heavily than moving voxels. We observe that the variance of moving voxels is usually larger than stationary voxels. Therefore, the variance of one voxel from first frame to previous frame indicates whether the voxel is stationary or moving. So, it is reasonable to use variance of the voxels to derive an estimation for ${\mathbf{\beta }}_{\text{t}}$. To suppress the impact of very early frames, only $\text{W}$ frames earlier than the frame being recovered are used to calculate ${\mathbf{\beta }}_{\text{t}}$. In addition, the value of ${\mathbf{\beta }}_{\text{t}}$ is normalized to (0, 1] to prevent amplifying the noise. For the first frame and because no variance is available, the CS-MRI [17] is used to reconstruct the image from under-sampled k-space data. Starting from the second frame, ${\mathbf{\beta }}_{\text{t}}$ is calculated using variance of the voxels, which will weight adaptively to catch the signal structure from a priori knowledge.

Given ${\mathbf{\beta }}_{\text{t}}$, Equation (8) is within the framework of the compressive sensing (CS). Among other requirements, the CS requires an L₁ solver. In this work, the Gradient Projection for Sparse Reconstruction with Barzilai and Borwein approach (GPSR-BB) [18] version 5.0 is used to recover the difference out of Equation (8). It is initialized with 0 and stops when the norm of the difference between two consecutive estimates divided by the norm of one of them falls below 10⁻⁵.

Methods

The proposed model is implemented in the CS framework and is verified by using a simulated phantom and cardiac MRI cine. In this work, the k-space is under-sampled using a Gaussian random pattern along the phase-encoding direction. Because our model reconstructs current frame by reconstructing the frame difference with previous frame, the reconstruction quality of the first frame will provide a good foundation for subsequent frames. During acquisition, a denser sampling pattern is implemented for the earlier frames and less data was collected in the later ones. Specifically, the first frame was under-sampled by a factor of two, the second frame by four and the rest of the frames by eight. And ${\mathbf{\beta }}_{\text{t}}$ is the normalized variance of $\text{W}=6$ frames earlier than the frame being recovered because $\text{W}=6$ yields better results overall. The state-of-the-art online method [9] and the offline method k–t FOCUSS [3] were compared with the proposed method, using the same under-sampling patterns. Theoretically, offline methods demonstrate a better performance than online ones because they use the information from all under-sampled data when recovering one frame. All simulations and reconstructions are implemented in MATLAB R2010b (MathWorks, Inc., Natick, MA).

Simulated phantom

A numerical phantom data-set was designed to reflect overall image structure in an axial MRI image of the heart. This data-set was generated using the same method as [9] by constructing a dynamic image series containing three pairs of concentric circles with a fixed radius with only changing grey level, a slowly oscillating radius and a rapidly oscillating radius, respectively. The circle with fixed radius simulates imaging background, and the slowly oscillating circle simulates respiration, while the rapidly oscillating circle simulates heartbeat. These image series are then transformed by the Fourier transform and are under-sampled randomly in the phase-encoding direction. This data-set contains 40 frames of image with the size of 128 × 128. For this simulated data-set, we set $\mathit{\alpha }=0.05$,$\mathit{\gamma }=0.01$ for Equation (8).

Cardiac MRI cine

For cardiac cine data, 25 frames of full k-space data were acquired using a 1.5T Philips scanner with a 345 × 270 mm FOV and the matrix size is 256 × 256. The data was obtained using a steady-state free precession sequence with a flip angle of 50 degrees and TR = 3.45 ms. This data-set covered one cardiac phase of the diastole and systole. In this work, the full k-space reconstruction is treated as ‘ground truth’. For this cardiac data-set, we set $\mathit{\alpha }=0.05$, $\mathit{\gamma }=0.02$ for Equation (8).

Results

Simulated phantom

Figure 3(a) shows the first frame of the simulated phantom. The circle at the bottom-right is stationary with only the grey level (concentration) changed; the one near the centre of the FOV is slowly oscillating, while the top-left is rapidly oscillating. We will inspect the area enclosed by dashed rectangle to verify how well the methods catch the rapid oscillation. For this simulated data-set, because the data is generated periodically and fits the assumption of k–t FOCUSS which enforces periodic motion, the k–t FOCUSS outperforms the Kalman and proposed method for most frames.

Figure 3. One frame of simulated phantom (a) with a rectangle ROI and the Gaussian random under-sampling pattern along phase-encoding direction (b).

Figure 4 shows results of one read-out position crossing the centre of a rapidly oscillating circle. All three methods capture overall dynamics of the motion, while all smooth out the minor grey level change at frames 6 and 31 marked by two arrows in Figure 4(b). The k–t FOCUSS demonstrates a better performance than the KF and proposed method as expected. However, k–t FOCUSS exhibits periodic pattern near frame 20 in recovered images, because it minimizes the L₁ norm of the x-f space using FOCUSS to enforce periodic pattern in x–t space. The absolute error maps show that proposed method slightly out-performs KF, especially in the area of bright inner circle (marked by arrows in Figure 4(d and e)).

Figure 4. The dashed line in (a) shows one read-out position across the rapidly oscillating circle, (b) is the truth of the x–t space at the read-out position, and (c,d,e) is the recovered x–t space(left half) and their absolute error maps(right half) of the k-t FOCUSS, KF and proposed method.

To evaluate how the proposed method catches slowly oscillating circle and stationary circle with only the grey level changed, we inspected another read-out position covering both slowly oscillating and stationary circles, as shown in Figure 5. From this figure, we can see that all three methods catch slowly oscillating circle well. However, for stationary circle with changing grey level, these methods fail to catch the variation. The k–t FOCUSS smoothens out almost all grey level changes, especially for frames 6 and 31 marked by two arrows in Figure 5(c). It can be seen from the absolute error map that proposed method achieved a smaller error than the KF, especially in areas marked by arrows in Figure 5(d and e).

Figure 5. The dashed line in (a) shows one read-out position across both the slowly oscillating and the stationary circles, (b) is the truth of the x–t space at the read-out position, and (c,d,e) are the recovered x–t space(left half) and their absolute error maps(right half) of k–t FOCUSS, KF and the proposed method.

To show the overall error of recovered images, Figure 6(a) shows the normalized mean squared error ($\text{NMSE}=\frac{||\mathbf{x}-{\mathbf{x}}_{true}{||}^2}{||{\mathbf{x}}_{true}{||}^2}$) of rapidly oscillating circle as marked by dashed rectangle in Figure 3(a). It can be seen that proposed method outperforms KF for all frames but 23 and 24. The k–t FOCUSS has a large NMSE for frame 31 because it cannot catch grey level changes of the inner bright circle as shown in the Figure 6(d).

Figure 6. The NMSE of the area covering the rapidly oscillating circle (a) and the truth, recovered images and their error maps of k–t FOCUSS, KF and proposed method for frames 6(b), 14(c), 31(d), respectively.

Cardiac MRI cine

Figure 7(a) shows one frame of the cardiac cine, in addition to the region of interest (ROI) covering both left and right ventricle cavities.

Figure 7. One frame of cardiac cine (a) with a rectangle ROI and the Gaussian random under-sampling pattern (b).

Figure 8(c) shows that k–t FOCUSS smoothens out some details (marked by black arrows) of the motion, because it enforces periodic pattern in its imaging model. This causes the k–t FOCUSS does not capture the motion of cardiac systole corresponding frames of 9–14. Visually, KF and proposed method had very similar results, except the areas marked by white arrows in Figure 8(d and e). In addition, the KF had a constant stripe in the position marked by the black arrow in Figure 8(d), while the proposed method had a similar pattern at the beginning of the imaging (corresponding to a cardiac diastole with a slow motion activity).This gradually disappears as time passes, as shown in Figure 8(e).This is because the KF updates weights using only the under-sampled Fourier transform and constant covariance matrices of the system and measurement noise, while proposed method corrects weights adaptively, using previous frames to incorporate the motion information.

Figure 8. (a) shows one read-out position across area of rapid cardiac deformation and (b) is the ground truth of the x–t space at the read-out position,(c,d,e) show the recovered x–t space and their error maps of the three methods.

Figure 9 shows the overall NMSE and SSIM (Structural SIMilarity) [20] of the ROI. It can be seen that the k–t FOCUSS had a smaller NMSE at the beginning and end of imaging (corresponding to a cardiac diastole), while KF and proposed method had smaller errors in the middle of the imaging (corresponding to a cardiac systole). This may also be observed in SSIM of Figure 9(e). Compared with KF, proposed method had a smaller NMSE for all frames but 13. However, the SSIM reveals that proposed method outperforms the KF gradually as time goes by. This verifies that the proposed model absorbs the motion information gradually to improve the reconstruction.

Figure 9. The NMSE (a) and SSIM (e) of the ROI and the ground truth, recovered images and their error maps of three methods for frames 8(b), 13(c), 24(d), respectively.

Figure 10 compares how three methods catch the motion of the heart. The k–t FOCUSS catches the motion better in the early and late stages of cardiac cycle, corresponding to slow motion of cardiac wall; however, this method offers a limited capability in depicting rapid contraction of the ventricle wall (see the 7–13 frames in Figure 10(b)). The Kalman and proposed method have a very similar error and follow the ground truth well (see Figure 10(c)). In comparison, without any restrictions in the periodic motion model set-up, proposed method can more accurately capture the motion activity over entire cardiac cycle (see Figure 10(d)).

Figure 10. The average signal intensity of a small region of size 3 × 3 marked in (a), (b,c,d) compares the three methods against the ground truth.

Discussion

In this work, a new online model was developed for rapid dynamic MRI using sub-sampling of the Cartesian k-space data. Previously, dynamic MRI has been investigated using the KF method; however, KF requires the system noise to be Gaussian, which is not suitable for modelling dynamic MRI. In dynamic MRI, the system noise is associated with the frame difference, that is, it contains the motion pattern that is usually not in the form of Gaussian distribution. With a basic assumption of measurement noise is time-invariant and white, the KF technique implements a weighted least square optimization to update the weighting matrix. Noticing that the difference of adjacent frames is sparse and the weighting scheme of KF cannot follow the motion of imaging target, we set up an elastic-net model using L₁ norm to promote sparsity and weighted L₂ norm to catch the motion adaptively.

In the KF method, the weighting matrix is determined by under-sampled Fourier transform, covariance of system and measurement noise, which are chosen to be constant all the time and are updated independently with state variables of the model. While in the proposed method, the weighting matrix is calculated according to the variance of each voxels. As a result, the voxels with larger variance will get reconstructed under sparsity constraints, while voxels with smaller variance will be dominated by previous frame. From comparative results shown above, we can see that proposed method offers a better performance in terms of depicting cardiac dynamics.

The proposed model is solved within the framework of compressive sensing, which benefits from several efficient algorithms. Nonetheless, the Cartesian trajectory enables concurrent reconstruction of all x–t images along read-out direction, which can further reduce the reconstruction time of a frame to only one column by allocating one unit of parallel computing hardware to each read-out position, thus making real-time imaging possible.

The limitation of the proposed model comes from three aspects. First, the model estimates current frame from previous frame using ${\mathbf{x}}_{\text{t}}={\mathbf{x}}_{\text{t}-1}+{\mathbf{u}}_{\text{t}}$, yielding the frame difference that is not Gaussian. If current frame can be better estimated from previous frames using ${\mathbf{x}}_{\text{t}}={\mathbf{Mx}}_{\text{t}-1}+{\mathbf{u}}_{\text{t}}$, where $\mathbf{M}$ is a matrix for motion estimation, to make the frame difference be nearly Gaussian, KF can be applied naturally. However, the motion estimation will incur other possibilities and complexity. Second, it is observed that, the current estimation of adaptive weights by calculating the variance is intuitive and straightforward; however, it may be not robust against the background noise. Third, the proposed model does not accurately represent the real cardiac motion and deformation. And the parameters in the model may contain parametric errors,[21,22] which may be suppressed by the model selection statistically.

Conclusion

This paper presented an online dynamic cardiac MRI model, which can be used to reconstruct image frames from partial acquisition of the k-space data by combining motion-adaptive, sparse constraints. Simulation results showed, compared with online KF method, the proposed method offered a better solution in depicting simulated phantom and cardiac dynamics. Future work will seek better methods for motion estimation and adding a prior structural information for higher spatial and temporal resolution. In addition, the efficiency of the proposed imaging method will be enhanced using hardware-based parallel-computing scheme, which is critical for dynamic MRI studies.

Acknowledgment

This work was supported by National Key Lab of Fundamental Science of Micro/Nano-Device[2014MS01] and Fundamental Research Funds for the Central Universities[CDJZR12090003]

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Fundamental Research Funds for the Central Universities [grant number CDJZR12090003].