{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 3: Normalize data to MNI template\n",
    "\n",
    "This example covers the normalization of data. Some people prefer to normalize the data during the preprocessing, just before smoothing. I prefer to do the 1st-level analysis completely in subject space and only normalize the contrasts for the 2nd-level analysis. But both approaches are fine.\n",
    "\n",
    "For the current example, we will take the computed 1st-level contrasts from the previous experiment (again once done with fwhm=4mm and fwhm=8mm) and normalize them into MNI-space. To show two different approaches, we will do the normalization once with ANTs and once with SPM."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "\n",
    "Before we can start with the ANTs example, we first need to download the already computed deformation field. The data can be found in the `derivatives/fmriprep` folder of the dataset and can be downloaded with the following `datalad` command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "datalad get -J 4 -d /data/ds000114 /data/ds000114/derivatives/fmriprep/sub-0[2345789]/anat/*h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** This might take a while, as datalad needs to download ~710MB of data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Alternatively: Prepare yourself\n",
    "We're using the precomputed warp field from [fmriprep](http://fmriprep.readthedocs.io), as this step otherwise would take up to 10 hours or more for all subjects to complete. If you're nonetheless interested in computing the warp parameters with ANTs yourself, without using [fmriprep](http://fmriprep.readthedocs.io), either check out the script [ANTS_registration.py](https://github.com/miykael/nipype_tutorial/blob/master/notebooks/scripts/ANTS_registration.py) or even quicker, use [RegistrationSynQuick](http://nipype.readthedocs.io/en/latest/interfaces/generated/interfaces.ants/registration.html#registrationsynquick), Nipype's implementation of `antsRegistrationSynQuick.sh`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalization with ANTs\n",
    "\n",
    "The normalization with ANTs requires that you first compute the transformation matrix that would bring the anatomical images of each subject into template space. Depending on your system this might take a few hours per subject. To facilitate this step, the transformation matrix is already computed for the T1 images.\n",
    "\n",
    "The data for it can be found under:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls /data/ds000114/derivatives/fmriprep/sub-*/anat/*h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Now, let's start with the ANTs normalization workflow!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports (ANTs)\n",
    "\n",
    "First, we need to import all the modules we later want to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os.path import join as opj\n",
    "from nipype import Workflow, Node, MapNode\n",
    "from nipype.interfaces.ants import ApplyTransforms\n",
    "from nipype.interfaces.utility import IdentityInterface\n",
    "from nipype.interfaces.io import SelectFiles, DataSink\n",
    "from nipype.interfaces.fsl import Info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experiment parameters (ANTs)\n",
    "\n",
    "It's always a good idea to specify all parameters that might change between experiments at the beginning of your script. And remember that we decided to run the group analysis without subject ``sub-01``, ``sub-06`` and ``sub-10`` because they are left-handed (see [this section](https://miykael.github.io/nipype_tutorial/notebooks/example_1stlevel.html#Special-case))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment_dir = '/output'\n",
    "output_dir = 'datasink'\n",
    "working_dir = 'workingdir'\n",
    "\n",
    "# list of subject identifiers (remember we use only right handed subjects)\n",
    "subject_list = ['02', '03', '04', '05', '07', '08', '09']\n",
    "\n",
    "# task name\n",
    "task_name = \"fingerfootlips\"\n",
    "\n",
    "# Smoothing widths used during preprocessing\n",
    "fwhm = [4, 8]\n",
    "\n",
    "# Template to normalize to\n",
    "template = '/data/ds000114/derivatives/fmriprep/mni_icbm152_nlin_asym_09c/1mm_T1.nii.gz'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note** if you're not using the corresponding docker image, than the **``template``** file might not be in your ``data``  directory. To get ``mni_icbm152_nlin_asym_09c``, either download it from this [website](https://files.osf.io/v1/resources/fvuh8/providers/osfstorage/580705089ad5a101f17944a9), unpack it and move it to ``/data/ds000114/derivatives/fmriprep/`` or run the following command in a cell:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```bash\n",
    "%%bash\n",
    "curl -L https://files.osf.io/v1/resources/fvuh8/providers/osfstorage/580705089ad5a101f17944a9 \\\n",
    "     -o /data/ds000114/derivatives/fmriprep/mni_icbm152_nlin_asym_09c.tar.gz\n",
    "     \n",
    "tar xf /data/ds000114/derivatives/fmriprep/mni_icbm152_nlin_asym_09c.tar.gz \\\n",
    "    -C /data/ds000114/derivatives/fmriprep/.\n",
    "    \n",
    "rm /data/ds000114/derivatives/fmriprep/mni_icbm152_nlin_asym_09c.tar.gz\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify Nodes (ANTs)\n",
    "\n",
    "Initiate all the different interfaces (represented as nodes) that you want to use in your workflow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply Transformation - applies the normalization matrix to contrast images\n",
    "apply2con = MapNode(ApplyTransforms(args='--float',\n",
    "                                    input_image_type=3,\n",
    "                                    interpolation='BSpline',\n",
    "                                    invert_transform_flags=[False],\n",
    "                                    num_threads=1,\n",
    "                                    reference_image=template,\n",
    "                                    terminal_output='file'),\n",
    "                    name='apply2con', iterfield=['input_image'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify input & output stream (ANTs)\n",
    "\n",
    "Specify where the input data can be found & where and how to save the output data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Infosource - a function free node to iterate over the list of subject names\n",
    "infosource = Node(IdentityInterface(fields=['subject_id', 'fwhm_id']),\n",
    "                  name=\"infosource\")\n",
    "infosource.iterables = [('subject_id', subject_list),\n",
    "                        ('fwhm_id', fwhm)]\n",
    "\n",
    "# SelectFiles - to grab the data (alternativ to DataGrabber)\n",
    "templates = {'con': opj(output_dir, '1stLevel',\n",
    "                        'sub-{subject_id}/fwhm-{fwhm_id}', '???_00??.nii'),\n",
    "             'transform': opj('/data/ds000114/derivatives/fmriprep/', 'sub-{subject_id}', 'anat',\n",
    "                              'sub-{subject_id}_t1w_space-mni152nlin2009casym_warp.h5')}\n",
    "selectfiles = Node(SelectFiles(templates,\n",
    "                               base_directory=experiment_dir,\n",
    "                               sort_filelist=True),\n",
    "                   name=\"selectfiles\")\n",
    "\n",
    "# Datasink - creates output folder for important outputs\n",
    "datasink = Node(DataSink(base_directory=experiment_dir,\n",
    "                         container=output_dir),\n",
    "                name=\"datasink\")\n",
    "\n",
    "# Use the following DataSink output substitutions\n",
    "substitutions = [('_subject_id_', 'sub-')]\n",
    "subjFolders = [('_fwhm_id_%ssub-%s' % (f, sub), 'sub-%s_fwhm%s' % (sub, f))\n",
    "               for f in fwhm\n",
    "               for sub in subject_list]\n",
    "subjFolders += [('_apply2con%s/' % (i), '') for i in range(9)] # number of contrast used in 1stlevel an.\n",
    "substitutions.extend(subjFolders)\n",
    "datasink.inputs.substitutions = substitutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify Workflow (ANTs)\n",
    "\n",
    "Create a workflow and connect the interface nodes and the I/O stream to each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initiation of the ANTs normalization workflow\n",
    "antsflow = Workflow(name='antsflow')\n",
    "antsflow.base_dir = opj(experiment_dir, working_dir)\n",
    "\n",
    "# Connect up the ANTs normalization components\n",
    "antsflow.connect([(infosource, selectfiles, [('subject_id', 'subject_id'),\n",
    "                                             ('fwhm_id', 'fwhm_id')]),\n",
    "                  (selectfiles, apply2con, [('con', 'input_image'),\n",
    "                                            ('transform', 'transforms')]),\n",
    "                  (apply2con, datasink, [('output_image', 'norm_ants.@con')]),\n",
    "                  ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the workflow (ANTs)\n",
    "\n",
    "It always helps to visualize your workflow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create ANTs normalization graph\n",
    "antsflow.write_graph(graph2use='colored', format='png', simple_form=True)\n",
    "\n",
    "# Visualize the graph\n",
    "from IPython.display import Image\n",
    "Image(filename=opj(antsflow.base_dir, 'antsflow', 'graph.png'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Run the Workflow (ANTs)\n",
    "\n",
    "Now that everything is ready, we can run the ANTs normalization workflow. Change ``n_procs`` to the number of jobs/cores you want to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "antsflow.run('MultiProc', plugin_args={'n_procs': 4})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normalization with SPM12\n",
    "\n",
    "The normalization with SPM12 is rather straightforward. The only thing we need to do is run the Normalize12 module. **So let's start!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports (SPM12)\n",
    "\n",
    "First, we need to import all the modules we later want to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os.path import join as opj\n",
    "from nipype.interfaces.spm import Normalize12\n",
    "from nipype.interfaces.utility import IdentityInterface\n",
    "from nipype.interfaces.io import SelectFiles, DataSink\n",
    "from nipype.algorithms.misc import Gunzip\n",
    "from nipype import Workflow, Node"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experiment parameters (SPM12)\n",
    "\n",
    "It's always a good idea to specify all parameters that might change between experiments at the beginning of your script. And remember that we decided to run the group analysis without subject ``sub-01``, ``sub-06`` and ``sub-10`` because they are left-handed (see [this section](https://miykael.github.io/nipype_tutorial/notebooks/example_1stlevel.html#Special-case))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment_dir = '/output'\n",
    "output_dir = 'datasink'\n",
    "working_dir = 'workingdir'\n",
    "\n",
    "# list of subject identifiers\n",
    "subject_list = ['02', '03', '04', '05', '07', '08', '09']\n",
    "\n",
    "# task name\n",
    "task_name = \"fingerfootlips\"\n",
    "\n",
    "# Smoothing withds used during preprocessing\n",
    "fwhm = [4, 8]\n",
    "\n",
    "template = '/opt/spm12-r7219/spm12_mcr/spm12/tpm/TPM.nii'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify Nodes (SPM12)\n",
    "\n",
    "Initiate all the different interfaces (represented as nodes) that you want to use in your workflow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gunzip - unzip the anatomical image\n",
    "gunzip = Node(Gunzip(), name=\"gunzip\")\n",
    "\n",
    "# Normalize - normalizes functional and structural images to the MNI template\n",
    "normalize = Node(Normalize12(jobtype='estwrite',\n",
    "                             tpm=template,\n",
    "                             write_voxel_sizes=[1, 1, 1]),\n",
    "                 name=\"normalize\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify input & output stream (SPM12)\n",
    "\n",
    "Specify where the input data can be found & where and how to save the output data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Infosource - a function free node to iterate over the list of subject names\n",
    "infosource = Node(IdentityInterface(fields=['subject_id', 'fwhm_id']),\n",
    "                  name=\"infosource\")\n",
    "infosource.iterables = [('subject_id', subject_list),\n",
    "                        ('fwhm_id', fwhm)]\n",
    "\n",
    "# SelectFiles - to grab the data (alternativ to DataGrabber)\n",
    "templates = {'con': opj(output_dir, '1stLevel',\n",
    "                        'sub-{subject_id}/fwhm-{fwhm_id}', '???_00??.nii'),\n",
    "             'anat': opj('/data/ds000114/derivatives', 'fmriprep', 'sub-{subject_id}',\n",
    "                         'anat', 'sub-{subject_id}_t1w_preproc.nii.gz')}\n",
    "\n",
    "selectfiles = Node(SelectFiles(templates,\n",
    "                               base_directory=experiment_dir,\n",
    "                               sort_filelist=True),\n",
    "                   name=\"selectfiles\")\n",
    "\n",
    "# Datasink - creates output folder for important outputs\n",
    "datasink = Node(DataSink(base_directory=experiment_dir,\n",
    "                         container=output_dir),\n",
    "                name=\"datasink\")\n",
    "\n",
    "# Use the following DataSink output substitutions\n",
    "substitutions = [('_subject_id_', 'sub-')]\n",
    "subjFolders = [('_fwhm_id_%ssub-%s' % (f, sub), 'sub-%s_fwhm%s' % (sub, f))\n",
    "               for f in fwhm\n",
    "               for sub in subject_list]\n",
    "substitutions.extend(subjFolders)\n",
    "datasink.inputs.substitutions = substitutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify Workflow (SPM12)\n",
    "\n",
    "Create a workflow and connect the interface nodes and the I/O stream to each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify Normalization-Workflow & Connect Nodes\n",
    "spmflow = Workflow(name='spmflow')\n",
    "spmflow.base_dir = opj(experiment_dir, working_dir)\n",
    "\n",
    "# Connect up SPM normalization components\n",
    "spmflow.connect([(infosource, selectfiles, [('subject_id', 'subject_id'),\n",
    "                                            ('fwhm_id', 'fwhm_id')]),\n",
    "                 (selectfiles, normalize, [('con', 'apply_to_files')]),\n",
    "                 (selectfiles, gunzip, [('anat', 'in_file')]),\n",
    "                 (gunzip, normalize, [('out_file', 'image_to_align')]),\n",
    "                 (normalize, datasink, [('normalized_files', 'norm_spm.@files'),\n",
    "                                        ('normalized_image', 'norm_spm.@image'),\n",
    "                                        ]),\n",
    "                 ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the workflow (SPM12)\n",
    "\n",
    "It always helps to visualize your workflow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create SPM normalization graph\n",
    "spmflow.write_graph(graph2use='colored', format='png', simple_form=True)\n",
    "\n",
    "# Visualize the graph\n",
    "from IPython.display import Image\n",
    "Image(filename=opj(spmflow.base_dir, 'spmflow', 'graph.png'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Run the Workflow (SPM12)\n",
    "\n",
    "Now that everything is ready, we can run the SPM normalization workflow. Change ``n_procs`` to the number of jobs/cores you want to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spmflow.run('MultiProc', plugin_args={'n_procs': 4})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison between ANTs and SPM normalization\n",
    "\n",
    "Now that we ran the normalization with ANTs and SPM, let us compare their output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.plotting import plot_stat_map\n",
    "%matplotlib inline\n",
    "anatimg = '/data/ds000114/derivatives/fmriprep/mni_icbm152_nlin_asym_09c/1mm_T1.nii.gz'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's compare the normalization of the **anatomical** images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_stat_map(\n",
    "    '/data/ds000114/derivatives/fmriprep/sub-02/anat/sub-02_t1w_space-mni152nlin2009casym_preproc.nii.gz',\n",
    "    title='anatomy - ANTs (normalized to ICBM152)', bg_img=anatimg,\n",
    "    threshold=200, display_mode='ortho', cut_coords=(-50, 0, -10));\n",
    "plot_stat_map(\n",
    "    '/output/datasink/norm_spm/sub-02_fwhm4/wsub-02_t1w_preproc.nii',\n",
    "    title='anatomy - SPM (normalized to SPM\\'s TPM)', bg_img=anatimg,\n",
    "    threshold=200, display_mode='ortho', cut_coords=(-50, 0, -10));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And what about the **contrast** images for **Finger > others**?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_stat_map(\n",
    "    '/output/datasink/norm_ants/sub-02_fwhm8/con_0005_trans.nii', title='contrast5 - fwhm=8 - ANTs',\n",
    "    bg_img=anatimg, threshold=2, vmax=5, display_mode='ortho', cut_coords=(-39, -37, 56));\n",
    "plot_stat_map(\n",
    "    '/output/datasink/norm_spm/sub-02_fwhm8/wcon_0005.nii', title='contrast5 - fwhm=8 - SPM',\n",
    "    bg_img=anatimg, threshold=2, vmax=5, display_mode='ortho', cut_coords=(-39, -37, 56));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.plotting import plot_glass_brain\n",
    "plot_glass_brain(\n",
    "    '/output/datasink/norm_ants/sub-02_fwhm8/con_0005_trans.nii', colorbar=True,\n",
    "    threshold=3, display_mode='lyrz', black_bg=True, vmax=6, title='contrast5 - fwhm=8 - ANTs')\n",
    "plot_glass_brain(\n",
    "    '/output/datasink/norm_spm/sub-02_fwhm8/wcon_0005.nii', colorbar=True,\n",
    "    threshold=3, display_mode='lyrz', black_bg=True, vmax=6, title='contrast5 - fwhm=8 - SPM');"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}