FCD_Tuber_Workflow.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook for Same same but different: a web-based deep learning application revealed classifying features for the histopathologic distinction of cortical malformations\n",
    "\n",
    "The following notebook contains part of our pipeline to preprocess tiles with a problem specific data-augmentation, to train a VGG16 CNN with Random-Rotate-Zoom datagenerator. Predict an unseen test-set, quality check the predictions with ROC and Confusion-Matrix and last but not least visualizing the testing results in form of Guided GradCAMs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib notebook\n",
    "import sys\n",
    "sys.path.append('/FCDIIb_TSC/utils')#(\"PATH to utils folder for modified Datagenerator for Random Rotate Zoom function\")\n",
    "\n",
    "#helper functions from helper.py, LRfinder and clr_callback\n",
    "from helper import as_keras_metric, looking_at_augmentation, plot_history, plot_roc, confusion_matrix, plot_confusion_matrix\n",
    "from keras_lr_finder import LRFinder\n",
    "from clr_callback import CyclicLR\n",
    "\n",
    "import numpy as np\n",
    "from glob import glob\n",
    "from pathlib import Path\n",
    "import os\n",
    "import cv2\n",
    "from matplotlib.pyplot import imshow\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import roc_curve, roc_auc_score, auc\n",
    "import pandas as pd\n",
    "import imgaug as ia\n",
    "from imgaug import augmenters as iaa\n",
    "\n",
    "import keras\n",
    "import tensorflow as tf\n",
    "import keras.backend as K\n",
    "from keras.models import Sequential, Model\n",
    "from keras import optimizers\n",
    "from keras import backend\n",
    "from keras import applications\n",
    "from keras.layers import Lambda, GlobalMaxPooling2D, GlobalAveragePooling2D, Dense, Input, Flatten, Conv2D, MaxPooling2D, BatchNormalization, Dropout, AveragePooling2D, Concatenate\n",
    "from keras.applications.xception import preprocess_input\n",
    "from keras.callbacks import ReduceLROnPlateau, TensorBoard, EarlyStopping, ModelCheckpoint, Callback, LambdaCallback, CSVLogger\n",
    "from keras.metrics import categorical_accuracy\n",
    "from keras.utils import to_categorical\n",
    "from image.image_data_generator import ImageDataGenerator\n",
    "from keras.preprocessing.image import array_to_img, img_to_array\n",
    "from image.directory_iterator import DirectoryIterator\n",
    "from image.dataframe_iterator import DataFrameIterator\n",
    "from tensorflow.python.framework import ops\n",
    "\n",
    "#custom dataaugmentation for the FCD/TSC classification from augmentation.py\n",
    "from augmentation import aug, get_seq\n",
    "\n",
    "#livestreaming auc \n",
    "auc_roc = as_keras_metric(tf.metrics.auc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#VGG-16 with BatchNorm + 2. Denselayer\n",
    "#VGG_base_nottrainable: only toplayer is trainable\n",
    "#VGG_base_trainable: whole model trainable\n",
    "def get_model_classif_VGG_base_nottrainable():\n",
    "    base_model_VGG = applications.VGG16(weights='imagenet', include_top=False, input_shape=(IMAGE_SIZE, IMAGE_SIZE, 3))\n",
    "    print('Model loaded.')\n",
    "\n",
    "    x = base_model_VGG.output\n",
    "    out = GlobalMaxPooling2D(name='1')(x)\n",
    "    x = BatchNormalization(name='5')(out)\n",
    "    x = Dropout(0.5, name='6')(x)\n",
    "    x = Dense(512, activation='relu',name='7') (x)\n",
    "    x = BatchNormalization(name='8')(x)\n",
    "    x = Dropout(0.2, name='9')(x)\n",
    "    x = Dense(256, activation='relu',name='10') (x)\n",
    "    x = BatchNormalization(name='11')(x)\n",
    "    \n",
    "    predictions = (Dense(2, activation='softmax'))(x)\n",
    "    model_VGG = Model(inputs=base_model_VGG.input, outputs=predictions)\n",
    "\n",
    "    adam = optimizers.Adam()\n",
    "    model_VGG.compile(loss='categorical_crossentropy',\n",
    "                  optimizer=adam,\n",
    "                  metrics=[categorical_accuracy])\n",
    "\n",
    "    model_VGG.summary()\n",
    "    return model_VGG\n",
    "\n",
    "def get_model_classif_VGG_base_trainable():\n",
    "    base_model_VGG = applications.VGG16(weights='imagenet', include_top=False, input_shape=(IMAGE_SIZE, IMAGE_SIZE, 3))\n",
    "    print('Model loaded.')\n",
    "  \n",
    "    x = base_model_VGG.output\n",
    "    out = GlobalMaxPooling2D(name='1')(x)\n",
    "    x = BatchNormalization(name='5')(out)\n",
    "    x = Dropout(0.5, name='6')(x)\n",
    "    x = Dense(512, activation='relu',name='7') (x)\n",
    "    x = BatchNormalization(name='8')(x)\n",
    "    x = Dropout(0.5, name='9')(x)\n",
    "    x = Dense(256, activation='relu',name='10') (x)\n",
    "    x = BatchNormalization(name='11')(x)\n",
    "    \n",
    "    predictions = (Dense(2, activation='softmax'))(x)\n",
    "    model_VGG = Model(inputs=base_model_VGG.input, outputs=predictions)\n",
    "\n",
    "    adam = optimizers.Adam()\n",
    "    model_VGG.compile(loss='categorical_crossentropy',\n",
    "                  optimizer=adam,\n",
    "                  metrics=[categorical_accuracy, auc_roc])\n",
    "\n",
    "    model_VGG.summary()\n",
    "    return model_VGG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IMAGE_SIZE = 100 #start with 100x100 upscale in 50x50 steps (150, 200, 250, 300)\n",
    "IMAGE_CHANNELS = 3\n",
    "train_batch_size = 128 \n",
    "val_batch_size = 64\n",
    "\n",
    "model_name = 'VGG_16'\n",
    "\n",
    "#base_path: path to base_dir containing training/validation and testing tiles\n",
    "base_path = ('/data/')\n",
    "if not os.path.exists(base_path):\n",
    "    os.makedirs(base_path)\n",
    "\n",
    "train_path = base_path + '/train_dir' #insert path to training images\n",
    "val_path = base_path + '/val_dir' #insert path to validation images\n",
    "test_path = base_path + '/test_dir' #insert path to test images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Initialize preprocessing and data-generators with included Random-Rotate-Zoom function (target_size=crop_size)\n",
    "\n",
    "train_datagen_aug = ImageDataGenerator(\n",
    "                rescale=1./255,\n",
    "                preprocessing_function=aug\n",
    "                )\n",
    "\n",
    "validation_datagen_aug = ImageDataGenerator(\n",
    "                rescale=1./255\n",
    "                )\n",
    "\n",
    "\n",
    "\n",
    "train_gen = DirectoryIterator(train_path, train_datagen_aug,\n",
    "                              target_size=(IMAGE_SIZE, IMAGE_SIZE),\n",
    "                              crop_size = (IMAGE_SIZE, IMAGE_SIZE),\n",
    "                              batch_size=train_batch_size,\n",
    "                              class_mode='categorical')\n",
    "\n",
    "val_gen = DirectoryIterator(val_path,\n",
    "                            validation_datagen_aug,\n",
    "                            target_size=(IMAGE_SIZE,IMAGE_SIZE),\n",
    "                            crop_size = (IMAGE_SIZE, IMAGE_SIZE),\n",
    "                            batch_size=val_batch_size,\n",
    "                            class_mode='categorical')\n",
    "\n",
    "\n",
    "num_train_samples = train_gen.samples\n",
    "num_val_samples = val_gen.samples\n",
    "print(num_train_samples)\n",
    "print(num_val_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "#helper function to visualize a 9x9 plot \n",
    "#of different training images with our dataaugmentation and Random-Rotate-Zoom\n",
    "looking_at_augmentation(train_gen, path = base_path, batchsize=9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize our defined VGG16\n",
    "model_VGG = get_model_classif_VGG_base_nottrainable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#LR_Finder from https://github.com/surmenok/keras_lr_finder\n",
    "#initialize LR-Finder for model while only training the toplayer\n",
    "lr_gen = DirectoryIterator(train_path,\n",
    "                            train_datagen_aug,\n",
    "                            target_size = (IMAGE_SIZE,IMAGE_SIZE),\n",
    "                            crop_size = (100, 100),\n",
    "                            batch_size = 10000,\n",
    "                            class_mode = 'categorical')\n",
    "\n",
    "x_lr_train,y_lr_train = lr_gen.next()\n",
    "print(x_lr_train.shape)\n",
    "lr_finder = LRFinder(model_VGG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "lr_finder.find(x_lr_train, y_lr_train, 0.000001, 0.1, batch_size=64, epochs=4)\n",
    "lr_finder.plot_loss(n_skip_beginning=1, n_skip_end=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#cyclic learning rate from https://github.com/bckenstler/CLR\n",
    "def get_callbacks_clr(name_weights):\n",
    "    ReduceLR = ReduceLROnPlateau(monitor='val_categorical_accuracy', factor=0.5, patience=5, verbose=1, mode='auto', cooldown=5, min_lr=0.00005)\n",
    "    earlystopping = EarlyStopping(monitor='val_categorical_accuracy', min_delta=0.001, patience=10, verbose=1, mode='auto')\n",
    "    checkpoint = ModelCheckpoint(name_weights, monitor='val_categorical_accuracy', verbose=1, save_best_only=True, save_weights_only=True, mode='max')\n",
    "    clr = CyclicLR(base_lr=0.0001, max_lr=0.01,\n",
    "                        step_size=num_train_samples/train_batch_size/2)\n",
    "    csv_logger = CSVLogger(filename = base_path + \"/model_history_{}_toplayer.csv\".format(model_name), append=True)\n",
    "    return [clr, earlystopping, ReduceLR, checkpoint, csv_logger]\n",
    "\n",
    "name_weights = (base_path + \"/{}_toplayer.h5\".format(model_name))\n",
    "callbacks_list = get_callbacks_clr(name_weights = name_weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "history = model_VGG.fit_generator(generator = train_gen, \n",
    "                                  steps_per_epoch=len(train_gen), \n",
    "                                  validation_data=val_gen,\n",
    "                                  validation_steps=len(val_gen),\n",
    "                                  epochs=40, verbose=1,\n",
    "                                  use_multiprocessing = False, workers=16,\n",
    "                                  callbacks=callbacks_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plotting and saving training for VGG_16_basenottrainable\n",
    "plot_history(history, modelname = model_name, path=base_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize our defined VGG16\n",
    "model_VGG = get_model_classif_VGG_base_trainable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize LR-Finder for whole model\n",
    "lr_gen = DirectoryIterator(train_path,\n",
    "                            train_datagen_aug,\n",
    "                            target_size = (IMAGE_SIZE,IMAGE_SIZE),\n",
    "                            crop_size = (100, 100),\n",
    "                            batch_size = 10000,\n",
    "                            class_mode = 'categorical')\n",
    "\n",
    "x_lr_train,y_lr_train = lr_gen.next()\n",
    "print(x_lr_train.shape)\n",
    "lr_finder = LRFinder(model_VGG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_finder.find(x_lr_train, y_lr_train, 0.000001, 0.1, batch_size=64, epochs=4)\n",
    "lr_finder.plot_loss(n_skip_beginning=1, n_skip_end=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import clr_callback\n",
    "def get_callbacks_clr(name_weights):\n",
    "    ReduceLR = ReduceLROnPlateau(monitor='val_categorical_accuracy', factor=0.5, patience=5, verbose=1, mode='auto', cooldown=5, min_lr=0.00005)\n",
    "    earlystopping = EarlyStopping(monitor='val_categorical_accuracy', min_delta=0.001, patience=10, verbose=1, mode='auto')\n",
    "    checkpoint = ModelCheckpoint(name_weights, monitor='val_categorical_accuracy', verbose=1, save_best_only=True, save_weights_only=True, mode='max')\n",
    "    clr = CyclicLR(base_lr=0.00001, max_lr=0.0001,\n",
    "                        step_size=num_train_samples/train_batch_size/2)\n",
    "    csv_logger = CSVLogger(filename = base_path + \"/model_history_{}_all.csv\".format(model_name), append=True)\n",
    "    return [clr, earlystopping, ReduceLR, checkpoint, csv_logger]\n",
    "\n",
    "name_weights = (base_path + \"/{}_all.h5\".format(model_name))\n",
    "callbacks_list = get_callbacks_clr(name_weights = name_weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load the model with only the trained toplayer to train the whole model with a small LR\n",
    "model_VGG.load_weights(base_path + \"/{}_toplayer.h5\".format(model_name)) \n",
    "history = model_VGG.fit_generator(generator = train_gen, \n",
    "                                  steps_per_epoch=len(train_gen), \n",
    "                                  validation_data=val_gen,\n",
    "                                  validation_steps=len(val_gen),\n",
    "                                  epochs=50, verbose=1,\n",
    "                                  use_multiprocessing = False, workers=16,\n",
    "                                  callbacks=callbacks_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plotting and saving training for VGG_16_basetrainable\n",
    "plot_history(history, modelname = model_name, path=base_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing on unseen data \n",
    "The following code predicts the unseen test-set. Shows ROC and Confusion matrix as a quality check and prepares the predictions to only calculate the most confident tiles for the calculation of Guided GradCAMs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size_gradcam= 1\n",
    "datagen_gradcam = ImageDataGenerator(rescale=1./255,\n",
    "                                  rotate_random_zoom_crop=True)\n",
    "\n",
    "gradcam_generator = datagen_gradcam.flow_from_directory(test_path, \n",
    "                                                  batch_size = batch_size_gradcam, \n",
    "                                                  shuffle=False,\n",
    "                                                  target_size=(IMAGE_SIZE, IMAGE_SIZE), \n",
    "                                                  crop_size = (IMAGE_SIZE, IMAGE_SIZE),\n",
    "                                                  class_mode='categorical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "gradcam_generator.reset()\n",
    "\n",
    "model_VGG = get_model_classif_VGG_base_nottrainable()\n",
    "model_VGG.load_weights(base_path + \"/{}_all.h5\".format(model_name))\n",
    "\n",
    "predictions = model_VGG.predict_generator(gradcam_generator, \n",
    "                                      steps=len(gradcam_generator), \n",
    "                                      workers=8, \n",
    "                                      use_multiprocessing=False, \n",
    "                                      verbose=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the true labels\n",
    "y_label = gradcam_generator.classes\n",
    "print(y_label)\n",
    "\n",
    "df_preds = pd.DataFrame(predictions, columns=['FCD', 'TSC'])\n",
    "print(df_preds.head())\n",
    "\n",
    "# Get the predicted labels as probabilities\n",
    "y_pred = df_preds['TSC']\n",
    "y_predictions = np.asarray(y_pred.round(3))\n",
    "print(y_predictions)\n",
    "\n",
    "#rounding predictions + make them int --> compare to label / heatmap\n",
    "y_pred = y_pred.round()\n",
    "y_pred_int = []\n",
    "for i in y_pred:\n",
    "    if i == 1.0:\n",
    "        y_pred_int.append(1)\n",
    "    else:\n",
    "        y_pred_int.append(0)\n",
    "\n",
    "y_pred_int = np.asarray(y_pred_int)\n",
    "print((y_pred_int.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "#view AUCROC and ConfusionMatrix to Quality-Check\n",
    "plot_roc(y_label, y_pred, modelname='try', path = base_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = confusion_matrix(y_label, y_pred.round())\n",
    "cm_plot_labels = ['FCD', 'TSC']\n",
    "\n",
    "plot_confusion_matrix(cm, cm_plot_labels, title='Confusion Matrix VGG')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filenames = gradcam_generator.filenames\n",
    "incorrect = np.where(y_pred!=y_label)[0]\n",
    "for i, incorrect in enumerate(incorrect):\n",
    "    print(incorrect)\n",
    "    print(filenames[incorrect])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Guided Grad CAMs\n",
    "The following code highlights those pixels which had a positive impact on the prediction of the class of interest red, while unimportant areas are marked grey.\n",
    "We adapted the code from https://github.com/jacobgil/keras-grad-cam to implement Guided Grad CAMs in Keras. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#All defenitions and functions used to calculate Guided Grad CAMs\n",
    "\n",
    "#helper functions for visualization\n",
    "def target_category_loss(x, category_index, nb_classes):\n",
    "    return tf.multiply(x, K.one_hot([category_index], nb_classes))\n",
    "\n",
    "def target_category_loss_output_shape(input_shape):\n",
    "    return input_shape\n",
    "\n",
    "def normalize(x):\n",
    "    # utility function to normalize a tensor by its L2 norm\n",
    "    return x / (K.sqrt(K.mean(K.square(x))) + 1e-5)\n",
    "\n",
    "def load_image(path):\n",
    "    img = image.load_img(path, target_size=(IMAGE_SIZE, IMAGE_SIZE))\n",
    "    x = image.img_to_array(img)\n",
    "    x = (x - np.min(x))/np.ptp(x)\n",
    "    x = np.expand_dims(x, axis=0)\n",
    "    return x\n",
    "\n",
    "#registering guidedrelu as a tensorflow gradient\n",
    "#the modified gradient only backpropagates gradients which had a positive activation during forward pass\n",
    "def register_gradient():\n",
    "    if \"GuidedBackProp\" not in ops._gradient_registry._registry:\n",
    "        @ops.RegisterGradient(\"GuidedBackProp\")\n",
    "        def _GuidedBackProp(op, grad):\n",
    "            dtype = op.inputs[0].dtype\n",
    "            return grad * tf.cast(grad > 0., dtype) * \\\n",
    "                tf.cast(op.inputs[0] > 0., dtype)\n",
    "\n",
    "\n",
    "#Generating the Saliency Map of the Image with respect to the prediction & the given model        \n",
    "def compile_saliency_function(model, activation_layer='block5_conv3'):\n",
    "    input_img = model.input\n",
    "    layer_dict = dict([(layer.name, layer) for layer in model.layers[1:]])\n",
    "    layer_output = layer_dict[activation_layer].output\n",
    "    max_output = K.max(layer_output, axis=3)\n",
    "    #Keras function returning the saliency map given an image input\n",
    "    saliency = K.gradients(K.sum(max_output), input_img)[0]\n",
    "    return K.function([input_img, K.learning_phase()], [saliency])\n",
    "\n",
    "#to implement guided backpropagation we need to duplicate the model with modified gradients & linear activation\n",
    "#replacing the model activations with linear activations\n",
    "#linear activation ensures that when computing derivatives of a class output, features which minimize the otehr classes arent encluded\n",
    "def modify_backprop(model, name, name_weights):\n",
    "    g = tf.get_default_graph()\n",
    "    with g.gradient_override_map({'Relu': name}):\n",
    "\n",
    "        # get layers that have an activation\n",
    "        layer_dict = [layer for layer in model.layers[1:]\n",
    "                      if hasattr(layer, 'activation')]\n",
    "\n",
    "        # replace relu activation\n",
    "        for layer in layer_dict:\n",
    "            if layer.activation == keras.activations.relu:\n",
    "                layer.activation = tf.nn.relu\n",
    "\n",
    "        # re-instanciate a new model\n",
    "        new_model = get_model_classif_VGG_base_trainable()\n",
    "        new_model.load_weights(name_weights)\n",
    "    return new_model\n",
    "\n",
    "def deprocess_image(x):\n",
    "    '''\n",
    "    Same normalization as in:\n",
    "    https://github.com/fchollet/keras/blob/master/examples/conv_filter_visualization.py\n",
    "    '''\n",
    "    if np.ndim(x) > 3:\n",
    "        x = np.squeeze(x)\n",
    "    # normalize tensor: center on 0., ensure std is 0.1\n",
    "    x -= x.mean()\n",
    "    x /= (x.std() + 1e-5)\n",
    "    x *= 0.1\n",
    "\n",
    "    # clip to [0, 1]\n",
    "    x += 0.5\n",
    "    x = np.clip(x, 0, 1)\n",
    "\n",
    "    # convert to RGB array\n",
    "    x *= 255\n",
    "    if K.image_dim_ordering() == 'th':\n",
    "        x = x.transpose((1, 2, 0))\n",
    "    x = np.clip(x, 0, 255).astype('uint8')\n",
    "    return x\n",
    "\n",
    "def grad_cam(input_model, image, category_index, layer_name):\n",
    "    #defining the layer & the position of our class on the prediction vector\n",
    "    nb_classes = 2\n",
    "    target_layer = lambda x: target_category_loss(x, category_index, nb_classes)\n",
    "    x = input_model.layers[-1].output\n",
    "    x = Lambda(target_layer, output_shape=target_category_loss_output_shape)(x)\n",
    "    model = keras.models.Model(input_model.layers[0].input, x)\n",
    "    \n",
    "    #getting the output of the layers with the predefined layername\n",
    "    loss = K.sum(model.layers[-1].output)\n",
    "    conv_output = [l for l in model.layers if l.name is layer_name][0].output\n",
    "    \n",
    "    #Gradient for our class with regard to the output-feature map of the predefined layer\n",
    "    grads = normalize(K.gradients(loss, conv_output)[0])\n",
    "    #Keras function taking image input and accessing the quantities of:\n",
    "    #\"grads\" (mean intensity gradient for a specific class) & \n",
    "    #\"conv_output\" (output featuremaps)\n",
    "    gradient_function = K.function([model.layers[0].input], [conv_output, grads])\n",
    "    output, grads_val = gradient_function([image])\n",
    "    #values of grads & conv_output as numpy arrays\n",
    "    output, grads_val = output[0, :], grads_val[0, :, :, :]\n",
    "    #normalising gradients via GlobalPooling\n",
    "    weights = np.mean(grads_val, axis = (0, 1))\n",
    "    #initializing array with right shape\n",
    "    cam = np.ones(output.shape[0 : 2], dtype = np.float32)\n",
    "    #multiplying each channel in the feature-map by \"how important it is\" with respect to the given class\n",
    "    for i, w in enumerate(weights):\n",
    "        cam += w * output[:, :, i]\n",
    "    \n",
    "    #formating for visualisation purposes: image, heatmap, GRAD CAM\n",
    "    cam = cv2.resize(cam, (IMAGE_SIZE, IMAGE_SIZE))\n",
    "    #sudo Relu\n",
    "    cam = np.maximum(cam, 0)\n",
    "    #Normalizing for Visualisation\n",
    "    heatmap = cam / np.max(cam)\n",
    "    heatmap = heatmap * 255\n",
    "    \n",
    "    \n",
    "    image = image[0, :]\n",
    "    image -= np.min(image)\n",
    "    image = np.minimum(image, 255)\n",
    "    image = 255 * image / np.max(image)\n",
    "\n",
    "    cam = cv2.applyColorMap(np.uint8(heatmap), cv2.COLORMAP_JET)\n",
    "    #cam[np.where(cam < 95)] = 0\n",
    "    #superimposed_img = cv2.addWeighted(image, 0.7, cam, 0.3, 0)\n",
    "    cam = np.float32(cam) + np.float32(image)\n",
    "    superimposed_img = 255 * cam / np.max(cam)\n",
    "    #print(np.uint8(cam))\n",
    "    return np.uint8(superimposed_img), heatmap\n",
    "\n",
    "def guided_grad_cam(model, preprocessed_input, heatmap, name_weights):\n",
    "    register_gradient()\n",
    "    guided_model = modify_backprop(model, 'GuidedBackProp', name_weights=name_weights)\n",
    "    saliency_fn = compile_saliency_function(guided_model)\n",
    "    saliency = saliency_fn([preprocessed_input, 0])\n",
    "    gradcam = saliency[0] * heatmap[..., np.newaxis]\n",
    "    ggc= deprocess_image(gradcam)\n",
    "    return ggc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot most confident FCD IIb Heatmaps + GradCAMs and save in new subdir\n",
    "mostconfident = np.where((y_label==0) & (y_label==y_pred_int))[0]\n",
    "print(len(mostconfident))\n",
    "\n",
    "gradcam_entity = 'FCD'\n",
    "\n",
    "GRAD_CAM_PATH_FCD = (base_path + '/{}/'.format(model_name))\n",
    "\n",
    "if not os.path.exists(GRAD_CAM_PATH_FCD):\n",
    "    os.makedirs(GRAD_CAM_PATH_FCD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k=0\n",
    "j=4\n",
    "\n",
    "fig, ax = plt.subplots(4, 4, figsize=(12,12), squeeze=False)\n",
    "fig.subplots_adjust(hspace=0.5, wspace=0.5)\n",
    "plt.suptitle('GRAD CAM: {}'.format(name_weights), fontsize=16)\n",
    "\n",
    "for n in range(len(mostconfident)):\n",
    "    print('Getting GRAD-CAMs for FCD {} - {}'.format(k,j))\n",
    "    for i , correct in enumerate(mostconfident[k:j]):\n",
    "        \n",
    "        preprocessed_input, y  = next(gradcam_generator)\n",
    "        im = (preprocessed_input - np.min(preprocessed_input))/np.ptp(preprocessed_input)\n",
    "        im = np.squeeze(im, axis=0)\n",
    "        print(im.shape)\n",
    "        predicted_class = np.argmax(y)\n",
    "        cam, heatmap = grad_cam(model_VGG, preprocessed_input, predicted_class, \"block5_conv3\")\n",
    "        ggc = guided_grad_cam(model_VGG, preprocessed_input, heatmap, name_weights=name_weights)\n",
    "        ax[0,i].set_title(filenames[correct], fontsize=6)\n",
    "        ax[0,i].imshow(im)\n",
    "\n",
    "        ax[1,i].set_title(\"Predicted: {}, Label: {}\".format(y_predictions[correct], y_label[correct]), fontsize=8)\n",
    "        map = ax[1,i].imshow(heatmap, interpolation=\"none\", cmap = 'rainbow')\n",
    "\n",
    "        ax[2,i].imshow(cv2.cvtColor(cam, cv2.COLOR_BGR2RGB))\n",
    "        ax[3,i].imshow(cv2.cvtColor(ggc, cv2.COLOR_BGR2RGB))\n",
    "\n",
    "        ax[0,0].set_ylabel('Original')\n",
    "        ax[1,0].set_ylabel('Heatmap')\n",
    "        ax[2,0].set_ylabel('GRAD CAM')\n",
    "        ax[3,0].set_ylabel('Guided GRAD CAM')\n",
    "    k+=8\n",
    "    j+=8\n",
    "    plt.show()\n",
    "    fig.savefig(GRAD_CAM_PATH_FCD + '/Mostconfident_FCDIIb_CAMs_{}-{}.png' .format(k,j), dpi=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot most confident TSC GradCAMs and save in new subdir\n",
    "correct_tsc = np.where((y_label==1) & (y_label==y_pred_int))[0]\n",
    "print(\"Found {} confident correct TSC labels: {}\" .format(len(correct_tsc), correct_tsc))\n",
    "\n",
    "gradcam_entity = 'TSC'\n",
    "\n",
    "GRAD_CAM_PATH_TSC = (base_path + '/{}/'.format(gradcam_entity))\n",
    "if not os.path.exists(GRAD_CAM_PATH_TSC):\n",
    "    os.makedirs(GRAD_CAM_PATH_TSC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "k=0\n",
    "j=4\n",
    "\n",
    "fig, ax = plt.subplots(4, 4, figsize=(12,12), squeeze=False)\n",
    "fig.subplots_adjust(hspace=0.5, wspace=0.5)\n",
    "plt.suptitle('GRAD CAM: {}'.format(name_weights), fontsize=16)\n",
    "\n",
    "for n in range(len(correct_tsc)):\n",
    "    print('Getting GRAD-CAMs for TSC {} - {}'.format(k,j))\n",
    "    for i , correct in enumerate(correct_tsc[k:j]):\n",
    "        preprocessed_input, y  = next(gradcam_generator)\n",
    "        im = (preprocessed_input - np.min(preprocessed_input))/np.ptp(preprocessed_input)\n",
    "        im = np.squeeze(im, axis=0)\n",
    "        print(im.shape)\n",
    "        predicted_class = np.argmax(y)\n",
    "        cam, heatmap = grad_cam(model_VGG, preprocessed_input, predicted_class, \"block5_conv3\")\n",
    "        ggc = guided_grad_cam(model_VGG, preprocessed_input, heatmap, name_weights=name_weights)\n",
    "        ax[0,i].set_title(filenames[correct], fontsize=6)\n",
    "        ax[0,i].imshow(im)\n",
    "\n",
    "        ax[1,i].set_title(\"Predicted: {}, Label: {}\".format(y_predictions[correct], y_label[correct]), fontsize=8)\n",
    "        map = ax[1,i].imshow(heatmap, interpolation=\"none\", cmap = 'rainbow')\n",
    "\n",
    "        ax[2,i].imshow(cv2.cvtColor(cam, cv2.COLOR_BGR2RGB))\n",
    "        ax[3,i].imshow(cv2.cvtColor(ggc, cv2.COLOR_BGR2RGB))\n",
    "\n",
    "        ax[0,0].set_ylabel('Original')\n",
    "        ax[1,0].set_ylabel('Heatmap')\n",
    "        ax[2,0].set_ylabel('GRAD CAM')\n",
    "        ax[3,0].set_ylabel('Guided GRAD CAM')\n",
    "    k+=4\n",
    "    j+=4\n",
    "    plt.show()\n",
    "    fig.savefig(GRAD_CAM_PATH_TSC + '/TSC_Mostconfident_CAMs_{}-{}.png' .format(k,j), dpi=100)"
   ]
  },
  {
   "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.6.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}