Journal of Radiology and Imaging
An International PeerReviewed Open Access Journal
ISSN 23998172
 Download PDF
 
 Download Citation
 
 Email a Colleague
 
 Share:

 Tweet

Journal of Radiology and Imaging
Volume 4, Issue 5, August 2020, Pages 30–39
Original researchOpen Access
Iterative versus noniterative image reconstruction methods for sparse magnetic resonance imaging
 ^{1 }Department of Radiology and Imaging Sciences, University of Utah, 729 Arapeen Drive, Salt Lake City, Utah, 84108, USA
 ^{2 }Department of Computer Science, Utah Valley University, 800 West University Parkway, Orem, UT 84058, USA
*Corresponding author: Gengsheng L. Zeng, Department of Radiology and Imaging Sciences, University of Utah, 729 Arapeen Drive, Salt Lake City, Utah, 84108, USA. Tel.: +1 (801) 5813918, +1 (801) 8634667; Fax: +1 (817) 5823918; Email: larry.zeng@hsc.utah.edu; larry.zeng@uvu.edu; larryzeng@live.com
Received 7 June 2020 Revised 25 July 2020 Accepted 30 July 2020 Published 3 August 2020
DOI: http://dx.doi.org/10.14312/23998172.20205
Copyright: © 2020 Zeng GL, et al. Published by NobleResearch Publishers. This is an openaccess article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution and reproduction in any medium, provided the original author and source are credited.
AbstractTop
Magnetic resonance imaging (MRI) using undersampled kspace data is a common method to shorten the imaging time. Iterative Bayesian algorithms are usually used for its image reconstruction. This paper compares an iterative Bayesian image reconstruction method that uses both spatial and temporal constraints and a noniterative reconstruction algorithm that does not use temporal constraints. Three patient studies are performed. It is interesting to notice that the images reconstructed by the iterative Bayesian algorithm may introduce more bias than the noniterative algorithm, even though the images provided by the iterative Bayesian algorithm look less noisy. The bias can be reduced by decreasing the influence of the temporal constraints.
Keywords: tomographic image reconstruction; undersampled measurements; fast MRI; analytics reconstruction
IntroductionTop
Realtime or dynamic imaging in magnetic resonance imaging (MRI) require fast data acquisition, and the kspace data acquired is incomplete [1, 2]. When the data is incomplete, the conventional analytical image reconstruction methods usually produce severe aliasing artifacts. Using undersampled data to reconstruct an image is sometimes referred to as a compressed sensing problem, which is commonly solved by constrained optimization. Iterative algorithms with constraints are able to produce images with less severe artifacts and with higher signaltonoise ratios [38]. The iterative methods in image reconstruction commonly have an objective function that contains spatial and/or temporal constraints. For example, if the image is sparse, the sparseness can be enforced by using the L_{0} norm or L_{1} norm constraints. If the image is piecewise constant, the total variation (TV) norm constraint can be used. In MRI, if the kspace is scanned using a nonCartesian sampling scheme, it is popular to first interpolate the nonCartesian kspace measurements into the Cartesian grid before reconstruction [911].
Nowadays deep learning methods are popular in almost all areas of this modern world and they are well accepted in MRI image reconstruction and denoising [1215]. Deep learning methods extract features from the training data, and their performance is training data dependent. It is not clear whether the selection of the training data affects the image bias, if the kspace data is incomplete and the image is reconstructed by a neural network.
Parallel imaging using multiple receiver coils in MRI is another way to achieve fast data acquisition [1618]. Each coil can measure a different region in the kspace. The combined data from all coils is equivalent to a complete data set. This paper’s focus is on incomplete data image reconstruction. Therefore, the patient studies in this paper only use the (incomplete) sparsedata obtained from only one coil.
For image reconstruction with an incomplete data set, some a priori information is enforced to supplement the missing data. As a priori information, for example, one could assume that the image changes slowly in time and the image is piecewise constant. Bayesian image reconstruction methods in general present images with higher signaltonoise ratios than the methods that do not use prior information.
There is one concern for all Bayesian methods. What happens if the a priori information is wrong? The Bayesian methods may thus cause bias. This paper compares an iterative Bayesian algorithm [19] and a noniterative nonBayesian algorithm [20] in terms of image bias using some MRI patient studies.
MethodsTop
Iterative algorithm with spatiotemporal constraints
In this paper, we consider sparse radial kspace sampling in MRI, in which the number of views is not sufficient. The sparseview data may cause streaking aliasing artifacts in the reconstructed images. Two algorithms are compared in this paper in terms of their performance for sparseview MRI.
The iterative Bayesian algorithm used in the comparison studies here was taken from reference [19] and is briefly outlined as follows. The aim of the iterative algorithm is to minimize an objective function that consists of one data fidelity term and two Bayesian terms. The objective function, C_{t}, can be expressed as
where m is the spatial domain image, F is the twodimensional (2D) Fourier transform, W is the masking matrix depending on the kspace sampling scheme, N is the total number of image pixels, and ε is a small positive constant. The mathematical symbols ∇_{x}, ∇_{y}, ∇_{t} are the gradients in the x, y, the time (t) directions, respectively. The image is represented in the xy Cartesian coordinate system. In the righthandside of Eq. (1), the first two terms use the L_{2} norm, while the third term uses the total variation (TV) norm. The control parameters α_{1} and α_{2} determine the balance among these three norms in the objective function (or cost function). In our implementation, α_{1} = 0.04, α_{1} = 0.006, ε = 0.00000001; these parameters were suggested in reference [19]. An iterative gradient descent algorithm is used to optimize the objective function, C_{t}, as defined in Eq. (1). The number of iterations is chosen to be 1000 as suggested in reference [19].
In this iterative algorithm, the sparse radial kspace measurements are first interpolated into a Cartesian grid, d. For a 2D image array m, d is a 2D complex matrix in the kspace. Since the kspace measurements are sparse, the majority of the elements in d is zero. Matrix W in Eq. (1) is a binary masking 2D matrix, with the same size as matrix d. The element of W is 1 if the corresponding element in d is measured. The element of W is 0 if the corresponding element in d is not measured. In Eq. (1), Fm denotes the 2D Fourier transform of the 2D image m. The notation WFm represents the elementbyelement multiplication between the mask matrix W and the 2D Fourier transform, Fm, of the image m.
The image m reconstructed by this iterative algorithm is complex, which consists of a real part mreal and an imaginary part mimaginary. The last step in the image reconstruction algorithm is the convert the two 2D image arrays to the norm image on the elementbyelement base. In other words, the final image is obtained as
Noniterative algorithms without temporal constraint
The noniterative algorithm used in our comparison studies is taken from reference [20] with a modification. This noniterative algorithm does not use any constraints in the temporal direction. In other words, the images at different time frames are assumed to be independent. The main idea of this noniterative algorithm is first to extend the sparseview kspace measurements to more synthetic views in order to reduce the undersampling aliasing artifacts and then to reconstruct the image with an analytical filtered backprojection (FBP) algorithm. The FBP algorithm is the most popular image reconstruction algorithm in computed tomography (CT).
Since the kspace data is complex, we process the real part and the imaginary part separately and obtain two images: the real part and the imaginary part. The final image is the norm image by combining the realpart image and the imaginaryimage as indicated by Eq. (2). This is only aspect that this algorithm differs from the one developed in reference [20], where the complex kspace data was first converted into norm data by combining the real part and imaginary part. We believe that it is more proper to construct the realpart image and imaginarypart image first and then to form the norm image by combining these two images as the final step. The current method is mathematically equivalent to use the complex measurements directly for image reconstruction and convert the complex image to a norm image at the end.
The outline of the noniterative algorithm is described as follows. Instead of interpolating the radial kspace measurements into a Cartesian grid, data on some unmeasured radial kspace lines is estimated by the deformation method. The deformation method is illustrated in Figure 1, where the two solid lines represent the locations where the measurements are available, and the two broken lines represent the locations where the measurements are not available and to be estimated. The values on line L1 are deformed to obtain the values on line L4, as indicated by the blue arrows in Figure 1, in the sense that for every value on L1 there is a corresponding close value on L4. The red dot on L1 and the red dot on L4 illustrate the corresponding values. Those blue arrows in Figure 1 form a deformation field. The name ‘deformation’ implies that L4 can be thought of being deformed from L1. The deformation method is to assign the ‘red dot’ value to lines L2 and L3 at the intersections with the deformation vector (i.e., the blue arrow). The deformation field (i.e., the blue arrows) can be estimated using a noniterative algorithm. For more details of this deformation method, please consult reference [20].
In this paper, the sparse MRI data contains 24 kspace radial lines. The extended MRI data contains 72 kspace lines. Between every pair of measured lines, measurements on two new lines are estimated, as shown in Figure 1.
After the extended radial measurements are obtained, an analytic FBP is carried out as following steps for real part and imaginary part of the data separately:
Step1: Linebyline filter the kspace measurements by a 1D transfer function
where is the absolute value of the frequency, that is, the distance to the center of the kspace. The parameter β determines the level of regularization, as suggested by reference [21].
Step 2: Linebyline take the 1D inverse Fourier transform.
Step 3: Perform backprojection.
Two images are obtained by application of the above three steps to the realpart and imaginarypart of the measurements. The final image is the norm image by combining these two images according to Eq. (2).
The transfer function (3) is a regularized version of the original ramp filter
as suggested by reference [21], aiming to find a minimum norm solution. This noniterative algorithm can have multiple versions. Some of the versions are used in this paper for comparison studies:
Version 1: As described above, the input data for ‘Step 1’ is the extended 72view measurements. The extended 72view measurements are obtained by the deformation method from the 24view raw measurements.
Version 2: The input data for ‘Step 1’ is the 24view raw measurements. The three steps are the same as presented above.
Version 3: This version is almost the same as Version 2, except that the filter transfer function for ‘Step 1’ is the ramp filter given by Eq. (4), instead of the one in Eq. (3).
Version 4: This version is almost the same as Version 3, except that the input data for ‘Step 1’ is the 72view ‘complete’ raw measurements. Note that these 72vew measurements are directly acquired from the MRI scanner; they are not extended measurements from the 24view raw measurements. The final image from Version 4 is treated as the gold standard for the results from other methods to compare.
MRI patient data
Actual patient cardiac perfusion MRI data is used in this paper for comparison studies. All studies were approved by our Institutional Review Board (IRB). All participating patients were properly consented. Patient identifiable information was removed before transferring to our research computers for processing.
Data acquisition was performed using a Siemens 3T Trio scanner. A phased array of multiple of coils was used during data acquisition; however, data from only one selected coil was used in image reconstruction. The scanner parameters for the radial acquisition were TR = 2.6 ms, TE = 1.1 ms, flip angle = 12^{o}, Gd dose = 0.03 mmol/kg, and slice thickness = 6 mm. Reconstruction pixel size was 1.8 × 1.8 mm^{2}. Each image was acquired in a 62 ms readout. The acquisition matrix size for an image frame was 256 × 72, and 75 sequential images were obtained at 75 different times. At each time frame, the kspace is sampled with 72 uniformly spaced radial lines over an angular range of 180^{o}.
For the comparison purposes of this paper, at each time frame we uniformly undersampled the 72 views into 24 views and call such obtained data as 24view raw data. The FBP reconstructed images with 72 views were treated as the gold standard. The image discrepancy between the gold standard and other images are evaluated in terms of rootmeansquareerror (RMSE). A median filter is applied to all images before the RMSEs were measured and displayed.
ResultsTop
Results are shown in Figures 2, 3 and 4, for patient 1, patient 2 and patient 3, respectively. For each patient, there are 5 sets of reconstructed images by the iterative Bayesian reconstruction algorithm and by the 3 versions of the noniterative reconstruction methodology. Each set consists of 75 sequential images. There would be too many images to show if all images were displayed. To make the paper more readable, 4 (instead of 75) images are shown for each set. Table 1 lists the rootmeansquareerrors (RMSEs) of images deviate from the gold standard.
It is clear from Table 1 to observe that directly using the 24view raw data gives the worst performance. The iterative algorithm and the noniterative algorithm have similar performance, with the noniterative algorithm’s results being slightly better.
72view noniterative (gold standard) β = 0 
24view noniterative (raw reconstruction) β = 0 
24view → 72view, noniterative, without temporal constraint β = 1 
24view, iterative Bayesian, with temporal constraint α_{1} = 0.04 

Patient 1  0 
6.7675 
3.6778 
4.1719 
Patient 2  0 
10.957 
5.0849 
5.3663 
Patient 3  0 
10.747 
4.9476 
5.3545 
All these algorithms have some parameters to adjust. It is interesting to observe that when the parameter α_{1} in the objective function (1) is reduced from 0.04 to 0.004, the RMSE is reduced. We conjecture that the RMSE for the iterative algorithm is caused by the bias. Reducing the influence of the temporal constraint in the objective function (1) reduces the bias. A Bayesian constraint can help the image look better, but at the same time, the Bayesian constraint can change the image value.
For the noniterative algorithm, a regularization factor β is introduced to encourage a minimum norm solution. When β = 0, there is no regularization. When β is large, the regularization is heavily enforced. In Figures 5∼7 and in Table 2, some other cases are shown. One case is the iterative reconstruction with a smaller α_{1} = 0.004. Another case is the noniterative reconstruction with heavier regularization using β = 2. The third case is the noniterative reconstruction with β = 1 and with 24view raw data, where the 24 views are not extended to 72 views.
From Figures 5∼7 and Table 2, we observe the following: For the iterative algorithm, reducing the temporal constraint helps reducing the bias. Increasing the regularization in the noniterative algorithm helps to reduce the noise. If the 24view data is not extended to 72data, the regularization helps denoising, but extension to 72view data can further reduce some artifacts.
72view noniterative (gold standard) β = 0 
24view noniterative (raw reconstruction) β = 1 
24view → 72view, noniterative, without temporal constraint β = 2 
24view, iterative Bayesian, with temporal constraint α_{1} = 0.04 

Patient 1  0 
4.2129 
3.4196 
3.9046 
Patient 2  0 
6.4403 
4.6963 
4.7018 
Patient 3  0 
6.2731 
4.5616 
4.6803 
The structural similarity index measure (SSIM) is another well accepted figureofmerit to find the similarity of two images; the definition details of SSIM are given in [22]. The counterparts of Tables 1 and 2 are Tables 3 and 4, respectively. The results of Tables 3 and 4 are consistent with those in Tables 1 and 2, indicating the superiority of the proposed noniterative algorithm. The optimal value for RMSE is 0, and the optimal value for SSIM is 1.
72view noniterative (gold standard) β = 0 
24view noniterative (raw reconstruction) β = 0 
24view → 72view, noniterative, without temporal constraint β = 1 
24view, iterative Bayesian, with temporal constraint α_{1} = 0.04 

Patient 1  1 
0.999635875413527 
0.999848361889880 
0.999807461269317 
Patient 2  1 
0.999413359233490 
0.999804906454934 
0.999768375404101 
Patient 3  1 
0.999421484562476 
0.999814647345931 
0.999769234887647 
72view noniterative (gold standard) β = 0 
24view noniterative (raw reconstruction) β = 1 
24view → 72view, noniterative, without temporal constraint β = 2 
24view, iterative Bayesian, with temporal constraint α_{1} = 0.04 

Patient 1  1 
0.999800444659817 
0.999849666451081 
0.999819971016024 
Patient 2  1 
0.999697329429499 
0.999803495341621 
0.999799850215280 
Patient 3  1 
0.999705810078143 
0.999813456704611 
0.999801863895034 
DiscussionTop
The iterative algorithm uses the TV norm for regularization in the image spatial domain. The TV norm also in favor of piecewise constant images and may cause staircase artifacts when encouraging sharp edges. On the other hand, the noniterative algorithm uses a lowpass filter regularize the image by controlling the noise. This lowpass filter was derived by encouraging a minimum norm solution [21]. We observe from Table 2 that by increasing β value from 1 to 2 the noise is reduced and the RMSE is also reduced. However, the SSIM is reduced for one patient and increases for the other two patients. Therefore, using a lowpass filter alone is not effective in reducing the aliasing artifacts caused by undersampled data. The noniterative deformation method is used to estimate some unmeasured data. Comparing Table 1 with Table 2 and Table 3 with Table 4, the extension of the 24view data to 72view data actually improves both the RMSE and SSIM. In other worlds, Version 1 is better than Version 2.
We must point out that the extended 72view data is different from the actually measured 72view data and contains estimation errors. Of course, there is no substitute for the real data. If the real data is not available, the estimated data is the second choice. There are numerous ways to estimate the unmeasured data. The iterative algorithm in this paper estimates the missing data from the measured data in different time frames adjacent to the current time, while the noniterative algorithm in this paper estimates the missing data from the measured data at the current time.
ConclusionTop
In dynamic MRI imaging, timeactivity curves are used to estimate tissue kinetic parameters, which have wide applications in clinical diagnoses. Accurate kinetic parameters depend on the unbiased timeactivity curves that are extracted from dynamic images. The dynamic images are reconstructed using sparse kspace data.
Two sparsedata image reconstruction methods are compared into this paper: an iterative Bayesian method that uses a temporal constraint and a noniterative method that does not have a temporal constraint. According to the rootmeansquareerror (RMSE) analysis and structural similarity index measure (SSIM) analysis, these two methods give similar accuracies as shown in Tables 1∼4. The noniterative method is much faster than the iterative method. In the iterative reconstruction algorithm, the 24view sparsedata is intrinsically extended to 72view data by using the temporal constraint as shown in the objective function (1). The ‘assisting’ 48view data is most likely different form the actual measured data. The influence of the ‘assisting’ 48view data may bring in image bias. This bias is the price we pay to reduce the data undersampling artifacts. After reducing the weighting parameter α_{1} from 0.04 to 0.004 in the iterative algorithm, the bias is reduced and the RMSE and SSIM are improved.
Our takehome message is that Bayesian constraints may make the images look better, but they may cause more image bias, deviating from the true solution. We must be aware of the tradeoff between looking good and accuracy.
Our future research includes investigations of deep learning methods for sparsedata MRI, in which the unmeasured data is supplemented by training data of other patients. The deep learning methods are effective but are not well understood, in terms of how the a priori information is presented and used. We will focus on whether the deep learning methods may introduce image bias when the kspace is not fully sampled.
Acknowledgements
The American Heart Association Precision Medicine Platform (https://precision.heart.org/) was used for data analysis.
Funding
This research is partially supported by American Heart Association under award number 18AJML34280074.
Conflicts of interest
Authors declare no conflict of interest.
A list of abbreviations
FBP: Filtered backprojection; MRI: Magnetic Resonance Imaging; 1D: One dimensional; RMSE: Root mean square error; SVD: Singular value decomposition; TV: Total variation.
Authors' contributions
All authors read and approved the final manuscript.
ReferencesTop
[1]Chu C, Liu G, Janowski M, Bulte JWM, Li S, et al. Realtime MRI guidance for reproducible hyperosmolar opening of the bloodbrain barrier in mice. Front. Neurol. 2018; 9:921.Article Pubmed
[2]Mukherjee RK, Chubb H, Roujol S, Razavi R, O'Neill MD. Advances in realtime MRIguided electrophysiology. Curr Cardiovasc Imaging Rep. 2019; 12(2):6.Article Pubmed
[3]Donoho DL. For most large underdetermined systems of linear equations the minimal 1norm solution is also the sparsest solution. Commun. Pure Appl Math. 2004; 59(6):797–829.Article
[4]Donoho DL. Compressed sensing. IEEE Trans Inf Theory. 2006; 52(4):1289–306.Article
[5]Candès EJ, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans Inf Theory. 2006; 52(2):489–509.Article
[6]Lustig M, Donoho DL, Pauly JM. Sparse MRI: the application of compressed sensing for rapid MR imaging. Magn Reson Med. 2007; 58(6):1182–1195.Article Pubmed
[7]Chen SW, Chao SC. A reweighted ℓ1minimization based compressed sensing for the spectral estimation of heart rate variability using the unevenly sampled data. PLoS One, 2014; 9(6):e99098.Article Pubmed
[8]Weizman L, Eldar YC, Bashat DB. Compressed sensing for longitudinal MRI: An adaptiveweighted approach. Med Phys. 2015; 42(9):5195.Article Pubmed
[9]Lingala SG, Toutios A, Toger J, Lim Y, Zhu Y, et al. Stateoftheart MRI protocol for comprehensive assessment of vocal tract structure and function. Proceedings of the Annual Conference of Interspeech, San Francisco, CA, September 812, 2016; 475–479.Article
[10]Kholmovski EG, Coulombe N, Silvernagel J, Angel N, Parker D, et al. Realtime MRIguided cardiac cryoablation: A feasibility study. J Cardiovasc Electrophysiol. 2016; 27(5):602–608.Article Pubmed
[11]Tian Y, Mendes J, Pedgaonkar A, Ibrahim M, Jensen L, et al. Feasibility of multipleview myocardial perfusion MRI using radial simultaneous multislice acquisitions. PLoS One, 2019; 14(2):e0211738.Article Pubmed
[12]Gibbons EK, Hodgson KK, Chaudhari AS, Richards LG, Majersik JJ, et al. Simultaneous NODDI and GFA parameter map generation from subsampled qspace imaging using deep learning. Magn Reson Med. 2019; 81(4):2399–2411.Article Pubmed
[13]Hyun CM, Kim HP, Lee SM, Lee S, Seo JK. Deep learning for undersampled MRI reconstruction. Phys Med Biol. 2018; 63(13):135007.Article Pubmed
[14]Jeelani H, Martin J, Vasquez F, Salerno M, Weller DS. Image quality affects deep learning reconstruction of MRI. 2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018), Washington, DC, 2018; 357–360.Article
[15]Yang G, Yu S, Dong H, Slabaugh G, Dragotti PL, et al. DAGAN: Deep dealiasing generative adversarial networks for fast compressed sensing MRI reconstruction. IEEE Trans Med Imaging. 2018; 37(6):1310–1321.Article Pubmed
[16]Sodickson DK, Manning WJ. Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays. Magn Reson Med. 1997; 38(4):591–603.Article Pubmed
[17]Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: sensitivity encoding for fast MRI. Magn Reson Med. 1999; 42(5):952–962.Article Pubmed
[18]Larkman DJ, Nunes RG. Parallel magnetic resonance imaging. Phys Med Biol. 2007; 52(7):R15–R55.Article
[19]Adluru G, McGann C, Speier P, Kholmovski EG, Shaaban A, et al. Acquisition and reconstruction of undersampled radial data for myocardial perfusion MRI. J Magn Reson Imaging. 2009; 29(2):466–473.Article Pubmed
[20]Zeng GL, DiBella RV. Noniterative image reconstruction from sparse magnetic resonance imaging radial data without priors. Vis Comput Ind Biomed Art. 2020; 3(1):9.Article Pubmed
[21]Zeng GL. A filtered backprojection MAP algorithm with nonuniform sampling and noise modeling. Med Phys. 2012; 39(4):2170–2178.Article Pubmed
[22]Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. Image quality assessment: from error visibility to structural similarity. IEEE Trans Image Process. 2004; 13(4):600–612.Article Pubmed