{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport sys\nimport warnings\nwarnings.simplefilter(action='ignore', category=FutureWarning)\n\n# Repository source: https://github.com/qubvel/efficientnet\nsys.path.append(os.path.abspath('../input/efficientnet/efficientnet-master/efficientnet-master/'))\nfrom efficientnet import EfficientNetB5","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-06-07T22:22:17.626380Z","iopub.execute_input":"2021-06-07T22:22:17.626673Z","iopub.status.idle":"2021-06-07T22:22:20.399278Z","shell.execute_reply.started":"2021-06-07T22:22:17.626626Z","shell.execute_reply":"2021-06-07T22:22:20.398469Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Standard dependencies\nimport cv2\nimport time\nimport scipy as sp\nimport numpy as np\nimport random as rn\nimport pandas as pd\nfrom tqdm import tqdm\nfrom PIL import Image\nfrom functools import partial\nimport matplotlib.pyplot as plt\n\n# Machine Learning\nimport tensorflow as tf\nimport keras\nfrom keras import initializers\nfrom keras import regularizers\nfrom keras import constraints\nfrom keras import backend as K\nfrom keras.activations import elu\nfrom keras.optimizers import Adam\nfrom keras.models import Sequential\nfrom keras.engine import Layer, InputSpec\nfrom keras.utils.generic_utils import get_custom_objects\nfrom keras.callbacks import Callback, EarlyStopping, ReduceLROnPlateau\nfrom keras.layers import Dense, Conv2D, Flatten, GlobalAveragePooling2D, Dropout\nfrom keras.preprocessing.image import ImageDataGenerator\nfrom sklearn.metrics import cohen_kappa_score\n\n# Path specifications\nKAGGLE_DIR = '../input/aptos2019-blindness-detection/'\nTRAIN_DF_PATH = KAGGLE_DIR + \"train.csv\"\nTRAIN_IMG_PATH = KAGGLE_DIR + \"train_images/\"\n\n# Set seed for reproducability\nseed = 1234\nrn.seed(seed)\nnp.random.seed(seed)\ntf.set_random_seed(seed)\nos.environ['PYTHONHASHSEED'] = str(seed)","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-06-07T22:22:20.860516Z","iopub.execute_input":"2021-06-07T22:22:20.860805Z","iopub.status.idle":"2021-06-07T22:22:21.383468Z","shell.execute_reply.started":"2021-06-07T22:22:20.860756Z","shell.execute_reply":"2021-06-07T22:22:21.382643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preparation <a id=\"2\"></a>","metadata":{}},{"cell_type":"markdown","source":"By examining the data we can readily see that we do not have that much data (± 700 samples per class). It is probably a good idea to use data augmentation to increase robustness of our model (See the modeling section).\n\nWe could also try to use additional data from previous competitions to increase performance. Although I do not implement this in the kernel, feel free to experiment with adding data. Additional data can be found in [this Kaggle dataset](https://www.kaggle.com/benjaminwarner/resized-2015-2019-blindness-detection-images) (± 35000 additional images).","metadata":{}},{"cell_type":"code","source":"print(\"Image IDs and Labels (TRAIN)\")\ntrain_df = pd.read_csv(TRAIN_DF_PATH)\n\n# Add extension to id_code\ntrain_df['id_code'] = train_df['id_code'] + \".png\"\nprint(f\"Training images: {train_df.shape[0]}\")\ndisplay(train_df.head())","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-06-07T22:22:23.944751Z","iopub.execute_input":"2021-06-07T22:22:23.945053Z","iopub.status.idle":"2021-06-07T22:22:23.999682Z","shell.execute_reply.started":"2021-06-07T22:22:23.945004Z","shell.execute_reply":"2021-06-07T22:22:23.999082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Specify image size\nIMG_WIDTH = 224 # 456\nIMG_HEIGHT = 224 # 456\nCHANNELS = 3","metadata":{"execution":{"iopub.status.busy":"2021-06-07T22:22:25.792900Z","iopub.execute_input":"2021-06-07T22:22:25.793197Z","iopub.status.idle":"2021-06-07T22:22:25.797226Z","shell.execute_reply.started":"2021-06-07T22:22:25.793146Z","shell.execute_reply":"2021-06-07T22:22:25.796459Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## EDA (Exploratory Data Analysis) <a id=\"4\"></a>","metadata":{}},{"cell_type":"markdown","source":"For EDA on image datasets I think one should at least examine the label distribution, the images before preprocessing and the images after preprocessing. Through examining these three aspects we can get a good sense of the problem. Note that the distribution on the test set can still vary wildly from the training data.","metadata":{}},{"cell_type":"code","source":"# Label distribution\ntrain_df['diagnosis'].value_counts().sort_index().plot(kind=\"bar\", \n                                                       figsize=(12,5), \n                                                       rot=0)\nplt.title(\"Label Distribution (Training Set)\", \n          weight='bold', \n          fontsize=18)\nplt.xticks(fontsize=15)\nplt.yticks(fontsize=15)\nplt.xlabel(\"Label\", fontsize=17)\nplt.ylabel(\"Frequency\", fontsize=17);","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-06-07T22:22:28.304181Z","iopub.execute_input":"2021-06-07T22:22:28.304492Z","iopub.status.idle":"2021-06-07T22:22:28.572598Z","shell.execute_reply.started":"2021-06-07T22:22:28.304438Z","shell.execute_reply":"2021-06-07T22:22:28.571805Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We will visualize a random image from every label to get a general sense of the distinctive features that seperate the classes. We will take this into account and try to enhance these features in our preprocessing. For these images there some to be increasingly more spots and stains on the retina as diabetic retinopathy worsens.","metadata":{}},{"cell_type":"code","source":"# Example from every label\nfig, ax = plt.subplots(1, 5, figsize=(15, 6))\nfor i in range(5):\n    sample = train_df[train_df['diagnosis'] == i].sample(1)\n    image_name = sample['id_code'].item()\n    X = cv2.imread(f\"{TRAIN_IMG_PATH}{image_name}\")\n    ax[i].set_title(f\"Image: {image_name}\\n Label = {sample['diagnosis'].item()}\", \n                    weight='bold', fontsize=10)\n    ax[i].axis('off')\n    ax[i].imshow(X);","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-06-07T22:22:29.965695Z","iopub.execute_input":"2021-06-07T22:22:29.965994Z","iopub.status.idle":"2021-06-07T22:22:31.708017Z","shell.execute_reply.started":"2021-06-07T22:22:29.965944Z","shell.execute_reply":"2021-06-07T22:22:31.707146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preprocessing <a id=\"5\"></a>","metadata":{}},{"cell_type":"markdown","source":"Here we will use the auto-cropping method with Ben's preprocessing as explained in [this kernel](https://www.kaggle.com/ratthachat/aptos-updatedv14-preprocessing-ben-s-cropping).","metadata":{}},{"cell_type":"code","source":"def crop_image_from_gray(img, tol=7):\n    \"\"\"\n    Applies masks to the orignal image and \n    returns the a preprocessed image with \n    3 channels\n    \n    :param img: A NumPy Array that will be cropped\n    :param tol: The tolerance used for masking\n    \n    :return: A NumPy array containing the cropped image\n    \"\"\"\n    # If for some reason we only have two channels\n    if img.ndim == 2:\n        mask = img > tol\n        return img[np.ix_(mask.any(1),mask.any(0))]\n    # If we have a normal RGB images\n    elif img.ndim == 3:\n        gray_img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)\n        mask = gray_img > tol\n        \n        check_shape = img[:,:,0][np.ix_(mask.any(1),mask.any(0))].shape[0]\n        if (check_shape == 0): # image is too dark so that we crop out everything,\n            return img # return original image\n        else:\n            img1=img[:,:,0][np.ix_(mask.any(1),mask.any(0))]\n            img2=img[:,:,1][np.ix_(mask.any(1),mask.any(0))]\n            img3=img[:,:,2][np.ix_(mask.any(1),mask.any(0))]\n            img = np.stack([img1,img2,img3],axis=-1)\n        return img\n\ndef preprocess_image(image, sigmaX=10):\n    \"\"\"\n    The whole preprocessing pipeline:\n    1. Read in image\n    2. Apply masks\n    3. Resize image to desired size\n    4. Add Gaussian noise to increase Robustness\n    \n    :param img: A NumPy Array that will be cropped\n    :param sigmaX: Value used for add GaussianBlur to the image\n    \n    :return: A NumPy array containing the preprocessed image\n    \"\"\"\n    image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n    image = crop_image_from_gray(image)\n    image = cv2.resize(image, (IMG_WIDTH, IMG_HEIGHT))\n    image = cv2.addWeighted (image,4, cv2.GaussianBlur(image, (0,0) ,sigmaX), -4, 128)\n    return image","metadata":{"execution":{"iopub.status.busy":"2021-06-07T22:22:32.036610Z","iopub.execute_input":"2021-06-07T22:22:32.036896Z","iopub.status.idle":"2021-06-07T22:22:32.049558Z","shell.execute_reply.started":"2021-06-07T22:22:32.036848Z","shell.execute_reply":"2021-06-07T22:22:32.048582Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"After preprocessing we have managed to enhance the distinctive features in the images. This will increase performance when we train our EfficientNet model.","metadata":{}},{"cell_type":"code","source":"# Example of preprocessed images from every label\nfig, ax = plt.subplots(1, 5, figsize=(15, 6))\nfor i in range(5):\n    sample = train_df[train_df['diagnosis'] == i].sample(1)\n    image_name = sample['id_code'].item()\n    X = preprocess_image(cv2.imread(f\"{TRAIN_IMG_PATH}{image_name}\"))\n    ax[i].set_title(f\"Image: {image_name}\\n Label = {sample['diagnosis'].item()}\", \n                    weight='bold', fontsize=10)\n    ax[i].axis('off')\n    ax[i].imshow(X);","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-06-07T22:22:33.505734Z","iopub.execute_input":"2021-06-07T22:22:33.506044Z","iopub.status.idle":"2021-06-07T22:22:35.505243Z","shell.execute_reply.started":"2021-06-07T22:22:33.505990Z","shell.execute_reply":"2021-06-07T22:22:35.504595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Modeling <a id=\"6\"></a>","metadata":{}},{"cell_type":"code","source":"train_df['diagnosis'] = train_df['diagnosis'].astype(str)","metadata":{"execution":{"iopub.status.busy":"2021-06-07T22:22:35.506975Z","iopub.execute_input":"2021-06-07T22:22:35.507419Z","iopub.status.idle":"2021-06-07T22:22:35.516076Z","shell.execute_reply.started":"2021-06-07T22:22:35.507219Z","shell.execute_reply":"2021-06-07T22:22:35.514890Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# We use a small batch size so we can handle large images easily\nBATCH_SIZE = 32\n\n# Add Image augmentation to our generator\ndatagen = ImageDataGenerator(zoom_range=0.15, rotation_range=120,\n                             horizontal_flip=True,\n                             vertical_flip=True,\n                             validation_split=0.15,\n                             shear_range=0.1,\n                             fill_mode='nearest',\n                             preprocessing_function=preprocess_image, \n                             rescale=1 / 255.)\n\n# Use the dataframe to define train and validation generators\ntrain_generator = datagen.flow_from_dataframe(train_df, \n                                                    x_col='id_code', \n                                                    y_col='diagnosis',\n                                                    directory = TRAIN_IMG_PATH,\n                                                    target_size=(IMG_WIDTH, IMG_HEIGHT),\n                                                    batch_size=BATCH_SIZE,\n                                                    color_mode = 'rgb',\n                                                    shuffle=True,\n                                                    class_mode=\"categorical\", #'sparse', \n                                                    subset='training',\n                                                    seed = 42)\n\nval_generator = datagen.flow_from_dataframe(train_df, \n                                                  x_col='id_code', \n                                                  y_col='diagnosis',\n                                                  directory = TRAIN_IMG_PATH,\n                                                  target_size=(IMG_WIDTH, IMG_HEIGHT),\n                                                  batch_size=BATCH_SIZE,\n                                                  color_mode = 'rgb',\n                                                  shuffle=False,\n                                                  class_mode=\"categorical\", #'sparse',\n                                                  subset='validation',\n                                                  seed = 42)","metadata":{"execution":{"iopub.status.busy":"2021-06-07T22:22:38.478942Z","iopub.execute_input":"2021-06-07T22:22:38.479234Z","iopub.status.idle":"2021-06-07T22:22:48.228834Z","shell.execute_reply.started":"2021-06-07T22:22:38.479185Z","shell.execute_reply":"2021-06-07T22:22:48.228021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t_x, t_y = next(train_generator)\n\nfig, m_axs = plt.subplots(2, 4, figsize = (16, 8))\nfor (c_x, c_y, c_ax) in zip(t_x, t_y, m_axs.flatten()):\n    c_ax.imshow(np.clip(c_x*127+127, 0, 255).astype(np.uint8))\n    c_ax.set_title('Severity {}'.format(c_y))\n    c_ax.axis('off')","metadata":{"execution":{"iopub.status.busy":"2021-06-07T22:22:48.230410Z","iopub.execute_input":"2021-06-07T22:22:48.230888Z","iopub.status.idle":"2021-06-07T22:22:54.198963Z","shell.execute_reply.started":"2021-06-07T22:22:48.230834Z","shell.execute_reply":"2021-06-07T22:22:54.193100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load in EfficientNetB5\nfrom efficientnet import EfficientNetB0\nfrom keras.applications import ResNet50\n\n# model 1 \neffnet = EfficientNetB0(weights=None,\n                        include_top=False,\n                        input_shape=(IMG_WIDTH, IMG_HEIGHT, CHANNELS))\neffnet.load_weights('../input/efficientnet-keras-weights-b0b5/efficientnet-b0_imagenet_1000_notop.h5')\n\n# model 2 \n# renet = ResNet50(weights = 'imagenet', include_top = False, input_shape=(IMG_WIDTH, IMG_HEIGHT, CHANNELS))","metadata":{"execution":{"iopub.status.busy":"2021-06-07T22:22:54.200668Z","iopub.execute_input":"2021-06-07T22:22:54.201084Z","iopub.status.idle":"2021-06-07T22:23:06.024470Z","shell.execute_reply.started":"2021-06-07T22:22:54.200904Z","shell.execute_reply":"2021-06-07T22:23:06.023716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import keras \n\ndef build_model():\n    model = Sequential()\n    model.add(effnet)\n    model.add(GlobalAveragePooling2D())\n    model.add(Dense(5, activation=\"softmax\"))\n    model.compile(loss='categorical_crossentropy', \n                  optimizer=keras.optimizers.Adam(lr=0.001), \n                  metrics=['acc'])\n    print(model.summary())\n    return model\n\n# Initialize model\nmodel = build_model()","metadata":{"execution":{"iopub.status.busy":"2021-06-07T22:23:06.028197Z","iopub.execute_input":"2021-06-07T22:23:06.028471Z","iopub.status.idle":"2021-06-07T22:23:09.562355Z","shell.execute_reply.started":"2021-06-07T22:23:06.028424Z","shell.execute_reply":"2021-06-07T22:23:09.561500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Monitor val loss to avoid overfitting and save best mode\ncheckpoint = keras.callbacks.ModelCheckpoint(\n    filepath='model.h5', \n    verbose=1, monitor='val_loss', \n    save_weights_only=True, save_best_only=True\n)   \nrlr = keras.callbacks.ReduceLROnPlateau(monitor='val_loss', \n                                        factor=0.5, patience=4, \n                                        verbose=1, mode='auto', \n                                        epsilon=0.0001)\n\n# Begin training\nmodel.fit_generator(train_generator,\n                    steps_per_epoch=train_generator.samples // BATCH_SIZE,\n                    epochs=10,\n                    validation_data=val_generator,\n                    validation_steps = val_generator.samples // BATCH_SIZE,\n                    callbacks=[checkpoint, rlr])","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-06-07T22:23:09.563678Z","iopub.execute_input":"2021-06-07T22:23:09.563967Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# load best saved weights\nmodel.load_weights('./model.h5')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history_df = pd.DataFrame(model.history.history)\nhistory_df[['loss', 'val_loss']].plot(figsize=(12,5))\n# plt.title(\"Loss\", fontsize=16, weight='bold')\nplt.xlabel(\"Epoch\")\nplt.ylabel(\"Loss\")\nhistory_df[['acc', 'val_acc']].plot(figsize=(12,5))\nplt.title(\"Accuracy\", fontsize=16, weight='bold')\nplt.xlabel(\"Epoch\")\nplt.ylabel(\"% Accuracy\");","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Evaluation <a id=\"7\"></a>\n\n- Confusion Matrix \n- Classification Report \n- Weighted Q. Kappa\n- ROC AUC Curve ","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import classification_report, confusion_matrix\nfrom sklearn.metrics import accuracy_score, classification_report\n\n#Confution Matrix and Classification Report\nval_generator.reset()\ny_pred = model.predict_generator(val_generator,\n                                 steps=val_generator.samples / BATCH_SIZE,\n                                 verbose = True)\ny_pred = np.argmax(y_pred, axis=1)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(val_generator.classes), len(y_pred) # both should be same ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"val_generator.classes[:5], y_pred[:5] # sanity check ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Confusion Matrix ","metadata":{}},{"cell_type":"code","source":"import seaborn as sns\nplt.figure(figsize=(10,5))\n\nconf_mat = confusion_matrix(val_generator.classes, y_pred)\nsns.heatmap(conf_mat, annot=True, fmt=\"d\", cbar = False, \n            cmap = plt.cm.Blues)\n\nplt.ylabel('True label')\nplt.xlabel('Predicted label')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Classification Report","metadata":{}},{"cell_type":"code","source":"# list(valid_generator.class_indices.keys())\ntarget = ['No DR', 'Mild', 'Moderate', 'Severe', 'Proliferative DR']\nprint('Classification Report')\n\nprint('Validation Acc: %2.2f%%' %(100*accuracy_score(val_generator.classes, y_pred)))\nprint(classification_report(val_generator.classes, y_pred, \n                            target_names = target))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Weighted Cohan Kappa","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import cohen_kappa_score\n\ncohen_kappa_score(val_generator.classes, \n                  y_pred, weights='quadratic')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## ROC AUC curve","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics import roc_curve, auc, roc_auc_score\nfrom sklearn.preprocessing import LabelBinarizer\n\nfig, c_ax = plt.subplots(1,1, figsize = (12, 10))\nsns.set(style=\"darkgrid\")\nall_labels = np.unique(train_df['diagnosis'])\n\ndef multiclass_roc_auc_score(y_test, y_pred, average=\"macro\"):\n    lb = LabelBinarizer()\n    lb.fit(y_test)\n    y_test = lb.transform(y_test)\n    y_pred = lb.transform(y_pred)\n\n    for (idx, c_label) in enumerate(all_labels):\n        fpr, tpr, thresholds = roc_curve(y_test[:,idx].astype(int), y_pred[:,idx])\n        c_ax.plot(fpr, tpr, label = '%s (AUC:%0.2f)'  % (c_label, auc(fpr, tpr)))\n    c_ax.plot(fpr, fpr, 'b-', label = 'Random Guessing')\n    return roc_auc_score(y_test, y_pred, average=average)\n\nprint('ROC AUC score:', multiclass_roc_auc_score(val_generator.classes, y_pred))\n\nc_ax.legend()\nc_ax.set_xlabel('False Positive Rate')\nc_ax.set_ylabel('True Positive Rate')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}