ENH: Base implementation for phase1/2 fieldmaps

Putting together the lessons learned in nipreps#30, leveraging nipreps#52 and nipreps#53
(unfolded from nipreps#30 too), and utilizing nipreps#50 and nipreps#51, this workflow adds
the phase difference map calculation, considering it one use-case of the
general phase-difference fieldmap workflow.

On top of this PR, we can continue the discussions held in nipreps#30.
Probably, we will want to address nipreps#23 the first - the magnitude
segmentation is sometimes really bad (e.g. see the phase1/2 unit test).

Another discussion arisen in nipreps#30 is the spatial smoothing of the
fieldmap (nipreps#22).

Finally, the plan is to revise this implementation and determine whether
the subtraction should happen before or after PRELUDE, and whether the
arctan2 route is more interesting.

sdcflows/interfaces/fmap.py

@@ -24,6 +24,43 @@
 LOGGER = logging.getLogger('nipype.interface')
 
 
+class _SubtractPhasesInputSpec(BaseInterfaceInputSpec):
+    in_phases = traits.List(File(exists=True), min=1, max=2,
+                            desc='input phase maps')
+    in_meta = traits.List(traits.Dict(), min=1, max=2,
+                          desc='metadata corresponding to the inputs')
+
+
+class _SubtractPhasesOutputSpec(TraitedSpec):
+    phase_diff = File(exists=True, desc='phase difference map')
+    metadata = traits.Dict(desc='output metadata')
+
+
+class SubtractPhases(SimpleInterface):
+    """Calculate a phase difference map."""
+
+    input_spec = _SubtractPhasesInputSpec
+    output_spec = _SubtractPhasesOutputSpec
+
+    def _run_interface(self, runtime):
+        if len(self.inputs.in_phases) != len(self.inputs.in_meta):
+            raise ValueError(
+                'Length of input phase-difference maps and metadata files '
+                'should match.')
+
+        if len(self.inputs.in_phases) == 1:
+            self._results['phase_diff'] = self.inputs.in_phases[0]
+            self._results['metadata'] = self.inputs.in_meta[0]
+            return runtime
+
+        self._results['phase_diff'], self._results['metadata'] = \
+            _subtract_phases(self.inputs.in_phases,
+                             self.inputs.in_meta,
+                             newpath=runtime.cwd)
+
+        return runtime
+
+
 class _FieldEnhanceInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc='input fieldmap')
     in_mask = File(exists=True, desc='brain mask')
@@ -647,3 +684,36 @@ def au2rads(in_file, newpath=None):
     out_file = fname_presuffix(in_file, suffix='_rads', newpath=newpath)
     nb.Nifti1Image(data, im.affine, hdr).to_filename(out_file)
     return out_file
+
+
+def _subtract_phases(in_phases, in_meta, newpath=None):
+    # Discard traits with copy(), so that pop() works.
+    in_meta = (in_meta[0].copy(), in_meta[1].copy())
+    echo_times = tuple([m.pop('EchoTime', None) for m in in_meta])
+    if not all(echo_times):
+        raise ValueError(
+            'One or more missing EchoTime metadata parameter '
+            'associated to one or more phase map(s).')
+
+    if echo_times[0] > echo_times[1]:
+        in_phases = (in_phases[1], in_phases[0])
+        in_meta = (in_meta[1], in_meta[0])
+        echo_times = (echo_times[1], echo_times[0])
+
+    in_phases_nii = [nb.load(ph) for ph in in_phases]
+    sub_data = in_phases_nii[1].get_fdata(dtype='float32') - \
+        in_phases_nii[0].get_fdata(dtype='float32')
+
+    new_meta = in_meta[1].copy()
+    new_meta.update(in_meta[0])
+    new_meta['EchoTime1'] = echo_times[0]
+    new_meta['EchoTime2'] = echo_times[1]
+
+    hdr = in_phases_nii[0].header.copy()
+    hdr.set_data_dtype(np.float32)
+    hdr.set_xyzt_units('mm')
+    nii = nb.Nifti1Image(sub_data, in_phases_nii[0].affine, hdr)
+    out_phdiff = fname_presuffix(in_phases[0], suffix='_phdiff',
+                                 newpath=newpath)
+    nii.to_filename(out_phdiff)
+    return out_phdiff, new_meta

sdcflows/workflows/phdiff.py

@@ -21,7 +21,7 @@
 from nipype.pipeline import engine as pe
 from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 
-from ..interfaces.fmap import Phasediff2Fieldmap, PhaseMap2rads
+from ..interfaces.fmap import Phasediff2Fieldmap, PhaseMap2rads, SubtractPhases
 from .gre import init_fmap_postproc_wf, init_magnitude_wf
 
 
@@ -107,6 +107,12 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
     # FSL PRELUDE will perform phase-unwrapping
     prelude = pe.MapNode(fsl.PRELUDE(), iterfield=['phase_file'], name='prelude')
 
+    calc_phdiff = pe.Node(SubtractPhases(), name='calc_phdiff',
+                          run_without_submitting=True)
+
+    workflow.connect([
+    ])
+
     fmap_postproc_wf = init_fmap_postproc_wf(omp_nthreads=omp_nthreads,
                                              fmap_bspline=False)
     compfmap = pe.Node(Phasediff2Fieldmap(), name='compfmap')
@@ -118,8 +124,11 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
                                  ('outputnode.fmap_mask', 'mask_file')]),
         (split, phmap2rads, [('map_file', 'in_file')]),
         (phmap2rads, prelude, [('out_file', 'phase_file')]),
-        (split, fmap_postproc_wf, [('meta', 'inputnode.metadata')]),
-        (prelude, fmap_postproc_wf, [('unwrapped_phase_file', 'inputnode.fmap')]),
+        (prelude, calc_phdiff, [('unwrapped_phase_file', 'in_phases')]),
+        (split, calc_phdiff, [('meta', 'in_meta')]),
+        (calc_phdiff, fmap_postproc_wf, [
+            ('phase_diff', 'inputnode.fmap'),
+            ('metadata', 'inputnode.metadata')]),
         (magnitude_wf, fmap_postproc_wf, [
             ('outputnode.fmap_mask', 'inputnode.fmap_mask'),
             ('outputnode.fmap_ref', 'inputnode.fmap_ref')]),

sdcflows/workflows/tests/test_phdiff.py

@@ -10,7 +10,7 @@
     'ds001600',
     'testdata',
 ])
-def test_workflow(bids_layouts, tmpdir, output_path, dataset, workdir):
+def test_phdiff(bids_layouts, tmpdir, output_path, dataset, workdir):
     """Test creation of the workflow."""
     tmpdir.chdir()
 
@@ -19,11 +19,56 @@ def test_workflow(bids_layouts, tmpdir, output_path, dataset, workdir):
     phdiff_wf = init_phdiff_wf(omp_nthreads=2)
     phdiff_wf.inputs.inputnode.magnitude = data.get(
         suffix=['magnitude1', 'magnitude2'],
-        acq='v4',
+        acquisition='v4',
         return_type='file',
         extension=['.nii', '.nii.gz'])
 
-    phdiff_files = data.get(suffix='phasediff', acq='v4',
+    phdiff_files = data.get(suffix='phasediff', acquisition='v4',
+                            extension=['.nii', '.nii.gz'])
+
+    phdiff_wf.inputs.inputnode.phasediff = [
+        (ph.path, ph.get_metadata()) for ph in phdiff_files]
+
+    if output_path:
+        from ...interfaces.reportlets import FieldmapReportlet
+        rep = pe.Node(FieldmapReportlet(reference_label='Magnitude'), 'simple_report')
+        rep.interface._always_run = True
+
+        dsink = pe.Node(DerivativesDataSink(
+            base_directory=str(output_path), keep_dtype=True), name='dsink')
+        dsink.interface.out_path_base = 'sdcflows'
+        dsink.inputs.source_file = phdiff_files[0].path
+
+        wf.connect([
+            (phdiff_wf, rep, [
+                ('outputnode.fmap', 'fieldmap'),
+                ('outputnode.fmap_ref', 'reference'),
+                ('outputnode.fmap_mask', 'mask')]),
+            (rep, dsink, [('out_report', 'in_file')]),
+        ])
+    else:
+        wf.add_nodes([phdiff_wf])
+
+    if workdir:
+        wf.base_dir = str(workdir)
+
+    wf.run()
+
+
+def test_phases(bids_layouts, tmpdir, output_path, workdir):
+    """Test creation of the workflow."""
+    tmpdir.chdir()
+
+    data = bids_layouts['ds001600']
+    wf = Workflow(name='phases_ds001600')
+    phdiff_wf = init_phdiff_wf(omp_nthreads=2)
+    phdiff_wf.inputs.inputnode.magnitude = data.get(
+        suffix=['magnitude1', 'magnitude2'],
+        acquisition='v2',
+        return_type='file',
+        extension=['.nii', '.nii.gz'])
+
+    phdiff_files = data.get(suffix=['phase1', 'phase2'], acquisition='v2',
                             extension=['.nii', '.nii.gz'])
 
     phdiff_wf.inputs.inputnode.phasediff = [