ENH: Add an activation count map interface

nipype/algorithms/stats.py

@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""
+Managing statistical maps
+"""
+from __future__ import (print_function, division, unicode_literals,
+                        absolute_import)
+import os
+import nibabel as nb
+import numpy as np
+
+from ..interfaces.base import (
+    BaseInterfaceInputSpec, TraitedSpec, SimpleInterface,
+    traits, InputMultiPath, File
+)
+from ..utils.filemanip import split_filename
+
+
+class ActivationCountInputSpec(BaseInterfaceInputSpec):
+    in_files = InputMultiPath(File(exists=True), mandatory=True,
+                              desc='input file, generally a list of z-stat maps')
+    threshold = traits.Float(
+        mandatory=True, desc='binarization threshold. E.g. a threshold of 1.65 '
+                             'corresponds to a two-sided Z-test of p<.10')
+
+
+class ActivationCountOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc='output activation count map')
+    acm_pos = File(exists=True, desc='positive activation count map')
+    acm_neg = File(exists=True, desc='negative activation count map')
+
+
+class ActivationCount(SimpleInterface):
+    """
+    Calculate a simple Activation Count Maps
+
+    Adapted from: https://github.com/poldracklab/CNP_task_analysis/\
+    blob/61c27f5992db9d8800884f8ffceb73e6957db8af/CNP_2nd_level_ACM.py
+    """
+    input_spec = ActivationCountInputSpec
+    output_spec = ActivationCountOutputSpec
+
+    def _run_interface(self, runtime):
+        allmaps = nb.concat_images(self.inputs.in_files).get_data()
+        acm_pos = np.mean(allmaps > self.inputs.threshold,
+                          axis=3, dtype=np.float32)
+        acm_neg = np.mean(allmaps < -1.0 * self.inputs.threshold,
+                          axis=3, dtype=np.float32)
+        acm_diff = acm_pos - acm_neg
+
+        template_fname = self.inputs.in_files[0]
+        ext = split_filename(template_fname)[2]
+        fname_fmt = os.path.join(runtime.cwd, 'acm_{}' + ext).format
+
+        self._results['out_file'] = fname_fmt('diff')
+        self._results['acm_pos'] = fname_fmt('pos')
+        self._results['acm_neg'] = fname_fmt('neg')
+
+        img = nb.load(template_fname)
+        img.__class__(acm_diff, img.affine, img.header).to_filename(
+            self._results['out_file'])
+        img.__class__(acm_pos, img.affine, img.header).to_filename(
+            self._results['acm_pos'])
+        img.__class__(acm_neg, img.affine, img.header).to_filename(
+            self._results['acm_neg'])
+
+        return runtime

nipype/algorithms/tests/test_auto_ActivationCount.py

@@ -0,0 +1,31 @@
+# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT
+from __future__ import unicode_literals
+from ..stats import ActivationCount
+
+
+def test_ActivationCount_inputs():
+    input_map = dict(
+        ignore_exception=dict(
+            deprecated='1.0.0',
+            nohash=True,
+            usedefault=True,
+        ),
+        in_files=dict(mandatory=True, ),
+        threshold=dict(mandatory=True, ),
+    )
+    inputs = ActivationCount.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_ActivationCount_outputs():
+    output_map = dict(
+        acm_neg=dict(),
+        acm_pos=dict(),
+        out_file=dict(),
+    )
+    outputs = ActivationCount.output_spec()
+
+    for key, metadata in list(output_map.items()):
+        for metakey, value in list(metadata.items()):
+            assert getattr(outputs.traits()[key], metakey) == value

nipype/algorithms/tests/test_stats.py

@@ -0,0 +1,45 @@
+# -*- coding: utf-8 -*-
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+
+import numpy as np
+import nibabel as nb
+from nipype.algorithms.stats import ActivationCount
+import pytest
+
+
+def test_ActivationCount(tmpdir):
+    tmpdir.chdir()
+    in_files = ['{:d}.nii'.format(i) for i in range(3)]
+    for fname in in_files:
+        nb.Nifti1Image(np.random.normal(size=(5, 5, 5)),
+                       np.eye(4)).to_filename(fname)
+
+    acm = ActivationCount(in_files=in_files, threshold=1.65)
+    res = acm.run()
+    diff = nb.load(res.outputs.out_file)
+    pos = nb.load(res.outputs.acm_pos)
+    neg = nb.load(res.outputs.acm_neg)
+    assert np.allclose(diff.get_data(), pos.get_data() - neg.get_data())
+
+
+@pytest.mark.parametrize("threshold, above_thresh", [
+    (1, 15.865),  # above one standard deviation (one side)
+    (2, 2.275),   # above two standard deviations (one side)
+    (3, 0.135)    # above three standard deviations (one side)
+])
+def test_ActivationCount_normaldistr(tmpdir, threshold, above_thresh):
+    tmpdir.chdir()
+    in_files = ['{:d}.nii'.format(i) for i in range(3)]
+    for fname in in_files:
+        nb.Nifti1Image(np.random.normal(size=(100, 100, 100)),
+                       np.eye(4)).to_filename(fname)
+
+    acm = ActivationCount(in_files=in_files, threshold=threshold)
+    res = acm.run()
+    pos = nb.load(res.outputs.acm_pos)
+    neg = nb.load(res.outputs.acm_neg)
+    assert np.isclose(pos.get_data().mean(),
+                      above_thresh * 1.e-2, rtol=0.1, atol=1.e-4)
+    assert np.isclose(neg.get_data().mean(),
+                      above_thresh * 1.e-2, rtol=0.1, atol=1.e-4)