Research Article  Open Access
Wei He, Yu Zhang, Junling Ding, Linman Zhao, "A Modified Phase Cycling Method for ComplexValued MRI Reconstruction", International Journal of Biomedical Imaging, vol. 2020, Article ID 8846220, 7 pages, 2020. https://doi.org/10.1155/2020/8846220
A Modified Phase Cycling Method for ComplexValued MRI Reconstruction
Abstract
The phase cycling method is a stateoftheart method to reconstruct complexvalued MR image. However, when it follows practical twodimensional (2D) subsampling Cartesian acquisition which is only enforcing random sampling in the phaseencoding direction, a number of artifacts in magnitude appear. A modified approach is proposed to remove these artifacts under practical MRI subsampling, by adding onedimensional total variation (TV) regularization into the phase cycling method to “preprocess” the magnitude component before its update. Furthermore, an operation used in SFISTA is employed to update the magnitude and phase images for better solutions. The results of the experiments show the ability of the proposed method to eliminate the ring artifacts and improve the magnitude reconstruction.
1. Introduction
In traditional magnetic resonance imaging (MRI), reconstructing an image requires the full Fourier domain (space) samples which are complex values and usually on a regular Cartesian grid. The full acquisition leads to an extremely timeconsuming scan process. Compressed sensing (CS) has proven to be effective to reconstruct highquality images from the incomplete set of space samples, accelerating the scan process [1–5]. And many researchers combine CS and parallel imaging techniques (CSPI) to further speed up the MRI acquisition [6–8]. Most approaches applying CS in MRI reconstruction [2–8] focus on the magnitude recovery, leaving the phase reconstruction less studied. However, MRI images are complexvalued, including magnitude and phase parts. And the phase of MRI signal also includes important information, such as field inhomogeneity or the velocity of blood flow [9], and can be utilized to assess field inhomogeneity and obtain clinically relevant physiological parameters [10]. Therefore, the phase map is of interest in many applications such as field map estimation [11] and phasecontrast imaging [12, 13]. If we can reconstruct accurate complexvalued MR image, i.e., its magnitude and phase maps simultaneously, from undersampled raw space data, it will make CS more generally applicable in MRI.
Zibetti and De Pierro [14] proposed a new penalty to separate the regularization of magnitude and phase, enforcing the magnitude to be sparse in finite differences domain and the phase to be smooth in firstorder roughness. This algorithm gives some better results than standard CS. But it is only applicable for firstorder differences operators in CS regularization, which is not the optimal choice. Moreover, in this method, the phase is still not reconstructed independently from the magnitude, because the parameter for the phase regularization term is constrained by the corresponding magnitude. It causes low SNR of the phase image in the locations having low magnitude. At last, the initialization for phase should be good when the firstorder phase difference is large, as the regularization term for phase is not convex here.
Some new regularization terms for phase component, exploiting the periodicity of phase, were introduced by Zhao et al. [15] to improve phase results. In this method, discontinuities are allowed to exist in phase, since the smoothness penalty for phase is modified. However, it has more computational complexity, e.g., about ten times, than the conventional CS.
The phase cycling method [16] proposed by Ong et al. is a stateoftheart algorithm for complexvalued MR image recovery. However, this approach follows the variable density Poissondisk plus partial Fourier (PDPF for short in this paper) subsampling scheme that is impractical for MRI hardware in 2D acquisition [17]. When random sampling is only imposed in the phaseencoding direction [18] (RSPe for short in this paper) on a 2D Cartesian grid, which is workable in practical MRI 2D acquisition, there are usually some noticeable artifacts in magnitude image, especially of brain data. The artifacts will be shown in Section 2.
In this paper, we utilize the smoothing property of total variation (TV) norm to remove the artifacts before each magnitude update under some sampling scheme and apply an update operation in the smoothingbased FISTA (SFISTA) [19] to improve the image recovery. The detailed procedure of this modified phase cycling method will be given in Section 2.
2. Materials and Methods
The cost function of the phase cycling method is where and are the magnitude and phase images, respectively, denotes the subsampled space measurements, is the system matrix for the complexvalued images, denotes the elementwise multiplication, is the elementwise exponential operator of the matrix , and and are the regularization functions for and , respectively. The first term in Equation (1) indicates the data fidelity, where is normally the product of the Fourier sampling operator and the sensitivity operator for multicoil data and reduces to a subsampling Fourier matrix for singlecoil data. For convenience, we call this data fidelity term .
The alternating minimization is performed to update the magnitude and phase images separately in this algorithm. And the proximal method [20, 21] is applied to solve each subproblem in the phase cycling method. The highlight of this method is that it shifts the phase wraps randomly over iterations to avoid significant error accumulating at the same location. It is because the phase wraps are artifacts in the initial solution and penalized by phase regularization in each iteration, resulting in errors to accumulate at the same position over iterations [16]. In this way, the magnitude and phase images are updated at iteration as follows: where is the proximal operator for function , and denote the stepsize for magnitude and phase , is a phase wrap picked randomly from the set of phase wraps generated by the initial solution with equal probability.
The phase cycling algorithm performs well based on PDPF undersampling scheme on a 2D Cartesian grid. Nevertheless, when it follows a practical RSPe subsampling scheme in 2D Cartesian acquisition, some artifacts pointed by the red arrows sometimes appear in magnitude image as shown in Figure 1. In addition, as the center of space data consists of low frequency signals, the central area is fullysampled in both undersampling schemes to improve the reconstruction quality and compute coil sensitivities for multicoil data.
Under the RSPe sampling pattern, the collected space data lack the whole echoes along several phaseencoding lines. Via inverse Fourier transform, the data on these missing lines are zeropadded, generating significant artifacts in the initial step in both magnitude and phase maps (see Figure 2). The artifacts are like arcs generally in the horizontal direction orthogonal to the phaseencoding coordinates. So far, our observation is that, under the RSPe sampling pattern, artifacts tend to happen in the initial step when the magnitude image has some very big values forming lines, arcs, or circles, which is common within brain images.
Since the phase cycling approach effectively averages the artifacts of phase wraps spatially in phase image over iterations by virtue of shifting the phase wraps randomly, the artifacts caused by the RSPe sampling pattern can be averaged incidentally. This can be deduced by comparing the phase maps under the RSPe sampling pattern in Figures 1 and 2. As phase values interact with magnitude values via calculating the gradient of the data fidelity term (see Equations (2) and (4)), the artifacts of magnitude reduce more or less eventually by means of delivering the abovementioned phase image with less or no artifacts to the data fidelity term. However, these errors in magnitude cannot be removed thoroughly in most situations.
In our method, before each update, onedimensional TV regularization is conducted on the magnitude images to eliminate the artifacts mentioned above. This is inspired by the piecewisesmooth nature of TV norm [22–24]. Since the artifacts appear to be oriented towards the coordinates orthogonal to the phaseencoding direction, the onedimensional TV regularization of magnitude shown in Equations (7) and (8) is enforced in the phaseencoding direction before each magnitude update in the phase cycling algorithm. where denotes the onedimensional TV norm in the phase encoding direction, denotes the magnitude after the onedimensional TV regularization for iteration , is the stepsize of the TV regularization, and is a positive regularization parameter.
Furthermore, inspired by SFISTA [19], the following generalized step called the smoothingbased proximal operator (SPO) in this paper replaces the proximal operator in the phase cycling method: where is the stepsize, denotes smooth approximation parameter, and represents the gradient of the data fidelity term with respect to the update variable .
The proposed method benefits from the smoothing effect of TV norm and the regularizationsolving capacity of SPO.
The pseudocode of the modified phase cycling method is shown as follows:

3. Results and Discussion
The proposed approach and the phase cycling method are implemented in MATLAB (MathWorks, Natick, MA) and run on a laptop with a 1.1 GHz Intel Core i7 with 6 multicores, and 16 GB memory. The stepsize for magnitude update is , and the stepsize for phase update is , where denotes the maximum eigenvalue calculated analytically. The maximum number of iterations is 100. The norm on the Daubechies 4 wavelet transform coefficients of magnitude is selected to be function unless otherwise specified. The smooth approximation parameter is set as 1 for the best performance of SPO. For each method, the regularization parameters are chosen for its optimal performance in terms of peaktosignal noise ratio (PSNR) and structural similarity index metric (SSIM) [25].
A fully sampled human brain dataset [16], with a 2D image size of and acquired by an 8channel head coil, is available in the software package at https://github.com/mikgroup/phase_cycling.git. regularization is imposed on the Daubechies 4 wavelet domain for the phase image. The positive TV regularization parameter is set to be 0.005. To get the measurements , data along 20% phaseencoding center lines are fully sampled and the other data are from 10% randomly chosen phaseencoding lines in the rest space. This undersampling pattern is shown in Figure 3.
Figure 4 shows the results of both methods and the proposed method without TV regularization on the brain dataset. To better display phase results, the background areas of the phase images were roughly masked out by thresholding the bottom 10% of the magnitude images. For the phase cycling algorithm, ring artifacts can be seen in the magnitude image. While, due to SPO, the proposed method without TV regularization reconstructs the magnitude image with less error compared to the phase cycling method, there are still some ring artifacts in the magnitude result. Furthermore, the proposed method removes the ring artifacts and promotes the quality of magnitude recovery significantly owe to the TV regularization and SPO. Comparing the magnitude images and their corresponding error maps of the proposed method with and without 1D TV regularization, it can be deduced that the 1D TV regularization can indeed decrease and even eliminate the ring artifacts. On the other hand, the magnitude image of the phase cycling method has more errors in general, which can be demonstrated by the contrast of the error maps. The three approaches perform comparably with each other in phase reconstruction.
Table 1 lists the PSNRs and SSIMs of the magnitude and phase maps by the proposed and phase cycling methods under the RSPe sampling pattern with a 30% sampling rate. Here, ten distinct RSPe sampling schemes with 30% samples are adopted to verify the robustness of the proposed method. The PSNR of the magnitude images by the proposed method improves about 1~3 dB, while the SSIM promotes about 0.01~0.02.

Another singlecoil brain dataset (see Dataset S1 in the Supplementary Material) is used to test the proposed algorithm. regularization is imposed on the Daubechies 6 wavelet domain for the phase image. . The data in the readout direction are always full sampled. And the random sampling only takes place in phaseencoding direction, with 20% center coordinates and 10% randomly selected coordinates elsewhere. Under this sampling pattern shown in Figure 5, Figure 6 displays the results of the proposed method with and without TV regularization and the phase cycling algorithm. The last method produces several artifacts pointed by arrows in the magnitude result. The artifact pointed by the yellow arrow may be mistaken as tissue in the brain. The proposed algorithm removes these artifacts distinctly thanks to the TV regularization, while without it the proposed method generates some ring artifacts. And the phase images reconstructed by three approaches are nearly the same.
Here, we still do ten experiments under ten different RSPe sampling schemes with 30% samples to calculate the PSNRs and SSIMs of magnitude images by the proposed and phase cycling methods. In terms of PSNR, the former results in , and the latter results in . And in terms of SSIM, the former results in , and the latter results in .
The value of the TV regularization parameter should be chosen cautiously. When it is too large, the magnitude image will blur and lose much detailed information, whereas too small renders this artifactelimination regularization ineffective. In our experience, the value between 0.001 and 0.005 is appropriate.
4. Conclusions
We add the onedimensional TV regularization of magnitude into the phase cycling algorithm, and use SPO to solve the magnitude and phase subproblems. The 1D TV regularization is like a preprocessing step for magnitude in each iteration. And the SPO improves the magnitude recovery to some extent. The modified algorithm can significantly reduce or totally eliminate the artifacts of magnitude by the phasecycling method under practical RSPe undersampling pattern, while the quality of the reconstructed phase image is comparable to the phase cycling method. However, the performance of the proposed algorithm is sensitive to the value of TV regularization parameter.
Data Availability
The brain data acquired by an 8channel head coil are from previously reported studies and datasets, which have been cited. The processed data are available in the software package at https://github.com/mikgroup/phase_cycling.git. The brain data used to support the findings of this study are included within the supplementary information file “Dataset S1.mat”.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This research was financially supported in part by the National Natural Science Foundation of China (61601396 and 31872704), in part by the Key Scientific Research Project of Colleges and Universities in Henan Province of China (19A520035).
Supplementary Materials
The brain data are acquired on 1.5 T GE Signa scanner (GE Healthcare, Waukesha, WI) with GRE sequence, supplied in a .mat file named “Dataset S1”. (Supplementary Materials)
References
 G. Adluru and E. V. R. DiBella, “Reordering for improved constrained reconstruction from undersampled kspace data,” International Journal of Biomedical Imaging, vol. 2008, Article ID 341684, 12 pages, 2008. View at: Publisher Site  Google Scholar
 M. Lustig, D. L. Donoho, and J. M. Pauly, “Sparse MRI: the application of compressed sensing for rapid MR imaging,” Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1182–1195, 2007. View at: Publisher Site  Google Scholar
 J. Huang, S. Zhang, and D. Metaxas, “Efficient MR image reconstruction for compressed MR imaging,” Medical Image Analysis, vol. 15, no. 5, pp. 670–679, 2011. View at: Publisher Site  Google Scholar
 Junfeng Yang, Yin Zhang, and Wotao Yin, “A fast alternating direction method for TVL1L2 signal reconstruction from partial Fourier data,” IEEE Journal of Selected Topics in Signal Processing, vol. 4, no. 2, pp. 288–297, 2010. View at: Publisher Site  Google Scholar
 S.S. Wang, J.B. Liu, Q.G. Liu et al., “Iterative feature refinement for accurate undersampled MR image reconstruction,” Physics in Medicine & Biology, vol. 61, no. 9, pp. 3291–3316, 2016. View at: Publisher Site  Google Scholar
 R. Otazo, D. Kim, L. Axel, and D. K. Sodickson, “Combination of compressed sensing and parallel imaging for highly accelerated firstpass cardiac perfusion MRI,” Magnetic Resonance in Medicine, vol. 64, no. 3, pp. 767–776, 2010. View at: Publisher Site  Google Scholar
 H. Chandarana, L. Feng, T. K. Block et al., “Freebreathing contrastenhanced multiphase MRI of the liver using a combination of compressed sensing, parallel imaging, and goldenangle radial sampling,” Investigative Radiology, vol. 48, no. 1, pp. 10–16, 2013. View at: Publisher Site  Google Scholar
 S. Wang, S. Tan, Y. Gao et al., “Learning jointsparse codes for calibrationfree parallel MR imaging,” IEEE Transactions on Medical Imaging, vol. 37, no. 1, pp. 251–261, 2018. View at: Publisher Site  Google Scholar
 S. Chavez, QingSan Xiang, and L. An, “Understanding phase maps in MRI: a new cutline phase unwrapping method,” IEEE Transactions on Medical Imaging, vol. 21, no. 8, pp. 966–977, 2002. View at: Publisher Site  Google Scholar
 W. Liu, X. Tang, Y. Ma, and J. Gao, “3D phase unwrapping using global expected phase as a reference: application to MRI global shimming,” Magnetic Resonance in Medicine, vol. 70, no. 1, pp. 160–168, 2013. View at: Publisher Site  Google Scholar
 A. K. Funai, J. A. Fessler, D. T. B. Yeo, V. T. Olafsson, and D. C. Noll, “Regularized field map estimation in MRI,” IEEE Transactions on Medical Imaging, vol. 27, no. 10, pp. 1484–1494, 2008. View at: Publisher Site  Google Scholar
 J. De Poorter, C. De Wagter, Y. De Deene, C. Thomsen, F. Ståhlberg, and E. Achten, “Noninvasive MRI thermometry with the proton resonance frequency (PRF) method: in vivo results in human muscle,” Magnetic Resonance in Medicine, vol. 33, no. 1, pp. 74–81, 1995. View at: Publisher Site  Google Scholar
 J. Nielsen and K. S. Nayak, “Referenceless phase velocity mapping using balanced SSFP,” Magnetic Resonance in Medicine, vol. 61, no. 5, pp. 1096–1102, 2009. View at: Publisher Site  Google Scholar
 M. V. W. Zibetti and A. R. De Pierro, “Separate magnitude and phase regularization in MRI with incomplete data: preliminary results,” in Proceeding of IEEE International Symposium on Biomedical Imaging: From Nano to Macro, pp. 736–739, Rotterdam, Netherlands, 1417 April 2010. View at: Google Scholar
 Feng Zhao, D. C. Noll, J. F. Nielsen, and J. A. Fessler, “Separate magnitude and phase regularization via compressed sensing,” IEEE Transactions on Medical Imaging, vol. 31, no. 9, pp. 1713–1723, 2012. View at: Publisher Site  Google Scholar
 F. Ong, J. Y. Cheng, and M. Lustig, “General phase regularized reconstruction using phase cycling,” Magnetic Resonance in Medicine, vol. 80, no. 1, pp. 112–125, 2018. View at: Publisher Site  Google Scholar
 M. Lustig, J. H. Lee, D. L. Donoho, and J. M. Pauly, “Faster imaging with randomly perturbed, undersampled spirals and L_1 reconstruction,” in Proceeding of the 13th Annual Meeting of ISMRM, Miami Beach, Florida, USA, 0713 May 2005. View at: Google Scholar
 M. V. W. Zibetti and A. R. De Pierro, “Improving compressive sensing in MRI with separate magnitude and phase priors,” Multidimensional Systems and Signal Processing, vol. 28, no. 4, pp. 1109–1131, 2017. View at: Publisher Site  Google Scholar
 Z. Tan, Y. C. Eldar, A. Beck, and A. Nehorai, “Smoothing and decomposition for analysis sparse recovery,” IEEE Transactions on Signal Processing, vol. 62, no. 7, pp. 1762–1774, 2014. View at: Publisher Site  Google Scholar
 P. L. Combettes and J. C. Pesquet, Proximal splitting methods in signal processing, Springer Optimization and Its Applications, Springer, New York, 2011.
 N. Parikh and S. Boyd, “Proximal algorithms,” Foundations and Trends in Optimization, vol. 1, no. 3, pp. 127–239, 2014. View at: Publisher Site  Google Scholar
 I. Daubechies, G. Teschke, and L. Vese, “Iteratively solving linear inverse problems under general convex constraints,” Inverse Problems and Imaging, vol. 1, no. 1, pp. 29–46, 2007. View at: Publisher Site  Google Scholar
 D. Strong and T. Chan, “Edgepreserving and scaledependent properties of total variation regularization,” Inverse Problems, vol. 19, no. 6, pp. 165–187, 2003. View at: Publisher Site  Google Scholar
 A. Beck and M. Teboulle, “Fast gradientbased algorithms for constrained total variation image denoising and deblurring problems,” IEEE Transactions on Image Processing, vol. 18, no. 11, pp. 2419–2434, 2009. View at: Publisher Site  Google Scholar
 Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600–612, 2004. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2020 Wei He et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.