Skip to content

Commit

Permalink
fix: revise total-readout-time, cleanup ees
Browse files Browse the repository at this point in the history 
This commit makes a deep revisiting of our total readout time implementation.
Most of the relevant details about this commit are included within the
documentation.

With thanks to @neurolabusc because this wouldn't have been possible without
his thorough work, to @parekhpravesh for the thread he opened at
Neurostars, to @Gilles86 for making us aware of these issues on the *fMRIPrep*
repo, and to @treanus for reporting #5.

Resolves: #5.
  • Loading branch information
@oesteban
oesteban committed Dec 1, 2020
1 parent 4f2d247 commit 4ecdf6e
Show file tree
Hide file tree
Showing 4 changed files with 223 additions and 192 deletions.
27 changes: 17 additions & 10 deletions sdcflows/fieldmaps.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -228,11 +228,12 @@ def __attrs_post_init__(self):
raise MetadataError(
f"Missing 'PhaseEncodingDirection' for <{self.path}>."
)
if not (
set(("TotalReadoutTime", "EffectiveEchoSpacing")).intersection(
self.metadata.keys()
)
):

from .utils.epimanip import get_trt

try:
get_trt(self.metadata, in_file=self.path)
except ValueError:
raise MetadataError(
f"Missing readout timing information for <{self.path}>."
)
Expand Down Expand Up @@ -380,10 +381,13 @@ def __attrs_post_init__(self):
raise ValueError("Insufficient sources to estimate a fieldmap.")

if not self.bids_id:
bids_ids = set([
f.metadata.get("B0FieldIdentifier")
for f in self.sources if f.metadata.get("B0FieldIdentifier")
])
bids_ids = set(
[
f.metadata.get("B0FieldIdentifier")
for f in self.sources
if f.metadata.get("B0FieldIdentifier")
]
)
if len(bids_ids) > 1:
raise ValueError(
f"Multiple ``B0FieldIdentifier`` set: <{', '.join(bids_ids)}>"
Expand Down Expand Up @@ -414,14 +418,17 @@ def get_workflow(self, **kwargs):
str(f.path) for f in self.sources if f.suffix.startswith("magnitude")
]
self._wf.inputs.inputnode.fieldmap = [
(str(f.path), f.metadata) for f in self.sources
(str(f.path), f.metadata)
for f in self.sources
if f.suffix in ("fieldmap", "phasediff", "phase2", "phase1")
]
elif self.method == EstimatorType.PEPOLAR:
from .workflows.fit.pepolar import init_topup_wf

self._wf = init_topup_wf(**kwargs)
elif self.method == EstimatorType.ANAT:
from .workflows.fit.syn import init_syn_sdc_wf

self._wf = init_syn_sdc_wf(**kwargs)

return self._wf
185 changes: 3 additions & 182 deletions sdcflows/interfaces/epi.py
View file
Original file line number Diff line number Diff line change
@@ -1,17 +1,4 @@
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""
Interfaces to deal with the various types of fieldmap sources.
.. testsetup::
>>> tmpdir = getfixture('tmpdir')
>>> tmp = tmpdir.chdir() # changing to a temporary directory
>>> nb.Nifti1Image(np.zeros((90, 90, 60)), None, None).to_filename(
... tmpdir.join('epi.nii.gz').strpath)
"""
"""Interfaces to deal with the various types of fieldmap sources."""

from nipype.interfaces.base import (
BaseInterfaceInputSpec,
Expand Down Expand Up @@ -40,177 +27,11 @@ class GetReadoutTime(SimpleInterface):
output_spec = _GetReadoutTimeOutputSpec

def _run_interface(self, runtime):
from ..utils.epimanip import get_trt

self._results["readout_time"] = get_trt(
self.inputs.metadata,
self.inputs.in_file if isdefined(self.inputs.in_file) else None,
)
self._results["pe_direction"] = self.inputs.metadata["PhaseEncodingDirection"]
return runtime


def get_trt(in_meta, in_file=None):
r"""
Extract the *total readout time* :math:`t_\text{RO}` from BIDS.
Calculate the *total readout time* for an input
:abbr:`EPI (echo-planar imaging)` scan.
There are several procedures to calculate the total
readout time. The basic one is that a ``TotalReadoutTime``
field is set in the JSON sidecar. The following examples
use an ``'epi.nii.gz'`` file-stub which has 90 pixels in the
j-axis encoding direction.
>>> meta = {'TotalReadoutTime': 0.05251}
>>> get_trt(meta)
0.05251
If the *effective echo spacing* :math:`t_\text{ees}`
(``EffectiveEchoSpacing`` BIDS field) is provided, then the
total readout time can be calculated reading the number
of voxels along the readout direction :math:`T_\text{ro}`
and the parallel acceleration factor of the EPI :math:`f_\text{acc}`.
.. math ::
T_\text{ro} = t_\text{ees} \, (N_\text{PE} / f_\text{acc} - 1)
>>> meta = {'EffectiveEchoSpacing': 0.00059,
... 'PhaseEncodingDirection': 'j-',
... 'ParallelReductionFactorInPlane': 2}
>>> get_trt(meta, in_file='epi.nii.gz')
0.05251
Some vendors, like Philips, store different parameter names:
>>> meta = {'WaterFatShift': 8.129,
... 'MagneticFieldStrength': 3,
... 'PhaseEncodingDirection': 'j-',
... 'ParallelReductionFactorInPlane': 2}
>>> get_trt(meta, in_file='epi.nii.gz')
0.018721183563864822
"""
import nibabel as nb

# Use case 1: TRT is defined
trt = in_meta.get("TotalReadoutTime", None)
if trt is not None:
return trt

# npe = N voxels PE direction
npe = nb.load(in_file).shape[_get_pe_index(in_meta)]

# Use case 2: TRT is defined
ees = in_meta.get("EffectiveEchoSpacing", None)
if ees is not None:
# Effective echo spacing means that acceleration factors have been accounted for.
return ees * (npe - 1)

# All other cases require the parallel acc and npe (N vox in PE dir)
acc = float(in_meta.get("ParallelReductionFactorInPlane", 1.0))
etl = npe // acc # effective train length
es = in_meta.get("EchoSpacing", None)
if es is not None:
return es * (etl - 1)

# Use case 3 (philips scans)
wfs = in_meta.get("WaterFatShift", None)
if wfs is not None:
fstrength = in_meta["MagneticFieldStrength"]
wfd_ppm = 3.4 # water-fat diff in ppm
g_ratio_mhz_t = 42.57 # gyromagnetic ratio for proton (1H) in MHz/T
wfs_hz = fstrength * wfd_ppm * g_ratio_mhz_t
return wfs / wfs_hz

raise ValueError("Unknown total-readout time specification")


def get_ees(in_meta, in_file=None):
r"""
Extract the *effective echo spacing* :math:`t_\text{ees}` from BIDS.
Calculate the *effective echo spacing* :math:`t_\text{ees}`
for an input :abbr:`EPI (echo-planar imaging)` scan.
There are several procedures to calculate the effective
echo spacing. The basic one is that an ``EffectiveEchoSpacing``
field is set in the JSON sidecar. The following examples
use an ``'epi.nii.gz'`` file-stub which has 90 pixels in the
j-axis encoding direction.
>>> meta = {'EffectiveEchoSpacing': 0.00059,
... 'PhaseEncodingDirection': 'j-'}
>>> get_ees(meta)
0.00059
If the *total readout time* :math:`T_\text{ro}` (``TotalReadoutTime``
BIDS field) is provided, then the effective echo spacing can be
calculated reading the number of voxels :math:`N_\text{PE}` along the
readout direction and the parallel acceleration
factor of the EPI
.. math ::
= T_\text{ro} \, (N_\text{PE} / f_\text{acc} - 1)^{-1}
where :math:`N_y` is the number of pixels along the phase-encoding direction
:math:`y`, and :math:`f_\text{acc}` is the parallel imaging acceleration factor
(:abbr:`GRAPPA (GeneRalized Autocalibrating Partial Parallel Acquisition)`,
:abbr:`ARC (Autocalibrating Reconstruction for Cartesian imaging)`, etc.).
>>> meta = {'TotalReadoutTime': 0.02596,
... 'PhaseEncodingDirection': 'j-',
... 'ParallelReductionFactorInPlane': 2}
>>> get_ees(meta, in_file='epi.nii.gz')
0.00059
Some vendors, like Philips, store different parameter names (see
http://dbic.dartmouth.edu/pipermail/mrusers/attachments/20141112/eb1d20e6/attachment.pdf
):
>>> meta = {'WaterFatShift': 8.129,
... 'MagneticFieldStrength': 3,
... 'PhaseEncodingDirection': 'j-',
... 'ParallelReductionFactorInPlane': 2}
>>> get_ees(meta, in_file='epi.nii.gz')
0.00041602630141921826
"""
import nibabel as nb
from sdcflows.interfaces.epi import _get_pe_index

# Use case 1: EES is defined
ees = in_meta.get("EffectiveEchoSpacing", None)
if ees is not None:
return ees

# All other cases require the parallel acc and npe (N vox in PE dir)
acc = float(in_meta.get("ParallelReductionFactorInPlane", 1.0))
npe = nb.load(in_file).shape[_get_pe_index(in_meta)]
etl = npe // acc

# Use case 2: TRT is defined
trt = in_meta.get("TotalReadoutTime", None)
if trt is not None:
return trt / (etl - 1)

# Use case 3 (philips scans)
wfs = in_meta.get("WaterFatShift", None)
if wfs is not None:
fstrength = in_meta["MagneticFieldStrength"]
wfd_ppm = 3.4 # water-fat diff in ppm
g_ratio_mhz_t = 42.57 # gyromagnetic ratio for proton (1H) in MHz/T
wfs_hz = fstrength * wfd_ppm * g_ratio_mhz_t
return wfs / (wfs_hz * etl)

raise ValueError("Unknown effective echo-spacing specification")


def _get_pe_index(meta):
pe = meta["PhaseEncodingDirection"]
try:
return {"i": 0, "j": 1, "k": 2}[pe[0]]
except KeyError:
raise RuntimeError('"%s" is an invalid PE string' % pe)
Loading

0 comments on commit 4ecdf6e

Please sign in to comment.