|Home | About | Journals | Submit | Contact Us | Français|
To improve myocardial perfusion MRI by reconstructing undersampled radial data with a spatiotemporal constrained reconstruction method (STCR).
The STCR method jointly reconstructs all of the time frames for each slice. In seven subjects at rest on a 3T scanner, the method was compared to a conventional(GRAPPA) Cartesian approach.
Increased slice coverage was obtained as compared to the Cartesian acquisitions. On average, 10 slices were obtained per heartbeat for radial acquisitions (8 of which are suitable for visual analysis and the remaining 2 slices may in theory be used for quantitative purposes), while 4 slices were obtained for the conventional Cartesian acquisitions. The new method was robust to interframe motion, unlike using Cartesian undersampling and STCR. STCR produced images with an image quality rating (1 for best and 5 for worst) of 1.7±0.5; the Cartesian images were rated 2.6±0.4, (p=0.0006). A mean improvement of 44% (±17%) in SNR and 46% (±22%) in CNR was observed for STCR.
The new radial data acquisition and reconstruction scheme for dynamic myocardial perfusion imaging is a promising approach to obtain significantly higher coverage and improved signal to noise ratios. Further testing during vasodilation in patients with coronary artery disease is warranted.
One of the major hurdles in the widespread clinical use of dynamic contrast enhanced myocardial perfusion imaging is the need to track the rapid kinetics of contrast agent uptake. The short data acquisition time available for each image often results in limited coverage of the heart, limited spatial and temporal resolution, or poor signal to noise ratio (SNR) images.
A number of approaches have been devised to address these challenges related to short acquisition time. UNFOLD (2,3), TSENSE (4) and k-t BLAST/k-t SENSE (5, 6), which exploit dynamic and/or parallel imaging, have all been applied to myocardial perfusion imaging. For example, the k-t SENSE method developed in (7) jointly reconstructs all of the time frames for each slice and was recently applied by Plein et al. (6) to actual undersampled data acquired from the scanner with a net acceleration of R=3.9. The acceleration was used to increase spatial resolution and the results were compared to the SENSE reconstruction method with R=2. It was found that dark rim artifacts (8) were reduced and also that the images had improved SNR with k-t SENSE. A significant limitation of the method is its sensitivity to respiratory motion (5, 6).
Another issue that arises for all of these types of speedup methods is that as readout time per image drops, a higher and higher percentage of the acquisition time is spent waiting for sufficient saturation recovery to allow reasonable SNR. One way to circumvent this issue is simply to acquire higher resolution images, so that the number of phase encodes or total readout time does not decrease significantly. Another approach is to image multiple slices after each saturation pulse. Two papers have used this type of approach. One of the methods acquired two slices after every saturation pulse, rather than one (4). The readouts for each line of the two slices were interleaved, so that it can take twice as long to reach the center of k-space, giving a longer effective saturation recovery time for each slice. This approach was implemented with TSENSE and acquired four slices per heartbeat at ~140 bpm.
In a similar vein, multiple slices were acquired sequentially after each saturation pulse in (1). The parallel imaging technique SENSE (R=2) was used and four short axis slices were acquired after a single saturation pulse each heartbeat. The slices had different SRTs and thus different contrasts and SNRs. SNR was also penalized by the use of the SENSE reconstruction. Good sensitivity/specificity of 88/82% was reported in 100 patients studied.
Constrained reconstruction techniques, also termed compressed sensing methods (9), have been used recently to accelerate myocardial perfusion imaging (10–13). The Temporally Constrained Reconstruction (TCR) method (11) was applied to Cartesian data undersampled in a variable density fashion with an acceleration of R=5. Good quality images were obtained from datasets with little motion, but artifacts arose when respiratory motion was present. Similar findings were reported in (13) with an L1 norm constraint in x-f space.
Here we propose a method that is relatively robust to respiratory motion and provides increased slice coverage by acquiring multiple slices of undersampled radial k-space data after every saturation pulse. The undersampled data are reconstructed with a spatiotemporal constrained reconstruction (STCR) method. Simulations and patient studies were performed to test the radial acquisition and reconstruction schemes. Results in human studies at 3T were compared to a standard (GRAPPA, R~1.7) Cartesian data acquisition.
In the standard turboFLASH saturation recovery sequence a single non-selective saturation pulse is applied and k-space data are acquired for each time frame for a given slice using slice-selective α (readout flip angle) pulses. This process is repeated for other slices in a given R-R interval. Although this results in acquiring more than one slice per heart beat, obtaining extended coverage is limited due to the time to perform multiple saturation pulses and due to the waiting time after every saturation pulse that is needed to obtain a reasonable signal before acquiring central k-space data.
To reduce the total preparation time and obtain increased coverage in a given R-R interval, we acquire multiple slices for every saturation pulse as shown in Figure 1. This results in different slices having different signal intensities and contrasts due to different saturation recovery times. As shown in Fig 1, the order of data acquisition for the slices is interleaved between the two saturation pulses in a given R-R interval. This reduces the cross-talk between physically adjacent slices that occurs with practical slice-select profiles. 24 rays in k-space equally spaced over 180° are acquired for a given time frame and a given slice. The 24 rays are acquired in four subsets; for example, for the first time frame the first set of 6 rays were acquired at angles 0, 4π/24, 8π/24,… the second set of 6 rays are acquired at 2π/24, 6π/24, … the third set of 6 rays are acquired at π/24, 5π/24, … and the fourth set of 6 rays are acquired at 3π/24, 7π/24, …Fig 2a shows this pattern. For subsequent time frames in the same slice, the set of 24 k-space rays are acquired in an interleaved fashion with a period of four, that is, for the first time frame k-space rays are acquired at angles 0, π/24, 2π/24, … for the second time frame data are acquired at angles π/96, 5π/96, … and so on.
Noise free simulations were performed using simple simulated left and right ventricular blood pools and myocardium to determine the effect of inconsistent k-space rays on the imaging method. 24 rays in k-space were generated in interleaved (in four subsets) fashion as shown in Fig. 2a. It is well known that interleaved acquisitions of rays help in obtaining better images than sequential acquisitions (14, 15). Radial lines in k-space were generated by taking into account the changing signal due to relaxation of magnetization after every readout α pulse (16, 17) according to equation .
In the above equation, n is the counter for the readout α pulses and TD is the time between the saturation pulse and the first readout pulse. The parameters used in the simulation for the pre-contrast case were M0=1000, T1 for the myocardium=1100 msec, T1 for the blood pools=1500 msec, TD=10 msec, SRT=40 msec, TR=2.5 msec, α =120, T2*=50 msec, TE=1.4 msec. The acquisitions were simulated to model the signal when contrast agent was present. T1 values for the myocardium and blood pools were modified according to equation  and using relaxivities given in (18) by Rohrer et al.
In equation  r is the relaxivity of Gd and its value is 5.5 mM−1s−1 for the blood and 4 mM−1s−1 the myocardium. T1pre is the pre-contrast T1 and T1obs is the observed T1 and [Gd] is the contrast agent concentration.
The radial lines in k-space created by the simulations were reconstructed with a spatially constrained reconstruction (SCR). The cost function for the constrained reconstruction method SCR consisted of a TV constraint only in the spatial dimension as there was no dynamic sequence of images. Although this is different from STCR, the fundamental idea of doing a TV constrained reconstruction is the same. It was shown by Candes et al. in (19) that a TV spatial constraint can reconstruct piecewise constant images faithfully from static undersampled radial data when there are no inconsistencies.
To illustrate the effect of respiratory motion that is often present over different time frames in the perfusion images, a fully sampled Cartesian k-space dataset was used that had complex motion of the heart and chest wall moving in different directions. The dataset was acquired using a standard turboFLASH sequence on a MAGNETOM Trio 3T (Siemens Healthcare AG) with TR = 1.8 msec, TE = 1 msec, flip angle = 12°, Gd dose = 0.025 mmol/kg, slice thickness = 6 mm, acquisition matrix 192 × 96. To simulate the undersampled radial approach proposed here, the data was undersampled offline using approximate radial lines on a Cartesian grid. These radial masks were rotated in an interleaved fashion for different time frames, as described in the pulse sequence section above.
Undersampled radial data with the proposed acquisition scheme was acquired for seven volunteers on the 3T Trio scanner. Six males and one female (age 61.3 (±22.4) years, weight 79.3 (±19.9) kg) were imaged. Four of the volunteers were normal, one had a previous coronary bypass, another had an artificial mitral valve, and one other had a history of cardiac problems. All studies were done specifically for this research and the patient data was acquired in accordance with the XXXX Institutional Review Board. The data was acquired using 9–12 coils selected from body and spine receiver arrays. The scanner parameters for the radial acquisition were TR = 2.6 msec, TE = 1.1 msec, flip angle = 12°, Gd dose = 0.03–0.04 mmol/kg, slice thickness = 6 mm. Reconstructed pixel size varied between 1.8 × 1.8 mm2 and 2.5 × 2.5 mm2. Each image was acquired in a ~62 msec readout, with radial field of view (FOV) ranging from 230 mm – 320 mm. The acquisition matrix size for an image frame was 256 × 24, including an oversampling factor of two in the readout direction. This corresponds to an acceleration factor of ~16 as compared to the Nyquist limit. But the image quality was good with little streaking when 96 projections were used with a standard Inverse Fourier Transform (IFT) reconstruction after gridding. So the current acquisition corresponds to an effective acceleration factor of four. Gradient timing errors were corrected using the method described in (20) by Speier et al. Resting perfusion data with radial undersampling was obtained in which patients were instructed to breathe shallowly.
Cartesian data was also obtained for all the seven patients for comparison. Gd doses for six of the patients were the same as that for the radial acquisition while a reduced dose of 0.015mmol/kg was used for one of the patients. The order of radial and Cartesian data acquisitions for different patients was random, with time between contrast injections on the order of 15 minutes. Cartesian data for four slices was acquired with R = 1.7 and reconstructed with GRAPPA. The parameters for the Cartesian scans were TR = 1.8 msec, TE = 1.1 msec, SRT = 100 msec, flip angle = 12°, slice thickness = 6.0 mm, acquisition matrix size before undersampling=192 × 116. Reconstructed pixel size varied between 1.9 × 1.9 mm2 and 2.1 × 2.1 mm2. FOV varied between 350 × 263 mm2 and 400 × 300 mm2.
After acquiring undersampled radial data for each slice, reconstruction was performed by iteratively minimizing a cost function (11). The cost function consists of a data fidelity term and total variation constraint terms, initially proposed by Rudin et al (21), in the temporal and spatial dimensions. The cost functional C is:
The first term in the above functional is the data fidelity term which quantifies the misfit between estimated solution and acquired data. W is the binary undersampling pattern that is used to obtain the undersampled data, F represents the 2D Fourier transform applied on each time frame in the dynamic sequence, is the estimated image data, is the acquired k-space data and ‖·‖2 is the L2 norm. The second term represents the temporal TV constraint in which α1 is the weighting factor for the constraint, t is the temporal gradient operator, ‖·‖1 is the L1 norm and ε is a small positive constant to avoid singularities in the derivative of the functional as shown by Acar et al in (22). The third term is the spatial TV constraint in which α2 is the weighting factor for the spatial TV constraint, x is the spatial gradient of the image in the x-direction, y is the spatial gradient of the image in y-direction. The TV constraints, in particular the temporal constraint, help in resolving the artifacts from undersampling while preserving sudden changes in the signals and improving SNR. Also, from a compressed sensing point of view, TV constraint exploits the implicit sparsity in the data (19), as the constraint corresponds to an L1 norm of finite differences.
Undersampled radial data was first interpolated onto Cartesian grid points that were within 0.5 units of a measured sample and eq.  was computed using the interpolated data . The value of α1 was chosen using the L-curve technique proposed by Hansen et al. (23) by using only the temporal constraint as in (11) in two datasets and was fixed for other datasets from different patients and from different coils. α2 was then chosen (by fixing the value of α1) as the value that minimized the total absolute error in the reconstruction as compared to full data reconstruction in simulated undersampled Cartesian data and was fixed for the acquired radial datasets. Based on the work in (11, 12) it is likely that the optimal weighting factors for the constraints do not change significantly for these types of perfusion datasets.
The constrained reconstructions obtained for acquisitions with varying Gd concentrations and flip angles were compared in terms of normalized reconstruction errors. The reconstruction error was computed as the normalized root mean square (RMS) error of the signal intensities between the reconstructed image (using the constrained reconstruction) from undersampled data and the original fully sampled image with signal corresponding to the acquisition of the middle (12th) α pulse readout and no data inconsistencies. The fully sampled and the reconstructed images for different cases were normalized to have the same signal intensity ranges, as the signal intensity is changing due to changing the [Gd] and flip angle parameters.
Undersampled patient data was reconstructed offline using STCR for each coil. The reconstructed data from multiple coils were combined in square root of sum of squares (SoS) fashion. Reconstructions from undersampled radial data were compared to the standard Cartesian acquisitions in terms of visual analysis. The reconstructions were scored independently by a cardiologist (XX) with eight years and a radiologist (XX) with three years of experience in reading cardiac perfusion MRI for image quality for both radial and Cartesian acquisitions on a scale from 1 to 5, with 1 being the best image quality. Two closely matching slices for each acquisition per patient were chosen and were scored based on evaluating the left ventricle.
The radial and Cartesian images were normalized to have the same signal intensity ranges and SNR was approximated for a post contrast time frame (visually picked when the signal is peak in the myocardium) as the ratio of the mean signal from a uniform region in the myocardium to the standard deviation of intensities from the same region. CNR was computed as (Myopost – Myopre)/σ, where Myopost is the mean signal in the myocardium in a post contrast time frame and Myopre is the mean signal in the myocardium in a pre contrast time frame. σ is the standard deviation of intensities from the region in the myocardium of the pre contrast time frame. While this scheme of approximating the SNR gives an insight to compare the quality of reconstruction, other methods can be used to create SNR maps reflecting spatially varying noise (24).
Figure 2 shows the results obtained for the simulated heart image. Figure 2b shows the image reconstructed using the constrained method from 24 k-space lines acquired in interleaved fashion (SRT=70msec and no contrast agent). Although not shown here the artifacts were worse with sequential acquisition of the k-space rays. Figure 2c shows the plot comparing normalized error in reconstruction from undersampled data with an SRT of 40 msec using the constrained method for different Gd concentrations and flip angles. The loge(normalized error) in the reconstruction is small, as evident in Figure 2b, and is reduced with increase in the Gd concentration and with increase in the flip angle.
Figure 3 shows an example from an artificially undersampled dataset with significant inter-frame respiratory motion already present in the images (approximate peak shift of 13 mm in column-direction and an approximate peak shift of 5 mm in row-direction over all of the time frames). Figure 3c shows the STCR reconstruction from R~4 data with Cartesian undersampling. The Cartesian data was undersampled in a variable density fashion in which 18 lines in the center of k-space were sampled for each time frame and the remaining data was sampled in a pseudo random fashion to give a net R value of four. The undersampling artifacts are not fully resolved in Fig. 3c. The RMS error for the reconstruction in Fig 3c is 47% higher than that for the reconstruction in Fig 3b as compared to Fig 3a. The mean signal intensity time curves in Figure 3d show that the radial undersampling case matches more closely with the full data reconstruction as compared to the Cartesian case.
Figure 4 shows the results obtained at one time frame for the undersampled radial acquisition and STCR method for a typical subject with rest perfusion at a heart rate of 80 bpm and shallow breathing. The figure shows the reconstructed images of all the 10 slices that were acquired in a single heartbeat during the perfusion sequence. Five slices were acquired after the first saturation pulse and five were acquired after the second saturation pulse. Slice numbers and SRTs (in msec) are shown on the images. The order of acquisition of the slices is interleaved between two saturation pulses and hence the signal intensities are increasing for slices 1, 3, 5, 7 and 9 and similarly for slices 2, 4, 6, 8, 10. A video showing the perfusion images for all of the reconstructed slices for another patient is presented online.
Figure 5 compares the radial reconstructions with the corresponding standard Cartesian acquisitions for closely matching slices for different patients. Patient number, the corresponding Cartesian acquisition and its matching radial slice are shown. Two slices per patient are shown. Image quality for the Cartesian and radial acquisitions is comparable except for the two slices from each radial acquisition which had the lowest saturation recovery time (not shown). A video comparing radial and Cartesian acquisitions for a closely matching slice for a patient is presented online.
The image quality mean scores for the radial and Cartesian acquisitions for the seven patients were significantly different (p=0.0006). The image quality score for the radial acquisitions was better than the Cartesian acquisitions, 1.7 (±0.5) vs. 2.6 (±0.4). A mean improvement of 44% (±17%) in SNR and 46% (±22%) in CNR was observed for the radial reconstructions as compared to the Cartesian ones.
An undersampled radial acquisition with 24 rays per image and a single saturation pulse for five slices was reconstructed with STCR to obtain increased coverage of the heart for myocardial perfusion. In this work we chose to increase the coverage by acquiring more short-axis slices but this can be traded for higher resolution (though resolution along each k-space ray is relatively high) or obtaining more slices with different orientations. The proposed acquisition results in different slices having different signal intensities and contrasts for a given time frame due to different saturation recovery times as implied by Fig. 1. The first slice that is acquired immediately after each saturation pulse, for example slice #1 and slice #2 in Fig. 1, have relatively poor image quality as compared to the other slices due to the low signal and more inconsistent projection data associated with their short saturation recovery times. These slices could be discarded or might be used to obtain a quantitative arterial input function (AIF) from the left ventricular blood pool. AIF from these short SRT slices is likely not saturated and has been used to obtain more accurate semi-quantitative and quantitative parameters (1, 25). The remaining eight slices can be used for both visual and quantitative analyses.
SRTs used here ranged from approximately 35–300 msec, with 90–300 msec being used in the visual assessment. Too long of a SRT could give inadequate T1-weightings. The range of SRTs used successfully in (1), from 42–410 msec, implies these SRTs are not too long.
Of all the slices from 7 patients, slices from two patients had the best image quality (ex. Radial slices in Fig 5 – Patient #3) and slices from the other five patients had slightly lower image quality (ex. Fig 4). The positioning and the choice of which coil data are used to combine in SoS fashion from multiple coils plays an important role in the reconstructed image quality as some of the coils may not have any signal from the heart and only contribute to increased streaking. This type of coil data can be discarded from the SoS combination. Here the coil reconstructions which only contributed to streaking were manually (based on visual inspection of signal present in the heart region after reconstruction) excluded from the SoS combination.
From Fig 2 we see that the normalized reconstruction error decreases with increasing Gd concentration and flip angle. This is interesting as one might guess that with increasing [Gd] the signal changes more rapidly and hence would lead to increased error (as the reconstruction error is related to inconsistencies between radial rays in k-space). But with the presence of even a small time delay (TD) between the saturation pulse and the first readout α pulse, signal inconsistencies are smaller than might be expected. Figure 6a shows the changing signal intensities for the blood pool region in the simulated heart image as 24 rays in k-space are acquired for the first slice and for different Gd concentrations. … The signal is increasing as the concentration is increased and also as radial rays in k-space are acquired. But when the signal curves for different Gd concentrations are normalized so that all the curves start from the same value as shown in Fig 6b, we see that the signal variation is decreasing with increasing Gd concentration. To quantify the amount of variation in the normalized curves, gradient of each curve was computed using finite differences and its L1 norm was computed. The L1 norms of the gradients for the curves are shown next to each curve in Figure 6b. The gradients are decreased with increasing [Gd]. Although not shown here, a similar trend in the variation of the normalized signal was observed with increasing the readout flip angle for a given Gd concentration. The source of this error in the reconstruction is only due to non-uniformity in k-space as these simulations are noise free.
In this study, the patients were instructed to breathe shallowly as opposed to holding their breath during the scan. Holding their breath for part of the scan often leads to much deeper breath later on resulting in large sudden changes in time. The radial acquisition scheme and reconstruction method were found to be much more robust to modest respiratory motion than similar undersampled Cartesian STCR methods as shown in Fig. 3. This robustness to respiratory motion is an important finding and is largely due to the radial sampling scheme in which the center of k-space, where there is most energy, is densely sampled in multiple directions and so the fidelity term helps in preserving sudden changes in time due to respiratory motion. Also, using an L1 as the norm in the temporal constraint helps in better preserving sudden changes in time due to motion as compared to using an L2 norm. When there are no sudden changes in time (no motion), using the L2 norm gave results nearly identical to those obtained using an L1 norm (26). Using a spatial constraint along with the temporal constraint further helped to reduce streaking and resulted in better reconstructions. This is different from the well known result of robustness of radial sampling scheme to motion artifacts that are caused due to motion of the object as k-space lines are acquired.
The STCR reconstruction method here was applied to each coil data separately and the resulting images were combined in a square root of sum of squares fashion. In preliminary work, we have found small improvements by using joint multi-coil reconstruction techniques that use coil sensitivity information in the constrained reconstruction.
In conclusion, undersampled radial acquisition with a spatiotemporal constrained reconstruction method can be used to improve slice coverage of the heart for perfusion imaging. The radial STCR method produced images with better or comparable quality to standard Cartesian acquisitions and also with increased coverage. The method is relatively robust to respiratory motion in the data and results in images with high SNR and CNR. More work is warranted to test the method in patients with ischemia during vasodilation.
Grant support: R01EB0006155 and the Ben B. and Iris M. Margolis Foundation