Single‐shot spiral diffusion‐weighted imaging at 7T using expanded encoding with compressed sensing

The expanded encoding model incorporates spatially‐ and time‐varying field perturbations for correction during reconstruction. To date, these reconstructions have used the conjugate gradient method with early stopping used as implicit regularization. However, this approach is likely suboptimal for low‐SNR cases like diffusion or high‐resolution MRI. Here, we investigate the extent that ℓ1$$ {\ell}_1 $$ ‐wavelet regularization, or equivalently compressed sensing (CS), combined with expanded encoding improves trade‐offs between spatial resolution, readout time and SNR for single‐shot spiral DWI at 7T. The reconstructions were performed using our open‐source graphics processing unit‐enabled reconstruction toolbox, “MatMRI,” that allows inclusion of the different components of the expanded encoding model, with or without CS.


INTRODUCTION
In an ideal setting, the Fourier transform effectively maps image space into k-space. The Fourier transform provides optimal conditioning of the reconstruction problem, and fast image reconstruction by its inverse Fourier transform. 1 However, in practice, artifacts occur when unwanted field perturbations that arise from eddy currents 2 and field inhomogeneity/drift are present. Long readout times and aggressive hardware usage, such as the use of diffusion gradients 3 exacerbate these issues due to increased phase errors and stronger eddy currents, respectively, which leads to artifacts like ghosting, blurring and geometric distortion. Notably, non-Cartesian k-space trajectories [4][5][6][7] and ultra-high-field imaging 5,6,8 are more susceptible to artifacts from non-ideal fields. Although there are several retrospective image-based methods to compensate for specific kinds of artifacts (e.g., diffusion gradient eddy-current correction 9,10 ), these solutions typically apply to Cartesian trajectories. The expanded encoding model [4][5][6]8,[11][12][13] incorporates static and dynamic field evolution to correct for unwanted field perturbations during the reconstruction process. A B 0 map allows for the compensation of static field offsets, whereas dynamic fields can originate from gradient coils, eddy currents, and physiological events. By monitoring a sequence in real-time with NMR field probes, the field dynamics can be characterized using spherical harmonics up to second-or third-order in space. Field monitoring has been shown to improve image quality for Cartesian EPI in diffusion 14 and functional acquisitions, 15 as well as with spiral acquisitions. [4][5][6][7][8] To date, the expanded encoding model has been incorporated in reconstructions via least-squares (LS) optimizations with no explicit regularization. Instead, reconstructions from previous works have utilized early stopping of conjugate gradient iterations as a form of implicit regularization. However, choosing the correct number of iterations requires manual investigation of image quality at different iterations, and in lower SNR sequences such as diffusion MRI, noise amplification from improper selection of the number of iterations can have severe effects on the quality of computed parameter maps. Furthermore, LS optimization is unlikely to accurately solve low SNR problems if not accompanied by a regularization term. Therefore, in this work, we investigate complementing the expanded encoding model with 1 regularization (i.e., compressed sensing [CS]) [16][17][18][19][20][21] to make reconstructions more robust to noise and to eliminate the need for manual tuning of the number of iterations: The data consistency term is based on the image x, the acquired k-space data from the scanner y, and the expanded encoding model A exp . The "expanded" version of the encoding model includes additional phase terms from static and dynamic fields that are provided by a field map and field monitoring, respectively. For more details we recommend reviewing Eq. (3) from Kasper et al. 5 The "normal" encoding model includes only the phase from coil sensitivity maps and the nominal k-space trajectory. The regularization term ||Wx|| 1 promotes sparsity in the reconstructed imagex, using wavelets W as the sparsifying transform, and is the regularization weighting. Although the determination of is not trivial, several approaches have been proposed for its automatic selection. [22][23][24] CS reconstructions with expanded encoding were implemented using the MatMRI package, which is the only public toolbox we know of for expanded encoding model reconstructions, with or without CS. Through in vivo data and simulations, we show that these regularized reconstructions improve image quality compared to early stopping for single-shot spiral diffusion MRI at 7T.

Reconstruction
All reconstructions were performed using our open-source MatMRI toolbox, where "Mat" signifies extensive usage of matrix-vector operations in the reconstruction process (https://doi.org/10.5281/zenodo.4495476). MatMRI is a graphics processing unit (GPU)-enabled reconstruction framework, programmed in Matlab (MathWorks), that can handle any k-space trajectory with the option to include a B 0 map, higher-order coefficients fitted from the monitored field dynamics, and 1 -norm regularization. Matlab was chosen as the coding environment due to its prevalence in the MRI community and its straightforward GPU implementation that requires no other libraries or compilation steps. MatMRI only uses standard Matlab libraries and toolboxes to eliminate the need for compilation of binary files, which simplifies installation and usage. We implemented 1 -norm regularization using the undecimated wavelet transform 25,26 W, and the determination of the regularization weighting is automatically determined from WA T y, similar to our earlier work. 22 Figure S1 shows a comparison between different methods to determine that were evaluated with different numbers of wavelet levels and spatial resolutions from in vivo data reconstructions. From these preliminary comparisons, we decided to perform experiments with one wavelet level and to define = ∕2, where is the SD of noise plus noise-like artifacts in the high-pass wavelet coefficients (see Figure S1 for more details).
Implicitly regularized LS was implemented using the conjugate gradient method 27 with early stopping of iterations and CS (i.e., wavelet regularized LS) was implemented using balanced fast iterative shrinkage-thresholding algorithm (bFISTA). 25,28 Figure S2 shows how the number of conjugate gradient iterations impacts the trade-off between artifacts and noise amplification. 29 From these preliminary comparisons, we chose to run the conjugate gradient method with 20 iterations for all experiments. Although the number of iterations for the conjugate gradient has been shown to vary according to the conditioning of the image reconstruction problem, 29 Figure S2 shows that 20 iterations provides good reconstruction quality in all reconstruction settings without significant image degradation or noise introduction. For bFISTA, iterations were stopped when the data consistency residual (i.e., ‖ ‖ A exp x − y ‖ ‖ 2 2 ) changed by less than 10 −4 % in consecutive iterations, or if the iterations reached a maximum of 1000. Given this stopping criteria, bFISTA reconstructions tended to stop between 200 and 600 iterations. It is important to mention that a large number of iterations for conjugate gradient could lead to noise amplification whereas for bFISTA it only ensures convergence to the global optimum. Regardless of the selection of the objective function, k-space filtering was performed during reconstruction by setting k-space points outside the trajectory to zero. The delay between the field probe and MRI data was determined using a model-based retrospective algorithm. 30 All experiments were run as 2D reconstruction problems on a workstation with an Intel i9-7900× processor and 128 GB random-access memory with GPU Nvidia GeForce GTX 1080 Ti with 11 GB of GDDR5X memory.

In vivo
This study was approved by the Institutional Review Board at Western University and informed consent was obtained prior to scanning. Scanning was performed on two healthy volunteers on a 7T head-only MRI scanner (Siemens MAGNETOM Terra Plus, Erlangen, Germany), with 80-mT/m gradient strength and a 400-T/m/s maximum slew rate. Concurrent field monitoring was performed using a radiofrequency coil (32-channel receive and 8-channel transmit) with an integrated 16-channel 19 F commercial field-probe system (Skope Clip-on Camera) to obtain field dynamics that were fit to second-order spherical harmonics. 31 The collected k-space data were coil compressed to 20 virtual coils to improve reconstruction speed, [32][33][34][35] and noise correlation between receivers was corrected using pre-whitening before reconstructions. 36 A Cartesian dual-echo gradient-echo acquisition was used to estimate the B 0 map (FOV = 240 × 240 mm 2 , spatial resolution = 1.5-mm isotropic, TE 1 /TE 2 = 4.08/5.10 ms). The B 0 map was interpolated to match the in-plane resolution of in vivo data. From the first echo, we estimated sensitivity coil maps using ESPIRiT. 37 We acquired 15 DTI protocols with spiral trajectories to evaluate reconstruction performance of diffusion-weighted images and computed fractional anisotropy (FA) maps. The scans consisted of combinations of acquisitions with acceleration factors (R) from 2× to 6× and in-plane spatial resolutions of 1.5, 1.3, and 1.1 mm-all other sequence parameters remained constant: slice thickness = 3 mm, number of slices = 10, TE/TR = 33/2500 ms, FOV = 192 × 192 mm 2 . The peak gradient strength for the diffusion modules were 75 mT/m, and 46 mT/m for the spiral readouts, and for both cases the maximum slew rate was 333 T/m/s. The DTI protocol used monopolar pulsed gradient spin-echo encoding using a b-value of 1000 s/mm 2 and six directions, plus one b = 0 s/mm 2 acquisition. The scan time for each DTI protocol was approximately 18 s. Following image reconstruction, the MRtrix3 38 package was used to estimate the diffusion tensor and obtain FA maps.

Simulations
To simulate reconstructions for comparisons with a known ground truth, the first echo from the Cartesian gradient-echo acquisition was first used as x in the forward model, y = A exp x, to simulate k-space data. The trajectories used in the forward model were acquired during the 1.5 mm in-plane resolution diffusion MRI scans described above. A slice near isocenter was used due to the presence of several high-contrast regions and B 0 variations. After applying the forward model to the Cartesian reference to obtain the simulated spiral k-space, complex Gaussian white noise was added in spiral k-space to generate noisy data with approximate SNR values of [5,10,20]. Figure S3 shows the ground truth, B 0 map, virtual receive-coil sensitivities, trajectory, and estimated higher-order coefficients for both 2-and 4-fold accelerated single-shot spirals supplied to A exp .
Reconstructions were performed with and without components of the expanded encoding model (B 0 map and second-order coefficients) to determine their impact on reconstruction quality. The four settings were reconstructed with both LS and CS methods to analyze the impact of adding wavelet regularization. These

F I G U R E 1
Blurring effect as a function of spatial resolution and readout time. Spatial resolution contributes to blurring due to partial volume effects whereas readout time contributes to blurring due to T 2 * decay. In each zoomed region, the left side shows the mean diffusion-weighted image, averaged over the six diffusion encoding directions, from least-squares (LS) reconstructions, and the right side shows the average diffusion-weighted image from compressed sensing (CS) reconstructions. In all panels, the right hemisphere of the averaged image from LS was reflected onto the left side to directly compare anatomy between LS and CS averaged images.
simulations were also used to investigate how CS regularization and the expanded encoding model act as complementary tools to improve image quality. Finally, we coil compressed physical coils to [8,16,32] virtual coils to analyze the trade-off between coil compression, reconstruction speed, and reconstruction quality.
To quantitatively evaluate reconstruction quality between reconstructions and the ground truth, we used the normalized RMS error (NRMSE). Figure 1 shows mean diffusion-weighted images, averaged over the six diffusion encoding directions, for the

F I G U R E 2
Fractional anisotropy (FA) maps obtained from reconstructions with different acceleration factors (R) and in-plane spatial resolutions. Least-squares (LS) reconstructions were also performed with compressed sensing (CS) regularization. In all panels, LS reconstructions from the right hemisphere were reflected onto the left side to directly compare anatomy between LS and CS reconstructions. Top-left reconstruction was performed with the highest signal-to-noise ratio (SNR) condition whereas the bottom-right reconstruction was performed with the lowest SNR condition. different spatial resolutions and 2× and 4× acceleration rates. Observable blurring occurs due to either low spatial resolution (partial volume effect) or long readout times (T 2 * decay), or both. Higher R and spatial resolutions mitigates blurring, but this comes at the expense of lowering SNR. Despite being an average across images, it is still possible to appreciate a slight noise level from LS reconstructions, which is noticeably reduced at 4× acceleration using CS reconstruction. Figure S4a shows the mean diffusion-weighted images acquired at 6× acceleration, where the appearance of artifacts suggests that the undersampling was too aggressive. Figure S4b shows the corresponding reconstruction times for processing a single slice on a DTI protocol (7 acquisitions). Figure 2 shows FA maps obtained from reconstructions with different R and in-plane spatial resolutions from one volunteer. CS-based reconstructions improve SNR for all cases when compared to their LS counterparts and as SNR decreases either by acceleration rate or spatial resolution, these differences become more apparent. Despite improvements from CS, at 6× acceleration both methods show an inability to preserve diffusion metric integrity.

F I G U R E 3
Comparison of fractional anisotropy (FA) maps from spiral acquisitions from the two subjects acquired at combinations of in-plane spatial resolutions and R of 1.5-mm and 4×, and 1.1-mm and 5×, which amount to comparable readout times of 13 ms and 16.5 ms, respectively. Image SNR allows for a qualitative comparison, despite images not being co-registered. Figure 3 shows a comparison between FA maps obtained from spirals with comparable readout times for each volunteer. SNR loss when simultaneously accelerating and increasing the spatial resolution is ameliorated using CS, whereby CS reconstruction demonstrates better preservation of FA features than the LS reconstruction at lower SNR cases. This is similarly illustrated for color FA maps in Figure S5, which shows how CS reconstruction improves fiber tract delineation and directional uniformity.  Figure 4 shows the effect of adding second-order phase terms in the expanded encoding model and CS regularization as a function of the R. When the B 0 map was not included in the encoding model, the NRMSE was higher than 25%; hence, this case was not included in the figure. At an R of 2×, the impact of the expanded encoding model can be appreciated given the noticeably lower

F I G U R E 4
Comparison between least-squares (LS) and compressed sensing (CS) reconstructions. Reconstruction quality, measured in terms of the normalized RMS error (NRMSE), is presented as a function of the R for the case of an SNR of 20 and 16 virtual coils. The "Expanded" case uses the full expanded encoding model, while the "Partial" case does not include the second-order phase terms.
NRMSE. The errors from omitting the second-order terms are reduced when moving to higher accelerations due to shorter readout durations (i.e., faster traversal through k-space), reaching a minimum at a factor of 4×. For higher accelerations than 4×, the predominant source of error is noise. At all R, the lowest NRMSE occurs for CS with the full expanded model. Figure S6 shows the reconstructed images for 2×, 4×, and 6× acceleration without, with the partial, and with the complete expanded encoding model. Figure 5 shows the reconstruction performance for LS and CS reconstructions for the complete expanded encoding model for different numbers of virtual coils. Figure 5A shows normalized root-mean-squared error (NRMSE) as a function of R, and Figure 5B shows NRMSE as a function of SNR. Both panels illustrate that both LS and CS reconstruction quality are proportional to the number of virtual coils used at either higher R or noise levels; however, 16 or more virtual coils only show incremental improvement. When comparing CS to LS, CS reduces NRMSE for all cases with a larger improvement for higher R.

DISCUSSION
In this work, we presented a regularized reconstruction framework that combines the expanded encoding model and CS. The qualitative in vivo experiments showed similar trends compared to the quantitative simulations. Together, they demonstrate that the expanded encoding model and CS regularization are complementary tools for improving reconstruction quality. On one hand, at low R and high spatial resolution, the main source of reconstruction error arises from field perturbations given a long readout time. In this acquisition regime, the expanded encoding model provides the greatest benefit during the reconstruction process. Despite improvements, long readout times exacerbate T 2 * blurring that is not accounted for in the expanded encoding model, making low R unfavorable for single-shot spirals. Additionally, long readout times can be detrimental from a field monitoring perspective given that excessively long readout times may be on the same time scale as the lifetime of 19 F field probe signals, which would impair probe phase measurements and subsequent estimation of spherical harmonic coefficients. On the other hand, at high R and at high spatial resolutions, the main source of error arises from SNR reduction. In this acquisition regime, CS reconstruction complements the expanded encoding model by using its denoising properties and by removing high-frequency artifacts. In fact, the reduction of aliasing artifacts for 6x acceleration suggests that CS also improves recovery of missing data points for single-shot spirals. When using the automatic selection of regularization weighting factor (see Figure S1), CS did not introduce additional blurring in reconstructed in vivo images that can occur with wavelet-based over-regularization. Although our experiments were limited to uniformly accelerated spirals (e.g., artifacts for 6× acceleration), it is expected that variable-density spirals or pseudo-randomized trajectories 39 could improve performance at higher R by better transforming aliasing into noise-like artifacts. While undersampled Cartesian trajectories like EPI have a point-spread-function with large coherent peaks that are not ideal for compressed sensing, they would likely still benefit from explicit wavelet regularization that favors solutions satisfying known tendencies of medical images (i.e., sparseness in wavelet domain), instead of implicit early stopping regularization that has poorly defined regularizing behavior.
In our experiments, coil compression was used to reduce memory requirements and reconstruction time. Results suggest that for varying SNR levels, higher R should be accompanied by a higher number of virtual coils (≥16) to preserve reconstruction quality, regardless of using LS or CS reconstruction. We recommend using 20-21 virtual coils for a 32-channel receiver as it provides a good trade-off between reconstruction time and image quality, and as shown in previous work. 40,41 Both LS and CS reconstructions have limitations. For LS reconstruction, the number of iterations must be determined before running the reconstruction. Although the number of iterations according to Figure S2 was quite consistent for both the reconstruction of simulated and in vivo data, the non-explicit nature of the regularization is a drawback. While it is important for LS to not overestimate the number of iterations as to avoid the global minimum of the objective function, which can amplify noise, the choice of iterations must still be large enough to correct aliasing artifacts. Therefore, too few or too many iterations will degrade image quality, requiring fine tuning of this parameter. In contrast, the global minimum is sought for CS, and adding extra iterations will never impair image quality. For CS reconstruction, the selection of used here was inspired by our previous work 22 but performed assuming a Rayleigh distribution on the set of wavelet coefficients for one level of the wavelet transform. The additional heuristically determined factor of 0.5 was used to improve reconstruction quality and to avoid blurriness from over-regularization. Despite the heuristic for the determination of the regularization weighting through a Rayleigh distribution, after its selection the reconstructions showed better performance than LS reconstructions in all cases for different subjects, spiral trajectories, R, spatial resolutions, virtual coils, and noise levels. Finally, CS requires more iterations than LS with early stopping, which leads to a factor of 2×-13× increase in reconstruction time ( Figure S4b). Future work will consider optimization alternatives to improve reconstruction speed (e.g., unrolled CS with deep learning 42 ) and improvements to computational hardware to explore sub-millimeter DWI, where the large matrix sizes will lead to even more lengthy reconstructions. 7,43

CONCLUSIONS
In this work, we have introduced a compressed sensing extension to the MatMRI toolbox, which was used to investigate the impact of CS on single-shot spiral diffusion MRI at 7T. Simulations and in vivo acquisitions showed improved reconstruction quality and diffusion metric integrity when combining the expanded encoding model with CS, with CS providing the largest improvements at higher R.

SUPPORTING INFORMATION
Additional supporting information may be found in the online version of the article at the publisher's website. (D) 0th to 2nd order dynamic spherical harmonic coefficients and second order concomitant fields monitored for diffusion-weighted spiral trajectory readouts acquired at R of 2 and 4. Figure S4. (A) Average of diffusion-weighted images as a function of acceleration and in-plane spatial resolution. The left hemisphere is the reconstruction from least-squares (LS) whereas the right hemisphere is the reconstruction from compressed sensing (CS). The left hemisphere was reflected into the right hemisphere to present the reconstructions from CS. In this set of images, it is clear to observe the signal-to-noise ratio (SNR) reduction due to either acceleration or in-plane spatial resolution. The arrows indicate the location of an artifact observed for the LS reconstruction with an R of 6×. (B) Reconstruction times for LS and CS methods to process a single slice on a diffusion tensor imaging (DTI) protocol (7 acquisitions total). Figure S5. Comparison of color FA maps from spiral acquisitions from the two subjects acquired at combinations of in-plane spatial resolutions and R of 1.5-mm and 4×, and 1.1-mm and 5×, which amount to comparable readout times of 13 and 16.5 ms respectively. Figure S6. Effect of both the encoding model and R for LS (left side) and CS (right side) reconstructions using a signal-to-noise ratio (SNR) of 10. In all panels, CS reconstructions from the left hemisphere were reflected onto the right side to directly compare anatomy between LS and CS reconstructions. Panels (A-C) are reconstructed with the standard encoding model (i.e., the B 0 map and 2nd order spherical harmonic terms are not included). Panels (D-F) are reconstructed with the standard encoding model and B 0 map, and panels (G-I) are reconstructed with the complete expanded encoding model. Peach arrows depict a region where aliasing artifacts appear for high accelerations, and olive arrows depict a region sensitive to simplifications in the encoding model.