Self‐calibrated subspace reconstruction for multidimensional MR fingerprinting for simultaneous relaxation and diffusion quantification

To propose a new reconstruction method for multidimensional MR fingerprinting (mdMRF) to address shading artifacts caused by physiological motion‐induced measurement errors without navigating or gating.


INTRODUCTION
MR fingerprinting (MRF) 1 is a novel framework for quantitative imaging, offering rapid and simultaneous multi-parameter mapping.Recently, multidimensional MRF (mdMRF), 2 has been developed for simultaneous relaxation and diffusion quantification.The sequence structure of mdMRF consists of multiple segments, each starting with a preparation module (T1, T2, or diffusion preparation), followed by continuous data acquisition, and ending with a wait time.This sequence design makes the signal sensitive to both relaxation and diffusion effects, allowing for simultaneous quantification of relaxation and diffusion.However, like conventional multi-shot (or segmented) diffusion MRI techniques, 3,4 mdMRF faces a challenge due to inter-segment phase variation caused by physiological motion, bulk motion, and eddy currents during diffusion encoding.Furthermore, the utilization of diffusion-preparation (DP) scheme 5,6 for diffusion encoding, in the absence of stabilizer gradients 7,8 to maintain high signal intensity, results in magnitude attenuation, as well as intra-segment phase variation, in addition to the above-mentioned inter-shot phase variation.The magnitude attenuation and phase variation can be viewed as measurement errors that deviate from intended signal evolutions, corrupting the diffusion-weighted images.As these measurement errors are not taken into account in the dictionary, shading artifacts arise in diffusion parametric maps, such as the ADC map.These shading artifacts typically manifest as uneven and/or overestimated ADC values. 2 To address such artifacts, the original implementation of mdMRF employed peripheral pulsation gating 9 as a prospective strategy to mitigate the impact of cardiac pulsation on diffusion encoding, thereby reducing measurement errors. 2 Despite its initial effectiveness, this approach comes with several limitations.First, it imposes restrictions on sequence design flexibility since each diffusion-prepared segment must align precisely with the cardiac cycles.Consequently, the wait times between segments are determined by peripheral pulsation gating, which may not be optimal from a sequence optimization point of view.Second, scan efficiency decreases due to lengthened wait times resulting from skipped peripheral pulsation gating.Third, reconstruction efficiency can be significantly reduced as the dictionary needs to be generated retrospectively after each specific scan.Last, the method's robustness decreases due to the potential occurrence of failed pulsation gating.This is because the optimal setting for peripheral pulsation gating is specific to each individual, and the peripheral pulsation gating does not account for other types of motion, such as cerebrospinal fluid pulsation and bulk motion.Motion compensation 10 can be used as an addition or alternative approach to prospective gating.However, this increases diffusion encoding duration time and leads to a decrease of signal intensity due to T2 decay. 11,12Therefore, a retrospective solution is highly desired to effectively address the shading artifacts in mdMRF, to enable scans without the need for prospective gating.
For conventional multi-shot diffusion MRI techniques, which typically utilize a diffusion-weighting (DW) scheme for diffusion encoding, followed by segmented EPI readout for data acquisition, 13,14 many reconstruction methods have been proposed to retrospectively address the issue of inter-shot phase variation (or inconsistency).These reconstruction methods can be divided into two categories based on the reconstruction models used.The first category is phase-corrected (PC) model-based reconstructions.][17][18] Typically, the explicit phase maps are estimated from low-resolution images acquired by an extra navigator 14,16 or self-navigator reconstructed from the densely sampled k-space center for each shot. 3,15,17,18The second category is low-rank (LR) model-based reconstructions, which is a relaxed and linear model implicitly characterizing inter-shot phase variation and exploring data correlation across all shots. 4,8,19Specifically, in LR subspace reconstruction, a subspace is pre-estimated from extra calibration data. 8While in LR matrix completion, the subspace and coefficient images are jointly estimated during reconstruction iterations. 4,19LR model-based reconstruction is considered more robust than PC model-based reconstruction because the former does not rely on the accurate estimation of explicit phase maps. 20However, all these reconstruction methods cannot be directly applied to mdMRF due to the following reasons: (1) mdMRF data involve not only phase variation but also magnitude attenuation, making PC model-based reconstructions inapplicable; (2) the structured LR method proposed in Ref. 4 and the locally LR (LLR) method 21 used in Ref. 19  rely on the assumption of spatially smooth phase modulations between shots.However, whether this assumption can be extended to mdMRF data involving both magnitude attenuation and phase variation has not been investigated; (3) the LR method proposed in Ref. 4 requires lifting data to construct a Hankel structured matrix followed by matrix completion, which is both computationand memory-consuming.It is particularly unsuitable for mdMRF, which employs time-resolved image reconstruction (i.e., an image is reconstructed for each timepoint associated with each TR); (4) the LR method used in Ref. 8 requires extra real-time calibration data to pre-train a temporal subspace, significantly reducing scan efficiency.
In this work, we propose a new reconstruction method, to retrospectively address shading artifacts in mdMRF scans without the need for prospective gating or navigator.The proposed method reconstructs aliasing-free, timeresolved, and high-resolution images, where the measurement errors are accurately represented in the reconstructed images.To the best of our knowledge, the proposed method is the first reconstruction method for multi-shot and/or multi-segment diffusion MRI data where both magnitude attenuation and phase variation are presented.The proposed method, termed self-calibrated subspace reconstruction, is essentially a data-driven LR model-based reconstruction but without the need for extra calibration data acquisition.It consists of two procedures: self-calibration and subspace reconstruction.In Procedure 1 (self-calibration), we reconstruct aliasing-free, time-resolved, and low-resolution images from a subset of imaging data extracted from the k-space center.The low-resolution images include non-corrupted images, as well as images corrupted by measurement errors.To achieve this, we perform segment-wise matrix completion, that is, matrix completion 22 for each segment (containing 96 images).To distinguish this strategy from global LR 23 and LLR, 21,24 we term it "temporally local matrix completion."Here, global LR is commonly used in dynamic imaging, 23 and LLR represents spatially local low rank.In each segment, the accumulated phase during diffusion encoding resulting from physiological motion, bulk motion, and eddy currents causes and determines the measurement errors within that segment. 7,25Since randomly occurring motions between segments are uncorrelated, the accumulated phases show no correlation between segments, resulting in uncorrelation of the measurement errors between segments.Consequently, this significantly weakens the global low-rank assumption.On the other hand, the measurement errors are caused and determined by the accumulated phase during diffusion encoding and constrained within each diffusion-encoded segment, without affecting the neighboring segments.As a result, the signals in each segment can be approximated by a temporally local low-rank subspace.Therefore, we employ temporally local matrix completion instead of temporally global matrix completion for this procedure.Notably, temporally local matrix completion does not require an assumption of spatially smooth phase modulations between segments, unlike the LR methods used in Refs 4,19.In Procedure 2 (subspace reconstruction), we perform temporally global subspace reconstruction. 26In this step, the temporal subspace is estimated from the low-resolution images reconstructed in the Procedure 1.It allows us to reconstruct aliasing-free, time-resolved, and high-resolution images, where the measurement errors are accurately represented in the reconstructed images, from the highly under-sampled imaging data.We choose temporally global subspace reconstruction for this procedure because it enables exploration of global data sharing, that is, exploring data correlation across all segments.This significantly enhances the quality of the reconstructed images, particularly for the SNR.After the reconstruction process, we apply a customized outlier detection algorithm to automatically detect and remove the corrupted segments caused by measurement errors.By removing those corrupted segments in the mapping step (pattern matching), we can generate artifact-free T1, T2, and ADC maps simultaneously.

Self-calibrated subspace reconstruction
A sequence diagram of the mdMRF is illustrated in Figure 1 in Ref. 2 For data acquisition, one k-space readout is acquired for each time point (associated with each TR), resulting in a highly under-sampled k-space.Spiral trajectories with golden-angle rotation across time points are employed for k-space sampling, creating an incoherent sampling pattern.
Figure 1 shows a flowchart of the proposed reconstruction method, consisting of two procedures split into five total steps.Procedure 1 (Steps 1-3) reconstructs aliasing-free, time-resolved, and low-resolution images.A set of data (d c ) is extracted from the k-space center (e.g., 2.5-fold resolution reduction relative to the original resolution) of the under-sampled imaging data (d u ) (Step 1), and then a sliding window (e.g., size = 5, temporal resolution = 30 ms) is used to form a denser under-sampled central k-space ( dc ) (Step 2).Temporally local matrix completion, which performs matrix completion for each segment separately, is used to reconstruct aliasing-free and low-resolution images ( mc ) from dc (Step 3).Next, Procedure 2 (Steps 4 and 5) reconstructs aliasing-free and high-resolution images.A temporally global subspace (V) is estimated by performing singular value decomposition (SVD) and truncation (e.g., within 1% approximation error) on the low-resolution images mc in k-t domain (Step 4).Finally, temporally global subspace reconstruction is used to reconstruct high-resolution images ( m), from highly under-sampled imaging data d u (Step 5).
In Procedure 1, the reconstruction model of temporally local matrix completion is: Flowchart of proposed self-calibrated reconstruction method, consisting of two mainprocedures split into five steps.Procedure 1 (self-calibration) reconstructs aliasing-free, time-resolved, and low-resolution images using temporally local matrix completion from a subset of imaging data extracted from the k-space center (Steps 1-3).Procedure 2 (subspace reconstruction) reconstructs aliasing-free, time-resolved, and high-resolution images using temporally global subspace reconstruction (Steps 4 and 5), where the temporal subspace is estimated from the low-resolution image reconstructed in Procedure 1.
where S c is the low-resolution coil sensitivity estimated from fully sampled data by combining dc along the time dimension, as in Refs 2,27, and F and  denote Fourier encoding in the spatial domain and under-sampling mask in k-t domain, respectively.m c is the low-resolution image series (x-t domain) to be reconstructed, and m c,s is a portion of m c corresponding to the s-th segment.N s is the number of segments and  l is the regularization parameter.|| ⋅ || * is the nuclear norm of a matrix, which is defined as the sum of all singular values of the matrix and used as a convex relaxation of the number of non-zero eigenvalues (i.e., rank), to enforce low rank of the matrix.Fast iterative shrinkage thresholding (FISTA) algorithm 28 is used to solve the problem.The reconstructed low-resolution images mc is used for estimation of the temporally global subspace V by performing SVD on mc in the following step.
In Procedure 2, the reconstruction model of temporally global subspace reconstruction is: where S is the high-resolution coil sensitivity estimated from fully sampled data by combining d u along the time dimension, as in Refs 2,27.U denotes the coefficient images to be reconstructed, that is, image series in spatial-SVD compressed temporal domain. 29Aliasing-free and high-resolution images m (x-t domain) can be finally formed by ÛV H .  denotes the additional sparsity constraint and  u is regularization parameter.Nonlinear conjugate gradient (NLCG) algorithm 30 was used to solve the problem.

Measurement errors detection and removal
Using the proposed reconstruction method, aliasing-free, high-resolution, and time-resolved images, including those corrupted by measurement errors (involving magnitude attenuation and phase variation), can be reconstructed.These corrupted segments, that is, the diffusion-encoded segments containing corrupted images, can be automatically detected using a customized outlier detection algorithm.Our customized outlier detection algorithm is based on the observation that phase variation is apparent in the first few images (∼20 images, from 0 to 120 ms after diffusion preparation) immediately following diffusion preparation in the corrupted segments.In contrast, phase errors are negligible in the rest images (∼76 images, from 120 to 576 ms after diffusion preparation) within those segments.

F I G U R E 2
Measurement error detection using our customized outlier detection algorithm.(A) For each diffusion encoded segment of the reconstructed images, two unit-magnitude exponential phase maps (m s 1 and m s e ) are extracted from the first and the last images, respectively.These maps are then utilized to calculate the energy of the phase map difference (E s ) in Step 1.The calculation formula is shown and the unit-magnitude exponential phase maps from a corrupted segment are showcased on the left side.Subsequently, as illustrated on the right side, the diffusion encoded segments are iteratively clustered into two groups: one for corrupted segments (Ω corr ) and the other for non-corrupted segments (Ω noncorr ), based on their E s using a threshold (E thresh ) in Step 2. The value of E thresh is determined and updated by multiplying the average value in Ω noncorr by a coefficient  (e.g., 2.5) in each iteration.The segments with phase difference energies higher than E thresh are clustered into the corrupted group Ω corr , while the segments with phase difference energies lower than E thresh are clustered into the non-corrupted group Ω noncorr .(B) After convergence (e.g., 10 iterations) of the algorithm, the phase difference energies of all segments are plotted with the outliers (corrupted segments) marked by red pentagrams.The threshold E thresh is indicated by a red dashed line on the plot.The x-axis represents the number of segments, while the y-axis displays the values of the normalized phase difference energies.
The customized outlier detection algorithm involves two steps.As shown in Figure 2A, in Step 1, two unit-magnitude exponential phase maps (m s 1 and m s e ) are extracted from the first and last images of each diffusion-encoded segment and used to calculate the energy of the phase map difference (E s ).Specifically, In Step 2, the diffusion-encoded segments are iteratively clustered into two groups: corrupted segments (Ω corr ) and non-corrupted segments (Ω noncorr ), based on their E s using a threshold (E thresh ), which is updated by multiplying the average value of E s in Ω noncorr by a coefficient  (2.5 was used in this study) in each iteration.The segments with phase difference energies higher than E thresh are clustered into the corrupted group Ω corr , while the segments with phase difference energies lower than E thresh are clustered into the non-corrupted group Ω noncorr .As showcased in Figure 2B, after convergence (10 iterations were used), the phase difference energies of all segments are plotted with the outliers (corrupted segments) marked by red pentagrams.The threshold E thresh is indicated by a red dashed line on the plot.The x-axis represents the number of segments, while the y-axis displays the values of the normalized phase difference energies.
Due to both magnitude attenuation and phase variation in the corrupted images, conventional phase correction alone is not applicable in mdMRF.A straightforward correction strategy -removal, was adopted, and the dictionary was truncated to exclude corrupted segments.Pattern matching 1 was then performed between the truncated dictionary and the corrected images to simultaneously generate T1, T2, and ADC maps, as well as the proton density (M0) map.

Simulation experiments
Given the complexities of realistic scenarios where potential sources such as physiological motion, bulk motion, eddy currents, and other system imperfections, we used images reconstructed from mdMRF data to simulate a set of 2D mdMRF images, consisting of 28 segments with 96 images per segment (2688 total images), with the following parameters: FOV 300 × 300 mm 2 ; matrix size 192 × 192; resolution 1.5 × 1.5 mm 2 ; slice thickness 5 mm.Measurement errors (i.e., magnitude attenuation and phase variation) occurred in six random segments.The forward process was simulated to acquire highly under-sampled (48-fold) and multi-channel k-space using nonuniform fast Fourier transform (NUFFT) and coil sensitivity encoding with a set of coil sensitivity maps of 32 channels.White Gaussian noise with a SD of 0.01 relative to the signal intensity in the k-space center (averaged in the time dimension), was added to the k-space data for both real and imaginary components.
In the comparison experiment for Procedure 1, temporally local matrix completion was compared to its alternative, temporally global matrix completion.For the latter, cases with various ranks (small, 5; intermediate, 15; large, 35) were used.A low resolution of 3.75 × 3.75 mm 2 (2.5-fold resolution reduction relative to the original resolution) was used and  l = 5 × 10 −3 for both methods.To evaluate their performance in subsequent steps, temporally global subspace reconstruction was fixed for Procedure 2 to reconstruct high-resolution images, and ADC maps were also quantified by pattern matching after measurement error detection and removal.
In the comparison experiment for Procedure 2, temporally global subspace reconstruction was compared to its alternative, temporally local subspace reconstruction, which estimates a temporal subspace from the low-resolution images (by performing SVD and truncation within 1% approximation error on the low-resolution images) and performs subspace reconstruction for each segment separately, while temporally local matrix completion was fixed for Procedure 1 to reconstruct the low-resolution images.To maintain noise distribution for clear visualization in the interest of comparison, the additional sparsity constraint  u was set to 0. The normalized RMS error (NRMSEs) relative to reference was calculated for the reconstructed images of all competing methods.

In vivo brain imaging
In vivo brain imaging was performed as part of an institutional review board (IRB) approved study.Five healthy volunteers were recruited following written informed consent.One volunteer (Subject 1) was scanned for a feasibility experiment and a scan efficiency experiment.The other four volunteers (Subjects 2-5) were scanned for a robustness experiment with different scanners and varied imaging settings (varied b-values for diffusion encoding).
In the feasibility experiment, Subject 1 was scanned on a 3T MAGNETOM Prisma scanner (Siemens Healthineers, Erlangen, Germany) equipped with a 32-channel head coil.Two mdMRF scans were performed (with and without peripheral pulsation gating, respectively), with the following parameters: TI (21 ms) for T1 preparation; TE (30, 50, 65 ms) for T2 preparation; b-values (300, 700, 1000 s/mm 2 ) with three encoding directions for diffusion preparation; SSFP with constant flip angles of 10 degrees for data acquisition; single-shot spiral trajectory with an under-sampling factor of 48 was used for each time point (associated with each TR, TR = 6 ms); FOV 300 × 300 mm 2 ; resolution 1.5 × 1.5 mm 2 ; slice thickness 5 mm.Twenty-eight segments, each with 96 images, were acquired.The scan time was 26 s for both scans.A conventional MRF scan 31 and a single-shot EPI-based diffusion MRI scan were performed to obtain reference relaxation (T1 and T2) and ADC maps, respectively.
We selected six bilateral brain regions of interest (ROIs) on each subject for quantitative analysis (see Figure S1 for details).In addition, Subject 1 was scanned for a scan efficiency experiment (see Supporting Information Section B for details).

Simulation experiments
In Figure 3, the first and second rows show the reconstructed low-resolution images using different methods

F I G U R E 3
Simulation experiment for Procedure 1 (self-calibration) of the proposed reconstruction method.Temporally local matrix completion (fourth column) is compared to its alternative, temporally global matrix completion in cases with various ranks (first-third columns).The low-resolution images reconstructed using these competing methods are shown (first and second rows).The high-resolution images reconstructed using temporally global subspace reconstruction are also shown (third and fourth rows).Both magnitude and phase (unit: radians) images are shown, as well as normalized RMS error (NRMSEs) of the competing methods relative to the reference (fifth column).The ADC maps (unit: 10 −6 mm 2 /s) quantified from these high-resolution images (after measurement error detection and removal) are also shown (fifth row), as well as the error maps (sixth row).The red solid boxes highlight the measurement errors (magnitude attenuation and phase variation), and the red arrows point out the shading artifacts, which appear as uneven and overestimated ADC values.
in Procedure 1, and the third and fourth rows show the high-resolution images reconstructed by temporally global subspace reconstruction, where the temporal subspace is estimated from the low-resolution images corresponding to the same column.The fifth row shows the quantified ADC maps from the high-resolution images corresponding to the same column, and the sixth row shows the error maps.
When a small rank (=5) is used, temporally global matrix completion fails to recover signal components associated with small singular values, such as measurement errors indicated by solid red boxes in the fifth column of the first and second rows.These signal components of measurement errors are not represented by the temporal subspace estimated from the low-resolution images and, consequently, are not recovered in the high-resolution images (third and fourth rows, first column).Note that the failed recovery of measurement errors does not imply an accurate reconstruction of the high-resolution images.The errors persist in the under-sampled k-space data and finally affect the high-resolution image reconstruction, leading to shading artifacts (highlighted by the red arrow, sixth row, first column).
When a large rank (=35) is used, the low-resolution images partly recover the measurement errors, but suffer from heavy aliasing artifacts (first and secnd rows, third column).These aliasing artifacts significantly impact the temporal subspace, leading to aliasing artifacts in the reconstructed high-resolution images (third and fourth rows, third column).The presence of these aliasing artifacts reduces the sensitivity and robustness of measurement error detection.Consequently, the ADC map continues to exhibit shading artifacts, as indicated by the red arrow (sixth row, third column).Moreover, these aliasing artifacts tend to make the shading artifacts severer in the ADC map.
When an intermediate rank (=15) was used, the reconstructed images demonstrated a balance between recovering measurement errors and reducing aliasing artifacts (second column).Note that fine-tuning an optimal case so that measurement errors are accurately recovered while aliasing artifacts are eliminated remains challenging.The ADC map still exhibits shading artifacts (highlighted by the red arrow, sixth row, second column).
In contrast, temporally local matrix completion successfully recovers signal components corresponding to small singular values in the low-resolution images (first and second rows, fourth column), as well as in the high-resolution images (third and fourth rows, fourth column), along with minimum NRMSEs.The measurement errors were automatically detected and removed by a customized outlier detection algorithm, leading to an artifact-free ADC map (fifth row, fourth column).
Figure 4 shows high-resolution image reconstructions (first row) using temporally global subspace reconstruction and its alternative (temporally local subspace reconstruction), along with quantified ADC maps and error maps relative to ADC reference (second row).The temporally global subspace reconstruction shows superior image quality with higher SNR (along with lower NRMSE) Simulation experiment for Procedure 2 (subspace reconstruction) of the proposed reconstruction method.Temporally global subspace reconstruction (first column) is compared to its alternative, temporally local subspace reconstruction (second column).Both magnitude and phase (unit: radians) of the reconstructed images are shown (first row), as well as normalized RMS error (NRMSEs) relative to the reference (third column).The apparent diffusion coefficient maps quantified from these reconstructed images (after measurement error detection and removal), as well as the error maps, are shown (second row).The red arrows point out the shading artifacts.
compared to temporally local subspace reconstruction.The high-SNR images provide robust error detection and thus artifact-free ADC maps.In contrast, the low-SNR images of temporally local subspace reconstruction reduce the sensitivity and robustness of measurement error detection, leading to shading artifacts in the quantified ADC map (highlighted by the red arrow, second row, fourth column).The measurement errors occur at random locations spatially and temporally, affecting only a few segments.As depicted in Figure 2B, the random occurrence of phase variation leads to high phase difference energies randomly distributed, clearly identified as outliers.

In vivo brain imaging
Figure 6 shows the T1, T2, and ADC maps generated from an mdMRF scan without peripheral pulsation gating, using the proposed reconstruction method.Before the correction, there were severe shading artifacts (overestimated ADC values) in the ADC map (third column, first row).After correction, shading artifacts are alleviated, resulting Reconstructed images by the proposed self-calibrated subspace reconstruction.(A) Images corresponding to four representative time points (the first and the last images of the non-corrupted and corrupted segments, respectively), as pointed out by the small green and red solid boxes in the signal evolution curve.Both magnitude (first row) and phase (second row, unit: degree) images are shown.(B) The signal evolution curve along time points (TR) corresponds to a voxel as pointed out by the blue arrow in the reconstructed images.The solid and hollow arrows point out the magnitude attenuation and phase variation in the corrupted segments, respectively.

F I G U R E 6
T1, T2, and ADC maps were acquired by mdMRF scan without peripheral pulsation gating, using the proposed self-calibrated subspace reconstruction.Before the correction, there are severe shading artifacts (overestimated ADC values) in the ADC map (first row).After correction, the artifacts are alleviated (second row).The reference is shown as well (third row).
in improved ADC accuracy (third column, second row) which is closer to the reference (third column, third row).For more details regarding the cause of shading artifacts, please see Figure S2.T1 and T2 maps are unaffected by the measurement errors and consistent before and after correction; however, T2 values in the mdMRF-generated T2 map are slightly higher compared to the reference.
Figure 7 shows the T1, T2, and ADC maps generated by an mdMRF scan with peripheral pulsation gating, using the proposed reconstruction method.Although prospective gating was used, there are still severe shading artifacts (overestimated ADC values) in the ADC map before correction (third column, first row) due to failed peripheral pulsation gating.After correction, the shading artifacts are alleviated (third column, second row) and the ADC values are closer to the reference (third column, third row).
Figure 8 shows the quantitative analysis for the six brain ROIs.In the non-corrected cases both with and without peripheral pulsation gating, the average ADC values (third row) are significantly biased (overestimated, blue bar), while in corrected cases the ADC values are more accurate (red bar) compared to the reference (yellow bar).The average T1 values (first row) are highly consistent before and after correction for mdMRF scans both with and without peripheral pulsation gating and are consistent with the reference.Similarly, the average T2 values (second row) are highly consistent before and after correction for mdMRF scans both with and without peripheral T1, T2, and ADC maps acquired by mdMRF scan with peripheral pulsation gating, using the proposed self-calibrated subspace reconstruction.Before the correction, there are severe shading artifacts (overestimated ADC values) in the ADC map (first row).After correction, the artifacts are alleviated (second row).The reference is shown as well (third row).pulsation gating, but slightly (∼10 ms) higher than the reference.
Figure 9 shows the ADC maps generated by mdMRF scans without peripheral pulsation gating using the proposed reconstruction method for multiple subjects, scanners, and with varied imaging settings (diffusion encoding b-values).In each ADC map, the red point indicates the shading artifact (uneven and overestimated ADC values) location with the corresponding ADC value shown beside.In each case, severe shading artifacts exist in the ADC map before correction (first row) with overestimated ADC values.Note that these shading artifacts typically exhibit spatial variation across different cases, indicating that the measurement errors are inconsistent and difficult to predict.After correction (second row), shading artifacts are alleviated in terms of visual inspection and measured ADC values.
Table 1 shows the average and standard deviation of ADC values on the six ROIs across five subjects (Subjects 1-5).Notably, the average ADC values quantified by non-corrected mdMRF are over-estimated when compared to the reference.The over-estimation percentages are as follows: WM-Frontal left 54.24%, WM-Frontal right 45.40%, GM-Putamen left 83.94%,GM-Putamen right 24.63%, WM-Parietal left 21.45%, WM-Parietal right 24.01%.The large difference in the over-estimation percentage of GM-Putamen between left and right implies asymmetrical shading artifacts on ADC maps quantified Graphics showing the average T1, T2, and ADC values for six regions of interest (ROIs) as highlighted in Figure S1.Each graphic shows the comparison before (blue bar) and after (orange bar) correction, for the cases without (non-gated) and with (gated) peripheral pulsation gating, as well as the reference (yellow bar).

F I G U R E 9
ADC maps generated by mdMRF without peripheral pulsation gating, using the proposed self-calibrated subspace reconstruction, among different subjects and under variant imaging settings (different scanners and b-value combinations for diffusion encoding).In each case, severe shading artifacts occur in the ADC map before correction (top row) and can be alleviated after correction (bottom row).The red solid point overlayed on the ADC map points out the location of shading artifacts that appear as overestimated ADC values (top row), or the case after correction (bottom row), with the ADC value shown beside.

T A B L E 1
Average and SD of ADC values on six ROIs across five subjects.S3, showing that mdMRF scans along with the proposed reconstruction method achieve a minimum scan time of less than 20 s per slice, generating artifact-free T1, T2, and ADC maps.

DISCUSSION
This work proposes a new reconstruction method for mdMRF that can retrospectively alleviate shading artifacts caused by physiological motion-induced measurement errors.The primary physiological motion involves brain pulsation through cardiac pulsation. 32The proposed method enables mdMRF scans without the need for prospective gating methods (e.g., peripheral pulsation gating).Aliasing-free, high-resolution, and time-resolved images, where the measurement errors are accurately represented, can be reconstructed.The high-quality reconstructions enable robust detection and removal of diffusion-weighted segments, leading to artifact-free T1, T2, and ADC maps simultaneously.The feasibility and robustness of the method have been demonstrated across various scanners and imaging parameters for five subjects, achieving a high scan efficiency of less than 20 s per slice.The original mdMRF implementation using peripheral pulsation gating has a tradeoff between scan efficiency and method robustness due to the delay time between pulsation trigger and diffusion preparation in each diffusion encoded segment.If the delay time is too long, it can lead to skipped gating for the next segment, increasing scan time.On the other hand, if it is too short, brain pulsations can impact diffusion preparations, causing measurement errors and shading artifacts.Figure 7 shows an example of failed peripheral pulsation gating resulting in shading artifacts despite gating being used.A balanced delay time is subject-specific and challenging to fine-tune.Additionally, peripheral pulsation gating cannot deal with other types of physiological motion (e.g., CSF pulsation) and bulk motion.The proposed method eliminates the need for peripheral pulsation gating and effectively address these issues.It also enables artifact-free ADC quantification in failed gating cases, as demonstrated in Figure 7.
The proposed method removes the measurement errors by making them visible or detectable and subsequently removing them in the mapping step.However, conventional reconstruction methods used in MRF are not applicable.NUFFT reconstruction 1,33 fails to detect these measurement errors due to severe aliasing.Dictionary-based low-rank subspace reconstruction 27,34,35 does not account for measurement errors in the dictionary, making the measurement errors invisible in reconstructed images.In contrast, the proposed reconstruction method can reconstruct images where the measurement errors are accurately represented in some corrupted segments.This is because it is a data-driven LR subspace reconstruction that can adaptively estimate the temporal subspace to represent both the signals and measurement errors.Notably, the proposed reconstruction method differs from other data-driven LR subspace reconstruction methods proposed for MRF 36 as it does not require a fully sampled k-space center, offering greater flexibility and scan efficiency.The reason we employ temporally local matrix completion instead of temporally global matrix completion is that the randomly occurred motions between segments are uncorrelated, resulting in uncorrelation of the measurement errors between segments.This significantly weakens the global low-rank assumption.Similar observations have been reported in other studies regarding only the inter-shot phase variation. 8,18,37On the other hand, the measurement errors are determined by the accumulated phase during diffusion encoding, 7,25 and are constrained within each diffusion-encoded segment, without affecting the neighboring segments.Thus, the signals in each segment can be approximated by a local low-rank subspace.Employing low-resolution images in Procedure 1 offers several benefits: (1) the small size of low-resolution data reduces the computation burden of matrix completion, resulting in approximately three-fold reduction in computation time compared to using original (high-resolution) data; (2) the non-Cartesian (spiral) sampling used by mdMRF naturally enables denser sampling in the central k-space, improving the conditions for matrix completion, thereby enhancing image quality and reducing the number of iterations.Note that the low-resolution images reconstructed in Procedure 1, as well as the high-resolution images reconstructed using the alternative method (temporally local subspace reconstruction) in Procedure 2, can also be used to identify measurement error-corrupted images.However, we choose to use the high-resolution images reconstructed in Procedure 2 by temporally global subspace reconstruction.This choice is based on their superior image quality, which ensures automatic and robust detection of measurement errors.
One major limitation of the proposed self-calibrated subspace reconstruction method is its still high computation burden, although the self-calibration strategy has been used.The offline reconstruction process takes ∼40 min, with Procedure 1 requiring around 30 min and Procedure 2 taking around 10 min.However, the use of GPU-based NUFFT implementation and spiral trajectories with denser central k-space sampling can help accelerate the reconstruction process.Another minor limitation is the compromised temporal resolution, which increases from 6 to 30 ms due to the use of a sliding window (size = 5) to ensure sufficient signal recovery and achieve aliasing-free image reconstruction in the self-calibration step using temporally local matrix completion.Although a sliding window may reduce quantification sensitivity, 38 a window size of five is sufficient for accurate T1, T2, and ADC mapping in mdMRF.T2 values generated by mdMRF are a little higher than the reference.This could be due to different signal modeling, as the reference scan, that is, conventional MRF, 31 only considers T1 and T2 in the dictionary simulation.The ADC values quantified by mdMRF in the corrected case are lower than the reference.This difference may be due to the different types of pulse sequences used (SSFP used in mdMRF and EPI used in reference).Additional physical and image priors such as signal sparsity in spatial and/or temporal dimensions, 39 LLR, 19,21,24 structured low rank, 40,41 patch-based low rank, 42 and virtual conjugate coils, 43,44 can be incorporated to further improve image reconstruction, especially to improve the condition in Procedure 1.Another limitation is the removal strategy used for measurement error correction, which sacrifices a part of the collected data.However, the randomized wait times help mitigate the motion's impact by introducing incoherence between motions and diffusion encodings, resulting in only a few segments being corrupted by measurement errors.Future optimization or randomization of preparation modules could further enhance this incoherence, although it has not been explored in this work.

CONCLUSIONS
This study presents a new reconstruction method for mdMRF that can retrospectively alleviate shading artifacts caused by measurement errors due to physiological motions etc.The proposed reconstruction method reconstructs aliasing-free, time-resolved, and high-resolution images, where the measurement errors are accurately represented in some corrupted segments.By adopting the proposed reconstruction method, mdMRF scans now eliminate the need for peripheral pulsation gating.This advancement brings several advantages, including the preservation of sequence design flexibility, enhanced method robustness, and improved scan and reconstruction efficiency.

Figure
Figure 5A shows the first and the last reconstructed images of two representative diffusion-encoded segments, one non-corrupted and the other corrupted by measurement errors.The corrupted segments exhibit measurement errors (both magnitude attenuation and phase variation) in the first ∼20 images (0-120 ms after diffusion preparation).The first and last images of the 20th segment are shown for example (top right, highlighted by red dashed