Susceptibility distortion correction

qsiprep.workflows.fieldmap.base.init_sdc_wf()

Hover over the image to flicker between the distorted and the corrected b=0 of each phase encoding direction; the contours are the subject’s T2w. Both directions end up on the same image, which is the sign that the correction worked.

Echo-planar images are stretched and compressed along the phase encoding axis where the magnetic field is inhomogeneous. Correcting this needs an estimate of the field, which can come from a fieldmap, from images acquired with opposite phase encoding, or from registration to an undistorted anatomical image. QSIPrep decides which one a series gets from the metadata (see Linking fieldmaps to DWI series), prefers reverse phase encoding over a GRE fieldmap when both are linked, and uses an anatomical reference only as a fallback unless --force sdc-anat-reference is given. The tool that applies the estimate depends on the head motion backend:

Estimated from	eddy	tortoise	shoreline
Reverse phase-encoded b=0 or DWI series	TOPUP, inside eddy; optionally refined by DRBUDDI	DRBUDDI	DRBUDDI
GRE fieldmap	Inside eddy	Applied after DIFFPREP	Applied after SHORELine
Synthetic b=0 (synb0)	TOPUP, inside eddy	T2Wreg	Not available
T2w (t2w)	Not available	T2Wreg	Not available
Inverted T1w (invt1w)	SyN	SyN	SyN

The correction that ran is named in the visual report, in the methods boilerplate, and in the EstimationMethod of the displacement map derivative (see Outputs).

TOPUP

TOPUP [1] estimates the field from b=0 images acquired with different phase encoding directions or readout times. From each distortion group of the correction unit, QSIPrep selects b=0 images, preferring those native to the series being corrected over borrowed ones, and passes them to TOPUP. The field is applied by eddy inside its own resampling, so it never exists as a standalone correction; for the displacement map derivative it is rebuilt from TOPUP’s field in Hz. TOPUP is only available with --hmc-method eddy.

DRBUDDI

qsiprep.workflows.fieldmap.drbuddi.init_drbuddi_wf()

DRBUDDI [2], from TORTOISE, estimates the field by registering the two phase encoding directions to each other after head motion correction. With reverse phase-encoded DWI series it uses the b=0 image and the fractional anisotropy image of each direction together, in a multi-modal registration; with a lone reverse b=0 (an epi fieldmap) it registers the b=0 images. A T2w image is added to the registration when one is available and t2w is not ignored; it is first rotated into the b=0 frame by an antsAI rotation search (see Registration to the anatomical), because DRBUDDI’s own rigid initialization cannot recover large rotations. Intensities are then adjusted with TORTOISE’s least-squares restoration (see Intensity modulation).

DRBUDDI runs after eddy, DIFFPREP or SHORELine, and processes exactly one matched pair of phase encoding directions per invocation. On the eddy backend it can either replace TOPUP (--sdc-method drbuddi) or refine it (--sdc-method topup+drbuddi), in which case the refinement is written as a separate desc-sdcrefinement displacement map.

The report also shows the fractional anisotropy image before and after the correction (hover to flicker). The FA is what DRBUDDI matched between the two directions, so distorted white matter should straighten out in it even where the b=0 is featureless.

(Source code, png, svg, pdf)

GRE fieldmaps

qsiprep.workflows.fieldmap.phdiff.init_phdiff_wf(), qsiprep.workflows.fieldmap.fmap.init_fmap_wf()

A phase-difference or two-phase fieldmap is unwrapped with ROMEO [3] as implemented in niimath [4], and converted to a field in Hz [5]; a Hz fieldmap is masked and median-filtered. The magnitude image is registered to the b=0 reference, and the field is turned into a displacement along the phase encoding axis using the series’ TotalReadoutTime and applied with ANTs. Since 26.1 none of this needs FSL, so GRE fieldmaps work in the FSL-free image as well.

With --hmc-method eddy, the field is handed to eddy and applied inside its model. With the other backends it is applied to the motion-corrected series at the final resampling, and an unwarped b=0 reference is computed for coregistration.

(Source code, png, svg, pdf)

Initializing DRBUDDI and T2Wreg from a GRE fieldmap

Registration-based correction infers the field by matching images, which is ambiguous where the field piles several voxels’ signal into one or drops it out. A GRE fieldmap measures the field directly. So when a series is corrected by DRBUDDI or T2Wreg and a GRE fieldmap also lists it, the GRE-derived warp is the registration’s starting point: for DRBUDDI, the initial transform of the primary direction and, negated, of the reverse direction; for T2Wreg, the initial transform. The initial warp is held fixed through the registration’s multi-resolution pyramid, so each stage estimates a residual correction on top of it rather than smoothing it away.

This happens whenever the linkage is present, on every backend that runs DRBUDDI (except when DRBUDDI only refines TOPUP) and on T2Wreg when the anatomical reference is forced. How to set up the sidecars is in Starting DRBUDDI or T2Wreg from a GRE fieldmap. The report’s distortion correction entry ends in (GRE-initialized) and the boilerplate describes the initialization.

SynB0

qsiprep.workflows.fieldmap.synb0.init_synb0_wf()

With --sdc-anat-reference synb0, a synthetic distortion-free b=0 image is generated from the T1w and the distorted b=0 with the SynB0-DISCO U-Net [6]. The synthetic image then stands in for a reverse phase-encoded acquisition: with eddy it joins the TOPUP inputs as a volume with zero readout time, and with DIFFPREP it is the T2Wreg target. The report shows the acquired and synthetic b=0 as a flicker pair, and the correction that TOPUP or T2Wreg estimated from them. Requires a T1w and a PhaseEncodingDirection on the DWI series.

The acquired b=0 and the synthetic one the U-Net produced from it and the T1w (hover to flicker). Nothing has been corrected yet; this view is for judging whether the synthetic target is credible. The contours are the white matter of the T1w.

The correction TOPUP estimated from that pair, before and after. This is the same session as the DRBUDDI example above, but with only one phase encoding direction and SynB0 as the reference.

The displacement map TOPUP estimated from the acquired and synthetic b=0. Compare it with the DRBUDDI map of the same session in Outputs: the field is recovered without a reverse phase-encoded acquisition, with the largest displacements in the same places.

T2Wreg

TORTOISE’s T2Wreg (--epi T2Wreg in DIFFPREP) nonlinearly registers the b=0 to the subject’s T2w image, or to the synthetic b=0 above, in the same run as head motion correction. The T2w is rotated into the b=0 frame by an antsAI rotation search first, for the same reason as with DRBUDDI. It is only available with --hmc-method tortoise. Its correction is applied without Jacobian modulation, as in TORTOISE; --force jacobian modulates it anyway.

SyN

qsiprep.workflows.fieldmap.syn.init_syn_sdc_wf()

With --sdc-anat-reference invt1w, the fieldmap-less method from fMRIPrep [7][8] registers the b=0 to the T1w with its intensity inverted, with the deformation constrained to the phase encoding axis and regularized by an average fieldmap template [9]. It runs after any head motion backend and is never selected by auto.

Intensity modulation

Correcting a spatial distortion moves signal between voxels, so the corrected image has to be rescaled by the local volume change, or regions the acquisition compressed stay too bright [10]. QSIPrep follows TORTOISE’s two implementations. Jacobian: each volume is multiplied by the derivative of the composed displacement along the phase encoding axis, used when only one phase encoding polarity was acquired (GRE and SyN fieldmaps, and single-polarity DIFFPREP runs). LSR: when both polarities exist, DRBUDDI computes the harmonic mean of the two corrected b=0 images and each direction’s volumes are scaled by the ratio of that reference to their own corrected b=0; this replaces the Jacobian entirely. Head motion, coregistration and normalization are never modulated: a rigid realignment is not a volume change, and modulating by a normalization warp would corrupt the signal.

Correction	Modulated by
Gradient nonlinearity	QSIPrep (Jacobian; not under LSR)
Susceptibility (TOPUP)	eddy, internally
Susceptibility (DRBUDDI)	QSIPrep (LSR)
Susceptibility (GRE)	eddy on its path, else QSIPrep
Susceptibility (SyN)	QSIPrep (Jacobian)
Susceptibility (T2Wreg)	none, unless --force jacobian
Eddy current (eddy)	eddy, internally
Eddy current (DIFFPREP)	QSIPrep (Jacobian; not under LSR)

--ignore jacobian disables only the modulation QSIPrep itself applies. eddy modulates whenever its resampling method is jac (the default); with "method": "lsr" in --eddy-config it does not modulate at all, and QSIPrep cannot retrofit that. The run then warns, and the gap is recorded in the UnmodulatedCorrections key of the Jacobian derivative’s sidecar. The weights QSIPrep applied are written as the desc-jacobian derivative described in Outputs.

References

[1]

Jesper L.R. Andersson, Stefan Skare, and John Ashburner. How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging. NeuroImage, 20(2):870–888, 2003. doi:10.1016/S1053-8119(03)00336-7.

[2]

M. Okan Irfanoglu, Pooja Modi, Amritha Nayak, Elizabeth B. Hutchinson, Joelle Sarlls, and Carlo Pierpaoli. Dr-buddi (diffeomorphic registration for blip-up blip-down diffusion imaging) method for correcting echo planar imaging distortions. NeuroImage, 106:284–299, 2015. doi:10.1016/j.neuroimage.2014.11.042.

[3]

Barbara Dymerska, Korbinian Eckstein, Beata Bachrata, Bernard Siow, Siegfried Trattnig, Karin Shmueli, and Simon Daniel Robinson. Phase unwrapping with a rapid opensource minimum spanning tree algorithm (ROMEO). Magnetic Resonance in Medicine, 85(4):2294–2308, 2021. doi:10.1002/mrm.28563.

[4]

Christopher Rorden, Matthew Webster, Chris Drake, Mark Jenkinson, Jonathan D. Clayden, Ningfei Li, and Taylor Hanayik. Niimath and fslmaths: replication as a method to enhance popular neuroimaging tools. Aperture Neuro, 2024. doi:10.52294/001c.94384.

[5]

Peter Jezzard and Robert S. Balaban. Correction for geometric distortion in echo planar images from B0 field variations. Magnetic Resonance in Medicine, 34(1):65–73, 1995. doi:10.1002/mrm.1910340111.

[6]

Kurt G. Schilling, Justin Blaber, Yuankai Huo, Allen Newton, Colin Hansen, Vishwesh Nath, Andrea T. Shafer, Owen Williams, Susan M. Resnick, Baxter Rogers, Adam W. Anderson, and Bennett A. Landman. Synthesized b0 for diffusion distortion correction (Synb0-DisCo). Magnetic Resonance Imaging, 64:62–70, 2019. doi:10.1016/j.mri.2019.05.008.

[7]

Sijia Wang, Daniel J. Peterson, J. C. Gatenby, Wenbin Li, Thomas J. Grabowski, and Tara M. Madhyastha. Evaluation of field map and nonlinear registration methods for correction of susceptibility artifacts in diffusion mri. Frontiers in Neuroinformatics, 2017. doi:10.3389/fninf.2017.00017.

[8]

Julia M. Huntenburg. Evaluating nonlinear coregistration of bold epi and t1w images. Master’s thesis, Freie Universität, Berlin, 2014.

[9]

Jeffrey Mark Treiber, Nathan S. White, Tyler Christian Steed, Hauke Bartsch, Dominic Holland, Nikdokht Farid, Carrie R. McDonald, Bob S. Carter, Anders Martin Dale, and Clark C. Chen. Characterization and correction of geometric distortions in 814 diffusion weighted images. PLOS ONE, 11(3):e0152472, 2016. doi:10.1371/journal.pone.0152472.

[10]

Guenther K. Rohde, Alan S. Barnett, Peter J. Basser, Stefano Marenco, and Carlo Pierpaoli. Comprehensive approach for correction of motion and distortion in diffusion-weighted mri. Magnetic Resonance in Medicine, 51(1):103–114, 2004. doi:10.1002/mrm.10677.