Magnetic Resonance in Medicine 73:523–535 (2015)

PROMISE: Parallel-Imaging and Compressed-Sensing Reconstruction of Multicontrast Imaging Using SharablE Information Enhao Gong,1,2 Feng Huang,3* Kui Ying,2,4 Wenchuan Wu,2 Shi Wang,4 and Chun Yuan2,5 INTRODUCTION

Purpose: A typical clinical MR examination includes multiple scans to acquire images with different contrasts for complementary diagnostic information. The multicontrast scheme requires long scanning time. The combination of partially parallel imaging and compressed sensing (CS-PPI) has been used to reconstruct accelerated scans. However, there are several unsolved problems in existing methods. The target of this work is to improve existing CS-PPI methods for multicontrast imaging, especially for two-dimensional imaging. Theory and Methods: If the same field of view is scanned in multicontrast imaging, there is significant amount of sharable information. It is proposed in this study to use manifold sharable information among multicontrast images to enhance CSPPI in a sequential way. Coil sensitivity information and structure based adaptive regularization, which were extracted from previously reconstructed images, were applied to enhance the following reconstructions. The proposed method is called Parallel-imaging and compressed-sensing Reconstruction Of Multicontrast Imaging using SharablE information (PROMISE). Results: Using L1-SPIRiT as a CS-PPI example, results on multicontrast brain and carotid scans demonstrated that lower error level and better detail preservation can be achieved by exploiting manifold sharable information. Besides, the privilege of PROMISE still exists while there is interscan motion. Conclusion: Using the sharable information among multicontrast images can enhance CS-PPI with tolerance to motions. C 2014 Wiley PeriodiMagn Reson Med 73:523–535, 2015. V cals, Inc. Key words: parallel imaging; compressed sensing; multicontrast; multichannel; sharable information

1 Magnetic Resonance System Research Lab, Department of Electrical Engineering, Stanford University, Stanford, California, USA. 2 Center for Biomedical Imaging Research, Department of Biomedical Engineering, Tsinghua University, Beijing, China. 3 Philips Research China, Shanghai, China. 4 Key Laboratory of Particle and Radiation Imaging, Department of Engineering Physics, Tsinghua University, Beijing, China. 5 Department of Radiology, University of Washington, Seattle, Washington, USA.

*Correspondence to: Feng Huang, Ph. D., Philips Research China, No.1, Building 10, Lane 888, Tianlin Road, Shanghai, China, 200233. E-mail: [email protected] Additional Supporting Information may be found in the online version of this article. Received 28 March 2013; revised 29 December 2013; accepted 2 January 2014 DOI 10.1002/mrm.25142 Published online 25 February 2014 in Wiley Online Library (wileyonlinelibrary. com). C 2014 Wiley Periodicals, Inc. V

MRI is a powerful and widely applied imaging modality in both clinical and research settings. One of the most important features of MRI is that it can provide images with multiple contrasts for complementary diagnostic information. Therefore, a typical clinical MR examination is composed of several scans to acquire images with different contrasts, such as T1 weighted (T1w), T2 weighted (T2w), and Proton density weight (PDw) images. The multicontrast scheme provides more information for diagnosis than single contrast imaging with longer acquisition, which may result in higher cost and higher sensitivity to motion artifacts. Therefore, reducing total acquisition time of multicontrast imaging is highly desired. Partially parallel imaging (PPI) (1,2), which exploits the data correlation among multiple coil elements for reconstruction with partially acquired data, has been widely used clinically for fast imaging. However, PPI techniques are limited by channel encoding capability and it reduces acquisition time at the cost of an increase of noise and artifact level. Compressed sensing (CS) (3) recovers undersampled data based on the assumption of image sparsity in certain transformed domain(s). Because these two categories of methods use different auxiliary information for reconstruction, the combination of these two methods, called CS-PPI in this work, can result in impressively better recovery than each individual method which have been shown in literatures (4–8). However, there are still several problems unsolved in CS-PPI. The most difficult problem is how to get accurate coil sensitivity information, which is essential for the accuracy of parallel imaging based reconstruction. The calculation of coil sensitivity uses either low resolution prescan (1) or self-calibration signal (2). Coil sensitivity estimation from prescan may suffer from prolonged acquisition time or misregistration due to patient motion. The self-calibration scheme needs sufficient autocalibration signal (ACS) for accurate coil sensitivity while the acquisition of ACS limits the net acceleration factor. Another problem is how to balance artifact-reduction and detail preservation using the predefined regularization parameter. An improper regularization may result in either residual noise/artifact or over-smoothed image. The scans for multiple contrasts share significant amount of redundant information, due to the same field of view (FOV) is imaged in the same system using the same radiofrequency coil. In this study, we propose to use manifold sharable information to tackle the existing

523

524

Gong et al.

problems and hence improve CS-PPI reconstruction for multicontrast imaging. In literatures, sharable information among multicontrast images has been used to improve either PPI or CS. In 2006, Larkman et al (9) proposed to use joint histogram entropy among images with different contrasts to enhance PPI. Recently, Li et al (10–12) proposed a method referred as “Correlation Imaging” that used the sharable sensitivity based information, called “correlation function,” to reconstruct the multicontrast images in the context of PPI. Bilgic et al (13) uses a Bayesian based method to exploit the similarity of sparsity among images with different contrasts in the context of compressed sensing (3). These previous methods do not fully take the advantages of the sharable information. Both Larkman’s method for PPI and Bilgic’s method for CS only uses the sharable anatomical information. Li’s method on correlation imaging only uses the similarity of coil sensitivity information. These previous methods also have some potential issues. Mutual information and Bayesian based method take high computational costs and Correlation Imaging mainly works with uniform sampling. Another potential problem of these existing techniques is the sensitivity to interscan motions. If there are interscan motions, the structural information will be mismatched between images with different contrasts. Strong constraints in the reconstruction model, on the similarity of structure [like in Larkman et al (9) and Bilgic et al (13)] or the coil sensitivity maps in the image domain (as in Li and Dumoulin) (10), could result in artifacts, which may mislead the diagnosis. Different from previous methods, in this study, it is proposed to use manifold sharable information to improve CS-PPI. Also, all constraints using the sharable information are more heuristic and softer than the existing methods to avoid side effects due to motions. Same as previous methods, this study also focus on multicontrast two-dimensional (2D) imaging because multicontrast 2D imaging is more intensively used clinically than 3D imaging (10) and can be more time consuming, The proposed method is called Parallel-imaging and compressed-sensing Reconstruction of Multicontrast Imaging using SharablE information (PROMISE). The rest of this study is organized as follows. In Theory section, how to use the sharable information to solve the existing problems in CS-PPI is introduced. The implementation details and the design of experiments are described in the Methods section. The validation of the method is presented in the Result section. In the Discussion section, the advantages of PROMISE and further generalization of PROMISE are included. Finally, we provide concluding remarks in the Conclusion section. THEORY

several properties of the scanned subject according to Boltzmann distribution (14). B 1 is the sensitivity map of the receive coil. F ð•Þ denotes a function of variables in parentheses, T1 and T2 are the T1 and T2 maps of the scanned subject, TR is the repetition time, and TE is the echo time. For T1 weighted, T2 weighted, and proton density weighted images, only TR and TE are different in acquisition. If the same location is scanned, these images share the same structural information in r, the same magnetic information B0 , B 1 , as well as T1 and T2 maps. In this work, it is proposed to take advantages of the set of sharable information, including coil sensitivity information and image structural information, to improve CS-PPI. The list above demonstrates a practical example of the set of sharable information in multicontrast fast imaging. For different applications, it is possible to use additional or alternative sets. Combination of PPI and CS The combination of PPI and CS is the trend of fast imaging because it can produce better images than those by each individual method (4–7,15). In this study, L1-SPIRiT (7), an iterative reconstruction method using L1 norm and k-space data self-consistency as regularization, is used as a specific example of the combination method. The objective function E(x) to minimize is defined as EðxÞ ¼ jjDx  kjj22 þ l0 jjGx  xjj22 þ l1 jjCðF 1 xÞjj1 þ l2 jjrðF 1 xÞjj1

[2]

where x is the reconstructed k-space of all channels, k is the partially acquired data, D is the operator for data undersampling, G is a general convolution operator for SPIRiT, C is the wavelet transform, and r is the gradient operator. The nonnegative regularization parameters l0 , l1 , and l2 , balance these four terms. The last two terms are sparsity constraints; they enforce the sparsity of the reconstructed image in these transformed domains. In Eq. [2], two sparsity constraints are used as an example. The number of sparsity constraints can be more or less than 2. F 1 ð•Þ is defined as below and iFFT is inverse fast Fourier transform vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u Nch uX 1 [3] F x¼t ðiFFTðxj ÞÞ2 j¼1

The reconstruction process is to minimize the objective function to solve x ¼ arg min EðxÞ. The accuracy of the CS PPI reconstruction using Eq. [2] relies on several factors, including the accuracy of the convolution kernels which implicitly uses the coil sensitivity and the regularization parameters. In the proposed method, sharable information among multicontrast images is used to improve these factors to solve the existing problems in L1-SPIRiT.

Sharable Information Using spin echo sequence as an example, the MRI signal S can be defined by the equation below (a simplified version): S¼

þ B 1 rFðB1 ; B0 Þð1

TR=T1

e

TE=T2

Þe

[1]

where r is the proton density which is a function of main magnetic density B0, transmit magnetic density Bþ 1 and

Sequential Reconstruction and Consideration of Acquisition Order The usage of the sharable information can be in either a sequential or a joint way. In the sequential way, the sharable information is extracted from previously reconstructed or acquired data and applied in the reconstruction of the following images. Images of each contrast are reconstructed

Improving CS-PPI Methods for Multicontrast Imaging

individually. In the joint way, the multicontrast images are reconstructed jointly with constraints on the similarity of sharable information (13). Currently, a clinical imaging protocol includes several scans processed sequentially and image of each scan is presented immediately after the acquisition. Hence, we propose to use the sequential way for the usage of sharable information because it follows the clinical convention. If a joint reconstruction scheme is used, the images of two or more scans need be reconstructed and presented together. The presentation, as well as the possible feedbacks from the acquired scans, has to be delayed. The clinical acceptance of joint-reconstruction is uncertain. Due to the sequential usage of sharable information, the order of acquisition should be decided by the speed and sensitivity to motion of these sequences to achieve higher total acceleration factor and more consistent image quality. For the first scan, there is no or very limited prior information. The faster and/or motion insensitive sequences should be used with low acceleration factor to provide sufficiently accurate sharable information without long acquisition time. After the reconstruction/acquisition of the first image, a set of full k-space with high resolution anatomical information is available. Hence, manifold sharable information (coil sensitivity maps and structural information) can be extracted as prior information for the following reconstructions with slower sequences at higher reduction factors. For example, T1w sequence is faster than T2w sequence in brain imaging. Hence, the T1w image can be acquired with low acceleration factor, say 1 to 2. Slower and motion sensitive T2w sequences are then acquired with higher acceleration factors, say 4 to 5, because extracted sharable information is available to enhance the reconstruction. Parallel-imaging and compressed sensing Reconstruction of Multicontrast Imaging Using SharablE Information (PROMISE) Based on the ideas above, a general framework, PROMISE, is proposed to exploit the sharable information among

FIG. 1. Flow chart of PROMISE.

525

images with different contrasts for fast imaging. Figure 1 shows the flowchart of the framework. Without loss of generalization, it is assumed that the examination is composed of three scans: Prescan for initial coil sensitivity maps, T1w scan, and T2w scan. As shown in Figure 1, the T1w image is acquired with equally spaced 1D-undersampling scheme at an acceleration factor of 2. The T2w image is acquired with higher acceleration factors with randomly 1D-undersampled trajectory. There are four major processes in the reconstruction: (i) reconstruction of the T1w with coil sensitivity maps from prescan; (ii) extraction of sharable information (implicit coil sensitivity information and spatially adaptive regularization parameters defined with structure information) from the reconstructed T1w image; (iii) partial reconstruction of the T2w image acquired using GRAPPA operator; and (iv) final reconstruction of the T2w image using a modified version of Eq. [2] with the information extracted from step 2. The implementation details are provided in the Methods section.

METHODS Extraction of Sharable Coil Sensitivity Information In this specific example, low-resolution sensitivity maps from prescan were used for T1w scan reconstruction. After the reconstruction of the T1w image, a full k-space of each individual channel was produced by projecting the sensitivity weighted image back to k-space. The 64 center k-space lines from the T1w image were used to calibrate convolution kernels for SPIRiT (16) and GRAPPA operator (17) to implicitly use the coil sensitivity information in k-space, which were shown to be more robust to interscan motion than sensitivity maps in the image domain. Because sufficient k-space data are used for the calibration, the convolution kernels are more accurate than those using a few self-calibration data. In addition, great numbers of ACS lines are no longer required for future scans.

526

Gong et al.

Wn ¼ ðpn þ dÞq ; 1  Wn  dq

Extraction of Sharable Structure Information as Spatially Adaptive Regularization

Wi ¼ ðpi þ dÞq ; 1  Wi  dq

The sharable structure information was extracted from previously acquired/recovered images using statistical estimation. The statistical probability of representing boundaries/edges of each component in sparsity transform, p, is computed using a Hidden Markov Tree based Expectation Maximization (EM) (18): p ¼ pðfS contributes to boundariesgjCm; QðCmÞÞ

[4]

where pðÞ is the operator for probability calculation, Cm is the wavelet transform of image m, S is the coefficient of the wavelet transforms, QðCmÞ is the superparameter characterizing the properties of the wavelet transform coefficients. The details of derivation are included in the Appendix, and more explanation can be found in Gong et al (19) and Duarte et al (20). To reduce the sensitivity to motion, the extracted edge information was not directly used to enforce the similarity of the edges between reference image and the to-be reconstructed image. Instead, it was used to define a spatially adaptive regularization parameter. The regularization parameter l (l1 , l2 , or both) in Eq. [2] balances data fidelity and the regularization. A larger l emphasizes regularization more, but may cause reduced sharpness and loss of fine structures; a smaller l emphasizes data fidelity more but may result in residual noise/artifacts. It has been shown in references (21– 23) that a spatially adaptive regularization based on prior image or noise distribution, rather than using a constant regularization parameter, can result in better balance between maintaining signal to noise ratio (SNR) and perserving boundaries. Based on similar but not the same idea, spatially adaptive regularization in the sparsity transform domains was used in this work, to preserve boundaries and fine structures while sufficiently remove noise/artifacts. The adaptive regularization implicitly and softly enforces the sparsity supports from the fuzzy structure information estimated from previously reconstructed images considering the possible existence of minute mismatches of structures. Hence, the objective function in Eq. [2] can be modified with estimated sharable information. EðxÞ ¼ jjDx  kjj22 þ l0 jjGx  xjj22 þ l1 jjWn CðF 1 xÞjj1 þ l2 jjWI rðF 1 xÞjj1 [5] Here the L1 norms of sparsity transforms are used as regularization with spatially adaptive weights Wn and WI , for wavelet and Total-Variance transform respectively. At locations without boundaries/edges or significant-value coefficients, a larger parameter is used to sufficiently suppress noise/artifacts; at locations near boundaries/edges or contains significant-value components, a smaller parameter is used to avoid blurring. The weighted Total-Variance transform is a nonquadratic analog to the smoothness based functional proposed in previous study (22), while the wavelet transform characterizes structures similarity robustly in multiple scales of resolutions. The weights are defined by a set of point-wise weighting functions (24).

[6]

pn and pi are the statistical probabilities of being boundaries/edges in wavelet and image space, which are computed using Eq. [3]. The parameter d is for overflow prevention; and parameter q is for controlling the order of regularization. Constant parameters, k1 and k2 in Eq. [2] and d and q in Eq. [4], are defined empirically (k1 and k2 are empirically set to be 0.001, d is set to be 1010 and q is 0.02).] PPI Enhanced Initialization for Data Fidelity Term Similar but not the same to the method proposed in (25), a generalized autocalibrating partial parallel acquisition (GRAPPA) operator (17) was applied to estimate a PPI based initialization for the recovery of partial k-space, which has the distance not larger than Dky to the acquired data. There are two reasons for this implementation. First, when random trajectory is used, GRAPPA operator is applicable because it does not require specific sampling patterns and is faster than other k-space methods due to the small convolution kernel. Second, the result of GRAPPA operator is accurate enough to be used as reference when the distance between the extrapolated data and the acquired data is not larger than Dky (26). Hence the result of GRAPPA operator was also used in the fidelity term in Eq. [5], in addition to initialization. Data Acquisition Three sets of 2D multicontrast images were acquired on a Philips 3 Tesla (T) system (Philips Healthcare, Best, Netherland). Two sets of the images are multislice axial brain images acquired with an eight-channel head coil (Invivo Corporation, Gainesville, FL). Each set was composed of three scans of brain: low-resolution T1w prescan, high resolution T1w and T2w scans. For the dataset 1, the volunteer was asked to keep still during the scan. To validate the proposed method with interscan rigid motions, the volunteer was asked to move head between the T1w and T2w scans when dataset 2 was acquired. The two T1w datasets used a fast field echo (FFE) sequence (FOV 230  230 mm2, repetition time (TR) 170 ms, echo time (TE) 3.9 ms, flip angle 80 , slice thickness 5 mm). The matrix sizes for the low and high resolution images are 64  64 and 256  256, respectively. The T2w dataset also used an FFE sequence (FOV 230  230 mm2, TR 3000 ms, TE 80 ms, flip angle 90 , slice thickness 5 mm). The matrix size is 256  256. For all datasets, the PE direction was left–right. The acquisition time for the 16 slices of T1w images and T2w images were 86 seconds and 108 s, respectively. To validate the sensitivity of the proposed method to nonrigid interscan motion, dataset 3 composed of two carotid MR scans was used. The two scans were 2D quadruple-inversion-recovery T1-weighted (QIR T1w) scan and 2D double-inversion-recovery T2-weighted (DIR T2w) scan. Both scans used turbo spin echo (TSE) sequence. A homemade 36-channel neurovascular coil (27), which provides 11 channels for carotid imaging

Improving CS-PPI Methods for Multicontrast Imaging

with two elements in phase encoding directions (anterior–posterior), was used in the experiments. The 16 axial slices are acquired (FOV 160160 mm2, voxel size 0.6  0.6 mm2, slice thickness 2 mm). The total acquisition time for QIR T1w scan (TR 800ms, TE 10ms, TSE factor 10) was 6 min 11 s and the total acquisition time for DIR T2w scan was 3 min 40 s (TR 4800 ms, TE 50 ms, TSE factor 12). Both datasets were fully acquired, and retrospectively undersampled in the experiments. For T1w data, uniform 1D undersampling along PE direction was used with reduction factor of 2. For T2w data, a random 1D undersampling scheme along PE direction was used for reconstruction schemes (shown in Table 1). The same 1D Variable Density random undersampling trajectory was used L1-SPIRiT, Joint Bayesian CS and PROMISE, which was generated based on point spread function which was proposed in literature (3) to optimize the incoherence of the undersampling pattern. Design of Experiments Four sets of experiments were designed to evaluate PROMISE. The first set of experiments was used to demonstrate the advantages of using each kind of sharable information. The set of extracted information, coil sensitivity information and structure based spatially adaptive regularization parameters, were added to a L1-SPIRiT based model one-by-one. In the second set of experiments, PROMISE was compared with existing methods listed in Table 1. Both sets of experiments used in vivo dataset 1. Two multicontrast brain imaging datasets with simulated and actual interscan motion, respectively, were used in the third set of experiments to validate the tolerance to rigid interscan motion. Dataset 1 was used to simulate the in-plane motion. One slice of composite T2w image was first synthesized from multichannel fully sampled T2w scan data. The composite T2w image was artificially displaced from the original position by twoby-two pixels (two pixels down along both FE direction and two pixels right along PE direction equivalent to approximately 3 mm displacement) and rotated by 5 counterclockwise around the center of slice. Then a sensitivity map is multiplied to the single image to synthesize a set of multichannel T2w image with the simulated interscan motion. Dataset 2 contains actual interscan

527

motion. Both in-plane and through-plane motions between the T1w and T2w scans can be clearly observed. In addition to the reconstruction using PROMISE, L1SPIRiT, and Joint Bayesian based CS was also used for comparison. The carotid dataset was used to test the robustness to nonrigid interscan motion in the fourth set of experiments. In this experiment, it is assumed that the T2w image which needs much shorter acquisition time was scanned first with full k-space. The T1w image was scanned with acceleration and reconstructed using sharable information extract from T2w scan. L1-SPIRiT and Joint Bayesian based CS were also used for comparison. Implementation of Existing Methods for Comparison For comparison, L1-SPIRiT and Joint Bayesian CS were implemented. Because the target of this work is to improve CS-PPI, specifically taking L1-SPIRiT as a basic CS-PPI framework, we implemented L1-SPIRiT according to the references (7,16), and the example codes (http:// www.eecs.berkeley.edu/mlustig/Software.html) were referred. A 7  7 convolution kernel, 20 ACS lines and Tikhonov regulation with weight of 0.01 were used for kernel calibration. The regularization parameter l was empirically optimized to be 0.001. The wavelet softthreshold parameter for L1-SPIRiT was empirically optimized to be around 0.002 to achieve a balance of aliasing reduction and block-artifact reduction. Nonlinear iterative reconstruction methods were used for L1-SPIRiT and PROMISE. The number of iterations follows the sampled code of L1-SPIRiT mentioned above to guarantee the minimum root mean square error (RMSE) was reached and the result with minimum RMSE was selected. Because Bayesian CS uses sharable structural information to improve CS, we also implemented it for comparison. The code was downloaded from online software package (http://web.mit.edu/berkin/www/software.html). We implemented the joint reconstruction function for complex image of each channel and then synthesized into one image using root-sum-of-square as in equation [3]. Image Quality Assessment For accuracy assessments, both the magnitude and the distribution of the error were measured. The images

Table 1 Reconstruction Methods Used in Comparison Acquisition Chosen method L1-SPIRiT (16,34) Joint Bayesian based CS (13) Proposed method: PROMISE

2D Multislice 1D trajectory Random undersampling Random undersampling

3D Isotropic 2D trajectory Poisson disk

Random undersampling

Poisson disk

Reconstruction

Sharable Information

PPI

CS

3

3

Sharable sensitiv- Sharable structure ity information information

3

3

3

3

3

3

528

reconstructed using the full k-space data were used as reference images. The normalized RMSE in a region-ofinterest (ROI) was used to evaluate the magnitude of error: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi .X X RMSE ¼ jIrecon  Iref j2 [7] jIref j2 where Iref is the reference image reconstructed with fully acquired k-space data and Irecon is the reconstruction with partially acquired k-space data. Error maps illustrating the spatial distribution of the differences in magnitude between the reconstruction and the reference were also presented. The ROI for error computing is defined by the image support thresholding out the background of the reference image. All the algorithms mentioned before were implemented using MATLAB (MathWorks Inc., Natick, MA) with Intel Core 2.4 GHz Duo CPU and 64-bit Operating System.

RESULTS Figure 2 shows one example of the extracted boundary information. Figures 2c and 2d are the extracted boundary information from T1w brain image in Figure 2a for the reconstruction of T2w image in Figure 2b. Figures 2c and 2d are in the wavelet and image domains, respectively. It can be clearly seen that the extracted bounda-

Gong et al.

ries are banded sparse, which are defined by the probability as boundaries. Figure 3 demonstrates the performance improvement of different components in PROMISE. By adding different component of manifold sharable information (shown in the table below the figures) one-by-one to conventional L1-SPIRiT, the reconstructed results at reduction factor of 5 are shown: L1-SPIRiT with random trajectory (Fig. 3b), L1-SPIRiT using shared sensitivity information from T1w data (Fig. 3c), L1-SPIRiT with spatially adaptive weight estimated from structure information (Fig. 3d) and PROMISE (with both sharable sensitivity information and structure information) (Fig. 3e). Corresponding error maps are five-times brightened and shown in Figure 3g–j. The detail of the methods and normalized RMSEs are shown in the table below. Figures 4 to 6 demonstrate the comparison of reconstructions using the 2D multislice in vivo brain dataset, which compares the proposed algorithm PROMISE with other algorithms listed in Table 1. Dataset 1 with net reduction factor 5 was used. Figure 4a is the reference image reconstructed with full k-space data. Figure 4b shows the 1D random undersampling trajectory used for all three methods (L1-SPIRiT, Joint Bayesian and the proposed method). The lower three rows of Figure 4 show the reconstructions (left) and the corresponding error maps (right) of L1-SPIRiT (Fig. 4c,d), Joint Bayesian based CS (Fig. 4e,f), and the proposed PROMISE (Fig.

FIG. 2. Demonstration of the extracted structural information. a,b: The T1w (a) and T2w (b) images. The statistical boundary information is extracted from the T1w image reconstructed with reduction factor of 2. c: The extracted boundary information in the wavelet domain. d: the extracted boundary information in the image domain.

Improving CS-PPI Methods for Multicontrast Imaging

529

FIG. 3. The improvements of T2w scan reconstruction by using manifold sharable information at reduction factor of 5. a: The reference image. b: L1-SPIRiT. c: L1-SPIRiT þshared coil sensitivity information (in k-space). d: L1-SPIRiT þshared structure information used for adaptive weight. e: L1-SPIRiT þshared coil sensitivity information þ shared structure information (PROMISE). f: the T1w image used for multicontrast information extraction. g–j: Corresponding difference maps which are five-times brightened. Different components of PROMISE are used as the table shows below and the resulting RMSE.

4g,h), respectively. The errors maps were brightened five times for better visibility. The values at the right-bottom of the error maps are the corresponding normalized RMSE. Figure 5 demonstrates the corresponding zoomed-in regions, which is defined in Figure 4a. Figure 6 shows the comparison of RMSE at different net reduction factors from 2 up to 6. In Figure 6a, the dotted line, dashed line, and solid lines are for L1-SPIRiT, Joint Bayesian CS, and PROMISE. In Figure 6b, reconstruction RMSEs while exploiting different sharable information are shown with labeling. The combinatory usage of sharable sensitivity and structure information significant improves the reconstruction. Figures 7–9 show the performance of PROMISE when there were different types of interscan motions. Dataset 1 with simulated in-plane motion, dataset 2 with both inplane and through-plane mismatches, and dataset 3 with nonrigid motion was used for Figure 7 to Figure 9, respectively. For comparison, the results of L1-SPIRiT and Joint Bayesian CS are also provided. Reduction factor of 4 was used in all three Figures. Both the reconstructed image and the corresponding error maps were presented in lower rows for L1-SPIRiT, Joint Bayesian CS and PROMISE respectively. All three methods shared the same 1D random undersampling trajectory. The error maps were brightened five times for better visibility. The normalized RMSEs were shown at the right bottom of each error maps.

DISCUSSION Advantages of Using Sharable Information This study proposes to sequentially use the sharable information to enhance fast multicontrast imaging and resolve currently existing problems in CS-PPI. As one specific example, sharable sensitivity and structure information from previous reconstructions was used to improve the reconstruction of multicontrast images using CS-PPI, specifically L1-SPIRiT. Figures 46 demonstrated that PROMISE resulted in lower RMSE and better structure preservation than existing methods when there was no interscan motion. Figures 7–9 demonstrated that using sharable information still resulted in better reconstruction when there were different kinds of interscan motions. PROMISE is superior to the existing methods because that PROMISE used manifold prior information shared in scans. In the specific implementation, both sensitivity and structure information was estimated and exploited in PROMISE. The contribution of individual kind and their combinations were shown in Figure 3. The normalized RMSE was reduced from 11.5% to 7.6% by using manifold sharable information, and each component (sharable sensitivity and structural information) of PROMISE has contributions to the improvement. As Figures 3–6 show, PROMISE consistently resulted in lower RMSE than L1-SPIRiT and Joint Bayesian CS for the

530

FIG. 4. Comparison to the existing algorithms (L1SPIRiT, Joint Bayesian based CS, and the proposed method PROMISE) for the reconstruction of simulated T2w brain image at net reduction factor of 5. The left column shows the image reconstructed with full k-space (a), L1-SPIRiT (c), Joint Bayesian based CS (e), and PROMISE (g). An ROI for zooming-in is shown in a. The same 1D random undersampling trajectory for all the methods is shown in b. The images in the right column are the corresponding error maps (d–h). The error maps were brightened five times.

same undersampled data, and the superiority increased with the net acceleration factor. Sensitivity to Interscan Motion As shown previously, the usage of sharable information can improve the reconstruction dramatically when there is no interscan motion. However, if there are interscan motions, then the coil sensitivity and the structure information would be misregistered. Models in the reconstruc-

Gong et al.

FIG. 5. The corresponding zoomed-in images from Figure 4. The zoomed in region is defined in Figure 4a. a: Corresponds to Figure 4a. b–g: Correspond to (c–h) in Figure 4. RMSEs in the zoomed-in ROI are shown at the right bottom at each zoomed-in error map.

tion strongly enforcing the similarity of misregistered information would result in artifacts (Figs. 8e,f and 9e,f). Other previously proposed algorithms working in the image domain also suffer the modeling error. PROMISE, however, uses two approaches to moderate the sensitivity to motion. The first one is to share implicit coil sensitivity information; the other one is to softly use the shared boundary information as spatially adaptive regularization parameters. In PROMISE, coil sensitivity is implicitly used as convolution kernels in GRAPPA operator and SPIRiT. These small size convolution kernels correspond to lowresolution coil sensitivity maps. Moderate motion does not change the sensitivity convolution kernel much. Therefore, adjacent slices or time frames can share the same set of

Improving CS-PPI Methods for Multicontrast Imaging

531

FIG. 6. a: Plot of RMSEs in image support of the T2w scan of dataset 1 at reduction factors from 2 up to 6. The lines represent the reconstructions using L1-SPIRiT (dotted), Joint Bayesian CS (dashed), and the proposed PROMISE (solid) respectively. b: Shows the plots of the proposed method PROMISE with different sharable information used (circle, sensitivity only; triangle, structure only; star, proposed usage of both sharable information), at different reduction factor from 2 up to 6. The ROI is defined by the image support extracted from reference T2w image.

convolution kernels. This feature has been used in literatures, such as zGRAPPA (28) and TGRAPPA (29). Second, PROMISE does not enforce that the multicontrast images explicitly share the same supports of sparsity or locations of boundaries. The shared structure information is only used to define regularization parameters. Similar ideas have been proposed to achieve robust recovery by using fuzzy prior information in weighted regularization (23). Hence, if there are interscan motions, even though the improper regularization parameter will result in nonoptimal balance of detail preservation and artifact reduction, no change of contrast or other artifacts will be introduced by using sharable information between scans. To further moderate this potential problem robustly, multilevel resolution approach is used in the extraction and implementation of structure information in both the wavelet domain and image domain. As shown in Figure 2, instead of just a single pixel width curve of deterministic binary values, an edge in image could result in bands of weighted pixels in multiple resolution levels. These two approaches help the proposed algorithm to robustly exploit the fuzzy structure information, and hence moderate motions in pixel level cannot change the result noticeably. Due to the usage of these two approaches, PROMISE consistently resulted in better reconstruction in our experiments than L1-SPIRiT and Joint Bayesian CS (Figs. 7–9) even when there were interscan motions. On the other hand, it is also true that the advantages of PROMISE over L1-SPIRiT, for exploiting structure information, will be gradually reduced with the increase of motion level. However, an efficient single step registration operation after the first few iteration of PROMISE can resolve the problem because PROMISE is robust for moderate interscan displacement as discussed before. The results for PROMISE with registration will be reported separately. Clinical Applicability Besides interscan motion, another issue of existing methods using sharable information is long reconstruc-

tion time (9,13). Compared with these existing methods, PROMISE has lower computational cost. PROMISE only has one additional step than L1-SPIRiT to extract sharable information. The extraction of boundary information took approximately 10 s in our implementation using Matlab. The remaining steps for the calibration of convolution kernels and optimizing function in Eqs. [2,5] took the same computational cost as L1-SPIRiT. When the target of the acquisition acceleration is to reduce motion artifacts, or breath-hold time, it is worthwhile to spend extra seconds for the information extraction between the scans. Moreover, the calculation can also be processed in parallel with the preparation steps for the following scans. Generalization of PROMISE This work focuses on 2D applications because multislice 2D imaging is widely used for multicontrast imaging clinically. Moreover the acquisition time could be even longer than 3D acquisition, as shown in Supplementary Table S1, which is available online. Faster 2D acquisition does not only save acquisition time, but also reduce motion artifacts, moderate geometry distortion artifacts in EPI acquisition (30) and improve the patient comfort. However, the concept of the proposed method can also be easily extended to 3D. One example of the application in 3D imaging was demonstrated in Supplementary Figure S1 in the supplemental material. The set of T1w and T2w images were acquired on a Philips 3T system with an eight-channel coil. Poisson disk distributed undersampling trajectory was applied for accelerated 3D sampling and the performance of L1-SPIRiT and PROMISE were compared. In PROMISE, the sensitivity kernels were computed based on T1w scan and the structural information was extracted using the proposed scheme in each plane of the 3D dataset. PROMISE resulted in lower RMSE (6.1%) than L1-SPIRiT (7.2%) at reduction factor of 12. Thus, the usage of sharable information helped the preservation of boundaries in 3D multicontrast imaging as well.

532

FIG. 7. Comparison of the sensitivity to rigid interscan motion (simulated in-plane motion of dataset 1) of L1-SPIRiT, Joint Bayesian, and PROMISE. The T2w image (a) contains simulated displacement from the T1w image (b). The T2w image was artificially translated by two pixels to the right and two pixels to the bottom and rotated by 5 degrees counterclockwise and then synthesized with sensitivity maps. Reconstructed T2w image with net reduction factor of 4 and the corresponding error maps are shown in c,d (L1-SPIRiT), e,f (Correlation Imaging), and g,h (PROMISE), respectively. The error maps are five times brightened and the RMSE in image support based on reference T2w image is labeled at the right-bottom of each error-map.

Limitation and Future Work This study reported the preliminary results of PROMISE. The major target of the study is to explain the advantages of using manifold sharable information (sensitivity, structure, etc.) in a sequential approach. More detail explanation of the boundary extraction step is provided in appendix.

Gong et al.

FIG. 8. Comparison of the sensitivity to rigid interscan motion (actual in-plane and out-of-plane motions in dataset 2) of L1-SPIRiT, Joint Bayesian CS, and the proposed PROMISE method. The fully sampled and reconstructed T2w image is shown in a. PROMISE exploits sharable information extracted from the T1w scan (b). Reconstructed image with reduction factor of 4 and the corresponding five times brightened error maps are shown in c,d (L1SPIRiT), e,f (Joint Bayesian CS), and g,h (PROMISE), respectively. RMSE in image support based on reference T2w image is labeled at the right bottom of each error-map.

The manuscript only provides a specific example set of sharable information between multicontrast scans. Modifications, especially the usage of additional sharable information, can be adopted for acquisition and reconstruction for different applications. For example, with previous reconstructed k-space as a reference, undersampling trajectory can probably be optimized (31,32) to improve the acquisition for multicontrast imaging with sharable information as well. To focus on the major target of this study, we will present a clinical applicable trajectory optimization scheme for multicontrast imaging in a separate work.

Improving CS-PPI Methods for Multicontrast Imaging

533

in the better balance of detail preservation and noise reduction than other tested methods (RMSE 7.6% for statistical method and 8.9% for Distance Transform approach when R ¼ 5). Due to the limitation of space, further discussions and experimental results have to be presented separately. CONCLUSIONS In this study, the idea of sequentially exploiting sharable information among multicontrast scans is introduced to resolve the challenges in multicontrast MRI reconstruction using Parallel Imaging and Compressed Sensing. A specific embodiment, PROMISE, shows significant and stable improvement for L1-SPIRiT and demonstrates advantages compared with previously proposed methods, when sharable information with or without interscan motion was used. APPENDIX Here, we present methods to quantify and determine the probabilities of representing boundaries for components in the wavelet domain (pn ) and image domain (pi ). 1. Extraction of Boundary Information in Wavelet Domain

FIG. 9. Comparison of the sensitivity to nonrigid interscan motion (actual in-plane motions and through-plane motions in carotid dataset) of L1-SPIRiT, Joint Bayesian, and the proposed PROMISE method. Reference T1w image (a) and reconstructed T2w image (b) with interscan motion are shown in the first row. PROMISE exploits sharable information from the T2w scan to reconstruct the undersampled T1w scan. A reconstructed T1w image with reduction factor of 4 and the corresponding five times brightened error maps are shown in c,d (L1-SPIRiT), e,f (Joint Bayesian CS) and g,h (PROMISE), respectively. The RMSE in image support based on the T2w image (from which sharable information is extracted) is labeled at the right bottom of each error-map. Zoomed-in images of the left artery, at the position indicated in the yellow square in a, are two times brightened and shown in the left bottom of each subfigure.

In this work, we proposed a statistical algorithm as a specific example to extract structural information in a sequential approach. There are other schemes for the same purpose, such as Distance Transform (33) of Canny Boundary operation results. Based on simulation and experiments, we found that the proposed method results

In the wavelet domain, boundaries are encoded as blocks of significant wavelet coefficients in multiple scales of resolutions. In contrast, piecewise constant components and additional noises are encoded as small wavelet coefficients in most scales and impulses in high resolution frequency. Thus, exploiting the properties of the wavelet domain can help to distinguish boundaries from background and noises. We implemented and improved the estimation method in the wavelet domain (24) to extract the statistical structure information. A Hidden Markov Tree (HMT) model was used to model the probability density function p of each wavelet coefficient un ¼ ðCmÞn as a Gaussian mixture density with a hidden binary state Sn , which in our application is either “Negligible” (Sn ¼ N) or “Significant” (Sn ¼ S). pð•Þ 2 ½0; 1; pðSn ¼ SÞ þ pðSn ¼ NÞ ¼ 1

[A1]

Extraction of structures at multiple levels of resolution can be achieved inherently due to the properties of wavelet transform (24). There are three properties of the Hidden Markov Quad-Tree model based on the properties of 2D wavelet transform. First, the values of wavelet coefficients obey a Gaussian mixture distribution. Second, the hidden states of coefficients in different scales in 2D wavelet quad-trees are correlated. The state that changes to daughter coefficient from a mother coefficient n can be formulated using a Transition Matrix Tn . Third, both the Gaussian mixture distribution and Transition Matrices between states are scale-dependent and can be represented by exponential functions of the scale J in wavelet-trees. These properties can be characterized by following equations. The mixture distribution function forms dual-Gaussian distribution with hidden state and scale dependent variance:

534

Gong et al.

f ðun Þ ¼ f ðun jSn ¼ SÞpðSn ¼ SÞ þ f ðun jSn ¼ NÞpðSn ¼ NÞ 2

f ðun jSn ¼ SÞ ¼ Nð0; sS;n Þ ¼ Nð0; CsS 4 2

JaS

f ðun jSn ¼ NÞ ¼ Nð0; sN;n Þ ¼ Nð0; CsN 4

Þ

J aN

Þ [A2]

The transition of the hidden states is modeled using transition matrix with scale dependent transition probability: 2

pðSn ¼ SÞ

4

3

2

pðSm ¼ SÞ

5 ¼Tn 4 pðSn ¼ NÞ

3

2

5¼4 pðSm ¼ NÞ

2

pðSm ¼ SÞ

4

pn S!S

pn N!S

S!N

N!N

pn

pn

3 5

ferent levels of resolution, MSTV can be generated using Eq. [A6]. MSTV shares similar properties as wavelet and can be used to extract structure information and estimate statistical regulation weights, as discussed in previous session. Then the estimated probability is used as regularization weights in total variance transform. In this way, we extracted boundary information in the image domain. In addition, the proposed method is supreme to conventional edge operators because it is efficient and able to correct false edges using multiresolution information.

3

REFERENCES

5

1. Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: sensitivity encoding for fast MRI. Magn Reson Med 1999;42:952– 962. 2. Griswold MA, Jakob PM, Heidemann RM, Mathias Nittka, Jellus V, Wang J, Kiefer B, Haase A. Generalized Autocalibrating Partially Parallel Acquisitions (GRAPPA). Magn Reson Med 2002;47:1202– 1210. 3. Lustig M, Donoho D, Pauly JM. Sparse MRI: the application of compressed sensing for rapid MR imaging. Magn Reson Med 2007;58: 1182–1195. 4. Block K, Uecker M, Frahm J. Undersampled radial MRI with multiple coils: Iterative image reconstruction using a total variation constraint. Magn Reson Med 2007;57:1086–1098. 5. Liang D, Liu B, Wang J, Ying L. Accelerating SENSE using compressed sensing. Magn Reson Med 2009;62:1574–1584. 6. King KF. Combining compressed sensing and parallel imaging. In Proceedings of the 16th Annual Meeting of ISMRM, Toronto, Canada, 2008. Abstract 1488. 7. Lustig M, Alley M, Vasanawala S, Donoho DL, Pauly JM. L1 SPIR-iT: autocalibrating parallel imaging compressed sensing. In Proceedings of the 17th Annual Meeting of ISMRM, Honolulu, Hawaii, 2009. Abstract 379. 8. Liu B, Sebert FM, Zou Y, Ying L. SparseSENSE: randomly-sampled parallel imaging using compressed sensing. In Proceedings of the 16th Annual Meeting of ISMRM, Toronto, Canada, 2008. Abstract 3154. 9. Larkman DJ, Batchelor PG, Atkinson D, Rueckert D, Hajnal JV. Beyond the g-Factor Limit in Sensitivity Encoding Using Joint Histogram Entropy. Magn Reson Med 2006;55:153–160. 10. Li Y, Dumoulin C. Correlation imaging for multiscan MRI with parallel data acquisition. Magn Reson Med 2012;68:2005–2017. 11. Li Y, Huang F, Lin W, Duensing GR, Dumoulin C. Correlation Imaging for Multi-Scan Acceleration in Clinical MRI. In Proceedings of the 20th Annual Meeting of ISMRM, Melbourne, Australia, 2012. Abstract 3347. 12. Li Y, Huang F, Lin W, Duensing GR, Dumoulin C. Correlation-based reconstruction using multiple calibration images. In Proceedings of the 20th Annual Meeting of ISMRM, Melbourne, Australia, 2012. Abstract 4253. 13. Bilgic B, Goyal VK, Adalsteinsson E. Multi-contrast reconstruction with Bayesian compressed sensing. Magn Reson Med 2011;66:1601– 1615. 14. Abragam A. The principles of nuclear magnetism. Oxford: Clarendon Press; 1983. 15. Liang D, King KF, Liu B, Ying L. Accelerating SENSE using distributed compressed sensing. In Proceedings of the 17th Annual Meeting of ISMRM, Honolulu, Hawaii, 2009. Abstract 377. 16. Lustig M, Pauly JM. SPIRiT: iterative self-consistent parallel imaging reconstruction from arbitrary k-space. Magn Reson Med 2010;64:457– 471. 17. Griswold MA, Blaimer M, Breuer F, Heidemann RM, Mueller M, Jakob PM. Parallel Magnetic Resonance Imaging Using the GRAPPA operator formalism. Magn Reson Med 2005;54:1553–1556. 18. Bilmes J. A gentle tutorial on the EM algorithm and its application to parameter estimation for Gaussian mixture and hidden Markov models 1997. Technical Report ICSI-TR-97-021 p. 19. Gong E, Wang X, Ying K, Wang S. Wavelet-domain Structure Based Compressed Sensing MRI Reconstruction Using a Hidden Markov

pðSm ¼ NÞ 2 S!S 3 pn N!S pn 5 Tn ¼ 4 S!N N!N pn pn 21 þ CAS 4JgS CAN 4JgN 64 6 ¼6 4 3  CAS 4J gs 1  CAN 4JgN 4

3 7 7 7 5

pn S!S þ pn S!N ¼ 1 pn N!S þ pn N!N ¼ 1 pn ! 2 ½0; 1 [A3] Thus, there are several superparameters of the HMT: Q ¼ faN ; aS ; gN ; gS ; CdN ; CdS ; CAN ; CAS g:

[A4]

Using a HMT based Expectation Maximization (EM) (18), we optimize the superparameters based on the measurement (wavelet transform of the reconstructed image m) and estimate the statistical probabilities as a representation of structure information in wavelet transform. pn ¼ pðSn ¼ SjCm; QðCmÞÞ:

[A5]

The details can be found in Gong et al (19) and Durate et al (20). 2. Extraction of Boundary Information in Image Domain To extract the boundary information in the image domain, we extend conventional TV transform to a multiple-scale TV transform. rXk ði; jÞ ¼ ðXk ði; jÞ  Xk ði; j  1Þ; Xk ði; jÞ  Xk ði  1; jÞÞ Xk ði; jÞ ¼ Xði; jÞ FLPk ði; jÞ [A6] Multiscale Total Variance (MSTV) is analog to a hard wavelet and is defined in [A6], in which Xði; jÞ is pixelwise value in the image and FLPk ði; jÞ is a low-pass filter in the image domain for level k. By using filters with dif-

Improving CS-PPI Methods for Multicontrast Imaging

20.

21.

22.

23.

24.

25.

26.

Tree Model. In Proceedings of the 20th Annual Meeting of ISMRM, Melbourne, Australia, 2012. Abstract 2251. Duarte M, Wakin M, Baraniuk R. Wavelet-domain compressive signal reconstruction using a hidden Markov tree model. In Proceedings of the ICASSP, Las Vegas, Nevada, USA, 2008. p 5137–5140. Huang F, Chen Y, Yin W, Lin W, Ye X, Guo W, Reykowski A. A Rapid and robust method for sensitivity encoding with sparsity constraints: self-feeding sparse SENSE. Magn Reson Med 2010;64:1078– 1088. Haldar JP, Hernando D, Song SK, Liang ZP. Anatomically constrained reconstruction from noisy data. Magn Reson Med 2008;59: 810–818. Haldar JP, Wedeen VJ, Nazamzadeh M, Dai G, Weiner MW, Schuff N, Liang ZP. Improved diffusion imaging through SNR-enhancing joint reconstruction. Magn Reson Med 2013;69:277–289. Crouse M, Nowak R, Baraniuk R. Wavelet-based statistical signal processing using hidden markov models. IEEE Trans Signal Process 1998;46:886–902. Lai P, Lustig M, Beatty P, Brau AC, Vasanawala SS, Alley M. Faster L1-SPIRiT reconstruction with parallel imaging initialization. In Proceedings of the ISMRM 3rd International Workshop Parallel MRI, 2009; Santa Cruz, California, USA. p 25. Lin W, Huang F, Li Y, Reykowski A. GRAPPA Operator for Wider Radial Bands (GROWL) with optimally regularized self-calibration. Magn Reson Med 2010;64:757–766.

535 27. Wang X, Li R, Hayes C, Balu N, Zhao X, Yuan C. A new designed 36channel neurovascular coil at 3T. In Proceedings of the 20th Annual Meeting of ISMRM, Melbourne, Australia, 2012. Abstract 2787. 28. Honal M, Bauer S, Ludwig U, Leupold J. Increasing Efficiency of Parallel Imaging for 2D Multislice Acquisitions. Magn Reson Med 2009; 61:1459–1470. 29. Felix A. Breuer, Peter Kellman, Mark A. Griswold, Jakob PM. Dynamic Autocalibrated Parallel Imaging Using Temporal GRAPPA (TGRAPPA). Magn Reson Med 2005;53:981–985. 30. Bammer R, Auer M, Keeling SL, Augustin M, Stables LA, Prokesch RW, Stollberger R, Moseley ME, Fazekas F. Diffusion tensor imaging using single-shot SENSE-EPI. Magn Reson Med 2002;48:128–136. 31. Seeger M, Nickisch H, Pohmann R, Scholkopf B. Optimization of kspace trajectories for compressed sensing by bayesian experimental design. Magn Reson Med 2010;63:116–126. 32. Xu D, Jacob M, Liang Z-P. Optimal sampling of k-space with Cartesian Grids for Parallel MR Imaging. In Proceedings of the 13th Annual Meeting of ISMRM, Miami, Florida, USA, 2005. Abstract 2450. 33. Borgefors G. Distance transformations in digital images. Computer vision, graphics, and image processing 1986;34:344–371. 34. Lustig M, Alley M, Vasanawala S, Donoho DL, Pauly JM. L1 SPIR-iT: autocalibrating parallel imaging compressed sensing. In Proceedings of the 17th Annual Meeting of ISMRM, Honolulu, Hawaii, USA, 2009. Abstract 379.