Recommended Implementation of Quantitative Susceptibility Mapping for Clinical Research in The Brain: A Consensus of the ISMRM Electro-Magnetic Tissue Properties Study Group

This article provides recommendations for implementing quantitative susceptibility mapping (QSM) for clinical brain research. It is a consensus of the ISMRM Electro-Magnetic Tissue Properties Study Group. While QSM technical development continues to advance rapidly, the current QSM methods have been demonstrated to be repeatable and reproducible for generating quantitative tissue magnetic susceptibility maps in the brain. However, the many QSM approaches available give rise to the need in the neuroimaging community for guidelines on implementation. This article describes relevant considerations and provides specific implementation recommendations for all steps in QSM data acquisition, processing, analysis, and presentation in scientific publications. We recommend that data be acquired using a monopolar 3D multi-echo GRE sequence, that phase images be saved and exported in DICOM format and unwrapped using an exact unwrapping approach. Multi-echo images should be combined before background removal, and a brain mask created using a brain extraction tool with the incorporation of phase-quality-based masking. Background fields should be removed within the brain mask using a technique based on SHARP or PDF, and the optimization approach to dipole inversion should be employed with a sparsity-based regularization. Susceptibility values should be measured relative to a specified reference, including the common reference region of whole brain as a region of interest in the analysis, and QSM results should be reported with – as a minimum – the acquisition and processing specifications listed in the last section of the article. These recommendations should facilitate clinical QSM research and lead to increased harmonization in data acquisition, analysis, and reporting.

While QSM techniques and their biomedical applications have been extensively reviewed 3,6,7,47-56,57(p31), 58,59 , and the QSM research community has held two challenges attempting to identify the best susceptibility calculation algorithm 60- 63 , there has been no community consensus or white paper on how to best perform brain QSM in the clinical setting. Some vendors have started implementing research-level QSM pipelines for their systems, nevertheless QSM is not yet a product on most magnetic resonance imaging (MRI) systems. On the other hand, demands for a robust "end-to-end" QSM recipe covering acquisition and processing through to presentation in scientific publications have been growing strongly from the neurological imaging community. For example, the need for a QSM consensus method was highlighted at the group discussion of the Recommendations represent a consensus among QSM experts and the community of experienced QSM users at the time of publication. For clarity, the committee decided that the emphasis of the paper would be on 3 tesla (T), the field strength most widely used in clinical brain research, and that general guidance would be provided for other field strengths (1.5 T and 7 T) where possible. Due to the focus on clinical applications, including clinical trials, the recommendations prioritize robustness and simplicity over state-of-the-art acquisitions and processing algorithms, and provide a starting point for application-specific improvements. Areas in which the committee could not yet reach a consensus (e.g., due to insufficient evidence) are indicated as such in the paper.

Organization
The article is organized into eight sections, covering data acquisition, image processing, image analysis, and presentation of QSM studies in scientific publications (Fig.1). Data acquisition is split into two sections: 1) pulse sequences and protocol, and 2) coil combination, saving, and exporting; image processing is split into four sections corresponding to four processing steps: phase unwrapping and echo combination, creation of masks, background field removal, and dipole inversion; the image analysis section focuses on the analysis of susceptibility maps; and the last section on presentation covers reporting in scientific publications. The contributions of the members of the QSM Consensus Organization Committee are listed in the Supplementary Materials I, Section S1.1 along with an overview of the history and approach of the initiative (Section S1.3). The Acknowledgement section lists all individuals who contributed significantly to the consensus recommendations in verbal or written form.
The paper is accompanied by harmonized pulse sequence protocols for several major platforms (current software versions in 2022) as well as sample code to perform the recommended processing steps (see Supplementary Materials II) using the Sepia toolbox 73 [available at https://github.com/kschan0214/sepia]. In addition, curated online resources are provided and will be updated as the field progresses. This information is detailed in the Supplementary Materials.

Background
Tissue magnetic susceptibility is directly related to tissue chemical composition. A substance gaining a magnetization opposing an externally applied magnetic field, e.g. calcium, is called diamagnetic. Diamagnetism is due to electron orbit perturbation. A substance gaining a magnetization in the same direction as the external field, e.g. iron, is called paramagnetic.
Paramagnetism is due to spin alignment of unpaired electrons, which have a magnetic moment 658 times greater than that of the hydrogen nucleus. Most biological tissues, like the human brain, have a susceptibility close to that of water, which is diamagnetic with a value of roughly -9 parts per million (ppm). In the brain, iron deposition makes tissue less diamagnetic than water, and calcium deposition and the presence of myelin make tissue more diamagnetic than water 4,11,15,[74][75][76] .
The static magnetic field of magnetic resonance imaging (MRI) scanners, 0 , affects electrons in tissue and, in linear materials, induces a magnetization proportional to 0 with the proportionality constant defined as magnetic susceptibility of the tissue. The induced tissue magnetization generates a field that is nonlocal and extends into the space surrounding the magnetization source, thereby inducing field inhomogeneities. These inhomogeneities can be described mathematically as a convolution of the susceptibility distribution with the unit dipole kernel [77][78][79] .
During MRI signal generation, water proton spins experience this tissue field. In particular, gradient-recalled echo (GRE) imaging is very sensitive to this tissue field. Consequently, there is complex destructive interfering or dephasing by the tissue field variation within a voxel seen as hypointensities or blooming susceptibility artifacts on the GRE magnitude images, which is known as T2*-weighted or susceptibility-weighted imaging [80][81][82] . The GRE phase images, traditionally discarded, represent the tissue field inhomogeneities averaged over the voxel. The process by which tissue magnetic susceptibility is computed from the magnetic field measurement is called quantitative susceptibility mapping (QSM) 83,84 . Several detailed review papers have been published about the technical foundations and clinical applications of QSM, and we refer the interested reader to these articles for a comprehensive overview of the technique 5

Pulse Sequences and Protocol
This section describes recommendations to robustly acquire data for QSM.

Overview
Most of the major MRI manufacturers, if not all of them, provide a radiofrequency (RF)-spoiled 3D multi-echo GRE pulse sequence from which one can obtain phase images for QSM in addition to the standard T2*-weighted magnitude images. The GRE sequence is probably the most elementary sequence in the MRI sequence tree, consisting in its simplest form of an RF excitation pulse followed by acquisition of a gradient recalled echo, which for the multi-echo variant is repeated at various echo times (TEs) before the next excitation pulse is applied. This sequence is often used for calibration steps, such as obtaining a field map for 0 shimming.
Optimizing a GRE sequence to maximize QSM contrast shares many of the principles for optimizing T2*-weighted contrast-to-noise ratio or T2* mapping 94 . We recommend the following principles for designing a QSM acquisition protocol: -Aim to set the last TE equal to at least the T2* value of the tissue of interest. For example, if targeting deep gray matter, the T2* of the putamen can be a good reference (55, 30 and 16 ms at 1.5, 3 and 7 T, respectively 95 ). The first echo time (TE1) should be as short as possible and the spacing between consecutive echoes (ΔTE) should be uniform throughout the echo train.
-Use the minimum repetition time (TR) and set the flip angle to the Ernst angle ( Ernst = cos −1 ( −TR/T 1 )) for the target region. Most deep gray matter structures have a T1 close to that of white matter, and white matter represents the largest volume fraction in the brain.
Hence, the Ernst angle may be calculated for white matter or deep gray matter with T1 = 650, 850 and 1220 ms at 1.5, 3 and 7 T respectively 96 .
-Use three or more monopolar echoes. While QSM can be achieved with one echo and two is the minimum number of echoes needed to separate the intrinsic transmit RF phase from the magnetic field-induced phase (see the Section 3 below), the use of a larger number of echoes (e.g., 5 echoes) will benefit the phase signal-to-noise ratio (SNR) for a range of tissues 97 . As phase SNR is maximal when TE = T2*, the use of multiple echoes ensures that high SNR field estimates are obtained for both short apparent T2* (venous blood or tissues close to air-tissue interfaces) and other tissues with longer T2* values.
-Use the minimum readout bandwidth (BW) that generates acceptable distortions. At 3T, 220 Hz/pixel is often sufficient (two-pixel fat-water shift). Such acquisitions negate the need to use fat suppression for brain applications. Note that as the BW is increased the ratio of the readout duration (1/BW) and echo spacing decreases (as a consequence of the required time to rewind the gradients), which results in a reduction in the SNR of the acquisition. The best trade-off depends on the system's gradient amplitude and slew rate.
-Use isotropic voxels of at most 1 mm to avoid susceptibility underestimation reported to occur at larger voxel sizes 98 .
-Use 3D acquisition instead of 2D acquisition to avoid potential slice-to-slice phase discontinuities in 2D phase maps that mainly occur in interleaved acquisitions and require additional processing or to avoid slice crosstalk in non-interleaved acquisitions 99 . Note that the 3D spatial coherence of the obtained phase/field maps is of great importance as QSM is ultimately a 3D spatial deconvolution process.
-Use a monopolar gradient readout (fly-back) to avoid geometric mismatch and eddy current related phase problems between even and odd echoes in bipolar acquisitions since these require additional correction in the pipeline used 100 .
-Consider using flow compensation when targeting vessels (e.g., oxygenation studies), but note that flow compensation is often only available and guaranteed for the first echo, while flow artifacts increase in later echoes. Flow compensation also reduces the number of echoes achievable as it increases the minimum TE, assuming a fixed BW. In comparison, flow compensation has a much smaller effect on the resulting QSM values than the choice of QSM processing pipeline 101 .
The following points focus on practical aspects of image orientation, FOV, and acceleration to achieve a whole brain data acquisition in approximately 6 minutes for clinical research on a 3 T system with as few as 8 receive channels using a standard 3D multi-echo RF-spoiled GRE sequence: -Patients should be scanned in the supine position with the head straight to minimize variability in white matter susceptibility related to myelin magnetic susceptibility anisotropy 102,103 and microstructural water compartmentation 104 .
-For whole brain acquisition (including cerebellum), it is recommended to use isotropic voxels 98 and prescribe the imaging volume with the readout direction along the anterior commissure -posterior commissure (AC-PC) line (oblique axial orientation), which reduces the number of phase encoding steps and, consequently, the total scan time 105 .
Setting the readout in this direction also restricts eye movement artifacts in the left-right direction, outside of the brain. Alternatively, sagittal acquisitions (readout along the headfoot direction) can be beneficial when the region of interest includes the brain stem but will require larger acceleration factors. The use of tilted imaging slab orientations (oblique axial or oblique sagittal) requires that the data header information (see Section 3 below) specifying the exact slab orientation is supplied to the subsequent processing pipeline, otherwise artifacts will occur in the resulting susceptibility map 106 .
-Parallel imaging and partial k-space coverage can be used to reduce acquisition time, as for any structural MRI acquisition. Below we give general guidelines for GRE brain imaging: -If accelerating in only one direction, use acceleration factors <3. If available, noninteger acceleration factors allow fine-tuning the acceleration factors to the available equipment.
-If available, elliptical k-space sampling (elliptical k-space shutter) and/or partial kspace acquisition (e.g., 7/8) with zero-filling can be used to reduce scan time but note that this can result in a reduction of the effective spatial resolution.
-If an RF coil with 32 or more channels is available, acceleration in both phase encoding directions by a factor of 2 (2x2) can be used, preferably with CAIPI patterns 107 .
-If available, similar or higher acceleration factors compared to conventional accelerated imaging may be obtained with incoherent sampling using compressed sensing reconstruction methods 108 or with deep learning reconstruction methods 109 .

Consensus Recommendation
The recommendation in this section is based on protocols currently available on commercial scanners that do not require a special research agreement with the scanner manufacturer: 2.a. 3D multi-echo RF-spoiled GRE with monopolar readout. Three or more echoes should be acquired and the TE range should include the T2* times of the target tissues.
2.c. Flip angle should be set to the Ernst angle for target tissues (e.g., white matter).
2.e. Resolution should be isotropic with a voxel edge length of at most 1 mm noninterpolated at 3T.
2.g. Use coil arrays with a large number of elements covering the whole brain.

Additional Considerations
Acquisitions for QSM reconstructions can be performed at all clinical field strengths with minimal image distortions using the protocol recommendations provided here, although higher field strengths provide several benefits. Susceptibility contrast and contrast-to-noise ratio (CNR) increase with the static magnetic field once relaxation time changes are considered 110 , and acquisition can be more time efficient at higher fields due to the shorter T2*. Our recommended protocols implicitly integrate a compensation for lower SNR at lower 0 through lower spatial resolution, while keeping the total acquisition times similar. While useful data can be obtained at all field strengths, finer anatomical details in certain applications, such as brain lesions in multiple sclerosis 111 , might be difficult to visualize at 1.5T.
Advanced sequences such as segmented 3D-EPI 112,113 or Wave-CAIPI 114 allow drastic decreases in acquisition time but are not available as commercial sequences on all vendors and may require a research agreement with the scanner manufacturer. Nevertheless, it should be noted that some of these methods may have constraints on the number of echoes that can be acquired and 3D-EPI sequences will have increased spatial distortions that may have to be corrected for and that may result in additional artifacts that subsequent processing steps such as background field removal need to account for.

Coil Combination, Saving, and Exporting
This section recommends effective and practical solutions for generating phase images that can be used for QSM and provides guidance on how to save and convert the format of the phase images in preparation for QSM analysis.

Overview
Modern MRI systems use phased array coils made up typically of 12, 32 or 64 elements, which provide higher SNR than a single birdcage coil and facilitate parallel imaging [115][116][117] . However, the images from individual coil elements are sensitive to only a part of the FOV and need to be combined to generate a single imagea process that requires consideration of the differences between the coil signals.
Neglecting phase wraps (see Section 4), the phase measured with a particular RF coil element c in a GRE sequence, , is given by: = B (TE)+ 0 . The first part, B (TE), is the phase shift caused by the deviation ΔB0 of the magnetic field from the uniform main magnetic field, 0 .
Neglecting non-linear effects (discussed in Section 9.1.6), B (TE) evolves linearly over time and is the term that is relevant for QSM. The second part, 0 , is the phase at TE=0, known as the phase offset (or initial phase) of coil element c. This contribution comprises effects that are common to all the RF receive coils in a phased array, such as the phase of the transmit RF field 1 + , the effects of tissue conductivity , gradient delays and eddy current effects, and contributions that are unique to each receive coil, such as the coil sensitivity. Coil-dependent phase offsets must be removed prior to a complex summation of the coil signals, as shown in Figure 2, to avoid destructive interference. Destructive interference leads to reduced SNR and unphysical phase wraps in regions of the image, which cause artifacts in QSM 118 . Complete destructive interference is often associated with phase wraps referred to as "open-ended fringe lines" or "phase singularities" (see Figure 3, left). These ill-behaved wraps do not represent isophase contours and cannot be fully removed by unwrapping (see Section 3 below).
All major 3T MR system vendors have effective solutions for generating phase images from RF array coils on their current software platforms. Most of these techniques remove individual coil sensitivities, which are estimated by referencing to a coil with a relatively homogenous sensitivity over the object (usually the body coil). The reference data is acquired in a separate, fast, automated measurement, and the coil sensitivity correction is carried out on complex data either in k-space 115 or image space 119 before extraction of the phase. To combine the signals from each coil element, we recommend using the available 'on-console' vendor solutions listed at the end of this section. These methods may not be available on systems running older software versions, in which case phase images can be reconstructed from 'raw' (k-space) data offline or, alternatively, phase images can be saved and exported separately for each coil element (e.g., as DICOM-format image data) to allow an appropriate coil combination method to be applied offline.
These offline solutions may require additional efforts in 1) dealing with the export and transfer of very large files or a large number of files, which can be problematic in a clinical setting and 2) obtaining research agreements with vendors for proprietary information on the raw data format and on special routines for performing the image reconstruction offline. Several offline possibilities It is important to make sure that phase images are scaled correctly and converted to the correct data type for subsequent analysis. The analysis pipeline supplied with this paper works with a range of data types and arbitrary phase scaling, but some phase unwrapping algorithms and the nonlinear complex fitting algorithm described in Section 4 require phase to be saved as floatingpoint numbers scaled between - and , and this rescaling and data type conversion needs to be performed by the user. Further, note that the sign convention is different for phase values on different vendor systems 125 and this needs to be checked and corrected where necessary so that relatively paramagnetic tissues (e.g., iron-rich deep-brain regions) have positive susceptibility values in the resulting susceptibility maps.

Consensus Recommendations
We provide recommendations for combining phase images from array coils and saving phase data for each of the major manufacturers, listing the software versions for which these solutions have been tested and whether a research agreement is required. Detailed step-by-step descriptions and solutions for older systems are provided in Supplementary Materials III.
3.a The recommended solution for saving phase images are, for • Canon: SPEEDER, a version of SENSE which is available from MPower version 2.3 onwards and allows phase images to be reconstructed through a vendor-provided service password.
• GE: ASSET, a SENSE-similar solution which reconstructs magnitude, phase, real, and imaginary images without a research key on platforms MR30 onward.
• Philips: SENSE, which provides well-combined phase images 116 without the need for a research key from software version 5 onwards.
• Siemens: "Adaptive-combined with prescan normalize" 126 , which is available from software version VE11 onward in the product GRE sequence.
• United Imaging: an inter-coil referencing and weighted correction approach which is available from software version v9 without the need for a research key.
3.b Exporting data: Data should be exported in (classic) DICOM format.
3.c Format conversion: if the analysis pipeline requires NIfTI data, DICOM data should be converted to NIfTI using DCM2NIIX.

Phase Unwrapping and Echo Combination
This section describes the methods used to resolve phase aliasing and calculate a field map from multi-echo GRE data.

Overview
MRI phase measurements are constrained to an interval of 2 and are, therefore, subject to phase wraps or phase aliasing artifacts, i.e., the measured phase = ( B (TE) + 0 ) mod 2 . Such phase wraps introduce a phase difference of an integer multiple of 2 between the measured phase, , and the true phase B (TE) + 0 . Phase wraps are usually visible as discontinuous phase jumps in the phase images ( Figs. 3 and 4). To a first-order approximation, B (TE) = 2 TE • Δ 0 , where γ is the proton gyromagnetic ratio. To obtain an accurate estimate of the field shift, Δ 0 , for QSM, both phase wraps and the phase offset 0 need to be removed from the measured phase, 54 . The coil combination methods recommended in the previous section remove the coil-specific contributions to the phase offset, 0 , but leave other non-B0-related contributions in 0 , common to all coils, which should be removed by the QSM processing pipeline. Note that the background field removal step (see Section 6 below) removes only harmonic fields (those which satisfy Laplace's equation) within the brain region and cannot completely remove 0 as it contains both harmonic and non-harmonic components 55 . Therefore, 0 must be explicitly removed for accurate QSM 127 .
Over the years, different phase unwrapping methods have been adapted, refined and applied to MR phase imaging; including time-domain unwrapping methods (with multi-echo acquisition) such as CAMPUS 128 and UMPIRE 129 , and spatial-domain unwrapping methods such as region-based PRELUDE 130 , SEGUE 131 , and SPUN 132 , path-based best-path unwrapping 133 and ROMEO 134 , and Laplacian unwrapping 135 . The Laplacian unwrapping method is robust and gives wrap-free phase results even with low SNR but can result in high-frequency errors that propagate into susceptibility maps that are hard to detect visually 136 . It is also noted that Laplacian unwrapping only gives an approximation of the underlying unwrapped phase, especially when using the commonly used Fourier-based implementation, while region-based and path-based unwrapping give quantitatively more accurate estimates of the unwrapped phase 54,134 . Regionbased and path-based methods are termed "exact unwrapping methods" below as they give the exact value of the unwrapped phase 134 . When comparing exact unwrapping methods to Laplacian unwrapping, unwrapping errors (e.g., in veins and hemorrhages) are observed to be smaller in the former, improving QSM quantification accuracy, e.g., for oxygenation estimation 137,138 .
Therefore, we recommend using exact unwrapping methods.
Multi-echo phase images, e.g., acquired using the recommended protocol (see Section 2), can be combined to achieve a more accurate estimate of the underlying field shift, Δ 0 , than can be obtained from single-echo phase images 97 . This is because combining multi-echo phase images can remove the phase offset contribution and give higher SNR in the estimated tissue field and susceptibility maps (Fig. 4). The optimum approach may depend on the application, but two echo combination methods have been widely used for QSMnonlinear complex data fitting 139 and weighted echo averaging 94 .
The nonlinear complex data fitting approach takes into account the Gaussian noise in the complex images 139 and estimates the field shift, Δ 0 , and phase offset, 0 , together as parameters from fitting the complex MR signal over multiple echoes, with the requirement of having acquired three or more echoes 139 . This approach usually needs spatial phase unwrapping to be performed after the fitting, i.e., on a scaled field-shift estimate, e.g., 2 • • 0 , with being the echo spacing, wrapped again between − and + , as it resolves phase wraps in the temporal dimension. Nonlinear complex data fitting is more robust than linear phase fitting against phase noise at long TEs and around large susceptibility sources, e.g., veins and hemorrhages.
If echoes are acquired over a useful range of TE (depending on T2* values of the tissues of interest, see Section 2), the weighted echo averaging approach gives higher SNR for estimating the field shift, 0 , than the complex data fitting approach, which estimates multiple parameters 94,140,141 . Unlike nonlinear complex data fitting, this approach needs the phase data at each TE to be unwrapped first. It also requires explicit removal of the phase offset through subtraction of the estimated phase offset from the phase measured at each TE (for example as in MCPC-3D-S and ASPIRE, where the phase offset can be estimated by extrapolating the linear phase evolution to zero echo time) 120 . Unwrapping errors at longer TEs in voxels with large field shifts and more pronounced noise can be reduced by the "template" unwrapping approach used in ROMEO, which performs path-based spatial unwrapping on an early echo and unwraps other echoes on the basis of the expected linear phase evolution 134 . This, combined with weighted echo averaging, reduces the effect of such errors in the estimated field map.
To improve QSM quality, the spatial noise map generated from nonlinear complex fitting 139 , or the phase "quality map" calculated in path-based unwrapping 134 can be used to mask out voxels with unreliable phase (see Section 5).
Most echo combination methods assume linear phase evolution over TE, ignoring nonlinear phase evolution due to microstructure-related compartmentalization effects or biased sampling of the sub-voxel field perturbation (see Section 9.1.6), flow, or signal dropout 142 .
Advanced modeling of the phase evolution over time may provide further information about tissue composition and microstructure 37,104,143-148 , but is beyond the scope of this paper.

4.a.
Use an exact phase unwrapping method.

4.b.
Perform echo combination before background field removal.

4.c.
The optimal pipeline for phase unwrapping and echo combination depends on the acquisition and application. We recommend using either nonlinear complex data fitting followed by spatial phase unwrapping, or weighted echo averaging after template phase unwrapping and explicit phase offset removal.

Creation of masks
This section provides recommendations on creating masks for background field removal (Section 6), dipole inversion (Section 7), and visualization (Section 9).

Overview
Masking is often overlooked when describing a QSM pipeline but is a crucial step 149 , particularly for background field removal (see Section 6). Masking refers to selecting a region of interest (ROI) within the whole field of view and applying a process or function only within this ROI. In QSM, field maps, Δ 0 , are masked primarily because most background field removal algorithms require a mask. In general, masks should cover the largest ROI possible to prevent exclusion of brain tissue with a sufficient signal-to-noise ratio to have reliable phase/field values. This is of special concern for studies of the cortex and the brainstem near the brain border or air-tissue interfaces.
Unreliable field map data is composed mostly of extremely noisy voxels resulting from phase noise in regions with very low MRI signal or rapid signal decay. The noise distribution in phase images (and hence, field maps) is generally non-Gaussian and depends on the local magnitude of the signal 150 . In practice, regions with very low MRI signal yield phase noise uniformly distributed throughout the whole -π to π range, obscuring any underlying phase contrast information 150 .
Masking is a binary operation. Voxels with mask values of 1 (or Boolean "true") are included in the selected ROI and voxels with mask values of zero (or Boolean "false") are excluded. Masks may be created by using heuristic thresholding operations on available subject images, including magnitude images, T2* maps, quality maps, or SNR maps. Masks created from differently thresholded images may also be joined or combined to exclude or include regions based on different criteria. In addition, segmentation algorithms may be based on pre-learned shapes or on the optimization of functionals [151][152][153][154][155][156][157] . In particular, the Brain Extraction Tool (BET) 158 from the FMRIB Software Library (FSL) is a widely-used method for brain masking (skull stripping), although it may fail when pathologies or injuries are present 159,160 . BET is a magnitude-image based algorithm that effectively removes non-brain tissues, air, and bone from magnitude images of the head. When acquiring multiple echoes, using the last-echo magnitude image for BET masking is robust to remove regions with rapid signal dropout 161 , which is undesirable if such regions are of interest. A more balanced approach with a larger ROI selection is achieved by using magnitude images combined across TEs (e.g., using sum of squares or weighted averaging). Alternatively, the magnitude image of the first echo can be used for brain extraction, with the use of a phase-quality map to further exclude voxels with unreliable phase values, as described below. Alternatives to BET include standard template-based brain-extraction 162 . Deep learning segmentation alternatives may also be considered, as this is a rapidly developing field [163][164][165][166] .
QSM is also vulnerable to errors and artifacts arising from unreliable phase data that may not be directly reflected in the corresponding magnitude data. These may be caused by coil combination errors, flow in vessels, and other factors. For this reason, it has been proposed to use phasebased quality maps in addition to magnitude-based masking to refine masking 161 . A straightforward method to obtain a phase quality map is to threshold the inverse of the noise map provided by the complex nonlinear multi-echo fitting algorithm (described in Section 4) at its mean value 139,161 . The thresholding at the mean value maintains an adequate number of voxels, as it is applied to the entire field of view (FOV), and the distribution of values exhibits bimodal characteristics. This approach effectively distinguishes between reliable and unreliable voxels, serving as a suitable initial approximation. However, modifying the threshold factor (e.g., by multiplying it by 1.2) may yield further enhancement in the results. Some exact phase unwrapping algorithms provide an alternative source of phase-based quality maps 134,167 . Setting a threshold for phase-based quality maps can help to identify voxels within the brain with unreliable phase values, and to provide a better estimation of the brain boundary 149,161,168 .
Most phase unwrapping and echo combination algorithms do not require masking 54 but suppressing extraneous data by masking can speed up some algorithms and improve their robustness. In contrast, almost all background field removal algorithms require masking to define the region of interest, outside which the susceptibility sources are classified as background field sources (see Section 6) 55 . Notable exceptions (i.e., background field methods that do not require masking) are recent deep learning approaches such as SHARQnet 169 and Total Field Inversion (which, however, requires a mask-like preconditioner) 170 . The performance of background field removal algorithms depends strongly on the mask and poor background field removal can negatively affect the quality of the reconstructed susceptibility maps 171,172 . Many dipole inversion algorithms use masks to exclude voxels with unreliable field values from the susceptibility computation or use masks for regularization 83,139,173 . Finally, susceptibility maps should be masked for display purposes to exclude streaks and spurious information outside the brain (see Section 7 on dipole inversion and Section 9 on presentation and publication).
Background field removal methods remove fields induced by all susceptibility sources outside the supplied brain mask. The accuracy of background field removal is lowest at the boundary of the ROI, such as the brain surface, and improves with increasing distance from the brain surface 55 .
Although further erosion of the mask after BFR is not explicitly required, it may be employed in specific cases to eliminate residual artifacts at the boundary. Also, many background removal algorithms are not able to recover a reliable local field over the whole ROI mask 55 , and further erosion is unavoidable. It should be noted that some background field removal algorithms (e.g., V-SHARP; see Section 6 below) will result in erosion of the brain mask and the resulting eroded mask should then be used in all subsequent operations and for display.
Holes inside the brain mask lead to the elimination of field contributions from the susceptibility sources within the holes during background field removal. Such holes can occur when thresholded (e.g., phase-based quality) maps are used to refine masking. For example, if a pathology (such as a hemorrhage or calcifications) creates unreliable phase data inside the brain, the affected region could be set to zero in the brain mask. Removing the field perturbations from susceptibility sources within the holes during the background elimination step will render these sources undetectable in the final susceptibility maps, which can be a significant problem in the clinical setting. Therefore, it is important that holes within the brain mask are filled before the mask is used for background field removal. Holes can be reintroduced in the mask used for dipole inversion as an effective way to prevent streaking artifacts from regions with unreliable phase values. This procedure has been included in some algorithms 139,174 . Since dipole inversion is a nonlocal operation, correct susceptibility values may be inferred inside relatively small holes excluded from the dipole inversion mask. To avoid streaking artifacts created by high contrast sources and pathologies, while preserving accurate susceptibility values inside the holes, some recent two-step approaches suggest performing reconstructions with and without holes and then merging both results 138,149,175,176 . This is useful to improve the accuracy of the reconstructions, and to characterize hemorrhages or calcifications.
Although some recent deep learning-based single-step QSM approaches have shown that explicit masking could be avoided 169,172,[177][178][179] , these methods require further study and validation to be considered for clinical applications.

Consensus Recommendations
The recommendations below are summarized in Figure 5, and differences between masks (Masks 1-4 in Fig. 5) are highlighted in Figure 6. 5.a. Create an initial brain mask (Mask 1) by applying a whole-brain segmentation tool (such as BET) to either the combined (sum of squares) or the first echo magnitude image. The goal of this initial mask is to remove air, skull and other tissues, while preserving cortical areas. Further refinement is performed in the following steps. 5.a. Create a mask of reliable phase values (Mask 2) by thresholding the phase quality map generated by the multi-echo combination method in Section 4. Multiply Mask 1 with Mask 2.
5.b. After multiplication, holes should be filled to obtain the mask to be used as an input to background field removal algorithms (Mask 3). 5.c. Holes from Mask 2 can be reintroduced to avoid streaking artifacts from unreliable phase data within the brain. For increased accuracy of susceptibility values inside pathological regions of low signal, e.g., hemorrhages and calcifications, mixing data from reconstructions with and without the holes can be performed. 5.d. The calculated susceptibility map should be multiplied by the mask used for background field removal (without holes; Mask 3) before display, reporting of susceptibility values or further analysis.

Background Field Removal
This section provides recommendations for the background field removal step.

Overview
In QSM, the background field is defined as the field generated by susceptibility sources outside a chosen ROI 55 , in our case the brain mask (see previous section). In the brain, the background fields are generated by the tissue and air surrounding the brain. The susceptibility difference between brain tissue and air is approximately 9 ppm 76 , which is almost two orders of magnitude larger than the naturally occurring susceptibility differences within the brain parenchyma.
Therefore, background fields can be significantly larger than the tissue field in the brain, but not always are. Certain pathologies, such as hemorrhages, can create a tissue field that is similar in magnitude locally. The term local field is also often used in the literature for fields generated by tissue within the ROI, but since the field is a nonlocal property, we use the term tissue field here.
Removal of the background field from the field map, ΔB0, allows focusing the inversion (see Section 7) on the spatial susceptibility variations located inside the ROI, which generate the socalled tissue field ΔBt ( Figure 7). When background fields are not completely removed from ΔB0, most dipole inversion methods will result in shadowing artifacts and/or experience a slow convergence rate.
Because of the spatial smoothness of the background field, spatial high-pass filtering has been a popular method to suppress background fields in the past. However, high-pass filtering also removes the low spatial frequency component, a major signal component, of the tissue field, which is not acceptable to QSM that requires quantitative accuracy of the corrected field maps.
Newer methods that directly exploit the harmonic function property of background fields have replaced heuristic filtering methods 55 . From Maxwell's equations, it can be derived that the background field is a harmonic field, i.e., it satisfies the Laplace equation within the ROI. A harmonic field is completely determined when it is known on the region boundary. In other words, the solution of the Laplace equation in a region with a given boundary condition is unique 55,180 .
The SHARP (Sophisticated Harmonic Artifact Reduction for Phase data) method 181 and variants thereof use the spherical mean value property of harmonic functions. This property implies that the average of a harmonic function over an arbitrary sphere centered at any location that fits within the region of interest is equal to the value of the harmonic function at that location. In practice, a radius is chosen that is somewhat large compared to the voxel size to overcome discretization effects. This means that the tissue field can only be computed for voxels that are at a distance The LBV (Laplacian Boundary Value) method 180 assumes that the field at the very boundary of the region of interest is entirely background field and determines a harmonic function that satisfies this boundary condition. An efficient algorithm has been introduced to solve this Laplacian boundary problem 180 . However, because LBV entirely relies on the field estimate at the mask boundary, it is sensitive to phase SNR at the boundary, rendering accurate masking particularly important. Often a small mask erosion is applied to remove low SNR voxels at the boundary.
While LBV can perform better than PDF and SHARP in some situations 55 , its performance has been observed to be highly dependent on the mask and on the quality of the field estimates at the mask boundary 189 .
Residual background fields can be suppressed by combining methods, such as applying additional polynomial fitting or V-SHARP after LBV. V-SHARP and PDF typically do not require additional polynomial fitting. B1 + related contributions in the field map, ΔB0 (see Section 4 above) are not removed by background field removal methods, but these are avoided when using multiecho data combined with field fitting as recommended above.
Susceptibility maps are dimensionless and are conventionally calculated and displayed in parts per million (ppm) (see Section 9.1.4). Assuming that the wrapped input phase was correctly scaled to (- to ) radians, the corresponding scaling can be done either before (on the tissue field map) or after (on the dipole inversion output) using the scale factor: ∆ (ppm) = ∆ (radians)•10 6 (radians•T −1 •s −1 )• 0 (T)•∆ (s) . When scaling is performed after, care has to be taken to adjust the default regularization parameter when using a total variation based dipole inversion method, as the regularization term scales linearly.

Consensus Recommendations
6.a. Use V-SHARP to achieve good results in many situations, as it is less sensitive to imperfections in brain masking. This comes at a cost of a one-voxel erosion of the brain mask used for dipole inversion (Mask 4 in Fig. 5) at the brain surface and reduced accuracy at the edge of the brain.
6.b. When whole brain mapping (including the cortex and superficial veins) is desired, use PDF. This method will be slightly more accurate throughout the brain. PDF requires a good brain mask.
6.c. Depending on the application, tissue field quality, i.e., the phase SNR especially near the boundary, must be balanced against mask erosion.

Additional Considerations
Single step 112,170,[190][191][192] or total field inversion methods fit the susceptibility directly to the field map, ΔB0, or even wrapped phase images. These are currently popular for applications of QSM outside the brain but are still under ongoing development to ensure robustness. Residual background field has been tackled for dipole inversion using weak harmonics modeling 136 . Finally, deep learning 52 has found application in background field removal as well and is the subject of ongoing development.

Dipole Inversion
This section provides recommendations for the field-to-susceptibility inversion step (Fig. 8), which derives from the tissue field map, ∆ (with background fields removed; see previous section), a map that is tissue magnetic susceptibility ( ) (up to a reference value 193,194 , see Section 8).

Overview
While susceptibility, ( ), is a local tissue property, the field is a summation of weighted contributions from the distribution of magnetic susceptibility in all space. Mathematically, this summation can be described as a convolution ( * ) of the susceptibility with the unit dipole kernel

77-79
: Convolution corresponds to multiplication in the spatial frequency domain, which facilitates its fast calculation and is used in most QSM implementations [77][78][79] to accelerate computations. The inversion step performs a deconvolution using the dipole kernel d(r), which reveals the local tissue susceptibility within the region of interest, ( ), from the background-corrected tissue field, ∆ ( ): ∆ ( ) * −1 ( ) = ( ). However, the dipole kernel value is zero at and very small near the cone surface of the magic angle (54. 7 0 ) relative to the direction of the main magnetic field, making this deconvolution a poorly conditioned inverse problem 75,84,[195][196][197] . The measured tissue field, ∆ , contains deviations from perfect dipole patterns, particularly in regions with small magnitude signal due to lack of water protons (field detectors) or rapid signal decay (largely caused by inhomogeneous fields). While the regions with most extreme deviations are usually eliminated through masking (see Section 5 above), remaining dipole deviations in the estimated field can cause deconvolution errors in the calculated susceptibility map reconstruction, manifesting as streaking and shadowing artifacts 53,198 . Additional information about the unknown susceptibility map, ( ), can be incorporated into the solution through regularization to suppress streaking and shadowing artifacts in the solution 25,59, 83,84,112,141,173,182,199,200 . An optimization approach for incorporating this additional information can be formulated according to Bayesian inference, which is the following minimization problem when approximating the noise in the field as Gaussian: Here the first term is the data fidelity term with spatially varying noise weighting and the second term, ( ), is the regularization term with as regularization strength 83,141 . The minimization problem is iteratively solved with the number of iterations determined by the desired convergence level. The optimal regularization strength ( ) depends on anatomy, susceptibility contrast, and SNR, and should be optimized to balance artifact suppression and image sharpness in each imaging protocol and application by varying using, e.g., the L-curve method 201,202 .
Various regularization strategies have been developed for the inverse problem in QSM 25,53,60, 83,84,136,141,173,182,192,197,[199][200][201][203][204][205][206][207] .  The simple sparse regularizer using L1 norm of the gradient (i.e., the total variation, TV) is a standard approach for brain QSM, as exemplified in MEDI, one of the most popular algorithms. where these methods yield better results than classical methods, the performance of these methods did not reach those of the best classical methods in the QSM challenges potentially due to generalization issues from limited training data 60,63 .
Some algorithms do not incorporate spatial constraints for suppressing streaking and shadowing artifacts but explicitly modify the dipole kernel instead 53,198 , for example, the thresholded k-space division 192,197,213 . Implicit regularizations based on the number of iterations may work, but these methods have limited denoising capabilities and may be less robust than the sparsity regularization optimization approach 59,205,214 .

Consensus Recommendation
7.a. Use an optimization approach for dipole inversion with a sparsity type regularization that is commonly used in compressed sensing 53 . Specific sparsity types include L1norm, total variation, and generalized total variation, which likely provide similar outcomes. Future algorithm developments and evaluations are needed to provide a more specific consensus on the sparsity type.

Additional Considerations
There may be streaking artifacts coming from strong susceptibility sources near borders and within the brain interior region. Major causes include the breakdown of the Gaussian noise assumption and other errors in the determined field 139 . These artifacts may be suppressed using methods such as masking out or reducing the weight of less trustworthy voxels in the optimization 139 . The border streaking can be removed by improving the brain mask 149,215,216 . The interior streaking can be reduced using techniques to improve convergence such as preconditioning 170,217 , and using in-painting techniques to compensate for field errors such as MERIT 139 and L1 data fidelity 174 .
There may be shadowing artifacts coming from residual background fields. This shadowing can be reduced by improving background field removal such as harmonic incompatibility removal 136,218 and by suppressing slowly varying spatial frequency components through regularization 219 or preconditioning 53,170 .

Analysis of Susceptibility Maps
This section provides recommendations for quality control and referencing of susceptibility maps, the quantification of susceptibility values, and the visualization of brain structures on susceptibility maps and derived contrasts in the context of clinical research performing group studies. The physical background for the consensus recommendations is briefly summarized. Possible tools to facilitate susceptibility quantification of brain structures and lesions as well as for group analyses are provided in Supplementary Materials I, Section S1.6. Figure 9 summarizes the recommendations of this section.  149,220 . Some of these tools may not be suitable for all possible QSM applications due to assumptions on patient cohorts of the implemented mask generation algorithms (see Section 5 above) or due to the need to adjust reconstruction parameters depending on the data (see Section 7 on dipole inversion). Most QSM applications still require multiple processing steps, which can result in error amplification/propagation or inconsistencies between steps, rendering QSM workflows prone to i) reconstruction artifacts (a list of common reconstruction artifacts is provided in Table 1)  and therefore 3D segmentation of reference regions is advisable to include a greater number of voxels. Consequently, whole-brain referencing is considered stable (largest possible mask) and reproducible (whole-brain mask readily available in all reconstruction pipelines).
The dependence of white matter apparent susceptibility on the fiber orientation with respect to the main magnetic field due to the geometry and complex microstructure of white matter fiber bundles 102,104,146,221,222 (see also Section 9.1.6) can be a source of additional variability when using a reference that includes white matter regions. Another challenge of referencing is that pathology or effects of age alter white and gray matter integrity, specifically myelination and brain iron levels, especially in deep gray nuclei [223][224][225][226][227][228] . In the case of widespread pathology when no ideal reference region exists, two reference regions could be used to evaluate if the choice of reference region affects the study results. If, for example, whole brain and CSF were used as reference leading to the same significant differences between patients and controls, the results can be assumes with greater confidence to originate from the presence of pathology, instead of being an artifact from susceptibility referencing 229 .
For local pathology, the use of contralateral or surrounding tissues as reference is an effective strategy to avoid introducing artificial susceptibility differences due to using a reference region affected by pathology. An accurate segmentation of ROIs is essential to uncover subtle changes in regional susceptibility values that might indicate pathology, or to establish normative values. While manual segmentation of regions by multiple expert readers is the gold standard for quantification of regional susceptibility values, this strategy is very time-consuming and therefore not feasible in larger studies. Many available automated neuroimaging segmentation tools are optimized for use with T1-weighted images or require T1-weighted input data 231 . However, when using these methods for the analysis of susceptibility maps, the segmentation and registration accuracy in many structures of interest (e.g., basal ganglia) can depend on T1 contrast 232 , which is also affected by tissue iron 233,234 , and the generally low visibility of some deep gray matter regions on T1-weighted images 235 (depending on sequence parameters). Consequently, ROI-based methods that rely solely on T1-weighted contrast may be biased and suffer from inaccuracies.
Previously, it has been shown that the use of a QSM or hybrid QSM-T1-weighted contrasts for template generation improves atlas and voxel-based analyses 30,31 . Therefore, using multicontrast segmentation can be considered the best approach to avoid template bias 236 . A list of recommended tools can be found in the Supplementary Materials I, Section S1.6. Partial volume effects might strongly affect susceptibility quantification both for voxel-based and ROI-based analyses, especially for small structures with relatively high susceptibility values such as veins 237 . This could be corrected for by eroding of ROIs 101 , only using high susceptibility voxels (in case of positive susceptibility) 238 or using a partial volume map for correction 239 .

Consensus Recommendations
8.a. When ROIs are affected by artifacts, exclude data by automated detection of outliers or outlier regions, use of image quality measures or visual inspection.

8.b. Ensure that analysis methods do not include voxels of the susceptibility map with
unreliable values, e.g., that lie outside of the eroded background field removal mask (see Section 5 above; Mask 4 without holes in Fig. 5).
8.c. Always reference susceptibility maps to an internal reference region before performing further analyses.
8.d. When choosing a reference region, consider the study design, influence of pathology, how pathology could bias the study findings and discuss accordingly. For widespread pathology, cross-checking results using two different reference regions (e.g., whole brain and CSF) can be considered safe to exclude bias.
8.f. Always include commonly used reference regions in your analysis and report mean and standard deviation in these regions along with those in other ROIs. responses collected for that item, respectively. The reporting of specific items was considered essential if there was unanimous consensus in reporting them among the authors ( = 1). Items that were not considered essential were assigned a "traffic light ranking" (green for 0.75 ≤ ≤ 1, orange for 0.5 ≤ ≤ 0.75, and red for ≤ 0.5). Standardized tables are provided to facilitate the reporting of a broad set of items. For the purpose of availability, unless there is limited space for particular journals, we recommend that items with S>0.5 be reported as other investigators may need this information. It is considered essential to report these items if they vary within the same study (e.g., if different scanners or different software releases are used within the same study).
Potential limitations and confounds should always be discussed. The last part of this section reviews some important aspects that should always be considered when interpreting and presenting QSM findings in scientific papers.

Acquisition hardware
Ideally, the acquisition hardware is described in one sentence reporting the scanner field strength, model, vendor, software release version, and type of coil used (including the number of channels). Table 3 provides an overview of consensus recommendations pertaining to acquisition hardware.

Acquisition sequence type and parameters
The QSM Consensus Organization Committee considered it as essential to indicate the acquisition sequence type (e.g. GRE, as recommended in this paper, or EPI; specify if the sequence is 3D or 2D) and several acquisition parameters including number of echoes, TEs, TR, flip angle, bandwidth, resolution and scan duration. Table 4 provides an overview of consensus recommendations pertaining to acquisition sequence type and parameters.

Reconstruction pipeline and analysis
It is considered essential to describe the toolbox and reconstruction pipeline and list the algorithms used. The numerical values of parameters used should be listed, even if they were the default parameters. Table 5

Sample paragraphs
Representative paragraphs describing data acquisition, processing and QSM calculation in a scientific paper are provided in the following. The description refers to the images shown in Figure   8. and have been referenced to the average susceptibility in the brain mask (imposed to zero by the adopted processing pipeline)." A complete description of the QSM calculation pipeline can be found in the Supplementary Materials II.

Interpretation of results
Potential limitations and confounds related to QSM should always be taken into account.
A potential confound that can affect the extraction of quantitative susceptibility values from MRI phase/frequency data arises from the fact that, even for a uniform voxel-averaged susceptibility distribution, the apparent field measured in a voxel depends on the subvoxel distribution and visibility of water protons (the sensors of the MRI signal) around susceptibility perturbers such as iron and myelin. This can lead to a phase shift resulting from the biased sampling of the fields generated by perturbers when sensor and perturber distributions spatially correlate and such correlation is anisotropic 243,244 . An example is water in and around myelinated fibers, whose anisotropic distribution leads to a fiber orientation-dependent shift in apparent frequency which can exceed 10 ppb 104,144,146,221,243,244 . In addition, substantial T1, MT, or T2* weighting may differentially affect water visibility in the various water compartments and render apparent frequency shifts dependent (in a non-linear manner) on TR or TE 142,143,145 . Especially above 3T, QSM values within and around fibers that run perpendicular to B0 should be interpreted with caution. For example, pathological changes in myelin structure but not myelin content in such fiber bundles may lead to QSM changes without actual changes in tissue susceptibility.
Another point of caution with interpretation of QSM are inaccuracies near the edge of the regions selected for the analysis (see Section 5). A notable example are areas near the surface of the brain, where phase data is unreliable (due to, e.g., the prevalence of paramagnetic blood in pial veins), or unavailable (due to the lack of signal in skull), or the tissue phase was partially removed in the background field removal step. Because of this, QSM values in some of cortical grey matter may be incorrect or have reduced spatial contrast and resolution.
Lastly, it should be realized that, when strong regularizations or prior information are used in QSM dipole inversion, potential smoothing and spatial resolution loss may occur or new features may be added 201,245,246 (see Section 7.1). Some anatomical detail, visible in phase or magnitude GRE data, may therefore be lost in the QSM.

Consensus Recommendations
9.a. Always report at least the essential information regarding the acquisition hardware (Table 3), acquisition sequence type and parameters (Table 4), reconstruction pipeline and analysis (Table 5). In summary, we recommend that data be acquired using a monopolar 3D multi-echo GRE sequence, that phase images be saved and exported in DICOM format and unwrapped using a quantitative approach. Echoes should be combined before background removal, and a brain mask created using a brain extraction tool with the incorporation of phase-quality based masking.
Background fields within the brain mask should be removed using a SHARP-based or PDF technique and the optimization approach to dipole inversion should be employed with a sparsity type regularization. Susceptibility values should be measured relative to a specified reference, including the common reference region of whole brain as a region of interest in the analysis, and QSM results should be reported withas a minimumthe acquisition and processing specifications listed in the final section.
The recommended steps for data acquisition, data preparation and post processing are intended to provide a uniform robust reference starting point for a brain-focused QSM study performed with a clinical scanner. Specialty applications such as the depiction of small structures might require spatial resolutions higher than recommended 247 . In this regard, limitations and further considerations are included in each section, but thorough testing of the processing pipeline is recommended before starting a large patient study.
We hope that the recommendations here will enable many medical research centers to perform comparable QSM studies on scanners from different vendors, and that the standardized acquisition protocols and the processing pipeline provided along with this article will facilitate these studies (see Supplementary Materials I, Section S1.5. and Supplementary Materials II). As more clinical QSM studies are performed, analyzed, and presented in scientific publications, and current and future technical innovations become mature, these QSM recommendations will need to be updated. CRediT authorship contribution statement  Figure 1: This consensus paper comprises eight sections. The first two sections cover image acquisition: 1) pulse sequences and protocol and 2) coil combination, saving, and exporting. The next four sections cover image processing: 3) phase unwrapping & echo combination, 4) creation of masks, 5) background field removal, and 6) dipole inversion. The last two sections cover analysis and presentation in scientific publications: 7) analysis of susceptibility maps, and 8) presentation and publication. The image output from each section is further detailed in the corresponding section.

Figure 2:
Steps in coil combination, saving, and exporting (illustrated for eight example coils from the 64-channel array). Each of the coils generates a phase image (left), which is modified by the coil sensitivity and other terms which make up the initial phase. The initial phase is removed using methods detailed in the text and referenced publications (center left) and phase images are combined in the manufacturer's reconstruction and saved for export in DICOM format (center right). QSM analysis software may require the DICOM data to be converted, offline, to NIfTI format (right). Images shown were acquired at 3 T (Siemens Healthineers, Erlangen, Germany; Prisma Fit, VE11C) with a head/neck 64 channel coil and the recommended multi-echo GRE sequence (TE1=5.25ms; echo spacing=5.83ms; 5 echoes) using monopolar readout. Additional details are reported in Section 9.1.5. The imaging data and a description of the acquisition protocol may be found in Supplementary Materials II.   Figure 2). More phase wraps can be observed at a later echo (bottom). B) Example unwrapped phase images (using ROMEO template phase unwrapping with MCPC3D-S phase offset correction). C) Total field estimation after echo combination using weighted echo averaging.  : By using the mask or reliable phase (Mask 2), the initial BET mask (Mask 1) can be further improved by removing unreliable phase data near the boundary (red). Mask 3 is used for background field removal (BFR). After BFR, Mask 3 may need to be further eroded depending on the output of the background field removal algorithms (eroded region shown in blue). This is used for visualization of the results and reporting. Unreliable phase data inside the brain can also be masked out for dipole inversion (holes in green, with the final Mask 4 used for dipole inversion in white).

Figure 7.
Process of background field removal estimates the background field component of the total field (first row; same images as shown on the right-hand side of Figure 4; unit is Hz) relative to a chosen region of interest (brain mask, third row) and subtracts it from the total field, resulting in the tissue field (fourth row; unit is Hz). The tissue field encodes the spatially varying susceptibility within the brain but is much smaller than the background field. This is illustrated by showing cross-sections (indicated by the dotted lines in the field images) in the total field (second row) and the tissue field (last row). The background field was calculated using the V-SHARP method. Figure 8. The process of dipole field inversion starts from the tissue field (first row, same images as shown in the bottom row of Figure 7; unit is Hz) and estimates the susceptibility map (second row; unit is ppm).

Figure 9:
Schematic for susceptibility map analysis in case of a study interested in susceptibility values of the putamen (blue ROI on the susceptibility map, ) and globus pallidus (red ROI). Data with streaking artifacts that affect the ROIs need to be excluded (or recalculated when applicable to all study data). ROI generation benefits from the inclusion of susceptibility contrast, e.g., by calculation of hybrid images (blue) or use of T1-weighted and susceptibility data (green). Susceptibility maps need to be referenced, then regional average susceptibility values ( ) can be computed from referenced susceptibility maps ( ). The shown susceptibility map without artifacts is the same as the one in Figure 8. Tables   Table 1. Reconstruction artifacts, possible sources, and strategies to identify, mitigate these artifacts and criteria to exclude the data.

Streaking and shadowing artifacts
Incorrect susceptibility values  · orientation dependence · might be affected by pathology, e.g. demyelination, gliosis, hemorrhage, atrophy, focal lesions · Relatively small region whole brain 229,236,236,257 · no extra mask required, brain mask from previous processing steps can be used · intrinsic for some methods · large region · might be affected by pathology and age (e.g. myelination, global demyelination, gliosis, iron accumulation, hemorrhage) · due to large WM fraction similar limitations as "white matter" above. Magnetic susceptibility values should always be reported in either ppm or ppb (parts-per-billion) and the reference region (see the Section 8) should be explicitly stated, even in the case the adopted method did implicit whole brain referencing.
When the reference region used in the study is not the whole-brain mask, its [mean ± std] susceptibility value when referenced to the whole-brain mask should be reported, to enable post-hoc rereferencing for meta-analyses. Generally, it should be discussed in the Discussion section how potential pathological changes within the reference region may have biased the study outcome.  , and -separation, which is an advanced susceptibility mapping method for susceptibility source separation (Shin et al., 2021).  group. In addition to the study group endorsement, the committee distributed an online form to provide the opportunity of endorsement to EMTP SG non-members, listed in Section S1.4 below.

S1.4 Individuals Endorsing the Consensus Who Are Not Members of the ISMRM EMTP Study Group
In the published paper, this section will include a list of all non-members who endorsed the manuscript after the study group endorsement.

S1.5 Suggested Protocols
Please note that the protocols below are simply suggestions (with timings rounded to the nearest millisecond) and that the parameters can, of course, be adjusted where needed, according to the principles in Section 2, to accommodate differences in scanner hardware and software.  This zip file contains all data and results that were produced by running all the scripts provided in the code directory.

1.5T Protocol
The following Data Preparation section provides information on all the pre-processing steps to prepare the unmodified DICOM images to the BIDS format data that is ready for QSM reconstruction in SEPIA. If the readers are interested in the QSM reconstruction and work on the "QSM_Consensus_Paper_example_Data_Result_Code.zip" file only, they may skip Section "Data Preparation".

Data Preparation
The scripts created for this section were tested on a Mac system (macOS 13.2) and a Linux system (CentOS 7). They were not tested on Windows systems and would require adaptations to work there.
Most imaging software in the field typically deals with images in Analyze or NIfTI format.
As such the raw data (imaging data after coil combination) provided has to be converted to this format in 4 steps: Step 1: Unzip the received data and reformat the directory structure Script: Preparation_01_rename_received_data.sh Step 2: Convert DICOM images into NIfTI format Dependency: dcm2niix (version 1.0.20220720)

Script: Preparation_02_convert_dicom2nii.sh
Step 3: Rename the files according to the BIDS format (Brain Imaging Data structure)

Dependency: Matlab R2016b onwards
The naming strategy is as follows: • Vendors are identified using the session tag: ses-<GE|PHILIPS|SIEMENS> • For GE and SIEMENS, different readout methods are identified using the acquisition tag: acq-<Bipolar|Monopolar>; • For PHILIPS, the normalisation method is also printed on the acquisition tag, i.e., This section describes all the QSM reconstruction settings that were used on the example data. All the methods and algorithm parameters mentioned were already specified in the SEPIA pipeline configuration files (sepia_<GE|PHILIPS|SIEMENS>_<Monopolar|Bipolar>_config.m), which can be found in the sub-directories of the code folder "QSM_Consensus_Paper_Example_Code/": "SEPIA_Pipeline_FANSI/" and "SEPIA_Pipeline_MEDI/". Here, we provide an overview of the main parameters of each of these pipelines (Tables S2.1-S2.4) for the readers' convenience.

Processing steps
Step 1: Preparation • (GE only) Phase data must be inverted before QSM reconstruction processing (i.e., phase = -phase), so that paramagnetic susceptibility gives a positive value while diamagnetic susceptibility gives a negative value, same as the data from other vendors. This step was performed with the option provided by SEPIA.
• Brain mask is obtained by using MEDI toolbox implementation of FSL's BET on the 1st echo magnitude image, using default setting -f 0.5 -g 0 • (Bipolar readout data only) Bipolar readout correction based on  using the implementation provided with SEPIA.
• Note that the relevant sequence parameters such as echo time and slice orientation are automatically derived from the data.
Step 2: Total field estimation and echo combination Step 3: Background field removal Step 4: Dipole inversion We demonstrate the dipole inversion steps with two recommended methods (FANSI and MEDI).
Step 4.1: FANSI dipole inversion This can be done by updating the "input" variable in the configuration file to the location of the input directory that contains all the essential data in your computer. Alternatively, if a graphical operation is preferred, the SEPIA pipeline configuration files can be imported to the SEPIA's GUI by using the "Load config" button on the bottom left of the GUI display and then select the configuration .m file. The GUI will then be updated to the specified methods and algorithm parameters according to the text in the configuration file. Readers can then specify the required input and output information on the "I/O" panel on the GUI. Figure S2.2: Susceptibility maps derived using the "SEPIA_Pipeline_FANSI" processing pipeline. Figure S2.3: Susceptibility maps derived using the "SEPIA_Pipeline_MEDI" processing pipeline.  (https://github.com/korbinian90/ROMEO), which will unwrap and also combine data over coils if there is a 5th dimension (x,y,z,echo,coil).

Example results
Features: Allows acceleration in the second phase-encode direction, not subject to the echo time constraints of ASPIRE.
Limitations: Requires export of separate channel data (a large number of files).
"Virtual Reference Coil": A coil combination method suitable for systems without a body coil, e.g. UHF (Parker et al. 2014).
Availability: All systems. Limitations: Requires export of separate channel data (a large number of files). Can fail in the cerebellum, for large objects or at field strengths above 7T. Needs Advanced User privileges for installation (online version).
"Multi-echo Coil Combination": A coil combination method for multi-echo data, suitable for 3T and 7T systems.
How to: In the C2P GRE sequence (customer/gre), set Coil Combination to Adaptive Combine.
Features: The sequence produces suitable magnitude and phase DICOM data directly on the scanner. Before channel combination, the phase of the first echo is subtracted from the phase of all echoes after which the channel phases are averaged to obtain the combined phase (Eq 13 in Bernstein et al., 1994). The magnitude is obtained using sum-of-squares.
Limitations: 1) The method only works for 2 or more echoes.
2) The image reconstruction method is memory intensive. Needs Advanced User privileges for installation.
Philips "SENSE or CS-SENSE": The product 3D FFE sequence allows reconstruction of coilcombined magnitude/phase/real/imaginary data using SENSE or CS-SENSE (compressed-sensing) Availability: SENSE is available in all systems with V4 and later software versions, CS-SENSE is available in some V5 systems How to: in the "Postproc" tab -> "Images" -> Select output of "M" and "P" for magnitude and phase output. For conversion with DCM2NII and DCM2NIIX, the "Philips precise scaling" parameter should be set to ON to avoid using the other rescaling factors provided (which only adjust relative pixel intensity but do not provide quantitative rescaled values).
When using the SWIp product sequence (e.g. the clinic, to get SWI images), save the magnitude and phase data using the "Delayed reconstruction" procedure (need to turn "Postproc" tab -> "Save raw data" to "yes"). On Release 5 of the software (R5), this does not require a research key and is performed by right-clicking on the exam card and selecting the delayed recon option.
Limitations: For V5, need to set "Postproc" tab -> "Images" -> "SWIp" to "no" to allow unfiltered phase output for QSM, otherwise have to use the "Delayed reconstruction". How to: The sequence allows approximating the recommended protocol (precise TE/TR will be scanner dependent). ASSET is required to obtain correct phase, as is the disabling of any image filter and 3D geometry correction. The default 2D gradient correction works fine. Use TE = minFull. Phase Image = OFF. Must appropriately set CVs rhfiesta, rhrcctrl and rhrcxres (to the acquired matrix size).
Limitations: precise TE/TR will be scanner dependent. Export of separate channel phase and magnitude images requires a research key.
"Product sequence": The product SWAN and MERGE sequences allow reconstructing mag/real/imag by following the steps indicated below. MERGE allows shorter TEs and TRs. To modify the necessary CVs, a research key is required.
Availability: SWAN and MERGE are commercial sequences provided by the vendor.
How to: The sequence allows approximating the recommended protocol (precise TE/TR will be scanner dependent). ASSET (or ARC) is required to obtain correct phase, and TE must be set to "minFull". It is mandatory to set Phase Image = OFF in the GUI. The following CVs must be set as follows (research key required): rhfiesta=0 to keep echoes separated; rhrcctrl=13 to produce also the real and imaginary parts.
For SWAN only, on DV25 onwards, the first echo can be imposed from the "Advanced" tab in the GUI and can be set as an integer value (CV16).
Limitations: No direct control on all TEs and TR.
Gradient nonlinearities on older systems (e.g. Signa HDx) lead to distortions when images are reconstructed with FFT from k-space data. These need to be corrected with the spherical harmonics of the gradient system (which are stored on the hard drive for GE).