enh: add a fieldmap report to replace old 3 reports.

Author: oesteban

sdcflows/interfaces/fmap.py

@@ -230,6 +230,63 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _FUGUEvsm2ANTSwarpInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True,
+                   desc='input displacements field map')
+    pe_dir = traits.Enum('i', 'i-', 'j', 'j-', 'k', 'k-',
+                         desc='phase-encoding axis')
+
+
+class _FUGUEvsm2ANTSwarpOutputSpec(TraitedSpec):
+    out_file = File(desc='the output warp field')
+    fieldmap = File(desc='field map in mm')
+
+
+class FUGUEvsm2ANTSwarp(SimpleInterface):
+    """Convert a voxel-shift-map to ants warp."""
+
+    _dtype = '<f4'
+    input_spec = _FUGUEvsm2ANTSwarpInputSpec
+    output_spec = _FUGUEvsm2ANTSwarpOutputSpec
+
+    def _run_interface(self, runtime):
+        phaseEncDim = {'i': 0, 'j': 1, 'k': 2}[self.inputs.pe_dir[0]]
+        phaseEncSign = [1.0, -1.0][len(self.inputs.pe_dir) != 2]
+
+        # Create new header
+        nii = nb.load(self.inputs.in_file)
+        hdr = nii.header.copy()
+        hdr.set_data_dtype(self._dtype)
+
+        # Get data, convert to mm
+        data = nii.get_fdata(dtype=self._dtype)
+        aff = np.diag([1.0, 1.0, -1.0])
+        if np.linalg.det(aff) < 0 and phaseEncDim != 0:
+            # Reverse direction since ITK is LPS
+            aff *= -1.0
+
+        aff = aff.dot(nii.affine[:3, :3])
+        data *= phaseEncSign * nii.header.get_zooms()[phaseEncDim]
+        self._results['fieldmap'] = fname_presuffix(
+            self.inputs.in_file, suffix='_units-mm_fieldmap', newpath=runtime.cwd)
+        nb.Nifti1Image(data, nii.affine, hdr).to_filename(self._results['fieldmap'])
+
+        # Compose a vector field
+        zeros = np.zeros_like(data, dtype=self._dtype)
+        field = [zeros, zeros]
+        field.insert(phaseEncDim, data)
+        field = np.stack(field, -1)
+
+        hdr.set_intent('vector', (), '')
+        # Write out
+        self._results['out_file'] = fname_presuffix(
+            self.inputs.in_file, suffix='_desc-field_sdcwarp', newpath=runtime.cwd)
+        nb.Nifti1Image(field[:, :, :, np.newaxis, :], nii.affine, hdr).to_filename(
+            self._results['out_file'])
+
+        return runtime
+
+
 def _despike2d(data, thres, neigh=None):
     """Despike axial slices, as done in FSL's ``epiunwarp``."""
     if neigh is None:

sdcflows/workflows/fmap.py

@@ -22,11 +22,10 @@
 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 ..interfaces.fmap import get_ees as _get_ees, FieldToRadS, FUGUEvsm2ANTSwarp
 from .gre import init_fmap_postproc_wf, init_magnitude_wf
 
 
@@ -111,7 +110,8 @@ def init_fmap_wf(omp_nthreads, fmap_bspline, name='fmap_wf'):
     return workflow
 
 
-def init_fmap2field_wf(omp_nthreads, debug, name='fmap2field_wf'):
+def init_fmap2field_wf(omp_nthreads, debug, name='fmap2field_wf',
+                       generate_report=True):
     """
     Convert the estimated fieldmap in Hz into a displacements field.
 
@@ -187,13 +187,13 @@ def init_fmap2field_wf(omp_nthreads, debug, name='fmap2field_wf'):
             'sdcflows', 'data/fmap-any_registration_testing.json')
 
     fmap2ref_reg = pe.Node(
-        ANTSRegistrationRPT(generate_report=True, from_file=ants_settings,
+        ANTSRegistrationRPT(generate_report=False, 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),
+        generate_report=False, dimension=3, interpolation='BSpline', float=True),
         name='fmap2ref_apply')
 
     fmap_mask2ref_apply = pe.Node(ANTSApplyTransformsRPT(
@@ -209,7 +209,7 @@ def init_fmap2field_wf(omp_nthreads, debug, name='fmap2field_wf'):
     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')
+    vsm2dfm = pe.Node(FUGUEvsm2ANTSwarp(), name='vsm2dfm')
 
     workflow.connect([
         (inputnode, fmap2ref_reg, [('fmap_ref', 'moving_image'),
@@ -228,11 +228,32 @@ def init_fmap2field_wf(omp_nthreads, debug, name='fmap2field_wf'):
             ('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')]),
+        (get_ees, gen_vsm, [('ees', 'dwell_time')]),
         (torads, gen_vsm, [('out_file', 'fmap_in_file')]),
         (vsm2dfm, outputnode, [('out_file', 'out_warp')]),
     ])
+
+    if generate_report:
+        from niworkflows.interfaces.bids import DerivativesDataSink
+        from ..interfaces.reportlets import FieldmapReportlet
+
+        fmap_rpt = pe.Node(FieldmapReportlet(
+            reference_label='EPI Reference',
+            moving_label='Magnitude', show='both'), name='fmap_rpt')
+        ds_report_sdc = pe.Node(
+            DerivativesDataSink(desc='fieldmap', suffix='bold'),
+            name='ds_report_fmap', mem_gb=0.01, run_without_submitting=True
+        )
+
+        workflow.connect([
+            (inputnode, fmap_rpt, [('in_reference', 'reference')]),
+            (fmap2ref_reg, fmap_rpt, [('warped_image', 'moving')]),
+            (fmap_mask2ref_apply, fmap_rpt, [('output_image', 'mask')]),
+            (vsm2dfm, fmap_rpt, [('fieldmap', 'fieldmap')]),
+            (fmap_rpt, ds_report_sdc, [('out_report', 'in_file')]),
+        ])
+
     return workflow
 
 

sdcflows/workflows/gre.py

@@ -84,7 +84,7 @@ def init_magnitude_wf(omp_nthreads, name='magnitude_wf'):
     return workflow
 
 
-def init_fmap_postproc_wf(omp_nthreads, fmap_bspline, median_kernel_size=3,
+def init_fmap_postproc_wf(omp_nthreads, fmap_bspline, median_kernel_size=5,
                           name='fmap_postproc_wf'):
     """
     Postprocess a B0 map estimated elsewhere.

sdcflows/workflows/outputs.py

@@ -13,12 +13,13 @@ def init_sdc_unwarp_report_wf(name='sdc_unwarp_report_wf', forcedsyn=False):
     This workflow generates and saves a reportlet showing the effect of fieldmap
     unwarping a BOLD image.
 
-    .. workflow::
-        :graph2use: orig
-        :simple_form: yes
+    Workflow Graph
+        .. workflow::
+            :graph2use: orig
+            :simple_form: yes
 
-        from sdcflows.workflows.outputs import init_sdc_unwarp_report_wf
-        wf = init_sdc_unwarp_report_wf()
+            from sdcflows.workflows.outputs import init_sdc_unwarp_report_wf
+            wf = init_sdc_unwarp_report_wf()
 
     Parameters
     ----------