hierarchical_covid_cxnet.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "hierarchical_covid_cxnet.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "machine_shape": "hm"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ig7Bd5HNTYdC",
        "colab_type": "text"
      },
      "source": [
        "# Hierarchical Classification with COVID-CXNet\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IDnf-zZ4DhLS",
        "colab_type": "text"
      },
      "source": [
        "## Importing libraries and run custom functions"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ep8uoYgw9yy0",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import glob, numpy as np, cv2, matplotlib.pyplot as plt, tensorflow as tf, copy\n",
        "\n",
        "from tensorflow.keras.preprocessing.image import load_img, img_to_array, ImageDataGenerator\n",
        "from tensorflow.keras import Model, callbacks as cb, optimizers\n",
        "from tensorflow.keras.models import load_model, model_from_json\n",
        "from tensorflow.keras.layers import Dense, Dropout, Lambda, Input, Conv2D, Conv2DTranspose, concatenate, MaxPooling2D\n",
        "from tensorflow.keras import Model\n",
        "\n",
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn.metrics import confusion_matrix\n",
        "from sklearn.utils import class_weight\n",
        "from sklearn.preprocessing import OneHotEncoder\n",
        "from sklearn.compose import ColumnTransformer\n",
        "\n",
        "from skimage.filters import rank\n",
        "from skimage.morphology import disk"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rQyswq-10erh",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def match_histogram(source, template):\n",
        "    \"\"\"\n",
        "    Adjust the pixel values of a gray-scale image such that its histogram\n",
        "    matches that of a target image\n",
        "    Arguments:\n",
        "    -----------\n",
        "        source: np.ndarray\n",
        "            Image to transform; the histogram is computed over the flattened\n",
        "            array\n",
        "        template: np.ndarray\n",
        "            Template image; can have different dimensions to source\n",
        "    Returns:\n",
        "    -----------\n",
        "        matched: np.ndarray\n",
        "            The transformed output image\n",
        "    \"\"\"\n",
        "\n",
        "    oldshape = source.shape\n",
        "    source = source.ravel()\n",
        "    template = template.ravel()\n",
        "\n",
        "    # get the set of unique pixel values and their corresponding indices and\n",
        "    # counts\n",
        "    s_values, bin_idx, s_counts = np.unique(source, return_inverse=True, return_counts=True)\n",
        "    t_values, t_counts = np.unique(template, return_counts=True)\n",
        "\n",
        "    # take the cumsum of the counts and normalize by the number of pixels to\n",
        "    # get the empirical cumulative distribution functions for the source and\n",
        "    # template images (maps pixel value --> quantile)\n",
        "    s_quantiles = np.cumsum(s_counts).astype(np.float64)\n",
        "    s_quantiles /= s_quantiles[-1]\n",
        "    t_quantiles = np.cumsum(t_counts).astype(np.float64)\n",
        "    t_quantiles /= t_quantiles[-1]\n",
        "\n",
        "    # interpolate linearly to find the pixel values in the template image\n",
        "    # that correspond most closely to the quantiles in the source image\n",
        "    interpolated_t_values = np.interp(s_quantiles, t_quantiles, t_values)\n",
        "    hist_matched_img = interpolated_t_values[bin_idx].reshape(oldshape)\n",
        "    result = np.array(hist_matched_img, dtype='uint8')\n",
        "\n",
        "    return result"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cFOnW0Jbg4Cm",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def subhist(image_pdf, minimum, maximum, normalize):\n",
        "    \"\"\"\n",
        "    Compute the subhistogram between [minimum, maximum] of a given histogram image_pdf\n",
        "    :param image_pdf: numpy.array\n",
        "    :param minimum: int\n",
        "    :param maximum: int\n",
        "    :param normalize: boolean\n",
        "    :return: numpy.array\n",
        "    \"\"\"\n",
        "    hi = np.zeros(shape=image_pdf.shape)\n",
        "    total = 0\n",
        "    for idx in range(minimum, maximum+1):\n",
        "        total += image_pdf[idx]\n",
        "        hi[idx] = image_pdf[idx]\n",
        "    if normalize:\n",
        "        for idx in range(minimum, maximum+1):\n",
        "            hi[idx] /= total\n",
        "    return hi\n",
        "\n",
        "\n",
        "def CDF(hist):\n",
        "    \"\"\"\n",
        "    Compute the CDF of the input histogram\n",
        "    :param hist: numpy.array()\n",
        "    :return: numpy.array()\n",
        "    \"\"\"\n",
        "    cdf = np.zeros(shape=hist.shape)\n",
        "    cdf[0] = hist[0]\n",
        "    for idx in range(1, len(hist)):\n",
        "        cdf[idx] = cdf[idx - 1] + hist[idx]\n",
        "    return cdf\n",
        "\n",
        "\n",
        "def BEASF(image, gamma):\n",
        "    \"\"\"\n",
        "    Compute the Bi-Histogram Equalization with Adaptive Sigmoid Functions algorithm (BEASF)\n",
        "    A python implementation of the original MATLAB code:\n",
        "    https://mathworks.com/matlabcentral/fileexchange/47517-beasf-image-enhancer-for-gray-scale-images\n",
        "    The algorithm is introduced by E. F. Arriaga-Garcia et al., in the research paper:\n",
        "    https://ieeexplore.ieee.org/document/6808563\n",
        "    :param image: numpy.ndarray\n",
        "    :param gamma: float [0, 1]\n",
        "    :return: numpy.ndarray\n",
        "    \"\"\"\n",
        "    m = int(np.mean(image, dtype=np.int32))\n",
        "    h = np.histogram(image, bins=256)[0] / (image.shape[0] * image.shape[1])\n",
        "    h_lower = subhist(image_pdf=h, minimum=0, maximum=m, normalize=True)\n",
        "    h_upper = subhist(image_pdf=h, minimum=m, maximum=255, normalize=True)\n",
        "\n",
        "    cdf_lower = CDF(hist=h_lower)\n",
        "    cdf_upper = CDF(hist=h_upper)\n",
        "\n",
        "    # Find x | CDF(x) = 0.5\n",
        "    half_low = 0\n",
        "    for idx in range(0, m+2):\n",
        "        if cdf_lower[idx] > 0.5:\n",
        "            half_low = idx\n",
        "            break\n",
        "    half_up = 0\n",
        "    for idx in range(m, 256):\n",
        "        if cdf_upper[idx + 1] > 0.5:\n",
        "            half_up = idx\n",
        "            break\n",
        "\n",
        "    # sigmoid CDF creation\n",
        "    tones_low = np.arange(0, m+1, 1)\n",
        "    x_low = 5.0 * (tones_low - half_low) / m  # shift & scale intensity x to place sigmoid [-2.5, 2.5]\n",
        "    s_low = 1 / (1 + np.exp(-gamma * x_low))  # lower sigmoid\n",
        "\n",
        "    tones_up = np.arange(m, 256, 1)\n",
        "    x_up = 5.0 * (tones_up - half_up) / (255 - m)  # shift & scale intensity x to place sigmoid [-2.5, 2.5]\n",
        "    s_up = 1 / (1 + np.exp(-gamma * x_up))  # upper sigmoid\n",
        "\n",
        "    mapping_vector = np.zeros(shape=(256,))\n",
        "    for idx in range(0, m+1):\n",
        "        mapping_vector[idx] = np.int32(m * s_low[idx])\n",
        "\n",
        "    minimum = mapping_vector[0]\n",
        "    maximum = mapping_vector[m]\n",
        "    for idx in range(0, m+1):\n",
        "        mapping_vector[idx] = np.int32((m / (maximum - minimum)) * (mapping_vector[idx] - minimum))\n",
        "    for idx in range(m+1, 256):\n",
        "        mapping_vector[idx] = np.int32(m + (255 - m) * s_up[idx - m - 1])\n",
        "\n",
        "    minimum = mapping_vector[m + 1]\n",
        "    maximum = mapping_vector[255]\n",
        "    for idx in range(m+1, 256):\n",
        "        mapping_vector[idx] = (255 - m) * (mapping_vector[idx] - minimum) / (maximum - minimum) + m\n",
        "\n",
        "    res = copy.deepcopy(image)\n",
        "    res[:, :] = mapping_vector[image[:, :]]\n",
        "    return res"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3lTao6TBTrSF",
        "colab_type": "text"
      },
      "source": [
        "## First level classification: pneumonia vs normal"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "R1EeVVDK-b_G",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# run directly with segmented images dataset\n",
        "X = np.load('/content/drive/My Drive/Datasets/covid_cxr_dataset/cxr_samples_segmented.npy')\n",
        "y = np.load('/content/drive/My Drive/Datasets/covid_cxr_dataset/cxr_labels_multiclass.npy')\n",
        "X = np.concatenate((X[:700], X[720:4370], X[5000:8650]), axis=0)\n",
        "y = np.concatenate((y[:700], y[720:4370], y[5000:8650]))\n",
        "# X = X[695:5000]\n",
        "# y = y[695:5000]\n",
        "\n",
        "# accumulate CAP and CP classes as pneumonia class\n",
        "y_binary = copy.copy(y)\n",
        "for label_idx in range(len(y_binary)):\n",
        "    if y_binary[label_idx] == 2:\n",
        "        y_binary[label_idx] = 1\n",
        "\n",
        "# compute num images in each class\n",
        "num_pneu_samples = 0\n",
        "num_normal_samples = 0\n",
        "for idx, img in enumerate(X):\n",
        "    if y_binary[idx] == 1:\n",
        "        num_pneu_samples += 1\n",
        "    else:\n",
        "        num_normal_samples += 1\n",
        "print('num pneumonia images:', num_pneu_samples)\n",
        "print('num normal images:', num_normal_samples)\n",
        "\n",
        "# perform histogram matching since majority of the images are from NIH-14 dataset\n",
        "base_image = X[0]\n",
        "X = np.array([match_histogram(source=image, template=base_image) for image in X])\n",
        "\n",
        "# apply image enhancements and concat with the original image\n",
        "X_beasf = np.array([BEASF(image=image, gamma=1.5) for image in X])\n",
        "X_clahe = np.array([cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8, 8)).apply(image) for image in X])\n",
        "X_clahe = np.array([np.expand_dims(a=image, axis=-1) for image in X_clahe])\n",
        "X = np.concatenate((X_beasf, X_clahe, X), axis=-1)\n",
        "\n",
        "# rescale images from [0, 255] to [0, 1]\n",
        "# X = np.array([X[idx] / 255. for idx in range(len(X))])\n",
        "\n",
        "print('number of total dataset images:', len(X))\n",
        "print('number of total dataset labels:', len(y))\n",
        "print('dataset shape:', X.shape)\n",
        "\n",
        "rnd_idx = np.random.choice(a=len(X), size=None)\n",
        "plt.imshow(X=X[rnd_idx].squeeze(), cmap='gray')\n",
        "plt.axis('off')\n",
        "plt.title(label='random image from dataset')\n",
        "plt.show()\n",
        "\n",
        "X_train, X_test, y_train, y_test = train_test_split(X, y_binary, test_size=0.2, \n",
        "                                                    random_state=18)\n",
        "\n",
        "# NOTE: sklearn train_test_split function copies the dataset, hence deleting initial data\n",
        "# variables will increase RAM space drastically (~ x2)\n",
        "del X, X_beasf, X_clahe"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rbGpi2wS_WZA",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# IMPORTANT: applying zoom and brightness at the same time will change image pixel range from\n",
        "# [0, 1] to [0, 255], hence we need to add rescaling again (despite doing it when importing\n",
        "# dataset). Otherwise, model will see all validation data (not undergone augmentation) as \n",
        "# black images giving a same prediction value for all\n",
        "augmenter = ImageDataGenerator(rotation_range=170, horizontal_flip=True, vertical_flip=True, \n",
        "                               zoom_range=[0.8, 1.5], brightness_range=[0.5, 1.3], rescale=1./255,\n",
        "                               width_shift_range=0.1, height_shift_range=0.1, fill_mode='constant')\n",
        "\n",
        "# NOTE: Keras ImageDataGenerator does not keep input image default ratio (as of v2.3.1)\n",
        "# NOTE: Keras ImageDataGenerator does not resize-and-padd when doing rotation (thus \n",
        "# loosing part of the image due to cropping while rotating)\n",
        "checkpoint = cb.ModelCheckpoint('/content/drive/My Drive/covid_model/eps={epoch:03d}_valLoss={val_loss:.4f}.hdf5',\n",
        "                                monitor='val_loss', verbose=1, save_best_only=True, mode='auto')\n",
        "cb_list = [checkpoint]"
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "959Qea5v7o1d",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# check image augmenter to see if it is working correctly and rescaling to [0, 1]\n",
        "idx = np.random.randint(low=0, high=len(X_train), size=None)\n",
        "print('max px value of input image:', np.max(X_train[idx]))\n",
        "\n",
        "augmented = augmenter.flow(X_train[idx:idx+1], batch_size=1).next()\n",
        "print('max px value of augmented image:', np.max(augmented))\n",
        "plt.subplot(1, 2, 1)\n",
        "plt.imshow(X_train[idx])\n",
        "plt.axis('off')\n",
        "plt.subplot(1, 2, 2)\n",
        "plt.imshow(augmented.squeeze())\n",
        "plt.axis('off')\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dL4rVvwQdpUh",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# prepare test-set for validating model at the end of each epoch\n",
        "# (training-set rescaling is done in the augmentation)\n",
        "X_test = np.array([X_test[idx] / 255. for idx in range(len(X_test))])"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "N3c2qTbc_wqi",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "backbone = load_model(\"/content/drive/My Drive/CheXNet_model.hdf5\")\n",
        "fc = Dense(units=10, activation='relu', name='fc')(backbone.layers[-2].output)\n",
        "do = Dropout(rate=0.2, name='dropout')(fc)\n",
        "output = Dense(units=1, activation='sigmoid', name='pred')(do)\n",
        "classifier = Model(inputs=backbone.input, outputs=output)\n",
        "classifier.trainable = True\n",
        "\n",
        "fast_adam = optimizers.Adam(learning_rate=0.0001)\n",
        "smoothened_bce = tf.losses.BinaryCrossentropy(label_smoothing=0.1)\n",
        "classifier.compile(optimizer=fast_adam, loss=smoothened_bce, metrics=['accuracy'])\n",
        "# classifier.summary()\n",
        "print('num pretrained network layers:', len(classifier.layers))\n",
        "del backbone"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "c1XyjExIGCFA",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "fine_tuning = classifier.fit(augmenter.flow(x=X_train, y=y_train, batch_size=16), \n",
        "                             steps_per_epoch=len(X_train)//16,\n",
        "                             callbacks=cb_list, epochs=20, verbose=1, \n",
        "                             validation_data=(X_test, y_test))\n",
        "\n",
        "fig = plt.figure(figsize=(16, 8))\n",
        "plt.subplot(1, 2, 1)\n",
        "plt.plot(fine_tuning.history['loss'], color='r', label='training_loss')\n",
        "plt.plot(fine_tuning.history['val_loss'], color='g', label='validation_loss')\n",
        "plt.legend()\n",
        "plt.subplot(1, 2, 2)\n",
        "plt.plot(fine_tuning.history['accuracy'], color='r', label='training_accuracy')\n",
        "plt.plot(fine_tuning.history['val_accuracy'], color='g', label='validation_accuracy')\n",
        "plt.legend()\n",
        "plt.show()\n",
        "# fig.savefig('/content/drive/My Drive/covid_model/fine_tuning.png')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "K8fKqbZ1XBNG",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "classifier.load_weights('/content/drive/My Drive/covid_model/L1_0.4425.hdf5')\n",
        "y_pred = classifier.predict(X_test)\n",
        "print('number of test-set images:', len(y_test))\n",
        "print(y_test)\n",
        "y_pred = np.round(np.reshape(a=y_pred, newshape=(1, -1)), decimals=2)[0]\n",
        "print(y_pred)\n",
        "y_pred_rnd = np.round(np.reshape(a=y_pred, newshape=(1, -1)))[0]\n",
        "cm = confusion_matrix(y_true=y_test, y_pred=y_pred_rnd)\n",
        "print('confusion matrix:')\n",
        "print(cm)\n",
        "print('test-set accuracy:', (cm[0][0] + cm[1][1])/np.sum(cm))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "37B4EuZanRPT",
        "colab_type": "text"
      },
      "source": [
        "## Second level Classification: CAP vs CP"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bRIOwgUUnWTx",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# run directly with segmented images dataset\n",
        "X = np.load('/content/drive/My Drive/Datasets/covid_cxr_dataset/cxr_samples_segmented.npy')\n",
        "y = np.load('/content/drive/My Drive/Datasets/covid_cxr_dataset/cxr_labels_multiclass.npy')\n",
        "X = np.concatenate((X[:700], X[720:4370], X[5000:8650]), axis=0)\n",
        "y = np.concatenate((y[:700], y[720:4370], y[5000:8650]))\n",
        "\n",
        "# accumulate CAP and CP classes as pneumonia class\n",
        "y_binary = copy.copy(y)\n",
        "for label_idx in range(len(y_binary)):\n",
        "    if y_binary[label_idx] == 1:\n",
        "        y_binary[label_idx] = 0\n",
        "    elif y_binary[label_idx] == 2:\n",
        "        y_binary[label_idx] = 1\n",
        "\n",
        "# compute num images in each class\n",
        "num_cap_samples = 0\n",
        "num_cp_samples = 0\n",
        "for idx, img in enumerate(X):\n",
        "    if y_binary[idx] == 0:\n",
        "        num_cap_samples += 1\n",
        "    else:\n",
        "        num_cp_samples += 1\n",
        "print('num CAP images:', num_cap_samples)\n",
        "print('num CP images:', num_cp_samples)\n",
        "\n",
        "# perform histogram matching since majority of the images are from NIH-14 dataset\n",
        "base_image = X[0]\n",
        "X = np.array([match_histogram(source=image, template=base_image) for image in X])\n",
        "\n",
        "# apply image enhancements and concat with the original image\n",
        "X_beasf = np.array([BEASF(image=image, gamma=1.5) for image in X])\n",
        "X_clahe = np.array([cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8, 8)).apply(image) for image in X])\n",
        "X_clahe = np.array([np.expand_dims(a=image, axis=-1) for image in X_clahe])\n",
        "X = np.concatenate((X_beasf, X_clahe, X), axis=-1)\n",
        "\n",
        "# rescale images from [0, 255] to [0, 1]\n",
        "# X = np.array([X[idx] / 255. for idx in range(len(X))])\n",
        "\n",
        "print('number of total dataset images:', len(X))\n",
        "print('number of total dataset labels:', len(y))\n",
        "print('dataset shape:', X.shape)\n",
        "\n",
        "rnd_idx = np.random.choice(a=len(X), size=None)\n",
        "plt.imshow(X=X[rnd_idx].squeeze(), cmap='gray')\n",
        "plt.axis('off')\n",
        "plt.title(label='random image from dataset')\n",
        "plt.show()\n",
        "\n",
        "X_train, X_test, y_train, y_test = train_test_split(X, y_binary, test_size=0.2, \n",
        "                                                    random_state=18)\n",
        "\n",
        "# NOTE: sklearn train_test_split function copies the dataset, hence deleting initial data\n",
        "# variables will increase RAM space drastically (~ x2)\n",
        "del X, X_beasf, X_clahe"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jgsmvnXEAki_",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# IMPORTANT: applying zoom and brightness at the same time will change image pixel range from\n",
        "# [0, 1] to [0, 255], hence we need to add rescaling again (despite doing it when importing\n",
        "# dataset). Otherwise, model will see all validation data (not undergone augmentation) as \n",
        "# black images giving a same prediction value for all\n",
        "augmenter = ImageDataGenerator(rotation_range=170, horizontal_flip=True, vertical_flip=True, \n",
        "                               zoom_range=[0.8, 1.5], brightness_range=[0.5, 1.3], rescale=1./255,\n",
        "                               width_shift_range=0.1, height_shift_range=0.1, fill_mode='constant')\n",
        "\n",
        "# NOTE: Keras ImageDataGenerator does not keep input image default ratio (as of v2.3.1)\n",
        "# NOTE: Keras ImageDataGenerator does not resize-and-padd when doing rotation (thus \n",
        "# loosing part of the image due to cropping while rotating)\n",
        "checkpoint = cb.ModelCheckpoint('/content/drive/My Drive/covid_model/eps={epoch:03d}_valLoss={val_loss:.4f}.hdf5',\n",
        "                                monitor='val_loss', verbose=1, save_best_only=True, mode='auto')\n",
        "cb_list = [checkpoint]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fsc6cOPwAoJH",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# prepare test-set for validating model at the end of each epoch\n",
        "# (training-set rescaling is done in the augmentation)\n",
        "X_test = np.array([X_test[idx] / 255. for idx in range(len(X_test))])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5GSmbwiqAtHO",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "outputId": "0384917d-39c7-4142-d3db-bd9225f7f5ae"
      },
      "source": [
        "backbone = load_model(\"/content/drive/My Drive/CheXNet_model.hdf5\")\n",
        "fc = Dense(units=10, activation='relu', name='fc')(backbone.layers[-2].output)\n",
        "do = Dropout(rate=0.2, name='dropout')(fc)\n",
        "output = Dense(units=1, activation='sigmoid', name='pred')(do)\n",
        "classifier = Model(inputs=backbone.input, outputs=output)\n",
        "classifier.trainable = True\n",
        "\n",
        "fast_adam = optimizers.Adam(learning_rate=0.0001)\n",
        "smoothened_bce = tf.losses.BinaryCrossentropy(label_smoothing=0.1)\n",
        "classifier.compile(optimizer=fast_adam, loss=smoothened_bce, metrics=['accuracy'])\n",
        "# classifier.summary()\n",
        "print('num pretrained network layers:', len(classifier.layers))\n",
        "del backbone"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "num pretrained network layers: 431\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Iuzw22a_Aw4Z",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "fine_tuning = classifier.fit(augmenter.flow(x=X_train, y=y_train, batch_size=16), \n",
        "                             steps_per_epoch=len(X_train)//16,\n",
        "                             callbacks=cb_list, epochs=40, verbose=1, \n",
        "                             validation_data=(X_test, y_test), initial_epoch=20)\n",
        "\n",
        "fig = plt.figure(figsize=(16, 8))\n",
        "plt.subplot(1, 2, 1)\n",
        "plt.plot(fine_tuning.history['loss'], color='r', label='training_loss')\n",
        "plt.plot(fine_tuning.history['val_loss'], color='g', label='validation_loss')\n",
        "plt.legend()\n",
        "plt.subplot(1, 2, 2)\n",
        "plt.plot(fine_tuning.history['accuracy'], color='r', label='training_accuracy')\n",
        "plt.plot(fine_tuning.history['val_accuracy'], color='g', label='validation_accuracy')\n",
        "plt.legend()\n",
        "plt.show()\n",
        "fig.savefig('/content/drive/My Drive/covid_model/fine_tuning.png')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iyN-YxSrBjwv",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# classifier.load_weights('/content/drive/My Drive/covid_model/eps=021_valLoss=0.2591.hdf5')\n",
        "y_pred = classifier.predict(X_test)\n",
        "print('number of test-set images:', len(y_test))\n",
        "print(y_test)\n",
        "y_pred = np.round(np.reshape(a=y_pred, newshape=(1, -1)), decimals=2)[0]\n",
        "print(y_pred)\n",
        "y_pred_rnd = np.round(np.reshape(a=y_pred, newshape=(1, -1)))[0]\n",
        "cm = confusion_matrix(y_true=y_test, y_pred=y_pred_rnd)\n",
        "print('confusion matrix:')\n",
        "print(cm)\n",
        "print('test-set accuracy:', (cm[0][0] + cm[1][1])/np.sum(cm))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TpLgn07gT8Kz",
        "colab_type": "text"
      },
      "source": [
        "## Heatmap generation for sample image"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Buon7zfGvU_-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class GradCAM:\n",
        "    def __init__(self, model, classIdx, layerName=None):\n",
        "        # store the model, the class index used to measure the class\n",
        "        # activation map, and the layer to be used when visualizing\n",
        "        # the class activation map\n",
        "        self.model = model\n",
        "        self.classIdx = classIdx\n",
        "        self.layerName = layerName\n",
        "        # if the layer name is None, attempt to automatically find\n",
        "        # the target output layer\n",
        "        if self.layerName is None:\n",
        "            self.layerName = self.find_target_layer()\n",
        "\n",
        "    def find_target_layer(self):\n",
        "        # attempt to find the final convolutional layer in the network\n",
        "        # by looping over the layers of the network in reverse order\n",
        "        for layer in reversed(self.model.layers):\n",
        "            # check to see if the layer has a 4D output\n",
        "            if len(layer.output_shape) == 4:\n",
        "                return layer.name\n",
        "        # otherwise, we could not find a 4D layer so the GradCAM\n",
        "        # algorithm cannot be applied\n",
        "        raise ValueError(\"Could not find 4D layer. Cannot apply GradCAM.\")\n",
        "\n",
        "    def compute_heatmap(self, image, eps=1e-8, normalize=True):\n",
        "        # construct our gradient model by supplying (1) the inputs\n",
        "        # to our pre-trained model, (2) the output of the (presumably)\n",
        "        # final 4D layer in the network, and (3) the output of the\n",
        "        # softmax activations from the model\n",
        "        gradModel = Model(inputs=[self.model.input], outputs=[self.model.get_layer(self.layerName).output,\n",
        "                                                              self.model.output])\n",
        "        # record operations for automatic differentiation\n",
        "        with tf.GradientTape() as tape:\n",
        "            # cast the image tensor to a float-32 data type, pass the\n",
        "            # image through the gradient model, and grab the loss\n",
        "            # associated with the specific class index\n",
        "            inputs = tf.cast(image, tf.float32)\n",
        "            inputs = np.expand_dims(inputs, axis=0)\n",
        "            (convOutputs, predictions) = gradModel(inputs)\n",
        "            loss = predictions[:, self.classIdx]\n",
        "        # use automatic differentiation to compute the gradients\n",
        "        grads = tape.gradient(loss, convOutputs)\n",
        "        # if gradients are too small (GradCAM is zero everywhere)\n",
        "        # equal to changing the value of 'eps' func arg\n",
        "        grads = grads / (grads.numpy().max() - grads.numpy().min())\n",
        "        # compute the guided gradients\n",
        "        castConvOutputs = tf.cast(convOutputs > 0, \"float32\")\n",
        "        castGrads = tf.cast(grads > 0, \"float32\")\n",
        "        guidedGrads = castConvOutputs * castGrads * grads\n",
        "        # the convolution and guided gradients have a batch dimension\n",
        "        # (which we don't need) so let's grab the volume itself and\n",
        "        # discard the batch\n",
        "        convOutputs = convOutputs[0]\n",
        "        guidedGrads = guidedGrads[0]\n",
        "        # compute the average of the gradient values, and using them\n",
        "        # as weights, compute the ponderation of the filters with\n",
        "        # respect to the weights\n",
        "        weights = tf.reduce_mean(guidedGrads, axis=(0, 1))\n",
        "        cam = tf.reduce_sum(tf.multiply(weights, convOutputs), axis=-1)\n",
        "        # grab the spatial dimensions of the input image and resize\n",
        "        # the output class activation map to match the input image\n",
        "        # dimensions\n",
        "        h, w = image.shape[:2]\n",
        "        heatmap = cv2.resize(cam.numpy(), (w, h))\n",
        "        print('avg heatmap value:', np.mean(heatmap))\n",
        "        print('max heatmap value:', np.max(heatmap))\n",
        "        # ignore certain values lower than a threshold to get sharper heatmaps\n",
        "        # heatmap[np.where(heatmap < 1)] = 0\n",
        "        if normalize:\n",
        "            # normalize the heatmap such that all values lie in the range\n",
        "            # [0, 1], scale the resulting values to the range [0, 255],\n",
        "            # and then convert to an unsigned 8-bit integer\n",
        "            numer = heatmap - np.min(heatmap)\n",
        "            denom = (heatmap.max() - heatmap.min()) + eps\n",
        "            heatmap = numer / denom\n",
        "            heatmap = (heatmap * 255)\n",
        "        return heatmap.astype(\"uint8\")\n",
        "\n",
        "    def overlay_heatmap(self, heatmap, image, alpha=0.5, colormap=cv2.COLORMAP_VIRIDIS):\n",
        "        # apply the supplied color map to the heatmap and then\n",
        "        # overlay the heatmap on the input image\n",
        "        heatmap = cv2.applyColorMap(heatmap, colormap)\n",
        "        # if the input image is grey-scale, convert it to 3-dim\n",
        "        if len(image.shape) == 2 or image.shape[-1] != 3:\n",
        "            image = cv2.cvtColor(src=image, code=cv2.COLOR_GRAY2RGB)\n",
        "        # if image px values are in [0, 1], upscale to [0, 255]\n",
        "        if np.max(image) <= 1.0:\n",
        "            image = image * 255.0\n",
        "        output = cv2.addWeighted(image.astype('uint8'), alpha, heatmap, 1 - alpha, 0)\n",
        "        return output"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XK4xqxeLtWIj",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# predictor = load_model('/content/drive/My Drive/COVID-CXNet_model.hdf5')\n",
        "predictor = classifier\n",
        "IMG_ID = np.random.randint(low=0, high=len(X_test), size=None)\n",
        "# IMG_ID = 841\n",
        "print('image index:', IMG_ID)\n",
        "\n",
        "\n",
        "# test_img = np.zeros(shape=(600, 600, 3))\n",
        "# test_img[140:460, 140:460] = X_test[IMG_ID]\n",
        "# test_img = cv2.resize(src=test_img, dsize=(320, 320))\n",
        "\n",
        "\n",
        "test_img = X_test[IMG_ID]\n",
        "plain_img = test_img[:, :, 2]\n",
        "temp_img = (plain_img * 255.).astype('uint8')\n",
        "\n",
        "labels = ['Normal', 'CAP', 'CP']\n",
        "pred = predictor.predict(np.expand_dims(test_img, axis=0))\n",
        "pred_label = labels[np.argmax(pred)]\n",
        "\n",
        "\n",
        "cam = GradCAM(model=predictor, classIdx=np.argmax(pred), layerName=None)\n",
        "heatmap = cam.compute_heatmap(image=test_img, normalize=True)\n",
        "overlaid_heatmap = cam.overlay_heatmap(heatmap=heatmap, image=temp_img, alpha=0.8, \n",
        "                                       colormap=cv2.COLORMAP_HSV)\n",
        "\n",
        "# overlaid_heatmap = overlaid_heatmap[75:245, 75:245, :]\n",
        "# test_img = test_img[75:245, 75:245]\n",
        "\n",
        "fig = plt.figure(num=2, figsize=(8, 8))\n",
        "plt.subplot(1, 2, 1)\n",
        "plt.imshow(overlaid_heatmap)\n",
        "plt.axis('off')\n",
        "plt.title('pred=%.4f    pred_label= %s' % (np.max(pred), pred_label))\n",
        "\n",
        "plt.subplot(1, 2, 2)\n",
        "plt.imshow(X_test[IMG_ID][:, :, 2], cmap='gray')\n",
        "plt.axis('off')\n",
        "plt.title('label=%s' % y_test[IMG_ID])\n",
        "fig.tight_layout()\n",
        "plt.show()\n",
        "# fig.savefig('/content/drive/My Drive/multi_%d.pdf' % IMG_ID, bbox_inches='tight', pad_inches=0.1)"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}