Merge pull request #2825 from richford/fsl-eddy-quad-interface

ENH: Add FSL `eddy_quad` interface

Author: effigies

.zenodo.json

@@ -589,6 +589,11 @@
       "affiliation": "MIT, HMS",
       "name": "Ghosh, Satrajit",
       "orcid": "0000-0002-5312-6729"
+    },
+    {
+      "affiliation": "University of Washington",
+      "name": "Richie-Halford, Adam",
+      "orcid": "0000-0001-9276-9084"
     }
   ],
   "keywords": [

nipype/interfaces/fsl/__init__.py

@@ -21,7 +21,7 @@
     WarpPointsFromStd, RobustFOV, CopyGeom, MotionOutliers)
 
 from .epi import (PrepareFieldmap, TOPUP, ApplyTOPUP, Eddy, EPIDeWarp, SigLoss,
-                  EddyCorrect, EpiReg)
+                  EddyCorrect, EpiReg, EddyQuad)
 from .dti import (BEDPOSTX, XFibres, DTIFit, ProbTrackX, ProbTrackX2, VecReg,
                   ProjThresh, FindTheBiggest, DistanceMap, TractSkeleton,
                   MakeDyadicVectors, BEDPOSTX5, XFibres5)

nipype/interfaces/fsl/epi.py

@@ -1232,3 +1232,193 @@ def _run_interface(self, runtime):
         if runtime.stderr:
             self.raise_exception(runtime)
         return runtime
+
+
+class EddyQuadInputSpec(FSLCommandInputSpec):
+    base_name = traits.Str(
+        'eddy_corrected',
+        usedefault=True,
+        argstr='%s',
+        desc=("Basename (including path) for EDDY output files, i.e., "
+              "corrected images and QC files"),
+        position=0,
+    )
+    idx_file = File(
+        exists=True,
+        mandatory=True,
+        argstr="--eddyIdx=%s",
+        desc=("File containing indices for all volumes into acquisition "
+              "parameters")
+    )
+    param_file = File(
+        exists=True,
+        mandatory=True,
+        argstr="--eddyParams=%s",
+        desc="File containing acquisition parameters"
+    )
+    mask_file = File(
+        exists=True,
+        mandatory=True,
+        argstr="--mask=%s",
+        desc="Binary mask file"
+    )
+    bval_file = File(
+        exists=True,
+        mandatory=True,
+        argstr="--bvals=%s",
+        desc="b-values file"
+    )
+    bvec_file = File(
+        exists=True,
+        argstr="--bvecs=%s",
+        desc=("b-vectors file - only used when <base_name>.eddy_residuals "
+              "file is present")
+    )
+    output_dir = traits.Str(
+        name_template='%s.qc',
+        name_source=['base_name'],
+        argstr='--output-dir=%s',
+        desc="Output directory - default = '<base_name>.qc'",
+    )
+    field = File(
+        exists=True,
+        argstr='--field=%s',
+        desc="TOPUP estimated field (in Hz)",
+    )
+    slice_spec = File(
+        exists=True,
+        argstr='--slspec=%s',
+        desc="Text file specifying slice/group acquisition",
+    )
+    verbose = traits.Bool(
+        argstr='--verbose',
+        desc="Display debug messages",
+    )
+
+
+class EddyQuadOutputSpec(TraitedSpec):
+    qc_json = File(
+        exists=True,
+        desc=("Single subject database containing quality metrics and data "
+              "info.")
+    )
+    qc_pdf = File(
+        exists=True,
+        desc="Single subject QC report."
+    )
+    avg_b_png = traits.List(
+        File(exists=True),
+        desc=("Image showing mid-sagittal, -coronal and -axial slices of "
+              "each averaged b-shell volume.")
+    )
+    avg_b0_pe_png = traits.List(
+        File(exists=True),
+        desc=("Image showing mid-sagittal, -coronal and -axial slices of "
+              "each averaged pe-direction b0 volume. Generated when using "
+              "the -f option.")
+    )
+    cnr_png = traits.List(
+        File(exists=True),
+        desc=("Image showing mid-sagittal, -coronal and -axial slices of "
+              "each b-shell CNR volume. Generated when CNR maps are "
+              "available.")
+    )
+    vdm_png = File(
+        exists=True,
+        desc=("Image showing mid-sagittal, -coronal and -axial slices of "
+              "the voxel displacement map. Generated when using the -f "
+              "option.")
+    )
+    residuals = File(
+        exists=True,
+        desc=("Text file containing the volume-wise mask-averaged squared "
+              "residuals. Generated when residual maps are available.")
+    )
+    clean_volumes = File(
+        exists=True,
+        desc=("Text file containing a list of clean volumes, based on "
+              "the eddy squared residuals. To generate a version of the "
+              "pre-processed dataset without outlier volumes, use: "
+              "`fslselectvols -i <eddy_corrected_data> -o "
+              "eddy_corrected_data_clean --vols=vols_no_outliers.txt`")
+    )
+
+
+class EddyQuad(FSLCommand):
+    """
+    Interface for FSL eddy_quad, a tool for generating single subject reports
+    and storing the quality assessment indices for each subject.
+    `User guide <https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddyqc/UsersGuide>`_
+
+    Examples
+    --------
+
+    >>> from nipype.interfaces.fsl import EddyQuad
+    >>> quad = EddyQuad()
+    >>> quad.inputs.base_name  = 'eddy_corrected'
+    >>> quad.inputs.idx_file   = 'epi_index.txt'
+    >>> quad.inputs.param_file = 'epi_acqp.txt'
+    >>> quad.inputs.mask_file  = 'epi_mask.nii'
+    >>> quad.inputs.bval_file  = 'bvals.scheme'
+    >>> quad.inputs.bvec_file  = 'bvecs.scheme'
+    >>> quad.inputs.output_dir = 'eddy_corrected.qc'
+    >>> quad.inputs.field      = 'fieldmap_phase_fslprepared.nii'
+    >>> quad.inputs.verbose    = True
+    >>> quad.cmdline
+    'eddy_quad eddy_corrected --bvals=bvals.scheme --bvecs=bvecs.scheme \
+--field=fieldmap_phase_fslprepared.nii --eddyIdx=epi_index.txt \
+--mask=epi_mask.nii --output-dir=eddy_corrected.qc --eddyParams=epi_acqp.txt \
+--verbose'
+    >>> res = quad.run() # doctest: +SKIP
+
+    """
+    _cmd = 'eddy_quad'
+    input_spec = EddyQuadInputSpec
+    output_spec = EddyQuadOutputSpec
+
+    def _list_outputs(self):
+        from glob import glob
+        outputs = self.output_spec().get()
+
+        # If the output directory isn't defined, the interface seems to use
+        # the default but not set its value in `self.inputs.output_dir`
+        if not isdefined(self.inputs.output_dir):
+            out_dir = os.path.abspath(os.path.basename(self.inputs.base_name) + '.qc.nii.gz')
+        else:
+            out_dir = os.path.abspath(self.inputs.output_dir)
+
+        outputs['qc_json'] = os.path.join(out_dir, 'qc.json')
+        outputs['qc_pdf'] = os.path.join(out_dir, 'qc.pdf')
+
+        # Grab all b* files here. This will also grab the b0_pe* files
+        # as well, but only if the field input was provided. So we'll remove
+        # them later in the next conditional.
+        outputs['avg_b_png'] = sorted(glob(
+            os.path.join(out_dir, 'avg_b*.png')
+        ))
+
+        if isdefined(self.inputs.field):
+            outputs['avg_b0_pe_png'] = sorted(glob(
+                os.path.join(out_dir, 'avg_b0_pe*.png')
+            ))
+
+            # The previous glob for `avg_b_png` also grabbed the
+            # `avg_b0_pe_png` files so we have to remove them
+            # from `avg_b_png`.
+            for fname in outputs['avg_b0_pe_png']:
+                outputs['avg_b_png'].remove(fname)
+
+            outputs['vdm_png'] = os.path.join(out_dir, 'vdm.png')
+
+        outputs['cnr_png'] = sorted(glob(os.path.join(out_dir, 'cnr*.png')))
+
+        residuals = os.path.join(out_dir, 'eddy_msr.txt')
+        if os.path.isfile(residuals):
+            outputs['residuals'] = residuals
+
+        clean_volumes = os.path.join(out_dir, 'vols_no_outliers.txt')
+        if os.path.isfile(clean_volumes):
+            outputs['clean_volumes'] = clean_volumes
+
+        return outputs
+