Merge pull request #107 from miykael/joss_review_feedback_2

Second round of updates in reaction to JOSS review

README.md

@@ -44,9 +44,9 @@ python setup.py install
 
 ## Usage
 
-`atlasreader` can either be run through the command line interface or directly
+AtlasReader can either be run through the command line interface or directly
 within Python. The commands to do so are rather straight forward. Let's say you
-want to apply `atlasreader` to a statistical image called
+want to apply AtlasReader to a statistical image called
 `file_name = 'stat_img.nii'`, and only want to keep clusters if they have more
 than 5 voxels:
 
@@ -64,7 +64,7 @@ atlasreader file_name 5
 
 ### Outputs
 
-After executing `atlasreader` on a given image, four kinds of outputs are generated:
+After executing AtlasReader on a given image, four kinds of outputs are generated:
 
 1. An **overview figure** that shows the results within the whole brain at once  
    ![Overview Figure](paper/fig_overview_figure.png)
@@ -90,7 +90,7 @@ After executing `atlasreader` on a given image, four kinds of outputs are genera
 
 ### Additional parameters
 
-`atlasreader` has many additional parameters that allow you to change the way
+`atlasreader.create_output` has many additional parameters that allow you to change the way
 the clusters are generated and what kind of outputs are generated:
 
 - **filename**: Niimg_like  
@@ -99,10 +99,14 @@ the clusters are generated and what kind of outputs are generated:
     Minimum number of contiguous voxels required to consider a cluster in `filename`
 - **atlas**: str or list, optional  
     Name of atlas(es) to consider for cluster analysis. ***Default***: `'default'`
-- **voxel_thresh**: int, optional  
-    Threshold to apply to `stat_img`. If a negative number is provided a
-    percentile threshold is used instead, where the percentile is
-    determined by the equation `100 - voxel_thresh`. ***Default***: `1.96`
+- **voxel_thresh**: float, optional
+    Threshold to apply to `stat_img`.  Use `direction` to specify the 
+    directionality of the threshold. If a negative number is provided a
+    percentile threshold is used instead, where the percentile is determined
+    by the equation `100 - voxel_thresh`. ***Default***: `1.96`
+- **direction**: str, optional
+    Specifies the direction in which `voxel_thresh` should be applied. Possible
+    values are `'both'`, `'pos'` or `'neg'`. ***Default***: `'both'`
 - **prob_thresh**: int, optional  
     Probability (percentage) threshold to apply to `atlas`, if it is
     probabilistic. ***Default***: `5`
@@ -115,8 +119,10 @@ the clusters are generated and what kind of outputs are generated:
     saved to the same folder as `filename`. ***Default***: `None`
 - **glass_plot_kws**: dict or None, optional  
     Additional keyword arguments to pass to `nilearn.plotting.plot_glass_brain`.
+    ***Default***: `None`
 - **stat_plot_kws**: dict or None, optional  
     Additional keyword arguments to pass to `nilearn.plotting.plot_stat_map`.
+    ***Default***: `None`
 
 For a more detailed explanation about the toolbox and the effect of the
 parameters above, see the [example notebook](https://github.com/miykael/atlasreader/blob/master/notebooks/atlasreader.ipynb).
@@ -141,8 +147,8 @@ about it!
 
 ## Licence
 
-`atlasreader` is licensed under the BSD-3 license; however, the atlases it uses 
+AtlasReader is licensed under the BSD-3 license; however, the atlases it uses 
 are separately licensed under more restrictive frameworks.
-By using `atlasreader`, you agree to abide by the license terms of the
+By using AtlasReader, you agree to abide by the license terms of the
 individual atlases. Information on these terms can be found online at:
 https://github.com/miykael/atlasreader/tree/master/atlasreader/data

atlasreader/__init__.py

@@ -4,10 +4,10 @@
 from atlasreader.atlasreader import create_output, get_statmap_info
 
 _LICENSE_MESSAGE = """\
-The Python package you are importing, `atlasreader`, is licensed under the
+The Python package you are importing, AtlasReader, is licensed under the
 BSD-3 license; however, the atlases it uses are separately licensed under more
 restrictive frameworks.
-By using `atlasreader`, you agree to abide by the license terms of the
+By using AtlasReader, you agree to abide by the license terms of the
 individual atlases. Information on these terms can be found online at:
 https://github.com/miykael/atlasreader/tree/master/atlasreader/data
 """

atlasreader/atlasreader.py

@@ -429,7 +429,7 @@ def read_atlas_cluster(atlastype, cluster, affine, prob_thresh=5):
             percentage[s] >= prob_thresh]
 
 
-def process_img(stat_img, cluster_extent, voxel_thresh=1.96):
+def process_img(stat_img, cluster_extent, voxel_thresh=1.96, direction='both'):
     """
     Parameters
     ----------
@@ -438,10 +438,13 @@ def process_img(stat_img, cluster_extent, voxel_thresh=1.96):
     cluster_extent : int
         Minimum number of voxels required to consider a cluster
     voxel_thresh : float, optional
-        Threshold to apply to `stat_img`. The same threshold is applied in both
-        directions, positive and negative. If a negative number is provided a
+        Threshold to apply to `stat_img`. Use `direction` to specify the
+        directionality of the threshold. If a negative number is provided a
         percentile threshold is used instead, where the percentile is
         determined by the equation `100 - voxel_thresh`. Default: 1.96
+    direction : str, optional
+        Specifies the direction in which `voxel_thresh` should be applied.
+        Possible values are 'both', 'pos' or 'neg'. Default: 'both'
 
     Returns
     -------
@@ -468,7 +471,15 @@ def process_img(stat_img, cluster_extent, voxel_thresh=1.96):
     # extract clusters
     min_region_size = cluster_extent * np.prod(thresh_img.header.get_zooms())
     clusters = []
-    for sign in ['pos', 'neg']:
+
+    if direction == 'both':
+        direction_list = ['pos', 'neg']
+    elif direction == 'pos':
+        direction_list = ['pos']
+    elif direction == 'neg':
+        direction_list = ['neg']
+
+    for sign in direction_list:
         # keep only data of given sign
         data = thresh_img.get_data().copy()
         data[(data < 0) if sign == 'pos' else (data > 0)] = 0
@@ -594,7 +605,8 @@ def get_cluster_data(clust_img, atlas='default', prob_thresh=5):
 
 
 def get_statmap_info(stat_img, cluster_extent, atlas='default',
-                     voxel_thresh=1.96, prob_thresh=5, min_distance=None):
+                     voxel_thresh=1.96, direction='both', prob_thresh=5,
+                     min_distance=None):
     """
     Extract peaks and cluster information from `clust_img` for `atlas`
 
@@ -608,10 +620,13 @@ def get_statmap_info(stat_img, cluster_extent, atlas='default',
     atlas : str or list, optional
         Name of atlas(es) to consider for cluster analysis. Default: 'default'
     voxel_thresh : float, optional
-        Threshold to apply to `stat_img`. The same threshold is applied in both
-        directions, positive and negative. If a negative number is provided a
+        Threshold to apply to `stat_img`. Use `direction` to specify the
+        directionality of the threshold. If a negative number is provided a
         percentile threshold is used instead, where the percentile is
         determined by the equation `100 - voxel_thresh`. Default: 1.96
+    direction : str, optional
+        Specifies the direction in which `voxel_thresh` should be applied.
+        Possible values are 'both', 'pos' or 'neg'. Default: 'both'
     prob_thresh : [0, 100] int, optional
         Probability (percentage) threshold to apply to `atlas`, if it is
         probabilistic. Default: 5
@@ -637,6 +652,7 @@ def get_statmap_info(stat_img, cluster_extent, atlas='default',
     # threshold + clusterize image
     clust_img = process_img(stat_img,
                             voxel_thresh=voxel_thresh,
+                            direction=direction,
                             cluster_extent=cluster_extent)
 
     clust_info, peaks_info = [], []
@@ -674,8 +690,8 @@ def get_statmap_info(stat_img, cluster_extent, atlas='default',
 
 
 def create_output(filename, cluster_extent, atlas='default', voxel_thresh=1.96,
-                  prob_thresh=5, min_distance=None, outdir=None,
-                  glass_plot_kws=None, stat_plot_kws=None):
+                  direction='both', prob_thresh=5, min_distance=None,
+                  outdir=None, glass_plot_kws=None, stat_plot_kws=None):
     """
     Performs full cluster / peak analysis on `filename`
 
@@ -697,10 +713,13 @@ def create_output(filename, cluster_extent, atlas='default', voxel_thresh=1.96,
     atlas : str or list, optional
         Name of atlas(es) to consider for cluster analysis. Default: 'default'
     voxel_thresh : float, optional
-        Threshold to apply to `stat_img`. The same threshold is applied in both
-        directions, positive and negative. If a negative number is provided a
+        Threshold to apply to `stat_img`. Use `direction` to specify the
+        directionality of the threshold. If a negative number is provided a
         percentile threshold is used instead, where the percentile is
         determined by the equation `100 - voxel_thresh`. Default: 1.96
+    direction : str, optional
+        Specifies the direction in which `voxel_thresh` should be applied.
+        Possible values are 'both', 'pos' or 'neg'. Default: 'both'
     prob_thresh : int, optional
         Probability (percentage) threshold to apply to `atlas`, if it is
         probabilistic. Default: 5
@@ -713,10 +732,10 @@ def create_output(filename, cluster_extent, atlas='default', voxel_thresh=1.96,
         saved to the same folder as `filename`. Default: None
     glass_plot_kws : dict or None, optional
         Additional keyword arguments to pass to
-        `nilearn.plotting.plot_glass_brain`.
+        `nilearn.plotting.plot_glass_brain`. Default: None
     stat_plot_kws : dict or None, optional
         Additional keyword arguments to pass to
-        `nilearn.plotting.plot_stat_map`.
+        `nilearn.plotting.plot_stat_map`. Default: None
     """
 
     # confirm input data is niimg_like to raise error as early as possible
@@ -744,6 +763,7 @@ def create_output(filename, cluster_extent, atlas='default', voxel_thresh=1.96,
     # generate stat map for plotting by collapsing all clusters into one image
     clust_img = process_img(stat_img,
                             voxel_thresh=voxel_thresh,
+                            direction=direction,
                             cluster_extent=cluster_extent)
     thresh_img = image.math_img('np.sum(img, axis=-1)', img=clust_img)
 
@@ -783,8 +803,8 @@ def create_output(filename, cluster_extent, atlas='default', voxel_thresh=1.96,
         # get cluster + peak information from image
         clust_frame, peaks_frame = get_statmap_info(
             stat_img, atlas=atlas, voxel_thresh=voxel_thresh,
-            cluster_extent=cluster_extent, prob_thresh=prob_thresh,
-            min_distance=min_distance)
+            direction=direction, cluster_extent=cluster_extent,
+            prob_thresh=prob_thresh, min_distance=min_distance)
 
         # write output .csv files
         clust_frame.to_csv(op.join(

atlasreader/cli.py

@@ -50,12 +50,17 @@ def _atlasreader_parser():
                         dest='voxel_thresh', metavar='threshold',
                         help='Value threshold that voxels in provided file '
                              'must surpass in order to be considered in '
-                             'cluster extraction. The same threshold value is '
-                             'applied in both directions, positive and '
-                             'negative. If a negative number is provided a '
-                             'percentile threshold is used instead, where the '
-                             'percentile is determined by the equation '
+                             'cluster extraction. Use `direction` to specify '
+                             'the directionallity of the threshold. If a '
+                             'negative number is provided a  percentile '
+                             'threshold is used instead, where the percentile '
+                             'is determined by the equation '
                              '`100 - voxel_thresh`. Default: 1.96')
+    parser.add_argument('-x', '--direction', type=str, default='both',
+                        dest='direction', metavar='direction',
+                        help='Specifies the direction in which the threshold '
+                             'should be applied. Possible values are '
+                             '\'both\', \'pos\' or \'neg\'. Default: \'both\'')
     parser.add_argument('-p', '--probability', type=_check_limit, default=5,
                         dest='prob_thresh', metavar='threshold',
                         help='Threshold to consider when using a '
@@ -91,6 +96,7 @@ def atlasreader_main():
     create_output(opts.filename,
                   atlas=check_atlases(opts.atlas),
                   voxel_thresh=opts.voxel_thresh,
+                  direction=opts.direction,
                   cluster_extent=opts.cluster_extent,
                   prob_thresh=opts.prob_thresh,
                   outdir=opts.outdir,

atlasreader/tests/test_atlasreader.py

@@ -4,7 +4,7 @@
 import nibabel as nb
 from nilearn.datasets import fetch_neurovault_motor_task
 import pytest
-import hashlib
+import pandas as pd
 
 STAT_IMG = fetch_neurovault_motor_task().images[0]
 EXAMPLE_COORDS = dict(
@@ -41,6 +41,20 @@
     ],
     bounding_shape=np.array([90, 90, 90])
 )
+EXPECTED_TABLES = dict(
+    cluster=np.array([[42,  -25,  58,  6.66003, 36936],
+                      [-36, -25,  55, -6.63604, 15012],
+                      [45,  -19,  16,  5.76538,  7722],
+                      [-15, -52, -26,  6.11673,  7101],
+                      [18,  -55, -23, -5.90086,  5184],
+                      [-36, -19,  19, -5.00687,   648]]),
+    peak=np.array([[42,  -25,  58,  7.94135, 36936],
+                   [-36, -25,  55, -7.94144, 15012],
+                   [45,  -19,  16,  7.94135,  7722],
+                   [-15, -52, -26,  7.94135,  7101],
+                   [18,  -55, -23, -7.94144,  5184],
+                   [-36, -19,  19, -6.21808,   648]])
+)
 
 
 def test_get_atlases():
@@ -105,6 +119,9 @@ def test_process_image():
     # check that defaults for processing image work
     img = atlasreader.process_img(stat_img, cluster_extent=20)
     assert isinstance(img, nb.Nifti1Image)
+    # check that one-sided thresholding works
+    img = atlasreader.process_img(stat_img, direction='neg', cluster_extent=20)
+    assert isinstance(img, nb.Nifti1Image)
     # check that negative voxel threshold works
     img = atlasreader.process_img(stat_img, cluster_extent=20,
                                   voxel_thresh=-10)
@@ -127,6 +144,7 @@ def test_create_output(tmpdir):
     # temporary output
     output_dir = tmpdir.mkdir('mni_test')
     atlasreader.create_output(STAT_IMG, cluster_extent=20,
+                              voxel_thresh=7,
                               atlas=['Harvard_Oxford'],
                               outdir=output_dir)
 
@@ -146,13 +164,15 @@ def test_plotting(tmpdir):
 
     # overwrite some default params
     atlasreader.create_output(STAT_IMG, cluster_extent=20,
+                              voxel_thresh=7,
                               atlas=['Harvard_Oxford'],
                               outdir=output_dir,
                               glass_plot_kws={'display_mode': 'ortho'},
                               stat_plot_kws={'black_bg': False})
 
     # add new parameter not already set by default
     atlasreader.create_output(STAT_IMG, cluster_extent=20,
+                              voxel_thresh=7,
                               atlas=['Harvard_Oxford'],
                               outdir=output_dir,
                               glass_plot_kws={'alpha': .4})
@@ -165,16 +185,11 @@ def test_table_output(tmpdir):
     # temporary output
     output_dir = tmpdir.mkdir('mni_test')
     atlasreader.create_output(STAT_IMG, cluster_extent=20,
-                              atlas=['AAL', 'desikan_killiany',
-                                     'Harvard_Oxford'],
+                              voxel_thresh=4, atlas='default',
                               outdir=output_dir)
 
     # test if output tables contain expected output
-    cluster_checksum = hashlib.md5(open(output_dir.join(
-        '{}_clusters.csv'.format(stat_img_name)), 'rb').read()).hexdigest()
-
-    peak_checksum = hashlib.md5(open(output_dir.join(
-        '{}_peaks.csv'.format(stat_img_name)), 'rb').read()).hexdigest()
-
-    assert cluster_checksum == '5d85805d58f8fbe22ef4e34ec75d8f42'
-    assert peak_checksum == 'f7cd664571413fe964eef0c45cd6f033'
+    df = pd.read_csv(output_dir.join('{}_clusters.csv'.format(stat_img_name)))
+    assert np.allclose(df[df.keys()[1:6]].values, EXPECTED_TABLES['cluster'])
+    df = pd.read_csv(output_dir.join('{}_peaks.csv'.format(stat_img_name)))
+    assert np.allclose(df[df.keys()[1:6]].values, EXPECTED_TABLES['peak'])

atlasreader/tests/test_cli.py

@@ -12,7 +12,7 @@ def test_cli(tmpdir):
     output_dir = tmpdir.mkdir('mni_test')
     ret = subprocess.run(['atlasreader',
                           '--atlas', 'harvard_oxford', 'aal',
-                          '--threshold', '1.96',
+                          '--threshold', '7.0',
                           '--probability', '5',
                           '--mindist', '20',
                           '--outdir', output_dir,