[ENH] Added afni.FWHMx interface

Author: Oscar Esteban

nipype/interfaces/afni/__init__.py

@@ -12,5 +12,5 @@
                          Fourier, Allineate, Maskave, SkullStrip, TCat, Fim,
                          BlurInMask, Autobox, TCorrMap, Bandpass, Retroicor,
                          TCorrelate, TCorr1D, BrickStat, ROIStats, AutoTcorrelate,
-                         AFNItoNIFTI, Eval, Means, Hist)
+                         AFNItoNIFTI, Eval, Means, Hist, FWHMx)
 from .svm import (SVMTest, SVMTrain)

nipype/interfaces/afni/preprocess.py

@@ -10,16 +10,18 @@
 """
 
 import os
+from sys import platform
 import os.path as op
 import re
-from warnings import warn
+import numpy as np
 
 from .base import AFNICommand, AFNICommandInputSpec, AFNICommandOutputSpec, Info, no_afni
-from ..base import CommandLineInputSpec, CommandLine, OutputMultiPath
+from ..base import CommandLineInputSpec, CommandLine
 from ..base import (Directory, TraitedSpec,
                     traits, isdefined, File, InputMultiPath, Undefined)
+from ...external.six import string_types
 from ...utils.filemanip import (load_json, save_json, split_filename)
-from ...utils.filemanip import fname_presuffix
+
 
 class To3DInputSpec(AFNICommandInputSpec):
     out_file = File(name_template="%s", desc='output image file name',
@@ -2160,3 +2162,200 @@ def _list_outputs(self):
         if not self.inputs.showhist:
             outputs['out_show'] = Undefined
         return outputs
+
+
+class FWHMxInputSpec(CommandLineInputSpec):
+    in_file = File(desc='input dataset', argstr='-input %s', mandatory=True, exists=True)
+    out_file = File(argstr='> %s', name_source='in_file', name_template='%s_fwhmx.out',
+                    position=-1, keep_extension=False, desc='output file')
+    out_subbricks = File(argstr='-out %s', name_source='in_file', name_template='%s_subbricks.out',
+                         keep_extension=False, desc='output file listing the subbricks FWHM')
+    mask = File(desc='use only voxels that are nonzero in mask', argstr='-mask %s', exists=True)
+    automask = traits.Bool(False, usedefault=True, argstr='-automask',
+                           desc='compute a mask from THIS dataset, a la 3dAutomask')
+    detrend = traits.Either(
+        traits.Bool(), traits.Int(), default=False, argstr='-detrend', xor=['demed'], usedefault=True,
+        desc='instead of demed (0th order detrending), detrend to the specified order.  If order '
+             'is not given, the program picks q=NT/30. -detrend disables -demed, and includes '
+             '-unif.')
+    demed = traits.Bool(
+        False, argstr='-demed', xorg=['detrend'],
+        desc='If the input dataset has more than one sub-brick (e.g., has a time axis), then '
+             'subtract the median of each voxel\'s time series before processing FWHM. This will '
+             'tend to remove intrinsic spatial structure and leave behind the noise.')
+    unif = traits.Bool(False, argstr='-unif',
+                       desc='If the input dataset has more than one sub-brick, then normalize each'
+                            ' voxel\'s time series to have the same MAD before processing FWHM.')
+    out_detrend = File(argstr='-detprefix %s', name_source='in_file', name_template='%s_detrend',
+                       keep_extension=False, desc='Save the detrended file into a dataset')
+    geom = traits.Bool(argstr='-geom', xor=['arith'],
+                       desc='if in_file has more than one sub-brick, compute the final estimate as'
+                            'the geometric mean of the individual sub-brick FWHM estimates')
+    arith = traits.Bool(argstr='-arith', xor=['geom'],
+                        desc='if in_file has more than one sub-brick, compute the final estimate as'
+                             'the arithmetic mean of the individual sub-brick FWHM estimates')
+    combine = traits.Bool(argstr='-combine', desc='combine the final measurements along each axis')
+    compat = traits.Bool(argstr='-compat', desc='be compatible with the older 3dFWHM')
+    acf = traits.Either(
+        traits.Bool(), File(exists=True), traits.Tuple(File(exists=True), traits.Float()),
+        argstr='-acf', desc='computes the spatial autocorrelation')
+
+
+class FWHMxOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc='output file')
+    out_subbricks = File(exists=True, desc='output file (subbricks)')
+    out_detrend = File(desc='output file, detrended')
+    fwhm = traits.Either(
+        traits.Tuple(traits.Float(), traits.Float(), traits.Float()),
+        traits.Tuple(traits.Float(), traits.Float(), traits.Float(), traits.Float()),
+        desc='FWHM along each axis')
+
+
+class FWHMx(CommandLine):
+    """
+    Unlike the older 3dFWHM, this program computes FWHMs for all sub-bricks
+    in the input dataset, each one separately.  The output for each one is
+    written to the file specified by '-out'.  The mean (arithmetic or geometric)
+    of all the FWHMs along each axis is written to stdout.  (A non-positive
+    output value indicates something bad happened; e.g., FWHM in z is meaningless
+    for a 2D dataset; the estimation method computed incoherent intermediate results.)
+
+    Examples
+    --------
+
+    >>> from nipype.interfaces import afni as afp
+    >>> fwhm = afp.FWHMx()
+    >>> fwhm.inputs.in_file = 'functional.nii'
+    >>> fwhm.cmdline
+    '3dFWHMx -input functional.nii -out functional_subbricks.out > functional_fwhmx.out'
+
+
+    (Classic) METHOD:
+
+      * Calculate ratio of variance of first differences to data variance.
+      * Should be the same as 3dFWHM for a 1-brick dataset.
+        (But the output format is simpler to use in a script.)
+
+
+    .. note:: IMPORTANT NOTE [AFNI > 16]
+
+      A completely new method for estimating and using noise smoothness values is
+      now available in 3dFWHMx and 3dClustSim. This method is implemented in the
+      '-acf' options to both programs.  'ACF' stands for (spatial) AutoCorrelation
+      Function, and it is estimated by calculating moments of differences out to
+      a larger radius than before.
+
+      Notably, real FMRI data does not actually have a Gaussian-shaped ACF, so the
+      estimated ACF is then fit (in 3dFWHMx) to a mixed model (Gaussian plus
+      mono-exponential) of the form
+
+        .. math::
+
+          ACF(r) = a * exp(-r*r/(2*b*b)) + (1-a)*exp(-r/c)
+
+
+      where :math:`r` is the radius, and :math:`a, b, c` are the fitted parameters.
+      The apparent FWHM from this model is usually somewhat larger in real data
+      than the FWHM estimated from just the nearest-neighbor differences used
+      in the 'classic' analysis.
+
+      The longer tails provided by the mono-exponential are also significant.
+      3dClustSim has also been modified to use the ACF model given above to generate
+      noise random fields.
+
+
+    .. note:: TL;DR or summary
+
+      The take-awaymessage is that the 'classic' 3dFWHMx and
+      3dClustSim analysis, using a pure Gaussian ACF, is not very correct for
+      FMRI data -- I cannot speak for PET or MEG data.
+
+
+    .. warning::
+
+      Do NOT use 3dFWHMx on the statistical results (e.g., '-bucket') from
+      3dDeconvolve or 3dREMLfit!!!  The function of 3dFWHMx is to estimate
+      the smoothness of the time series NOISE, not of the statistics. This
+      proscription is especially true if you plan to use 3dClustSim next!!
+
+
+    .. note:: Recommendations
+
+      * For FMRI statistical purposes, you DO NOT want the FWHM to reflect
+        the spatial structure of the underlying anatomy.  Rather, you want
+        the FWHM to reflect the spatial structure of the noise.  This means
+        that the input dataset should not have anatomical (spatial) structure.
+      * One good form of input is the output of '3dDeconvolve -errts', which is
+        the dataset of residuals left over after the GLM fitted signal model is
+        subtracted out from each voxel's time series.
+      * If you don't want to go to that much trouble, use '-detrend' to approximately
+        subtract out the anatomical spatial structure, OR use the output of 3dDetrend
+        for the same purpose.
+      * If you do not use '-detrend', the program attempts to find non-zero spatial
+        structure in the input, and will print a warning message if it is detected.
+
+
+    .. note:: Notes on -demend
+
+      * I recommend this option, and it is not the default only for historical
+        compatibility reasons.  It may become the default someday.
+      * It is already the default in program 3dBlurToFWHM. This is the same detrending
+        as done in 3dDespike; using 2*q+3 basis functions for q > 0.
+      * If you don't use '-detrend', the program now [Aug 2010] checks if a large number
+        of voxels are have significant nonzero means. If so, the program will print a
+        warning message suggesting the use of '-detrend', since inherent spatial
+        structure in the image will bias the estimation of the FWHM of the image time
+        series NOISE (which is usually the point of using 3dFWHMx).
+
+
+    """
+    _cmd = '3dFWHMx'
+    input_spec = FWHMxInputSpec
+    output_spec = FWHMxOutputSpec
+
+    def _parse_inputs(self, skip=None):
+        if not self.inputs.detrend:
+            if skip is None:
+                skip = []
+            skip += ['out_detrend']
+        return super(FWHMx, self)._parse_inputs(skip=skip)
+
+    def _format_arg(self, name, trait_spec, value):
+        if name == 'detrend':
+            if isinstance(value, bool):
+                if value:
+                    return trait_spec.argstr
+                else:
+                    return None
+            elif isinstance(value, int):
+                return trait_spec.argstr + ' %d' % value
+
+        if name == 'acf':
+            if isinstance(value, tuple):
+                return trait_spec.argstr + ' %s %f' % value
+            elif isinstance(value, string_types):
+                return trait_spec.argstr + ' ' + value
+        return super(FWHMx, self)._format_arg(name, trait_spec, value)
+
+    def _run_interface(self, runtime):
+        if platform == 'darwin':
+            # http://afni.nimh.nih.gov/afni/community/board/read.php?1,145346,145347#msg-145347
+            runtime.environ['DYLD_FALLBACK_LIBRARY_PATH'] = '/usr/local/afni/'
+
+        return super(FWHMx, self)._run_interface(runtime)
+
+
+    def _list_outputs(self):
+        outputs = super(FWHMx, self)._list_outputs()
+
+        if self.inputs.detrend:
+            fname, ext = op.splitext(self.inputs.in_file)
+            if '.gz' in ext:
+                _, ext2 = op.splitext(fname)
+                ext = ext2 + ext
+            outputs['out_detrend'] += ext
+        else:
+            outputs['out_detrend'] = Undefined
+
+        outputs['fwhm'] = tuple(np.loadtxt(outputs['out_file']))  #pylint: disable=E1101
+        return outputs

nipype/interfaces/afni/tests/test_auto_FWHMx.py

@@ -0,0 +1,80 @@
+# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT
+from ....testing import assert_equal
+from ..preprocess import FWHMx
+
+
+def test_FWHMx_inputs():
+    input_map = dict(acf=dict(argstr='-acf',
+    ),
+    args=dict(argstr='%s',
+    ),
+    arith=dict(argstr='-arith',
+    xor=['geom'],
+    ),
+    automask=dict(argstr='-automask',
+    usedefault=True,
+    ),
+    combine=dict(argstr='-combine',
+    ),
+    compat=dict(argstr='-compat',
+    ),
+    demed=dict(argstr='-demed',
+    xorg=['detrend'],
+    ),
+    detrend=dict(argstr='-detrend',
+    usedefault=True,
+    xor=['demed'],
+    ),
+    environ=dict(nohash=True,
+    usedefault=True,
+    ),
+    geom=dict(argstr='-geom',
+    xor=['arith'],
+    ),
+    ignore_exception=dict(nohash=True,
+    usedefault=True,
+    ),
+    in_file=dict(argstr='-input %s',
+    mandatory=True,
+    ),
+    mask=dict(argstr='-mask %s',
+    ),
+    out_detrend=dict(argstr='-detprefix %s',
+    keep_extension=False,
+    name_source='in_file',
+    name_template='%s_detrend',
+    ),
+    out_file=dict(argstr='> %s',
+    keep_extension=False,
+    name_source='in_file',
+    name_template='%s_fwhmx.out',
+    position=-1,
+    ),
+    out_subbricks=dict(argstr='-out %s',
+    keep_extension=False,
+    name_source='in_file',
+    name_template='%s_subbricks.out',
+    ),
+    terminal_output=dict(nohash=True,
+    ),
+    unif=dict(argstr='-unif',
+    ),
+    )
+    inputs = FWHMx.input_spec()
+
+    for key, metadata in list(input_map.items()):
+        for metakey, value in list(metadata.items()):
+            yield assert_equal, getattr(inputs.traits()[key], metakey), value
+
+
+def test_FWHMx_outputs():
+    output_map = dict(fwhm=dict(),
+    out_detrend=dict(),
+    out_file=dict(),
+    out_subbricks=dict(),
+    )
+    outputs = FWHMx.output_spec()
+
+    for key, metadata in list(output_map.items()):
+        for metakey, value in list(metadata.items()):
+            yield assert_equal, getattr(outputs.traits()[key], metakey), value