ENH: Refactor fieldmap-unwarping flows, more homogeneous interface

This PR closes #19, reorganizing the unwarping tools for fieldmap-based
solutions:

  - [x] Separate the displacements field generation from actual
        unwarping (i.e., applying the nonlinear transform).
        By doing this, we are not only addressing #19 directly, we are
        also readying the ground for tackling #21 since now all the
        unwarping paths (pepolar, fieldmaps/phases/phasediff, or syn)
        have more consistent interfaces.
  - [x] Update the root workflow to reflect these changes.
  - [x] Unloaded ``unwarp.py`` of the workflow to write reports, which
        has been moved to a new module called ``outputs.py`` following
        some of the latest *fMRIPrep* non-written standards.
  - [x] Polished some minor documentation and stylistic issues within
        the scope of the proposed changes.

Author: oesteban

sdcflows/workflows/base.py

@@ -151,6 +151,7 @@ def init_sdc_estimate_wf(fmaps, epi_meta, omp_nthreads=1, debug=False, ignore=No
 
     # FIELDMAP path
     elif 'fieldmap' in fmaps or 'phasediff' in fmaps:
+        from .fmap import init_fmap2field_wf
         from .unwarp import init_sdc_unwarp_wf
 
         # SyN works without this metadata
@@ -179,21 +180,28 @@ def init_sdc_estimate_wf(fmaps, epi_meta, omp_nthreads=1, debug=False, ignore=No
                 m for m, _ in fmaps['phasediff']['magnitude']]
             fmap_wf.inputs.inputnode.phasediff = fmaps['phasediff']['phases']
 
+        fmap2field_wf = init_fmap2field_wf(omp_nthreads=omp_nthreads, debug=debug)
+        fmap2field_wf.inputs.inputnode.metadata = epi_meta
+
         sdc_unwarp_wf = init_sdc_unwarp_wf(
             omp_nthreads=omp_nthreads,
             debug=debug,
             name='sdc_unwarp_wf')
-        sdc_unwarp_wf.inputs.inputnode.metadata = epi_meta
 
         workflow.connect([
+            (inputnode, fmap2field_wf, [
+                ('epi_file', 'inputnode.in_reference'),
+                ('epi_brain', 'inputnode.in_reference_brain')]),
             (inputnode, sdc_unwarp_wf, [
                 ('epi_file', 'inputnode.in_reference'),
-                ('epi_brain', 'inputnode.in_reference_brain'),
-                ('epi_mask', 'inputnode.in_mask')]),
-            (fmap_wf, sdc_unwarp_wf, [
+                ('epi_brain', 'inputnode.in_reference_brain')]),
+            (fmap_wf, fmap2field_wf, [
                 ('outputnode.fmap', 'inputnode.fmap'),
                 ('outputnode.fmap_ref', 'inputnode.fmap_ref'),
                 ('outputnode.fmap_mask', 'inputnode.fmap_mask')]),
+            (fmap2field_wf, sdc_unwarp_wf, [
+                ('outputnode.out_warp', 'inputnode.in_warp')]),
+
         ])
     elif not only_syn:
         raise ValueError('Fieldmaps of types %s are not supported' %

sdcflows/workflows/fmap.py

@@ -8,21 +8,25 @@
 Direct B0 mapping sequences
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 When the fieldmap is directly measured with a prescribed sequence (such as
-:abbr:`SE (spiral echo)`), we only need to calculate the corresponding B-Spline
-coefficients to adapt the fieldmap to the TOPUP tool.
+:abbr:`SE (spiral echo)`), we only need to calculate the corresponding
+displacements field that accounts for the distortions.
 This procedure is described with more detail `here
 <https://cni.stanford.edu/wiki/GE_Processing#Fieldmaps>`__.
 
 This corresponds to `this section of the BIDS specification
 <https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#a-real-fieldmap-image>`__.
 
 """
+import pkg_resources as pkgr
 
 from nipype.pipeline import engine as pe
 from nipype.interfaces import utility as niu, fsl
 from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+from niworkflows.interfaces import itk
 from niworkflows.interfaces.images import IntraModalMerge
+from niworkflows.interfaces.registration import ANTSApplyTransformsRPT, ANTSRegistrationRPT
 
+from ..interfaces.fmap import get_ees as _get_ees, FieldToRadS
 from .gre import init_fmap_postproc_wf, init_magnitude_wf
 
 
@@ -71,6 +75,10 @@ def init_fmap_wf(omp_nthreads, fmap_bspline, name='fmap_wf'):
 
     """
     workflow = Workflow(name=name)
+    workflow.__desc__ = """\
+A B0-nonuniformity map (or *fieldmap*) was directly measured with an MRI scheme
+designed with that purpose (typically, a spiral pulse sequence).
+"""
     inputnode = pe.Node(niu.IdentityInterface(
         fields=['magnitude', 'fieldmap']), name='inputnode')
     outputnode = pe.Node(niu.IdentityInterface(fields=['fmap', 'fmap_ref', 'fmap_mask']),
@@ -101,3 +109,138 @@ def init_fmap_wf(omp_nthreads, fmap_bspline, name='fmap_wf'):
         (fmap_postproc_wf, outputnode, [('outputnode.out_fmap', 'fmap')]),
     ])
     return workflow
+
+
+def init_fmap2field_wf(omp_nthreads, debug, name='fmap2field_wf'):
+    """
+    Convert the estimated fieldmap in Hz into a displacements field.
+
+    This workflow takes in a fieldmap and calculates the corresponding
+    displacements field (in other words, an ANTs-compatible warp file).
+
+    Workflow Graph
+        .. workflow ::
+            :graph2use: orig
+            :simple_form: yes
+
+            from sdcflows.workflows.fmap import init_fmap2field_wf
+            wf = init_fmap2field_wf(omp_nthreads=8,
+                                    debug=False)
+
+    Parameters
+    ----------
+    omp_nthreads : int
+        Maximum number of threads an individual process may use.
+    debug : bool
+        Run fast configurations of registrations.
+    name : str
+        Unique name of this workflow.
+
+    Inputs
+    ------
+    in_reference
+        the reference image
+    in_reference_brain
+        the reference image (skull-stripped)
+    metadata
+        metadata associated to the ``in_reference`` EPI input
+    fmap
+        the fieldmap in Hz
+    fmap_ref
+        the reference (anatomical) image corresponding to ``fmap``
+    fmap_mask
+        a brain mask corresponding to ``fmap``
+
+
+    Outputs
+    -------
+    out_reference
+        the ``in_reference`` after unwarping
+    out_reference_brain
+        the ``in_reference`` after unwarping and skullstripping
+    out_warp
+        the corresponding :abbr:`DFM (displacements field map)` compatible with
+        ANTs
+    out_jacobian
+        the jacobian of the field (for drop-out alleviation)
+    out_mask
+        mask of the unwarped input file
+
+    """
+    workflow = Workflow(name=name)
+    workflow.__desc__ = """\
+The *fieldmap* was then co-registered to the target EPI (echo-planar imaging)
+reference run and converted to a displacements field map (amenable to registration
+tools such as ANTs) with FSL's `fugue` and other *SDCflows* tools.
+"""
+    inputnode = pe.Node(niu.IdentityInterface(
+        fields=['in_reference', 'in_reference_brain', 'metadata',
+                'fmap_ref', 'fmap_mask', 'fmap']), name='inputnode')
+    outputnode = pe.Node(niu.IdentityInterface(
+        fields=['out_warp']), name='outputnode')
+
+    # Register the reference of the fieldmap to the reference
+    # of the target image (the one that shall be corrected)
+    ants_settings = pkgr.resource_filename('sdcflows', 'data/fmap-any_registration.json')
+    if debug:
+        ants_settings = pkgr.resource_filename(
+            'sdcflows', 'data/fmap-any_registration_testing.json')
+
+    fmap2ref_reg = pe.Node(
+        ANTSRegistrationRPT(generate_report=True, from_file=ants_settings,
+                            output_inverse_warped_image=True, output_warped_image=True),
+        name='fmap2ref_reg', n_procs=omp_nthreads)
+
+    # Map the VSM into the EPI space
+    fmap2ref_apply = pe.Node(ANTSApplyTransformsRPT(
+        generate_report=True, dimension=3, interpolation='BSpline', float=True),
+        name='fmap2ref_apply')
+
+    fmap_mask2ref_apply = pe.Node(ANTSApplyTransformsRPT(
+        generate_report=False, dimension=3, interpolation='MultiLabel',
+        float=True),
+        name='fmap_mask2ref_apply')
+
+    # Fieldmap to rads and then to voxels (VSM - voxel shift map)
+    torads = pe.Node(FieldToRadS(fmap_range=0.5), name='torads')
+
+    get_ees = pe.Node(niu.Function(function=_get_ees, output_names=['ees']), name='get_ees')
+
+    gen_vsm = pe.Node(fsl.FUGUE(save_unmasked_shift=True), name='gen_vsm')
+    # Convert the VSM into a DFM (displacements field map)
+    # or: FUGUE shift to ANTS warping.
+    vsm2dfm = pe.Node(itk.FUGUEvsm2ANTSwarp(), name='vsm2dfm')
+
+    workflow.connect([
+        (inputnode, fmap2ref_reg, [('fmap_ref', 'moving_image'),
+                                   ('in_reference_brain', 'fixed_image')]),
+        (inputnode, fmap2ref_apply, [('fmap', 'input_image'),
+                                     ('in_reference', 'reference_image')]),
+        (inputnode, fmap_mask2ref_apply, [('in_reference', 'reference_image'),
+                                          ('fmap_mask', 'input_image')]),
+        (inputnode, get_ees, [('in_reference', 'in_file'),
+                              ('metadata', 'in_meta')]),
+        (inputnode, gen_vsm, [(('metadata', _get_pedir_fugue), 'unwarp_direction')]),
+        (inputnode, vsm2dfm, [(('metadata', _get_pedir_bids), 'pe_dir')]),
+        (fmap2ref_reg, fmap2ref_apply, [
+            ('composite_transform', 'transforms')]),
+        (fmap2ref_reg, fmap_mask2ref_apply, [
+            ('composite_transform', 'transforms')]),
+        (fmap2ref_apply, torads, [('output_image', 'in_file')]),
+        (fmap_mask2ref_apply, gen_vsm, [('output_image', 'mask_file')]),
+        (get_ees, gen_vsm, [('ees', 'dwell_time')]),
+        (gen_vsm, vsm2dfm, [('shift_out_file', 'in_file')]),
+        (torads, gen_vsm, [('out_file', 'fmap_in_file')]),
+        (vsm2dfm, outputnode, [('out_file', 'out_warp')]),
+    ])
+    return workflow
+
+
+# Helper functions
+# ------------------------------------------------------------
+def _get_pedir_bids(in_dict):
+    return in_dict['PhaseEncodingDirection']
+
+
+def _get_pedir_fugue(in_dict):
+    return in_dict['PhaseEncodingDirection'].replace('i', 'x').replace('j', 'y').replace('k', 'z')

sdcflows/workflows/outputs.py

@@ -0,0 +1,80 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""Writing out outputs."""
+from nipype.pipeline import engine as pe
+from nipype.interfaces import utility as niu
+from niworkflows.interfaces.bids import DerivativesDataSink
+
+
+def init_sdc_unwarp_report_wf(name='sdc_unwarp_report_wf', forcedsyn=False):
+    """
+    Save a reportlet showing how SDC unwarping performed.
+
+    This workflow generates and saves a reportlet showing the effect of fieldmap
+    unwarping a BOLD image.
+
+    .. workflow::
+        :graph2use: orig
+        :simple_form: yes
+
+        from sdcflows.workflows.outputs import init_sdc_unwarp_report_wf
+        wf = init_sdc_unwarp_report_wf()
+
+    Parameters
+    ----------
+    name : str, optional
+        Workflow name (default: ``sdc_unwarp_report_wf``)
+    forcedsyn : bool, optional
+        Whether SyN-SDC was forced.
+
+    Inputs
+    ------
+    in_pre
+        Reference image, before unwarping
+    in_post
+        Reference image, after unwarping
+    in_seg
+        Segmentation of preprocessed structural image, including
+        gray-matter (GM), white-matter (WM) and cerebrospinal fluid (CSF)
+    in_xfm
+        Affine transform from T1 space to BOLD space (ITK format)
+
+    """
+    from niworkflows.interfaces import SimpleBeforeAfter
+    from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms
+    from niworkflows.interfaces.images import extract_wm
+
+    DEFAULT_MEMORY_MIN_GB = 0.01
+
+    workflow = pe.Workflow(name=name)
+
+    inputnode = pe.Node(niu.IdentityInterface(
+        fields=['in_pre', 'in_post', 'in_seg', 'in_xfm']), name='inputnode')
+
+    map_seg = pe.Node(ApplyTransforms(
+        dimension=3, float=True, interpolation='MultiLabel'),
+        name='map_seg', mem_gb=0.3)
+
+    sel_wm = pe.Node(niu.Function(function=extract_wm), name='sel_wm',
+                     mem_gb=DEFAULT_MEMORY_MIN_GB)
+
+    bold_rpt = pe.Node(SimpleBeforeAfter(), name='bold_rpt',
+                       mem_gb=0.1)
+    ds_report_sdc = pe.Node(
+        DerivativesDataSink(desc='sdc' if not forcedsyn else 'forcedsyn',
+                            suffix='bold'), name='ds_report_sdc',
+        mem_gb=DEFAULT_MEMORY_MIN_GB, run_without_submitting=True
+    )
+
+    workflow.connect([
+        (inputnode, bold_rpt, [('in_post', 'after'),
+                               ('in_pre', 'before')]),
+        (bold_rpt, ds_report_sdc, [('out_report', 'in_file')]),
+        (inputnode, map_seg, [('in_post', 'reference_image'),
+                              ('in_seg', 'input_image'),
+                              ('in_xfm', 'transforms')]),
+        (map_seg, sel_wm, [('output_image', 'in_seg')]),
+        (sel_wm, bold_rpt, [('out', 'wm_seg')]),
+    ])
+
+    return workflow

sdcflows/workflows/pepolar.py

@@ -103,9 +103,9 @@ def init_pepolar_unwarp_wf(omp_nthreads=1, matched_pe=False,
     """
     workflow = Workflow(name=name)
     workflow.__desc__ = """\
-A deformation field to correct for susceptibility distortions was estimated
-based on two echo-planar imaging (EPI) references with opposing phase-encoding
-directions, using `3dQwarp` @afni (AFNI {afni_ver}).
+A B0-nonuniformity map (or *fieldmap*) was estimated based on two (or more)
+echo-planar imaging (EPI) references with opposing phase-encoding
+directions, with `3dQwarp` @afni (AFNI {afni_ver}).
 """.format(afni_ver=''.join(['%02d' % v for v in afni.Info().version() or []]))
 
     inputnode = pe.Node(niu.IdentityInterface(

sdcflows/workflows/phdiff.py

@@ -33,7 +33,7 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
     magnitude images corresponding to two or more :abbr:`GRE (Gradient Echo sequence)`
     acquisitions.
     The most delicate bit of this workflow is the phase-unwrapping process: phase maps
-    are clipped in the range :math:`[0 \dotsb 2 \cdot \pi )`.
+    are clipped in the range :math:`[0 \dotsb 2\pi )`.
     To find the integer number of offsets that make a region continously smooth with
     its neighbour, FSL PRELUDE is run [Jenkinson2003]_.
     FSL PRELUDE takes wrapped maps in the range 0 to 6.28, `as per the user guide
@@ -83,11 +83,11 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
     """
     workflow = Workflow(name=name)
     workflow.__desc__ = """\
-A deformation field to correct for susceptibility distortions was estimated
-based on a field map that was co-registered to the EPI (echo-planar imaging) reference
-run, using a custom workflow of *SDCFlows* derived from D. Greve's `epidewarp.fsl`
-[script](http://www.nmr.mgh.harvard.edu/~greve/fbirn/b0/epidewarp.fsl) and
-further improvements of HCP Pipelines [@hcppipelines].
+A B0-nonuniformity map (or *fieldmap*) was estimated based on a phase-difference map
+calculated with a dual-echo GRE (gradient-recall echo) sequence, processed with a
+custom workflow of *SDCFlows* inspired by the
+[`epidewarp.fsl` script](http://www.nmr.mgh.harvard.edu/~greve/fbirn/b0/epidewarp.fsl)
+and further improvements in HCP Pipelines [@hcppipelines].
 """
 
     inputnode = pe.Node(niu.IdentityInterface(fields=['magnitude', 'phasediff']),