A correction algorithm for improved magnetic field monitoring with distal field probes

A significant source of artifacts in MRI are field fluctuations. Field monitoring is a new technology that allows measurement of field dynamics during a scan via “field probes,” which can be used to improve image reconstruction. Ideally, probes are located within the volume where gradients produce nominally linear field patterns. However, in some situations probes must be located far from isocenter where rapid field variation can arise, leading to erroneous field‐monitoring characterizations and images. This work aimed to develop an algorithm that improves the robustness of field dynamics in these situations.


INTRODUCTION
MRI has become the imaging modality of choice for the diagnosis of many pathologies, largely because of its noninvasiveness, excellent soft-tissue contrast, and ability to produce high-quality anatomical and functional images.However, reconstructed images can often suffer from varying degrees of artifacts and imperfections resulting from magnetic field perturbations that manifest during scan acquisitions.6][7] The ability to reduce artifacts depends on how well field perturbations can be corrected or accounted for, which is particularly important to improve diagnostic efficacy using diffusion imaging, which generates strong eddy currents.Field monitoring (FM) is a technique that uses field probes, [8][9][10][11][12] which individually measure NMR signals and accrue phase according to the local magnetic field.This has demonstrated great utility for determining dynamic field perturbations that can be incorporated in image reconstructions to correct artifacts.One particular commercially available FM system consists of 16 NMR field probes and can measure perturbations up to third order in space, 13,14 which is advantageous in correcting higher-order perturbations that commonly originate from complex field evolution such as aforementioned eddy currents 15 and patient motion. 16o accurately characterize patient-induced and time-dependent effects, monitoring must be performed concurrently with patient scanning, as opposed to sequential monitoring, which characterizes field dynamics in the empty scanner and uses the measurements to correct images from a subsequent identical acquisition. 170][21] The latter approach is advantageous, as it enables more streamlined setup workflows and designs that jointly consider both probe and RF geometry.
For concurrent field monitoring, field probes are distributed around the anatomy of interest and, ideally, are situated well inside the diameter spherical volume (DSV), where the main magnet and gradient fields are nominally homogenous and vary almost linearly, respectively.As such, field-probe phase accrual would follow near-linear behavior as a function of position in the DSV, and approximating these phase values using a spherical harmonics expansion can accurately capture the field patterns within the imaging volume.In instances in which field probes may reside near or outside the DSV where linearity and uniformity errors exceed acceptable limits, such as when using high-performance gradient coils with limited DSVs, when imaging over the large FOVs required for whole-body imaging, or when peripheral equipment prevents probes from being placed near the subject, accrued phase may deviate more from linearity, leading to solid harmonic coefficient fitting errors when determining field dynamics and a reduction in the quality of higher-order field correction during image reconstruction.
In this work, we propose a simple three-component correction algorithm that diminishes solid harmonic coefficient fitting errors introduced when probes are outside the DSV of the gradient coil and/or main magnet. 22In particular, this work features a head-only gradient system in which gradient nonlinearity is reported as being less than 7% at a DSV of 22 cm.We evaluate the algorithm by using our previously developed RF head coil with integrated field probes. 23In this RF coil, probes were located as close to isocenter as possible, pursuant to design conditions; however, the physical constraints of the RF coil required all probes to be placed in the nonlinear region of the head-only gradient coil.We therefore used this architecture as an exemplar to evaluate the performance of our correction algorithm.We demonstrate the benefits of incorporating higher-order coefficients by comparing the quality of concurrently monitored field dynamics and resulting brain images when using conventional fitting and our developed fitting methods, for first through third spatial order fits.The performance of the algorithm is demonstrated for both single-shot spiral and EPI acquisitions.We also compare computed field dynamics and reconstructed images informed by different levels of algorithm correction to assess the efficacy of each correction component.Using a phantom, algorithm performance is assessed quantitatively for the field probe coil probe arrangement, along with other combinations of probe positionings.Finally, we assess the effect of the correction algorithm on sequentially monitored field dynamics, where all probes were located within the DSV, and thus were not heavily influenced by higher-order phase variations.The correction algorithm has been made publicly available to promote dissemination. 24

Theory
As discussed by Barmet et al. 8 and Wilm et al., 13 the spatially and temporally varying magnetic field magnitude that encodes the MR signal during an acquisition can be expressed in terms of the static field and spatially smooth dynamic field components, which can be expanded using relatively few basis functions-most notably using the set of N L real-valued solid harmonics.The magnetic field evolution B(r,t) for a field probe located at position r and at time t can then be approximated as 8 where B ref (r) denotes the local magnetic field offset relative to the field that the spectrometer is tuned to for each probe; b l denotes the first out of N L solid harmonic, which can be categorized by spatial order, and c l is the corresponding field coefficient.The time integrated phase evolution is then described by where  is the respective gyromagnetic ratio, and k l now denotes the phase coefficients.For the first-order linear harmonics (e.g., b l [r] = x), k l become the traditional k-space coefficients.
To estimate (r,t) from field probe raw data, two preprocessing steps are first performed: (1) The phase of the raw data is extracted and unwrapped, and (2) γB ref (r)t is subtracted from the phase.The value of B ref for each probe is estimated using a calibration scan with no gradients applied.
For an array of N P probes, Eq. ( 2) can be cast into vector-matrix notation 8 as follows: with where P is referred to as the probing matrix, whose entries represent the solid harmonic basis functions up to the desired spatial order that are sampled by the probe positions.For a probe array consisting of 16 probes, P can be filled from zeroth up to the third spatial order (N L = 16).Probe positions can be determined using a separate calibration scan by acquiring probe FIDs when individual x-, y-, or z-gradients are on, assuming linear gradients are applied.The value of γ P denotes the gyromagnetic ratio of the probe sample, and γ is the gyromagnetic ratio of the nuclei being imaged.An approximate solution to Eq. ( 3) can be written as where P + = ( P T P ) −1 P T denotes the Moore-Penrose pseudoinverse of P. Therefore, to obtain the phase coefficients during a scan acquisition, the phase time courses of each probe are measured during the scan, and Eq. ( 4) is carried out independently at each point in time.
While not explicitly included in these equations, the calculation of solid harmonic coefficients also considers the phase contributions from second-order concomitant field terms, which require special consideration because they depend on the amplitude of the first-order linear harmonics. 25Accordingly, the expected phase from concomitant fields is computed from the first-order k terms, which is then subtracted from  P (t) before re-evaluating Eq. ( 4).This process is done iteratively to arrive at a final solution for k(t) and the time-varying concomitant fields. 14,17For notational simplicity, henceforth P and k(t) are considered to include the concomitant field basis functions and coefficients, respectively, and Eq. ( 4) is considered to carry out this iterative procedure.

Correction algorithm
The proposed correction algorithm considers modifications to how k(t) is calculated.The fully corrected fitting algorithm consists of three components: stepwise-fitting, weighted least squares (WLS) optimization, and weighted phase residual (WPR) removal, which are described in more detail subsequently.

Stepwise fitting correction
Conventionally, the coefficients k(t) are calculated simultaneously using Eq. ( 4) by incorporating all basis functions in P. When the probes are near isocenter, a simultaneous higher-order fit is able to accurately represent the time-varying magnetic field, as validated by improved image reconstruction. 15,23However, in the presence of substantial higher-order phase variation for distant field probes that may be outside the DSV (likely higher than the third-order harmonics that typical fits use), a simultaneous higher-order fit fails to accurately capture the zeroth-order and first-order behaviors of the gradient field, because the very high orders that the model is not equipped to fit are projected to lower orders.To understand this in one dimension (x) for a gradient exhibiting nonlinearity, if many probes sample the phase in the nonlinear region, a second-order fit up to x 2 inherits significant higher-order behavior and does an inadequate job at fitting the true linear behavior of the gradient field near isocenter (Figure 1A,B).To avoid this underfitting of zeroth-order and first-order terms, we propose solving k(t) one order at a time, which acts as an implicit regularization that strongly favors lower order terms.The zeroth-order and first-order terms are computed together in the first step, as is traditionally done for first-order fits, as calculation of the zeroth-order term requires knowledge of the probe positions (first-order solid harmonics), to not overfit the zeroth-order solution, except in the very ideal case in which probes are completely symmetric about isocenter.
After calculating k(t) for each order, the residual phase is calculated as per Eq. ( 5), and is used in the calculation of coefficients for subsequent orders: In the one-dimensional case, performing stepwise correction such that the basis functions (1, x) are used to calculate the zeroth-order and first-order coefficients followed by using x 2 in a subsequent step, instead of using all the basis functions at once for fitting, leads to a fit that much more closely aligns with the ground-truth field (Figure 1C).

WLS correction
To further suppress the effects from high-order phase variation, Eq. ( 4) was modified to include WLS optimization as follows: where W is a diagonal matrix of weights.Here, we choose that the weight for the jth probe depends inversely on the squared Euclidean probe distance from isocenter r j , because, by definition, the higher orders of phase variation that should be suppressed increase in magnitude with F I G U R E 1 (A) One-dimensional probe fitting scenario illustrating how fitting a phase profile with probes outside the linear region becomes more accurate when employing the proposed correction strategies.For the arbitrary direction x, the basis functions selected to fit the phase were (1, x, x 2 ).(B, C) "Conventional fit" indicates that all three basis functions were simultaneously used in the fit (B), whereas "stepwise" first fit (1, x) and then x 2 in a subsequent step (C).(D) Zoom-in shows how the step + weighted least squares (WLS) fit moves toward the inner probe.Notably, weighted phase residuals (WPR) only significantly adjusts the phase for the probes farthest from isocenter.(E) This is highlighted in the zoom-in, where an arrow indicates the large shift in the newly sampled phase by a distal probe.(F) The full correction scheme.
distance from isocenter: where the numerator normalizes the entries of W to be between 0 and 1. Building off the one-dimensional model, the incorporation of WLS makes the fitting more accurate by enabling probes affected less by nonlinearity to contribute more strongly to the fit (Figure 1D).

WPR removal
The third correction assumes that any residual phase after fitting using the first two corrections is primarily attributed to errors from unaccounted-for higher orders, especially for the probes far from isocenter.Weighted residual phases are thus removed from the initial probe phases, and the entire fitting procedure is repeated a specified number of times.The residuals that get subtracted are weighted using weights of min(W)/W.This causes probes far from isocenter to have a larger correction than probes near isocenter.Heuristically, we saw the best results when using two to three iterations for this phase-fitting process.As shown in Figure 1E, the act of removing residual phase from the probes can be thought of as displacing the sampled phase values such that the phase distribution is of a lower order in space, which in turn leads to a more accurate fitting of the lower-order coefficients.
A flowchart outlining how phase coefficients are calculated from probe phase data with the inclusion of the correction steps is depicted in Figure 2.

Algorithm access
The correction algorithm has been implemented in the MatMRI toolbox that is publicly available at https:// doi.org/10.5281/zenodo.4495476and https://gitlab.com/cfmm/matlab/matmri, 24,26 where the FM k-coefficients are computed offline using the raw field-probe data files collected with a commercial Clip-On Camera (Skope, Zurich, Switzerland) during an FM acquisition.

MR acquisition
Two healthy volunteers and a silicone oil phantom were scanned on a 7T head-only MRI scanner (Siemens MAG-NETOM Terra Plus) operating on the VE12U software platform, at Western University's Center for Functional and Metabolic Mapping (80 mT/m gradient strength and 400 T/m/s maximum slew rate).This study was approved by the institutional review board at Western University, and informed consent was obtained before scanning.Diffusion-weighted single-shot spiral acquisitions were performed with parallel imaging acceleration rates of 2, 3, and 4. The imaging parameters were as Flowchart outlining the calculation of phase coefficients from probe phase data while using the proposed correction steps: stepwise correction, weighted least squares, and weighted phase residuals.The values of k m (t) and P m contain all the coefficients and spherical harmonic bases, respectively, of the mth order.For m = 1, the zeroth-order and first-order harmonics (l = {0,1,2,3}) are included in k m (t) and P m ; for m = 2, the second-order harmonics (l = {4,5,6,7,8}) are included, and so on, for higher orders.The weighted vector W is calculated as per Eq. ( 7) and is reshaped into a diagonal matrix when multiplied by P m .Contribution of accrued phase from concomitant field terms occurs during the "calculate phase coefficients" step.A Cartesian EPI diffusion-weighted acquisition was also acquired using the same six-direction diffusion scheme and the following imaging parameters: acceleration rate = 2, partial Fourier = 6/8, FOV = 192 × 192 mm 2 , in-plane resolution = 1.5 × 1.5 mm 2 , slice thickness = 3 mm, number of slices = 10, TE/TR = 45/2500 ms, and flip angle = 70 • .
A Cartesian dual-echo gradient-echo acquisition was used to estimate B 0 maps for inclusion in a model-based reconstruction to correct for static off-resonance effects.The imaging parameters were as follows: FOV = 240 × 240 mm 2 , spatial resolution = 1.5-mm isotropic, and TE 1 /TE 2 = 4.08/5.10ms.From the first echo, sensitivity coil maps were estimated using ESPIRiT. 27

2.3
Field monitoring

Concurrent field monitoring
FM was performed concurrently using a Skope Clip-On Camera consisting of 16 transmit/receive 19 F field probes that are integrated into a 8-channel transmit, 32-channel receive RF head coil, 23 for both in vivo and phantom acquisitions.The mean distance between field probes and isocenter as reported in the CAD design is 13.5 cm (range: 11.9-14.9cm), with probes being 8%-35% outside the DSV.

FM probe array
To acquire more probe locations and enable higher than third-order fits for phantom investigations, the scanner bed holding the field probe coil was translated in 1-cm increments along the z-direction, and subsequent calibrations and field monitoring measurements were performed.This was repeated over 14 positions.

Sequential FM
Sequential FM of an identical diffusion-weighted single-shot spiral acquisition was also performed before integration of the field probes into the RF coil.In this arrangement, all probes were located within the DSV, having an average distance from isocenter of 9.0 cm (range: 8.2-9.9 cm).

Image reconstruction
Image reconstruction was performed in MATLAB using the in-house-developed MatMRI toolbox 24,26 that uses an iterative expanded encoding model-based reconstruction. 13All images were reconstructed using 20 iterations with the conjugate gradient method and coil compression [28][29][30][31] to 20 virtual coils was performed to improve reconstruction speed.Synchronization delay between the MRI and field-probe data was corrected using an automatic, retrospective algorithm. 32Noise correlation between receivers was corrected using prewhitening before any reconstructions. 33Before reconstruction, field dynamics were also adjusted to account for vendor-Maxwell corrections and eddy current compensation that are not measured by the field probes, using methods described previously by the authors 23 and reproduced here with notable differences expounded upon.Namely, additional B 0 and linear Maxwell terms are produced, given the asymmetric gradient coil design, 25 and so the MRI scanner corrects for the linear terms automatically by applying dynamic gradients to counteract them and corrects for the B 0 terms by applying the appropriate demodulation phase. 34The linear correction was presumed to be measured by the field probes, whereas the DC correction is invisible to the field probe system and was thus reversed during reconstruction using Maxwell term predictions based on the first-order field probe measurements.Similarly, the scanner performs B 0 eddy current compensation using similar demodulation, which is also invisible to the field probes and needs to be accounted for in measured field dynamics retrospectively.This is done by estimating the zeroth-order corrections via offline simulation of the full acquisition (including diffusion gradients) and reversing them during reconstruction in a similar manner to the DC Maxwell terms.Following image reconstruction, fractional anisotropy (FA) maps of the brain were computed using the MRtrix3 package. 355 Data analysis

k-Coefficient comparison
To assess changes in field dynamics for all three fitting orders with the addition of the correction algorithm, first-order, second-order, and third-order k-coefficients were calculated the conventional way and with the correction algorithm, and their time courses were plotted.Difference profiles were also calculated and plotted.Fourth-order-fit k-coefficients were also calculated by including 90 probe positions from the experiment described in Section 2.3.2 (positions ranging from 7 cm to 15 cm from isocenter).

Qualitative image comparison
In vivo images were reconstructed using concurrently acquired field dynamics, and non-DWI, DWI, and FA map quality was qualitatively assessed.Similar analysis was performed for both spiral and EPI trajectories.Difference images were calculated with respect to the images that used fully corrected second-order field dynamics, which appeared to provide the best image quality.Additionally, to evaluate each correction component's contribution to overall correction, image reconstructions using second-order fits of spiral trajectories were performed using no correction algorithm, only stepwise correction, stepwise+WLS, and all three correction techniques (full correction), followed by qualitative inspection of images and FA maps.Difference images were calculated with respect to the images informed by fully corrected field dynamics.Images and FA maps from third-order sequential monitoring were also compared to evaluate correction performance when probes were situated well within the DSV.Finally, to assess residual blurring, image quality for the different accelerated spirals was compared.

Quantitative comparison
Performance was assessed quantitatively by calculating the normalized RMS error (NRMSE) of reconstructed phantom DWI with a ground truth.Given the very low diffusivity of the silicone oil phantom, diffusion-weighted acquisitions exhibited negligible diffusion weighting. 36ccordingly, the non-DWI was used as a ground truth.The DWI reconstructions used concurrently monitored field dynamics, whereas the ground-truth image reconstructions used third-order fully corrected field dynamics, calculated from a virtual probe array consisting of 25 probes located primarily near the edge of the DSV, acquired as described in Section 2.3.2.The average probe distance from isocenter was 10.3 cm (range: 6.9-12.1 cm).The NRMSE of DWI was computed relative to the ground-truth GT , where S DWI and S GT are the DWI and ground-truth signal values, respectively, where the summation was across all voxels, and then averaged over all slices and diffusion directions.Similar to the qualitative analysis, average NRMSE was compared for images reconstructed using first-order through third-order k-coefficients calculated with the conventional method and with the full correction algorithm.Furthermore, average NRMSE was compared for images using subsets of the correction algorithm.To investigate correction effectiveness in various probe scenarios, several other virtual probe arrangements were investigated: 16 probes on the boundary of the DSV (mean: 9.9 cm, range: 6.9-12.0cm), probes that were moved outward by approximately 15% further than the original field probe coil positions (mean: 15.9 cm range: 13.5-17.9cm), and nonoptimal probe spacing where a probe from the boundary DSV case was substituted for a probe position farther from isocenter.Finally, to assess the effect of an overdetermined system on third-order fits, 25 probes outside the DSV were used to fit the field dynamics and compared with the traditional 16 probes via NRMSE.

RESULTS
Third-order k-coefficients calculated using the full correction method (stepwise + WLS + WPR), which were acquired from a spiral DWI concurrently monitored trajectory (Figure 3A), showed significant zeroth-order to third-order trajectory differences when compared with conventionally calculated k-coefficients (Figure 3B).
Applying the stepwise correction appeared to induce the largest difference in the profiles out of the three correction strategies.When using a second-order fit, the differences in these terms is smaller yet still apparent.Similar to third-order fit, stepwise appears to induce the largest change in coefficients.First-order fits showed small yet observable differences in trajectories when adding WLS and WPR correction.Note that stepwise correction is not applicable for first-order fits.Comparable results were observed for all basis functions (Figure S1A).Notably, results from the full set of basis functions consistently showed larger contributions with conventional fitting for terms having z dependency, especially those having z 2 and z 3 dependency, such as k 6 and k 12 , respectively (Figure S1B).This observation also applied to fourth-order terms (Figure S2).Also, the deviations varied significantly depending on DWI gradient direction, likely due to diffusion gradient-induced eddy currents (Figure S3).Images and FA maps reconstructed with and without k-coefficient correction, for first-order through third-order fits of concurrently monitored FM data, show degrading image quality when using higher-order conventional fits (Figure 4).When correcting the k-coefficients, image quality is improved, particularly for FA maps.This is evidenced by the reduction in blurring and artifactually elevated FA in the gray matter and CSF for higher-order fits.FA map zoom-ins show how k-coefficient correction enables the accurate use of second-order fits, which leads to improved image quality when compared with images (A) Fully corrected concurrently monitored k-coefficients generated with a third-order fit from a diffusion-weighted single-shot spiral acquisition with an acceleration rate of 2. (B) Difference between the fully corrected k-coefficients and conventionally calculated coefficients (black), coefficients corrected using only stepwise correction (red), and coefficients corrected using stepwise + weighted least squares (WLS) correction (blue), for maximum order of the fits ranging from first through third order.Sample basis function is displayed for each order.Insets added in the first-order plots highlight the observable differences in the profiles from 15 ms onward.using first-order fits.This is illustrated more effectively in Video S1, which shows apparent DWI blurring reduction and improved FA map integrity when transitioning from first-order corrected images to second-order corrected images.The zoom-ins and video also enable effective comparisons of images reconstructed using first-order fits, which show clear improvements in quality when correcting the k-coefficients.Furthermore, DWI and FA maps informed by corrected third-order fits were slightly blurrier and noisier, respectively, than those informed by corrected second-order fits.Similar trends were observed for images and FA maps (Figure 5 and Video S2), as well as for k-coefficient difference profiles (Video S2) from the diffusion-weighted EPI acquisition.
In Figure 6, images and FA maps reconstructed using second-order concurrently monitored k-coefficient data corrected using only stepwise correction appeared to provide most of the correction benefit to the images.Despite this, individually adding WLS and WPR corrections further improved image quality and integrity of FA maps, which is more clearly shown in the zoom-ins, as well as in Video S3.
An example showing average DWI and ground-truth reconstructed phantom images can be seen in Figure S4.

F I G U R E 4
Non-diffusion-weighted images (non-DWI) (b = 0 s/mm 2 , left), single-direction DWI (b = 1000 s/mm 2 , middle), and fractional anisotropy (FA) maps (right) reconstructed from single-shot spiral acquisitions, acceleration rate (R) = 2 using first-order through third-order concurrently monitored field dynamics.Reconstructions incorporated conventionally calculated k-coefficients and corrected k-coefficients.Difference images (left hemisphere) were calculated with respect to images that used fully corrected second-order field dynamics.Zoom-ins highlight improved FA map integrity when using corrected first-order fits, which are further improved using corrected second-order fits.

F I G U R E 5
Non-DWI (b = 0 s/mm 2 , left), single-direction DWI (b = 1000 s/mm 2 , middle), and fractional anisotropy (FA) maps (right) reconstructed from single-shot EPI acquisitions (R = 2) using first-order through third-order concurrently monitored field dynamics.Reconstructions incorporated conventionally calculated k-coefficients and corrected k-coefficients.Difference images (left hemisphere) were calculated with respect to images that used fully corrected second-order field dynamics.Zoom-ins show improved FA map integrity when using corrected first-order fits, which are further improved using corrected second-order fits.

F I G U R E 6
Non-DWI (b = 0 s/mm 2 , top), single-direction DWI (b = 1000 s/mm 2 , middle), and fractional anisotropy maps (bottom) reconstructed from single-shot spiral acquisitions (R = 2) using second-order concurrently monitored field dynamics.Comparison of image quality for reconstructions incorporating conventional, stepwise-only corrected, stepwise + weighted least squares (WLS)-corrected, and fully corrected (stepwise + WLS + weighted phase residuals k-coefficients.Difference images calculated with respect to fully corrected images are shown in the left hemisphere of each image.DWI blurring reduction and improved FA map integrity is seen as the various components of the correction are added. Figure 7A illustrates a reduction in NRMSE values for the reconstructions when using the full correction strategy to fit the k-coefficients as opposed to the conventional fit, with the second-order fit providing the lowest NRMSE.Third-order conventional fit presented a significantly higher NRMSE than first and second order.When investigating the effect of correction strategy on NRMSE, the first two correction strategies provided observable reductions in NRMSE compared with the use of conventional fitting, with the reduction due to WPR correction only being incremental (Figure 7B).Interestingly, however, the DWI images become more consistent between diffusion directions, as evidenced by the smaller SD after WPR is added to the algorithm.
When 16 probes were positioned on the boundary of the DSV, the NRMSE was significantly lower than the value for the 16-probe arrangement in the field probe coil, yet the correction algorithm still considerably reduced the NRMSE.When these probes were moved outward by approximately 15% further than the original field probe coil positions, the NRMSE increased further.Correction resulted in a decreased NRMSE, but the NRMSE did not reach the values listed when the probes were closer to isocenter (Figure 8A).When only one probe was located (A,B) Comparison of normalized RMS error (NRMSE) between DWI and non-DWI ground-truth phantom images for concurrently monitored single-shot spiral acquisitions, averaged over all diffusion directions and slices, when performing first-order through third-order fits using conventionally calculated and fully corrected k-coefficients (A) and when performing second-order fits incorporating conventional, stepwise-only corrected, stepwise + weighted least squares (WLS)-corrected, and fully corrected (stepwise + WLS + weighted phase residuals) k-coefficients (B).Error bars represent the SD of the mean across diffusion directions.

F I G U R E 8
(A-C) Average normalized RMS error (NRMSE) for phantom DWI reconstructed using conventionally calculated and fully corrected k-coefficients with third-order fits, for mean probe distances of 15.9 cm, 13.5 cm, and 9.9 cm (A), when all probes are uniformly situated on the diameter spherical volume (DSV) boundary (mean distance = 9.9 cm) compared with when a single probe is substituted for a probe located well outside the DSV (B), and probe arrays outside the DSV using 16 probes (blue) and 25 probes (red) (C).

F I G U R E 9
Comparisons of concurrent and sequential field dynamic correction.Non-DWI (b = 0 s/mm 2 , top), single-direction DWI (b = 1000 s/mm 2 , middle), and fractional anisotropy maps (bottom) reconstructed using (left to right) concurrently monitored third-order conventional k-coefficients, concurrently monitored third-order corrected k-coefficients, sequentially monitored third-order corrected k-coefficients, and sequentially monitored third-order conventional k-coefficients.outside the DSV, a large increase in NRMSE was observed, which was reduced with correction (Figure 8B).For probes located outside the DSV, when performing a third-order fit with 25 probes versus the traditional 16 probes, the NRMSE of reconstructed images was significantly lower when 25 probes were used to conventionally calculate the coefficients.Full correction still provided a slight reduction in NRMSE in this case (Figure 8C).
Correction of sequentially monitored third-order k-coefficients showed small differences in k-coefficients (Figure S5), with the degree of difference being far smaller in magnitude compared with the concurrent case.This is also illustrated via image reconstructions and FA maps in Figure 9, which shows negligible differences in image quality when using corrected sequentially monitored k-coefficients, as opposed to the substantial image-quality improvements observed for corrected concurrently monitored k-coefficients.
Accelerating the spiral acquisitions by 3-fold and 4-fold effectively minimized the blurring observed in 2-fold accelerated non-DWI and DWI (Figure 10), with 4-fold being the least blurry.

DISCUSSION
In this work, we present a solution for correcting computed k-coefficients from field monitoring data for situations in which higher-order phase variations significantly affect spherical harmonics fitting, such as the case that arises when field probes are placed outside the homogenous/linear region of the main magnet and/or gradient coil.
When implementing a third-order fit, the conventional inverse problem that solves all orders simultaneously introduces large errors for fitted terms, as evidenced Comparison of average DWI reconstructed from single-shot spiral acquisitions using acceleration factors of 2, 3, and 4 (left to right, respectively).Fully corrected second-order concurrent monitored field dynamics were used in the expanded encoding model.Zoom-ins are shown under each respective image.Reduction in blurring can be seen as the acceleration rate is increased.by the large deviations in k-coefficients, significant image artifacts, and large NRMSE value (Figures 3, 4, 7A).The large sensitivity to error is due to the problem being a balanced system of equations with an equal number of equations (i.e., probes) to unknowns (i.e., spherical harmonic k-coefficients), which causes the probes far from isocenter in regions of high nonlinearity to have a large erroneous effect on the solution.To improve the conditioning of the probing matrix, more probes than basis functions are required.For instance, when using a probe array with 25 probes and performing a third-order fit, conventional fitting was shown to perform nearly as well as when using 16 probes with the full correction strategy (Figure 8C).However, creating a large probe array from multiple scans is a cumbersome process that would be required every time a new sequence is performed, and ultimately is not applicable to concurrent field monitoring.Larger errors were exhibited by terms having stronger z-dependency, which is consistent with the rapidly changing fields along the z-direction that are expected for a head-only scanner and gradient coil (for both eddy currents and applied gradients).In other words, spherical harmonic orders greater than 3 are likely required for probes far from isocenter, but there are insufficient probes for this type of model.Figure S2 illustrates phase accrued at an arbitrary position approximately 14.2 cm from isocenter for terms up to the fourth order using data from 90 probes.When using the conventional fitting method, most second-order to fourth-order terms showed apparent phase accrual that is reduced with application of the full correction fit.Even following correction, nonnegligible fourth-order phase accrual remained.When the probes are close to isocenter, these problematic high-order terms that scale with the fourth or higher powers of radius become negligible compared with the lower orders.These k-coefficient errors become less apparent as the fit order is reduced, likely due to the system becoming overdetermined, but are still present and possible to mitigate with the proposed fitting algorithm (Figure 3).Although the stepwise correction is the most effective of the three strategies used, additional image-quality improvement is observed with the addition of each of the WLS and WPR corrections (Figures 6, 7B and Video S3).For optimal performance, we suggest using WPR correction up to fits that are one order lower than the allowable fit (i.e., for a 16-probe system that allows up to third-order fits, this correction is expected to perform well up to second order).This is because the WPR correction uses leftover residual phase in its correction, but fitting 16 probes to third order will leave behind little to no residual phase to use in the WPR correction in subsequent iterations.To be consistent with the correction techniques used, WPR correction was used in correcting third-order fits in this work, given that there was little difference in image quality when and when not using WPR correction for third-order fits.Altogether, the use of these corrections enables the production of images of higher quality when using higher-order fits compared with using a first order fit alone (Figures 4, 5, and Video S1).However, in the event that a first-order fit is chosen by the user for image reconstruction, WLS and WPR corrections still provide better quality than when not using the correction strategies (Figures 4, 5, and Video S1).
Interestingly, images were of slightly lower quality overall when using corrected third-order k-coefficients compared with corrected second-order k-coefficients.This is likely because all three components of the proposed algorithm capitalize on overdetermined solutions, and a third-order solution has balanced numbers of equations and unknowns for 16 probes.In contrast, the improvements in image quality for second-order corrected fitting versus both first-order corrected fitting and second-order conventional fitting (Figure 4) suggest that the second-order terms are more accurately estimated via the proposed methods.The larger distortions observed for DWI when using the conventional method to calculate coefficients compared with non-DWI is likely due to diffusion gradient-induced eddy currents, which likely have strong spatial variation along the z-direction, given the short spatial extent of the head-only gradient coils and possibly the cryostat of the head-only MRI that the eddy currents primarily reside in.
Blurring appears to be present in some images even after the correction, which may be a combination of uncorrected phase errors and T 2 * blurring.However, a modest acceleration rate of 2 was chosen for this work to better demonstrate the effects of the corrections, because errors in k-coefficients are exacerbated by slower traversal through k-space (similar to distortions/artifacts from uncorrected B 0 inhomogeneity).When using acceleration rates of 3 or 4, which have been shown to have low g-factors for single-shot spiral acquisitions, 37 this blurring is greatly reduced, likely given the reduced effects of T 2 * blurring and uncorrected phase errors at these shorter readouts (Figure 10).
When implementing WLS correction, different weighting schemes were tested, including weightings based solely on probe distance in the z-direction, probe signal magnitude, and probe signal decay rate.However, net distance from isocenter was qualitatively most effective in consistently improving image quality.The poor performance for distance along z, even though it is the higher-order modes with large z-terms that are most affected by the corrections, is due to the fact that probes that are distant in both z and x/y likely reside in regions with greater inhomogeneity than those distant only along z.Interestingly, full removal of the most distant probe(s) consistently reduced image quality, suggesting that even the farthest probes still have useful information to provide during the fits.Nevertheless, there may be more optimal choices for probe weightings.
The negligible difference in k-coefficients (Figure S5) and image quality (Figure 9) when correcting sequential FM data was expected, given that all probes were within the vendor-defined DSV of the scanner; accordingly, image reconstructions using sequentially monitored data were already of high quality.These results also suggest that the algorithm is suitable for routine use in sequential or concurrent monitoring environments that do not experience large phase errors, as no reconstruction errors were introduced by the algorithm.Also, if all the probes can be positioned within the DSV, fitting to the highest order possible would likely be appropriate, given the lack of errors expected.Accordingly, in this situation the WPR correction could be omitted, as it is expected to have little effect on the fit when the system is not overdetermined.
The proposed algorithm requires very few modifications to the original workflow, and negligible computation time is added.The duration of the conventional fit on our hardware for a single slice takes approximately 180 ms.The inclusion of the correction strategies' stepwise fitting, WLS, and WPR correction on average adds 5 ms, 3 ms, and 45 ms of computation time, respectively.The proposed approach is also effective in correcting different trajectories, as shown for a diffusion-weighted EPI trajectory of the same resolution and a comparable readout time (Figure 5 and Video S2).The k-coefficient errors present in the EPI acquisition translate to different artifacts than for the spiral acquisition, with bulk shifts in the phase-encoding direction, compression, and ghosting being observed.This translates to significantly corrupted FA maps for higher-order uncorrected fits, which are corrected using the proposed correction techniques.Notably, similar correction improvements were observed for the second volunteer (Figure S6).
A clear guideline for when the correction algorithm can provide improvements in image quality remains unclear.However, a rule of thumb from personal findings is that if there is a zeroth-order phase difference of over 0.1 radians when fitting the k-coefficients with a conventional first-order fit versus a third-order fit, then the higher order fits likely corrupt the lower order coefficients.This suggests a potential regime of complex field behavior that may benefit from the algorithm.
Using the data from multiple scans described in Section 2.3.2,other probe arrangements were explored to assess algorithm effectiveness in different situations.The small improvement in image quality when correcting probe data on the boundary of the DSV suggests that higher-order behavior is less influential but can still affect probes at this distance, and that the algorithm may be beneficial in such a case and/or when nonoptimal probe spacing is exhibited.Additionally, when probes were moved outward by approximately 15% further than the original field probe coil positions, as expected, the NRMSE was higher than the field probe coil arrangement when using conventional coefficient fitting due to even more problematic higher-order effects.The reduction in NRMSE implies that the algorithm may still provide a benefit when probes uniformly extend even as far as 45% beyond the DSV on average, but evidently the performance declines.Finally, when only one probe was located substantially outside the DSV, a surprisingly large increase in NRMSE was observed, suggesting that the probing matrix becomes even more poorly conditioned (Figure 8B).Fitting using the correction algorithm once again showed a large improvement in image quality.This may be useful when one or several probes have to be positioned further than the rest of the array due to design constraints.
It is important to note that while the algorithm may be advantageous in the last two scenarios described, moving the probes too far away may also cause rapid signal decay from dephasing, causing phase and resulting coefficients to be severely corrupted, leading to other image errors.While the performance of the algorithm has been investigated in our custom-made RF coil and with various positions and numbers of probes, it has only been tested in limited conditions, namely, a head-only gradient coil.Additionally, the probe arrays constructed from multiple scans with differing bed positions had unoptimized probe positions, which may have compromised calculated coefficients more than optimally designed arrangements.Finally, it was difficult to access uniform probe locations closer than the edge of the DSV because of the restricted movement of the field probe coil to along only the z-direction.
A fourth correction step was also investigated.Because calibrations rely on linear gradient assumptions to calculate probe positions, large errors in calculated probe positions can be expected when probes experience large phase deviations from linearity, such as when placed outside the DSV.Therefore, this correction method used vendor-provided gradient nonlinearity characterizations 38,39 to adjust the calibrated probe positions, and the x, y, and z k-coefficient basis functions supplied to P when fitting the phase coefficients.However, it was observed that while this correction independently provided an improvement in image quality, the benefit was negligible when combined with the other correction techniques and was not deemed necessary, given the additional computation time required to correct the probe positions (Figure S7), and due to the fact that the three-step algorithm can remain vendor-agnostic and be made entirely open-source, as it does not rely on proprietary vendor information.The poor performance of this approach was likely because it does not correct for eddy current modes with rapid spatial variations far from isocenter.Also, the vendor characterization may only be valid within the DSV of the gradient coils.Although it may be possible to calibrate the field nonlinearity of the gradients and eddy current modes outside the DSV, it would likely require a lengthy calibration process.
Another potential correction approach could be to introduce a tailored regularization term in Eq. ( 4) that penalizes large contributions of z-dependent spherical harmonic terms.However, this approach would require tuning of the regularization strength for the different basis functions and would require knowledge of the spherical harmonic basis functions that are more prone to error.In contrast, the proposed approach makes no assumptions about the coil/scanner geometry (e.g., perhaps coils with x or y asymmetry may also benefit from it 40 ) and does not require tuning of parameters.That said, regularization could also potentially be used in combination with the proposed algorithm, or the proposed algorithm could be used to identify problematic basis functions that could be penalized with regularization, which are promising avenues for future work.

CONCLUSIONS
A correction algorithm that successfully improves the reliability of phase-corrupted FM data is presented.This encourages the use of FM systems to achieve higher-quality images even when probes may lie outside the imaging volume, including complex FM environments such as concurrent FM in head-only scanners and imaging of anatomy distant from isocenter with large FOVs.
are publicly available in the MatMRI toolbox: https:// doi.org/10.5281/zenodo.4495476and https://gitlab.com/cfmm/matlab/matmri. 25A demo calculating phase coefficients from raw probe data acquired from a spiral acquisition with and without the proposed corrections is included.
single-shot spiral acquisitions (R = 2) using first-order through third-order concurrently monitored field dynamics.Reconstructions incorporated conventionally calculated k-coefficients and corrected k-coefficients using the proposed full correction strategy.Video S1.Video comparison of first-order and second-order corrections.Improved DWI and fractional anisotropy (FA) map quality are observed when using corrected first-order field dynamics compared with conventional first order.This is further improved when using corrected second-order fits.
Video S2.Video shows difference profiles in second-order fully corrected k-coefficients and conventionally calculated k-coefficients (black), coefficients corrected using only stepwise correction (red), and coefficients corrected using stepwise + weighted least squares (WLS) correction (blue), for concurrently monitored k-coefficients from non-diffusion-weighted (top row) and diffusion-weighted (bottom row) EPI acquisitions (R = 2).Non-DWI (b = 0 s/mm 2 , left), single-direction DWI (b = 1000 s/mm 2 , middle), and fractional anisotropy (FA) maps (right) reconstructed from EPI acquisitions using concurrently monitored third-order through first-order conventionally calculated and fully corrected k-coefficients, respectively.Fitting errors in EPI are shown to primarily translate as bulk image shifts in the phase-encoding direction, which lead to erroneous FA map quality.
Video S3.Video comparison of DWI and fractional anisotropy (FA) map quality when using conventional fitting, subsets, and full use of the correction algorithm when performing second-order fits of concurrently monitored field dynamics.Stepwise correction provides the most correction benefit, followed by incremental improvement with each successive correction.

Figure S7 .
(A) Non-DWI (b = 0 s/mm 2 , left), single-direction DWI (b = 1000 s/mm 2 , middle), and FA maps (right) reconstructed from single-shot spiral acquisitions using third-order concurrently monitored field dynamics.Shown is a comparison of image quality for images incorporating conventionally calculated (upper left quadrant), gradient nonlinearity (GNL) corrected (bottom left quadrant), fully corrected stepwise, weighted least squares (WLS), and weighted phase residual (WPR) correction (top right quadrant), and fully corrected plus GNL-corrected (bottom-right quadrant) k-coefficients.Difference images calculated between images in each column.(B) Average normalized RMS error (NRMSE) values between DWI and ground-truth phantom images for DWI reconstructed when including k-coefficients as calculated using the methods described in (A).By itself, GNL correction improves image quality, but not sufficiently.Applying GNL correction in addition to the conventional three-component correction algorithm provides no additional improvement in quality beyond what is offered by the three-component correction strategy.