Iterative Restoration of the Fringe Phase (REFRASE) for QSM

In quantitative susceptibility mapping (QSM), reconstructed results can be critically biased by misinterpreted or missing phase data near the edges of the brain support originating from the non-local relationship between field and susceptibility. These data either have to be excluded or corrected before further processing can take place. To address this, our iterative restoration of the fringe phase (REFRASE) approach simultaneously enhances the accuracy of multi-echo phase data QSM maps and the extent of the area available for evaluation. Data loss caused by strong local phase gradients near the surface of the brain support is recovered within the original phase data using harmonic and dipole-based fields extrapolated from a robust support region toward an extended brain mask. Over several iterations, phase data are rectified prior to the application of further QSM processing steps. The concept is successfully validated on numerical phantoms and brain scans from a cohort of volunteers. The increased extent of the mask and improved numerical stability within the segmented globus pallidus confirm the efficacy of the presented method in comparison to traditional evaluation.


INTRODUCTION
Since emerging as a novel magnetic resonance imaging (MRI) contrast less than a decade ago, quantitative susceptibility mapping (QSM) has rapidly evolved in terms of evaluation strategies and applications, exploiting the magnetic properties of tissue within the brain based on the phase information (Schenck, 1996;Reichenbach, 2012;Duyn and Schenck, 2017;Eskreis-Winkler et al., 2017). The foundation of any successful phase processing scheme is, however, defined by the quality of the phase data and, as a direct consequence, by the extent of the exploitable data support (Elkady et al., 2016). Being an ill-conditioned, spatially convolved problem, QSM requires a continuous volume for evaluation. As three-dimensional deconvolution is the central step of QSM, even small areas of incorrect phase or field information at the rim of the volume of interest, typically defined by a mask, will strongly propagate errors into the entire susceptibility map. This can obscure or impair the contrast of small anatomic details (e.g., Fortier and Levesque, 2018).
Other than a small number of approaches that aim to directly evaluate the complex signal (e.g., Liu et al., 2013;Wang and Liu, 2015), most algorithms initially unwrap the acquired phase in both the space and time domains or only in the time domain (e.g., de Rochefort et al., 2010;Liu et al., 2011b;Schweser et al., 2011Schweser et al., , 2016Sun and Wilman, 2015). Unwrapping is, in general, only possible when a continuous phase function describes the true data. Furthermore, the signal-to-noise ratio (SNR) should be high enough to represent this function and the true phase change per voxel or time point must remain below π, fulfilling the Nyquist-Shannon criterion. Susceptibility contrast between brain tissue and surrounding areas, such as blood vessels, bone, or air, can create strong local field distortions (Schenck, 1996;Liu et al., 2015;Dixon, 2018;Fortier and Levesque, 2018), These distortions can generate undersampling in the time domain, whenever the relative field shift is strong enough, and in the spatial domain, when the phase gradients between voxels exceed the Nyquist criterion. Areas within the measured volume that are affected by any of the aforementioned impairments will be referred to as erroneous phase regions (EPRs) in the following and have to be excluded from the evaluation area (EA) by masking or a weighting step prior to data processing. The focus of this study is on EPRs that do not contain extreme susceptibility shifts or signal voids and are only undersampled due to strong gradients of external origin.
In healthy subjects, strong sources of field distortions, known as background fields, tend to reside outside of the brain and their influences on QSM are identified and addressed by tailor-made background field removal strategies (e.g., Shmueli et al., 2009;Liu et al., 2011a;Schweser et al., 2011;Sun and Wilman, 2013;Lindemeyer et al., 2015). Most background field removal approaches determine the nature of the background fields based on physical principles such as harmonic distortions (e.g., Schweser et al., 2011;Wu et al., 2012;Topfer et al., 2015;Özbay et al., 2017) or dipole fields (e.g., de Rochefort et al., 2010;Liu et al., 2011a). Only a few approaches (such as Neelavalli et al., 2009;Wharton et al., 2010) exist that quantify the value and position of actual dipole sources, rather than a pseudo-susceptibility distribution, as this step is not required for the correction of the EA. Some approaches, such as SHARP (Schweser et al., 2011), even reduce the spatial coverage of the EA, depending on the applied filter kernel. In a related approach, Wen et al. (2014) make use of iterative optimisation to improve the accuracy of the background field fitting. Laplacian approaches, such as that from Zhou et al. (2014), remove background influences and, at the same time, can potentially unwrap the raw phase. Nevertheless, their reliability is uncertain, especially when facing strong local contrast and undersampling (Fortier and Levesque, 2018). An extensive discussion on the underlying theory of various approaches for background field removal has been published by Schweser et al. (2017).
In contrast to correcting estimated field maps after unwrapping, our novel approach aims to retrospectively correct the phase within the EPRs in an iterative manner based on the obtained background field correction. Previously published iterative approaches, such as Buch et al. (2014), have successfully modelled air and bone contributions to improve the quality and the coverage of susceptibility maps. In contrast, we address arbitrary sources of field distortion without the need for discrete source localisation. The presented algorithm, REstoration of the FRinge phASE (REFRASE), makes use of the employed background field correction within the EA to obtain additional information on global background fields. In the present study, we apply this novel concept to the harmonic and dipole field corrections obtained with MUltistage BAckground FIeld REmoval (MUBAFIRE) (Lindemeyer et al., 2015). These corrections are modified to produce approximative background field maps for the EPR based only on data from within the EA. With this information, the measured phase can be iteratively corrected before further processing takes place, ideally recovering data from the EPRs and hence extending the overall EA.
The difference between field mapping applied to a dataset with and without prior phase correction is illustrated in Figure 1. Erroneous unwrap results are triggered when the measured phase is directly employed, biasing the estimate of the background field. Hence, the estimated local field strongly deviates from the true local field for pixels beyond the EPR (Figure 1, pixel 18). This makes evaluation of the local contrast (blurry black spot in the image), e.g., for QSM, impossible.
Since QSM maps are based on a non-local deconvolution of the corrected field map, structures near the surface of the EA and contrast in its centre suffer from partial data coverage (missing data) or regions with incorrectly estimated field values. By offering partial data recovery, our REFRASE method, first outlined in the form of a prototype in Lindemeyer et al. (2017), enables QSM with a greater coverage and enhances the reconstruction quality of QSM maps within the EA. Here, the feasibility and quality of REFRASE are assessed with numerical samples and via the homogeneity within segmented anatomical structures in a cohort of in vivo subjects. i]| ≤ π, describing the Nyquist limit, holds for every pair (i, i + 1). This inversion is commonly known as the unwrapping operation and will be denoted by ϕ

MATERIALS AND
For multidimensional datasets, this condition must hold for at least one path between two locations to render these mutually unwrappable. The set can only be entirely unwrapped if such a path exists for every pair of voxels within a (multidimensional) set of voxels. Depending on the unwrapping strategy, additional conditions might apply.

Data Masking
Whilst the true phase change is well-known in synthetic data, it remains unknown in measured data as it is concealed by the acquisition process. Generally speaking, even a measured phase dataset, in which phase wraps are spaced wider than the Nyquist limit of two voxels, can originate in a phase profile that was measured with significant undersampling, meaning that the actual phase wrap distance is much smaller than the limit (e.g., Dagher et al., 2014). Consequently, the resulting unwrap would significantly underestimate the true slope. An MRI measurement of a human subject acquired with reasonable parametrisation normally contains a region where unwrapping in time is possible and where all true phase gradients lie below the Nyquist limit. With increasing distance from this region, which typically resides in the geometric centre of the imaged volume and near the magnet's isocentre, field, and thus phase homogeneity, decreases due to global and local perturbations. Additionally, the SNR drops rapidly within strong field perturbations or when leaving the support of a continuous object.
Data fidelity has to be assessed based on the observed phase. If the slope of a function exceeds the Nyquist limit at a certain point, its gradient will first approach π/Voxel, leading to a strong dephasing of the local neighbourhood of the observed point at position [ijk] in the spatial domain. This can be evaluated by the so-called local coherence (LC) metric as described by Witoszynskyj et al. (2009): The presence of noise will further diminish Q LC , making the LC a good estimator for data validity and hence for delineating the EA. In order to avoid non-random noise bias in the near neighbourhood of 27 voxels, the LC can be smoothed, e.g., by a narrow-width Gaussian, e.g., with σ = 2 voxels, prior to thresholding. The thresholding limit, Q LC ≥ Q LC,min , is initially determined heuristically. In order to achieve reasonable masking, one has to consider the employed static magnetic field strength, voxel dimensioning, sequence parametrisation, and the evaluated echo time. Furthermore, one must also consider acquisition-specific artefacts such as fold-ins, regions of low signal or phase poles. For a given experimental setup, Q LC,min can be optimised for best result fidelity based on processing a cohort of measurements with several parametrisations. TE be a function along a discrete path in space, indexed by i = 1, . . . , n. Without loss of generality, we neglect phase offsets in this description (⇒ ϕ[i, 0] : = 0). Let further b be composed of a dominant function b bg and a local characteristic b local , with 2πt

Field Separation
Omitting the path index, this gives the commonly used separation (e.g., Deistung et al., 2017) It is clear that: While b local can be unwrapped in a local neighbourhood, b and b bg may exceed the Nyquist limit. Further,φ(t TE ) = mod 2π (2πt TE · b) is the measured phase andb bg is an estimate for b bg . The local phase estimate can be derived from the difference: . (4) Now,φ local can be unwrapped, leading the field mapb local : Reinserting the previously subtracted background field estimate, the true field distribution is recovered. Including (2), (4), and (5), this leads to the comprehensive term: A region can be entirely unwrapped ifb bg can be estimated everywhere within and if each pair of member points can be connected via a discrete path r[i], indexed by i = [1, . . . , n], that lies inside the region.

REFRASE Optimisation Concept
Equation (6) only holds in regions where the argument in the numerator fulfils the Nyquist criterion, i.e., whereb bg is a sufficient approximation. The b bg within EPRs can only be estimated based on background field removal within the EA, where data fidelity is provided by prior phase coherence thresholding. This motivates the use of an iterative correction scheme working to grow the reliable EA in each step. For phase data,φ, of a typical multiple-echo gradient-echo measurement, let m max be the maximum achievable brain mask and m EA ⊆ m max be a valid EA that can be fully unwrapped. A contiguous EA can be segmented via the process of determining the local phase coherence Q LC as in Equation (1) and through the removal of disjoint regions. To compensate for the low phase SNR at early echo times, a Gaussian filter (σ = 2 voxels) is applied to the first echo before field mapping. Subsequently, the field map calculated fromφ is expected to be valid within m EA and allows the background field to be estimated.
Following these steps, harmonic and dipole-shaped fields are determined by MUBAFIRE (Lindemeyer et al., 2015). The harmonic result is composed from a superposition of orthonormalised solid spherical harmonic functions. As theb SSH is only dependent on the harmonic order and the coefficients, it can be calculated within the expanded region m max . Dipole fields are usually estimated via a minimisation term within the EA and originate from its complement, not(m EA ) (e.g., de Rochefort et al., 2010). An additional condition, preventing field sources within m max , is introduced, resulting in the minimisation expression: where B 0 is the static magnetic field, d is the dipole kernel (de Rochefort et al., 2010), χ ext is the external susceptibility distribution and * is the spatial convolution operation. The Tikhonov regularisation suppresses susceptibility sources inside m max . In order to avoid strong streaking artefacts at the volume surfaces, morphologic dilation can be used such that m max extends m EA by at least one voxel. Using the previously determined χ ext , the externally generated dipole background field for m max ,b DIP is calculated usingb Summing up both components, the background phase derived fromb bg =b SSH +b DIP is subtracted from the original phase using Equation (4). The processing pipeline is illustrated in In order to avoid any bias in the results due to additional magnitude conditioning, spatial priors are not used. With each step, the raw phase is extended by the part of the EPR where Q LC fulfils the threshold condition for "healthy data" following background phase subtraction (see Figures 6, 10). In accordance with the convergence observations made in the first numeric test runs, a total of five iterations were used.

QSM Processing Without REFRASE Optimisation
The following steps were performed for conventional processing:

Simulation
A simplified numerical model, containing a sphere-shaped m max that mimics the observable brain support, is used to provide a proof-of-principle for REFRASE at ultra-high field (7 T). The generated phantom with dimensions of 128 × 128 × 128 voxels represents a simplified human head and neck and contains airfilled cavities and a blood-filled bubble for the demonstration of strong field distortions. Furthermore, bone is introduced as a spherical shell inside the head structure. Typical susceptibility values, roughly adapted to values from Schenck (1996), Hopkins and Wehrli (1997), and Spees et al. (2001), were chosen for all tissue types: χ Air = 0.36 ppm, χ Skull = −0.9 ppm, χ BloodBubble = −0.7 ppm, χ Tissue = −9.0 ppm, including three spherical-shaped regions in the central brain imitating deep grey matter, χ R 1 = χ Tissue + 0.2 ppm (hereafter serving as the control region), χ R 2 = χ Tissue + 0.25 ppm, and χ R 3 = χ Tissue + 0.3 ppm.
The field map was simulated based on fifth order solid spherical harmonic functions (SSHs) using randomised coefficients. Gaussian random values with zero mean were assigned with standard deviation ranges of 1, 1, 2.5 · 10 −3 , 1.25 · 10 −4 , 1.25 · 10 −7 , and 1.25 · 10 −8 for the zeroth to the fifth order, respectively. Field distortions generated by the numeric susceptibility model were included using dipole convolution (Marques and Bowtell, 2005) With a simplified, uniform T 2 = 80 ms, T * 2 was calculated as: within the brain, and 0 elsewhere. Based on the simulated field and the T * 2 -distribution, the true phase and the noise-free measurement signal were calculated based on and Measurement conditions were simulated by adding complex Gaussian noise with a relative standard deviation of σ noise = 0.01 to the signal, resulting inS(t TE ). The signal was simulated for t TE = [4, 16, 28, 40, 52] ms, with M 0 set to 1 wherever non-air structures were present. Figure 3 illustrates the individual steps of the simulation procedure. Conventional evaluation and REFRASE analysis with five iteration steps was conducted on the phase data obtained using a mask that contains the maximum extent of the simulated brain structure. This was repeated for different configurations of the Q LC,min threshold between 0.4 and 0.9 applied to the phase of the second echo. Susceptibility reconstruction was performed for each iteration step and each configuration using λ = 0.03 (Tikhonov) and µ = 0.001 (gradient regularisation). These regularisation values were determined heuristically, based on manual optimisation of image contrast vs. artefact level.
The resulting masks were then compared to the full brain mask of the numerical model, representing the maximum achievable coverage via the relative difference in mask voxel count: The mean susceptibility value, χ(R 1 ), within control region R 1 was determined in the reconstructed maps and its heterogeneity was assessed via the observed standard deviation, σ χ (R 1 ). A set of 10 numerical instances was simulated in a Monte Carlo-like approach and analysed. Each instance featured randomised parameters such as: positions of the air-filled cavities, the blood-bubble and differing coefficients of the harmonic background field characteristics.

Measurements
A 3D multi-echo gradient echo sequence with monopolar readout, parametrised for two different studies, was acquired at 1 mm isotropic resolution, FA = 11 • and TR = 18 ms with a spatially-selective pulse. Having obtained written informed consent, measurements on a cohort of ten volunteers were conducted using a 4 T MRI scanner with a single-channel birdcage coil. In total, three echoes were recorded at 2.5, 6.0, and 12.0 ms (with BW = [500, 500, 160] Hz/Px respectively), but only the first and last echo were processed in this study. A maximum coverage brain mask, m max , was derived using bet2 (FSL Toolkit, Smith, 2002). For each volunteer, bilateral regions within the globus pallidus of both hemispheres were manually segmented via ITKSnap Yushkevich et al. (2006) based on the susceptibility maps obtained. The segmentation is illustrated in Figure 4. All data underwent conventional evaluation as well as REFRASE analysis with five iterations. Regularisation parameters for QSM reconstruction were identical to the simulation. This was repeated for Q LC,min threshold values between 0.4 and 0.9, applied to the phase of the second echo. Here, n rel was calculated with reference to m max , and σ χ and χ were determined within the control region.
To demonstrate the merits of REFRASE one in vivo dataset was processed with STI Suite 3 (University of California, Berkeley; https://people.eecs.berkeley.edu/~chunlei.liu/software. html). V-SHARP (Li et al., 2011) was selected for background removal and StarQSM (Wei et al., 2015) was selected for susceptibility reconstruction. Three different configurations were computed: (a) using a classic FSL-bet2 brain mask, (b) using an eroded FSL-bet2 brain mask, and (c) using the REFRASEpreprocessed phase and mask. Figure 5 shows a sample slice of the susceptibility maps obtained using conventional and REFRASE processing. The conventional case shown in the top left, and particularly the first iteration step with the three lowest Q LC,min , clearly demonstrates how faulty or uncorrected field areas close to the rim of the mask distort the susceptibility result of the entire map, causing streaking artefacts. The homogeneity of the map increases with the iteration count and with higher threshold levels Q LC,min . However, it can be seen that very high thresholds result in reduced mask coverage. Figure 6 demonstrates how the rim of the mask changes between conventional evaluation and each REFRASE iteration. Whereas the conventional mode naturally has full mask coverage, REFRASE approaches this solution with each iteration, while excluding noise-contaminated areas.

Simulation
The behaviour of the relative mask difference and the susceptibility results are illustrated in Figure 7. The first row shows n rel . For all Monte Carlo cases, it can be seen that it significantly increases in the first iteration step but then decreases and approaches a constant level closer to zero, e.g., full coverage, with every further step. For Q LC,min ≥ 0.7, the achievable constant difference minimum grows rapidly with the threshold level, implying a lower achievable coverage of m max . The absolute standard deviation within the control region R 1 , plotted in the second row, approaches a comparable minimum for all cases. However, a finite minimum is to be expected due to the artificially introduced noise level σ noise . The convergence only seems comparably slow for very low Q LC,min = 0.4 and the standard deviation in iteration five lies slightly above the results for the higher Q LC,min . The relative standard deviation in comparison to the conventional approach, shown in the third row, supports this observation and makes the behaviour clearer for all Monte Carlo cases. With the exception of the first iteration of a few Monte Carlo cases, the standard deviation lies below that of the conventional analysis. As the heterogeneity should only be biased by noise, REFRASE generates the more precise estimate. The last row, illustrating the mean susceptibility value within R 1 , shows that the true value of χ(R 1 ) = 0.02 ppm is well-approximated for Q LC,min ≥ 0.6. For Q LC,min = 0.5, the convergence is slightly slower, and for Q LC,min = 0.4, the control region R 1 is not well approximated after five iterations. Therefore, threshold values between 0.6 and 0.7 are a good trade-off between mask coverage and quality.
In order to gain a global insight into the distribution of reconstruction errors, an additional Root Mean Square Error (RMSE) analysis was performed for two different areas: Global RMSE, covering the full evaluation mask, m EA , and rim RMSE, defining a the rim of m EA within six voxel widths from the mask surface. Since the presented QSM analysis and the data are reference-less, global constant offsets between result and simulated ground truth maps have been subtracted prior to RMSE calculation.
The results of the RMSE analysis are displayed in Figure 8, showing a sample rim mask image and RMSE plots of the rim and global masks for all ten simulation cases. Figure 9 shows susceptibility maps for one volunteer in a sample slice computed with the conventional approach and REFRASE. As observed in the numerical sample, the measurement in the first iteration also clearly shows the influence of faulty field areas at rim of the mask, distorting the susceptibility result of the entire map and the image centre by causing streaking artefacts as well as over-and underestimated areas. The homogeneity of the result increases with iteration count and with Q LC,min , while the mask difference appears to decrease. Figure 10 demonstrates the changes near the rim of the mask between conventional evaluation and REFRASE with each iteration. Again, the results are comparable to the synthetic data. The REFRASE iterations approach m max with each iteration, while noise-contaminated areas remain excluded. The relative mask difference and the susceptibility results are illustrated in Figure 11. While n rel again exhibits a peak in the first iteration of REFRASE, it decreases more slowly with each iteration compared to the synthetic data. This is potentially due to the higher complexity of the data. Nevertheless, a similar trend is observed. As was also the case in the simulation, a more rapid increase of minimum achievable n rel appears from Q LC,min ≥ 0.7 onward. In all cases, the absolute, as well as the relative, standard deviation within the segmented basal ganglia approach a value that is clearly lower than the conventional analysis. Leaving aside the first iteration, better homogeneity in the control region is confirmed for the acquired measurements from the second iteration onwards. Slow convergence is only observed for Q LC,min = 0.4, similar to the numeric samples. The mean susceptibility, or rather the value it is converging to, naturally differs due to the individual physical properties of the volunteers. Even for the lowest Q LC , it shows convergence after four to five iterations. For Q LC,min ≥ 0.6, the mean value curve already reaches steady behaviour after the second or third iteration. Figure 12 shows the results of the QSM processing comparison with STI Suite and REFRASE using our custom QSM reconstruction. The extent of the EA and the improvement of heterogeneity within the control region is also shown.

DISCUSSION AND CONCLUSIONS
A comparison between the behaviour of the relative mask difference, mean susceptibility, and homogeneity of the simulation with the in vivo measurements shows that the employed numeric model is a simplified, yet adequate representation of disturbances observed within actual measurements.
The results of the simulation confirm the functionality of REFRASE by showing clear convergence toward the known susceptibility values as well as demonstrating improved performance for the intra-region homogeneity within the control region, compared to the conventional analysis. The RMSE analysis of the rim shows better performance for all settings after FIGURE 6 | Illustration of the evaluation mask for conventional (raw) processing and REFRASE analysis on the synthetic data (coronal slice). The purple frame indicates the masked region shown on the right-hand side. For REFRASE, a threshold of Q LC = 0.6 was chosen. Each coloured rim indicates the iteration step, shown in masked view on the right-hand side. The phase of the fourth echo is shown in the background.
five REFRASE iterations. However, while the global RMSE shows similar behaviour for strict thresholds, it increases significantly with iteration count for very weak thresholds (Q LC,min = 0.4 or 0.5). This was to be expected as very low thresholds inevitably include unwanted EPRs. These noise-contaminated parts cannot be corrected for by REFRASE since they are not described by the physical model used for QSM and, hence, propagate into the entire map during QSM calculation. This explains the poorer convergence and homogeneity within the control region, and, consequently, here the results are closer to the conventional analysis.
In contrast, a very high Q LC,min would provide excellent results within the central control regions, at the cost of substantial mask shrinkage. Nevertheless, even in this scenario, a convergent behaviour develops after a few iterations. Partial compensation for the mask shrinkage may be achieved by employing iteration counts much higher than those tested in the present study. However, results imply that when thresholding is too strict, regions with steep phase gradients will be irrevocably excluded. In summary, the highest Q LC,min values are not suitable for high mask coverage and should only be chosen if the rim of the mask provides extremely poor phase quality and if the evaluation focusses on structures in the centre of the brain. Very low settings FIGURE 7 | Susceptibility and mask results for the simulation. One individual colour is used for each Monte Carlo case. From the top row to bottom: relative mask difference between true mask and phase masking result (n rel ); heterogeneity (standard deviation, σ χ ) of susceptibility within reference area; heterogeneity normalised to conventional result without REFRASE; mean susceptibility value inside of reference area (χ). Each graph shows the conventional calculation without REFRASE as 0th point and the iteration result of REFRASE.  of Q LC,min ≤ 0.5, on the other hand, are even less desirable as they tend to increase the reconstruction error.
The measurement-based results relating to mask difference and the intra-region homogeneity of the susceptibility maps obtained confirm the threshold values estimated in the simulation. This implies that even under realistic conditions, the proposed REFRASE scheme successfully improves stability within the control region and mask coverage of the results. The standard deviation in the control region has its lower limit at the tissue-specific heterogeneity. Thus, considering the physical principles of the processing steps applied, the decrease of σ χ with iteration count is not due to a loss of contrast. Rather, it has to be FIGURE 10 | Illustration of the evaluation mask for conventional (raw) processing and REFRASE analysis on the measurement data (sagittal slice). The purple frame indicates the masked region shown on the right-hand side. For REFRASE, a threshold of Q LC,min = 0.6 was chosen. Each coloured rim indicates the iteration step, shown in masked view on the right-hand side. The phase of the second echo is shown in the background.
implied by the larger spatial data support and by the exclusion of voxels with false field information, which would otherwise lead to artefact-induced heterogeneity.
As most other commonly used QSM background removal strategies at least partially shrink the mask to avoid edge effects, the correction of EPRs by conventional processing is not possible. This leads to either artefact contamination or, in case the EPRs are excluded in advance, implies a massive reduction of the m EA . While neither method can revert noise-induced phase falsification, Nyquist-exceeding phase gradients can be successfully recovered by using REFRASE (see Figure 10). This offers a valuable advantage, not only for dealing with typical artefacts generated by e.g., the nasal cavities but also for QSM applications in patients with brain tumours or traumatic injuries, especially near the surface of the brain. In both cases, bleeding or local aggregation of blood generate a susceptibility interface due to the iron content in the blood, which is much stronger than any usual intra-brain contrast and makes traditional evaluation in the surrounding areas difficult. While REFRASE is not able to recover the bleeding itself, it makes imaging the surrounding areas possible. However, care must be taken when applying REFRASE in subjects with bleeding inside of the brain rather than at the rim. Depending on the chosen parametrisation, such regions might be excluded by REFRASE, because they exhibit FIGURE 11 | Susceptibility and mask results from in vivo measurements. One individual colour is used for each subject. From the top row to bottom: relative mask difference, n rel , between a mask from the classic approach and the REFRASE result; heterogeneity (standard deviation, σ χ ) of susceptibility within reference area; heterogeneity relative to conventional result without REFRASE; mean susceptibility value (χ) inside reference area. Each graph shows the conventional calculation without REFRASE as 0th point and the iteration result of REFRASE. strong local phase gradients. Consequently, the investigation of a segmentation or preprocessing step to separate such regions from the brain support and from the exterior is an objective for further development.
To the best of our knowledge, only Topfer et al. (2015) and Buch et al. (2014) have presented rim phase correction approaches remotely comparable to REFRASE. Topfer et al. (2015) reported recovery limitations of one or two voxels closest to the mask edge as a result of employing first and second order derivatives. Even though our approach is not limited by a discrete kernel size or derivatives at a voxel level, the minimisation for an external susceptibility distribution, as employed in MUBAFIRE, causes overcompensation near the very edges of the mask and hence, REFRASE can, in its current form, also be affected by overcompensation in single voxels close to the outer rim. This issue was also reported by Liu et al. (2011a) in the original publication of the projection onto dipole fields (PDF) approach. The authors proposed the use of a dilated brain mask as W T , pushing the external susceptibility distribution further away from the boundary. Since this is identical to the REFRASE strategy, by using an EA smaller than m max , the scenario is much less likely to occur than in traditional PDF. While Buch et al. (2014) successfully demonstrated the modelling of distinct anatomic structures, such as the sinuses, and recovery of their effect on neighbouring tissue, the present approach represents a less source-specific, yet, more general means to recover phase data. A very recent publication from Buch et al. (2019) discusses the estimation of phase and susceptibility in the dural sinuses. Among other features, this approach utilises Taylor expansions of the phase to overcome signal voids. This nicely complements the rim phase recovery abilities of REFRASE by continuing the phase recovery beyond the brain surface.
To exclude additional bias introduced by the magnitude and to emphasise the advantages of REFRASE in a non-preconditioned scenario, spatial priors were deliberately excluded from all QSM analyses in this study. As a result, the QSM maps obtained are not as homogeneous and rich in contrast as common showcase maps with enforced contrast. Furthermore, as acquisition was via a single-channel birdcage, improvements in SNR by per-channel phase computation were not possible in the present study. However, our results also imply improvements for spatial-prior based QSM reconstruction.
The applied background identification steps, namely harmonic fields and external dipole sources from MUBAFIRE, were shown to perform very well on the processed datasets. Nevertheless, the newly introduced REFRASE principle can be applied with any other background removal strategy, as long as its field estimate can be estimated in or extrapolated to the maximum mask, m max . Alternatively, REFRASE-preprocessed phase can be fed into other QSM pipelines.
Finally, the QSM-Comparison showcase in Figure 12 outlines the importance of the phase quality at the mask rim for QSM reconstruction. Even a sophisticated QSM reconstruction performed on the inside of a maximised brain mask without restoration or exclusion of the rim phase, as in method (a), shows a heavy impact on the resulting QSM maps, even in central regions, as shown in subplot (f). While the homogeneity of the maps and within the control region generated with method (b) (eroded mask and V-SHARP/StarQSM) is significantly improved, the extent of the mask is clearly reduced. Evidently, method (c) (STI Suite with REFRASE preprocessing) and REFRASE with custom QSM produce the most robust maps with low control region heterogeneity, while keeping the mask coverage as high as possible.
It is important to note that the presented method used phase data from a birdcage coil and does not address the recombination of phase from multichannel receive arrays. However, REFRASE can be easily applied directly after an appropriate preprocessing step that combines the channels.
In conclusion, we have demonstrated a novel iterative data preparation method, REFRASE, which achieves significant improvements in terms of the quality of the susceptibility maps obtained and mask coverage for QSM. Both factors are significant for clinicians, as the quantitative precision of QSM is increased and the mask coverage becomes more comparable to that of other standard MRI contrasts. Furthermore, as ultrahigh field MRI becomes more accessible, the method will evolve to its full potential, compensating for stronger local magnetic field gradients. Finally, the presented method is very valuable for low-resolution field mapping approaches that suffer from strong intra-voxel gradients. The free choice of the backgroundestimation method renders REFRASE a universal tool for arbitrary phase contrast applications.

DATA AVAILABILITY STATEMENT
The in vivo datasets for this manuscript are not publicly available due to ethical restrictions. Requests to access the datasets should be directed to the corresponding author: Prof. N. Jon Shah. The segmented susceptibility maps of the Monte Carlo set as well as the according result maps of conventional and REFRASE processing are available on FigShare (https://doi.org/10.6084/m9. figshare.12769853). Further, for the first Monte Carlo case, the maps of all simulation steps are included.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethics Committee of the Medical Faculty of the Rheinisch-Westfaelische Technische Hochschule Aachen (RWTH Aachen University, University Hospital, Aachen, Germany). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
JL conceived the presented idea, developed the theory and performed the computations, and wrote the manuscript. WW, AS, and JL planned, designed, and performed the measurements. NS provided equipment and infrastructure. JL, WW, AS, and NS discussed the method and results. All authors contributed to the final manuscript.