From 4eca2c3f1792c98ae0b5e389c405654d4b3771d5 Mon Sep 17 00:00:00 2001
From: Matthew Brett <matthew.brett@gmail.com>
Date: Thu, 13 Feb 2014 22:27:32 -0800
Subject: [PATCH] NF: tsdiffana option to save results volumes

Add option for tsdiffana to be able to save the results volumes created
during the analysis.
---
 nipy/algorithms/diagnostics/commands.py       | 45 +++++++++--
 .../diagnostics/tests/test_commands.py        | 60 ++++++++++++--
 nipy/tests/test_scripts.py                    | 29 +++++--
 scripts/nipy_tsdiffana                        | 79 +++++++++++++++++--
 4 files changed, 189 insertions(+), 24 deletions(-)

diff --git a/nipy/algorithms/diagnostics/commands.py b/nipy/algorithms/diagnostics/commands.py
index da6a28709f..8bbc906fcb 100644
--- a/nipy/algorithms/diagnostics/commands.py
+++ b/nipy/algorithms/diagnostics/commands.py
@@ -11,13 +11,17 @@
 The implementation here accepts the ``args`` object from ``argparse`` and does
 the work.
 """
+from os.path import split as psplit, join as pjoin
+
+import numpy as np
 
 from nibabel import AnalyzeHeader
+from nibabel.filename_parser import splitext_addext
 
 import nipy
 
 from .tsdiffplot import plot_tsdiffs
-from .timediff import time_slice_diffs, time_slice_diffs_image
+from .timediff import time_slice_diffs_image
 
 
 def parse_fname_axes(img_fname, time_axis, slice_axis):
@@ -90,24 +94,51 @@ def tsdiffana(args):
     Parameters
     ----------
     args : object
-        object with (None or ``str``) attributes:
+        object with attributes
 
-        * filename : 4D image filename
-        * out_file: graphics file to write to instead of leaving graphics on
-          screen
-        * time_axis : name or number of time axis in `filename`
-        * slice_axis : name or number of slice axis in `filename`
+        * filename : str - 4D image filename
+        * out_file : str - graphics file to write to instead of leaving
+          graphics on screen
+        * time_axis : str - name or number of time axis in `filename`
+        * slice_axis : str - name or number of slice axis in `filename`
+        * write_results : bool - if True, write images and plots to files
+        * out_path : None or str - path to which to write results
+        * out_fname_label : None or filename - suffix of output results files
 
     Returns
     -------
     axes : Matplotlib axes
        Axes on which we have done the plots.
     """
+    if not args.out_file is None and args.write_results:
+        raise ValueError("Cannot have OUT_FILE and WRITE_RESULTS options "
+                         "together")
     img, time_axis, slice_axis = parse_fname_axes(args.filename,
                                                   args.time_axis,
                                                   args.slice_axis)
     results = time_slice_diffs_image(img, time_axis, slice_axis)
     axes = plot_tsdiffs(results)
+    if args.out_file is None and not args.write_results:
+        # interactive mode
+        return axes
     if not args.out_file is None:
+        # plot only mode
         axes[0].figure.savefig(args.out_file)
+        return axes
+    # plot and images mode
+    froot, ext, addext = splitext_addext(args.filename)
+    fpath, fbase = psplit(froot)
+    fpath = fpath if args.out_path is None else args.out_path
+    fbase = fbase if args.out_fname_label is None else args.out_fname_label
+    axes[0].figure.savefig(pjoin(fpath, 'tsdiff_' + fbase + '.png'))
+    # Save image volumes
+    for key, prefix in (('slice_diff2_max_vol', 'dv2_max_'),
+                        ('diff2_mean_vol', 'dv2_mean_')):
+        fname = pjoin(fpath, prefix + fbase + ext + addext)
+        nipy.save_image(results[key], fname)
+    # Save time courses into npz
+    np.savez(pjoin(fpath, 'tsdiff_' + fbase + '.npz'),
+             volume_means=results['volume_means'],
+             slice_mean_diff2=results['slice_mean_diff2'],
+            )
     return axes
diff --git a/nipy/algorithms/diagnostics/tests/test_commands.py b/nipy/algorithms/diagnostics/tests/test_commands.py
index 415e62fa76..4babed82d7 100644
--- a/nipy/algorithms/diagnostics/tests/test_commands.py
+++ b/nipy/algorithms/diagnostics/tests/test_commands.py
@@ -1,7 +1,9 @@
 """ Testing diagnostics.command module
 """
 
+import os
 from os.path import dirname, join as pjoin, isfile
+import shutil
 
 import numpy as np
 
@@ -11,9 +13,14 @@
 
 from nipy import load_image
 from ..commands import parse_fname_axes, tsdiffana
+from ..timediff import time_slice_diffs_image
 
-from numpy.testing import (assert_almost_equal,
-                           assert_array_equal)
+from numpy.testing import (assert_almost_equal, assert_array_equal, decorators)
+
+from nibabel.optpkg import optional_package
+
+matplotlib, HAVE_MPL, _ = optional_package('matplotlib')
+needs_mpl = decorators.skipif(not HAVE_MPL, "Test needs matplotlib")
 
 from nose import SkipTest
 from nose.tools import (assert_true, assert_false, assert_raises,
@@ -102,20 +109,24 @@ def check_axes(axes, img_shape, time_axis, slice_axis):
     assert_equal(len(axes), 4)
     # First x axis is time point differences
     assert_array_equal(axes[0].xaxis.get_data_interval(),
-                        [0, img_shape[time_axis]-2])
+                       [0, img_shape[time_axis]-2])
     # Last x axis is over slices
     assert_array_equal(axes[-1].xaxis.get_data_interval(),
-                        [0, img_shape[slice_axis]-1])
+                       [0, img_shape[slice_axis]-1])
 
 
+@needs_mpl
 def test_tsdiffana():
     # Test tsdiffana command
     args = Args()
     img = load_image(funcfile)
-    with InTemporaryDirectory():
+    with InTemporaryDirectory() as tmpdir:
         args.filename = funcfile
         args.time_axis = None
         args.slice_axis = None
+        args.write_results = False
+        args.out_path = None
+        args.out_fname_label = None
         args.out_file = 'test.png'
         check_axes(tsdiffana(args), img.shape, -1, -2)
         assert_true(isfile('test.png'))
@@ -131,3 +142,42 @@ def test_tsdiffana():
         check_axes(tsdiffana(args), img.shape, 0, -2)
         args.slice_axis = 't'
         check_axes(tsdiffana(args), img.shape, 0, -1)
+        # Check absolute path works
+        args.slice_axis = 'j'
+        args.time_axis = 't'
+        args.out_file = pjoin(tmpdir, 'test_again.png')
+        check_axes(tsdiffana(args), img.shape, -1, -3)
+        # Check that --out-images incompatible with --out-file
+        args.write_results=True
+        assert_raises(ValueError, tsdiffana, args)
+        args.out_file=None
+        # Copy the functional file to a temporary writeable directory
+        os.mkdir('mydata')
+        tmp_funcfile = pjoin(tmpdir, 'mydata', 'myfunc.nii.gz')
+        shutil.copy(funcfile, tmp_funcfile)
+        args.filename = tmp_funcfile
+        # Check write-results generates expected images
+        check_axes(tsdiffana(args), img.shape, -1, -3)
+        assert_true(isfile(pjoin('mydata', 'tsdiff_myfunc.png')))
+        max_img = load_image(pjoin('mydata', 'dv2_max_myfunc.nii.gz'))
+        assert_equal(max_img.shape, img.shape[:-1])
+        mean_img = load_image(pjoin('mydata', 'dv2_max_myfunc.nii.gz'))
+        assert_equal(mean_img.shape, img.shape[:-1])
+        exp_results = time_slice_diffs_image(img, 't', 'j')
+        saved_results = np.load(pjoin('mydata', 'tsdiff_myfunc.npz'))
+        for key in ('volume_means', 'slice_mean_diff2'):
+            assert_array_equal(exp_results[key], saved_results[key])
+        # That we can change out-path
+        os.mkdir('myresults')
+        args.out_path = 'myresults'
+        check_axes(tsdiffana(args), img.shape, -1, -3)
+        assert_true(isfile(pjoin('myresults', 'tsdiff_myfunc.png')))
+        max_img = load_image(pjoin('myresults', 'dv2_max_myfunc.nii.gz'))
+        assert_equal(max_img.shape, img.shape[:-1])
+        # And out-fname-label
+        args.out_fname_label = 'vr2'
+        check_axes(tsdiffana(args), img.shape, -1, -3)
+        assert_true(isfile(pjoin('myresults', 'tsdiff_vr2.png')))
+        max_img = load_image(pjoin('myresults', 'dv2_max_vr2.nii.gz'))
+        assert_equal(max_img.shape, img.shape[:-1])
+        del max_img, mean_img
diff --git a/nipy/tests/test_scripts.py b/nipy/tests/test_scripts.py
index 193179ea8c..7b83e5d9b6 100644
--- a/nipy/tests/test_scripts.py
+++ b/nipy/tests/test_scripts.py
@@ -21,7 +21,7 @@
 from nipy import load_image, save_image
 from nipy.core.api import rollimg
 
-from nose.tools import assert_true, assert_false, assert_equal
+from nose.tools import assert_true, assert_false, assert_equal, assert_raises
 
 from ..testing import funcfile
 from numpy.testing import decorators, assert_almost_equal
@@ -123,13 +123,32 @@ def test_nipy_tsdiffana():
     cmd_template = 'nipy_tsdiffana "{0}" --out-file="{1}"'
     with InTemporaryDirectory():
         for i, cmd in enumerate((cmd_template,
-                                 cmd_template + ' --time=axis=0',
-                                 cmd_template + ' --slice=axis=0',
-                                 cmd_template + ' --slice=axis=0 ' +
+                                 cmd_template + ' --time-axis=0',
+                                 cmd_template + ' --slice-axis=0',
+                                 cmd_template + ' --slice-axis=0 ' +
                                  '--time-axis=1')):
             out_png = 'ts_out{0}.png'.format(i)
-            run_command(cmd_template.format(funcfile, out_png))
+            run_command(cmd.format(funcfile, out_png))
             assert_true(isfile(out_png))
+    # Out-file and write-results incompatible
+    assert_raises(RuntimeError,
+                  run_command,
+                  cmd_template.format(funcfile, out_png) + '--write-results')
+    # Can save images
+    cmd_root = 'nipy_tsdiffana "{0}" '.format(funcfile)
+    with InTemporaryDirectory():
+        os.mkdir('myresults')
+        run_command(cmd_root + '--out-path=myresults --write-results')
+        assert_true(isfile(pjoin('myresults', 'tsdiff_functional.png')))
+        assert_true(isfile(pjoin('myresults', 'tsdiff_functional.npz')))
+        assert_true(isfile(pjoin('myresults', 'dv2_max_functional.nii.gz')))
+        assert_true(isfile(pjoin('myresults', 'dv2_mean_functional.nii.gz')))
+        run_command(cmd_root + '--out-path=myresults --write-results '
+                    '--out-fname-label=vr2')
+        assert_true(isfile(pjoin('myresults', 'tsdiff_vr2.png')))
+        assert_true(isfile(pjoin('myresults', 'tsdiff_vr2.npz')))
+        assert_true(isfile(pjoin('myresults', 'dv2_max_vr2.nii.gz')))
+        assert_true(isfile(pjoin('myresults', 'dv2_mean_vr2.nii.gz')))
 
 
 @script_test
diff --git a/scripts/nipy_tsdiffana b/scripts/nipy_tsdiffana
index 4fce3a9394..189734427f 100755
--- a/scripts/nipy_tsdiffana
+++ b/scripts/nipy_tsdiffana
@@ -1,27 +1,92 @@
 #!/usr/bin/env python
 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
 # vi: set ft=python sts=4 ts=4 sw=4 et:
-''' Analyze, plot time series difference metrics'''
+DESCRIP = ' Analyze, plot time series difference metrics'
+EPILOG = \
+'''nipy_tsdiffana runs the time series diference algorithm over a 4D
+image volume, often and FMRI volume.
 
-import nipy.externals.argparse as argparse
+It works in one of three modes:
+
+* interactive : the time series difference plot appears on screen.  This is the
+  default mode
+* non-interactive, plot only : write time series difference plot to graphic
+  file. Use the "--out-file=<myfilename>" option to activate this mode
+* non-interactive, write plot, images and variables : write plot to file, and
+  write generated diagnostic images and variables to files as well. Use the
+  "--write-results" flag to activate this option. The generated filenames come
+  from the results of the "--out-path" and "--out-fname-label" options (see
+  help).
+
+Write-results option, generated files
+-------------------------------------
+
+When doing the time point analysis, we will make a difference volume between
+each time point and the next time point in the series.  If we have T volumes
+then there will be (T-1) difference volumes. Call the vector of difference
+volumes DV and the first difference volume DV[0].  So DV[0] results from
+subtraction of the second volume in the 4D input image from the first volume in
+the 4D input image.  The element-wise squared values from DV[0] is *DV2[0]*.
+
+The following images will be generated. <ext> is the input filename extension
+(e.g. '.nii'):
+
+* "dv2_max_<label><ext>" : 3D image volume, where each slice S is slice from
+  all of DV2[0] (slice S) throudh DV2[T-1] (slice S) that has the maximum
+  summed squared values.  This volume gives an idea of the worst (highest
+  difference) slices across the whole time series.
+* "dv2_mean_<label><ext>" : the mean of all DV2 volumes DV2[0] .. DV[T-1]
+  across the volume (time) dimension.  Higher voxel values in this volume mean
+  that time-point to time point differences tended to be high in this voxel.
+
+We also write the mean signal at each time point, and the mean squared
+difference between each slice in time, as variables to a 'npz' file named
+"tsdiff_<label>.npz"
+
+The filenames for the outputs are of the form
+<out-path>/<some_prefix><label><file-ext> where <out-path> is the path
+specified by the --out-path option, or the path of the input filename;
+<some_prefix> is one of the standard prefixes above, <label> is given by
+--out-label, or by the filename of the input image (with path and extension
+removed), and <file-ext> is '.png' for graphics, or the extension of the input
+filename for volume images.  For example, specifying only the input filename
+``/some/path/fname.img`` will generate filenames of the form
+``/some/path/tsdiff_fname.png, /some/path/dv2_max_fname.img`` etc.
+'''
+
+from nipy.externals.argparse import (ArgumentParser,
+                                     RawDescriptionHelpFormatter)
 
 
 def main():
-    # create the parser
-    parser = argparse.ArgumentParser()
+    parser = ArgumentParser(description=DESCRIP,
+                            epilog=EPILOG,
+                            formatter_class=RawDescriptionHelpFormatter)
     # add the arguments
+    parser.add_argument('filename', type=str,
+                        help='4D image filename')
     parser.add_argument('--out-file', type=str,
                         help='graphics file to write to instead '
                         'of leaving image on screen')
-    parser.add_argument('filename', type=str,
-                        help='4D image filename')
+    parser.add_argument('--write-results', action='store_true',
+                        help='if specified, write diagnostic images and '
+                        'analysis variables, plot to OUT_PATH. Mutually '
+                        'incompatible with OUT_FILE')
+    parser.add_argument('--out-path', type=str,
+                        help='path for output image files (default from '
+                        'FILENAME path')
+    parser.add_argument('--out-fname-label', type=str,
+                        help='mid part of output image / plot filenames')
     parser.add_argument('--time-axis', type=str, default='t',
                         help='Image axis for time')
     parser.add_argument('--slice-axis', type=str, default=None,
                         help='Image axis for slice')
     # parse the command line
     args = parser.parse_args()
-    show = args.out_file is None
+    if not args.out_file is None and args.write_results == True:
+        raise RuntimeError("OUT_FILE option not compatible with WRITE_RESULTS"
+                           " option")
+    show = args.out_file is None and args.write_results == False
     if not show:
         import matplotlib
         matplotlib.use('Agg')