Methods and implementation
==========================
*SDCFlows* defines a clear :abbr:`API (application programming interface)` that divides
the problem of susceptibility distortions (SD) into two stages:
#. **Estimation**:
the MRI acquisitions in the protocol for :abbr:`SD (susceptibility distortions)` are
discovered and preprocessed to estimate a map of B\ :sub:`0` non-uniformity in Hz (:math:`\Delta B_0`).
The theory behind these distortions is well described in the literature ([Jezzard1995]_, [Hutton2002]_),
and further discussed below (see a summary in :numref:`fig-1`).
*SDCFlows* builds on freely-available software to implement three major strategies for estimating
:math:`\Delta B_0` (Eq. :math:`\eqref{eq:fieldmap-1}`).
These strategies are described below, and implemented within :py:mod:`sdcflows.workflows.fit`\ .
#. **Application**:
the B-Spline basis coefficients used to represent the estimated :math:`\Delta B_0` map mapped into the
target :abbr:`EPI (echo-planar imaging)` 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 :abbr:`PE (phase-encoding)` will be proportional to :math:`\Delta B_0 \cdot T_\text{ro}`,
where :math:`T_\text{ro}` is the *total readout time* of the target :abbr:`EPI (echo-planar imaging)` (:numref:`fig-1`).
The implementation of these workflows is found in the submodule :py:mod:`sdcflows.workflows.apply`\ .
.. _fig-1:
.. figure:: _static/sdcflows-OHBM21-fig1.png
:width: 100%
:align: center
Susceptibility distortions in a nutshell
.. admonition:: BIDS Specification
See the section `Echo-planar imaging and *B0* mapping
`__.
Fieldmap estimation: theory and methods
---------------------------------------
The displacement suffered by every voxel along the :abbr:`PE (phase-encoding)` direction
can be derived from Eq. (2) of [Hutton2002]_:
.. math::
\Delta_\text{PE} (i, j, k) = \gamma \cdot \Delta B_0 (i, j, k) \cdot T_\text{ro},
\label{eq:fieldmap-1}\tag{1}
where
:math:`\Delta_\text{PE} (i, j, k)` is the *voxel-shift map* (VSM) along the :abbr:`PE (phase-encoding)` direction,
:math:`\gamma` is the gyromagnetic ratio of the H proton in Hz/T
(:math:`\gamma = 42.576 \cdot 10^6 \, \text{Hz} \cdot \text{T}^\text{-1}`),
:math:`\Delta B_0 (i, j, k)` is the *fieldmap variation* in T (Tesla), and
:math:`T_\text{ro}` is the readout time of one slice of the :abbr:`EPI (echo-planar imaging)` dataset
we want to correct for distortions.
Let :math:`V` represent the *fieldmap* in Hz (or equivalently,
*voxel-shift-velocity map*, as Hz are equivalent to voxels/s), with
:math:`V(i,j,k) = \gamma \cdot \Delta B_0 (i, j, k)`, then, introducing
the voxel zoom along the phase-encoding direction, :math:`s_\text{PE}`,
we obtain the nonzero component of the associated displacements field
:math:`\Delta D_\text{PE} (i, j, k)` that unwarps the target :abbr:`EPI (echo-planar imaging)` dataset:
.. math::
\Delta D_\text{PE} (i, j, k) = V(i, j, k) \cdot T_\text{ro} \cdot s_\text{PE}.
\label{eq:fieldmap-2}\tag{2}
.. _sdc_direct_b0 :
Direct B0 mapping sequences
~~~~~~~~~~~~~~~~~~~~~~~~~~~
.. admonition:: BIDS Specification
See the section `Types of fieldmaps - Case 3: Direct field mapping
`__
in the BIDS specification.
Some MR schemes such as :abbr:`SEI (spiral-echo imaging)` can directly
reconstruct an estimate of *the fieldmap in Hz*, :math:`V(i,j,k)`.
These *fieldmaps* are described with more detail `here
`__.
.. _sdc_phasediff :
Phase-difference B0 estimation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
.. admonition:: 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, :math:`\Delta B_0 (i, j, k)`, that is necessary to obtain
:math:`\Delta_\text{PE} (i, j, k)` in Eq. :math:`\eqref{eq:fieldmap-1}` can be
calculated from two subsequent :abbr:`GRE (Gradient-Recalled Echo)` echoes,
via Eq. (1) of [Hutton2002]_:
.. math::
\Delta B_0 (i, j, k) = \frac{\Delta \Theta (i, j, k)}{2\pi \cdot \gamma \, \Delta\text{TE}},
\label{eq:fieldmap-3}\tag{3}
where
:math:`\Delta \Theta (i, j, k)` is the phase-difference map in radians,
and :math:`\Delta\text{TE}` is the elapsed time between the two GRE echoes.
For simplicity, the «*voxel-shift-velocity map*» :math:`V(i,j,k)`, which we
can introduce in Eq. :math:`\eqref{eq:fieldmap-2}` to directly obtain
the displacements field, can be obtained as:
.. math::
V(i, j, k) = \frac{\Delta \Theta (i, j, k)}{2\pi \cdot \Delta\text{TE}}.
\label{eq:fieldmap-4}\tag{4}
This calculation is further complicated by the fact that :math:`\Theta_i`
(and therefore, :math:`\Delta \Theta`) are clipped (or *wrapped*) within
the range :math:`[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]_).
.. _sdc_pepolar :
:abbr:`PEPOLAR (Phase Encoding POLARity)` techniques
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
.. admonition:: 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 :abbr:`PE (phase-encoding)` polarity is reversed.
*SDCFlows* integrates two implementations based on FSL ``topup`` [Andersson2003]_,
and AFNI ``3dQwarp`` [Cox1997]_.
By default, FSL ``topup`` will be used.
.. _sdc_fieldmapless :
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 :math:`B_0`
inhomogeneity map, researchers resort to nonlinear registration to an
«*anatomically correct*» map of the same individual (normally acquired with
:abbr:`T1w (T1-weighted)`, or :abbr:`T2w (T2-weighted)` 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 :abbr:`T1w (T1-weighted)` image and
a reference :abbr:`EPI (Echo-Planar Imaging)` 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 :math:`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 :abbr:`EPI (echo-planar imaging)` sources) are prepared for
registration, including a linear pre-alignment of both, calculation of a 3D :abbr:`EPI (echo-planar imaging)` *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
:py:func:`~sdcflows.workflows.fit.syn.init_syn_preprocessing_wf`.
Second, the outputs of the preprocessing workflow are fed into
:py:func:`~sdcflows.workflows.fit.syn.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 :abbr:`EPI (echo-planar imaging)`
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 :abbr:`PE (phase-encoding)` 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
:abbr:`EPI (echo-planar imaging)` given as input into its *unwarped* mapping.
The estimated transform is then converted into a :math:`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 :abbr:`EPI (echo-planar imaging)` data [Esteban2016]_.
Estimation tooling
~~~~~~~~~~~~~~~~~~
The workflows provided by :py:mod:`sdcflows.workflows.fit` make use of several utilities.
The cornerstone of these tools is the fieldmap representation with B-Splines
(:py:mod:`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
----------------------------
:py:mod:`sdcflows.workflows.apply` contains workflows to project fieldmaps represented by B-Spline
basis into the space of the target :abbr:`EPI (echo-planar imaging)` data.
Discovering fieldmaps in a BIDS dataset
---------------------------------------
To ease the implementation of higher-level pipelines integrating :abbr:`SDC (susceptibility distortion correction)`,
*SDCFlows* provides :py:func:`sdcflows.utils.wrangler.find_estimators`.
Explicit specification with ``B0FieldIdentifier``
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
.. admonition:: 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``
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
.. admonition:: 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 :abbr:`EPI (echo-planar imaging)` 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 :py:func:`~sdcflows.utils.wrangler.find_estimators`
will attempt to find :abbr:`BOLD (blood-oxygen level-dependent)` or :abbr:`DWI (diffusion-weighted imaging)`
instances within single sessions that are consistent in :abbr:`PE (phase-encoding)` 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
:abbr:`BOLD (blood-oxygen level-dependent)` or :abbr:`DWI (diffusion-weighted imaging)` 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
:py:func:`~sdcflows.utils.wrangler.find_estimators` utility.
References
----------
.. [Jezzard1995] 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 `__.
.. [Hutton2002] Hutton et al., (2002) Image Distortion Correction in fMRI: A Quantitative
Evaluation, NeuroImage 16(1):217-240. doi:`10.1006/nimg.2001.1054
`__.
.. [Jenkinson2003] Jenkinson, M. (2003) Fast, automated, N-dimensional phase-unwrapping
algorithm. MRM 49(1):193-197. doi:`10.1002/mrm.10354
`__.
.. [Andersson2003] 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 `__.
.. [Cox1997] 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
`__.
.. [Studholme2000] 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
`__.
.. [Treiber2016] 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
`_.
.. [Wang2017] 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
`_.
.. [Huntenburg2014] Huntenburg, J. M. (2014) `Evaluating Nonlinear
Coregistration of BOLD EPI and T1w Images
`__,
Berlin: Master Thesis, Freie Universität.
.. [Zaitsev2004] 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 `__.
.. [Esteban2016] 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 `__.