gen_error_maps.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% Imports\n",
    "import os\n",
    "import shutil\n",
    "import nibabel as nib\n",
    "import numpy as np\n",
    "import random\n",
    "import pickle\n",
    "import matplotlib.pyplot as plt\n",
    "from glob import glob\n",
    "from scipy import ndimage\n",
    "from nilearn.image import resample_to_img, resample_img\n",
    "from nilearn.masking import compute_background_mask, compute_epi_mask\n",
    "from nilearn.plotting import plot_roi, plot_epi\n",
    "from scipy.spatial.distance import directed_hausdorff\n",
    "from nipype.algorithms.metrics import Distance\n",
    "from sklearn.metrics import roc_curve, auc, f1_score, precision_recall_curve, average_precision_score\n",
    "from scipy import interp\n",
    "from itertools import chain\n",
    "from scipy.ndimage.morphology import binary_dilation, binary_erosion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set working directory\n",
    "os.chdir('/home/uziel/DISS')\n",
    "\n",
    "# Set root of models to be post-processed\n",
    "root = \"./milestones_6\"\n",
    "\n",
    "# Set path for error maps\n",
    "#error_map_path = './data_processed/error_maps'\n",
    "error_map_path = './data_processed/error_maps_R'\n",
    "# Load best performing model variant\n",
    "model_variant = 'DM_V0_R_[0-4]' # this must be the best performing model\n",
    "tmp = model_variant.split('_')\n",
    "if len(tmp) == 3:\n",
    "    model_name = tmp[1]\n",
    "elif len(tmp) == 4:\n",
    "    model_name = tmp[1] + '_' + tmp[2]\n",
    "else:\n",
    "    model_name = tmp[1] + '_' + tmp[2] + '_' + tmp[3]\n",
    "\n",
    "# Load all trained models (k-folds) of model_variant\n",
    "trained_models = sorted(glob(os.path.join(root, model_variant)))\n",
    "\n",
    "# If data augmentation, set to (number of clones + 1), otherwise 1\n",
    "n_subject_copies = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**GET ERROR MAPS**\n",
    "\n",
    "Compute error maps for the best performing model. Get probability maps for class 0 and 1 for all validation subjects across all folds (should sum up to 43). Obtain error maps by substracting each subject's probability map from the ground truth. Save error maps. For each trained model, create files listing error maps corresponding to the training subjects (and validation?)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create folder in data for error maps\n",
    "if os.path.exists(error_map_path): shutil.rmtree(error_map_path)\n",
    "os.mkdir(error_map_path)\n",
    "\n",
    "# Get all probability maps for class 0 and 1 from model variant X (43)\n",
    "prob_maps_class_0 = []\n",
    "prob_maps_class_1 = []\n",
    "for model in trained_models:\n",
    "    root_1 = os.path.join(model, 'output/predictions/testSession/predictions')\n",
    "\n",
    "    # Load probability maps of background\n",
    "    prob_maps_class_0 += (glob(os.path.join(root_1, '*ProbMapClass0.nii.gz')))\n",
    "    # Load probability maps of foreground\n",
    "    prob_maps_class_1 += (glob(os.path.join(root_1, '*ProbMapClass1.nii.gz')))\n",
    "\n",
    "prob_maps_class_0 = sorted(prob_maps_class_0)\n",
    "prob_maps_class_1 = sorted(prob_maps_class_1)\n",
    "\n",
    "# Get all ground truth labels for all training subjects (43)\n",
    "root_2 = './data_processed/ISLES2017/training'\n",
    "subject_labels = sorted([y\n",
    "                         for x in os.walk(root_2)\n",
    "                         for y in glob(os.path.join(x[0], '*OT*.nii.gz'))\n",
    "                         if 'clone' not in y\n",
    "                        ])\n",
    "\n",
    "# Subject code in prediction files comes from MTT channel\n",
    "subject_mtt = sorted([y\n",
    "                      for x in os.walk(root_2)\n",
    "                      for y in glob(os.path.join(x[0], '*MTT*.nii.gz'))\n",
    "                      if 'clone' not in y\n",
    "                     ])\n",
    "\n",
    "# Compute error maps\n",
    "for i in range(len(subject_mtt)):\n",
    "    # Load label\n",
    "    label = nib.load(subject_labels[i])\n",
    "    \n",
    "    # Get subject code\n",
    "    code = os.path.basename(subject_mtt[i]).split('.')[-3]\n",
    "    \n",
    "    # Get probability maps of subject code\n",
    "    pmap_0 = [m for m in prob_maps_class_0 if code in m][0]\n",
    "    pmap_1 = [m for m in prob_maps_class_1 if code in m][0]\n",
    "    \n",
    "    pmap_0_img = nib.load(pmap_0)\n",
    "    pmap_1_img = nib.load(pmap_1)\n",
    "    \n",
    "    # Compute square error map\n",
    "    emap_0 = ((label.get_data() == 0).astype(int) - pmap_0_img.get_data())**2\n",
    "    emap_1 = (label.get_data() - pmap_1_img.get_data())**2\n",
    "    \n",
    "    # Normalize\n",
    "    emap_0 = (emap_0 - np.mean(emap_0))/np.std(emap_0)\n",
    "    emap_1 = (emap_1 - np.mean(emap_1))/np.std(emap_1)\n",
    "    \n",
    "    # Save error maps\n",
    "    nib.save(nib.Nifti1Image(emap_0, pmap_0_img.affine),\n",
    "             os.path.join(error_map_path, 'EMAP.0.' + code + '.nii.gz'))\n",
    "    nib.save(nib.Nifti1Image(emap_1, pmap_1_img.affine),\n",
    "             os.path.join(error_map_path, 'EMAP.1.' + code + '.nii.gz'))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**CREATE CONFIGURATION FILES FOR WEIGHTED MAPS**\n",
    "\n",
    "weightedMapsForSamplingEachCategoryTrain = [\"./weightMapsForeground.cfg\", \"./weightMapsBackground.cfg\"]\n",
    "#weightedMapsForSamplingEachCategoryVal = [\"./validation/weightMapsForeground.cfg\", \"./validation/weightMapsBackground.cfg\"]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load error map paths\n",
    "emaps_0 = sorted([x for x in os.listdir(error_map_path) if 'EMAP.0.' in x])\n",
    "emaps_1 = sorted([x for x in os.listdir(error_map_path) if 'EMAP.1.' in x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def data_to_file(data, path):\n",
    "    out = open(path, \"w\")\n",
    "    for line in data:\n",
    "        print >> out, line\n",
    "    out.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for model in trained_models:\n",
    "    # Set new model path\n",
    "    old_model_path = './ischleseg/deepmedic/versions/' + os.path.basename(model)\n",
    "    new_model_path = './ischleseg/deepmedic/versions/DM_V2_R_' + model.split('_')[-1]\n",
    "    if not os.path.exists(new_model_path): shutil.copytree(old_model_path, new_model_path)\n",
    "    \n",
    "    # Load train and validation MTT channels\n",
    "    trainChannels_MTT = [line.rstrip('\\n')\n",
    "                                for line in open(os.path.join(model,\n",
    "                                                              'configFiles/train/trainChannels_MTT.cfg'))]\n",
    "    validationChannels_MTT = [line.rstrip('\\n')\n",
    "                                     for line in open(os.path.join(model,\n",
    "                                                                   'configFiles/validation/validationChannels_MTT.cfg'))]\n",
    "    \n",
    "    # Get subject codes for train and validation\n",
    "    train_codes = [x.split('.')[-3] for x in trainChannels_MTT]\n",
    "    indexes = np.unique(train_codes, return_index=True)[1]\n",
    "    train_codes = [train_codes[index] for index in sorted(indexes)]\n",
    "    val_codes = [x.split('.')[-3] for x in validationChannels_MTT]\n",
    "    \n",
    "    # Get train weight maps\n",
    "    train_weightMapsBackground = []\n",
    "    train_weightMapsForeground = []\n",
    "    for t_code in train_codes:\n",
    "        # n_subject_copies = (clones + 1)\n",
    "        for _ in range(n_subject_copies):\n",
    "            train_weightMapsBackground.append([x for x in emaps_1 if x.split('.')[-3] in t_code][0])\n",
    "            train_weightMapsForeground.append([x for x in emaps_0 if x.split('.')[-3] in t_code][0])\n",
    "    \n",
    "    # Save train weight maps\n",
    "    data_to_file(train_weightMapsBackground,\n",
    "                 os.path.join(new_model_path, 'configFiles/train/weightMapsBackground.cfg'))\n",
    "    data_to_file(train_weightMapsForeground,\n",
    "                os.path.join(new_model_path, 'configFiles/train/weightMapsForeground.cfg'))\n",
    "    \n",
    "    \n",
    "    # Get val weight maps    \n",
    "    val_weightMapsBackground = []\n",
    "    val_weightMapsForeground = []\n",
    "    for t_code in train_codes:        \n",
    "        val_weightMapsBackground.append([x for x in emaps_1 if x.split('.')[-3] in t_code][0])\n",
    "        val_weightMapsForeground.append([x for x in emaps_0 if x.split('.')[-3] in t_code][0])\n",
    "\n",
    "    # Save train weight maps\n",
    "    data_to_file(val_weightMapsBackground,\n",
    "                os.path.join(new_model_path, 'configFiles/validation/weightMapsBackground.cfg'))\n",
    "    data_to_file(val_weightMapsForeground,\n",
    "                os.path.join(new_model_path, 'configFiles/validation/weightMapsForeground.cfg'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}