Merge remote-tracking branch 'upstream/master' into enh/ReviseResourceProfiler

oesteban · oesteban · commit 8a887d33faf1 · 2017-10-03T11:35:06.000-07:00
diff --git a/nipype/algorithms/rapidart.py b/nipype/algorithms/rapidart.py
@@ -99,6 +99,29 @@ def _calc_norm(mc, use_differences, source, brain_pts=None):
 
     """
 
+    affines = [_get_affine_matrix(mc[i, :], source)
+               for i in range(mc.shape[0])]
+    return _calc_norm_affine(affines, use_differences, brain_pts)
+
+
+def _calc_norm_affine(affines, use_differences, brain_pts=None):
+    """Calculates the maximum overall displacement of the midpoints
+    of the faces of a cube due to translation and rotation.
+
+    Parameters
+    ----------
+    affines : list of [4 x 4] affine matrices
+    use_differences : boolean
+    brain_pts : [4 x n_points] of coordinates
+
+    Returns
+    -------
+
+    norm : at each time point
+    displacement : euclidean distance (mm) of displacement at each coordinate
+
+    """
+
     if brain_pts is None:
         respos = np.diag([70, 70, 75])
         resneg = np.diag([-70, -110, -45])
@@ -107,49 +130,34 @@ def _calc_norm(mc, use_differences, source, brain_pts=None):
     else:
         all_pts = brain_pts
     n_pts = all_pts.size - all_pts.shape[1]
-    newpos = np.zeros((mc.shape[0], n_pts))
+    newpos = np.zeros((len(affines), n_pts))
     if brain_pts is not None:
-        displacement = np.zeros((mc.shape[0], int(n_pts / 3)))
-    for i in range(mc.shape[0]):
-        affine = _get_affine_matrix(mc[i, :], source)
-        newpos[i, :] = np.dot(affine,
-                              all_pts)[0:3, :].ravel()
+        displacement = np.zeros((len(affines), int(n_pts / 3)))
+    for i, affine in enumerate(affines):
+        newpos[i, :] = np.dot(affine, all_pts)[0:3, :].ravel()
         if brain_pts is not None:
-            displacement[i, :] = \
-                np.sqrt(np.sum(np.power(np.reshape(newpos[i, :],
-                                                   (3, all_pts.shape[1])) -
-                                        all_pts[0:3, :],
-                                        2),
-                               axis=0))
+            displacement[i, :] = np.sqrt(np.sum(
+                np.power(np.reshape(newpos[i, :],
+                                    (3, all_pts.shape[1])) - all_pts[0:3, :],
+                         2),
+                axis=0))
     # np.savez('displacement.npz', newpos=newpos, pts=all_pts)
-    normdata = np.zeros(mc.shape[0])
+    normdata = np.zeros(len(affines))
     if use_differences:
         newpos = np.concatenate((np.zeros((1, n_pts)),
                                  np.diff(newpos, n=1, axis=0)), axis=0)
         for i in range(newpos.shape[0]):
             normdata[i] = \
-                np.max(np.sqrt(np.sum(np.reshape(np.power(np.abs(newpos[i, :]), 2),
-                                                 (3, all_pts.shape[1])), axis=0)))
+                np.max(np.sqrt(np.sum(
+                    np.reshape(np.power(np.abs(newpos[i, :]), 2),
+                               (3, all_pts.shape[1])),
+                    axis=0)))
     else:
         newpos = np.abs(signal.detrend(newpos, axis=0, type='constant'))
         normdata = np.sqrt(np.mean(np.power(newpos, 2), axis=1))
     return normdata, displacement
 
 
-def _nanmean(a, axis=None):
-    """Return the mean excluding items that are nan
-
-    >>> a = [1, 2, np.nan]
-    >>> _nanmean(a)
-    1.5
-
-    """
-    if axis:
-        return np.nansum(a, axis) / np.sum(1 - np.isnan(a), axis)
-    else:
-        return np.nansum(a) / np.sum(1 - np.isnan(a))
-
-
 class ArtifactDetectInputSpec(BaseInterfaceInputSpec):
     realigned_files = InputMultiPath(File(exists=True),
                                      desc="Names of realigned functional data files",
@@ -376,11 +384,11 @@ def _detect_outliers_core(self, imgfile, motionfile, runidx, cwd=None):
                     vol = data[:, :, :, t0]
                     # Use an SPM like approach
                     mask_tmp = vol > \
-                        (_nanmean(vol) / self.inputs.global_threshold)
+                        (np.nanmean(vol) / self.inputs.global_threshold)
                     mask = mask * mask_tmp
                 for t0 in range(timepoints):
                     vol = data[:, :, :, t0]
-                    g[t0] = _nanmean(vol[mask])
+                    g[t0] = np.nanmean(vol[mask])
                 if len(find_indices(mask)) < (np.prod((x, y, z)) / 10):
                     intersect_mask = False
                     g = np.zeros((timepoints, 1))
@@ -390,7 +398,7 @@ def _detect_outliers_core(self, imgfile, motionfile, runidx, cwd=None):
                 for t0 in range(timepoints):
                     vol = data[:, :, :, t0]
                     mask_tmp = vol > \
-                        (_nanmean(vol) / self.inputs.global_threshold)
+                        (np.nanmean(vol) / self.inputs.global_threshold)
                     mask[:, :, :, t0] = mask_tmp
                     g[t0] = np.nansum(vol * mask_tmp) / np.nansum(mask_tmp)
         elif masktype == 'file':  # uses a mask image to determine intensity
@@ -400,15 +408,15 @@ def _detect_outliers_core(self, imgfile, motionfile, runidx, cwd=None):
             mask = mask > 0.5
             for t0 in range(timepoints):
                 vol = data[:, :, :, t0]
-                g[t0] = _nanmean(vol[mask])
+                g[t0] = np.nanmean(vol[mask])
         elif masktype == 'thresh':  # uses a fixed signal threshold
             for t0 in range(timepoints):
                 vol = data[:, :, :, t0]
                 mask = vol > self.inputs.mask_threshold
-                g[t0] = _nanmean(vol[mask])
+                g[t0] = np.nanmean(vol[mask])
         else:
             mask = np.ones((x, y, z))
-            g = _nanmean(data[mask > 0, :], 1)
+            g = np.nanmean(data[mask > 0, :], 1)
 
         # compute normalized intensity values
         gz = signal.detrend(g, axis=0)  # detrend the signal
diff --git a/nipype/interfaces/afni/__init__.py b/nipype/interfaces/afni/__init__.py
@@ -20,7 +20,7 @@
                          TShift, Volreg, Warp, QwarpPlusMinus, Qwarp)
 from .svm import (SVMTest, SVMTrain)
 from .utils import (ABoverlap, AFNItoNIFTI, Autobox, Axialize, BrickStat,
-                    Bucket, Calc, Cat, CatMatvec, Copy, Dot,
+                    Bucket, Calc, Cat, CatMatvec, CenterMass, Copy, Dot,
                     Edge3, Eval, FWHMx, MaskTool, Merge, Notes, NwarpApply,
                     OneDToolPy,
                     Refit, Resample, TCat, TCatSubBrick, TStat, To3D, Unifize,
diff --git a/nipype/interfaces/afni/tests/test_auto_CenterMass.py b/nipype/interfaces/afni/tests/test_auto_CenterMass.py
@@ -0,0 +1,59 @@
+# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT
+from __future__ import unicode_literals
+from ..utils import CenterMass
+
+
+def test_CenterMass_inputs():
+    input_map = dict(all_rois=dict(argstr='-all_rois',
+    ),
+    args=dict(argstr='%s',
+    ),
+    automask=dict(argstr='-automask',
+    ),
+    cm_file=dict(argstr='> %s',
+    descr='File to write center of mass to',
+    hash_files=False,
+    keep_extension=False,
+    name_source='in_file',
+    name_template='%s_cm.out',
+    position=-1,
+    ),
+    environ=dict(nohash=True,
+    usedefault=True,
+    ),
+    ignore_exception=dict(nohash=True,
+    usedefault=True,
+    ),
+    in_file=dict(argstr='%s',
+    copyfile=True,
+    mandatory=True,
+    position=-2,
+    ),
+    local_ijk=dict(argstr='-local_ijk',
+    ),
+    mask_file=dict(argstr='-mask %s',
+    ),
+    roi_vals=dict(argstr='-roi_vals %s',
+    ),
+    set_cm=dict(argstr='-set %f %f %f',
+    ),
+    terminal_output=dict(nohash=True,
+    ),
+    )
+    inputs = CenterMass.input_spec()
+
+    for key, metadata in list(input_map.items()):
+        for metakey, value in list(metadata.items()):
+            assert getattr(inputs.traits()[key], metakey) == value
+
+
+def test_CenterMass_outputs():
+    output_map = dict(cm=dict(),
+    cm_file=dict(),
+    out_file=dict(),
+    )
+    outputs = CenterMass.output_spec()
+
+    for key, metadata in list(output_map.items()):
+        for metakey, value in list(metadata.items()):
+            assert getattr(outputs.traits()[key], metakey) == value
diff --git a/nipype/interfaces/afni/utils.py b/nipype/interfaces/afni/utils.py
@@ -642,6 +642,106 @@ def _format_arg(self, name, spec, value):
             return spec.argstr%(' '.join([i[0]+' -'+i[1] for i in value]))
         return super(CatMatvec, self)._format_arg(name, spec, value)
 
+
+class CenterMassInputSpec(CommandLineInputSpec):
+    in_file = File(
+        desc='input file to 3dCM',
+        argstr='%s',
+        position=-2,
+        mandatory=True,
+        exists=True,
+        copyfile=True)
+    cm_file = File(
+         name_source='in_file',
+         name_template='%s_cm.out',
+         hash_files=False,
+         keep_extension=False,
+         descr="File to write center of mass to",
+         argstr="> %s",
+         position=-1)
+    mask_file = File(
+        desc='Only voxels with nonzero values in the provided mask will be '
+             'averaged.',
+        argstr='-mask %s',
+        exists=True)
+    automask = traits.Bool(
+        desc='Generate the mask automatically',
+        argstr='-automask')
+    set_cm = traits.Tuple(
+        (traits.Float(), traits.Float(), traits.Float()),
+        desc='After computing the center of mass, set the origin fields in '
+             'the header so that the center of mass will be at (x,y,z) in '
+             'DICOM coords.',
+        argstr='-set %f %f %f')
+    local_ijk = traits.Bool(
+        desc='Output values as (i,j,k) in local orienation',
+        argstr='-local_ijk')
+    roi_vals = traits.List(
+        traits.Int,
+        desc='Compute center of mass for each blob with voxel value of v0, '
+             'v1, v2, etc. This option is handy for getting ROI centers of '
+             'mass.',
+        argstr='-roi_vals %s')
+    all_rois = traits.Bool(
+        desc='Don\'t bother listing the values of ROIs you want: The program '
+             'will find all of them and produce a full list',
+        argstr='-all_rois')
+
+
+class CenterMassOutputSpec(TraitedSpec):
+    out_file = File(
+        exists=True,
+        desc='output file')
+    cm_file = File(
+        desc='file with the center of mass coordinates')
+    cm = traits.Either(
+        traits.Tuple(traits.Float(), traits.Float(), traits.Float()),
+        traits.List(traits.Tuple(traits.Float(), traits.Float(),
+                                 traits.Float())),
+        desc='center of mass')
+
+
+class CenterMass(AFNICommandBase):
+    """Computes center of mass using 3dCM command
+
+    .. note::
+
+      By default, the output is (x,y,z) values in DICOM coordinates. But
+      as of Dec, 2016, there are now command line switches for other options.
+
+
+    For complete details, see the `3dCM Documentation.
+    <https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dCM.html>`_
+
+    Examples
+    ========
+
+    >>> from nipype.interfaces import afni
+    >>> cm = afni.CenterMass()
+    >>> cm.inputs.in_file = 'structural.nii'
+    >>> cm.inputs.cm_file = 'cm.txt'
+    >>> cm.inputs.roi_vals = [2, 10]
+    >>> cm.cmdline  # doctest: +ALLOW_UNICODE
+    '3dCM -roi_vals 2 10 structural.nii > cm.txt'
+    >>> res = 3dcm.run()  # doctest: +SKIP
+    """
+
+    _cmd = '3dCM'
+    input_spec = CenterMassInputSpec
+    output_spec = CenterMassOutputSpec
+
+    def _list_outputs(self):
+        outputs = super(CenterMass, self)._list_outputs()
+        outputs['out_file'] = os.path.abspath(self.inputs.in_file)
+        outputs['cm_file'] = os.path.abspath(self.inputs.cm_file)
+        sout = np.loadtxt(outputs['cm_file'])  # pylint: disable=E1101
+        if len(sout) > 1:
+            outputs['cm'] = [tuple(s) for s in sout]
+        else:
+            outputs['cm'] = tuple(sout)
+        return outputs
+
+
 class CopyInputSpec(AFNICommandInputSpec):
     in_file = File(
         desc='input file to 3dcopy',
