#!/usr/bin/env python
# coding: utf-8

# # AtlasReader
# 
# [`atlasreader`](https://github.com/miykael/atlasreader) is a Python interface that generates coordinate tables and region labels from statistical MRI images.
# 
# This notebook will show you how to install AtlasReader, how to run the toolbox and what the different parameters allow you to do.
# 
# 
# ## Installing AtlasReader
# 
# Provided you have `pip` at your disposal, installing AtlasReader is as simple as this:
# 
# ```bash
# pip install atlasreader
# ```
# 
# If you want to install AtlasReader directly from the source code, use the following code:
# 
# ```bash
# git clone https://github.com/miykael/atlasreader.git
# cd atlasreader
# python setup.py install
# ```
# 
# ## Using AtlasReader
# 
# 
# ### License terms of AtlasReader
# 
# Before showcasing the toolbox, we want to point to the particular **license issue** that you need to consider when you are using AtlasReader.
# 
# AtlasReader itself is licensed under the BSD-3 license, but the included atlases that are used to lookup the coordinate labels are under their own license terms. For this reason, we ask you to inform yourself about the atlas specific license terms under [https://github.com/miykael/atlasreader/tree/master/atlasreader/data](https://github.com/miykael/atlasreader/tree/master/atlasreader/data) and to cite the corresponding papers should you use AtlasReader in a publication.
# 
# We've included the following information message after every first module import, to make sure that users are aware about this issue:
# 
#     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
#     individual atlases. Information on these terms can be found online at:
#     https://github.com/miykael/atlasreader/tree/master/atlasreader/data
# 
# 
# ### Citing atlasreader
# 
# If you use AtlasReader in your publication, please cite the following paper:
# 
#     Notter M. P., Gale D., Herholz P., Markello R. D., Notter-Bielser M.-L., & Whitaker K. (2019). AtlasReader: A Python package to generate coordinate tables, region labels, and informative figures from statistical MRI images. *Journal of Open Source Software, 4(34), 1257*, [https://doi.org/10.21105/joss.01257](https://doi.org/10.21105/joss.01257)
# 
# 

# ### Prepare an example file
# 
# AtlasReader can either be run through the command line interface, or directly within Python. To showcase both of those examples, let's first import relevant plotting functions for this notebook.

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
from nilearn.plotting import plot_glass_brain
from IPython.display import Image


# Next, let's grab an example stat map. We'll use `nilearn` to grab a motor task stat map from [neurovault](https://neurovault.org/) and to threshold it at a value of 3 to create some clusters. Note that you'll need the most up to date version of `nilearn` in order to run `fetch_neurovault_motor_task()`.

# In[2]:


# Get a stat map from neurovault using nilearn
from nilearn.datasets import fetch_neurovault_motor_task

motor_images = fetch_neurovault_motor_task()
stat_img = motor_images.images[0]


# In[3]:


# Threshold image at a value of 3 to create clusters
from nilearn.image import threshold_img

stat_img = threshold_img(stat_img, threshold=3)


# In[4]:


# Save the thresholded image into a NIfTI file
file_name = "stat_img.nii.gz"
stat_img.to_filename(file_name)


# And what does our example data look like?

# In[6]:


plot_glass_brain(
    file_name,
    black_bg=True,
    colorbar=True,
    plot_abs=False,
    display_mode="lyrz",
    title="Finger tapping task",
);


# Ok, so we can see that there are a few particular clusters, if we threshold the data at a value of 3. But how big are they and in which regions are they?
# 
# To answer these questions, you can use AtlasReader. It gives you the opportunity to better understand where the peaks of those clusters are, over which regions the clusters extend, and much more. So, how is it done?

# ## Calling AtlasReader from Python
# 
# If you want to run AtlasReader directly in Python, just import the `create_output()` function, and you're good to go:

# In[7]:


from atlasreader import create_output


# Great, so now we can run AtlasReader. The only thing that we have to provide is our **stat map** (either the `file_name` or the image object) and the **cluster extent**, which specifies the minimum number of contiguous voxels required for a cluster to be considered for the analysis. For this example, let's set the `cluster_extent` to 5.

# In[8]:


create_output(file_name, cluster_extent=5, direction="both")


# After execution, we get four different kind of files:
# 
# - An **overview figure** that shows the results within the whole brain at once
# - For **each cluster**, an **informative figure** showing the sagittal, coronal and trasnveral plane centered at this center of the cluster
# - A **csv file** containing relevant information about the **peak** of each cluster
# - A **csv file** containing relevant information about each **cluster**
# 
# So, let's take a closer look at those outputs:

# ### Overview figure
# 
# The overview figure shows all clusters through the brain, plotted on a glass brain plot. The name of file is the name of provided statistical image, but with the file ending `.png`

# In[9]:


Image("stat_img.png")


# ### Informative cluster figure
# 
# AtlasReader creates for each cluster an informative figure, showing the cluster center in the three anatomical planes: sagittal, coronal and trasnveral. The name of those files all end in `_clusterXY.png`, i.e. `_cluster01.png`, `_cluster02.png`, ...:

# In[10]:


Image("stat_img_cluster01.png")


# In[11]:


Image("stat_img_cluster03.png")


# ### Peak table
# 
# The csv file ending with `_peaks.csv` contains the location of each cluster's peak, it's signal value at this location, the cluster extend (in mm, not in number of voxels) and the atlas correspondence of the peak:

# In[12]:


import pandas as pd

pd.read_csv("stat_img_peaks.csv")


# ### Cluster table
# 
# The csv file ending with `_clusters.csv` contains the location of each cluster's peak, the mean value within the cluster, the cluster extend (in mm, not in number of voxels), as well as the membership of each cluster, given a particular atlas.

# In[13]:


pd.read_csv("stat_img_clusters.csv")


# For example, in the csv table shown above, we know that 60.17% of the second cluster is in the left postcentral and 26.32% is in the left precentral area, according to the AAL atlas.

# ## Additional parameters
# 
# As mentioned before, `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. So, let's have a closer look:
# 
# - **filename**: Niimg_like  
#     A 3D statistical image.
# - **cluster_extent**: int  
#     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**: 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`
# - **min_distance**: float, optional  
#     Specifies the minimum distance (in mm) required between sub-peaks in a
#     cluster. If None, sub-peaks will not be examined and only the primary
#     cluster peak will be reported. ***Default***: `None`
# - **outdir**: str or None, optional  
#     Path to desired output directory. If None, generated files will be
#     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`.
# - **stat_plot_kws**: dict or None, optional  
#     Additional keyword arguments to pass to `nilearn.plotting.plot_stat_map`.
#     
# The first two, we've already used in the previous example, but let's take a closer look at the others.

# ### `atlas` and `voxel_thresh`
# 
# First, let's change the **atlas** and the **voxel_thresh** parameters:

# In[14]:


create_output(
    file_name, cluster_extent=5, atlas=["destrieux", "marsatlas"], voxel_thresh=6
)


# As you can see, the usage of a higher voxel threshold created smaller clusters and caused the clusters with below threshold voxels to disappear.

# In[15]:


Image("stat_img.png")


# Additionally, the csv files now contain the information about the `destrieux` and `marsatlas`.

# In[16]:


pd.read_csv("stat_img_clusters.csv")


# **Note**: The three atlases from the previous example, `aal`, `desikan_killiany`, `harvard_oxford` are the **default** atlases, which can either be specified with `atlas='default'`, `atlas=['aal', 'desikan_killiany', 'harvard_oxford']` or by just omitting the `atlas` parameter.

# If a negative number is provided for `voxel_thresh`, a percentile threshold is used instead, where the percentile is determined by the equation `100 - voxel_thresh`. So, let's threshold the data at 99%.

# In[17]:


create_output(file_name, cluster_extent=5, voxel_thresh=-1)


# In[18]:


Image("stat_img.png")


# **Note:** This percentage term is misleading if the statistical images contains negative values. For example, 99% in the above example, is the 99 percentile of all values between -7.9 to 7.9.

# In[19]:


pd.read_csv("stat_img_clusters.csv")


# ### `direction`
# 
# The **direction** parameter allows you to specify the directionality of the thresholding. By default, the directionality is `'both'`, which means that the `voxel_thresh` value is applied to positive and negative values. Changing `direction` to `'pos'` or `'neg'` allows you to apply the threshold only to one direction and ignore everything else in the other direction.

# In[20]:


create_output(
    file_name, cluster_extent=5, atlas=["aal"], voxel_thresh=4, direction="both"
)


# In[21]:


Image("stat_img.png")


# In[22]:


create_output(file_name, cluster_extent=5, atlas=["aal"], voxel_thresh=4, direction="pos")


# In[23]:


Image("stat_img.png")


# In[24]:


create_output(file_name, cluster_extent=5, atlas=["aal"], voxel_thresh=4, direction="neg")


# In[25]:


Image("stat_img.png")


# This thresholding directionality is of course also preserved in the tables.

# In[26]:


pd.read_csv("stat_img_clusters.csv")


# ### `prob_thresh`
# 
# The **prob_thresh** parameter can be used to restrict the information in the **cluster table** to a given probability value. This will only influence the information in the table, but won't change the cluster size or any other outcome.
# 
# So let's take again an example from before and set the probability threshold to 33%

# In[27]:


create_output(file_name, cluster_extent=5, prob_thresh=33.0)


# As you can see, the **cluster** table has now fewer entries in the atlas column, as atlas corresponances below 33% where omitted. Be careful, if this threshold is set to high, no atlas corresponance will survive. In which case the entry is replaces by `NaN`.

# In[28]:


pd.read_csv("stat_img_clusters.csv")


# ### `min_distance`
# 
# By default, AtlasReader only shows information of the main peak in the **peak table**. But some clusters might be rather big with many different peaks. If you also want to extract peak information of other peaks within a given cluster, you can use the `min_distance` parameter. This parameter specifies the minimum distance required between sub-peaks in a cluster.
# 
# To better illustrate this example, let's focus only on the big four cluster in our dataset (i.e. we will set `cluster_extent=100`).

# In[29]:


create_output(file_name, cluster_extent=100, atlas="aal", min_distance=10)


# In[30]:


Image("stat_img.png")


# Setting the minimum distance between peaks to 10, means that the **peak table** will now include multiple peaks withini a cluster, if they have a euclidean distance of at least 10mm.

# In[31]:


pd.read_csv("stat_img_peaks.csv")


# **Note**: The usage of the `min_distance` parameter has only an effect on the **peak table**, and will not influence any other output file.

# ### `outdir`
# 
# The **outdir** parameter allows you to store the output files of AtlasReader at a particular location. If none is specified, then AtlasReader will store the results at the same location, as the input image.
# 
# As an example, let's store the output at `/tmp/example`:

# In[32]:


create_output(file_name, cluster_extent=100, outdir="/tmp/example")


# In[33]:


ls /tmp/example/


# ### `glass_plot_kws` and `stat_plot_kws`
# 
# The last two parameters allow you to adjust the style of the output plots. Let's say that we want to tweak the default plotting style to have a white background, or adjust the colormap used to convey the clusters. This customization is done separately for the glass brain plot and the cluster plots using the `_kws` arguments in `create_output`. These parameters are dictionaries that contain keyword arguments for the `nilearn` functions used to create these plots (see the keyword arguments for [plot_glass_brain](http://nilearn.github.io/modules/generated/nilearn.plotting.plot_glass_brain.html#nilearn.plotting.plot_glass_brain) and [plot_stat_map](http://nilearn.github.io/modules/generated/nilearn.plotting.plot_stat_map.html#nilearn.plotting.plot_stat_map)).
# 
# For example:

# In[34]:


create_output(
    file_name,
    cluster_extent=100,
    glass_plot_kws={"black_bg": False, "vmax": 20, "colorbar": False},
    stat_plot_kws={"black_bg": False, "title": None},
)


# This sets the background color to white, sets the maximum value of the colorbar to 20 and removes the colorbar from the overview plot.

# In[35]:


Image("stat_img.png")


# The cluster specific plots also have now a white background and the image title is removed from the figure.

# In[36]:


# Title is now removed for all cluster plots
Image("stat_img_cluster01.png")


# ## Calling AtlasReader from the Command Line
# 
# This is all super exciting! But as promised before, AtlasReader can also be run directly from the command line. Assuming you installed it via `pip`.
# 
# Let's begin with running the `help` argument so we can look at the input parameters.

# In[37]:


get_ipython().run_cell_magic('bash', '', 'atlasreader -h \n')


# **Perfect!** As you can see, first, we have again the information about the license terms, followed by the description of the software. The parameters that you can specify are the same input parameters as if we would run AtlasReader within Python.
# 
# So let's try it out. Let's say we want to set the **probability threshold** to `20%` to remove small probabilities, use `all` **atlases** to extract information from, **threshold** the image at `3` and only keep clusters with more than `10` voxels:

# In[38]:


get_ipython().run_cell_magic('bash', '', 'atlasreader -a all -t 3 -p 20 stat_img.nii.gz 10\n')


# The overview plot and cluster specific plots look similar to what we saw before. Great!

# In[39]:


Image("stat_img.png")


# In[40]:


Image("stat_img_cluster01.png")


# But because we set `atlas` to `all` we now have much more information in the **peak** and **cluster tables** (you need to scroll to the right to see all of them).

# In[41]:


pd.read_csv("stat_img_peaks.csv")


# In[42]:


pd.read_csv("stat_img_clusters.csv")


# ## `queryatlas`
# 
# Last but certainly not least: If you install AtlasReader via `pip`, then you will also have access to the `queryatlas` tool. This command line interface is based on the main functions in AtlasReader and allows you to get the **atlas information** at a particular location.
# 
# Let's for example say, we want to know which atlas corresponance the MNI coordinate 32.3, -14.8, -18.5 has:

# In[43]:


get_ipython().run_cell_magic('bash', '', 'queryatlas 32.3 -14.8 -18.5\n')


# As before, you first get again the license term information, followed by the atlas name and the corresponding label at the given location.
# 
# As before, `queryatlas` gives you also some optional parameters. Use the `-h` flag to see them:

# In[44]:


get_ipython().run_cell_magic('bash', '', 'queryatlas -h\n')