#!/usr/bin/env python
# -*- coding: utf-8 -*-
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""
Phase-difference B0 estimation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The field inhomogeneity inside the scanner (fieldmap) is proportional to the
phase drift between two subsequent :abbr:`GRE (gradient recall echo)`
sequence.
Fieldmap preprocessing workflow for fieldmap data structure
8.9.1 in BIDS 1.0.0: one phase diff and at least one magnitude image
"""
from niworkflows.nipype.interfaces import ants, fsl, utility as niu
from niworkflows.nipype.pipeline import engine as pe
# Note that demean_image imports from nipype
from niworkflows.nipype.workflows.dmri.fsl.utils import siemens2rads, demean_image, \
cleanup_edge_pipeline
from niworkflows.interfaces.masks import BETRPT
from ...interfaces import ReadSidecarJSON, IntraModalMerge, DerivativesDataSink
[docs]def init_phdiff_wf(reportlets_dir, omp_nthreads, name='phdiff_wf'):
"""
Estimates the fieldmap using a phase-difference image and one or more
magnitude images corresponding to two or more :abbr:`GRE (Gradient Echo sequence)`
acquisitions. The `original code was taken from nipype
<https://github.com/nipy/nipype/blob/master/nipype/workflows/dmri/fsl/artifacts.py#L514>`_.
.. workflow ::
:graph2use: orig
:simple_form: yes
from fmriprep.workflows.fieldmap.phdiff import init_phdiff_wf
wf = init_phdiff_wf(reportlets_dir='.', omp_nthreads=1)
Outputs::
outputnode.fmap_ref - The average magnitude image, skull-stripped
outputnode.fmap_mask - The brain mask applied to the fieldmap
outputnode.fmap - The estimated fieldmap in Hz
"""
inputnode = pe.Node(niu.IdentityInterface(fields=['magnitude', 'phasediff']),
name='inputnode')
outputnode = pe.Node(niu.IdentityInterface(
fields=['fmap', 'fmap_ref', 'fmap_mask']), name='outputnode')
def _pick1st(inlist):
return inlist[0]
# Read phasediff echo times
meta = pe.Node(ReadSidecarJSON(), name='meta', mem_gb=0.01, run_without_submitting=True)
dte = pe.Node(niu.Function(function=_delta_te), name='dte', mem_gb=0.01)
# Merge input magnitude images
magmrg = pe.Node(IntraModalMerge(), name='magmrg')
# de-gradient the fields ("bias/illumination artifact")
n4 = pe.Node(ants.N4BiasFieldCorrection(
dimension=3, copy_header=True, num_threads=omp_nthreads), name='n4', n_procs=omp_nthreads)
bet = pe.Node(BETRPT(generate_report=True, frac=0.6, mask=True),
name='bet')
ds_fmap_mask = pe.Node(DerivativesDataSink(
base_directory=reportlets_dir, suffix='fmap_mask'), name='ds_fmap_mask',
mem_gb=0.01, run_without_submitting=True)
# uses mask from bet; outputs a mask
# dilate = pe.Node(fsl.maths.MathsCommand(
# nan2zeros=True, args='-kernel sphere 5 -dilM'), name='MskDilate')
# phase diff -> radians
pha2rads = pe.Node(niu.Function(function=siemens2rads), name='pha2rads')
# FSL PRELUDE will perform phase-unwrapping
prelude = pe.Node(fsl.PRELUDE(), name='prelude')
denoise = pe.Node(fsl.SpatialFilter(operation='median', kernel_shape='sphere',
kernel_size=3), name='denoise')
demean = pe.Node(niu.Function(function=demean_image), name='demean')
cleanup_wf = cleanup_edge_pipeline(name="cleanup_wf")
compfmap = pe.Node(niu.Function(function=phdiff2fmap), name='compfmap')
# The phdiff2fmap interface is equivalent to:
# rad2rsec (using rads2radsec from nipype.workflows.dmri.fsl.utils)
# pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True), name='ComputeFieldmapFUGUE')
# rsec2hz (divide by 2pi)
workflow = pe.Workflow(name=name)
workflow.connect([
(inputnode, meta, [('phasediff', 'in_file')]),
(inputnode, magmrg, [('magnitude', 'in_files')]),
(magmrg, n4, [('out_avg', 'input_image')]),
(n4, prelude, [('output_image', 'magnitude_file')]),
(n4, bet, [('output_image', 'in_file')]),
(bet, prelude, [('mask_file', 'mask_file')]),
(inputnode, pha2rads, [('phasediff', 'in_file')]),
(pha2rads, prelude, [('out', 'phase_file')]),
(meta, dte, [('out_dict', 'in_values')]),
(dte, compfmap, [('out', 'delta_te')]),
(prelude, denoise, [('unwrapped_phase_file', 'in_file')]),
(denoise, demean, [('out_file', 'in_file')]),
(demean, cleanup_wf, [('out', 'inputnode.in_file')]),
(bet, cleanup_wf, [('mask_file', 'inputnode.in_mask')]),
(cleanup_wf, compfmap, [('outputnode.out_file', 'in_file')]),
(compfmap, outputnode, [('out', 'fmap')]),
(bet, outputnode, [('mask_file', 'fmap_mask'),
('out_file', 'fmap_ref')]),
(inputnode, ds_fmap_mask, [('phasediff', 'source_file')]),
(bet, ds_fmap_mask, [('out_report', 'in_file')]),
])
return workflow
# ------------------------------------------------------
# Helper functions
# ------------------------------------------------------
[docs]def phdiff2fmap(in_file, delta_te, out_file=None):
r"""
Converts the input phase-difference map into a fieldmap in Hz,
using the eq. (1) of [Hutton2002]_:
.. math::
\Delta B_0 (\text{T}^{-1}) = \frac{\Delta \Theta}{2\pi\gamma \Delta\text{TE}}
In this case, we do not take into account the gyromagnetic ratio of the
proton (:math:`\gamma`), since it will be applied inside TOPUP:
.. math::
\Delta B_0 (\text{Hz}) = \frac{\Delta \Theta}{2\pi \Delta\text{TE}}
.. [Hutton2002] Hutton et al., Image Distortion Correction in fMRI: A Quantitative
Evaluation, NeuroImage 16(1):217-240, 2002. doi:`10.1006/nimg.2001.1054
<http://dx.doi.org/10.1006/nimg.2001.1054>`_.
"""
import numpy as np
import nibabel as nb
import os.path as op
import math
# GYROMAG_RATIO_H_PROTON_MHZ = 42.576
if out_file is None:
fname, fext = op.splitext(op.basename(in_file))
if fext == '.gz':
fname, _ = op.splitext(fname)
out_file = op.abspath('./%s_fmap.nii.gz' % fname)
image = nb.load(in_file)
data = (image.get_data().astype(np.float32) / (2. * math.pi * delta_te))
nb.Nifti1Image(data, image.affine, image.header).to_filename(out_file)
return out_file
def _delta_te(in_values, te1=None, te2=None):
if isinstance(in_values, float):
te2 = in_values
te1 = 0.
if isinstance(in_values, dict):
te1 = in_values.get('EchoTime1')
te2 = in_values.get('EchoTime2')
if not all((te1, te2)):
te2 = in_values.get('EchoTimeDifference')
te1 = 0
if isinstance(in_values, list):
te2, te1 = in_values
if isinstance(te1, list):
te1 = te1[1]
if isinstance(te2, list):
te2 = te2[1]
if te1 is None or te2 is None:
raise RuntimeError(
'No echo time information found')
return abs(float(te2)-float(te1))