APIR4EMC: Autocalibrated parallel imaging reconstruction for extended multi-contrast imaging.

PURPOSE
To improve image quality of multi-contrast imaging with the proposed Autocalibrated Parallel Imaging Reconstruction for Extended Multi-Contrast Imaging (APIR4EMC).


METHODS
APIR4EMC reconstructs multi-contrast images in an autocalibrated parallel imaging reconstruction framework by adding contrasts as virtual coils. Compensation of signal evolution along the echo train of different contrasts is performed to improve signal prediction for missing samples. As a proof of concept, we performed prospectively accelerated phantom and in-vivo brain acquisitions with T1, T1-fat saturated (Fatsat), T2, PD, and FLAIR contrasts. The k-space sampling patterns of these acquisitions were jointly optimized. Images were jointly reconstructed with the proposed APIR4EMC method as well as individually with GRAPPA. Root mean square error (RMSE) to fully sampled reference images and g-factor maps were computed for both methods in the phantom experiment. Visual evaluation was performed in the in-vivo experiment.


RESULTS
Compared to GRAPPA, APIR4EMC reduced artifacts and improved SNR of the reconstructed images in the phantom acquisitions. Quantitatively, APIR4EMC substantially reduced noise amplification (g-factor) as well as RMSE compared to GRAPPA. Signal evolution compensation reduced artifacts. In the in-vivo experiments, 1 mm3 isotropic 3D images with contrasts of T1, T1-Fatsat, T2, PD, and FLAIR were acquired in as little as 7.5 min with the acceleration factor of 9. Reconstruction quality was consistent with the phantom results.


CONCLUSION
Compared to single contrast reconstruction with GRAPPA, APIR4EMC reduces artifacts and noise amplification in accelerated multi-contrast imaging.


Introduction
In magnetic resonance imaging (MRI), multiple contrasts, like T1, T2, proton-density (PD), and FLAIR weighted images, are routinely acquired in clinical practice, as these images provide complementary information for diagnosis [1]. However, acquiring multiple images comes at the cost of prolonged scan time, and thus lower patient comfort and throughput, and more potential for motion artifacts. Therefore, reducing scan time for multi-contrast MRI is desired [2][3][4][5][6][7][8].
Parallel imaging [9][10][11][12][13] is one common method for accelerated imaging by reducing the number of acquired samples in k-space. While more advanced methods exist [14][15][16][17], the basic reconstruction method of parallel imaging either recovers the full k-space (GRAPPA [11] and ARC [12]) or unwraps the aliased image (SENSE [9]) by exploiting the sensitivity profile of a multi-channel coil set. Compressed sensing [18,19], another successful image acquisition acceleration approach, reconstructs the full image from incoherently subsampled acquisition under the assumed sparsity of the image in a specific domain, such as gradient or wavelet [19]. Compressed sensing has been extended to multi-contrast scenarios by performing joint reconstruction, exploiting similarities of spatial structures across different contrasts by assuming the same set of support in the sparsity domain [2][3][4][5][6]8]. By regularizing on both joint and individual sparsity, [20] alleviated the feature leakage across contrasts in compressed sensing method. Alternatively, multiple contrasts can be added as virtual coils in parallel imaging to exploit the correlation among contrasts. The phase independent GRASE reconstruction [21,22] and APIR4GRASE [23] assume spin echo and gradient echo as virtual coils for autocalibrated reconstruction. They exploit the correlation between the different echo types as the sensitivity encoding of virtual coil channels, and reconstruct different contrast weighted images. The JVC-GRAPPA [24] further considers each echo in an echo train as virtual coil for a joint parallel imaging and this improves the reconstructed image quality substantially. Nonetheless, the integrated parallel imaging with the generally used multi-contrast weightings, e.g., T1, T2, or PD weighting, acquired in separate acquisitions has not been explored and validated yet. A recent extension of JVC-GRAPPA [25], using variational neural networks, achieved good image quality with 16fold acceleration for multi-contrast images, although fully sampled data was needed for training of the networks.
Inspired by JVC-GRAPPA [24] and APIR4GRASE [23], we propose the Autocalibrated Parallel Imaging Reconstruction (APIR) for Extended Multi-Contrast (APIR4EMC) method to jointly reconstruct multicontrast acquisitions with optimized k-space patterns. The main challenge in autocalibrated parallel imaging reconstruction of multi-contrast acquisitions is introduced by the different signal evolutions along the echo train for different contrast weightings. In APIR4EMC we propose to compensate these differences by exploiting additionally acquired zerophase echo trains to stabilize the signal evolution. Thanks to this signal stabilization, APIR4EMC goes beyond JVC-GRAPPA [24] and APIR4GRASE [23] by enabling joint parallel imaging reconstruction for widely used multi-contrast imaging protocols with separate weighted acquisitions. In APIR4EMC reconstruction, the k-space of each contrast is regularly subsampled, with a small fully sampled autocalibration signal (ACS) region in k-space center. The sampling patterns of different contrasts can be shifted among each other to allow additional spatial encoding. For a given subsampling factor, we derive an optimal k-space sampling pattern by exhaustive search across all possible patterns. By combining signal stabilization with the optimized k-space pattern, APIR4EMC aims to improve the quality of the reconstructed image in accelerated multi-contrast MRI acquisitions.
We demonstrate APIR4EMC on a combination of commonly used structural imaging contrasts: T1, T1 with fat saturation (T1-Fatsat), T2, PD, and FLAIR weighted images. We acquired all contrast weightings with a variable flip angle 3D fast spin echo sequence (FSE) [26][27][28] with appropriate spin preparation. Compared to the conventional single contrast GRAPPA, APIR4EMC achieved substantially improved image quality in terms of artifacts and noise amplification.

Overview
APIR4EMC regards multiple contrasts as virtual coils, as illustrated with an example in Fig. 1, in a way similar to APIR4GRASE [23] and JVC-GRAPPA [24]. The autocalibrated parallel imaging reconstruction is performed to interpolate unsampled signals for each contrast using signals of all contrasts as in GRAPPA [11]. The image of every contrast is reconstructed using root sum of squares of channel images.
The k-space of each contrast is subsampled regularly, possibly with CAIPIRINHA-like patterns [29]. The k-space center is fully sampled as the ACS region. The k-space pattern for different contrasts is identical but possibly shifted in the phase encoding (PE) directions.

Contrast preparation
From an FSE sequence acquisition, with different TRs, TEs, and/or preparation pulses, signals with multiple contrast weightings can be acquired. Specifically, as a proof-of-concept, we acquire and reconstruct the images for multiple contrast weightings of T1, T1 with fat saturation (T1-Fatsat), T2, PD, and FLAIR, using 3D FSE [26] with variable flip angle [27,28].
In practice, in the acquisition of PD and T2, we use the first half of each echo train for PD and the second half for T2. By distributing each echo train part radially in k-space [30], its contrast is dominated by the contrast at the start of that part, which can be PD or T2 weighted with a proper protocol [31]. Together with T1 and T1-Fatsat, they follow a radial encoding trajectory [30] for the view ordering. For FLAIR, to create a T2 weighted contrast, a linear encoding trajectory is used [30,32].

Stabilization of signal evolution
When performing autocalibrated parallel imaging reconstruction from multi-contrast acquisitions, artifacts can arise due to the signal evolution differences among contrasts along the echo train. The GRAPPA convolution kernel is computed from the ACS region that is filled by a small part of the echo train and only represents the signal evolution along this small part. However, the kernel which encodes correlation among contrasts is used to predict unsampled signals over the whole k-space, where the signal evolution becomes very different and thus the contrast correlation is also different. This can introduce bias in predicting. To prevent this, in APIR4EMC, we propose to compensate the differences by acquiring an additional zero-phase encoded echo train as reference, and dividing each echo of each contrast by the norm of the corresponding reference echo before the reconstruction. Precisely, echo n in every echo train in the k-space is computed as where S zero− phase, n is a vector of all (complex) data acquired by the multichannel coil in the n-th echo of the echo train of length N with zero phase encoded gradient amplitudes, and S orginal, n the originally acquired echo n in every echo train in the k-space. To limit noise amplification, 1/3 2 of the median energy in the zero-phase encoded echoes is added to the denominator to limit the maximum compensation factor to 3 times the maximum compensation factor needed to compensate the 50% of echoes with the largest energy.
With the stabilization of signal evolution, the apparent blurring induced by the decay along the echo train can be reduced as well.
To reduce Gibbs ringing and noise amplification due to the echo train stabilization, while preserving informative structures, a Fermi filter [33] (kT = 0.05, μ = 0.8) is applied in each k-space before reconstruction.

k-space pattern optimization
Similar to the k-space allocation with different echo types in APIR4GRASE [23], the k-spaces of different contrasts are acquired with the same pattern, whereas it is possible to have relative shifts in PE directions between each other. To maximize the acceleration capability of APIR4EMC, we identify an near-optimal k-space pattern using a pattern optimization method similar to the one in APIR4GRASE [23]. For this pattern optimization, we acquire full k-spaces of all contrasts by the same protocols with the same field of view (FOV) as the APIR4EMC acquisition, but in a lower resolution to achieve feasible scan time. All possible subsampled patterns for specific subsampling factors are retrospectively constructed and reconstructed by APIR4EMC. For each pattern we evaluated the root mean square error (RMSE) of the reconstructed images with respect to the reference images reconstructed from the full acquisitions. The pattern with the lowest RMSE is selected. This pattern is expected to be close to optimal, even when resolution or object changes, as the kernel depends mostly on the acquisition coil and distance between k-space samples (1/FOV).

k-space pattern optimization
To design the optimal k-space pattern in a realistic setting, MRI images of a volunteer were acquired after obtaining the Institutional Review Board approval and informed consent. The images were acquired with a 3T clinical scanner (Discovery MR750, General Electric Medical Systems, Waukesha, WI) and an eight-channel birdcage-like receive brain coil (8HRBRAIN). The scan parameters of the T1, T1-Fatsat, T2, PD, and FLAIR contrasts are given in Table 1, with the only difference that these scans were fully sampled with a Matrix Size = 112 × 112 × 88, which resulted in a resolution of around 2 × 2 × 2 mm 3 .
All possible regular patterns with the subsampling factor of 6, i.e. 2 × 3 or 3 × 2, and of 9, i.e. 3 × 3, in PE1×PE2 directions, were retrospectively constructed from the full k-space with an ACS region size of 25 × 25 in PE1×PE2 directions. With subsampling factor of 9, which exceeds the ability of GRAPPA using an eight-channel brain coil, we were able to test if APIR4EMC can still perform a reasonable reconstruction with the additional signals from other contrasts. To reduce computation time, k-space was truncated in the FE direction, preserving the center 25 positions. Each of these datasets was reconstructed by APIR4EMC. The GRAPPA convolution kernel size in APIR4EMC was set to cover a neighborhood patch in k-space with the size of [5,7] in [FE, PE1, PE2] directions, which includes both sampled and unsampled positions. Tikhonov regularization [34] was used in the computation of the kernel. The regularization strength was selected as a value trying to suppress noise amplification while avoiding significant artifacts through visual inspection of the reconstructed image.
The RMSE of the reconstructed image with respect to the reference image was computed for every retrospectively constructed k-space pattern. In the computation of RMSE, the images were normalized such that the highest intensity of all images was equal to 1. The reference image was computed as inverse Fourier transform of the fully sampled kspace followed by root sum of squares channel combination. For each subsampling factor, the k-space pattern leading to the minimum averaged RMSE of all contrasts was selected as the optimal pattern.

Qualitative and quantitative evaluation with phantom
A phantom experiment was conducted to evaluate APIR4EMC both qualitatively and quantitatively. The q-MRI ISMRM system phantom [35] was acquired prospectively with the selected optimal k-space patterns (as shown in Fig. 2) for both subsampling factors of 6 and 9 in the same contrasts and the same scanner and receive coil as in the in-vivo acquisition for the k-space pattern optimization. Scan parameters were the same as well except FOV = 224 × 224 × 224 mm 3 and Matrix Size = Besides APIR4EMC reconstruction, GRAPPA reconstruction was also performed on the k-space of each individual contrast for comparison. The GRAPPA convolution kernel size was set to cover the same size neighborhood as in the previous experiment for k-space pattern optimization, and Tikhonov regularization was used as well for the kernel computation. For both methods, the reconstructions without stabilization of signal evolution, and with stabilization but without Fermi filtering, were also performed to inspect their effects on the reconstructed image. Without filtering and stabilization the reconstruction method reduces to APIR4GRASE [23].
The RMSE of the reconstructed image with respect to the fully sampled reference image was computed for both methods. To allow for a fair comparison, the reference images for the RMSE computation were reconstructed with the same pre-processing options (i.e., with or without stabilization and Fermi filtering) as the accelerated images. Since the stabilization and Fermi filtering may affect the energy of the image, the images were divided by the L2 norm of the corresponding reference image, to make the RMSE values of different pre-processing options comparable. The g-factor map was also computed for both APIR4EMC and GRAPPA with the pseudo multiple replica method [36] with 100 iterations by adding Gaussian white noise to the acquired kspace. The magnitude level of the simulated noise was estimated by the standard deviation of a border region of the k-space.

In-vivo validation with high resolution 3D isotropic multi-contrast brain imaging
We validated APIR4EMC for in-vivo acquisitions with high resolution 3D isotropic multi-contrast brain imaging. For this, in the same scanning session, two separate prospective acquisitions were performed with the selected optimal k-space patterns. The scans had an acceleration factor of 6 and 9, respectively. The scan parameters are shown in Table 1. The effective scan time was 10.6 min and 7.5 min with subsampling factors of 6 and 9, respectively. APIR4EMC and GRAPPA reconstructions were performed. The reconstruction settings, e.g. the size of the GRAPPA convolution kernel, and the employment of regularization, were the same as the phantom experiment. For both methods, the reconstructions without stabilization, and with stabilization but without Fermi filtering were performed as well to see their effects on the in-vivo images.
The image quality was compared visually for contrasts, artifacts and the SNR level between the GRAPPA images and the APIR4EMC images. Only visual evaluation was performed for the in-vivo experiments, since quantitative evaluation would require a fully sampled image as reference. In our experiment, a fully sampled 3D acquisition with matrix size 224 × 224 × 178 would need an effective scan time of around 74 min. Together with the prospectively subsampled acquisitions, that is prohibitively long for a volunteer scan. Additionally, we expect artifacts caused by subject motion during/between such long scans to dominate the error measure. Therefore, we consider reliable quantitative in-vivo evaluation not feasible.
To estimate how much influence each of the five contrasts had on the reconstruction of each contrast, we performed additional reconstruction experiments for APIR4EMC. While still using the same k-space pattern and the same computed GRAPPA convolution kernel, we set the k-space of one of the five contrasts to zeros, and kept the rest the same as the original to reconstruct the images. Each contrast was reconstructed with the leave-one-contrast-out method to get rid of the contribution from this left-out contrast. We compared the image quality visually between the leave-one-contrast-out reconstructions and the original reconstructions. For the subsampling factor of 9, the average RMSE was 2.45 × 10 − 2 and the optimal RMSE attained was 2.20 × 10 − 2 . Fig. 3 shows the reconstructed images with RMSEs and g-factor maps for both APIR4EMC and GRAPPA. To have a fair comparison, the shown GRAPPA reconstruction is with stabilization of signal modulation and with Fermi filtering as well. Compared to the GRAPPA images, which suffered from substantial artifacts, APIR4EMC shows lower noise amplification and less severe artifacts in each contrast image. Fig. 4 presents the zero-phase encoded signal evolution for the five contrasts where obvious differences in modulation along the echo trains of different contrasts can be observed. Fig. 5 shows the reconstructed images for both methods with and without stabilization of signal modulation with acceleration factor of 6. For the reconstructions with stabilization of signal modulation, images with and without Fermi filtering are shown as well. The images with acceleration factor of 9 are shown in Fig. S1 in the supplementary material. As can be observed in Fig. 5, by normalizing the signal evolution to a consistent magnitude, APIR4EMC does not suffer from the scaling differences along the echo train among different contrasts and the artifacts are reduced, while for GRAPPA the difference is negligible, though the apparent blurring is reduced. With further Fermi filtering, the excessive Gibbs ringing and the noise are reduced while the valid details of the images are not affected.

Phantom experiments
The RMSEs of reconstructed images by both methods with and without stabilization of signal modulation are shown in Table 2. Compared to GRAPPA, APIR4EMC achieved lower RMSE, which is more obvious in acceleration factor of 9. In APIR4EMC reconstructions, by including stabilization and Fermi filtering, the best overall RMSE was achieved among others, although in two contrasts (T1-Fatsat and PD) with acceleration factor of 9 it was slightly suboptimal. Fig. 6 shows one axial slice of the reconstructed 3D images for both APIR4EMC and GRAPPA using the selected optimal patterns in both subsampling factors of 6 and 9. Images of the coronal and sagittal slices are shown in Fig. S2 and S3, respectively, in the supplementary material. The GRAPPA images in Fig. 6 are with stabilization of signal modulation Fig. 3. One axial slice (the PE1 × PE2 plane) of the reconstructed phantom images (odd columns) of all contrasts by GRAPPA and APIR4EMC and their corresponding g-factor maps (even columns) with the optimized patterns for subsampling factors of 6 and 9. The g-factor intensities range linearly from 0 to 10. The reference images with the full acquisition are shown in the last row. Note that the stripes through the small spheres in (e) are probably due to blurring caused by modulation of the high intensities in the spheres with the severe signal decay along echo train in the PE1 (vertical) direction. and with Fermi filtering as well. As can be seen, the GRAPPA images overall are noisier than the APIR4EMC image with the same subsampling factors. With subsampling factor increased to 9, the GRAPPA images become very noisy and this is substantially improved in the APIR4EMC images. Both reconstructions are with proper contrasts as intended to achieve, although the contrast to noise ratio (CNR) of GRAPPA is lower than that of APIR4EMC. As pointed out by the red arrows in Fig. 6, strong aliasing artifacts in the brain region of GRAPPA reconstruction with acceleration factor of 9 are present and are substantially reduced in APIR4EMC reconstruction with the same acceleration factor. For APIR4EMC, the images with acceleration factor of 6 show more artifacts than with factor of 9. This is expected to be primarily due to motion, especially in the relatively long acquisition of the FLAIR contrast with acceleration factor of 6. Overall in APIR4EMC images, PD and FLAIR images show more artifacts than others. Besides the effects of motion, it may also be partially due to over-regularization since the contrasts are inherently different in SNR and a constant regularization strength might be suboptimal for some of the contrasts.

In-vivo experiments
The reconstructed images for both methods with and without stabilization of signal modulation with acceleration factors of 6 and 9 are shown in Fig. S4 and S5, respectively, in the supplementary material. APIR4EMC with signal stabilization substantially improves image quality over GRAPPA without stabilization, showing much less apparent blurring, less severe artifacts, and lower noise amplification, especially with the high acceleration factor 9.
The leave-one-contrast-out reconstructions with acceleration factors of 6 and 9 are shown in Fig. 7 and the supplementary Fig. S6, respectively. Without the k-space of the reconstructed contrast itself, it is not able to reconstruct the image. However, it becomes clear that the other contrasts contribute (or rather, compensate for) artifact-like image content. Indeed, when leaving out one of the other contrasts, the artifacts become more severe than in the original reconstruction. This observation is consistent between acceleration factors of 6 and 9. With the higher acceleration 9 the artifacts and the noise are overall stronger.

Discussion
This paper presents an autocalibrated parallel imaging reconstruction for extended multi-contrast imaging (APIR4EMC) method to improve image quality for multi-contrast imaging by exploiting intrinsic signal correlation among different contrasts with optimized k-space patterns.
Using the same scan time, APIR4EMC improves the image quality of the multi-contrast imaging compared to GRAPPA reconstruction on a single contrast in terms of both noise amplification and artifacts. Compared to JVC-GRAPPA [24] and APIR4GRASE [23], which reconstruct k-spaces of different echo times or echo types in one echo train, thanks to the stabilization and Fermi filtering, APIR4EMC explores the possibilities to reconstruct more generally different contrasts with separate contrast preparation and different signal evolutions. As Similar to other methods [23,24,37] which essentially employ the virtual coil concept to exploit additional signal information for reconstruction, the computation time of APIR4EMC also increases with more  Table 2 RMSE (×10 − 2 ) for phantom images with acceleration factors of 6 and 9, reconstructed using different methods.
contrasts acting as virtual coils. In the current implementation the time complexity of APIR4EMC with regard to the number of contrasts (c) is O (c 4 ). With an optimized implementation and a slight redefinition of the kernel neighborhood, probably resulting in equal reconstruction quality, the time complexity of APIR4EMC could be reduced to O(c 3 ).
Using the entire echo train for imaging, like in standard FSE, was enabled by stabilizing the signal evolution with the zero-phase encoded echo train. Thanks to this stabilization, artifacts were reduced in APIR4EMC. This is probably due to reduction of bias in the predictions made with the kernel that combines all contrasts and is estimated in the ACS region. The ACS is acquired by only a small part of the echo train. Without compensation, the signal modulations change the correlation between contrasts in other parts of the k-space. With stabilization, artifacts are reduced.
The stabilization could also alter the appearance of the image. Due to the T2 decay, the FSE using a center-out view ordering suppresses the high frequency signals in the periphery of the k-space which are acquired at the end of the echo trains, which would lead to a smoothed appearance of the image. With stabilization, the suppressed high frequency signals are reverted to the same magnitude level as other frequency components, thus avoiding the smoothed appearance. However, the noise is amplified substantially by this procedure as well. While other, more advanced techniques to address this have been proposed [38,39], a straightforward Fermi filtering as used in our work already proved to substantially suppress the noise amplification and improve SNR. Finally, the signal stabilization does not change the image contrast, because with a center-out view ordering the signal stabilization has a minor effect in the center of k-space (which has the largest impact on the perceived contrast).
In this paper, FSE was used to acquire multi-contrast images. One interesting question would be if APIR4EMC also works with other sequences, e.g., gradient echo based sequences which are routinely used in clinical practice as well. With gradient echo sequences, the echo is formed with the rephasing gradient instead of the refocusing RF pulse, in which case dephasing due to magnetic field inhomogeneities, tissue susceptibility, or chemical shifts are not rephased in the gradient echo. The phase evolution therefore is substantially different from the one with spin echo sequences, due to which APIR4EMC may need careful further investigation for working with both sequences.
A limitation of the proposed approach is that the separate scans of multiple contrasts are vulnerable to motion which can lower the image quality of the reconstruction. Any displacement during the entire scan time of all contrasts can lead to motion artifacts. To compensate the motion, the proposed method can potentially be extended with a retrospective motion compensation method as in [40] which compensates motion by maximizing the signal correlation among coils in kspace using the ACS region. Another factor that can cause suboptimal reconstruction of multiple contrasts is the regularization strength applied in the convolution kernel estimation. With different acquisition protocols, different contrasts are inherently different in SNR and may require different optimal regularization strength. As such the constant regularization strength used in this work may be suboptimal for some of As an extension to GRAPPA, with virtual coils integrated in the reconstruction, the current APIR4EMC framework is still a linear method using least squares fitting to compute the prediction kernel coefficients and subsequent (linear) convolution to predict the unsampled k-space positions. However, the signals among the different contrasts might be related in non-linear ways. In the current approach such nonlinear relations cannot be recovered. This may put a limit on the acceleration capability of the proposed method. Methods with the ability to approximate non-linear relations, such as APIR-Net [41] or RAKI [42], may improve the image quality at high(er) acceleration factors at the expense of a more complicated and potentially less stable reconstruction performance.

Conclusions
By including multi-contrast data into the reconstruction process, APIR4EMC is able to reduce artifacts and noise amplification for each contrast compared to GRAPPA with the same subsampling factor. Thus, APIR4EMC enables increased acceleration capability for multi-contrast acquisitions and this was validated by in-vivo multi-contrast brain imaging.  Fig. 7. The in-vivo images of all contrasts with acceleration of 6 reconstructed by APIR4EMC where in each row the sampled k-space data of one contrast was set to zero when computing unsampled k-space positions of all contrasts. These images should be compared to the last row of Fig. 6, where none of the data was set to zero. For each reconstructed image, one axial slice (the PE1×PE2 plane) is shown.

Author statement
Curation,Writing -Reviewing and Editing,Supervision.