Iterative static field map estimation for off-resonance correction in non-Cartesian susceptibility weighted imaging.

PURPOSE: Patient-induced inhomogeneities in the magnetic field cause distortions and blurring during acquisitions with long readouts such as in susceptibility-weighted imaging (SWI). Most correction methods require collecting an additional ΔB0$$ \Delta {\mathrm{B}}_0 $$ field map to remove these artifacts. THEORY: The static ΔB0$$ \Delta {\mathrm{B}}_0 $$ field map can be approximated with an acceptable error directly from a single echo acquisition in SWI. The main component of the observed phase is linearly related to ΔB0$$ \Delta {\mathrm{B}}_0 $$ and the echo time (TE), and the relative impact of non- ΔB0$$ \Delta {\mathrm{B}}_0 $$ terms becomes insignificant with TE$$ \mathrm{TE} $$ >20 ms at 3 T for a well-tuned system.
METHODS: The main step is to combine and unfold the multi-channel phase maps wrapped many times, and several competing algorithms are compared for this purpose. Four in vivo brain data sets collected using the recently proposed 3D spreading projection algorithm for rapid k-space sampling (SPARKLING) readouts are used to assess the proposed method.
RESULTS: The estimated 3D field maps generated with a 0.6 mm isotropic spatial resolution provide overall similar off-resonance corrections compared to reference corrections based on an external ΔB0$$ \Delta {\mathrm{B}}_0 $$ acquisitions, and even improved for 2 of 4 individuals. Although a small estimation error is expected, no aftermath was observed in the proposed corrections, whereas degradations were observed in the references.
CONCLUSION: A static ΔB0$$ \Delta {\mathrm{B}}_0 $$ field map estimation method was proposed to take advantage of acquisitions with long echo times, and outperformed the reference technique based on an external field map. The difference can be attributed to an inherent robustness to mismatches between volumes and external ΔB0$$ \Delta {\mathrm{B}}_0 $$ maps, and diverse other sources investigated.

1 INTRODUCTION

Susceptibility‐weighted imaging (SWI) is commonly used in high‐resolution brain venography or traumatic brain injuries for its sensitivity to blood products and calcium. The magnetic susceptibility of those compounds contributes to local field distortions visible on filtered phase maps. Long echo times TEs (e.g., ms at 3 T) are typically used to enhance the susceptibility contribution, resulting in compromises on slice thickness to reach clinically acceptable scan times. Many parallel imaging and compressed‐sensing (CS) methods , , , , , , have been proposed over the last 2 decades to accelerate MRI acquisitions and non‐Cartesian sampling patterns , have recently gained popularity. In particular, the trajectories based on the spreading projection algorithm for rapid k‐space sampling (SPARKLING) in 2D and 3D imaging , can reach acceleration factor (AF) superior to 15 in scan times compared to fully sampled Cartesian imaging in high resolution (0.6 mm) isotropic brain imaging by taking advantage of all degrees of freedom offered by modern MR scanners and long observation windows (i.e., readouts). However, this experimental setting also results in off‐resonance artifacts amplification that causes geometric distortions and image blurring. These artifacts are mainly caused by patient‐induced static field inhomogeneities that are especially pronounced near air‐tissue interfaces (e.g., close to the oral cavity and the ear canals). Non‐Cartesian sampling patterns tend to be more sensitive to inhomogeneities causing local k‐space inconsistencies over the different gradient directions. Therefore, specific attention must be paid to reduce signal and detail dropouts in these regions.

Different approaches exist to compensate for off‐resonance artifacts during the acquisitions or later on in image reconstruction or post‐processing. The default technique integrated in all scanners is the spherical harmonic shimming , that allows the partial mitigation of patient‐induced static inhomogeneities using a set of shim coils. A field map is acquired and minimized with respect to the gradients that are expanded over a spherical harmonic basis (generally up to the second or third order). This calibration step is performed quickly using a low‐resolution map because the technique does not allow precise correction. Although this approach provides a significant improvement on image quality, a complementary technique is required in the most demanding setups such as in EPI , or in non‐Cartesian imaging, that makes use of long readouts. , , , , Recent works have been proposed to improve this method using more efficient coil designs that accurately fit the expected inhomogeneity profiles using complementary data. Those experimental results permit to achieve improved correction performances, but remain limited by the number of shimming coils to be deployed to reach the theoretical limits. For that reason, alternative retrospective approaches are instrumental in correcting off‐resonance effects once the k‐space data has been collected.

In the particular case of Cartesian EPI, a well‐known technique , consists in inverting the gradient direction at every time frame to emphasize the geometric distortions and therefore, obtain a deformation map that can then be directly corrected. This technique is, however, difficult to transpose to non‐Cartesian readouts (e.g., spirals, rosette, SPARKLING) , , , , , , because the gradient directions are strongly varying along multiple dimensions simultaneously. In such context, a more generic and well‐established method , , , , consists in compensating the undesired gradient by modifying the Fourier operator to integrate the prior knowledge of a field map. This technique can be applied to both Cartesian and non‐Cartesian data at the cost of much longer (e.g., 15‐fold) image reconstruction times.

For multi‐echo sequences, the field map can be directly derived from the acquisitions. For single echo sequences, the constraint of knowing the field map has been considered from different perspectives in the literature , , , , , to avoid additional medium‐to‐high resolution multi‐echo acquisitions that would counterbalance the goal of accelerating the original scan. A first solution is to integrate the multiple echoes into the desired sequences, , however, they represent an additional constraint on the pulse sequence design that is not compatible with some protocols because of their long readouts. This holds in particular for multishot spiral, segmented EPI, or SPARKLING imaging. Post‐processing approaches have been developed, notably by the teams that contributed to the correction method. , One of the most popular approaches over the years initiated by Sutton et al is to solve a non‐convex optimization problem involving the extended signal equation model during image reconstruction. Recent developments , rely on tight regularizing constraints to make the solution locally smooth by projecting the gradients onto a set of functions (e.g., 2D Gaussians). Compared to the original implementation, this approach guarantees a successful convergence to a local minimizer. However, the field maps shown in Patzig et al are low resolution because of the involved smoothing constraints. Another possibility lies in the simulation of the field map , with mask‐sensitive binary anatomy models, however, proposing an accurate mask is challenging without even considering the model limitations. Because the stated non‐Cartesian SWI problem arises from long readouts centered over long TEs, in this paper, we propose to take advantage of the phase properties in such context. Indeed the contribution to the phase maps grows linearly with and becomes significantly dominating, whereas most other contributions remain constant or decay. , , The disproportion reaches a point where it can be assumed that the phase is equivalent to the field map with minimal error.

In this work, a pipeline is proposed to convert the phase map extracted either from a spherical stack of 2D SPARKLING trajectories or from full 3D SPARKLING ones into a high‐resolution 3D field map that is subsequently used to correct the SWI acquisition for inhomogeneities. These SPARKLING trajectories were specifically designed for SWI acquisitions with different acceleration factors [as defined later in Equation (10)]. The iterative procedure requires notably to address the problems of coil combination , , , , and phase unwrapping , , , for which 2 and 4 algorithms, respectively, were jointly confronted. On 4 healthy volunteers at 3 T, the field map estimates in the brain are directly compared to the individually matched acquired references and the benefit of this internal estimation is subsequently assessed on SWI brain images after offline image reconstruction.

2 THEORY

The simplest model to recover a complex‐valued MR image from an acquired complex‐valued signal over time window at the grid position from k‐space trajectory is:

The above equation corresponds to an adjoint Fourier transform of the raw data. However, this model does not account for the different perturbations of the signal that occur during acquisition. A common extension that allows the recovery of a partially corrected image for the undesired effects of field inhomogeneities reads as follows:

The additional gradient term within the exponential assumes the knowledge of the field map. The most common method to find it relies on the measurements of phase maps. Indeed, the observed phase of an MR image at position can be linked to the field map as follows:

and where , , and are the phase contributions emanating from major sources: inhomogeneities in the and fields for the first 2 terms, and field induced by eddy currents for the third one. The term includes the effect of minor sources such as heating and physiological motion. A common approximation is to consider related to the through a linear relationship , :

The common approach to estimate a map is to perform acquisitions with multiple TEs, and consider other phase contributions as constant with respect to . For example, considering 2 TEs such that < , we obtain:

This formulation implies that must be recovered from by solving Equation (3). Although many unwrapping algorithms have been proposed in the literature, , , , , convenient analytical choices are preferred to reduce the observed wrapping, and therefore, simplify or even avoid this step, such as applying a division over the complex‐valued volumes acquired with close TEs. In contrast, our approach consists in taking advantage of acquisitions with long TEs to propose a high‐resolution map based on the phase image from a single echo SWI acquisition through the following approximation: considering that other contributions are almost negligible for long TEs at high magnetic field for well‐tuned MR scanners. , , , , Indeed the absolute error in the previous approximation can be decomposed into individually measurable contributions:

The main component in Equation (8) is because it includes both the radiofrequency pulse inhomogeneities and the coil sensitivities. An optimal coil combination with respect to the signal‐to‐noise ratio can be computed when the coil sensitivity maps are known. Similarly to the map, these sensitivity maps require an additional acquisition and are therefore, usually unknown. Various methods exist to estimate them from the acquired multi‐channel volume. , , After recombination, a constant phase shift remains undetermined depending on the reference coil. This term can be simply removed if we assume that a 0‐order shimming is applied and matches the acquisition box because no constant term is expected any longer. Both coil combination algorithms and 0‐order shimming hypothesis are considered in later experiments.

The additional fields created by the radiofrequency pulse inhomogeneities and the eddy currents , , through the MR scanner and the coil are generally limited for well‐designed and calibrated hardware. , An experiment on oil phantom proposed in Supporting Information Section S1 suggests an error inferior to 2 Hz from the eddy currents contribution in the proposed experimental setup. However, the relative contribution of both sources increases with the field strength and may have a significant impact on the estimation at very high magnetic fields. , Overall, the related error at 3 T is expected to be inferior to 10 Hz , , on average, whereas off‐resonance frequencies can reach up to 300 Hz.

Based on the literature, other sources such as heating and physiological motion can be neglected for static contributions. Note that the contribution from heating follows a similar relationship with respect to but is still considered negligible:

where is the Larmor frequency and the frequency shift per degree Celsius, whereas refers to the temperature variation at voxel .

Stronger motion such as bulk motion during the sequence acquisition could greatly impede the estimation. This issue is considered to be addressed separately, because it would also hinder the overall quality of the images. Different methods exist to correct motion artifacts, , however, the current clinical practice is to simply repeat the exam. Note that avoiding the map acquisition already allows us to ignore potential inter‐acquisitions mismatches caused by such motions. Overall, the main contributions to the estimation error are the radiofrequency pulse inhomogeneities, coils sensitivity, and eddy currents, , along with an undetermined constant phase term supposedly close to 0.

Finally, the contribution from magnetic susceptibility of tissues has to be specifically preserved in the case of SWI. Indeed the susceptibility information, used to enhance diverse contrast sensitivities, originates from minor inhomogeneities , and therefore, the estimated field map should exclude those contributions to avoid canceling them during the off‐resonance correction. Because that information is usually extracted from high frequencies using filters, similar filters can be used over the field map to remove this information from the correction step.

3 METHODS

Our contribution is 3‐fold: first, different algorithms for coil combination and phase unwrapping is jointly studied. Second, the best combination is integrated into an iterative pipeline aiming to estimate the field map. Third, this map is used to perform off‐resonance corrections on the acquired k‐space data. The performance of the later is evaluated by comparisons with corrections obtained from external field maps.

3.1 Coil combination and phase unwrapping

A particular attention has been paid to the above‐mentioned critical steps, namely the coil combination and phase unwrapping. Multiple algorithms have been jointly evaluated to reach an accurate and reliable estimation.

In regards to the coil combination, the technique developed by Parker et al referred here to as virtual coil combination (VCC) appears as the state of the art among the phase‐preserving techniques. Another method introduced by Walsh et al, referred here to as adaptive coil combination (ACC), has been improved contemporarily by Inati et al to specifically preserve phase during combination. Both algorithms were internally implemented on graphical processing units (GPU) in Python using the CuPy (https://github.com/cupy/cupy) module, with specific modifications to the ACC method detailed in Supporting Information Section S2.

For phase unwrapping, 4 algorithms , , , were considered to compare different approaches. Only the method from Ghiglia et al was implemented internally on GPU, with noticeable changes described in Supporting Information Section S3 to extend it to 3D imaging. The other approaches were tested on CPU thanks to different open‐source implementations: scikit‐image (https://github.com/scikit‐image/scikit‐image) in Python for Herráez et al, the original implementation of ROMEO (https://github.com/korbinian90/ROMEO) in Julia by Dymerska et al, and PyMRT (https://github.com/norok2/pymrt) in Python for Schofield et al. The method from Ghiglia et al was the only one observed to benefit from masking the phase before application, so the others were not masked by default. Additionally, the method developed by Dymerska et al was used along with the magnitude information, as recommended by the authors.

A qualitative evaluation of the estimated maps by visual inspection after 1 iteration was performed over 4 volumes. The main criterion was to avoid anomalies in the regions and outside the skull where the phase is locally inconsistent because of a void signal. The second criterion was to obtain off‐resonance values as close as possible to the acquired map used as reference.

field map estimation

The map estimation pipeline is represented in Figure 1. First, the non‐Cartesian multi‐coil k‐space data was compressed using the principal component analysis‐based method proposed by Buehrer et al. The channel dimension is reduced to 30 components for image reconstruction, and to 5 components for the map estimation. An adjoint non‐uniform fast Fourier transform (NUFFT) with a pre‐computed density compensation was then applied to the 5 compressed channels to produce multi‐channel complex‐valued volumes. The different channels were combined using the ACC method and the resulting magnitude image was used to produce a brain mask, and then combined with magnitude and phase images.

FIGURE 1

Pipeline for iterative estimation followed by SWI reconstruction. The pipeline represents the overall process of estimation (A,B) and the reconstruction with inhomogeneity correction (C). A simplified version is displayed in the upper left corner for clarification purpose. In block A, the undersampled acquired k‐space is transformed into a complex‐valued volume, and then converted in block B into the required field map. Blocks A and B can be looped over for a few iterations to ensure a robust estimation (blue arrows) by adding residual estimations from each step (orange‐circled blue plus sign). The produced field map is then used along with the undersampled k‐space in block C to obtain the corrected SWI volume

At this stage, the observed phase remains wrapped over the (−, ] domain. Phase unwrapping was performed using the algorithm proposed in Ghiglia and Romero, which consists in solving a Poisson equation weighted by the previously computed magnitude image. Finally, the phase maps were scaled by following Equation (7) and a low‐pass filter was applied using a Hanning window over the central third of the Fourier domain. This filtering step is specific to SWI and is carried out to preserve the susceptibility contributions in the high frequencies during the correction. This point is detailed in Supporting Information Section S6. Other contrasts relying on phase information might also require specific map processing.

Therefore, a first estimation of the map is computed. For highly accelerated acquisitions (e.g., AF >15 as defined below in Equation [10]), the signal‐to‐noise ratio in the regions can be too low to provide an accurate field map. Therefore, the process should be iterated to improve the first estimation as shown in Figure 1. In the next iterations, the regular NUFFT operator is replaced by a correcting pseudo‐NUFFT operator, combined with the previous estimation. This step actually consists of interpolating the usually intractable gradient compensation through a weighted sum of regular NUFFT operators. The histogram‐based coefficients were used as weights with a number of interpolators to maintain an approximation root‐mean‐square error (RMSE) below . This way, some signal is recovered in regions where some of the inhomogeneities are being removed. The updated map estimate provides information about previously missed inhomogeneities, and is added to the previous one(s) until convergence. Hence, the resulting map is also corrected for geometric distortions caused by inhomogeneities. In the following experiments, 3 additional iterations were run to reach convergence, based on experiments available in Supporting Information Section S5.

3.3 Data acquisition and reconstruction

A total of 4 SWI volumes were acquired with non‐Cartesian 3D gradient echo sequences, each on a different healthy volunteer at 3 T (Magnetom PrismaFIT, Siemens Healthcare, Erlangen, Germany) with a 64‐channel head/neck coil array. Local and national ethical committees have approved the protocol, and written consent was obtained from the volunteers. Two different sampling patterns were used with distinct acceleration factors defined as:

with the in‐plane resolution, the number of slices and the number of spokes. The recently proposed 3D spherical stack of SPARKLING was used for 2 volumes (referred to as 1 and 2) with AF = 10, and its recent extension, namely full 3D SPARKLING, was considered for the other 2 (referred to as 3 and 4) with AF = 20. All acquisitions were performed with the following parameters: a 0.6 mm isotropic resolution, a field‐of‐view of 24 cm in‐plane ( = 384) over 12.5 cm , an observation time of ms centered around an ms, a repetition time ms. For the spherical stack of SPARKLING acquisitions (1, 2), we used a dwell time s, and a number of spokes = 8192 that resulted in an acquisition time of 5 min. For the full 3D SPARKLING (3, 4), a smaller dwell time s was used to balance the smaller number of spokes = 4096 that resulted in a shorter acquisition of 2 min 30 s. For comparison purposes, an additional reference map was acquired with a 2D gradient echo sequence using the following parameters: an acquisition time of 2 min 43 s (no acceleration), same field‐of‐view, a 2 mm isotropic resolution, 2 ms and ms. Those TEs allow for an excursion of 203 Hz, resulting in 1 phase wrap present in all references and unwrapped using the method proposed by Herráez et al.

MR image reconstruction was performed offline using the pysap‐mri (https://github.com/CEA‐COSMIC/pysap‐mri) software, , (Figure 1C) which implements 3D self‐calibrated compressed sensing reconstruction (‐norm regularization in the wavelet domain uses symlet 8 and 3 scales of decomposition). The correction was performed using the approximation described in Fessler et al, with the number of interpolators chosen as mentioned previously. Overall, the volumes were reconstructed considering 3 competing strategies: without correction, with correction based on the acquired field map used as a reference, and with correction based on the estimated field map. Further details about the offline reconstruction and correction are given in Supporting Information Section S4.

Finally, SWI specific processing as described in Haake et al was applied. The low frequencies were extracted by applying a Hanning window over the central third of the k‐space, before being removed from the phase image to obtain a high frequency map, subsequently normalized to produce a continuous mask. The magnitude image was multiplied 5 times by the mask, and a minimum‐intensity projection (mIP) was computed using a thickness of 8 mm. All post‐processing was run on a 2560 cores Quadro P5000 GPU and 16 GB of GDDR5 VRAM (NVIDIA, Santa Clara, California, USA).

4 RESULTS

4.1 Coil combination and phase unwrapping

Different algorithms were considered in the estimation pipeline for the coil combination and unwrapping tasks, and results are shown in Figures 2 and 3, respectively. The 2 coil combination algorithms compared in Figure 2, namely VCC (Figure 2B) and ACC (Figure 2C), obtain similar performances for AF = 10 (row 1), whereas an overall degradation can be observed at AF = 20 (row 2). The front region close to the sinuses in (Figure 2B2) with VCC is dominated by noise, making difficult to estimate even visually the unwrapped phase values. On the other hand, the ACC method produces locally consistent results (Figure 2C2) with only some blurring in the regions associated with the highest off‐resonance frequencies, therefore, obtaining more reliable phase maps.

FIGURE 2

Comparison between coil combination algorithms from 1‐step reconstruction. Two algorithms are compared one another to combine the 5 channels used during the map estimation process over volunteer 1 (1) with AF = 10 and volunteer 3 (2) with AF = 20. The magnitude (A), displayed to help identify the artifacts, was considered identical in both cases. The phases obtained using the VCC algorithm (B) and the ACC algorithm (C) are shown to compare the degradation of the information caused by field inhomogeneities in both cases

FIGURE 3

One‐step estimations using different coil combination and phase unwrapping methods. Different algorithms are compared one another to estimate the maps over 2 volumes from the phase maps shown in Figure 2: top rows (1,2) are obtained from volunteer 1 with AF = 10, and bottom rows (3, 4) are obtained from volunteer 3 with AF = 20. The virtual coil combination (1,3) and the adaptive coil combination (2,4) techniques are used to obtain the wrapped phase maps. Different unwrapping algorithms are then used in all 4 situations to produce maps (B)‐(E) close to those acquired 1 (A), Ghiglia et al (B), Herràez et al (C), Dymerska et al (D), and Schofield et al (E). Note that all field maps were masked post‐estimation

The phase maps from Figure 2 are then unwrapped using different algorithms and processed to yield the estimates shown in Figure 3, respectively, using VCC (rows 1,3) and ACC (rows 2,4) methods on both volunteers 1 (rows 1,2) and 3 (rows 3,4). Most unwrapping algorithms are consistent with the reference at AF = 10 (rows 1,2) except for Schofield et al's method (Figure 3E). However at AF = 20 (rows 3,4), no unwrapping algorithm managed to obtain results close to the reference maps in 1 step. Particularly, frequency values seem to drop significantly in key regions pointed on (Figure 3C3) and (Figure 3D3‐4), using respectively, Herráez et al's and Dymerska et al's methods. Such behavior can be detrimental to the additional iteration steps. For the other algorithms the field map remains locally smooth, however, its values are underestimated as illustrated in Figure 3B4,C4. An improved estimation is achieved in the inner region in Figure 3C4 with the method proposed by Herráez et al, balanced by extreme frequencies outside the skull on the left side, also observed in Figure 3C1,D1. Overall, the approach proposed by Ghiglia et al, shown in Figure 3B, is preferred because of the absence of detrimental patterns, particularly when combined with the ACC method , shown in rows (2,4).

The first estimation step described previously took ∼40 s to compute on average, with the following details: adjoint NUFFT computation in 16 s, coil combination in 5 s, phase unwrapping in 8 s, and masking, filtering, and scaling in 11 s. Additional steps may last between 3 and 5 min depending on the number of interpolators used for correction, which tends to increase at each step. Note that the correcting pseudo‐NUFFT operator can be parallelized with more resources to achieve a duration for additional steps close to 40 s as well. Overall, the estimation time contributes marginally to the 6 to 8 h required to perform the state‐of‐the‐art correction.

The intermediate estimation steps are represented in Figure 4. For both volumes, the low‐frequency phase is progressively canceled over the iterations, whereas the map converges. The resulting maps are shown for all volumes in Figure 5 along with the acquired reference maps. The phase images obtained after 1 correction step with the acquired maps are shown on Figure 5 row 4. For each volume, an important phase component remains, particularly in the bucco‐nasal regions. This observed phase is wrapped from negative to positive values, almost reaching an entire 2 cycle that corresponds to an inappropriate over‐correction of 50 Hz [Equation (5)]. A remaining phase component from other sources is indeed expected in Equation (4), but those regions are specifically associated with inhomogeneities and the values are too large to correspond to or contributions. On the other hand, the estimated maps produce spatially homogeneous phase images (see row 6 in Figure 5). More details are available in Supporting Information Section S7.

FIGURE 4

Iterative estimations using the selected pipeline. The 4‐iteration estimation process is illustrated from columns (A) to (D) over 2 volumes to show the evolution of their map estimates: top rows (1‐3) are obtained from volunteer 1 with AF = 10, and bottom rows (4‐6) are obtained from volunteer 3 with AF = 20. Rows 1 and 4 correspond to the magnitude images, rows 2 and 5 to the wrapped phase images, and rows 3 and 6 to the iterative estimated field maps obtained from (1,2) and (4,5), respectively

FIGURE 5

Comparison of the field maps and the resulting phase images. All 4 volumes are displayed from left to right in columns (A) to (D), respectively. Rows 1 and 2 correspond to the uncorrected magnitude and phase images, row 3 corresponds to the acquired field maps used to produce the corrected phase images (row 4) after 1 iteration, and row 5 corresponds to the estimated field maps used to produce the corrected phase images (row 6) after 1 iteration

A comparison between the estimated and acquired maps is proposed for all volumes in Figure 6. As shown in the superimposed histograms, the estimated frequencies tend to be lower than the acquired ones, therefore, confirming that references might be over‐estimated as observed in Figure 5. Additionally, this is quantitatively assessed through the linear regression slopes around 0.9. The distributions are, however, similar as shown by the Pearson correlation coefficients . The RMSE varies from 6.69 Hz for the first volume to 10.80 Hz for the third volume. A noticeable shift of 7.27 Hz is observed for the third volume (originating from the acquired map) although it has the highest value.

FIGURE 6

Comparison between acquired and estimated maps. The acquired and estimated field maps are confronted one another through superimposed histograms (A) and linear regressions (B) over the spherical stacks of SPARKLING acquisitions (1,2) with AF = 10 and the full 3D SPARKLING acquisitions (3,4) with AF = 20

A global specific absorption rate (SAR) of 3% was measured during examinations. Considering a homogeneous temperature rise, the heating contribution to the phase can be roughly estimated around 5 mHz using Equation (9) (with = 127.74 MHz, = −0.01 ppm/°C). This confirms that the temperature contribution to the phase variation can be neglected for the proposed sequence.

artifacts correction

The volumes with SWI processing are displayed in Figure 7 for AF = 10 and Figure 8 for AF = 20 with emphasis on the bucco‐nasal region. Overall, the less affected volume by the correction is 2, whereas the more impacted is 1. This striking difference cannot be explained by acquisition parameters, which were similar for both. Besides the main artifact, the shape of the anterior region for volume 1 has indeed become sharper. Stronger artifacts are observed in volumes 3 and 4 along with a larger correction. Quantitative comparisons are available in Supporting Information Section S7. Overall, the corrections are improving the image quality for all volumes, although some regions cannot be recovered.

FIGURE 7

Comparison of correction with acquired and estimated maps for AF = 10. The different volumes obtained from the spherical stack of SPARKLING acquisitions on volunteers 1 and 2 with AF = 10, SWI processing and a 8 mm minimum intensity projection (mIP) are displayed on the left column (A) and right column (B), respectively, according to the following conventions: without correction in row 1, corrected using acquired field maps in row 2, corrected using the estimated field maps in row 3. Different artifact and correction details are pointed out using the following color‐coded arrows: artifacted (red), mitigated (orange), and corrected (green). The dotted lines overlaid to the frontal lobe in the axial views of column 1 are used to project the edges from (A1) onto (A2) and (B3) for comparison purpose

FIGURE 8

Comparison of correction with acquired and estimated maps for AF = 20. The different volumes obtained from the full 3D SPARKLING acquisitions on volunteers 3 and 4 with AF = 20, SWI processing and a 8 mm minimum intensity projection (mIP) are displayed on the left column (A) and right column (B), respectively, according to the following conventions: without correction in row (1), corrected using acquired field maps in row 2, corrected using the estimated field maps in row (3). Different artifact and correction details are pointed out using the following color‐coded arrows: artifacted (red), mitigated (orange), and corrected (green)

Although the corrections obtained with the reference and proposed maps are similar, some noticeable differences persist. On volume 1 (Figure 7A), the right part of the brain visible on the axial and coronal views (Figure 7A2) suffers from degradation only present in the reference correction using the external map. The anterior region previously mentioned recovers a wider area shown with dotted lines when using the estimated map. Almost no difference is visible between the 2 corrections on volumes 2 and 4 except when comparing the axial slices (Figure 7B,C). Finally, both volumes 1 and 3 demonstrate a definite advantage of using the proposed field map estimation technique compared to using the external field map.

5 DISCUSSION

An iterative pipeline was established by comparing 2‐coil combination and 4 phase unwrapping algorithms to produce a robust field map estimate. The coil combination method from Walsh et al, further improved by Inati et al, along with the phase unwrapping method developed by Ghiglia et al and embedded in an iterative pipeline permit to obtain a stable estimation that empirically converged in few iterations. The resulting maps were observed to be highly correlated with the acquired maps used as reference. Improved corrections were generally observed with the proposed method over the references, whereas none of the expected sources of error seemed to impact the correction.

The differences observed between the acquired and estimated maps have a few distinct causes. The first one to consider is the quality of the acquired map. The number of TEs and their proximity in time can influence the measurements and therefore, the computation of the map using Equation (6). Although 2 close echoes facilitate the phase unwrapping task, as mentioned in the Theory section, a third and more distant echo can be used to strengthen the linear regression for peak values, which might otherwise be less accurate. , The residual phase shown in Figure 5D suggests an over‐correction and hence, that the acquired field maps are all over‐evaluated. Therefore, the difference observed between acquired and estimated field maps should be considered carefully. Additionally, a mismatch between the external and estimated field maps could also be caused by inter‐scan motion, but it would not explain the overall lower values in the proposed estimation. The resolution was also suspected to impact the correction, but experiments proposed in Supporting Information Section S6 showed almost no difference when downsampling the map estimates. Overall, the external field maps could be improved to provide a better quality assessment of the map estimation algorithm.

The second element to take into account is the modeling error involved in Equation (8). Besides the contributions from the and fields, the hypothesis of zero‐order shimming stated in the Theory section was not systematically respected. Indeed, as observed in Figure 6, the linear regressions yield non‐null constant frequency shifts. For 3 volumes, a variation of <2 Hz is observed, whereas almost 8 Hz are reached with volume 3. However, this difference does not seem to imbalance the correction, because this volume still presents clear improvements using the proposed estimation over the reference.

Finally, the last possible explanation for the observed difference is the methods used in the estimation pipeline. The strongly accelerated acquisitions (AF = 20) showed phase degradation near the bucco‐nasal region, but still converged to a similar quality of map estimates as those obtained at AF = 10. However, as pointed out in Robinson et al, the selected method by Ghiglia et al is a non‐exact method that removes some spatial components of the phase. It could cause the lower estimated frequencies observed in Figure 6, but would not explain that phase residuals are present in the references rather than in the proposed correction in Figure 5. An iterative solution allowed us to favor the robustness of the unwrapping algorithm over its accuracy, but methods that are more recent could be explored to obtain a better compromise in fewer iterations. ,

In terms of image quality, 2 volumes ( 1, 3) under study showed a clear advantage in using the proposed method. In contrast, for the other 2 (2, 4), almost no difference was observed with the references (Figures 7 and 8). Volume 4 represents the worst estimation case, because it also has the lowest correlation coefficient between the acquired and estimated maps (Figure 6). On volume 1, the degradation brought by the reference correction on the right side of the brain (Figure 7A2) could be explained by unexpected inter‐scan movements causing a mismatch between the map and the SWI volume. The proposed method overcomes such issue as the map is internally extracted from the same data set. Besides the partial correction with both the acquired and estimated maps, no spatial degradation that could be attributed to the pipeline or the and fields has been observed. The above‐mentioned sources of error seem less impactful than the gain issuing from high resolution, robustness to motion, and reduced acquisition time.

6 CONCLUSIONS

MR imaging contrasts based on long readouts tend to suffer from inhomogeneities and therefore, from signal dropout, especially in non‐Cartesian acquisitions used to reach shorter scan times. Moreover, long TEs such as in SWI facilitate the estimation of a field map that can be further used for artifact correction as proposed in this paper. The proposed approach actually relies on the internal phase measurements during SWI acquisition and on some tenable assumptions concerning its main contributions to produce a high resolution and robust‐to‐motion field map in a few minutes only. This approach can be used a posteriori during the image reconstruction step and can deliver equivalent to improved corrections compared to the reference, which requires itself additional scan time. The estimation duration could be reduced in the future by improving the different pipeline steps with more advanced coil combination or phase unwrapping algorithms. Last, future work could focus on integrating optimization‐based estimation methods to account for bulk motion, or on extending this work from static imaging to EPI, similarly to recent works by Dymerska et al.

CONFLICT OF INTEREST

Guillaume Daval‐Frérot, Aurélien Massire, Boris Mailhe and Mariappan Nadar are all employed by Siemens Healthcare.

Section S1 Phase contributions from eddy currents. Additional experiment proposed to estimate the phase contributions of eddy currents during the main experiments.

Section S2 Coil combination. Detail of the diverse modification of the ACC algorithm used during estimations and reconstructions.

Section S3 Phase unwrapping. Detail of the diverse modification of the phase unwrapping algorithm used during estimations and reconstructions.

Section S4 Image reconstruction and ΔB0 correction. Detailed reconstruction and correction techniques used for both estimation and reconstruction pipelines.

Section S5 Convergence of the ΔB0 estimation. Additional reconstructions and details showing the empirical convergence of the estimation algorithm.

Section S6 ΔB0 map resolution impact on estimation and correction. Additional experiment proposed to show the impact of reduced field map resolution post‐estimation on correction quality.

Section S7 Additional quantitative results on ΔB0 estimation and correction. Diverse additional results proposed to complete the main discussion.

Figure S1 Phase contributions from the field on oil phantom. The contribution from eddy currents was estimated from phase maps in 4 different contexts: without (A) or with (B) a pre‐pulse gradient of 10 mT/m, and with either a body coil array (1) or a 64‐channel head/neck coil array (2). The voxels of minimum (blue arrows) and maximum (red arrows) values of the phase are automatically pointed out

Figure S2 Evolution of 4 map estimations and corresponding corrections over 15 steps. The estimation (A B) and consequent correction (C,D) are monitored during the estimation process over 15 steps, using an average L2 norm and SSIM score, respectively. The evolution is first observed by comparing each step with the fixed initial steps (A,C), corresponding to 0‐filled maps for the estimation (A) and the non‐corrected volume for the correction (C). The relative progression is shown by comparing each step with the next one (B,D). The results are displayed for each volunteer and each score, and a vertical dotted line shows the number of steps used in our proposition.

Figure S3 Comparison of different volumes corrected using multiple map resolutions. The map estimated with a resolution of 0.6 mm isotropic from acquisition 3 and used to produce the 3D SWI volume (A) was downsampled to 4 mm, 2 mm, and 1 mm to obtain on row (1) the corrected volumes (C), (D) and (E), respectively. The original volume without correction (B) is shown for reference. The absolute value of the differences between (B‐E) and (A) are displayed on row (2), multiplied by a factor 5.

Figure S4 Phase images from corrected volumes with or without map filtering. The impact of filtering during the estimation is shown over the 5‐channel 1‐step volumes obtained during the estimation process for volunteer 3.

Figure S5 Bland–Altman plots comparing the estimated and acquired maps. For each volunteer (1 to 4), the differences between estimated and acquired maps are shown according to the average between both values, for only 1 every 10 voxels for clarity purpose. The average and standard deviation (SD) of the differences are shown with plain and dotted red lines, respectively.

Table S1 Scores for different SWI volumes corrected using multiple map resolutions. Different scores were computed to compare the SWI volume obtained from the map estimated with a resolution of 0.6 mm isotropic from acquisition 3 and the SWI volumes obtained by downsampling the same map to 4 mm, 2 mm, and 1 mm.

Table S2 Scores for the SWI volumes produced using estimated maps. Different scores were computed for each volunteer to compare the corrected SWI volume obtained from the estimated map, as described in the main experiments, and the SWI volumes without correction and corrected using an acquired map.

Click here for additional data file.