{"cells":[{"metadata":{},"cell_type":"markdown","source":"# References\n\nHeatmap code is adapted from: https://www.kaggle.com/ratthachat/aptos-augmentation-visualize-diabetic-retinopathy/\n\nPreprocessing filter from: https://github.com/btgraham/SparseConvNet/tree/kaggle_Diabetic_Retinopathy_competition\n\nIdea for this kind of visualisation from: http://gradcam.cloudcv.org/"},{"metadata":{},"cell_type":"markdown","source":"# Introduction and Explanation\n\nRoughly speaking, this kernel is a way of visualising what a trained model is paying attention to when it classifies an image. \n\nWhen a neural network is being trained, we use the chain rule to find the extent to which each of its parameters contributes to the network's error rate. \n\nHere we apply the same basic principle, with two differences:\n\n1. We are only interested in the model's final convolutional layer prior to its first fully-connected layer. That is, we are interested in the final point in the model at which information is 'spatially encoded'.\n\n2. We use the chain rule to find the extent to which each parameter in this layer contributes to the network's final prediction, rather than its error rate. \n\nThe final convolutional layer will contain multiple filters; the gradients across each of these filters are averaged and negative values clamped to zero. The resulting 'heatmap' is overlaid on the input image to input regions that the model has identified as important. "},{"metadata":{},"cell_type":"markdown","source":"# Implementation"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport random\nimport numpy as np\nimport pandas as pd\nimport tensorflow as tf\n\nimport cv2\nimport matplotlib.pyplot as plt\n\nfrom efficientnet import *\n\nimport keras.backend as K\nfrom keras import layers, models\nfrom keras.applications import DenseNet121, MobileNetV2\nfrom keras.callbacks import Callback, ModelCheckpoint\nfrom keras.models import Model, load_model, Sequential\nfrom keras.optimizers import Adam\nfrom keras.preprocessing.image import ImageDataGenerator\nfrom PIL import Image\nfrom sklearn.metrics import cohen_kappa_score\nfrom sklearn.model_selection import train_test_split   \nfrom keras.preprocessing.image import apply_channel_shift\nfrom scipy.ndimage.filters import convolve, gaussian_filter, sobel","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will define a ModelInfo class that contains each model's backbone (eg, MobileNetV2), its pretrained weights, and the name of its last (convolutional) activation layer, which is where the visualisation information comes from. This makes it easy to switch between different models for comparison. \n\nWe also load the images and csv to visualise from. "},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"K.clear_session()\n\n# Weights for each model\nMOBILENETV2_WEIGHTS = \"../input/mobilenetv2-2015-2019-messidor-full-cropped/run9.h5\"\nEFFICIENTNETB0_WEIGHTS = \"../input/effcientnetb0-weights/weights.h5\"\nEFFICIENTNETB1_WEIGHTS = \"../input/efficentnetb1-weights/weights.h5\"\n\n# Layers to get visualisation information from\nMOBILENETV2_LAYER = 'out_relu'\nEFFICIENTNETB0_LAYER = 'swish_49'\nEFFICIENTNETB1_LAYER = \"swish_118\"\n\nPATH_TO_IMAGE_ARRAY_224 = \"../input/image-resizer-224-224-crop/2019_Cropped.npy\"\nPATH_TO_IMAGE_ARRAY_240 = \"../input/image-resizer-240-240-crop/2019_Cropped.npy\"\nPATH_TO_CSV = \"../input/aptos2019-blindness-detection/train.csv\"\n\ndf = pd.read_csv(PATH_TO_CSV)\nimages_224 = np.load(PATH_TO_IMAGE_ARRAY_224)\nimages_240 = np.load(PATH_TO_IMAGE_ARRAY_240)\n\n# Backbones for each model\neffNetB0 = EfficientNetB0(\n    weights = None,\n    include_top = False,\n    input_shape = (None, None, 3)\n)\n\neffNetB1 = EfficientNetB1(\n    weights = None,\n    include_top = False,\n    input_shape = (None, None, 3)\n)\n\nmobilenet = MobileNetV2(\n    weights = None,\n    include_top = False,\n    input_shape = (None,None, 3)\n)\n\nclass ModelInfo():\n    def __init__(self, backbone, weights, last_layer, images):\n        self.backbone = backbone\n        self.weights = weights\n        self.last_layer = last_layer\n        self.images = images\n\nmobileNetInfo = ModelInfo(mobilenet, MOBILENETV2_WEIGHTS, MOBILENETV2_LAYER, images_224)\nefficientNetB0Info = ModelInfo(effNetB0, EFFICIENTNETB0_WEIGHTS, EFFICIENTNETB0_LAYER, images_224)\nefficientNetB1Info = ModelInfo(effNetB1, EFFICIENTNETB1_WEIGHTS, EFFICIENTNETB1_LAYER, images_240)\n\n# Only need to change this line\ncurrentModel = mobileNetInfo\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To get the necessary visualisation information from the model, it needs to be parsed from Sequential form to Functional form, as per the following helper method:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_visualisation_model(backbone, weights):\n    \n#     K.clear_session()\n    \n    GAP_layer = layers.GlobalAveragePooling2D()\n    drop_layer = layers.Dropout(0.5)\n    dense_layer = layers.Dense(4, activation='sigmoid')\n\n    model = Sequential()\n    model.add(backbone)\n    model.add(GAP_layer)\n    model.add(drop_layer)\n    model.add(dense_layer)\n    \n    model.load_weights(weights)\n\n    base_model = backbone\n    x = GAP_layer(base_model.layers[-1].output)\n    x = drop_layer(x)\n    final_output = dense_layer(x)\n    model = Model(base_model.layers[0].input, final_output)\n    return model\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Build the functional model:"},{"metadata":{"trusted":true},"cell_type":"code","source":"model = build_visualisation_model(\n    currentModel.backbone,\n    currentModel.weights\n)\n\nprint(model.summary())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Code for building the heatmap:"},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\nimg - the image to visualise (preprocessed)\nmodel0 - the functional model\nlayer_name - the name of the last conv layer [out_relu, I think]\nviz_img - the image that the heatmap will be overlaid on\n    (May wish to use a preprocessed image for this - seems to show up better)      \n\"\"\"\ndef gen_heatmap_img(img, model0, layer_name='last_conv_layer',viz_img=None,orig_img=None,\n    img_position=-1):\n    \n    # Promote the image to an array of length 1, then pass the array to the model to predict\n    preds_raw = model0.predict(img[np.newaxis])\n    # Convert the prediction to categorical form\n    preds = preds_raw > 0.5 \n    class_idx = (preds.astype(int).sum(axis=1) - 1)[0]\n\n    # Because we are using multilabel encoding, the output tensor will consist of multiple \n    # nodes. Note that prediction = 0 implies class_idx = -1 = 3. So the whole output layer\n    # will be examined for 0-level images. \n    class_output_tensor = model0.output[:, class_idx]\n    \n    viz_layer = model0.get_layer(layer_name)\n    # Returns the derivatives of class_output_tensor with respect to viz_layer.output\n    # Essentially, identifies the extent to which each image region at the final convolution\n    # layer contributed to the final classification. \n    grads = K.gradients(\n                        class_output_tensor ,\n                        viz_layer.output\n                        )[0] \n    \n    # Average the gradients by image region (there are multiple convolution matrices for\n    # each region)\n    pooled_grads=K.mean(grads,axis=(0,1,2))\n\n    # Return the pooled gradients and raw output for each pixel in the input image (?)\n    iterate=K.function([model0.input],[pooled_grads, viz_layer.output[0]]) \n    pooled_grad_value, viz_layer_out_value = iterate([img[np.newaxis]])\n\n    # Multiply each output value at the final convolution layer by the extent to which\n    # it contributed to the final classification\n    for i in range(pooled_grad_value.shape[0]):\n        viz_layer_out_value[:,:,i] *= pooled_grad_value[i]\n    # Average the result for each region of the image\n    heatmap = np.mean(viz_layer_out_value, axis=-1)\n    # Clamp all negative values to 0\n    heatmap = np.maximum(heatmap,0)\n    # Normalise the heatmap\n    heatmap /= np.max(heatmap)\n\n    # Standardise size of visualisation image and heatmap\n    viz_img=cv2.resize(viz_img,(img.shape[1],img.shape[0]))\n    heatmap=cv2.resize(heatmap,(viz_img.shape[1],viz_img.shape[0]))\n    # Apply preset colour map from OpenCV to the heatmap\n    # For COLORMAP_JET, low = red --> high = blue\n    heatmap_color = cv2.applyColorMap(np.uint8(heatmap*255), cv2.COLORMAP_JET)/255\n    # Overlay the visualisation image and the heatmap, with different degrees of transparency\n    heated_img = heatmap_color*0.3 + viz_img*0.5\n    \n    display_images(img, viz_img, heatmap_color, heated_img, preds_raw, img_position)\n    \n    plt.show()\n    return heated_img\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Code for displaying the heatmap:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Source: https://github.com/btgraham/SparseConvNet/tree/kaggle_Diabetic_Retinopathy_competition\ndef gaussian_preprocess(img):\n    img=cv2.addWeighted(img,4, cv2.GaussianBlur( img , (0,0) ,  10) ,-4 ,128)\n    return img/255\n\ndef display_images(img, viz_img, heatmap_color, heated_img, preds_raw, img_position):\n    fig, axs = plt.subplots(1, 4, constrained_layout=True, figsize=(10, 4))\n    \n    image_id = \"Not found\" if img_position == -1 else df.iat[img_position, 0]\n\n    preds_raw = preds_raw[0, :]\n    preds = []\n    for i, pred in enumerate(preds_raw):\n        preds.append('%.3f' % preds_raw[i])\n    preds = '[%s]' % ', '.join(map(str, preds))\n\n    predicted = sum(preds_raw > 0.5)\n\n    actual = \"NA\" if 'diagnosis' not in df.columns else df.iat[img_position, 1] \n\n    fig.suptitle(f\"Image: {image_id} \\n\\n Output: {preds} \\n\\n Predicted: {predicted}     Actual: {actual}\")\n\n    for axis in axs:\n        axis.get_xaxis().set_visible(False)\n        axis.get_yaxis().set_visible(False)\n    \n    axs[0].imshow(img)\n    axs[0].set_title(\"Original\")\n\n    axs[1].imshow(viz_img)\n    axs[1].set_title(\"Processed\")\n\n    axs[2].imshow(heatmap_color)\n    axs[2].set_title(\"Heatmap\")\n\n    axs[3].imshow(heated_img)\n    axs[3].set_title(\"Overlay\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Results"},{"metadata":{},"cell_type":"markdown","source":"## Visualising the first 5 images\n\nThe 'Processed' images are produced by creating a blurred version of the original image, then subtracting the blurred image from the original image. This technique was used by the winner of a previous diabetic retinopathy detection competition on Kaggle (https://github.com/btgraham/SparseConvNet/tree/kaggle_Diabetic_Retinopathy_competition), but we found that it did not improve our results in this competition. It is used here because it makes small exudates more visible to the naked eye. Images used by the model to make predictions are simply resized and cropped. "},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(5):\n    gen_heatmap_img(\n    currentModel.images[i, :, :, :], \n    model, \n    layer_name=currentModel.last_layer, \n    viz_img=gaussian_preprocess(currentModel.images[i, :, :, :]),\n    img_position=i\n    )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Observations (MobileNetV2): **\n\n* Small exudates are correctly identified in Image 1\n* Larger exudates and 'cotton wool spots' are correctly identified in Image 2\n* The model appears to be focusing on the shape of Image 3 - particularly the 'notch' in the top-right corner. This is not ideal, since these features have nothing to do with diabetic retinopathy. Also, it has failed to identify the microaneurysm in the bottom-centre.\n* Image 4 has not noticed anything, which is probably what we want for an image without any sign of DR. \n* I'm not sure what's being noticed in Image 5: maybe a high density of blood vessels? Nonetheless the image is still diagnosed correctly. \n"},{"metadata":{},"cell_type":"markdown","source":"**Observations (EfficientNetB0) **\n\n* Seems to do a better job of identifying small features\n* Definitely more vulnerable to overfitting 0s based on shape"},{"metadata":{},"cell_type":"markdown","source":"# General Observations\n\nThere is a persistent problem where the model learns to classify level-0 (no DR) images based on their shape, rather than the contents of the image. This is probably a consequence of bias in the training data, where most images with DR were either horizontally or vertically cropped, and most images without DR were full size. We have tried to reduce this tendency by cropping the images to a uniform standard in preprocessing, by using data augmentation and by finding additional data to train on, but none of these have fully eliminated the problem. \n\nThe higher-capacity the model, the more vulnerable it is to this kind of overfitting. "},{"metadata":{},"cell_type":"markdown","source":"# Visualisations by diagnosis"},{"metadata":{"trusted":true},"cell_type":"code","source":"diagnosis = 2\nnum_to_display = 5\n\nfor i in range(len(currentModel.images)):\n    if num_to_display == 0: break\n    if df.iat[i, 1] == diagnosis:\n        num_to_display -= 1\n        gen_heatmap_img(\n            currentModel.images[i, :, :, :], \n            model, \n            layer_name=currentModel.last_layer, \n            viz_img=gaussian_preprocess(currentModel.images[i, :, :, :]),\n            img_position=i\n        )\n        \n        ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Observations:\n\n-The 5th picture for diagnosis=3 appears to show an eye that has been treated with [laser coagulation](https://en.wikipedia.org/wiki/Laser_coagulation)"},{"metadata":{},"cell_type":"markdown","source":"# Visualisations by correctness (or lack thereof)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# 0 = correct prediction; 4 = was 0, predicted 4 (or vice versa)\ndegree_of_error = 4\nnum_to_display = 5\n\nfor i in range(len(currentModel.images)):\n    if num_to_display == 0: break\n    img = currentModel.images[i, :, :, :]\n    img = np.expand_dims(img, axis=0)\n    y_pred = (model.predict(img) > 0.5).sum(axis=1)\n    if abs(y_pred - df.iat[i, 1]) == degree_of_error: \n        num_to_display -= 1\n        gen_heatmap_img(\n            currentModel.images[i, :, :, :], \n            model, \n            layer_name=currentModel.last_layer, \n            viz_img=gaussian_preprocess(currentModel.images[i, :, :, :]),\n            img_position=i\n        )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Observations:\n\n-The model seems to struggle with identifying 'big' exudates and haemorrhages, possibly because there aren't many examples in the dataset\n-When they are identified, they're not being translated to the appropriate label. I wonder if this indicates a problem with multilabel encoding. "},{"metadata":{},"cell_type":"markdown","source":"# Visualising Random Images\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"num_to_display = 5\n\nfor i in range(num_to_display):\n    if num_to_display == 0: break\n    n = random.randrange(len(currentModel.images))\n#     img = currentModel.images[n, :, :, :]\n#     img = np.expand_dims(img, axis=0)\n    num_to_display -= 1\n    gen_heatmap_img(\n        currentModel.images[n, :, :, :], \n        model, \n        layer_name=currentModel.last_layer, \n        viz_img=gaussian_preprocess(currentModel.images[n, :, :, :]),\n        img_position=n\n    )","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}