# Methods and implementation#

*SDCFlows* defines a clear API that divides
the problem of susceptibility distortions (SD) into two stages:

**Estimation**: the MRI acquisitions in the protocol for SD are discovered and preprocessed to estimate a map of B_{0}non-uniformity in Hz (\(\Delta B_0\)). The theory behind these distortions is well described in the literature ([Jezzard1995], [Hutton2002]), and further discussed below (see a summary in Fig. 1).*SDCFlows*builds on freely-available software to implement three major strategies for estimating \(\Delta B_0\) (Eq. \(\eqref{eq:fieldmap-1}\)). These strategies are described below, and implemented within`sdcflows.workflows.fit`

.**Application**: the B-Spline basis coefficients used to represent the estimated \(\Delta B_0\) map mapped into the target EPI scan to be corrected, and a displacement field in NIfTI format that is compatible with ANTs is interpolated from the B-Spline basis. The voxel location error along the PE will be proportional to \(\Delta B_0 \cdot T_\text{ro}\), where \(T_\text{ro}\) is the*total readout time*of the target EPI (Fig. 1). The implementation of these workflows is found in the submodule`sdcflows.workflows.apply`

.

BIDS Specification

See the section Echo-planar imaging and *B0* mapping.

## Fieldmap estimation: theory and methods#

The displacement suffered by every voxel along the PE direction can be derived from Eq. (2) of [Hutton2002]:

where
\(\Delta_\text{PE} (i, j, k)\) is the *voxel-shift map* (VSM) along the PE direction,
\(\gamma\) is the gyromagnetic ratio of the H proton in Hz/T
(\(\gamma = 42.576 \cdot 10^6 \, \text{Hz} \cdot \text{T}^\text{-1}\)),
\(\Delta B_0 (i, j, k)\) is the *fieldmap variation* in T (Tesla), and
\(T_\text{ro}\) is the readout time of one slice of the EPI dataset
we want to correct for distortions.

Let \(V\) represent the *fieldmap* in Hz (or equivalently,
*voxel-shift-velocity map*, as Hz are equivalent to voxels/s), with
\(V(i,j,k) = \gamma \cdot \Delta B_0 (i, j, k)\), then, introducing
the voxel zoom along the phase-encoding direction, \(s_\text{PE}\),
we obtain the nonzero component of the associated displacements field
\(\Delta D_\text{PE} (i, j, k)\) that unwarps the target EPI dataset:

### Direct B0 mapping sequences#

BIDS Specification

See the section Types of fieldmaps - Case 3: Direct field mapping in the BIDS specification.

Some MR schemes such as SEI can directly
reconstruct an estimate of *the fieldmap in Hz*, \(V(i,j,k)\).
These *fieldmaps* are described with more detail here.

### Phase-difference B0 estimation#

BIDS Specification

See the section Types of fieldmaps - Case 2: Two phase maps and two magnitude images in the BIDS specification.

Some scanners produce one `phasediff`

map, where the drift between the two echos has
already been calculated, see the section
Types of fieldmaps - Case 1: Phase-difference map and at least one magnitude image.

The fieldmap variation in T, \(\Delta B_0 (i, j, k)\), that is necessary to obtain \(\Delta_\text{PE} (i, j, k)\) in Eq. \(\eqref{eq:fieldmap-1}\) can be calculated from two subsequent GRE echoes, via Eq. (1) of [Hutton2002]:

where \(\Delta \Theta (i, j, k)\) is the phase-difference map in radians, and \(\Delta\text{TE}\) is the elapsed time between the two GRE echoes.

For simplicity, the «*voxel-shift-velocity map*» \(V(i,j,k)\), which we
can introduce in Eq. \(\eqref{eq:fieldmap-2}\) to directly obtain
the displacements field, can be obtained as:

This calculation is further complicated by the fact that \(\Theta_i\)
(and therefore, \(\Delta \Theta\)) are clipped (or *wrapped*) within
the range \([0 \dotsb 2\pi )\).
It is necessary to find the integer number of offsets that make a region
continuously smooth with its neighbors (*phase-unwrapping*, [Jenkinson2003]).

### PEPOLAR techniques#

BIDS Specification

See the section Types of fieldmaps - Case 4: Multiple phase encoded directions (“pepolar”).

Alternatively, it is possible to estimate the field by exploiting the symmetry of the
distortion when the PE polarity is reversed.
*SDCFlows* integrates two implementations based on FSL `topup`

[Andersson2003],
and AFNI `3dQwarp`

[Cox1997].
By default, FSL `topup`

will be used.

### Fieldmap-less approaches#

Many studies acquired (especially with legacy MRI protocols) do not have any
information to estimate susceptibility-derived distortions.
In the absence of data with the specific purpose of estimating the \(B_0\)
inhomogeneity map, researchers resort to nonlinear registration to an
«*anatomically correct*» map of the same individual (normally acquired with
T1w, or T2w sequences).
One of the most prominent proposals of this approach is found in [Studholme2000].

*SDCFlows* includes an (experimental) procedure, based on nonlinear image registration
with ANTs’ symmetric normalization (SyN) technique.
This workflow takes a skull-stripped T1w image and
a reference EPI image, and estimates a field of nonlinear
displacements that accounts for susceptibility-derived distortions.
To more accurately estimate the warping on typically distorted regions, this
implementation uses an average \(B_0\) mapping described in [Treiber2016].
The implementation is a variation on those developed in [Huntenburg2014] and
[Wang2017].

The process is divided in two steps.
First, the two images to be aligned (anatomical and one or more EPI sources) are prepared for
registration, including a linear pre-alignment of both, calculation of a 3D EPI *reference* map,
intensity/histogram enhancement, and *deobliquing* (meaning, for images where the physical
coordinates axes and the data array axes are not aligned, the physical coordinates are
rotated to align with the data array).
Such a preprocessing is implemented in
`init_syn_preprocessing_wf()`

.
Second, the outputs of the preprocessing workflow are fed into
`init_syn_sdc_wf()`

,
which executes the nonlinear, SyN registration.
To aid the *Mattes* mutual information cost function, the registration scheme is set up
in *multi-channel* mode, and laplacian-filtered derivatives of both anatomical and EPI
reference are introduced as a second registration channel.
The optimization gradients of the registration process are weighted, so that deformations
effectively possible only along the PE axis.
Given that ANTs’ registration framework performs on physical coordinates, it is necessary
that input images are not *oblique*.
The anatomical image is used as *fixed image*, and therefore, the registration process
estimates the transformation function from *unwarped* (anatomically *correct*) coordinates
to *distorted* coordinates.
If fed into `antsApplyTransforms`

, the resulting transform will effectively *unwarp* a distorted
EPI given as input into its *unwarped* mapping.
The estimated transform is then converted into a \(B_0\) fieldmap in Hz, which can be
stored within the derivatives folder.

Danger

This procedure is experimental, and the outcomes should be scrutinized one-by-one and used with caution. Feedback will be enthusiastically received.

### Other (unsupported) approaches#

There exist some alternative options to estimate the fieldmap, such as mapping the point-spread-function [Zaitsev2004], or by means of nonlinear registration of brain surfaces onto the distorted EPI data [Esteban2016].

### Estimation tooling#

The workflows provided by `sdcflows.workflows.fit`

make use of several utilities.
The cornerstone of these tools is the fieldmap representation with B-Splines
(`sdcflows.interfaces.bspline`

).
B-Splines are well-suited to plausibly smooth the fieldmap and provide a compact
representation of the field with fewer parameters.
This representation is also more accurate in the case the images that were used for estimation
are not aligned with the target images to be corrected because the fieldmap is not directly
interpolated in the projection, but rather, the position of the coefficients in space is
updated and then the fieldmap reconstructed.

## Unwarping the distorted data#

`sdcflows.workflows.apply`

contains workflows to project fieldmaps represented by B-Spline
basis into the space of the target EPI data.

## Discovering fieldmaps in a BIDS dataset#

To ease the implementation of higher-level pipelines integrating SDC,
*SDCFlows* provides `sdcflows.utils.wrangler.find_estimators()`

.

### Explicit specification with `B0FieldIdentifier`

#

BIDS Specification

See the section Expressing the MR protocol intent for fieldmaps - Using *B0FieldIdentifier* metadata.

If one or more `B0FieldIdentifier`

(s) are set within the input metadata (following BIDS’ specifications),
then corresponding estimators will be built from the available input data.

### Implicit specification with `IntendedFor`

#

BIDS Specification

See the section Expressing the MR protocol intent for fieldmaps - Using *IntendedFor* metadata.

In the case no `B0FieldIdentifier`

(s) are defined, then *SDCFlows* will try to identify as many fieldmap
estimators as possible within the dataset following a set of heuristics.
Then, those estimators may be linked to target EPI data by interpreting the
`IntendedFor`

field if available.

### Fieldmap-less by registration to a T1-weighted image#

If none of the two previous options yielded any workable estimation strategy, and provided that
the argument `fmapless`

is set to `True`

, then `find_estimators()`

will attempt to find BOLD or DWI
instances within single sessions that are consistent in PE direction and
*total readout time*, assuming they have been acquired with the same shimming settings.

If one or more anatomical images are found, and if the search for candidate BOLD or DWI data yields results, then corresponding fieldmap-less estimators are set up.

It is possible to force the fieldmap-less estimation by passing `force_fmapless=True`

to the
`find_estimators()`

utility.

## References#

Jezzard, P. & Balaban, R. S. (1995) Correction for geometric distortion in echo planar images from B0 field variations. Magn. Reson. Med. 34:65–73. doi:10.1002/mrm.1910340111.

Hutton et al., (2002) Image Distortion Correction in fMRI: A Quantitative Evaluation, NeuroImage 16(1):217-240. doi:10.1006/nimg.2001.1054.

Jenkinson, M. (2003) Fast, automated, N-dimensional phase-unwrapping algorithm. MRM 49(1):193-197. doi:10.1002/mrm.10354.

Andersson, J. (2003) How to correct susceptibility distortions in spin-echo echo-planar images: application to diffusion tensor imaging. NeuroImage 20:870–888. doi:10.1016/s1053-8119(03)00336-7.

Cox, R. (1997) Software tools for analysis and visualization of fMRI data. NMR Biomed. 10:171–178, doi:10.1002/(sici)1099-1492(199706/08)10:4/5%3C171::aid-nbm453%3E3.0.co;2-l.

Studholme et al. (2000) Accurate alignment of functional EPI data to anatomical MRI using a physics-based distortion model, IEEE Trans Med Imag 19(11):1115-1127, 2000, doi: 10.1109/42.896788.

Treiber, J. M. et al. (2016) Characterization and Correction of Geometric Distortions in 814 Diffusion Weighted Images, PLoS ONE 11(3): e0152472. doi:10.1371/journal.pone.0152472.

Wang S, et al. (2017) Evaluation of Field Map and Nonlinear Registration Methods for Correction of Susceptibility Artifacts in Diffusion MRI. Front. Neuroinform. 11:17. doi:10.3389/fninf.2017.00017.

Huntenburg, J. M. (2014) Evaluating Nonlinear Coregistration of BOLD EPI and T1w Images, Berlin: Master Thesis, Freie Universität.

Zaitsev, M. (2004) Point spread function mapping with parallel imaging techniques and high acceleration factors: Fast, robust, and flexible method for echo-planar imaging distortion correction, MRM 52(5):1156-1166. doi:10.1002/mrm.20261.

Esteban, O. (2016) Surface-driven registration method for the structure-informed segmentation of diffusion MR images. NeuroImage 139:450-461. doi:10.1016/j.neuroimage.2016.05.011.