High‐dimensionality undersampled patch‐based reconstruction (HD‐PROST) for accelerated multi‐contrast MRI

Purpose To develop a new high‐dimensionality undersampled patch‐based reconstruction (HD‐PROST) for highly accelerated 2D and 3D multi‐contrast MRI. Methods HD‐PROST jointly reconstructs multi‐contrast MR images by exploiting the highly redundant information, on a local and non‐local scale, and the strong correlation shared between the multiple contrast images. This is achieved by enforcing multi‐dimensional low‐rank in the undersampled images. 2D magnetic resonance fingerprinting (MRF) phantom and in vivo brain acquisitions were performed to evaluate the performance of HD‐PROST for highly accelerated simultaneous T1 and T2 mapping. Additional in vivo experiments for reconstructing multiple undersampled 3D magnetization transfer (MT)‐weighted images were conducted to illustrate the impact of HD‐PROST for high‐resolution multi‐contrast 3D imaging. Results In the 2D MRF phantom study, HD‐PROST provided accurate and precise estimation of the T1 and T2 values in comparison to gold standard spin echo acquisitions. HD‐PROST achieved good quality maps for the in vivo 2D MRF experiments in comparison to conventional low‐rank inversion reconstruction. T1 and T2 values of white matter and gray matter were in good agreement with those reported in the literature for MRF acquisitions with reduced number of time point images (500 time point images, ~2.5 s scan time). For in vivo MT‐weighted 3D acquisitions (6 different contrasts), HD‐PROST achieved similar image quality than the fully sampled reference image for an undersampling factor of 6.5‐fold. Conclusion HD‐PROST enables multi‐contrast 2D and 3D MR images in a short acquisition time without compromising image quality. Ultimately, this technique may increase the potential of conventional parameter mapping.


| INTRODUCTION
In MRI, multiple contrasts are exploited to extract clinically relevant tissue parameters and pathological tissue changes. These multiple contrasts are achieved using different imaging sequences and preparation pulses. Multi-contrast acquisitions also find important applications in parameter mapping (e.g., T 1 and T 2 mapping) and magnetic resonance fingerprinting (MRF). 1,2 However, these acquisitions lead to long scan times because multiple images with different contrasts need to be acquired, making parameter imaging more sensitive to physiological motion. [3][4][5][6] Parallel imaging (PI), [7][8][9][10][11] compressed sensing (CS), 12,13 as well as the combination of both undersampled reconstruction techniques 14,15 have been proposed to overcome the long scan times associated with multi-contrast imaging and parameter mapping. PI can accelerate multi-contrast imaging by undersampling each individual image and exploiting the information provided by multiple coil arrays, yet at a SNR penalty generally marked for high acceleration factors. Sparse CS alone has been shown to cope with the problem of undersampling through the use of random or pseudo-random sampling patterns and efficient regularized reconstructions that make the assumption that the multi-contrast images share common and sparse information in a specific domain. [16][17][18][19][20][21] Even though these strategies have achieved acceleration factors that have not previously been possible to attain with parallel imaging alone, CS-based techniques still suffer from residual aliasing artifacts for high acceleration factors, which compromise the diagnostic value of the reconstructed multi-contrast images.
Recently, novel techniques that exploit the strong anatomical correlations observed in the contrast dimension (or parameter dimension) on a global or local scale have been proposed. Indeed, the nature of signal evolution in multicontrast acquisitions exhibits a low-rank structure in the contrast dimension that can be exploited to further reduce scan times. 17,[22][23][24] These types of reconstruction techniques, also known as the globally (GLR) or locally low-rank (LLR) methods, 25 have been efficiently used in many applications such as T 2 mapping 26 or dynamic contrast enhanced MRI. 27 More recently, high-order tensor decomposition techniques, exploiting global correlation, have been efficiently used to allow for highly accelerated multi-dimensional cardiac MRI acquisitions. 28,29 Although those techniques have shown promise for motion-resolved quantitative cardiac imaging by efficiently solving a global low-rank tensor decomposition, they do not exploit the strong non-local correlations between neighboring patches.
Motivated by the LLR techniques that exploit localized correlations in the contrast dimension, patch-based image reconstructions exploiting non-local spatial redundancies and low-rank matrix structures have been introduced for single-contrast MRI reconstruction to lead to even sparser representation. 30,31 By modeling the similarity of image patches through block-matching, low-rank representation and filtering, 2D, 32 and 3D 33 patch-based reconstructions have been shown to outperform conventional CS reconstructions by recovering better image details and edges and exhibiting better overall image quality.
In this study, we present a new reconstruction technique for highly accelerated 2D and 3D multi-channel multicontrast MRI that combines the promising performances of patch-based reconstructions and the potential of low-rank image reconstruction through higher-order tensor decomposition. The proposed high-dimensionality undersampled patch-based reconstruction (HD-PROST) technique is first applied to accelerated 2D radial MRF, for various acceleration factors, where a high degree of inherent redundancy can be exploited locally, non-locally, and through the contrast dimension. In a second application, HD-PROST is used to acquire multiple undersampled high-resolution 3D Cartesian magnetization transfer contrast (MTC) images with several MT weightings in a reduced scan time.

| THEORY
The framework presented hereafter jointly reconstructs multi-channel multi-contrast images from undersampled 2D or 3D MR acquisitions. This is achieved by generalizing our previously proposed PROST technique 33 to high dimensional imaging. A description of the proposed HD-PROST reconstruction is presented, followed by the description of 2 multicontrast applications (2D radial and 3D Cartesian) where high-dimensionality can be exploited to reduce acquisition time, which is often a key factor for clinical translation.

| High-dimensionality undersampled patch-based reconstruction (HD-PROST)
Let X ∈ ℂ M x × M y × M z × L be the multi-contrast complex images that we seek to reconstruct, where M x , M y and M z are the number of voxels in the x, y and z spatial directions, and L is the number of contrast-weighted images. The corresponding complex receive-coil sensitivity maps for the N c channels are denoted as The joint multi-contrast undersampled reconstruction can be combined with parallel imaging and cast as the following inverse problem where A is the undersampling operator that acquires k-space data for each contrast-weighted image, F denotes the Fourier transform operator and ‖ ⋅ ‖ F is the Frobenius norm.
(1) argmin Mathematically, this inverse problem is ill-posed, in the sense that the exact solution might not exist or not be unique, making precise recovery of X hardly possible, and prior assumptions on the unknown solution X have to be considered. The principle behind HD-PROST reconstruction assumes that a multi-contrast image X can be expressed as a high-order low-rank representation on a patch scale, with respect to an appropriately chosen patch selection operator. The recovery problem can be formulated as the following constrained optimization on the high-order low-rank tensor  where p is the nonnegative sparsity-promoting regularization parameter and ‖ ⋅ ‖ * is the nuclear norm that enforces multi-dimensional low-rank on a multi-contrast patch scale. The patch selection operator P p (⋅) forms a 3D tensor from a patch centered at pixel p from a set of multi-contrast images (see optimization 2 below). Now considering the constraint  p = P p (X), and the encoding operator E = AFS, we can form the unconstrained Lagrangian of Equation 2 by linearly combining the constraint and cost function 31,33 where b is the Lagrange multiplier, and > 0 is the penalty parameter. Equation 3 can be efficiently solved through operator-splitting via alternating direction method of multipliers (ADMM). 34 ADMM simplifies the optimization process by alternating the minimization with respect to the multicontrast set of images X (optimization 1) and the high-order tensor  (optimization 2) followed by an update of the augmented multiplier b, and repeating these 3 steps until a convergence criterion is satisfied.

| Optimization 1: joint MR reconstruction update
The first sub-problem is a joint multi-contrast MR reconstruction that incorporates the denoised tensor  (obtained at the end of optimization 2) as prior information in a parallel imaging fashion to obtain X Equation 4 corresponds to a standard iterative SENSE reconstruction with Tikhonov regularization, where the solution X can be efficiently computed using the Conjugate Gradient 35 algorithm.

| Optimization 2: high order singular value decomposition (HOSVD)-based denoising
Considering the variable  p = P p (X) + b p , the second subproblem minimizes with respect to the high-order tensor  and is given by X denotes multiple MR images with different contrasts. Several observations can be made about X: (1) on a local scale, voxels at a specific location for a given contrast exhibit similar intensity to their nearest neighbors (within a patch), (2) on a non-local scale, images for a given contrast contain selfrepeating patterns (measured as patch similarity within a neighborhood), and (3) on a contrast scale, common structures and features are shared across multiple contrast images. Motivated by these observations, the proposed joint multi-channel multicontrast problem can be cast as a multi-dimensional low-rank reconstruction. Bearing this in mind, equation 5 can be solved on a multi-contrast patch level. The construction of the highorder tensor  is performed as a 3-step process: Step 1 -Similar overlapping patches in X + b are grouped together to form a third-order tensor: considering a 3D + L reference patch of size N x × N y × N z × L, we build a high dimensional tensor  p ∈ ℂ N × K × L of K − 1 similar 3D + L patches, with N = N x × N y × N z (see Figure 1, "unfolding" and "tensor stacking"). A fixed local window is used for the patch search, whereas the contrast signature remains unchanged. Along this line, the proposed reconstruction can exploit as much of the contrast and spatial correlations as possible.
Step 2 -The tensor  p exhibits a strong low multilinear rank structure and can therefore be compressed into a tensor of smaller size (i.e., the core tensor) through tensor decomposition (see Supporting Information Table S1 and Figure 1, "High-Order Tensor Decomposition"). The dominant components of the core tensor can be extracted by computing a complex-valued higher-order singular value decomposition (HOSVD) 36,37 and by only keeping the largest (given by the thresholding parameter 2 p ) multilinear singular vectors and high-order singular values. This step effectively acts as a high-order denoising process where the small discarded coefficients mainly reflect contributions from noise and noise-like artifacts.
Step 3 -The denoised tensor  p is then rearranged to form the denoised patches. Steps 1-3 are repeated over all patches in the image in a sliding window fashion. Because a single patch might belong to several groups in step 1, the final denoised multi-contrast complex-valued images  are obtained by averaging ( Figure 1, "Aggregation") the different estimates.
The solution  to this optimization problem is a denoised version of  that is incorporated in the optimization 1 as prior knowledge, as described before. The Lagrangian multiplier b is then updated and optimizations 1 and 2 are processed iteratively to improve the quality of the reconstructed images. In the spirit of reproducible research, codes and examples for the proposed HD-PROST technique are made available at http://www.kclcardiacmr.com/downloads/.
The generalized reconstruction framework described before considers 2D or 3D Cartesian multi-contrast acquisitions (as the 3D undersampled Cartesian multi MT-weighted acquisitions considered in this study). Slight modifications in the reconstruction process are required for the accelerated non-Cartesian 2D MRF application considered in this study and will be described in the next section.

| HD-PROST for accelerated 2D radial parameter mapping with MRF
MRF 1 is a novel quantitative MRI approach that allows the simultaneous acquisition of multi-parametric maps (e.g., T 1 , T 2 , M 0 ) in a single efficient scan. Conventional MRF sequences acquire in the order of thousands or more highly undersampled time point images by pseudo-randomly collecting the MR data in a continuous fashion with timevarying acquisition parameters (e.g., repetition time, flip angle). The spatial and temporal incoherencies provide a unique signal evolution (or fingerprint) for each tissue. These unique fingerprints can be matched, through pattern matching, to a pre-generated MRF dictionary representative of the MRF sequence, and whose atoms are composed of simulated signal evolution curves. This matching process is performed on a voxel-by-voxel basis to identify the underlying tissue properties and generate quantitative parameter maps. The highly undersampled pseudo-random MRF acquisition results in a high level of noise and aliasing in the reconstructed time point images. Several iterative techniques have been recently proposed to improve the reconstruction quality of each time point image. [38][39][40][41][42] Zhao et al 38 proposed to enforce low-rank and subspace modeling in the temporal dimension to reconstruct high-quality time point images. Assländer et al 39 recently introduced a low-rank ADMM reconstruction technique to temporally compress the time point images, resulting in a reduced number of singular value images. The reconstruction of the temporally compressed images is faster and better posed than reconstructing each time point image separately. 39 This temporal compression operator U r is obtained through compression of the MRF dictionary at an appropriate rank r. F I G U R E 1 Flowchart of the optimization 2 of the proposed high-dimensionality patch-based reconstruction (HD-PROST). Denoising of multi-contrast images is performed using 2D (respectively 3D) block matching, which groups similar 2D (respectively 3D) patches in the multicontrast images. Similar patches are then unfolded together in a simple 2D matrix. A third-order tensor  p is formed by stacking the unfolded patches in the contrast dimension. The high-order tensor of size N × K × L admits a low multilinear rank approximation and can be compressed, through tensor decomposition, by truncating the multilinear singular vectors that correspond to small multilinear singular values. The outputs of this step are the denoised multi-contrast images that are then used in the joint MR reconstruction process (optimization 1) as prior knowledge. An overview of the algorithm is provided in Supporting Information Table S1 | 3709 BUSTIN eT al.
Because of the multi-contrast nature of MRF, HD-PROST can be used to explicitly exploit the local, non-local and contrast information of the temporally compressed images by integrating the compression operator into the encoding operator in Equation 3 as follows

| METHODS
The proposed HD-PROST reconstruction was evaluated on accelerated radial 2D MRF phantom and in vivo brain acquisitions, and on accelerated Cartesian 3D magnetization transfer imaging with varying MT-weighting in in vivo brain data. The 2 applications are described in detail below along with imaging and reconstruction parameters. Written informed consent was obtained from all subjects before undergoing MRI scans and the study was approved by the Institutional Review Board.

| Accelerated 2D magnetic resonance
fingerprinting MRF acquisitions were performed on a 1.5 T Ingenia MR system (Philips, Best, the Netherlands) equipped with a 15-element head coil.

| Phantom and in vivo experiments
A 2D MRF acquisition was performed on a standardized (T1MES) T 1 /T 2 phantom containing 9 agarose-based tubes with different T 1 and T 2 combinations (range, T 1 : 255 ms to 1489 ms, T 2 : 44 ms to 243 ms). 43 Relevant scan parameters included: balanced steady-state free precession radial sequence, TE = 2 ms, fixed TR = 4.4 ms, FOV = 160 × 160 mm 2 , in-plane resolution = 1 × 1 mm 2 , slice thickness = 8 mm, bandwidth = 723.4 Hz/pixel. Only 1 radial spoke was acquired at each time point (resulting in an acceleration factor of ~251 with respect to a fully sampled radial acquisition). A total of 2000 time points were acquired in 10 s. A flip angle (FA) pattern similar to the one proposed in Assländer et al 44 for optimized T 1 /T 2 mapping was used and is shown in Supporting Information Figure S1. This RF pattern, which has been shown to be optimal in a Cramér-Rao lower bound sense, consists of intrinsic repetitive loops that offers the advantage to lengthen the scan time by simple concatenation. The experiments consisted of undersampling the acquired data by keeping only Reference T 1 and T 2 times for each vial were obtained from gold standard spin echo (SE) acquisitions. For T 1 values, an inversion-recovery SE (IRSE) sequence was used with 8 inversion times from 25 ms to 3200 ms with TR = 10s, TE = 14.75 ms. For T 2 values, the SE sequence was performed with 8 TEs from 10 ms to 640 ms. T 1 and T 2 values were obtained by mono-exponential curve fitting.
Single-slice 2D MRF brain data were acquired in 5 healthy subjects (4 men, mean age: 32 years; range: 28-37 years) using the same scan parameters as in the phantom experiments.

| Image reconstruction
For both phantom and in vivo 2D MRF experiments, data was temporally compressed with r = 10, leading to only 10 singular value images to reconstruct (i.e., in this study, L = 10 and M z = 1).
HD-PROST reconstruction was implemented using the algorithm described in Supporting Information Table S2 and performed offline on a workstation with a 16-core Dual Intel Xeon Processor (23 GHz, 256 GB RAM). The joint MR reconstruction step (optimization 1) was implemented in MATLAB (v7.1, The MathWorks, Natick, MA) and the multi-contrast patch-based denoising step (optimization 2) in C. Coil sensitivity maps were estimated using the eigenvaluebased approach ESPIRiT. 45 The encoding operator E MRF was implemented using the nonuniform fast Fourier transform. 46 The tolerance of the conjugate gradient was set to CG eps = 1e −4 and a maximum number of CG iter = 15 iterations was chosen as stopping criterion. The regularization parameter , which balances the contribution of the prior term (obtained at the end of optimization 2) and the data fidelity term, was set to 5e −3 .
The proposed high-order patch-based denoising strategy was implemented as described in Supporting Information Table S1. The performance of the proposed strategy relies on the optimal selection of several parameters. The patch size, which controls the degree of local image features, was set to N = 7 × 7. We set the search window radius around each pixel to 20 and restricted the number of similar patches selected to K = 20 to form a third-order tensor  p of size 49 × 20 × 10. The l 2 distance was chosen as measure of patch similarity and was defined as To save computational complexity, a sliding-window approach was performed with a patch offset of 3 pixels at each image dimension. The performance of HD-PROST was assessed on several data sets (not reported here) by comparing the quality of the reconstructions with several regularization parameters (the same was used for all patches: p = for all p). The optimal value was shown to be proportional to the number of MRF measurements and was set to = −1e −3 × n + 0.4 for each decomposition, with (6) E MRF = AU r FS n being the number of MRF radial spokes. The joint MR reconstruction and denoising steps were iteratively interleaved and the reconstruction was terminated after 5 ADMM iterations.
The proposed HD-PROST reconstruction for 2D MRF was compared to the low-rank inversion (LRI) reconstruction 24,38 with r = 10 and using 10 conjugate gradient iterations, which were seen to be enough for convergence.

| Dictionary generation and pattern recognition
The MRF dictionary was generated using the extended phase graphs (EPG) formalism. 47 The dictionary was calculated for a T 1 in the range of ([50:10:1400, 1430:30:1600, 1700:100:2200, 2400:200:3000] ms) and T 2 in the range of ([5:2:80, 85:5:150, 160:10:300, 330:30:600] ms). Slice profile was simulated for each RF pulse using 51 isochromats distributed along the slice selection direction and was included in the dictionary generation to correct for profile imperfections. 48 Template-matching between fingerprints and dictionary were performed using the inner product as in Ma et al. 1

| Acquisition
A 3D accelerated MTC experiment was performed to evaluate the proposed HD-PROST reconstruction on 3D Cartesian acquisitions with multiple MT-weighted images. In vivo brain acquisitions were performed on 3 healthy subjects (1 man, age range: 24-30 y) on a 1.5 T MR scanner (Magnetom Aera, Siemens Healthcare, Erlangen, Germany) equipped with a 20-channel head coil. Acquisitions consisted of 1 reference image without magnetization preparation and 5 images with different MT preparations (i.e., in this study, L = 6 and M z > 1).
A prototype 3D Cartesian variable-density trajectory was integrated in the sequence to allow for fast acquisition of multiple MT-weighted images. The Cartesian trajectory with spiral profile order 33,49 samples the k y -k z phaseencoding plane following approximate spiral interleaves on the Cartesian grid with variable density along each spiral arm and with 2 successive spiral interleaves being rotated by the golden ratio. A golden angle rotation between different contrast acquisitions was incorporated here (shifted VD-CASPR) to introduce incoherently distributed aliasing artifacts along the contrast dimension and noise-like artifacts in the spatial dimension, which is beneficial from a CS and low-rank point of view. 50 The MT weighting was achieved by applying a train of sinc-shaped, off-resonance RF pulses before image acquisition with the following parameters: MT off-resonance frequency (ΔF) = 3 kHz, 20 MT pulse repetitions, MT bandwidth = 401 Hz/pixel. Relevant scan parameters included: 3D gradient echo sequence, axial orientation, FOV = 230 × 230 × 160 mm 3 , nominal resolution 1 × 1 × 2 mm 3 , FA = 15°, TE = 1.78 ms, TR = 4.06 ms, receiver bandwidth = 925 Hz/pixel, 32 readouts per spiral interleave. Six measurements were acquired with different MT pulse flip angles ( MT = [0 • , 160 • , 320 • , 480 • , 640 • , 800 • ]) with a 5-s pause between them. Acquisitions were performed with an acceleration factor of 6.5-fold for each weighted image. The total scan time to acquire the 6 measurements was 13:18 [min:s]. A fully sampled acquisition of the 6 measurements at this resolution would take more than 1 h. Therefore, for comparison purposes, an additional fully sampled acquisition was performed only for the reference image ( MT = 0 • ). The total scan time for this single-contrast fully sampled acquisition was 12:57 [min:s].
The proposed HD-PROST reconstruction was compared with 2 well-established state-of-the-art reconstruction techniques. The first technique is LLR, proposed by T. Zhang et al 26 for accelerating MR parameter mapping. LLR exploits the redundancy in the contrast dimension on local image regions in an iterative low-rank framework. LLR was implemented using our ADMM framework by replacing the patch-based denoising step by the low-rank thresholding. This allows for fair comparisons because the same optimization was used, and only the manner in which the denoising is performed was modified. The rank threshold LLR was fixed and set to 5% of the highest singular value. Because the acquired MT-weighted data was fully sampled in the readout direction, the MR reconstruction step was accelerated for both LLR and HD-PROST reconstructions by computing a 1D inverse FFT and considering multiple separable 2D reconstruction problems independently.
The second technique is an iterative CS reconstruction with spatial Wavelet sparsity constraint as described in Lustig et al 12 Figure 2 shows T 1 and T 2 values for the 2D MRF phantom experiments with 2000, 1000, and 500 time points in comparison to the gold standard IRSE and SE acquisitions for both LRI and HD-PROST reconstructions. T 1 values obtained from both strategies were in good agreement with the IRSE acquisition even for reconstructions with 500 time points, with an excellent linear relationship with the reference T 1 values (goodness-to-fit R 2 > 0.98). T 2 accuracy was also preserved with the proposed reconstruction with a slight T 2 degradation observed for long T 2 values and high acceleration for both reconstructions. Figure 3 depicts the precision of T 1 and T 2 values, as characterized by the SD (aggregated based on the variance of each vial). An increase in precision was observed for both T 1 /T 2 values using the proposed HD-PROST reconstruction compared with LRI even for reconstructions with 500 time points, corresponding to 2.5s scan time. Corresponding T 1 and T 2 maps are shown in Supporting Information Figure S2. From the above analysis, it follows that 500 MRF time points or less might be sufficient and suitable for accurate and precise in vivo T 1 /T 2 maps acquisitions in <2.5 s. Figure 4 depicts the first 4 2D MRF singular images from the reference LRI and the proposed HD-PROST reconstruction for 1 representative subject reconstructed with 1000 time points. A clear superior image quality can be observed on the HD-PROST singular images with a sharp and clear delineation of the brain structures. A high level of streaking artifacts and noise can be seen on the last singular value components (e.g., singular images 3 and 4) with LRI, whereas HD-PROST not only produces images with considerably less noise but is also able to recover small structures that were lost below the noise level with LRI ( Figure 4, yellow arrows). T 1 and T 2 maps are displayed in Figures 5 and 6 for 2 subjects and 3 different measurement lengths (2000, 1000, and 500 time points) for both LRI and HD-PROST reconstructions.

| In vivo study
The reconstructed maps from 1 additional subject are shown in Supporting Information Figure S3. A number of F I G U R E 2 Phantom results for the 2D accelerated MRF using low-rank inversion (LRI) and the proposed HD-PROST reconstructions. Plots are comparing the mean T 1 (A) and T 2 (B) values derived from 2000, 1000, and 500 time points, with conventional inversion-recovery spin-echo (IRSE) and spin-echo (SE) acquisitions (identity lines). T 1 and T 2 accuracies are preserved with the 2 strategies, with a slight bias observed for long T 2 s at high accelerations for both methods. The mean values were obtained from ROIs drawn around each phantom vial. Abbreviations: LRI, lowrank inversion; HD-PROST, high-dimensionality undersampled patch-based reconstruction interesting observations can be made. Reducing the number of measurements tends to blur the T 1 maps with LRI whereas the T 2 maps suffer from noise amplification, showing an overall noisier appearance. Conversely, by enforcing low-rank in the local, non-local and contrast dimension, HD-PROST reconstruction delivers higher image quality, recovering sharpness for T 1 and reducing the noise for T 2 . The improvement is more pronounced for the 500 time points acquisition (2.5 s scan time). In vivo T 1 and T 2 relaxation times measured in regions of interest in the white and gray matters with LRI and the proposed HD-PROST are shown in Table 1. Both reconstructions converged to very comparable values that are in good agreement with values obtained from the literature for T 1 . Moreover, the proposed HD-PROST reconstruction tends to lower the SDs of T 1 and T 2 times, which is in accordance with the noise reduction seen in the quantitative

F I G U R E 4 Reconstructed first 4 MRF singular images with low-rank inversion (LRI) (A) and the proposed HD-PROST (B) in in vivo brain
experiments in a representative subject acquired with 1000 time points. A clear improvement in image quality and image sharpness can be observed on the HD-PROST reconstruction with considerable reduction of noise and streaking artifacts, particularly for the last singular images maps. Note that the T 2 relaxation times for both techniques are slightly biased and depart from the literature values. This may be partly explained by the fact that B 1 imperfections 52 as well as other sources of bias such as magnetization transfer 53 and diffusion-weighting 54 were not considered in the proposed study. The average reconstruction time for 2D MRF with HD-PROST was ~10 min per data set. Additional comparisons with single-contrast PROST reconstruction (i.e., reconstructing each singular image independently) and with a global low-rank tensor decomposition (in the spirit of cardiac multitasking) 28,29 are provided in Supporting Information Figure S4 Figure 7C, showing excellent agreement between HD-PROST and the fully sampled reference. Six different undersampled MT-weighted images were acquired in 13 min 18s, whereas the fully sampled acquisition of a single contrast took 12 min 57 s. Figure 8 compares HD-PROST to conventional CS reconstruction from a 6.5-fold acceleration. Comparisons with zero-filling and LLR reconstructions are provided in Supporting Information Figures S5 and S6. As expected, zero-filling exhibits a low image quality with apparent aliasing artifacts and blurring. Exploiting contrast redundancy through local image regions with LLR improves the overall image quality and enables the recovery of small structures, particularly for low-contrast images (e.g., MT = 800 • ), while the apparent noise is still large. By contrast, CS reconstruction with spatial regularization is able to recover images with reduced level of noise but fails to recover small structures for low contrast images  Table 1) can be seen in Supporting Information Figures S7 and S8. The average computation time for 3D HD-PROST reconstruction was ~27 min for all 6 contrasts in the acquisitions performed in this study.

| DISCUSSION
HD-PROST reconstruction enables accelerated acquisition of 2D or 3D multi-contrast MR images by exploiting the high local and non-local redundancies and the similarities between the multi-contrast images through a high-order low-rank tensor approximation.
The proposed technique was applied to accelerated non-Cartesian 2D MRF and accelerated Cartesian 3D MTC imaging to enable undersampling factors that go beyond the limit of traditional PI and CS reconstructions (i.e., ~2.5 s acquisition for 2D MRF and 6.5-fold acceleration for 3D MTC), while removing residual aliasing artifacts. Phantom experiments in accelerated 2D MRF were carried out to investigate the impact of rapid acquisition (i.e., reduced number of time point images) on accuracy and precision of T 1 and T 2 relaxation times. High agreement with reference T 1 /T 2 values was observed using HD-PROST, even for high accelerations, with increased precision compared to conventional LRI reconstruction.
For in vivo MRF, streaking artifacts and noise amplification often propagated in the T 1 maps with LRI reconstruction, whereas blurring was observed on the T 2 maps for high acceleration factors. HD-PROST achieved improved sharpness and F I G U R E 7 Three-dimensional reconstruction of a MT-weighted 6.5-fold undersampled brain data in a healthy subject (subject 1). HD-PROST reconstruction (B) is compared to the fully sampled acquisition (A) for the reference image only ( MT = 0 • ). Line profiles going through a structure with sharp edges are shown in (C). HD-PROST is able to recover high fidelity 3D images and retrieve sharp edges in agreement with the fully sampled acquisition. Six different undersampled MT-weighted images were acquired in 13 min 18 s, whereas the fully sampled acquisition of a single contrast took 12 min 57 s F I G U R E 8 6.5-fold accelerated 3D MT-weighted images for 6 different contrasts from 1 representative subject (subject 1) reconstructed with compressed-sensing (CS) and the proposed HD-PROST reconstruction. Fine anatomic structures can be efficiently retrieved with HD-PROST as shown by the arrows. See Supporting Information Figure S5 for the visualization of the whole axial images and Supporting Information Figure S6 for comparisons with zero-filling and locally low-rank reconstructions reduced noise level in comparison to the low-rank inversion reconstruction, especially for acquisitions with reduced number of time points. Nevertheless, a systemic underestimation of the T 2 values, previously reported in MRF literature, was observed in the in vivo study. This finding may be partly explained by the fact that B 1 imperfection, 52 magnetization transfer, 53 and diffusion-weighting 54 were not considered in this MRF study and could lead to inaccurate T 2 measurements.
HD-PROST has a modular design, which allows for its straightforward extension to 3D or n-D imaging by simple patch vectorization. In line with the previous 2D MRF study, accelerated 3D MTC using HD-PROST showed improved image quality over conventional CS and low-rank reconstructions for an acceleration factor of 6.5, with visual quality comparable to the fully sampled acquisition. High denoising performance was achieved because of the existence of multiple MT-weighted images of the same object with varying contrasts, leading to high redundancy that can be exploited by HD-PROST. The pseudo-random sampling, given by the proposed shifted VD-CASPR, causes aliasing artifacts that spread incoherently in the contrast dimension and exhibits noise-like perturbations at the image scale, providing an excellent basis for HD-PROST reconstruction. This study was only performed on a small number of subjects and further evaluations on larger cohorts are needed. Nevertheless, this proof of concept suggests an opportunity for high-resolution quantitative magnetization transfer imaging in a clinically feasible scan time.
The efficient multithreaded implementation of the highorder patch-based denoising allowed for fast image denoising of large data sets (e.g., in the order of 200 s for a 3D data set with a matrix size of 200 × 256 × 104 × 6). Further speedups could be achieved to reach clinically acceptable runtimes by implementing the joint MR optimization step on multiple GPUs 55 and using coil compression algorithms. 56 HD-PROST imposes low-rank in the complex domain and therefore captures the possible cross-correlation observed between the real and imaginary components, allowing for accurate and faithful reconstruction of both phase and magnitude. Our framework makes use of ADMM to decouple the main optimization problem into 2 simpler sub-problems that have straightforward solutions. Although most of the noise and undersampling artifacts can be efficiently removed after the first iteration, aliasing may still exist depending on the quality of the input images. This behavior mainly stems from the fact that corrupted images can negatively affect the block matching step, resulting in a sub-optimal grouping. Therefore, several ADMM iterations (5 in this study) are needed to achieve good image quality reconstructions (see Supporting Information Figure S9).
The technique proposed in this article can potentially change conventional multi-contrast imaging by making efficient use of the rich and redundant information available locally and temporally. Two applications were introduced in this study, nonetheless HD-PROST stays generic and should be easily extendable to many MR applications where multiple contrasts are involved, such as conventional T 1 and T 2 mapping, perfusion imaging, 57 4D flow MRI, 58 or low SNR applications such as arterial spin labeling. 59

| CONCLUSION
We present a new framework, termed HD-PROST, for efficient reconstruction of undersampled multi-channel multi-contrast MR images. HD-PROST aims at achieving high image quality by exploiting the high local and nonlocal redundancies, and the similarities between the multicontrast images through a high-dimensionality low-rank tensor decomposition. HD-PROST was validated in accelerated 2D MRF to generate precise T 1 and T 2 maps in ~2.5 s without affecting T 1 /T 2 accuracy. For accelerated multiple 3D MT-weighted acquisitions, HD-PROST can recover high quality images, comparable to a fully sampled acquisition, in a clinically reasonable time frame. The straightforward, yet efficient, application of HD-PROST to 2D and 3D multi-contrast data sets provides several opportunities for future research, particularly in areas where high-dimensionality is likely to increase in importance.