notebooks/warp_individual_to_templates.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Register Individual Macaque MRI Data to Templates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Authors:** Nikoloz Sirmpilatze (nsirmpilatze@dpz.eu) & Chris Klink (c.klink@nin.knaw.nl)                \n",
    "\n",
    "**Last updated:** May 11th, 2020     \n",
    "\n",
    "**Requirements:**    \n",
    "* _python_ >= 3.7\n",
    "* _nipype_ >= 1.2.0\n",
    "* _nilearn_ >= 0.5.2\n",
    "    * Used only for visualisation\n",
    "* _nibabel_ >= 2.3.3\n",
    "* _joblib_ >= 0.14.1\n",
    "    * Used only for parallel processing, not necessary when registrations are done serially\n",
    "* _ANTs_ >= 2.4.0\n",
    "    * _antsRegistration_,  _antsApplyTransforms_ and _antsAverageImages_ need to be in your path as executables    \n",
    "\n",
    "**Citation**: Sirmpilatze, Nikoloz and Klink, P. Christiaan (2020). RheMAP: Non-linear warps between common rhesus macaque brain templates (Version 1.2)[Data set]. Zenodo. https://doi.org/10.5281/zenodo.3786357    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following templates are used: \n",
    "1. [NMT v1.2](https://github.com/jms290/NMT)\n",
    "2. [NMT v1.3](https://github.com/jms290/NMT)\n",
    "3. [D99](https://afni.nimh.nih.gov/Macaque)\n",
    "4. [INIA19](https://www.nitrc.org/projects/inia19/https://www.nitrc.org/projects/inia19/)\n",
    "5. [MNI macaque](http://www.bic.mni.mcgill.ca/ServicesAtlases/Macaque)\n",
    "6. [Yerkes19](https://github.com/Washington-University/NHPPipelines)   \n",
    "7. [NMT v2.0 (4 templates)](https://github.com/jms290/NMT) (NB! at the time of writing, NMTv2.0 was not online yet)      \n",
    "\n",
    "**NB!** We do not provide copies of the actual templates (licenses often forbids redistribution), but instead suggest you follow the links above and get them at the source. We do offer the warp files and warped templates that will be produced by this workflow. They can be downloaded from [Zenodo](https://doi.org/10.5281/zenodo.3786357). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to emulate this code, you can consider setting up your templates in the following folder structure:    \n",
    "\n",
    "|--- RheMAP   \n",
    "&emsp; |--- notebooks     \n",
    "&emsp; |--- templates     \n",
    "&emsp; &emsp; |--- D99     \n",
    "&emsp; &emsp; &emsp; |--- D99_atlas_1.2a.nii.gz     \n",
    "&emsp; &emsp; &emsp; |--- D99_atlas_1.2a_in_MNI.nii.gz     \n",
    "&emsp; &emsp; &emsp; |--- D99_template.nii.gz      \n",
    "&emsp; &emsp; |--- INIA   \n",
    "&emsp; &emsp; &emsp; |--- inia19-t1-brain_truncated.nii.gz      \n",
    "&emsp; &emsp; |--- MNI     \n",
    "&emsp; &emsp; &emsp; |--- macaque_25_model-MNI_brain.nii.gz      \n",
    "&emsp; &emsp; |--- NMT     \n",
    "&emsp; &emsp; &emsp; |--- NMT_v1.2     \n",
    "&emsp; &emsp; &emsp; &emsp; |--- NMT_SS.nii.gz      \n",
    "&emsp; &emsp; &emsp; |--- NMT_v1.3     \n",
    "&emsp; &emsp; &emsp; &emsp; |--- NMT_SS.nii.gz     \n",
    "&emsp; &emsp; |--- YRK      \n",
    "&emsp; &emsp; &emsp; |--- MacaqueYerkes19_T1w_0.5mm_brain.nii.gz       \n",
    "&emsp; |--- subjects     \n",
    "&emsp; &emsp; |--- SUBJ-01     \n",
    "&emsp; &emsp; &emsp; |--- SUBJ-01_T1.nii.gz\n",
    "&emsp; &emsp; |--- SUBJ-02      \n",
    "&emsp; &emsp; &emsp; |--- SUBJ-01_T1.nii.gz     \n",
    "\n",
    "After downloading the warp files and warped templates from [Zenodo](https://zenodo.org/record/3776856#.XqqfI3UzZjE), we suggest you include them like this:      \n",
    "\n",
    "|--- RheMAP   \n",
    "&emsp; |--- notebooks     \n",
    "&emsp; |--- templates     \n",
    "&emsp; |--- warps       \n",
    "&emsp; &emsp; |--- final     \n",
    "&emsp; &emsp; &emsp; |--- D99_to_INIA_CompositeWarp.nii.gz     \n",
    "&emsp; &emsp; &emsp; |--- D99_to_MNI_CompositeWarp.nii.gz      \n",
    "&emsp; &emsp; &emsp; |--- etc   \n",
    "&emsp; &emsp; |--- linear     \n",
    "&emsp; &emsp; &emsp; |--- D99_to_INIA_affine_0GenericAffine.mat     \n",
    "&emsp; &emsp; &emsp; |--- D99_to_MNI_affine_0GenericAffine.mat     \n",
    "&emsp; &emsp; &emsp; |--- etc   \n",
    "&emsp; &emsp; |--- nonlinear     \n",
    "&emsp; &emsp; &emsp; |--- D99_to_INIA_1InverseWarp.nii.gz     \n",
    "&emsp; &emsp; &emsp; |--- D99_to_INIA_1Warp.nii.gz     \n",
    "&emsp; &emsp; &emsp; |--- D99_to_MNI_1InverseWarp.nii.gz     \n",
    "&emsp; &emsp; &emsp; |--- D99_to_MNI_1Warp.nii.gz  \n",
    "&emsp; &emsp; &emsp; |--- etc   \n",
    "&emsp; |--- warped_templates       \n",
    "&emsp; &emsp; |--- final     \n",
    "&emsp; &emsp; &emsp; |--- D99_in_INIA_composite.nii.gz     \n",
    "&emsp; &emsp; &emsp; |--- D99_in_MNI_composite.nii.gz       \n",
    "&emsp; &emsp; &emsp; |--- etc   \n",
    "&emsp; &emsp; |--- linear     \n",
    "&emsp; &emsp; &emsp; |--- D99_in_INIA_linear.nii.gz    \n",
    "&emsp; &emsp; &emsp; |--- D99_in_MNI_linear.nii.gz          \n",
    "&emsp; &emsp; &emsp; |--- etc   \n",
    "&emsp; &emsp; |--- nonlinear     \n",
    "&emsp; &emsp; &emsp; |--- D99_in_INIA_linear+SyN.nii.gz     \n",
    "&emsp; &emsp; &emsp; |--- D99_in_MNI_linear+SyN.nii.gz     \n",
    "&emsp; &emsp; &emsp; |--- etc  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the warp files and warped templates you could of course reconstruct the original templates with something like the following.     \n",
    "\n",
    "On the command line:     \n",
    "```bash\n",
    "antsApplyTransforms -i <TEMPLATE1_in_TEMPLATE2_composite.nii.gz> \\\n",
    "                    -r <TEMPLATE2_in_TEMPLATE1.nii.gz> \\        \n",
    "                    -o <RECONSTRUCTED_ORIGINAL_TEMPLATE1.nii.gz> \\\n",
    "                    -t [<TEMPLATE1_to_TEMPLATE2_CompositeWarp>,1] \\\n",
    "                    -n Linear \\\n",
    "                    -d 3\n",
    "```    \n",
    "\n",
    "In NiPype:     \n",
    "```python\n",
    "import nipype.interfaces.ants as ants    \n",
    "ants.ApplyTransforms(\n",
    "            input_image=<TEMPLATE1_in_TEMPLATE2_composite.nii.gz>,\n",
    "            reference_image=<TEMPLATE2_in_TEMPLATE1.nii.gz>,        \n",
    "            output_image=<RECONSTRUCTED_ORIGINAL_TEMPLATE1.nii.gz>,\n",
    "            transforms=<TEMPLATE1_to_TEMPLATE2_CompositeWarp>,\n",
    "            invert_transform_flags=True,\n",
    "            interpolation='Linear',\n",
    "            dimension=3)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 0: Preparations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0a. Import required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import time\n",
    "import shutil as sh\n",
    "import nibabel as nb\n",
    "\n",
    "from itertools import combinations\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import nipype.interfaces.fsl as fsl # nipype interface for FSL\n",
    "import nipype.interfaces.ants as ants # nipype interface for ANTs\n",
    "from nilearn.plotting import plot_anat # Plotting function from nilearn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0b. Parallel or serial    \n",
    "Registration can be a time consuming process. If you have the resources, you may opt to parallelize some of the steps. You can configure that here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "do_parallel = True # Set to True for parallel processing\n",
    "par_numjobs = 5     # How many parallel processes? This is only used when do_parallel is True\n",
    "if do_parallel is True:\n",
    "    from joblib import Parallel, delayed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0c. Define relative paths to files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "THe preferred way is to use the skull-stripped isotropic volumetric images for registration.\n",
    "* For **NMT**, **D99**, and **YRK** the provided skull-stripped brains were used\n",
    "* In **INIA** the brain stem extends further down the spinal cord compared to the other temlplates, so the braisn stem was truncated at a level similar to the others\n",
    "* No skull-stripped image was provided with **MNI**, so the brain was segmented semi-manually using ITK-SNAP\n",
    "* For **YRK** we used the version provided together with [NHPPipelines](https://github.com/Washington-University/NHPPipelines)\n",
    "\n",
    "If you want to do skullstripping of your single sibject within this notebook, you should also have a full-head template and a corresponding brainmask. Most templates provide these.     \n",
    "\n",
    "First let's specify the template paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ============================================\n",
    "# NB! THIS COULD BE DIFFERENT FOR EVERY USER\n",
    "# ============================================\n",
    "BASE_path     = os.path.dirname(os.getcwd()) + '/'   # repo base folder\n",
    "\n",
    "# these follow the directory structure as outline above\n",
    "TEMPLATE_path = BASE_path + 'templates/'                   # templates base folder\n",
    "NMTv12_path   = TEMPLATE_path + 'NMT/NMT_v1.2/'\n",
    "NMTv13_path   = TEMPLATE_path + 'NMT/NMT_v1.3/'\n",
    "NMTv20_path   = TEMPLATE_path + 'NMT/NMT_v2.0/'\n",
    "D99_path      = TEMPLATE_path + 'D99/'\n",
    "INIA_path     = TEMPLATE_path + 'INIA/'\n",
    "MNI_path      = TEMPLATE_path + 'MNI/'\n",
    "YRK_path      = TEMPLATE_path + 'YRK/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NMTv12_brain = NMTv12_path + 'NMT_SS.nii.gz'\n",
    "NMTv13_brain = NMTv13_path + 'NMT_SS.nii.gz'\n",
    "NMTv20_acpc_sym_brain = NMTv20_path + 'NMT_acpc_sym_2.0_SS.nii.gz'\n",
    "NMTv20_acpc_asym_brain = NMTv20_path + 'NMT_acpc_asym_2.0_SS.nii.gz'\n",
    "NMTv20_stereo_sym_brain = NMTv20_path + 'NMT_stereo_sym_2.0_SS.nii.gz'\n",
    "NMTv20_stereo_asym_brain = NMTv20_path + 'NMT_stereo_asym_2.0_SS.nii.gz'\n",
    "D99_brain    = D99_path + 'D99_template.nii.gz'\n",
    "INIA_brain   = INIA_path + 'inia19-t1-brain_truncated.nii.gz'\n",
    "MNI_brain    = MNI_path + 'macaque_25_model-MNI_brain.nii.gz'\n",
    "YRK_brain    = YRK_path + 'MacaqueYerkes19_T1w_0.5mm_brain.nii.gz'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now create a list of only the templates you want to calculate warps for"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # all\n",
    "# temp_names   = ['NMTv1.2', 'NMTv1.3',\n",
    "#                 'NMTv2.0_acpc_sym','NMTv2.0_acpc_asym',\n",
    "#                 'NMTv2.0_stereo_sym','NMTv2.0_stereo_asym',\n",
    "#                 'D99', 'INIA', 'MNI', 'YRK']\n",
    "# temp_brains  = [NMTv12_brain, NMTv13_brain, \n",
    "#                 NMTv20_acpc_sym_brain, NMTv20_acpc_asym_brain,\n",
    "#                 NMTv20_stereo_sym_brain,NMTv20_stereo_asym_brain,\n",
    "#                 D99_brain, INIA_brain, MNI_brain, YRK_brain]\n",
    "\n",
    "# selected\n",
    "temp_names   = ['NMTv1.3','NMTv2.0_stereo_sym']\n",
    "temp_brains  = [NMTv13_brain,NMTv20_stereo_sym_brain]\n",
    "\n",
    "# strength of contrast in the plot_anat function (\"dim\")\n",
    "temp_contrasts = [0]*len(temp_brains) # create a list of zeros because this should exist\n",
    "#temp_contrasts = [-1, 0, 0, 0, 0, 0, 0, -1, -0.5, 0] # specify to improve display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the individual anatomical(s) you would like to warp to a template space. Ideally you would use a skull-stripped version with the highest resolution and contrast you have available (consider averaging multiple scans). There are various ways to perform skullstripping on your data (see the [PRIME-RE website](https://prime-re.github.io/). You can also opt for the RheMAP built-in solution and let us do the skull-stripping for you, but this is a fast-over-precision soution and your mileage may vary. The RheMAP method requires that you have a template that includes a full-head (non-skull-stripped) scan and a brainmask volume (or a brain volume you can change into a mask by binarizing, e.g. with `fslmaths BRAIN.nii.gz -bin BRANMASK.nii.gz` (but you'll have to do this manually first)     \n",
    "\n",
    "In short the RheMAP method does the following:     \n",
    "- Linearly register the full head single-subject scan to a full-head template scan\n",
    "- Transform the template brainmask to single-subject space with the inverse tranform\n",
    "- Apply the transformed mask to skullstrip the single subject data    \n",
    "\n",
    "You can either skullstrip all subject in a run of the notebook or none, not some of them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SUBJECT_path = BASE_path + 'subjects/' \n",
    "SUBJ_01_path = SUBJECT_path + 'Danny/'\n",
    "SUBJ_02_path = SUBJECT_path + 'Eddy/'\n",
    "subj_paths = [SUBJ_01_path, SUBJ_02_path]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DoSkullStrip = True\n",
    "if DoSkullStrip:\n",
    "    SUBJ_01_head = SUBJ_01_path + 'Danny.nii.gz'\n",
    "    SUBJ_02_head = SUBJ_02_path + 'Eddy2019.nii.gz'\n",
    "    REF_head = NMTv13_path + 'NMT.nii.gz'\n",
    "    REF_brainmask = NMTv13_path + 'NMT_brainmask.nii.gz'\n",
    "    subj_heads = [SUBJ_01_head, SUBJ_02_head]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Concatenate the templates and subject into one list for easy handling (within this notebook we will later exclude template-to-tamplate warping)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_names = subj_names + temp_names\n",
    "all_brains = subj_brains + temp_brains "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0d. Skullstripping (optional)    \n",
    "If you are not doing skullstripping with the RheMAP method, the brain files linked in the next cell should already exist.       \n",
    "If skullstripping is performed, they are generated (and overwritten if they already exist).    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SUBJ_01_brain = SUBJ_01_path + 'Danny_brain.nii.gz'\n",
    "SUBJ_02_brain = SUBJ_02_path + 'Eddy2019_brain.nii.gz'\n",
    "\n",
    "subj_names   = ['Danny', 'Eddy2019']\n",
    "subj_brains  = [SUBJ_01_brain, SUBJ_02_brain]\n",
    "\n",
    "subj_contrasts = [0]*len(subj_brains) # create a list of zeros because this should exist\n",
    "all_contrasts = subj_contrasts + temp_contrasts "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def skullstrip_subject(subj_name, subj_head, subj_brain, subj_path,\n",
    "                       ref_head, ref_brainmask):\n",
    "\n",
    "    \"\"\"Skullstrips an individual anatomical by coregistering with an existing mask\n",
    "\n",
    "    Parameters\n",
    "    ----------   \n",
    "    subj_name: string\n",
    "        name of the subject being skullstrippedor 'target' is on this list\n",
    "    subj_head: 3D NIFTI image file of subject head\n",
    "    subj_brain: Name of a 3D NIFTI image file that will be the skullstripped result\n",
    "    ref_head: 3D NIFTI image file of reference head\n",
    "    \n",
    "    Return\n",
    "    ------\n",
    "    Skullstripped 3D NIFTI\n",
    "    \"\"\"\n",
    "\n",
    "    # linear registration \n",
    "    rgstr = ants.registration.Registration(\n",
    "               dimension=3,\n",
    "               float=False,\n",
    "               fixed_image=ref_head,\n",
    "               moving_image=subj_head,\n",
    "               output_transform_prefix='{0}_to_ref_affine_'.format(subj_name),\n",
    "               output_warped_image='{0}_in_ref_linear.nii.gz'.format(subj_name),\n",
    "               initial_moving_transform_com=0,\n",
    "               winsorize_lower_quantile=0.05,\n",
    "               winsorize_upper_quantile=0.95,\n",
    "               interpolation='Linear',\n",
    "               use_histogram_matching=[True],\n",
    "               transforms=['Affine'],\n",
    "               transform_parameters=[(0.1,)],\n",
    "               metric=['MI'],\n",
    "               metric_weight=[1],\n",
    "               radius_or_number_of_bins=[32],\n",
    "               sampling_strategy=['Regular'],\n",
    "               sampling_percentage=[0.2],\n",
    "               number_of_iterations=[[1000, 500, 250, 100]],\n",
    "               convergence_threshold=[1e-6],\n",
    "               convergence_window_size=[10],\n",
    "               shrink_factors=[[8, 4, 2, 1]],\n",
    "               smoothing_sigmas=[[3, 2, 1, 0]],\n",
    "               sigma_units=['vox'])\n",
    "    \n",
    "    rgstr.run()\n",
    "    \n",
    "    # apply inverse transform to mask\n",
    "    invreg_mask = ants.ApplyTransforms(\n",
    "            input_image=ref_brainmask,\n",
    "            reference_image=subj_head,        \n",
    "            output_image='{0}_brainmask.nii.gz'.format(subj_name),\n",
    "            transforms='{0}_to_ref_affine_0GenericAffine.mat'.format(subj_name),\n",
    "            invert_transform_flags=True,\n",
    "            interpolation='Linear',\n",
    "            dimension=3)\n",
    "    \n",
    "    invreg_mask.run()\n",
    "    \n",
    "    # dilate mask a bit to be sure we don't cut into brain\n",
    "    dil_mask = fsl.maths.DilateImage(in_file='{0}_brainmask.nii.gz'.format(subj_name),operation='mean',\n",
    "                                     out_file='{0}_brainmask.nii.gz'.format(subj_name))\n",
    "    dil_mask.run()\n",
    "    \n",
    "    # mask\n",
    "    apply_mask = fsl.maths.ApplyMask(\n",
    "            in_file=subj_head,\n",
    "            mask_file='{0}_brainmask.nii.gz'.format(subj_name),        \n",
    "            out_file=subj_brain,\n",
    "            )\n",
    "    \n",
    "    apply_mask.run() \n",
    "    \n",
    "    # clean up\n",
    "    os.remove('{0}_to_ref_affine_0GenericAffine.mat'.format(subj_name))\n",
    "    os.remove('{0}_in_ref_linear.nii.gz'.format(subj_name))\n",
    "    os.rename(os.getcwd() + '/' + '{0}_brainmask.nii.gz'.format(subj_name),\n",
    "              subj_path + '/{0}_brainmask.nii.gz'.format(subj_name));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if DoSkullStrip:\n",
    "    if do_parallel: # parallel\n",
    "        print('\\nPerforming skullstripping in parallel...')\n",
    "        start_time = time.time()\n",
    "        Parallel(n_jobs=par_numjobs, backend='multiprocessing')(\n",
    "            delayed(skullstrip_subject)\n",
    "            (subj_name, subj_head, subj_brain, subj_path, REF_head, REF_brainmask) \n",
    "            for subj_name, subj_head, subj_brain, subj_path in zip(subj_names, subj_heads, subj_brains, subj_paths))\n",
    "        end_time = time.time()\n",
    "        print('Skullstripping completed in  {0:.2f} minutes'.format((end_time - start_time)/60))\n",
    "        \n",
    "    else: # serial\n",
    "        for subj_name, subj_head, subj_brain, subj_path in zip(subj_names, subj_heads, subj_brains, subj_paths):\n",
    "            print('Skullstripping (serial) ' + subj_name) \n",
    "            start_time = time.time()\n",
    "            skullstrip_subject(subj_name, subj_head, subj_brain, subj_path, REF_head, REF_brainmask)\n",
    "            end_time = time.time()\n",
    "            print('Skullstripping completed in  {0:.2f} minutes'.format((end_time - start_time)/60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0e. Plot individual subject brains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, brain,con in zip(subj_names, subj_brains, subj_contrasts):\n",
    "    plot_anat(brain, display_mode='ortho', title=name, \n",
    "              draw_cross=False, annotate=True, dim=con);\n",
    "    plt.savefig(name + '_brain.png'); # save png\n",
    "    plt.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Linear registration\n",
    "* This is performed using [ANTs](http://stnava.github.io/ANTs/) **Affine** registration (12 degrees of freedom: 3 translations + 3 rotations + 3 scalings + 3 shears)\n",
    "* The _antsRegistration_ function is called through its [nipype wrapper](https://nipype.readthedocs.io/en/latest/interfaces/generated/interfaces.ants/registration.html#registration).\n",
    "* [Anatomy of an antsRegistration call](https://github.com/ANTsX/ANTs/wiki/Anatomy-of-an-antsRegistration-call) does a good job of explaining what the various parameters mean."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1a. Define function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def register_linear(moving, target,\n",
    "                    moving_prefix, target_prefix,\n",
    "                    subj_names,include_sub2sub=False):\n",
    "    \"\"\"Linearly registers moving image to target image\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    moving: 3D NIFTI image file\n",
    "    target: 3D NIFTI image file\n",
    "    moving_prefix: string\n",
    "        one of the template/subj names \n",
    "        is used in the naming of output files\n",
    "    target_prefix: string\n",
    "        one of the template/subj names \n",
    "        is used in the naming of output files\n",
    "    subj_names: list\n",
    "        list of ssubjects to include in execution\n",
    "        will only run registration if 'moving' or 'target' is on this list\n",
    "    include_sub2sub: boolean\n",
    "        Boolean that determines whether warps between subjects are calculated\n",
    "        Default is false (not calculated)\n",
    "    \n",
    "    Return\n",
    "    ------\n",
    "    Linear transformation matrix\n",
    "        file named as '{moving_prefix}_to_{target_prefix}_0GenericAffine.mat'\n",
    "    Transformed image\n",
    "        moving image transformed into target image space\n",
    "        file named as: {moving_prefix}_in_{target_prefix}_linear.nii.gz\n",
    "    \"\"\"\n",
    "\n",
    "    # Check whether this combinations should be done\n",
    "    DoFunc=False\n",
    "    if include_sub2sub and (moving_prefix in subj_names or target_prefix in subj_names):\n",
    "        DoFunc=True\n",
    "    if not include_sub2sub and ((moving_prefix in subj_names and target_prefix not in subj_names) or (moving_prefix not in subj_names and target_prefix in subj_names)):\n",
    "        DoFunc=True\n",
    "           \n",
    "    # run it   \n",
    "    if DoFunc: \n",
    "        rgstr = ants.registration.Registration(\n",
    "               dimension=3,\n",
    "               float=False,\n",
    "               fixed_image=target,\n",
    "               moving_image=moving,\n",
    "               output_transform_prefix='{0}_to_{1}_affine_'.format(moving_prefix, target_prefix),\n",
    "               output_warped_image='{0}_in_{1}_linear.nii.gz'.format(moving_prefix, target_prefix),\n",
    "               initial_moving_transform_com=0,\n",
    "               winsorize_lower_quantile=0.05,\n",
    "               winsorize_upper_quantile=0.95,\n",
    "               interpolation='Linear',\n",
    "               use_histogram_matching=[True],\n",
    "               transforms=['Affine'],\n",
    "               transform_parameters=[(0.1,)],\n",
    "               metric=['MI'],\n",
    "               metric_weight=[1],\n",
    "               radius_or_number_of_bins=[32],\n",
    "               sampling_strategy=['Regular'],\n",
    "               sampling_percentage=[0.2],\n",
    "               number_of_iterations=[[1000, 500, 250, 100]],\n",
    "               convergence_threshold=[1e-6],\n",
    "               convergence_window_size=[10],\n",
    "               shrink_factors=[[8, 4, 2, 1]],\n",
    "               smoothing_sigmas=[[3, 2, 1, 0]],\n",
    "               sigma_units=['vox'])\n",
    "    \n",
    "        rgstr.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1b. Apply function to register all unique template pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if do_parallel: # parallel\n",
    "    print('\\nPerforming linear registration in parallel...')\n",
    "    start_time = time.time()\n",
    "    Parallel(n_jobs=par_numjobs, backend='multiprocessing')(\n",
    "        delayed(register_linear)\n",
    "        (brains[0], brains[1], names[0],names[1], subj_names, include_sub2sub=False) \n",
    "        for names, brains in zip(combinations(all_names, 2), combinations(all_brains, 2)))\n",
    "    end_time = time.time()\n",
    "    print('Linear registration completed in  {0:.2f} minutes'.format((end_time - start_time)/60))\n",
    "        \n",
    "else: # serial\n",
    "    for names, brains in zip(combinations(all_names, 2), combinations(all_brains, 2)):\n",
    "        \n",
    "        # explicit here for clarity\n",
    "        mov_pref = names[0]\n",
    "        targ_pref = names[1]\n",
    "        mov_brain = brains[0]\n",
    "        targ_brain = brains[1]\n",
    "    \n",
    "        print('\\nLinearly registering {0} to {1} template...'.format(mov_pref, targ_pref))\n",
    "        start_time = time.time()\n",
    "        register_linear(mov_brain, targ_brain, mov_pref, targ_pref, subj_names, include_sub2sub=False)\n",
    "        end_time = time.time()\n",
    "        print('Linear registration completed in  {0:.2f} minutes'.format((end_time - start_time)/60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1c. Visualize linear registration results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for names, brains in zip(combinations(all_names, 2), combinations(all_brains, 2)):\n",
    "\n",
    "    mov_pref = names[0]\n",
    "    targ_pref = names[1]\n",
    "    mov_brain = brains[0]\n",
    "    targ_brain = brains[1]\n",
    "    targ_index = all_names.index(targ_pref)\n",
    "\n",
    "    # only do this for exiting pairs\n",
    "    if os.path.isfile(os.getcwd() + '/{0}_in_{1}_linear.nii.gz'.format(mov_pref, targ_pref)):      \n",
    "        display = plot_anat(targ_brain, display_mode='ortho',\n",
    "                            title='{0} edges on {1} volume'.format(mov_pref, targ_pref),\n",
    "                            draw_cross=False, annotate=False, dim=all_contrasts[targ_index])\n",
    "        brain = '{0}_in_{1}_linear.nii.gz'.format(mov_pref, targ_pref)\n",
    "        display.add_edges(brain)\n",
    "        plt.savefig('Linear_{0}_on_{1}.png'.format(mov_pref, targ_pref))\n",
    "        plt.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Non-linear registration\n",
    "* Symmetric Diffeomorphic ([SyN](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2276735/)) registration is performed using [ANTs](http://stnava.github.io/ANTs/).\n",
    "* The starting point for this step is the transformed image from the previous (linear) step\n",
    "* The _antsRegistration_ function is called through its [nipype wrapper](https://nipype.readthedocs.io/en/latest/interfaces/generated/interfaces.ants/registration.html#registration).\n",
    "* [Anatomy of an antsRegistration call](https://github.com/ANTsX/ANTs/wiki/Anatomy-of-an-antsRegistration-call) does a good job of explaining what the various parameters mean."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2a. Define function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def register_SyN(moving, target,\n",
    "                 moving_prefix, target_prefix,\n",
    "                 subj_names,include_sub2sub=False):\n",
    "    \"\"\"Non-linearly registers moving image to target image using ANTs SyN\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    Same as register_linear function\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    Non-linear forward warp\n",
    "        file named as '{moving_prefix}_to_{target_prefix}_1Warp.nii.gz'\n",
    "    Non-linear backward warp (inverse of the above)\n",
    "        file named as '{moving_prefix}_to_{target_prefix}_1InverseWarp.nii.gz'\n",
    "    Transformed image\n",
    "        moving image transformed into target image space\n",
    "        file named as: {moving_prefix}_in_{target_prefix}_linear+SyN.nii.gz\n",
    "    \"\"\"\n",
    "    \n",
    "    # Check whether this combinations should be done\n",
    "    DoFunc=False\n",
    "    if include_sub2sub and (moving_prefix in subj_names or target_prefix in subj_names):\n",
    "        DoFunc=True\n",
    "    if not include_sub2sub and ((moving_prefix in subj_names and target_prefix not in subj_names) or (moving_prefix not in subj_names and target_prefix in subj_names)):\n",
    "        DoFunc=True\n",
    "           \n",
    "    # run it   \n",
    "    if DoFunc: \n",
    "        rgstr = ants.registration.Registration(\n",
    "               dimension=3,\n",
    "               float=False,\n",
    "               fixed_image=target,\n",
    "               moving_image=moving,\n",
    "               output_transform_prefix='{0}_to_{1}_'.format(moving_prefix, target_prefix),\n",
    "               output_warped_image='{0}_in_{1}_linear+SyN.nii.gz'.format(moving_prefix, target_prefix),\n",
    "               initial_moving_transform_com=0,\n",
    "               winsorize_lower_quantile=0.05,\n",
    "               winsorize_upper_quantile=0.95,\n",
    "               interpolation='Linear',\n",
    "               use_histogram_matching=[True],\n",
    "               transforms=['SyN'],\n",
    "               transform_parameters=[(0.1,)],\n",
    "               metric=['CC'],\n",
    "               metric_weight=[1],\n",
    "               radius_or_number_of_bins=[4],\n",
    "               sampling_strategy=['Regular'],\n",
    "               sampling_percentage=[0.2],\n",
    "               number_of_iterations=[[500, 200, 50]],\n",
    "               convergence_threshold=[1e-6],\n",
    "               convergence_window_size=[10],\n",
    "               shrink_factors=[[8, 4, 2]],\n",
    "               smoothing_sigmas=[[3, 2, 1]],\n",
    "               sigma_units=['vox'])\n",
    "\n",
    "        rgstr.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2b. Apply function to register all unique template pairs    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if do_parallel: # parallel\n",
    "    print('\\nPerforming non-linear registration in parallel...')\n",
    "    start_time = time.time()\n",
    "    Parallel(n_jobs=par_numjobs, backend='multiprocessing')(\n",
    "        delayed(register_SyN)\n",
    "        ('{0}_in_{1}_linear.nii.gz'.format(names[0], names[1]), brains[1], names[0], names[1], subj_names, include_sub2sub=False) \n",
    "        for names, brains in zip(combinations(all_names, 2), combinations(all_brains, 2)))\n",
    "    end_time = time.time()\n",
    "    print('Non-linear registration completed in  {0:.2f} minutes'.format((end_time - start_time)/60))\n",
    "    \n",
    "else: # serial\n",
    "    for names, brains in zip(combinations(all_names, 2), combinations(all_brains, 2)):\n",
    "\n",
    "        # explicit here for clarity\n",
    "        mov_pref = names[0]\n",
    "        targ_pref = names[1]\n",
    "        mov_brain = '{0}_in_{1}_linear.nii.gz'.format(mov_pref, targ_pref)\n",
    "        targ_brain = brains[1]\n",
    "    \n",
    "        print('\\nNon-linearly registering {0} to {1} template...'.format(mov_pref, targ_pref))\n",
    "        start_time = time.time()\n",
    "        register_SyN(mov_brain, targ_brain, mov_pref, targ_pref, subj_names, include_sub2sub=False)\n",
    "        end_time = time.time()\n",
    "        print('Non-linear registration completed in  {0:.2f} minutes'.format((end_time - start_time)/60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2c: Visualize non-linear registration results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for names, brains in zip(combinations(all_names, 2), combinations(all_brains, 2)):\n",
    "\n",
    "    mov_pref = names[0]\n",
    "    targ_pref = names[1]\n",
    "    targ_brain = brains[1]\n",
    "    targ_index = all_names.index(targ_pref)\n",
    "\n",
    "        # only do this for exiting pairs\n",
    "    if os.path.isfile(os.getcwd() + '/{0}_in_{1}_linear+SyN.nii.gz'.format(mov_pref, targ_pref)):\n",
    "        display = plot_anat(targ_brain, display_mode='ortho',\n",
    "                            title='{0} edges on {1} volume'.format(mov_pref, targ_pref),\n",
    "                            draw_cross=False, annotate=False, dim=all_contrasts[targ_index])\n",
    "        brain = '{0}_in_{1}_linear+SyN.nii.gz'.format(mov_pref, targ_pref)\n",
    "        display.add_edges(brain)\n",
    "        plt.savefig('Nonlinear_{0}_on_{1}.png'.format(mov_pref, targ_pref))\n",
    "        plt.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "## Step 3: Derive composite warps by combining linear and non-linear steps\n",
    "* This is done with the _antsApplyTransforms_ function, using its [nipype wrapper](https://nipype.readthedocs.io/en/latest/interfaces/generated/interfaces.ants/resampling.html#applytransforms).\n",
    "* **Caution:** the order with which transforms are passed is akin to linear algebra notation, so first transform is applied last\n",
    "* Both forward (templateA-to-templateB) and backward (templateB-to-templateA) warps are computed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3a. Define function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def combine_transforms(moving, target,\n",
    "                       moving_prefix, target_prefix,\n",
    "                       subj_names,include_sub2sub=False):\n",
    "\n",
    "    \"\"\"Combines linear transform with non-linear warp to generate a composite warp\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    Same as register_linear function\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    Composite forward warp\n",
    "        file named as '{moving_prefix}_to_{target_prefix}_CompositeWarp.nii.gz'\n",
    "    Composite backward warp (inverse of the above)\n",
    "        file named as '{target_prefix}_to_{moving_prefix}_CompositeWarp.nii.gz'\n",
    "    \"\"\"\n",
    "    \n",
    "    # Check whether this combinations should be done\n",
    "    DoFunc=False\n",
    "    if include_sub2sub and (moving_prefix in subj_names or target_prefix in subj_names):\n",
    "        DoFunc=True\n",
    "    if not include_sub2sub and ((moving_prefix in subj_names and target_prefix not in subj_names) or (moving_prefix not in subj_names and target_prefix in subj_names)):\n",
    "        DoFunc=True\n",
    "           \n",
    "    # run it   \n",
    "    if DoFunc: \n",
    "        forward = ants.ApplyTransforms(\n",
    "                input_image=moving,\n",
    "                output_image='{0}_to_{1}_CompositeWarp.nii.gz'.format(moving_prefix, target_prefix),\n",
    "                reference_image=target,\n",
    "                transforms=['{0}_to_{1}_1Warp.nii.gz'.format(moving_prefix, target_prefix),\n",
    "                            '{0}_to_{1}_affine_0GenericAffine.mat'.format(moving_prefix, target_prefix)],\n",
    "                dimension=3,\n",
    "                print_out_composite_warp_file=True)\n",
    "        forward.run()\n",
    "\n",
    "        backward = ants.ApplyTransforms(\n",
    "                input_image=target,\n",
    "                output_image='{1}_to_{0}_CompositeWarp.nii.gz'.format(moving_prefix, target_prefix),\n",
    "                reference_image=moving,\n",
    "                transforms=['{0}_to_{1}_affine_0GenericAffine.mat'.format(moving_prefix, target_prefix),\n",
    "                            '{0}_to_{1}_1InverseWarp.nii.gz'.format(moving_prefix, target_prefix)],\n",
    "                invert_transform_flags=[True, False], # invert the linear transform\n",
    "                dimension=3,\n",
    "                print_out_composite_warp_file=True)\n",
    "\n",
    "        backward.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3b. Apply function to generate all composite warp pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if do_parallel: # parallel\n",
    "    print('\\nCreating composite warps in parallel...')\n",
    "    start_time = time.time()\n",
    "    Parallel(n_jobs=par_numjobs, backend='multiprocessing')(\n",
    "        delayed(combine_transforms)\n",
    "        (brains[0], brains[1], names[0], names[1], subj_names, include_sub2sub=False) \n",
    "        for names, brains in zip(combinations(all_names, 2), combinations(all_brains, 2)))\n",
    "    end_time = time.time()\n",
    "    print('Composite warps computed in {0:.2f} minutes'.format((end_time - start_time)/60))\n",
    "    \n",
    "else: # serial\n",
    "    for names, brains in zip(combinations(all_names, 2), combinations(all_brains, 2)):\n",
    "\n",
    "        # make explicit here\n",
    "        mov_pref = names[0]\n",
    "        targ_pref = names[1]\n",
    "        mov_brain = brains[0]\n",
    "        targ_brain = brains[1]\n",
    "\n",
    "        print('\\nCreating composite warps between {0} and {1} templates...'.format(mov_pref, targ_pref))\n",
    "        start_time = time.time()\n",
    "        combine_transforms(mov_brain, targ_brain, mov_pref, targ_pref, subj_names, include_sub2sub=False)\n",
    "        end_time = time.time()\n",
    "        print('Composite warps computed in {0:.2f} minutes'.format((end_time - start_time)/60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3c. Visually assess the derived composite warp pairs\n",
    "* First use the derived composite warps to transform brains in both directions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for names, brains in zip(combinations(all_names, 2), combinations(all_brains, 2)):\n",
    "\n",
    "    mov_pref = names[0]\n",
    "    targ_pref = names[1]\n",
    "    mov_brain = brains[0]\n",
    "    targ_brain = brains[1]\n",
    "    \n",
    "    if os.path.isfile(os.getcwd() + '/{0}_to_{1}_CompositeWarp.nii.gz'.format(mov_pref, targ_pref)):\n",
    "        forward = ants.ApplyTransforms(\n",
    "                input_image=mov_brain,\n",
    "                output_image='{0}_in_{1}_composite.nii.gz'.format(mov_pref, targ_pref),\n",
    "                reference_image=targ_brain,\n",
    "                transforms=['{0}_to_{1}_CompositeWarp.nii.gz'.format(mov_pref, targ_pref)],\n",
    "                dimension=3)\n",
    "        backward = ants.ApplyTransforms(\n",
    "                input_image=targ_brain,\n",
    "                output_image='{1}_in_{0}_composite.nii.gz'.format(mov_pref, targ_pref),\n",
    "                reference_image=mov_brain,\n",
    "                transforms=['{1}_to_{0}_CompositeWarp.nii.gz'.format(mov_pref, targ_pref)],\n",
    "                dimension=3)\n",
    "\n",
    "        print('\\nTransforming {0} into {1} space...'.format(mov_pref, targ_pref))\n",
    "        start_time = time.time()\n",
    "        forward.run()\n",
    "        end_time = time.time()\n",
    "        print('Transformation complete in {0:.2f} minutes'.format((end_time - start_time)/60))\n",
    "        print('Transforming {1} into {0} space...'.format(mov_pref, targ_pref))\n",
    "        start_time = time.time()\n",
    "        backward.run()\n",
    "        end_time = time.time()\n",
    "        print('Transformation complete in {0:.2f} minutes'.format((end_time - start_time)/60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Then, overlay the transformed (by composite warps) brains on target brains for visualization\n",
    "* two plots (forward/backward) will be generated for each template pair"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for names, brains in zip(combinations(all_names, 2), combinations(all_brains, 2)):\n",
    "\n",
    "    mov_pref = names[0]\n",
    "    targ_pref = names[1]\n",
    "    mov_brain = brains[0]\n",
    "    targ_brain = brains[1]\n",
    "    mov_index = all_names.index(mov_pref)\n",
    "    targ_index = all_names.index(targ_pref)\n",
    "\n",
    "    if os.path.isfile(os.getcwd() + '/{0}_in_{1}_composite.nii.gz'.format(mov_pref, targ_pref)):\n",
    "        display = plot_anat(targ_brain, display_mode='ortho',\n",
    "                            title='{0} edges on {1} volume'.format(mov_pref, targ_pref),\n",
    "                            draw_cross=False, annotate=False, dim=all_contrasts[targ_index])\n",
    "        brain = '{0}_in_{1}_composite.nii.gz'.format(mov_pref, targ_pref)\n",
    "        display.add_edges(brain)\n",
    "        plt.savefig('Composite_{0}_on_{1}.png'.format(mov_pref, targ_pref))\n",
    "        plt.draw()\n",
    "\n",
    "        display = plot_anat(mov_brain, display_mode='ortho',\n",
    "                            title='{1} edges on {0} volume'.format(mov_pref, targ_pref),\n",
    "                            draw_cross=False, annotate=False, dim=all_contrasts[mov_index])\n",
    "        brain = '{1}_in_{0}_composite.nii.gz'.format(mov_pref, targ_pref)\n",
    "        display.add_edges(brain)\n",
    "        plt.savefig('Composite_{1}_on_{0}.png'.format(mov_pref, targ_pref))\n",
    "        plt.draw()\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Clean-up\n",
    "* Create subfolders for linear and non-linear warps, and for the warped volumes\n",
    "* Put all files in sensible subfolders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.makedirs(SUBJECT_path + \"warps/linear\" , exist_ok=True)\n",
    "os.makedirs(SUBJECT_path + \"warps/nonlinear\" , exist_ok=True)\n",
    "os.makedirs(SUBJECT_path + \"warps/final\" , exist_ok=True)\n",
    "\n",
    "[os.rename(os.getcwd() + '/' + f, SUBJECT_path + \"warps/linear/\" + f) for f in glob.glob('*_affine_*Affine.mat')];\n",
    "[os.rename(os.getcwd() + '/' + f, SUBJECT_path + \"warps/nonlinear/\" + f) for f in glob.glob('*_to_*_1Warp.nii.gz')];\n",
    "[os.rename(os.getcwd() + '/' + f, SUBJECT_path + \"warps/nonlinear/\" + f) for f in glob.glob('*_to_*_1InverseWarp.nii.gz')];\n",
    "[os.rename(os.getcwd() + '/' + f, SUBJECT_path + \"warps/final/\" + f) for f in glob.glob('*_to_*_CompositeWarp.nii.gz')];\n",
    "[os.remove(f) for f in glob.glob('*Affine.mat')]; # spurious linear warps (generated by nonlinear)\n",
    "\n",
    "os.makedirs(SUBJECT_path + \"warped_templates/linear\" , exist_ok=True)\n",
    "os.makedirs(SUBJECT_path + \"warped_templates/nonlinear\" , exist_ok=True)\n",
    "os.makedirs(SUBJECT_path + \"warped_templates/final\" , exist_ok=True)\n",
    "\n",
    "[os.rename(os.getcwd() + '/' +  f, SUBJECT_path + \"warped_templates/linear/\" + f) for f in glob.glob('*_in_*linear.nii.gz')];\n",
    "[os.rename(os.getcwd() + '/' +  f, SUBJECT_path + \"warped_templates/nonlinear/\" + f) for f in glob.glob('*_in_*linear+SyN.nii.gz')];\n",
    "[os.rename(os.getcwd() + '/' +  f, SUBJECT_path + \"warped_templates/final/\" + f) for f in glob.glob('*_in_*composite.nii.gz')];\n",
    "\n",
    "os.makedirs(SUBJECT_path + \"images/brains\" , exist_ok=True)\n",
    "os.makedirs(SUBJECT_path + \"images/linear_reg\" , exist_ok=True)\n",
    "os.makedirs(SUBJECT_path + \"images/nonlinear_reg\" , exist_ok=True)\n",
    "os.makedirs(SUBJECT_path + \"images/all_warp_pairs\" , exist_ok=True)\n",
    "\n",
    "[os.rename(os.getcwd() + '/' + f, SUBJECT_path + \"images/brains/\" + f) for f in glob.glob('*_brain.png')];\n",
    "[os.rename(os.getcwd() + '/' + f, SUBJECT_path + \"images/linear_reg/\" + f) for f in glob.glob('Linear*.png')];\n",
    "[os.rename(os.getcwd() + '/' + f, SUBJECT_path + \"images/nonlinear_reg/\" + f) for f in glob.glob('Nonlinear*.png')];\n",
    "[os.rename(os.getcwd() + '/' + f, SUBJECT_path + \"images/all_warp_pairs/\" + f) for f in glob.glob('Composite*.png')];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.6"
  },
  "toc-showcode": false,
  "toc-showmarkdowntxt": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}