fMRI_FigGen.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate figures for fMRI report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load necessary items\n",
    "\n",
    "%pylab inline\n",
    "np.set_printoptions(precision=2, suppress=True) # Set numpy to print only 2 decimal digits for neatness\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import math\n",
    "import ast\n",
    "\n",
    "#import pylab as plt\n",
    "import nibabel as nb\n",
    "\n",
    "import os, os.path\n",
    "import json\n",
    "from os.path import join as opj\n",
    "import pprint\n",
    "\n",
    "from scipy.io.matlab import loadmat\n",
    "\n",
    "from nilearn import image as nli\n",
    "from nilearn.plotting import plot_stat_map\n",
    "from nilearn.plotting import plot_epi\n",
    "from nilearn.plotting import plot_anat\n",
    "from nilearn.plotting import plot_glass_brain\n",
    "from nilearn.plotting import plot_stat_map\n",
    "from nilearn.image import resample_to_img\n",
    "from nilearn.image import math_img\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "from IPython.display import SVG, Image, Markdown, HTML, clear_output\n",
    "\n",
    "import shutil\n",
    "\n",
    "import csv\n",
    "import operator\n",
    "\n",
    "# First level\n",
    "import pandas as pd\n",
    "from nipype.interfaces.base import Bunch\n",
    "\n",
    "import sys\n",
    "from PIL import Image as pimage\n",
    "from IPython.display import Image\n",
    "\n",
    "from cairosvg import svg2png"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define some parameters for the pipeline and the functions\n",
    "\n",
    "wDir = '/data/Patients'\n",
    "infoDir = '/data/Patients/SubjectInfoFiles'\n",
    "oDir = '/data/output'\n",
    "iDir = '/data/Patients/SubjectInfoFiles'\n",
    "patDcmRoot = '/data/Patients/tmp_dicom'\n",
    "bidsSource = '/data/BIDS_source'\n",
    "\n",
    "# For the LI curves:\n",
    "aalNii = '/opt/spm12/toolbox/aal_for_SPM12/atlas/AAL.nii'\n",
    "aalLabels = '/data/masks/ROI_MNI_V4.csv'\n",
    "aalImg = nb.load(aalNii)\n",
    "aalROIs = pd.read_csv(aalLabels, delimiter='\\t')\n",
    "\n",
    "brainRegions = ['Precentral', 'Frontal', 'Rolandic', 'Supp_Motor', 'Olfactory', 'Rectus', 'Insula', 'Cingulum', \n",
    "           'Hippocampus', 'ParaHippocampal', 'Amygdala', 'Calcarine', 'Cuneus', 'Lingual', 'Occipital', 'Fusiform', \n",
    "           'Postcentral', 'Parietal', 'SupraMarginal', 'Angular', 'Precuneus', 'Paracentral_Lobule', 'Caudate', 'Putamen',\n",
    "           'Pallidum', 'Thalamus', 'Heschl', 'Temporal', 'Cerebelum', 'Vermis']\n",
    "\n",
    "enconewDirlist = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pause():\n",
    "    programPause = input()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def printmd(string, color=None):\n",
    "    colorstr = \"<span style='color:{}'>{}</span>\".format(color, string)\n",
    "    fontstr = \"<font\"\n",
    "    display(Markdown(colorstr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select subjects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subject selection:\n",
    "def subjects4report():\n",
    "    \"\"\"\n",
    "    \n",
    "    Screens all entrys with sub-<#> in wDir for completion of datatreatment\n",
    "    \n",
    "    returns: patientStatList = [subdir, [patientinfo = dcm tags], {figures info per task}]\n",
    "        \n",
    "    \"\"\"\n",
    "    \n",
    "    # Here we define the list of subjects in wDir: subList\n",
    "    subList = []\n",
    "    for item in os.listdir(infoDir):\n",
    "        if item.endswith('_log.json'):\n",
    "            subList.append(item)\n",
    "    \n",
    "    # Create a list 'patientList' with for every patient the following elements:\n",
    "    # [[ sub-<accession#>, {dcmtags dict}, [list of 1stLevel tasks]], [idem], ....]\n",
    "    \n",
    "    patientStatList = []\n",
    "            \n",
    "    for i, sub in enumerate(subList):\n",
    "        subInfoList = []\n",
    "        subID = sub.split('_')[0]\n",
    "        jsonFile = opj(infoDir, sub)\n",
    "        with open(jsonFile, \"r\") as jFile:\n",
    "            jData = json.load(jFile)\n",
    "        if jData[\"Statistics\"] != ['None']:\n",
    "            subInfoList.append(subID)\n",
    "            subInfoList.append(jData[\"Subject DCM tags\"])\n",
    "            subInfoList.append(jData[\"Statistics\"])\n",
    "            patientStatList.append(subInfoList)\n",
    "\n",
    "    return patientStatList"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# List subjects without report\n",
    "def subjectsWithoutReport():\n",
    "    \"\"\"\n",
    "    Generates a list of subjects without report as from the /data/Patients/SubjectInfoFiles/subjects.tsv\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    subjectsWOR = []\n",
    "    \n",
    "    with open(opj(wDir, 'SubjectInfoFiles', 'subjects.tsv'), 'rt') as tsv_file:\n",
    "        reader = csv.reader(tsv_file, delimiter='\\t')\n",
    "        lines = list(reader)\n",
    "\n",
    "    figDict = {}\n",
    "    for line in lines:\n",
    "        #print('%s -> %s' %(line[2], line[-2]))\n",
    "        figDict[line[2]] = line[-2]\n",
    "\n",
    "    item = 0\n",
    "    for subject in sorted(lines):\n",
    "        subIDnum = subject[2]\n",
    "        if subject[-2] == 'N':\n",
    "            item += 1\n",
    "            print('%i: sub-%s - > %s, %s' % (item, subject[2], subject[0], subject[5]))\n",
    "            subjectsWOR.append(subject)\n",
    "    \n",
    "    return subjectsWOR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hasNumbers(inputString):\n",
    "    return bool(re.search(r'\\d', inputString))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def select_subject4figures(pList):\n",
    "    \"\"\"\n",
    "    If no subject is selected, the procedure will be aborted\n",
    "    A list of available subjects and their finished 1st-level tasks will be shown,\n",
    "    \n",
    "    : returns: subID and firstlevelList of selected\n",
    "        \n",
    "    \"\"\"\n",
    "    \n",
    "    for i, subject in enumerate(subjectsWOR):\n",
    "        subName = subjectsWOR[i][0]\n",
    "        subDir = 'sub-' + subjectsWOR[i][2]\n",
    "        subParadigms = subjectsWOR[i][5]\n",
    "        print('%i = %s (%s) with %s' %(i+1, subName, subDir, subParadigms))\n",
    "        \n",
    "    subSelectNumber = int(eval(input('Select a subject to generate figures by entering the number: ')))\n",
    "    \n",
    "    subID = 'sub-' + pList[subSelectNumber-1][2]\n",
    "    firstlevelList = ast.literal_eval(pList[subSelectNumber-1][5])\n",
    "    \n",
    "    #for i in range(len(firstlevelList)):\n",
    "    #    firstlevelList[i] = str(firstlevelList[i])\n",
    "    \n",
    "    return subID, firstlevelList"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### show_slices, plot_art"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_slices(slices):\n",
    "    \"\"\" Function to display row of image slices \"\"\"\n",
    "    fig, axes = plt.subplots(1, len(slices))\n",
    "    for i, slice in enumerate(slices):\n",
    "        axes[i].imshow(slice.T, cmap=\"gray\", origin=\"lower\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_art(task):\n",
    "    \"\"\" \n",
    "    Function to plot the SVG file to screen\n",
    "    subID, subDir, outputDir and subDatasinkDir are defined\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    fileID = subID + '_task-' + task + '_bold'\n",
    "    fName = opj(subDatasinkDir, 'preproc', 'plot.r' + fileID + '.svg')\n",
    "    if os.path.isfile(fName) is True:\n",
    "        print('Artefact anaysis result for ' + task.upper())\n",
    "        display(SVG(filename = fName))\n",
    "    else:\n",
    "        print(fName + \" doesn't exist!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot_detrended_voxel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_detrended_voxel():\n",
    "    \"\"\"\n",
    "    Procedure to select a relevant voxel to plot trendy curves,\n",
    "    followed by plot of curves in that voxel for all tasks\n",
    "    \n",
    "    subID, subDir, outputDir and subDatasinkDir are defined\n",
    "        \n",
    "    \"\"\"\n",
    "    \n",
    "    task = firstlevelList[0]\n",
    "    \n",
    "    fileID = subID + '_task-' + task + '_bold'\n",
    "    \n",
    "    output = opj(outputDir, 'work_preproc_' + task)\n",
    "    try:\n",
    "        mc = nb.load(opj(output, 'applywarp_' + task, 'ra' + fileID + '_flirt.nii'))\n",
    "    except:\n",
    "        mc = nb.load(opj(output, 'applywarp_' + task, 'r' + fileID + '_flirt.nii'))\n",
    "    \n",
    "    if os.path.isdir(opj(output, 'susan_' + task)):\n",
    "        smooth = nb.load(opj(output, 'susan_' + task,  'smooth', 'mapflow', '_smooth0', 'ra' + fileID + '_flirt_smooth.nii.gz'))\n",
    "    elif os.path.isdir(opj(output, 'smooth_' + task)):\n",
    "        try:\n",
    "            smooth = nb.load(opj(output, 'smooth_' + task, 'sra' + fileID + '_flirt.nii'))\n",
    "        except:\n",
    "            smooth = nb.load(opj(output, 'smooth_' + task, 'sr' + fileID + '_flirt.nii'))\n",
    "                  \n",
    "    detrendedFile = nb.load(opj(subDatasinkDir, 'preproc', 'detrend-' + task + '.nii.gz'))\n",
    "    detrended = detrendedFile.get_data()\n",
    "    detrShape = detrended.shape\n",
    "    midX = detrShape[0]//2\n",
    "    midY = detrShape[1]//2\n",
    "    midZ = detrShape[2]//2\n",
    "        \n",
    "    slice0 = detrended[midX, :, :, 1]\n",
    "    slice1 = detrended[:, midY, :, 1]\n",
    "    slice2 = detrended[:, :, midZ, 1]\n",
    "    show_slices([slice0, slice1, slice2])\n",
    "    plt.suptitle('Center slices for detrended image for paradigm ' + task)\n",
    "    plt.pause(0.05)\n",
    "        \n",
    "    # try to select an axial slice in the brain for displaying\n",
    "    zSlice = int(input('Select a z coordinate for displaying a brain slice:'))\n",
    "        \n",
    "    slice2 = detrended[:, :, zSlice, 1]\n",
    "    show_slices([slice0, slice1, slice2])\n",
    "\n",
    "    plt.suptitle(\"Selected axial slice.\")\n",
    "    plt.grid(True)\n",
    "    plt.pause(0.05)\n",
    "    \n",
    "    x = int(input('Select a valid x coordinate for display of time signal in detrend.nii.gz:'))\n",
    "    y = int(input('Select a valid y coordinate for display of time signal in detrend.nii.gz:'))\n",
    "    z = zSlice\n",
    "    \n",
    "    for task in firstlevelList:\n",
    "        output = opj(outputDir, 'work_preproc_' + task)\n",
    "        fileID = subID + '_task-' + task + '_bold'\n",
    "        try:\n",
    "            mc = nb.load(opj(output, 'applywarp_' + task, 'ra' + fileID + '_flirt.nii'))\n",
    "        except:\n",
    "            mc = nb.load(opj(output, 'applywarp_' + task, 'r' + fileID + '_flirt.nii'))\n",
    "                  \n",
    "        if os.path.isdir(opj(output, 'susan_' + task)):\n",
    "            smooth = nb.load(opj(output, 'susan_' + task,  'smooth', 'mapflow', '_smooth0', 'ra' + fileID + '_flirt_smooth.nii.gz'))\n",
    "        elif os.path.isdir(opj(output, 'smooth_' + task)):\n",
    "            try:\n",
    "                smooth = nb.load(opj(output, 'smooth_' + task, 'sra' + fileID + '_flirt.nii'))\n",
    "            except:\n",
    "                smooth = nb.load(opj(output, 'smooth_' + task, 'sr' + fileID + '_flirt.nii'))\n",
    "        \n",
    "        detrendedFile = nb.load(opj(subDatasinkDir, 'preproc', 'detrend-' + task + '.nii.gz'))\n",
    "        detrended = detrendedFile.get_data()\n",
    "        fig = plt.figure(figsize=(12, 4))\n",
    "        plt.suptitle('Timeline of voxel ' + '[' + str(x) + ',' + str(y) + ',' + str(z) +'] for ' + task.upper())\n",
    "        plt.plot(mc.get_data()[x, y, z, :])\n",
    "        plt.plot(smooth.get_data()[x, y, z, :])\n",
    "        plt.plot(detrendedFile.get_data()[x, y, z, :])\n",
    "        plt.legend(['motion_corrected', 'smoothed', 'detrended']);\n",
    "        plt.savefig(opj(subDatasinkDir, 'figures', 'detrended_curve_' + task + '_loc_' + str(x) + '_' + str(y) + '_' + str(z) + '.png'))\n",
    "        plt.pause(0.05)\n",
    "    \n",
    "    voxel = (x, y, z)\n",
    "    \n",
    "    return voxel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### file_len, select_task"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for finding the number of outliers, we count the number of lines\n",
    "def file_len(fName):\n",
    "    with open(fName) as f:\n",
    "        for i, l in enumerate(f):\n",
    "            pass\n",
    "    try:\n",
    "        return i + 1\n",
    "    except:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def select_task():\n",
    "    \"\"\"\n",
    "    Function to select the task for reporting SPM results\n",
    "    the variable 'firstlevelList' must be defined\n",
    "    \n",
    "    returns variable 'task'\n",
    "    \"\"\"\n",
    "    \n",
    "    print('Task(s) for evaluation: %s' % firstlevelList)\n",
    "    \n",
    "    task = input('Choose a task? ')\n",
    "    if task in firstlevelList:\n",
    "        print(task)\n",
    "    else:\n",
    "        print('Please select one in the list!')\n",
    "        raise\n",
    "        \n",
    "    return task"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### contrastList_gen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def contrastList_gen():\n",
    "    '''\n",
    "    function to define a list with all the contrasts\n",
    "    firstlevelList is defined\n",
    "        \n",
    "    '''\n",
    "\n",
    "    contrastList = []\n",
    "    subjectInfoList = []\n",
    "\n",
    "    for i, task in enumerate(firstlevelList):\n",
    "        eventFile = opj(subDir, 'func', subID + '_task-' + task + '_events.tsv')\n",
    "        trialInfo = pd.read_table(eventFile)\n",
    "        #trialInfo = pd.read_table(opj(bidsSource, task + '_seconds.tsv'))\n",
    "\n",
    "        conditions = []\n",
    "        onsets = []\n",
    "        durations = []\n",
    "\n",
    "        for group in trialInfo.groupby('trial_type'):\n",
    "            conditions.append(group[0])\n",
    "            onsets.append(list(group[1].onset))\n",
    "            durations.append(group[1].duration.tolist())\n",
    "\n",
    "        subjectInfo = [Bunch(conditions=conditions,\n",
    "                          onsets=onsets,\n",
    "                          durations=durations,\n",
    "                          )]\n",
    "\n",
    "        subjectInfoList.append('subject_info_' + task)\n",
    "        subjectInfoList[i] = subjectInfo\n",
    "\n",
    "        conditionNames = []\n",
    "        if task == 'enco_new':\n",
    "            conditionNames = [conditions[1], conditions[0]]\n",
    "        else:\n",
    "            conditionNames = [conditions[0], conditions[1]]\n",
    "\n",
    "        # Contrasts\n",
    "        cont01 = ['average', 'T', conditionNames, [1/2., 1/2.]]\n",
    "        cont02 = [conditionNames[0], 'T', conditionNames, [1, 0]]\n",
    "        cont03 = [conditionNames[1],         'T', conditionNames, [0, 1]]\n",
    "        cont04 = [conditionNames[0] + ' >> ' + conditionNames[1], 'T', conditionNames, [1, -1]]\n",
    "        cont05 = [conditionNames[1] + ' >> ' + conditionNames[0], 'T', conditionNames, [-1, 1]]\n",
    "\n",
    "        contrastList.append('contrastList_' + task)\n",
    "        contrastList[i] = [cont01, cont02, cont03, cont04, cont05]\n",
    "        \n",
    "    return contrastList"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### report_activation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_enconew_info():\n",
    "    enconewDirlist = [x for x in os.listdir(opj(subDatasinkDir, '1stLevel')) if x.startswith('enconew')]\n",
    "    return enconewDirlist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def report_activation(task):\n",
    "    \"\"\"\n",
    "    print some nice figures to aggregate in a report\n",
    "    will print automagically:\n",
    "    - the design matrix, \n",
    "    - the anatomy in NMI space, \n",
    "    - a glass brain with the activation \n",
    "        \n",
    "    \"\"\"\n",
    "    \n",
    "    #try:\n",
    "    #    fileID = subID + '_task-' + task.split('_')[0] + '_bold'\n",
    "    #    taskOrig = task.split('_')[0]\n",
    "    #except:\n",
    "    #    fileID = subID + '_task-' + task + '_bold'\n",
    "    \n",
    "    fileID = subID + '_task-' + task + '_bold'\n",
    "    \n",
    "    numOutliers = file_len(opj(subDatasinkDir, 'preproc', 'art.r' + fileID + '_outliers.txt'))\n",
    "    print('The number of outliers is: %s' % numOutliers)\n",
    "    print('For a nicer figure, define the aspect ratio for the plot of the Design Matrix')\n",
    "    print('[default = 1 (square), lower is vertical rectangle, higher is horizontal rectangle]')\n",
    "    aspectInput = input('Aspect ratio: (number between 0.5 and 2.5) default [1.0]: ')\n",
    "    if aspectInput == '':\n",
    "        aspectRatio = 1.0\n",
    "    else:\n",
    "        try:\n",
    "            aspectRatio = float(aspectInput)\n",
    "            if 0.5 <= aspectRatio <= 2.5:\n",
    "                pass\n",
    "            else:\n",
    "                print('Number out of range, taking default [1.0].')\n",
    "        except:\n",
    "            aspectRatio = 1.0\n",
    "            print('Not a number, taking default [1.0].')\n",
    "\n",
    "    # Using scipy's loadmat function we can access SPM.mat\n",
    "    spmmat = loadmat(opj(subDatasinkDir, '1stLevel', task, 'SPM.mat'),\n",
    "                     struct_as_record=False)\n",
    "\n",
    "    designMatrix = spmmat['SPM'][0][0].xX[0][0].X\n",
    "    names = [i[0] for i in spmmat['SPM'][0][0].xX[0][0].name[0]]\n",
    "\n",
    "    normedDesign = designMatrix / np.abs(designMatrix).max(axis=0)\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(int(round(aspectRatio * 8)), 8))\n",
    "    plt.imshow(normedDesign, aspect='auto', cmap='gray', interpolation='none')\n",
    "    plt.suptitle('Design matrix for ' + task.upper())\n",
    "    ax.set_ylabel('Volume id')\n",
    "    ax.set_xticks(np.arange(len(names)))\n",
    "    ax.set_xticklabels(names, rotation=90)\n",
    "    plt.show()\n",
    "    #fig.savefig(opj(subDatasinkDir, 'figures', 'DesignMatrix_' + taskOrig + '.png'));\n",
    "    fig.savefig(opj(subDatasinkDir, 'figures', 'DesignMatrix_' + task + '.png'));\n",
    "    \n",
    "    anatImg = opj(subDatasinkDir, 'preproc', 'w' + subID +'_T1w.nii')\n",
    "    #epiImg4D = opj(subDatasinkDir, 'preproc', 'wdetrend-' + taskOrig + '.nii')\n",
    "    epiImg4D = opj(subDatasinkDir, 'preproc', 'wdetrend-' + task + '.nii')\n",
    "    epiImg = nli.index_img(epiImg4D, 1)\n",
    "    \n",
    "    if 'enco' in task:\n",
    "        contrast = 'con_0004.nii'\n",
    "    else:\n",
    "        contrast = 'con_0005.nii'\n",
    "    \n",
    "    if task == 'enconew' and len(enconewDirlist) > 1:\n",
    "        for enTask in enconewDirlist:\n",
    "            conImg = opj(subDatasinkDir, '1stLevel', enTask, contrast)\n",
    "            plot_glass_brain(conImg, threshold=0.5, title = enTask.upper() + ': Activation on GlassBrain')\n",
    "            posConImg = math_img('img1*img2', img1 = conImg, img2 = math_img('img>0', img = conImg))\n",
    "            plot_stat_map(posConImg, threshold=0.25, colorbar=False, title = enTask.upper() + ' on standard brain')\n",
    "            # Save the images\n",
    "            plot_glass_brain(conImg, threshold=0.5, title = enTask.upper() + ': Activation on GlassBrain', output_file = opj(subDatasinkDir, 'figures', 'GlassBrain_' + enTask + '.png'))\n",
    "            plot_stat_map(posConImg, threshold=0.25, colorbar=False, title = enTask.upper() + ' on standard brain', output_file = opj(subDatasinkDir, 'figures',enTask.upper() + '_on_standard_brain.png'))\n",
    "    else:\n",
    "        conImg = opj(subDatasinkDir, '1stLevel', task, contrast)\n",
    "        plot_glass_brain(conImg, threshold=0.5, title = task.upper() + ': Activation on GlassBrain')\n",
    "        posConImg = math_img('img1*img2', img1 = conImg, img2 = math_img('img>0', img = conImg))\n",
    "        plot_stat_map(posConImg, threshold=0.25, colorbar=False, title = task.upper() + ' on standard brain')\n",
    "        # Save the images\n",
    "        plot_glass_brain(conImg, threshold=0.5, title = task.upper() + ': Activation on GlassBrain', output_file = opj(subDatasinkDir, 'figures', 'GlassBrain_' + task + '.png'))\n",
    "        plot_stat_map(posConImg, threshold=0.25, colorbar=False, title = task.upper() + ' on standard brain', output_file = opj(subDatasinkDir, 'figures', task.upper() + '_on_standard_brain.png'))\n",
    "    \n",
    "    #if 'enco' not in task:\n",
    "    #    conImg = opj(subDatasinkDir, '1stLevel', task, 'con_0005.nii')\n",
    "    #else:\n",
    "    #    conImg = opj(subDatasinkDir, '1stLevel', task, 'con_0004.nii')\n",
    "\n",
    "    #plot_glass_brain(conImg, threshold=0.5, title = task.upper() + ': Activation on GlassBrain')\n",
    "    \n",
    "    #posConImg = math_img('img1*img2', img1 = conImg,\n",
    "    #                                     img2 = math_img('img>0', img = conImg))\n",
    "    \n",
    "    #plot_stat_map(posConImg, threshold=0.25, colorbar=False, title = task.upper() + ' on standard brain')\n",
    "    \n",
    "    # Save the images\n",
    "    #plot_glass_brain(conImg, threshold=0.5, title = task.upper() + ': Activation on GlassBrain', output_file = opj(subDatasinkDir, 'figures', 'GlassBrain_' + task + '.png'))\n",
    "    #plot_stat_map(posConImg, threshold=0.25, colorbar=False, title = task.upper() + ' on standard brain', output_file = opj(subDatasinkDir, 'figures', task.upper() + '_on_standard_brain.png'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plotparam_spmcontrast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotparam_spmcontrast(task):\n",
    "    \"\"\"\n",
    "    Print a page of SPM overlayed on anatomy or EPI image\n",
    "    \n",
    "    returns plotParam, a dictionary with plotting parameters\n",
    "        \n",
    "    \"\"\"\n",
    "        \n",
    "    # dictionary of paramaters for plotting\n",
    "    plotParam = {}\n",
    "    \n",
    "    #try:\n",
    "    #    taskNew = task.split('_')[0]\n",
    "    #except:\n",
    "    #    taskNew = task\n",
    "        \n",
    "    taskNew = task\n",
    "    \n",
    "    # Show a list of contrasts\n",
    "    print('')\n",
    "    print('List of spm contrasts for task ' + taskNew.upper())\n",
    "    for i, item in enumerate(contrastList[firstlevelList.index(taskNew)]):\n",
    "        if item[1][-1] == 'T':\n",
    "            print('  ' + str(i+1) + '. spmT_000' + str(i+1) + '.nii = ' + str(item[:2]) + str(item[-1]))\n",
    "        else:\n",
    "            print('  ' + str(i+1) + '. spmF_000' + str(i+1) + '.nii = ' + str(item[:2]))\n",
    "            \n",
    "    spmNum = int(input('Select a contrast by choosing the corresponding number: '))\n",
    "    print('')\n",
    "    \n",
    "    if contrastList[firstlevelList.index(taskNew)][spmNum-1][1] == 'T':\n",
    "        plotParam['c2p'] = 'spmT_000' + str(spmNum)\n",
    "        thrshld = input('Lower T cut-off value: [default = 3.5]')\n",
    "        if thrshld == '':\n",
    "            plotParam['thrshld'] = 3.5\n",
    "        else:\n",
    "            try:\n",
    "                plotParam['thrshld'] = float(thrshld)\n",
    "            except:\n",
    "                print('Wrong entry, taking 3.5 for T-threshold.')\n",
    "        maxVal = input('Highest T value: [default = None]')\n",
    "        if maxVal == '':\n",
    "            plotParam['maxv'] = None\n",
    "        else:\n",
    "            try:\n",
    "                plotParam['maxv'] = int(maxVal)\n",
    "            except:\n",
    "                print('Wrong entry, taking \"None\" for maxv.')\n",
    "    else:\n",
    "        plotParam['c2p'] = 'spmF_000' + str(spmNum)\n",
    "        thrshld = input('Lower F cut-off value: [default = 20]')\n",
    "        if thrshld == '':\n",
    "            plotParam['thrshld'] = 20\n",
    "        else:\n",
    "            try:\n",
    "                plotParam['thrshld'] = float(thrshld)\n",
    "            except:\n",
    "                print('Wrong entry, taking 20 for F-threshold.')\n",
    "        maxVal = input('Highest T value: [default = None]')\n",
    "        if maxVal == '':\n",
    "            plotParam['maxv'] = 200\n",
    "        else:\n",
    "            try:\n",
    "                plotParam['maxv'] = int(maxVal)\n",
    "            except:\n",
    "                print('Wrong entry, taking 200 for maxv.')\n",
    "    \n",
    "    anatImg = opj(subDatasinkDir, 'preproc', 'w' + subID +'_T1w.nii')\n",
    "    epiImg4D = opj(subDatasinkDir, 'preproc', 'wdetrend-' + taskNew + '.nii')\n",
    "    epiImg = nli.index_img(epiImg4D, 1)\n",
    "       \n",
    "    # inquire which direction to plot, sag, cor or ax\n",
    "    print(\"\")\n",
    "    plotDirection = input('Select a plotting direction (sagittal = x, coronal = y, axial = z): ')\n",
    "    if plotDirection.lower() == 'x' or plotDirection.lower().startswith('sag'):\n",
    "        plotParam['first'] = -60\n",
    "        plotParam['last'] = 60\n",
    "        plotParam['increment'] = int(math.ceil((plotParam['last']-plotParam['first'])/float(25)))\n",
    "        plotParam['title'] = task.upper() + ' L>R'\n",
    "        print('selected sagittal')\n",
    "        plotParam['displ_mod'] = 'x'\n",
    "    elif plotDirection.lower() == 'y' or plotDirection.lower().startswith('cor'):\n",
    "        plotParam['first'] = -80\n",
    "        plotParam['last'] = 30\n",
    "        plotParam['increment'] = int(math.ceil((plotParam['last']-plotParam['first'])/float(25)))\n",
    "        plotParam['title'] = task.upper()\n",
    "        print('selected coronal')\n",
    "        plotParam['displ_mod'] = 'y'\n",
    "    elif plotDirection.lower() == 'z' or plotDirection.lower().startswith('ax'):\n",
    "        plotParam['first'] = -50\n",
    "        plotParam['last'] = 75\n",
    "        plotParam['increment'] = int(math.ceil((plotParam['last']-plotParam['first'])/float(25)))\n",
    "        plotParam['title'] = task.upper()\n",
    "        print('selected axial')\n",
    "        plotParam['displ_mod'] = 'z'\n",
    "    \n",
    "    print(opj(subDatasinkDir, '1stLevel', task, plotParam['c2p'] + '.nii'))\n",
    "    \n",
    "    #plotParam['image'] = anatimg\n",
    "    backgroundIma = input('Background anatomical or EPI [default = anatomical]?: ')\n",
    "\n",
    "    if backgroundIma.lower().startswith('e'):\n",
    "        plotParam['image'] = epiImg\n",
    "        plotParam['dimming'] = 0\n",
    "    else:\n",
    "        plotParam['image'] = anatImg\n",
    "        plotParam['dimming'] = -1\n",
    "    \n",
    "    return plotParam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### aggregate_images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def aggregate_images(imageList, plotParam):\n",
    "    \"\"\"\n",
    "    Function for stacking images vertically and saving in an aggregated image file\n",
    "    \n",
    "    saves aggr_image.<extention> in the dir of the images\n",
    "    \n",
    "    returns the name of the new image including full path\n",
    "    \n",
    "    \"\"\"\n",
    "   \n",
    "    plotContrast = plotParam['c2p']\n",
    "    \n",
    "    extention = imageList[0].split('.')[-1]\n",
    "\n",
    "    images = list(map(pimage.open, imageList))\n",
    "    widths, heights = list(zip(*(i.size for i in images)))\n",
    "    \n",
    "    maxWidth = max(widths)\n",
    "    totalHeight = sum(heights)\n",
    "    \n",
    "    newImg = pimage.new('RGB', (maxWidth, totalHeight))\n",
    "    \n",
    "    yOffset = 0\n",
    "    for im in images:\n",
    "        newImg.paste(im, (0, yOffset))\n",
    "        yOffset += im.size[1]\n",
    "        \n",
    "    newImgName = imageList[0].split('_1.' + extention)[0] + '_' + plotContrast + '.' + extention\n",
    "    \n",
    "    newImg.save(opj(figDir, newImgName))\n",
    "    \n",
    "    print('New Image %s saved.' % newImgName)\n",
    "    print('Image size = %s' % str(newImg.size))\n",
    "    \n",
    "    return newImgName"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plot_fMRI_page"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_fMRI_page(plotParam, task, print2screen='y'):\n",
    "    \"\"\"\n",
    "    Inside function to plot to screen or to file\n",
    "    \n",
    "    print2screen is 'y' or 'n' in wich case safe to file\n",
    "    \n",
    "    returns outfileList, a list with the paths to the created images\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    #plotParam['title'] = task.upper()\n",
    "    \n",
    "    try:\n",
    "        if 't1w' in plotParam['image'].lower():\n",
    "            bckgrnd = 'anat'\n",
    "    except:\n",
    "        bckgrnd = 'epi'\n",
    "        \n",
    "    if plotParam['displ_mod'] == 'x': modPlt = 'sag'\n",
    "    if plotParam['displ_mod'] == 'y': modPlt = 'cor'\n",
    "    if plotParam['displ_mod'] == 'z': modPlt = 'tra'\n",
    "        \n",
    "    outfileList = []\n",
    "     \n",
    "    index = 1\n",
    "\n",
    "    for n in range(plotParam['first'], plotParam['last'], (5*plotParam['increment'])):    \n",
    "\n",
    "        coordinates = []\n",
    "\n",
    "        for i in range(n, (n + (5*plotParam['increment'])), plotParam['increment']):\n",
    "            coordinates.append(i)\n",
    "\n",
    "        coord = tuple(coordinates)\n",
    "        \n",
    "        # threshold the contrast to only see the positive values\n",
    "        \n",
    "        img2Thrhld = opj(subDatasinkDir, '1stLevel', task, plotParam['c2p'] + '.nii')\n",
    "        posSPMT = math_img('img1*img2', img1 = img2Thrhld,\n",
    "                                         img2 = math_img('img>0', img = img2Thrhld))\n",
    "\n",
    "        if print2screen.lower().startswith('y'):\n",
    "            plot_stat_map(posSPMT, \n",
    "                            title=plotParam['title'] + ' : T>' + str(plotParam['thrshld']), vmax = plotParam['maxv'], bg_img=plotParam['image'], \n",
    "                            threshold=plotParam['thrshld'], display_mode=plotParam['displ_mod'], cut_coords=coord, colorbar=False, \n",
    "                            dim=plotParam['dimming'], output_file = None);\n",
    "            plt.show()\n",
    "        elif print2screen.lower().startswith('n'):\n",
    "            outFile = (opj(subDatasinkDir, 'figures', 'overlay_' + bckgrnd + '_' + task.upper() + '_' + modPlt + '_' + str(index) + '.png'))\n",
    "            print('File will be saved to %s' % outFile)\n",
    "            plot_stat_map(posSPMT, \n",
    "                            title=plotParam['title'] + ' : T>' + str(plotParam['thrshld']), vmax = plotParam['maxv'], bg_img=plotParam['image'], colorbar=False, \n",
    "                            threshold=plotParam['thrshld'], display_mode=plotParam['displ_mod'], cut_coords=coord, \n",
    "                            dim=plotParam['dimming'], output_file = outFile);\n",
    "            plt.show()\n",
    "            outfileList.append(outFile)\n",
    "\n",
    "        index = index + 1\n",
    "\n",
    "        n = n + (5*plotParam['increment'])\n",
    "                \n",
    "    return outfileList"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### show_spmresult"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_spmresult(task):\n",
    "    \"\"\"\n",
    "    function to print the spm on background to the output\n",
    "    with the possibility to change parameters\n",
    "    \n",
    "    input: task, should be in de firstlevelList\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    print('Generating output for %s' % task)\n",
    "    #to continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def human_size(bytes, units=[' bytes','KB','MB','GB','TB', 'PB', 'EB']):\n",
    "    \"\"\" \n",
    "    Returns a human readable string reprentation of bytes by using >>, a bitwise operator \n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    return str(bytes) + units[0] if bytes < 1024 else human_size(bytes>>10, units[1:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LI curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def selectROIs():\n",
    "    \"\"\"\n",
    "    Procedure to select brain regions from AAL\n",
    "    \n",
    "    returns ROIList\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    ROIList = []\n",
    "    \n",
    "    print(brainRegions)\n",
    "    region = input('\\nType a brain region to add to ROI: <return to stop>\\n')\n",
    "\n",
    "    while region != '':\n",
    "        if region in brainRegions and region not in ROIList:\n",
    "            ROIList.append(region)\n",
    "            region = input('\\nType a brain region to add to ROI: <return to stop>\\n')\n",
    "        else:\n",
    "            print('\\nOeps, typo or already selected. Try again!\\n')\n",
    "            region = input('\\nType a brain region to add to ROI: <return to stop>\\n')\n",
    "        \n",
    "    print(ROIList)\n",
    "            \n",
    "    return ROIList"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def createROIs(ROIList):\n",
    "    \"\"\"\n",
    "    Function to create composite left and right ROIs from the regions in ROIList\n",
    "\n",
    "    returns leftROI and rightROI\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    print('\\nCreating the composite ROIs from this list: %s\\n' %str(ROIList))\n",
    "    AAL = np.array(aalImg.dataobj)\n",
    "    leftROI = np.zeros(AAL.shape)\n",
    "    rightROI = np.zeros(AAL.shape)\n",
    "    \n",
    "    for region in ROIList:\n",
    "        leftROIs = aalROIs.loc[(aalROIs.Name.str.contains(region)) & (aalROIs.Name.str.contains('_L'))]\n",
    "        rightROIs = aalROIs.loc[(aalROIs.Name.str.contains(region)) & (aalROIs.Name.str.contains('_R'))]\n",
    "\n",
    "        leftROIList = list(leftROIs.loc[:,'Number'])\n",
    "        for element in leftROIList:\n",
    "            leftROI += (AAL == element).astype(int)\n",
    "\n",
    "        rightROIList = list(rightROIs.loc[:,'Number'])\n",
    "        for element in rightROIList:\n",
    "            rightROI += (AAL == element).astype(int)\n",
    "    \n",
    "    return leftROI, rightROI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def LI_curves(leftROI, rightROI, spmT, fN = 0.1):\n",
    "    \"\"\"\n",
    "    Procedure to create Laterality Index (LI) array's to be plotted\n",
    "    LeftROI and rightROI are the ROI's for which the LI-curves have to be calculated\n",
    "    spmT is the nifti file of the contrast\n",
    "    fN is the finese for calculating the LI series (best leave unchanged)\n",
    "    \n",
    "    returns aIndexDict, with 2 arrays with keys:\n",
    "        aIndexNum (sum of voxels above theshold)\n",
    "        aIndexWei (sum of of voxels * T-value above threshold)\n",
    "        \n",
    "    \"\"\"\n",
    "    # resample the spmT image to the AAL template (loaded at start-up)\n",
    "    resampled_spmT = resample_to_img(spmT, aalImg)\n",
    "    \n",
    "    # import image as numpy array (3D)\n",
    "    resArr = np.array(resampled_spmT.dataobj)\n",
    "    \n",
    "    # calculate the masked spmT arrays\n",
    "    leftROI_masked_spmT = leftROI * resArr\n",
    "    rightROI_masked_spmT = rightROI * resArr\n",
    "        \n",
    "    # define a dict with 2 arrays\n",
    "    aIndexDict = {}\n",
    "    aIndexDict['aIndexNum'] = []\n",
    "    aIndexDict['aIndexWei'] = []\n",
    "    \n",
    "    # define the max T-value of the spmT array\n",
    "    Tmax = max(np.max(leftROI_masked_spmT), np.max(rightROI_masked_spmT))//1.5 # not to have to many 1.0's in LI series\n",
    "    \n",
    "    # create the LI series\n",
    "    for tH in range(0, int((Tmax // 1)/fN)):\n",
    "        leftROI_thrshld = leftROI_masked_spmT[leftROI_masked_spmT>(tH*fN)]\n",
    "        rightROI_thrshld = rightROI_masked_spmT[rightROI_masked_spmT>(tH*fN)]\n",
    "\n",
    "        leftNum = leftROI_thrshld.shape[0]\n",
    "        rightNum = rightROI_thrshld.shape[0]\n",
    "\n",
    "        leftW = sum(leftROI_thrshld)\n",
    "        rightW = sum(rightROI_thrshld)\n",
    "\n",
    "\n",
    "        aIndexNum = (leftNum - rightNum) / (leftNum + rightNum)\n",
    "        aIndexW = (leftW - rightW) / (leftW + rightW)\n",
    "\n",
    "\n",
    "        aIndexDict['aIndexNum'].append(aIndexNum)\n",
    "        aIndexDict['aIndexWei'].append(aIndexW)\n",
    "        \n",
    "    return aIndexDict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fig Gen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fMRI_fig(task):\n",
    "    \"\"\"\n",
    "    Function to generate an SPM with correct T cut-off\n",
    "    Select T value, ax, cor or sag and coordinate limits\n",
    "    \n",
    "    \"\"\"\n",
    "        \n",
    "    plotParam = plotparam_spmcontrast(task)\n",
    "    outfileList = plot_fMRI_page(plotParam, task, print2screen='y')        \n",
    "    \n",
    "    answer1 = input('\\nChange the threshold? y/[n]')\n",
    "    \n",
    "    while answer1.upper().startswith('Y'):\n",
    "        newT = input('Value = ')\n",
    "        plotParam['thrshld'] = float(newT)\n",
    "        outfileList = plot_fMRI_page(plotParam, task, print2screen='y')\n",
    "        answer1 = input('\\nChange the threshold? y/[n]')\n",
    "    \n",
    "    answer2 = input('\\nChange the range of the plot? y/[n]')\n",
    "    \n",
    "    while answer2.upper().startswith('Y'):\n",
    "        first = input('New value for first slice (now %s): ' % plotParam['first'])\n",
    "        plotParam['first'] = int(first)\n",
    "        last = input('New value for last slice (now %s): ' % plotParam['last'])\n",
    "        plotParam['last'] = int(last)\n",
    "        increment = input('New value for increment (now %s): ' % plotParam['increment'])\n",
    "        plotParam['increment'] = int(increment)        \n",
    "        outfileList = plot_fMRI_page(plotParam, task, print2screen='y')\n",
    "        answer2 = input('\\nChange the range of the plot? y/[n]')\n",
    "        \n",
    "    answer3 = input('Do you want to save this result: y/[n]')\n",
    "    \n",
    "    if answer3.upper().startswith('Y'):\n",
    "        outfileList = plot_fMRI_page(plotParam, task, print2screen='n')\n",
    "        newImgName = aggregate_images(outfileList, plotParam)\n",
    "        print(newImgName)\n",
    "\n",
    "        # delete the individual overlay files\n",
    "        for ff in outfileList:\n",
    "            os.remove(ff)\n",
    "\n",
    "        try:\n",
    "            if 'T1w' in plotParam['image']:\n",
    "                title = 'Overlay of T1 anatomy and ' + plotParam['c2p'] + ' for ' + task.upper()\n",
    "                print('')\n",
    "                print('')\n",
    "                printmd(title, color = 'blue')    \n",
    "        except:\n",
    "            title = 'Overlay of EPI and ' + plotParam['c2p'] + ' for ' + task.upper()\n",
    "            print('')\n",
    "            print('')\n",
    "            printmd(title, color = 'blue')\n",
    "\n",
    "        im = pimage.open(newImgName)\n",
    "        plt.figure(figsize = (24,16))\n",
    "        plt.axis('off')\n",
    "        plt.imshow(im)\n",
    "\n",
    "    else:\n",
    "        print('Figure not saved!')\n",
    "    \n",
    "    return plotParam, outfileList"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def LI_fig(task):\n",
    "    \"\"\"\n",
    "    Function to generate activation laterality curves\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    spmT = opj(subDatasinkDir, '1stLevel', task, plotParam['c2p'] + '.nii')\n",
    "    tHold = plotParam['thrshld']\n",
    "    \n",
    "    answer1 = input('\\nLI Curves? y/[n]')\n",
    "    \n",
    "    while answer1.upper().startswith('Y'):\n",
    "        if task.startswith('enco'):\n",
    "            #ROIList = ['Hippocampus', 'ParaHippocampal', 'Fusiform', 'Amygdala']\n",
    "            ROIList = selectROIs()\n",
    "        else:\n",
    "            ROIList = selectROIs() \n",
    "        \n",
    "        leftROI, rightROI = createROIs(ROIList)\n",
    "        aIndexDict = LI_curves(leftROI, rightROI, spmT, fN=0.1)\n",
    "        \n",
    "        region = ''.join(ROIList)\n",
    "        figName = opj(subDatasinkDir, 'figures', 'LI_' + task.upper() + '_' + region + '.png')\n",
    "        \n",
    "        fN = 0.1\n",
    "        numTrace = np.asarray(aIndexDict['aIndexNum'])\n",
    "        weiTrace = np.asarray(aIndexDict['aIndexWei'])\n",
    "\n",
    "        X = [x*fN for x in range(0, len(weiTrace))]\n",
    "        \n",
    "        fig = plt.figure()\n",
    "        fig.set_size_inches(7, 5)        \n",
    "        \n",
    "        plt.plot(X, weiTrace, label = 'Weighted')\n",
    "        plt.plot(X, numTrace, label = 'Numbers')\n",
    "        plt.legend(loc='upper left')\n",
    "        xlabel('T-value')\n",
    "        ylabel('Assymmetry Index')\n",
    "        plt.title(task.upper()+': '+str(ROIList) , fontdict=None, loc='center', pad=None)\n",
    "        plt.axvline(x=tHold)\n",
    "        #plt.show()\n",
    "        plt.savefig(figName, format='png', dpi=300)\n",
    "        plt.show()\n",
    "        \n",
    "        answer2 = input('Save figure? y/[n]')\n",
    "\n",
    "        if answer2.upper().startswith('N') or answer2 == '': \n",
    "            os.remove(figName)\n",
    "            print('LI figure deleted.')\n",
    "        \n",
    "        answer1 = input('LI Curves? y/[n]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def svgtopng(task):\n",
    "    \"\"\"\n",
    "    Convert the error volume svg to png format\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    svgFile = opj(subDatasinkDir, 'preproc', 'plot.r' + subID + '_task-' + task + '_bold.svg')\n",
    "    pngFile = opj(subDatasinkDir, 'preproc', 'plot.r' + subID + '_task-' + task + '_bold.png')\n",
    "    \n",
    "    svg2png(url=svgFile, write_to=pngFile)\n",
    "    \n",
    "    print('converted %s to %s' % (svgFile, pngFile))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_subjects():\n",
    "    \"\"\"\n",
    "    function to display a tabel of the subjects in /data/Patients/SubjectInfoFiles/subjects.tsv\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "    print(' #   Subject                                  Adrema          Accession  Date Exam  Referral                                           fMRI tasks                                                       DTI     F-map     Figs    Report')\n",
    "    print('-'*220)\n",
    "    with open(opj(iDir, 'subjects.tsv'), 'rt') as tsvFile:\n",
    "        reader = csv.reader(tsvFile, delimiter='\\t')\n",
    "        #sortedLines = sorted(reader, key=operator.itemgetter(3))\n",
    "        idx = 0\n",
    "        for row in reader:\n",
    "            idx = idx + 1\n",
    "            row = [idx] + row\n",
    "            print(' {:<3} {:<40} {:<15} {:<10} {:<10} {:<50} {:<65} {:<8} {:<8} {:<8} {:<8}'.format(*row))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_subjects():\n",
    "    \"\"\"\n",
    "    Function to update certain parameters in /data/Patients/SubjectInfoFiles/subjects.tsv\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    show_subjects()\n",
    "    \n",
    "    with open(opj(wDir, 'SubjectInfoFiles', 'subjects.tsv'), 'rt') as tsv_file:\n",
    "        reader = csv.reader(tsv_file, delimiter='\\t')\n",
    "        lines = list(reader)\n",
    "\n",
    "    rows = len(lines) \n",
    "    \n",
    "    try:\n",
    "        answer = int(input('\\nWelke lijn wilt u aanpassen?:'))-1\n",
    "        print(lines[answer])\n",
    "        \n",
    "        param2Change = {'Referral':4, 'Figs':8, 'Report':9}\n",
    "\n",
    "        if answer in range(rows) and answer != 0:\n",
    "            \n",
    "            item = input('\\nWelke waarde wil je veranderen? (kies uit: %s): ' % (', '.join(str(i) for i in list(param2Change.keys()))))\n",
    "            if item in param2Change.keys():\n",
    "                position = param2Change[item]\n",
    "                print(\"\\nHuidige waarde variable '%s' = %s\" % (item, lines[answer][position]))\n",
    "                newvalue = input(\"\\nNieuwe waarde voor variabele '%s' = \" % item)\n",
    "                lines[answer][position] = newvalue\n",
    "                \n",
    "                print(lines[answer])\n",
    "                \n",
    "                with open(opj(wDir, 'SubjectInfoFiles', 'newsubjects.tsv'), 'w') as tsv_file:\n",
    "                    writer = csv.writer(tsv_file, delimiter='\\t')\n",
    "                    writer.writerows(lines)\n",
    "                \n",
    "                shutil.copy(opj(wDir, 'SubjectInfoFiles', 'subjects.tsv'), opj(wDir, 'SubjectInfoFiles', 'subjects.tsv.backup'))\n",
    "                shutil.copy(opj(wDir, 'SubjectInfoFiles', 'newsubjects.tsv'), opj(wDir, 'SubjectInfoFiles', 'subjects.tsv'))\n",
    "                os.remove(opj(wDir, 'SubjectInfoFiles', 'newsubjects.tsv'))\n",
    "                \n",
    "                clear_output()\n",
    "                \n",
    "                print('Dit is the aangepaste subjects.tsv file:\\n')\n",
    "                show_subjects()\n",
    "                \n",
    "            else:\n",
    "                print('Deze parameter bestaat niet!')\n",
    "\n",
    "        else:\n",
    "            print('Deze rij bestaat niet (of is hoofding)!')\n",
    "\n",
    "    except:\n",
    "        print('Fout! Duid een lijn aan met een nummer!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "# Start generating figures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "patientStatList = subjects4report()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subjectsWOR = subjectsWithoutReport()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subID, firstlevelList = select_subject4figures(subjectsWOR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define some useful variables for further processing\n",
    "subDir = opj(wDir, subID)\n",
    "outputDir = opj(oDir, subID)\n",
    "subDatasinkDir = opj(oDir, 'results_' + subID + '-N')\n",
    "figDir = opj(subDatasinkDir, 'figures')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a directory for the figures\n",
    "dir2Create = subDatasinkDir + '/figures'\n",
    "if os.path.isdir(dir2Create):\n",
    "    print('Directory %s exists! Nothing done.' % dir2Create)\n",
    "else:\n",
    "    os.makedirs(dir2Create)\n",
    "    print('Directory %s created!' % dir2Create)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contrastList = contrastList_gen()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "firstlevelList# = ['wgen', 'read', 'encoold', 'enconew']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the realignment plots and save copies to disc under the ../results sections per task\n",
    "\n",
    "for task in firstlevelList:\n",
    "    fileID = subID + '_task-' + task + '_bold'\n",
    "    par = np.loadtxt(opj(subDatasinkDir, 'preproc', 'rp_' + fileID + '.txt'))\n",
    "    fig, axes = plt.subplots(2, 1, figsize=(15, 5))\n",
    "    plt.suptitle('Realignment plot of ' + task.upper())\n",
    "    axes[0].set_ylabel('rotation (radians)')\n",
    "    axes[0].plot(par[0:, :3])\n",
    "    axes[1].plot(par[0:, 3:])\n",
    "    axes[1].set_xlabel('time (scans)')\n",
    "    axes[1].set_ylabel('translation (mm)');\n",
    "    plt.savefig('realignmentplot.png')\n",
    "    copy2Dir = subDatasinkDir + '/figures'\n",
    "    shutil.copy('realignmentplot.png', copy2Dir + '/realignmentplot_' + task + '.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the artefact detection timelines\n",
    "for task in firstlevelList:\n",
    "    plot_art(task)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for task in firstlevelList:\n",
    "    fileID = subID + '_task-' + task + '_bold'\n",
    "    numOutliers = file_len(opj(subDatasinkDir, 'preproc', 'art.r' + fileID + '_outliers.txt'))\n",
    "    print('Number of artefacted volumes in ' + task + ' = ' + str(numOutliers))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voxel = plot_detrended_voxel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select task"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "task = select_task()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "report_activation(task)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    svgtopng(task)\n",
    "except:\n",
    "    print('No error volumes svg file for %s' % task)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotParam, outfileList = fMRI_fig(task)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LI_fig(task)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate SPM directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Put all files necessary for SPM results into a directory SPM under <oDir>/results_sub_ID/figures\n",
    "SPMdir = opj(subDatasinkDir, 'figures', 'SPM')\n",
    "\n",
    "if os.path.isdir(SPMdir):\n",
    "    print('Directory %s exists! Nothing done.' % SPMdir)\n",
    "else:\n",
    "    os.makedirs(SPMdir)\n",
    "    print('Directory %s created!' % SPMdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# copy the wdetrend en wanat files\n",
    "for f in os.listdir(opj(subDatasinkDir, 'preproc')):\n",
    "    if f.startswith('w'):\n",
    "        shutil.copy(opj(subDatasinkDir, 'preproc', f), opj(SPMdir, f))\n",
    "\n",
    "# copy the SPM12 files\n",
    "for task in firstlevelList:\n",
    "    #print(task)\n",
    "    # copy the task SPM.mat files and rename\n",
    "    shutil.copy(opj(subDatasinkDir, '1stLevel', task, 'SPM.mat'), opj(SPMdir, 'SPM_' + task + '.mat'))\n",
    "\n",
    "    # copy the spm_#####.nii files corresponding to the saved figures and rename\n",
    "    for f in os.listdir(opj(subDatasinkDir, 'figures')):\n",
    "        if task.upper() in f and f.startswith('overlay_'):\n",
    "            nbr = f.split('.')[0].split('_')[-1]\n",
    "            cntrst = f.split('.')[0].split('_')[-2]\n",
    "            shutil.copy(opj(subDatasinkDir, '1stLevel', task, cntrst + '_' + nbr + '.nii'), opj(SPMdir, cntrst + '_' + nbr + '_' + task + '.nii'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Files in SPM dir: ')\n",
    "for i, f in enumerate(os.listdir(SPMdir)):\n",
    "    fbytes = os.path.getsize(opj(SPMdir, f))\n",
    "    fsize = human_size(fbytes, units=[' bytes','KB','MB','GB','TB', 'PB', 'EB'])\n",
    "    print('%i. %s - %s' % (i+1, f, fsize))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Update subjects.tsv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "update_subjects()"
   ]
  },
  {
   "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.12"
  },
  "toc-autonumbering": true,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}