{"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":"markdown","source":"# Team 45 Notebook\n\nHello everyone, and welcome to our Notebook !  \nIt is a summary of the various experimentations we tried during this challenge.\n\n\nWe experimented with differents methods :\n# Table of Contents\n1. [Data visualisation](#data-visualisation)\n2. [Data preparation](#data-preparation)\n3. [Preprocessing](#preprocessing)\n4. [Models](#models)\n5. [Loss function](#loss-function)\n6. [Training](#training)\n7. [Validation](#validation)\n8. [Best model overall: baseline tuning](#baseline-tuning)","metadata":{}},{"cell_type":"markdown","source":"## Imports \nFirst, we import all the library needed for our work.  \nWe have decided to use Tensorflow rather than PyTorch, because of our better understanding of this package. ","metadata":{}},{"cell_type":"code","source":"# Globals imports\nimport os\nimport json\nimport csv\nimport random\nimport pathlib\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom os import listdir\n\n# Pandas for the submission\nimport pandas as pd\n\n# Tensorflow 2.0\nimport tensorflow as tf\nprint(\"Tensorflow version :\", tf.__version__)\n\n# for our models\nfrom scipy.ndimage.measurements import label\nfrom sklearn.cluster import KMeans, MeanShift, DBSCAN, AgglomerativeClustering\n\n# Image processing\nimport cv2\nfrom PIL import Image\nimport skimage\nimport skimage.segmentation as seg\nimport skimage.filters as filters\nimport skimage.draw as draw\nimport skimage.color as color\nfrom skimage.color import rgb2gray\nfrom skimage import io\n# morphology for Veins visualization\nfrom skimage import morphology\nfrom skimage import measure\n\n\nfrom scipy import ndimage\n\n\n# keras imports for our models\nimport keras\nfrom keras.utils import to_categorical\nimport keras.backend as K\n\nfrom tensorflow.keras import layers, models","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-05-23T14:56:58.595856Z","iopub.execute_input":"2021-05-23T14:56:58.596198Z","iopub.status.idle":"2021-05-23T14:56:58.606923Z","shell.execute_reply.started":"2021-05-23T14:56:58.596165Z","shell.execute_reply":"2021-05-23T14:56:58.605856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next, we define the different paths of the 3 datasets: train, test and validation.","metadata":{}},{"cell_type":"code","source":"# PATHS\n\ndata_path = \"../input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data\"\ntrain_path = data_path + \"/train\"\ntest_path = data_path + \"/test\"\nval_path = data_path + \"/val\"","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:56:58.608434Z","iopub.execute_input":"2021-05-23T14:56:58.608897Z","iopub.status.idle":"2021-05-23T14:56:58.622015Z","shell.execute_reply.started":"2021-05-23T14:56:58.608859Z","shell.execute_reply":"2021-05-23T14:56:58.621230Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Visualization <a name=\"data-visualisation\"></a>\n\nThis part is here to show you how we try to extract features from the images, by visualizing veins, cup etc...  ","metadata":{}},{"cell_type":"markdown","source":"### Veins visualization","metadata":{}},{"cell_type":"code","source":"\"\"\"\nReturn a mask containing the veins detected from image\nBackground should contain the background mask\n\"\"\"\ndef veins(image, background):\n    \n    #1st threshold\n    threshold = filters.threshold_local(image,block_size=91) \n    image = image <= threshold\n    \n    #Blur the result to get ride of high frequency details\n    sigma = 4\n    image = skimage.filters.gaussian(image, sigma=(sigma, sigma), truncate=3.5)\n    \n    #2nd threshold, and small imperfections removal\n    image = image > 0.62\n    image = morphology.remove_small_objects(image, min_size=100_000)\n    \n    #Remove artifacts outside the eye\n    image = np.logical_and(image, background)\n    \n    return image","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:56:58.623605Z","iopub.execute_input":"2021-05-23T14:56:58.624145Z","iopub.status.idle":"2021-05-23T14:56:58.633344Z","shell.execute_reply.started":"2021-05-23T14:56:58.624107Z","shell.execute_reply":"2021-05-23T14:56:58.632559Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nReturn a mask containing the nerve detected from image\nBackground should contain the background mask\n\"\"\"\ndef nerve(image, background):\n\n    #Blur the image to get ride of high frequency details\n    sigma = 10\n    image = skimage.filters.gaussian(image, sigma=(sigma, sigma), truncate=3.5)\n    \n    #Threshold the nerve\n    threshold = filters.threshold_local(image, block_size=1001, offset=-20/100) \n    image = image > threshold\n    \n    #Identify the resulting areas, and select the bigger one\n    labels = measure.label(image, background=0)\n    props = measure.regionprops(labels)\n    m = -1\n    i = 0\n    for p in props:\n        if p['Area']>m and p['Label']!=0:\n            m = p['Area']\n            i = p['Label']\n    image = labels == i\n    \n    #Remove artifacts outside the eye\n    image = np.logical_and(image, background)\n    \n    return image","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:56:58.637021Z","iopub.execute_input":"2021-05-23T14:56:58.637401Z","iopub.status.idle":"2021-05-23T14:56:58.644831Z","shell.execute_reply.started":"2021-05-23T14:56:58.637374Z","shell.execute_reply":"2021-05-23T14:56:58.643912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nReturn a mask containing everything but the background of image\n\"\"\"\ndef background(image):\n    image = image > 0\n    return image","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:56:58.648899Z","iopub.execute_input":"2021-05-23T14:56:58.649158Z","iopub.status.idle":"2021-05-23T14:56:58.657718Z","shell.execute_reply.started":"2021-05-23T14:56:58.649136Z","shell.execute_reply":"2021-05-23T14:56:58.656903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nRead the image from path, and returns it resized\n\"\"\"\ndef read_original(path):\n    image = cv2.imread(path)\n    return cv2.resize(image, (500, 500), interpolation = cv2.INTER_NEAREST).astype('uint8')","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:56:58.659079Z","iopub.execute_input":"2021-05-23T14:56:58.659661Z","iopub.status.idle":"2021-05-23T14:56:58.667418Z","shell.execute_reply.started":"2021-05-23T14:56:58.659586Z","shell.execute_reply":"2021-05-23T14:56:58.666657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nReturn various vizualizations of the image in path\n--> The original image, a black and white version, and another black and white version with equalized histograms\n--> The three channels (B,G,R)\n\"\"\"\ndef visu_channels(path):\n    \n    #Read the image\n    image = cv2.imread(path)\n    #Resize it\n    image = cv2.resize(image, (400, 400), interpolation = cv2.INTER_NEAREST).astype('uint8')\n    #Extract the three channels\n    b_image, g_image, r_image = cv2.split(image)\n    #Copy the orignal image, to show it later\n    image_src = np.copy(image)\n    \n    #Compute the first line of the result (the original image, a black and white version, and another black and white version with equalized histograms)\n    def visu_channels_aux1(image):\n        bw_image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)\n        eq_image = cv2.equalizeHist(bw_image)\n        bw_image, eq_image = cv2.merge((bw_image, bw_image, bw_image)), cv2.merge((eq_image, eq_image, eq_image))\n        return np.concatenate((image, bw_image, eq_image), axis=1)\n    \n    #Compute the second line of the result (the three channels (B,G,R))\n    def visu_channels_aux2(b_image, g_image, r_image):\n        b_image, g_image, r_image = cv2.merge((b_image, b_image, b_image)), cv2.merge((g_image, g_image, g_image)), cv2.merge((r_image, r_image, r_image))\n        b_image, g_image, r_image = cv2.resize(b_image, (400, 400), interpolation = cv2.INTER_NEAREST).astype('uint8'), cv2.resize(g_image, (400, 400), interpolation = cv2.INTER_NEAREST).astype('uint8'), cv2.resize(r_image, (400, 400), interpolation = cv2.INTER_NEAREST).astype('uint8')\n        return np.concatenate((b_image, g_image, r_image), axis=1)\n    \n    #Return the result\n    return np.concatenate((visu_channels_aux1(image), visu_channels_aux2(b_image, g_image, r_image)), axis=0)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:56:58.728740Z","iopub.execute_input":"2021-05-23T14:56:58.728982Z","iopub.status.idle":"2021-05-23T14:56:58.738530Z","shell.execute_reply.started":"2021-05-23T14:56:58.728959Z","shell.execute_reply":"2021-05-23T14:56:58.737584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nReturn an image where the background of the eye, the veins, and the optic nerve are segmented.\n\"\"\"\ndef get_segmented_image(path):\n    \n    #Read the image\n    image = cv2.imread(path)\n    #Convert it to black and white\n    bw_image = rgb2gray(image)\n    #Extract the three channels\n    b_image, g_image, r_image = cv2.split(image)\n    \n    #Compute the various masks\n    b = background(b_image)\n    v = veins(g_image, b)\n    c = nerve(r_image, b)\n    \n    #Merge the mask to create a visual representation\n    w,h = bw_image.shape\n    image = np.zeros(shape=(w, h, 3))\n    for i in range(w):\n        for j in range(h):\n            if v[i][j]:\n                image[i][j] = np.array([255, 0, 0])\n            elif c[i][j]:\n                image[i][j] = np.array([0, 0, 255])\n            elif b[i][j]:\n                image[i][j] = np.array([200, 200, 200])\n            else:\n                image[i][j] = np.array([0, 0, 0])\n    \n    #Return the original image and the segmented version, side by side\n    image = cv2.resize(image, (500, 500), interpolation = cv2.INTER_NEAREST).astype('uint8')\n    image_src = read_original(path)\n    return np.concatenate((image_src, np.zeros(shape=(500, 50, 3), dtype=np.uint8), image), axis=1)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:56:58.740040Z","iopub.execute_input":"2021-05-23T14:56:58.740463Z","iopub.status.idle":"2021-05-23T14:56:58.752670Z","shell.execute_reply.started":"2021-05-23T14:56:58.740430Z","shell.execute_reply":"2021-05-23T14:56:58.751842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path = os.path.join(data_path, 'train', 'images', 'g0001.jpg')","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:56:58.754470Z","iopub.execute_input":"2021-05-23T14:56:58.755002Z","iopub.status.idle":"2021-05-23T14:56:58.766606Z","shell.execute_reply.started":"2021-05-23T14:56:58.754966Z","shell.execute_reply":"2021-05-23T14:56:58.765853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize = (18,12))\nplt.imshow(cv2.cvtColor(visu_channels(path), cv2.COLOR_BGR2RGB))","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:56:58.768269Z","iopub.execute_input":"2021-05-23T14:56:58.768642Z","iopub.status.idle":"2021-05-23T14:56:59.352825Z","shell.execute_reply.started":"2021-05-23T14:56:58.768607Z","shell.execute_reply":"2021-05-23T14:56:59.351828Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This image allows us to see the same image, differently:\n* 1st line: The original image, a black and white version, and another black and white version with equalized histograms\n* 2nd line: The three channels (B,G,R)","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize = (18,9))\nplt.imshow(cv2.cvtColor(get_segmented_image(path), cv2.COLOR_BGR2RGB))","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:56:59.354193Z","iopub.execute_input":"2021-05-23T14:56:59.354536Z","iopub.status.idle":"2021-05-23T14:57:30.333051Z","shell.execute_reply.started":"2021-05-23T14:56:59.354501Z","shell.execute_reply":"2021-05-23T14:57:30.332112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"On the left is the original image, and on the right is the segmented version, with:\n* In grey the background of the image\n* In red the optic nerve\n* In blue the veins","metadata":{}},{"cell_type":"markdown","source":"### Global segmentation\nWe try to segment the image by using various algorithms, e.g. Kmeans. <p>","metadata":{}},{"cell_type":"code","source":"\"\"\"\nAuxiliary function that generates an array with the colors used for the visualisation.\n\"\"\"\ndef gen_colors(n=10):\n    L = [0, 128, 255] if n>8 else [0, 255]\n    return np.array([[b,g,r] for r in L for g in L for b in L])","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:57:30.334435Z","iopub.execute_input":"2021-05-23T14:57:30.334820Z","iopub.status.idle":"2021-05-23T14:57:30.340998Z","shell.execute_reply.started":"2021-05-23T14:57:30.334782Z","shell.execute_reply":"2021-05-23T14:57:30.339776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Various hyperparameters for the segmentation\n\nSRC_DIRECTORY   = f'{data_path}/train/images/'\n\nNB_CLUSTERS     = 6\nBANDWIDTH       = 20\nEPSILON         = 10\nMIN_SAMPLES     = 100\n\nCOLORS          = gen_colors()\nSIZE            = 150\nSIZE_FACTOR     = 1/20\nBW_FACTOR       = 0","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:57:30.342422Z","iopub.execute_input":"2021-05-23T14:57:30.342870Z","iopub.status.idle":"2021-05-23T14:57:30.354673Z","shell.execute_reply.started":"2021-05-23T14:57:30.342834Z","shell.execute_reply":"2021-05-23T14:57:30.353856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Although we do not really need 6 classes, this number was chosen because:\n* A smaller number of classes will not distinguish the OC from the OD\n* A bigger number of classes just adds noise\n\n<p> Other parameters were chosen to get similar results.","metadata":{}},{"cell_type":"code","source":"\"\"\"\nConvert a 2D array (coordinates) with 5D points (RGB + coordinates) into a 1D array with 3D points (RGB)\n\"\"\"\ndef pixels_list(X):\n    L = []\n    for i in range(len(X)):\n        for j in range(len(X[0])):\n            b,g,r,gr = X[i][j]\n            L.append([i*SIZE_FACTOR,j*SIZE_FACTOR,b,g,r,gr*BW_FACTOR])\n    return np.array(L)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:57:30.357302Z","iopub.execute_input":"2021-05-23T14:57:30.357661Z","iopub.status.idle":"2021-05-23T14:57:30.365566Z","shell.execute_reply.started":"2021-05-23T14:57:30.357625Z","shell.execute_reply":"2021-05-23T14:57:30.364787Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nReturns the clustered result from 'algo' applied on the image in 'img_path'\n\"\"\"\ndef cluster(img_path, algo):\n    \n    img = cv2.imread(img_path)\n    img = cv2.resize(img, (SIZE, SIZE), interpolation = cv2.INTER_AREA)\n    \n    img_src = img.copy()\n    img_dst = img.copy()\n    \n    gr_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n    gr_img = cv2.equalizeHist(gr_img)\n    b_img, g_img, r_img = cv2.split(img)\n    img = cv2.merge((b_img, g_img, r_img, gr_img))\n    \n    if algo==\"KMeans\":\n        clustering = KMeans(n_clusters=NB_CLUSTERS, random_state=0).fit(pixels_list(img))\n    elif algo==\"MeanShift\":\n        clustering = MeanShift(bandwidth=BANDWIDTH).fit(pixels_list(img))\n    elif algo==\"DBSCAN\":\n        clustering = DBSCAN(eps=EPSILON, min_samples=MIN_SAMPLES, n_jobs=-1).fit(pixels_list(img))\n    else:\n        clustering = AgglomerativeClustering(n_clusters = NB_CLUSTERS).fit(pixels_list(img))\n    \n    labels = list(clustering.labels_)[::-1]\n    for i in range(len(img_dst)):\n        for j in range(len(img_dst[0])):\n            img_dst[i][j] = COLORS[labels.pop()]\n    return img_src, img_dst","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:57:30.368921Z","iopub.execute_input":"2021-05-23T14:57:30.369187Z","iopub.status.idle":"2021-05-23T14:57:30.379443Z","shell.execute_reply.started":"2021-05-23T14:57:30.369163Z","shell.execute_reply":"2021-05-23T14:57:30.378568Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Image segmented with the KMeans algorithm:","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize = (14,7))\nx, y   = cluster(SRC_DIRECTORY + 'g0001.jpg', \"KMeans\")\nresult = cv2.resize(np.concatenate((x, y), axis=1), (600, 300), interpolation = cv2.INTER_AREA)\nplt.imshow(result)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:57:30.380585Z","iopub.execute_input":"2021-05-23T14:57:30.381103Z","iopub.status.idle":"2021-05-23T14:57:31.515134Z","shell.execute_reply.started":"2021-05-23T14:57:30.381068Z","shell.execute_reply":"2021-05-23T14:57:31.514349Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Image segmented with the MeanShift algorithm:","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize = (14,7))\nx, y   = cluster(SRC_DIRECTORY + 'g0001.jpg', \"MeanShift\")\nresult = cv2.resize(np.concatenate((x, y), axis=1), (600, 300), interpolation = cv2.INTER_AREA)\nplt.imshow(result)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T14:57:31.516447Z","iopub.execute_input":"2021-05-23T14:57:31.516798Z","iopub.status.idle":"2021-05-23T15:03:40.584529Z","shell.execute_reply.started":"2021-05-23T14:57:31.516761Z","shell.execute_reply":"2021-05-23T15:03:40.583560Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Image segmented with the DBSCAN algorithm:","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize = (14,7))\nx, y   = cluster(SRC_DIRECTORY + 'g0001.jpg', \"DBSCAN\")\nresult = cv2.resize(np.concatenate((x, y), axis=1), (600, 300), interpolation = cv2.INTER_AREA)\nplt.imshow(result)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T15:03:40.586024Z","iopub.execute_input":"2021-05-23T15:03:40.586384Z","iopub.status.idle":"2021-05-23T15:03:43.896391Z","shell.execute_reply.started":"2021-05-23T15:03:40.586347Z","shell.execute_reply":"2021-05-23T15:03:43.893626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Image segmented with the AgglomerativeClustering algorithm:","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize = (14,7))\nx, y   = cluster(SRC_DIRECTORY + 'g0001.jpg', \"AgglomerativeClustering\")\nresult = cv2.resize(np.concatenate((x, y), axis=1), (600, 300), interpolation = cv2.INTER_AREA)\nplt.imshow(result)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T15:03:43.897650Z","iopub.execute_input":"2021-05-23T15:03:43.898001Z","iopub.status.idle":"2021-05-23T15:04:03.655105Z","shell.execute_reply.started":"2021-05-23T15:03:43.897965Z","shell.execute_reply":"2021-05-23T15:04:03.654304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Conclusion\nThose methods worked quite well. However, the segmentation method of the baseline still had a better accuracy. Indeed, the images segmented using our various methods are always somehow noisy. The result is therefore not exploitable.","metadata":{}},{"cell_type":"markdown","source":"## Data Preparation <a name=\"data-preparation\"></a>\n\n### Region of interest\nWe have decided to use algorithm in order to find the region of interest of each images. \nWe think that, by cropping around the cup, the models must better predict the cup size.  \nThen, we will able to use the ROI everywhere, instead of full images, in each methods.","metadata":{}},{"cell_type":"code","source":"ROI_OFFSET = 300 # the size of the cropped ROI \ndef get_roi(image):\n    b,g,r = cv2.split(image)\n    \n    # First we apply a Gaussian Blur on the green channel of the image    \n    g = cv2.GaussianBlur(g,(15,15),0)\n    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(15,15))\n    g = ndimage.grey_opening(g,structure=kernel)\n    \n    # Finally we apply the minMaxLoc function to crop the optic nerve which is very bright on the image\n    (minVal, maxVal, minLoc, maxLoc) = cv2.minMaxLoc(g)\n    x0 = max(int(maxLoc[0])- ROI_OFFSET,0)\n    y0 = max(int(maxLoc[1])- ROI_OFFSET,0)\n    x1 = max(int(maxLoc[0])+ ROI_OFFSET,0)\n    y1 = max(int(maxLoc[1])+ ROI_OFFSET,0)\n    return (y0, y1,x0, x1)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T15:04:03.656308Z","iopub.execute_input":"2021-05-23T15:04:03.656637Z","iopub.status.idle":"2021-05-23T15:04:03.663966Z","shell.execute_reply.started":"2021-05-23T15:04:03.656601Z","shell.execute_reply":"2021-05-23T15:04:03.663175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"example_image = cv2.imread(train_path + '/images/g0001.jpg')\nexample_colored_image = cv2.cvtColor(example_image, cv2.COLOR_BGR2RGB)\n\n\nplt.subplot(1,3,1)\nplt.imshow(example_colored_image)\nplt.title('Original image')\n\nplt.subplot(1,3,2)\nb,g,r = cv2.split(example_colored_image)\ng = cv2.GaussianBlur(g,(15,15),0)\nplt.imshow(g)\nplt.title('Channel division')\n\nplt.subplot(1,3,3)\n(y0, y1,x0, x1) = get_roi(example_colored_image)\nROI_image = example_colored_image[y0:y1, x0:x1]\nplt.imshow(ROI_image)\nplt.title(\"after MaxLoc\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2021-05-23T15:04:03.665391Z","iopub.execute_input":"2021-05-23T15:04:03.666049Z","iopub.status.idle":"2021-05-23T15:04:08.532841Z","shell.execute_reply.started":"2021-05-23T15:04:03.666009Z","shell.execute_reply":"2021-05-23T15:04:08.531821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Darkside of roi cropping** : our algorithm takes around 3 ans 4 seconds to process each images : the load time of each dataset will be impacted...\n","metadata":{}},{"cell_type":"markdown","source":"## Dataset Pre-processing <a name=\"pre-processing\"></a>\n\nWe will use the ROI method for each dataset, in order to load cropped dataset in memory  \n**Remarks**: We will not use the informations of the fovea position, available in the index.json, which become false by cropping the image","metadata":{}},{"cell_type":"code","source":"def load_images(path,split):\n    cdr_list = []\n\n    k = 0\n    length = len(listdir(path))\n    for filename in listdir(path):\n        \n        print(f'Loading {split} image {k}/{length}...', end='\\r')\n        k += 1\n        image = cv2.imread(path + \"/\" + filename)\n        colored_image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)        \n        (x0, x1, y0, y1) = get_roi(colored_image)\n        ROI_image = colored_image[x0:x1, y0:y1]        \n        resized_image = cv2.resize(ROI_image, (256,256), interpolation = cv2.INTER_AREA)                             \n        cdr_list.append(resized_image)\n                        \n    return(cdr_list)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T15:04:08.534189Z","iopub.execute_input":"2021-05-23T15:04:08.534582Z","iopub.status.idle":"2021-05-23T15:04:08.542045Z","shell.execute_reply.started":"2021-05-23T15:04:08.534546Z","shell.execute_reply":"2021-05-23T15:04:08.541133Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_images = load_images(train_path + \"/images\",\"train\")\nval_images = load_images(val_path + \"/images\",\"val\")\ntest_images = load_images(test_path + \"/images\",\"test\")","metadata":{"execution":{"iopub.status.busy":"2021-05-23T15:04:08.543706Z","iopub.execute_input":"2021-05-23T15:04:08.544164Z","iopub.status.idle":"2021-05-23T16:02:58.725770Z","shell.execute_reply.started":"2021-05-23T15:04:08.544128Z","shell.execute_reply":"2021-05-23T16:02:58.724873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This process takes around 1 hour to load all datasets...\n\nSecondly, we extract the labels from the index.json file :","metadata":{}},{"cell_type":"code","source":"# Labels :\n\ndef load_json(split):\n    assert split in ['train', 'val', 'test']\n    with open(os.path.join(data_path, split, 'index.json'), 'r') as file: \n        data = json.load(file)\n    labels = list(map(lambda x: x['Label'], data.values()))\n    print(f\"\\033[95m Labels for split {split}\")\n    print(f\"Number of positive :\", sum(labels), '/', len(labels), '\\033[0m' )\n    return labels\n\ntrain_target = load_json('train')\nval_target = load_json('val')\n","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:02:58.727176Z","iopub.execute_input":"2021-05-23T16:02:58.727546Z","iopub.status.idle":"2021-05-23T16:02:58.740327Z","shell.execute_reply.started":"2021-05-23T16:02:58.727506Z","shell.execute_reply":"2021-05-23T16:02:58.739375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Classes imbalance\n\nOne of the main parameters of this challenge is unbalanced data. This makes it difficult for some models to learn, like CNNs.","metadata":{}},{"cell_type":"code","source":"# ones = np.count_nonzero(np.array(train_target) == 1)\n# zeros = np.count_nonzero(np.array(train_target) == 0)\nones = sum(train_target)\nzeros = len(train_target) - ones\n\nclasses = ['Non glaucomatous','Glaucomatous']\n\nplt.bar(classes,[zeros,ones],color='blue',edgecolor='black')\nplt.xticks(classes)\n\nplt.bar([0,1],[zeros,ones])\nplt.xlabel('Class', fontsize=16)\nplt.ylabel('Count', fontsize=16)\nplt.title(\"Diagram of imbalance of the two classes\")\nplt.show()\n\nprint(\"Class imbalance : \",sum(train_target)/len(train_target)*100, \"% of glaucomatous\")\n","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:02:58.741983Z","iopub.execute_input":"2021-05-23T16:02:58.742571Z","iopub.status.idle":"2021-05-23T16:02:58.855057Z","shell.execute_reply.started":"2021-05-23T16:02:58.742533Z","shell.execute_reply":"2021-05-23T16:02:58.854121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Images visualization\n\nLets see our data : how do our images look after cropping and labelizing ? ","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(10,10))\nfor i in range(25, 50):\n    plt.subplot(5,5,i-25+1)\n    plt.xticks([])\n    plt.yticks([])\n    plt.grid(False)\n    plt.imshow(train_images[i], cmap=plt.cm.binary)\n    plt.xlabel('Glaucomatous' if (train_target[i] == 1) else 'Non Glaucomatous')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:02:58.856370Z","iopub.execute_input":"2021-05-23T16:02:58.856711Z","iopub.status.idle":"2021-05-23T16:03:00.510300Z","shell.execute_reply.started":"2021-05-23T16:02:58.856676Z","shell.execute_reply":"2021-05-23T16:03:00.509146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Augmentation <a name=\"data-augmentation\"></a>","metadata":{}},{"cell_type":"markdown","source":"We curently have 400 images in the training set. We have decided to increase this number of data, by adding transformed images in the dataset.","metadata":{}},{"cell_type":"code","source":"def flip(originalImage):\n    flipBoth = cv2.flip(originalImage, -1)\n    return flipBoth","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:00.517332Z","iopub.execute_input":"2021-05-23T16:03:00.517681Z","iopub.status.idle":"2021-05-23T16:03:00.527398Z","shell.execute_reply.started":"2021-05-23T16:03:00.517646Z","shell.execute_reply":"2021-05-23T16:03:00.526203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"example_image = cv2.imread(train_path + '/images/g0001.jpg')\nexample_colored_image = cv2.cvtColor(example_image, cv2.COLOR_BGR2RGB)\n\n\nplt.subplot(1,3,1)\nplt.imshow(example_colored_image)\nplt.title('Original image')\n\nplt.subplot(1,3,2)\nflipped = flip(example_colored_image)\nplt.imshow(flipped)\nplt.title('Image flipped')\n\nplt.subplot(1,3,3)\n(y0, y1,x0, x1) = get_roi(flipped)\nROI_image = flipped[y0:y1, x0:x1]\nplt.imshow(ROI_image)\nplt.title(\"flipped & cropped\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:00.530354Z","iopub.execute_input":"2021-05-23T16:03:00.530833Z","iopub.status.idle":"2021-05-23T16:03:05.410500Z","shell.execute_reply.started":"2021-05-23T16:03:00.530798Z","shell.execute_reply":"2021-05-23T16:03:05.409572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So let's do it for each of the images, and increase the number of data. We also use a treshold, in order to not duplicate all images.","metadata":{}},{"cell_type":"code","source":"def augment_set(images_set,targets, treshold=0.5):\n    for i in range(len(images_set)):\n        u = random.random()\n        if(u <= treshold):\n            images_set.append(flip(images_set[i]))\n            targets.append(targets[i])\n    return images_set,targets\n\ntrain_images, train_target = augment_set(train_images,train_target)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:05.411862Z","iopub.execute_input":"2021-05-23T16:03:05.412202Z","iopub.status.idle":"2021-05-23T16:03:05.428695Z","shell.execute_reply.started":"2021-05-23T16:03:05.412166Z","shell.execute_reply":"2021-05-23T16:03:05.427946Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#convert lists into arrays\ntrain_images = np.array(train_images)\nval_images = np.array(val_images)\ntest_images = np.array(test_images)\n\ntrain_target = np.array(train_target)\nval_target = np.array(val_target)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:05.430024Z","iopub.execute_input":"2021-05-23T16:03:05.430432Z","iopub.status.idle":"2021-05-23T16:03:05.526303Z","shell.execute_reply.started":"2021-05-23T16:03:05.430391Z","shell.execute_reply":"2021-05-23T16:03:05.525429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#targets one-hot encoding\ny_train = to_categorical(np.array(train_target), num_classes=2)\ny_val = to_categorical(np.array(val_target), num_classes=2)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:05.527653Z","iopub.execute_input":"2021-05-23T16:03:05.528011Z","iopub.status.idle":"2021-05-23T16:03:05.532910Z","shell.execute_reply.started":"2021-05-23T16:03:05.527975Z","shell.execute_reply":"2021-05-23T16:03:05.532116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Models<a name=\"models\"></a>\n\n### Custom CNN\n\nIn a first time, we have tried to use our own CNN, which have the advantage to be simple.","metadata":{}},{"cell_type":"code","source":"model = models.Sequential()\n\n# CNN\n\nmodel.add(layers.Conv2D(256, (3, 3), activation='relu', input_shape=(256, 256, 3)))\nmodel.add(layers.MaxPooling2D((2, 2)))\nmodel.add(layers.Conv2D(128, (3, 3), activation='relu'))\nmodel.add(layers.MaxPooling2D((2, 2)))\nmodel.add(layers.Conv2D(64, (3, 3), activation='relu'))\nmodel.add(layers.MaxPooling2D((2, 2)))\nmodel.add(layers.Conv2D(32, (3, 3), activation='relu'))\nmodel.add(layers.MaxPooling2D((2, 2)))\nmodel.add(layers.Conv2D(16, (3, 3), activation='relu'))\nmodel.add(layers.MaxPooling2D((2, 2)))\nmodel.add(layers.Conv2D(8, (3, 3), activation='relu'))\nmodel.add(layers.Dense(2))","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:05.534094Z","iopub.execute_input":"2021-05-23T16:03:05.534586Z","iopub.status.idle":"2021-05-23T16:03:05.621565Z","shell.execute_reply.started":"2021-05-23T16:03:05.534550Z","shell.execute_reply":"2021-05-23T16:03:05.620749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.summary()","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:05.622723Z","iopub.execute_input":"2021-05-23T16:03:05.623066Z","iopub.status.idle":"2021-05-23T16:03:05.634959Z","shell.execute_reply.started":"2021-05-23T16:03:05.623031Z","shell.execute_reply":"2021-05-23T16:03:05.633800Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### ResNet50\n\nWe have decided to use a more powerful and more complex convolutional neural network: ResNet50 and InceptionV3.","metadata":{}},{"cell_type":"code","source":"from tensorflow.keras.applications import ResNet50\n\ncategories = 2\nmodel = models.Sequential()\n\n# The fully connected top layer of ResNet50 is not to added in this model\nmodel.add(ResNet50(include_top = False, pooling = 'avg', weights = 'imagenet'))\n\n# All inputs and outputs are connected to neurons (Dense Layers)\nmodel.add(layers.Dense(categories, activation = \"softmax\"))\nmodel.layers[0].trainable = False","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:05.636313Z","iopub.execute_input":"2021-05-23T16:03:05.636654Z","iopub.status.idle":"2021-05-23T16:03:07.472900Z","shell.execute_reply.started":"2021-05-23T16:03:05.636622Z","shell.execute_reply":"2021-05-23T16:03:07.472076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### InceptionV3\nWe also try an InceptionV3 model","metadata":{}},{"cell_type":"markdown","source":"##### Custom loss function and metrics\n\nTo fight the targets imbalance, we used a custom f1_loss function on one-hot encoded targets. We use a lot of different metrics to keep track of false positives and false negatives ratios.","metadata":{}},{"cell_type":"code","source":"# custom loss functions\n\ndef f1_loss(y_true, y_pred):\n    \n    tp = K.sum(K.cast(y_true*y_pred, 'float'), axis=0)\n    tn = K.sum(K.cast((1-y_true)*(1-y_pred), 'float'), axis=0)\n    fp = K.sum(K.cast((1-y_true)*y_pred, 'float'), axis=0)\n    fn = K.sum(K.cast(y_true*(1-y_pred), 'float'), axis=0)\n\n    p = tp / (tp + fp + K.epsilon())\n    r = tp / (tp + fn + K.epsilon())\n\n    f1 = 2*p*r / (p+r+K.epsilon())\n    f1 = tf.where(tf.math.is_nan(f1), tf.zeros_like(f1), f1)\n    return 1 - K.mean(f1)\n\nMETRICS = [\n      keras.metrics.TruePositives(name='tp'),\n      keras.metrics.FalsePositives(name='fp'),\n      keras.metrics.TrueNegatives(name='tn'),\n      keras.metrics.FalseNegatives(name='fn'), \n      keras.metrics.BinaryAccuracy(name='accuracy'),\n      keras.metrics.Precision(name='precision'),\n      keras.metrics.Recall(name='recall'),\n      keras.metrics.AUC(name='auc'),\n      f1_loss\n]\n\nones = sum(y_train)\nzeros = len(y_train) - ones\n\nCLASS_WEIGHT = {0: zeros,1: ones}","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:07.474257Z","iopub.execute_input":"2021-05-23T16:03:07.474585Z","iopub.status.idle":"2021-05-23T16:03:07.514286Z","shell.execute_reply.started":"2021-05-23T16:03:07.474551Z","shell.execute_reply":"2021-05-23T16:03:07.513451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##### Model definition","metadata":{}},{"cell_type":"code","source":"from tensorflow.keras.applications.inception_v3 import InceptionV3\nfrom tensorflow.keras.optimizers import RMSprop\nfrom keras.layers import Flatten, Dense","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:07.515503Z","iopub.execute_input":"2021-05-23T16:03:07.515837Z","iopub.status.idle":"2021-05-23T16:03:07.520198Z","shell.execute_reply.started":"2021-05-23T16:03:07.515802Z","shell.execute_reply":"2021-05-23T16:03:07.519246Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"HEIGHT = 256\nWIDTH = 256\n\nBATCH_SIZE = 8\nclass_list = [\"0\", \"1\"]\nFC_LAYERS = [1024, 512, 256]\ndropout = 0.5\nNUM_EPOCHS = 100\nBATCH_SIZE = 8","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:07.521622Z","iopub.execute_input":"2021-05-23T16:03:07.521995Z","iopub.status.idle":"2021-05-23T16:03:07.532472Z","shell.execute_reply.started":"2021-05-23T16:03:07.521957Z","shell.execute_reply":"2021-05-23T16:03:07.531563Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Build the InceptionV3 model","metadata":{}},{"cell_type":"code","source":"def build_model(base_model, dropout, fc_layers, num_classes):\n    for layer in base_model.layers:\n        layer.trainable = False\n\n    x = base_model.output\n    x = Flatten()(x)\n    for fc in fc_layers:\n        print(fc)\n        x = layers.Dense(fc, activation='relu')(x)\n        x = layers.Dropout(dropout)(x)\n    preditions = layers.Dense(num_classes, activation='softmax')(x)\n    finetune_model = models.Model(inputs = base_model.input, outputs = preditions)\n    return finetune_model","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:07.533695Z","iopub.execute_input":"2021-05-23T16:03:07.534116Z","iopub.status.idle":"2021-05-23T16:03:07.542295Z","shell.execute_reply.started":"2021-05-23T16:03:07.534082Z","shell.execute_reply":"2021-05-23T16:03:07.541503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"base_model_2 = InceptionV3(weights = 'imagenet',\n                       include_top = False,\n                       input_shape = (256, 256, 3))","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:07.543477Z","iopub.execute_input":"2021-05-23T16:03:07.543985Z","iopub.status.idle":"2021-05-23T16:03:09.362563Z","shell.execute_reply.started":"2021-05-23T16:03:07.543927Z","shell.execute_reply":"2021-05-23T16:03:09.361684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inception_model = build_model(base_model_2, dropout = dropout, fc_layers = FC_LAYERS, num_classes = 2)\ninception_model.compile(optimizer = \"adam\", loss=f1_loss, metrics=METRICS)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:09.364073Z","iopub.execute_input":"2021-05-23T16:03:09.364425Z","iopub.status.idle":"2021-05-23T16:03:09.439181Z","shell.execute_reply.started":"2021-05-23T16:03:09.364385Z","shell.execute_reply":"2021-05-23T16:03:09.438231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Multi-step CNN","metadata":{}},{"cell_type":"code","source":"def CNN_1():\n    \n    model = models.Sequential([\n        \n        layers.Input(shape=(IMG_SIZE, IMG_SIZE, 3)),\n    \n        layers.Conv2D(filters=8,\n               kernel_size=5,\n               padding='same',\n               activation='relu',\n               kernel_regularizer=regularizers.l2(REGULARIZATION_FACTOR_CONV),\n               bias_regularizer=regularizers.l2(REGULARIZATION_FACTOR_CONV)),\n        layers.BatchNormalization(),\n        layers.Dropout(DROPOUT_RATE),\n        \n        layers.Conv2D(filters=16,\n               kernel_size=5,\n               padding='same',\n               activation='relu',\n               kernel_regularizer=regularizers.l2(REGULARIZATION_FACTOR_CONV),\n               bias_regularizer=regularizers.l2(REGULARIZATION_FACTOR_CONV)),\n        layers.BatchNormalization(),\n        layers.Dropout(DROPOUT_RATE),\n        \n        layers.Conv2D(filters=2,\n               kernel_size=5,\n               padding='same',\n               activation='relu',\n               kernel_regularizer=regularizers.l2(REGULARIZATION_FACTOR_CONV),\n               bias_regularizer=regularizers.l2(REGULARIZATION_FACTOR_CONV)),\n        \n    ])\n    \n    model.compile(optimizer='adam',\n                  loss=keras.losses.MeanSquaredError())\n    \n    return model","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:09.440555Z","iopub.execute_input":"2021-05-23T16:03:09.440902Z","iopub.status.idle":"2021-05-23T16:03:09.448873Z","shell.execute_reply.started":"2021-05-23T16:03:09.440863Z","shell.execute_reply":"2021-05-23T16:03:09.447920Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def CNN2():\n        \n    model = models.Sequential([\n        \n        layers.Input(shape=(IMG_SIZE, IMG_SIZE, 2)),\n    \n        layers.Conv2D(filters=8,\n               kernel_size=5,\n               padding='same',\n               activation='relu',\n               kernel_regularizer=regularizers.l2(REGULARIZATION_FACTOR_CONV),\n               bias_regularizer=regularizers.l2(REGULARIZATION_FACTOR_CONV)),\n        layers.BatchNormalization(),\n        layers.Dropout(DROPOUT_RATE),\n        \n        layers.Conv2D(filters=16,\n               kernel_size=5,\n               padding='same',\n               activation='relu',\n               kernel_regularizer=regularizers.l2(REGULARIZATION_FACTOR_CONV),\n               bias_regularizer=regularizers.l2(REGULARIZATION_FACTOR_CONV)),\n        layers.BatchNormalization(),\n        layers.MaxPooling2D(),\n        layers.Dropout(DROPOUT_RATE),\n        \n        layers.Flatten(),\n        \n        layers.Dense(64,\n              activation='relu',\n              kernel_regularizer=regularizers.l2(REGULARIZATION_FACTOR_DENSE),\n              bias_regularizer=regularizers.l2(REGULARIZATION_FACTOR_DENSE)),\n        \n        layers.Dense(1,\n              activation='softmax',\n              kernel_regularizer=regularizers.l2(REGULARIZATION_FACTOR_DENSE),\n              bias_regularizer=regularizers.l2(REGULARIZATION_FACTOR_DENSE))\n        \n    ])\n    \n    model.compile(optimizer='adam',\n                  loss=keras.losses.BinaryCrossentropy(),\n                  metrics=[keras.metrics.BinaryAccuracy()])\n    \n    return model","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:09.450490Z","iopub.execute_input":"2021-05-23T16:03:09.450987Z","iopub.status.idle":"2021-05-23T16:03:09.463388Z","shell.execute_reply.started":"2021-05-23T16:03:09.450943Z","shell.execute_reply":"2021-05-23T16:03:09.462450Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def copy_layer(layer):\n    config = layer.get_config()\n    weights = layer.get_weights()\n    cloned_layer = type(layer).from_config(config)\n    cloned_layer.build(layer.input_shape)\n    cloned_layer.set_weights(weights)\n    return cloned_layer","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:09.464802Z","iopub.execute_input":"2021-05-23T16:03:09.465286Z","iopub.status.idle":"2021-05-23T16:03:09.476365Z","shell.execute_reply.started":"2021-05-23T16:03:09.465249Z","shell.execute_reply":"2021-05-23T16:03:09.475568Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def concat_models(model_1, model_2):\n    \n    model = models.Sequential()\n    model.add(layers.Input(shape=(IMG_SIZE, IMG_SIZE, 3)))\n    \n    for layer in model_1.layers:\n        model.add(copy_layer(layer))\n    for layer in model_2.layers:\n        model.add(copy_layer(layer))\n        \n    model.compile(optimizer='adam',\n                  loss=keras.losses.BinaryCrossentropy(),\n                  metrics=[keras.metrics.BinaryAccuracy()])\n    \n    return model","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:09.477451Z","iopub.execute_input":"2021-05-23T16:03:09.477914Z","iopub.status.idle":"2021-05-23T16:03:09.485411Z","shell.execute_reply.started":"2021-05-23T16:03:09.477883Z","shell.execute_reply":"2021-05-23T16:03:09.484620Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nTrain the first part of the network, then the second one, and then fine-tune the result by training them together\n\"\"\"\n\n# x, seg, y = create_dataset()\n# w         = get_balanced_weights(y)\n# \n# model_1 = CNN_1()\n# history_1 = model_1.fit(\n#     x                 = x,\n#     y                 = seg,\n#     validation_split  = 0.1,\n#     batch_size        = BATCH_SIZE,\n#     epochs            = EPOCHS[0]\n# )\n# \n# model_2 = CNN_2()\n# history_2 = model_2.fit(\n#     x                 = seg,\n#     y                 = y,\n#     class_weight      = w,\n#     validation_split  = 0.1,\n#     batch_size        = BATCH_SIZE,\n#     epochs            = EPOCHS[1]\n# )\n# \n# model_3 = concat_models(model_1, model_2)\n# history_3 = model_3.fit(\n#     x                 = x,\n#     y                 = y,\n#     class_weight      = w,\n#     validation_split  = 0.1,\n#     batch_size        = BATCH_SIZE,\n#     epochs            = EPOCHS[2]\n# )","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:09.486725Z","iopub.execute_input":"2021-05-23T16:03:09.487274Z","iopub.status.idle":"2021-05-23T16:03:09.499931Z","shell.execute_reply.started":"2021-05-23T16:03:09.487237Z","shell.execute_reply":"2021-05-23T16:03:09.499113Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training <a name=\"training\"></a>\n\nWe have decided to fit the ResNet50 model, and try to analyze this result","metadata":{}},{"cell_type":"code","source":"model.compile(optimizer='adam',\n              loss='binary_crossentropy',\n              metrics=METRICS)\n\n#history = model.fit(train_images, y_train, epochs=50, batch_size=32, validation_data=(val_images, y_val), class_weight=CLASS_WEIGHT, verbose=1)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:09.501305Z","iopub.execute_input":"2021-05-23T16:03:09.501689Z","iopub.status.idle":"2021-05-23T16:03:09.518118Z","shell.execute_reply.started":"2021-05-23T16:03:09.501653Z","shell.execute_reply":"2021-05-23T16:03:09.517359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Best model overall : Baseline Tuning <a name=\"baseline-tuning\"></a>\n\nThis part is only here to show you how we have scored **0.98572** on the leaderboard. We have tuned the baseline by adding preprocessing, as explained in the notebook. We have added a ROI crop, and a data augmentation of 40 per cent.\n### ROI + Data Augmentation on RefugeDataset class","metadata":{}},{"cell_type":"code","source":"import os\nimport json\nimport csv\nimport random\nimport pickle\nimport cv2\nfrom scipy import ndimage\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport numba as nb\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torch.optim as optim\nimport torchvision.transforms as transforms\n\nfrom PIL import Image\nfrom torch.utils.data import Dataset, DataLoader\nfrom scipy.ndimage.measurements import label\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.metrics import roc_auc_score, roc_curve","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:09.519097Z","iopub.execute_input":"2021-05-23T16:03:09.519447Z","iopub.status.idle":"2021-05-23T16:03:11.698625Z","shell.execute_reply.started":"2021-05-23T16:03:09.519411Z","shell.execute_reply":"2021-05-23T16:03:11.697797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"root_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:11.699901Z","iopub.execute_input":"2021-05-23T16:03:11.700243Z","iopub.status.idle":"2021-05-23T16:03:11.704588Z","shell.execute_reply.started":"2021-05-23T16:03:11.700193Z","shell.execute_reply":"2021-05-23T16:03:11.703415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Augmentation example","metadata":{}},{"cell_type":"code","source":"with open(os.path.join(root_dir, 'train', 'index.json')) as f:\n    index = json.load(f)\nimg_name = os.path.join(root_dir, 'train', 'images', index['0']['ImgName'])\nimg = np.array(Image.open(img_name).convert('RGB'))\n(x0, x1, y0, y1) = get_roi(img)\nimg_roi = img[x0 : x1, y0 : y1]\nimg_flipped = img_roi.copy()[::-1].copy()\n# Segmented\nseg_name = os.path.join(root_dir, 'train', 'gts', index['0']['ImgName'].split('.')[0]+'.bmp')\nseg = np.array(Image.open(seg_name)).copy()\nseg_roi = seg[x0 : x1, y0 : y1]\nseg_flipped = seg_roi.copy()[::-1].copy()\n\nplt.subplot(2,3,1)\nplt.imshow(img, cmap=plt.cm.binary)\nplt.xlabel('Input image')\n\nplt.subplot(2,3,2)\nplt.imshow(img_roi, cmap=plt.cm.binary)\nplt.xlabel('Image centered on ROI')\n\nplt.subplot(2,3,3)\nplt.imshow(img_flipped, cmap=plt.cm.binary)\nplt.xlabel('Image flipped')\n\nplt.subplot(2,3,4)\nplt.imshow(seg, cmap=plt.cm.binary)\nplt.xlabel('Segmentation')\n\nplt.subplot(2,3,5)\nplt.imshow(seg_roi, cmap=plt.cm.binary)\nplt.xlabel('Segmentation ROI')\n\nplt.subplot(2,3,6)\nplt.imshow(seg_flipped, cmap=plt.cm.binary)\nplt.xlabel('Segmentation ROI flipped')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:11.706007Z","iopub.execute_input":"2021-05-23T16:03:11.706487Z","iopub.status.idle":"2021-05-23T16:03:16.647302Z","shell.execute_reply.started":"2021-05-23T16:03:11.706451Z","shell.execute_reply":"2021-05-23T16:03:16.646498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Dataset class + ROI","metadata":{}},{"cell_type":"code","source":"class RefugeDataset(Dataset):\n\n    def __init__(self, root_dir, split='train', output_size=(256,256)):\n        # Define attributes\n        self.output_size = output_size\n        self.root_dir = root_dir\n        self.split = split\n        TRESHOLD = 0.05\n\n        # Load data index\n        with open(os.path.join(self.root_dir, self.split, 'index.json')) as f:\n            self.index = json.load(f)\n        self.images = []\n        self.augmented = []\n        self.roi = []\n        self.augmented_index = []\n        self.length = len(self.index)\n\n        for k in range(self.length):\n            print('Loading {} image {}/{}...'.format(split, k, self.length), end='\\r')\n            img_name = os.path.join(self.root_dir, self.split, 'images', self.index[str(k)]['ImgName'])\n            img = np.array(Image.open(img_name).convert('RGB'))\n            roi = get_roi(img)\n            self.roi.append(roi)\n            img = img[roi[0]:roi[1], roi[2] : roi[3]]\n            img = transforms.functional.to_tensor(img)\n            img = transforms.functional.resize(img, self.output_size, interpolation=Image.BILINEAR)\n            self.images.append(img)\n        print(f'{len(self.augmented)} images augmented ')\n        self.images = self.images + self.augmented\n        self.augmented = []\n        \n            \n        # Load ground truth for 'train' and 'val' sets\n        if split != 'test':\n            self.segs = []\n            self.augmented_segs = []\n            for k in range(self.length):\n                print('Loading {} segmentation {}/{}...'.format(split, k, self.length), end='\\r')\n                seg_name = os.path.join(self.root_dir, self.split, 'gts', self.index[str(k)]['ImgName'].split('.')[0]+'.bmp')\n                seg = np.array(Image.open(seg_name)).copy()\n                roi = self.roi[k]\n                seg = seg[roi[0]:roi[1], roi[2] : roi[3]]\n                seg = 255. - seg\n                od = (seg>=127.).astype(np.float32)\n                oc = (seg>=250.).astype(np.float32)\n                od = torch.from_numpy(od[None,:,:])\n                oc = torch.from_numpy(oc[None,:,:])\n                od = transforms.functional.resize(od, self.output_size, interpolation=Image.NEAREST)\n                oc = transforms.functional.resize(oc, self.output_size, interpolation=Image.NEAREST)\n                seg = torch.cat([od, oc], dim=0)\n                self.segs.append(seg)\n            self.segs = self.segs + self.augmented_segs\n            self.augmented_segs = []\n        self.augmented_index = []\n        print('sizes : ', len(self.index), \"labels and\", len(self.images), 'images')\n\n        print('Succesfully loaded {} dataset.'.format(split) + ' '*50)\n    \n    def plot(self, items):\n        plt.figure(figsize=(10,10))\n        for i in range(25, 50):\n            plt.subplot(5,5,i-25+1)\n            plt.xticks([])\n            plt.yticks([])\n            plt.grid(False)\n            plt.imshow(items[i], cmap=plt.cm.binary)\n            plt.xlabel(f'Img n°{i}')\n        plt.show()\n    def __len__(self):\n        return len(self.index)\n\n    def __getitem__(self, idx):\n        # Image\n        img = self.images[idx]\n    \n        # Return only images for 'test' set\n        if self.split == 'test':\n            return img\n        \n        # Else, images and ground truth\n        else:\n            # Label\n            lab = torch.tensor(self.index[str(idx)]['Label'], dtype=torch.float32)\n            \n            # Segmentation masks\n            seg = self.segs[idx]\n\n            # Fovea localization\n            f_x = self.index[str(idx)]['Fovea_X']\n            f_y = self.index[str(idx)]['Fovea_Y']\n            fov = torch.FloatTensor([f_x, f_y])\n        \n            return img, lab, seg, fov, self.index[str(idx)]['ImgName']","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:16.648568Z","iopub.execute_input":"2021-05-23T16:03:16.648907Z","iopub.status.idle":"2021-05-23T16:03:16.669491Z","shell.execute_reply.started":"2021-05-23T16:03:16.648871Z","shell.execute_reply":"2021-05-23T16:03:16.668461Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Metrics","metadata":{}},{"cell_type":"code","source":"EPS = 1e-7\n\ndef compute_dice_coef(input, target):\n    '''\n    Compute dice score metric.\n    '''\n    batch_size = input.shape[0]\n    return sum([dice_coef_sample(input[k,:,:], target[k,:,:]) for k in range(batch_size)])/batch_size\n\ndef dice_coef_sample(input, target):\n    iflat = input.contiguous().view(-1)\n    tflat = target.contiguous().view(-1)\n    intersection = (iflat * tflat).sum()\n    return (2. * intersection) / (iflat.sum() + tflat.sum())\n\n\ndef vertical_diameter(binary_segmentation):\n    '''\n    Get the vertical diameter from a binary segmentation.\n    The vertical diameter is defined as the \"fattest\" area of the binary_segmentation parameter.\n    '''\n\n    # get the sum of the pixels in the vertical axis\n    vertical_axis_diameter = np.sum(binary_segmentation, axis=1)\n\n    # pick the maximum value\n    diameter = np.max(vertical_axis_diameter, axis=1)\n\n    # return it\n    return diameter\n\n\n\ndef vertical_cup_to_disc_ratio(od, oc):\n    '''\n    Compute the vertical cup-to-disc ratio from a given labelling map.\n    '''\n    # compute the cup diameter\n    cup_diameter = vertical_diameter(oc)\n    # compute the disc diameter\n    disc_diameter = vertical_diameter(od)\n\n    return cup_diameter / (disc_diameter + EPS)\n\ndef compute_vCDR_error(pred_od, pred_oc, gt_od, gt_oc):\n    '''\n    Compute vCDR prediction error, along with predicted vCDR and ground truth vCDR.\n    '''\n    pred_vCDR = vertical_cup_to_disc_ratio(pred_od, pred_oc)\n    gt_vCDR = vertical_cup_to_disc_ratio(gt_od, gt_oc)\n    vCDR_err = np.mean(np.abs(gt_vCDR - pred_vCDR))\n    return vCDR_err, pred_vCDR, gt_vCDR\n\n\ndef classif_eval(classif_preds, classif_gts):\n    '''\n    Compute AUC classification score.\n    '''\n    auc = roc_auc_score(classif_gts, classif_preds)\n    return auc\n\n\ndef fov_error(pred_fov, gt_fov):\n    '''\n    Fovea localization error metric (mean root squared error).\n    '''\n    err = np.sqrt(np.sum((gt_fov-pred_fov)**2, axis=1)).mean()\n    return err","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:16.670931Z","iopub.execute_input":"2021-05-23T16:03:16.671520Z","iopub.status.idle":"2021-05-23T16:03:16.684662Z","shell.execute_reply.started":"2021-05-23T16:03:16.671483Z","shell.execute_reply":"2021-05-23T16:03:16.683698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Post-processing functions","metadata":{}},{"cell_type":"code","source":"def refine_seg(pred):\n    '''\n    Only retain the biggest connected component of a segmentation map.\n    '''\n    np_pred = pred.numpy()\n        \n    largest_ccs = []\n    for i in range(np_pred.shape[0]):\n        labeled, ncomponents = label(np_pred[i,:,:])\n        bincounts = np.bincount(labeled.flat)[1:]\n        if len(bincounts) == 0:\n            largest_cc = labeled == 0\n        else:\n            largest_cc = labeled == np.argmax(bincounts)+1\n        largest_cc = torch.tensor(largest_cc, dtype=torch.float32)\n        largest_ccs.append(largest_cc)\n    largest_ccs = torch.stack(largest_ccs)\n    \n    return largest_ccs","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:16.685907Z","iopub.execute_input":"2021-05-23T16:03:16.686268Z","iopub.status.idle":"2021-05-23T16:03:16.697269Z","shell.execute_reply.started":"2021-05-23T16:03:16.686232Z","shell.execute_reply":"2021-05-23T16:03:16.696458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Network","metadata":{}},{"cell_type":"code","source":"class UNet(nn.Module):\n    def __init__(self, n_channels=3, n_classes=2):\n        super(UNet, self).__init__()\n        self.n_channels = n_channels\n        self.n_classes = n_classes\n        self.epoch = 0\n\n        self.inc = DoubleConv(n_channels, 64)\n        self.down1 = Down(64, 128)\n        self.down2 = Down(128, 256)\n        self.down3 = Down(256, 512)\n        factor = 2 \n        self.down4 = Down(512, 1024 // factor)\n        self.up1 = Up(1024, 512 // factor)\n        self.up2 = Up(512, 256 // factor)\n        self.up3 = Up(256, 128 // factor)\n        self.up4 = Up(128, 64)\n        self.output_layer = OutConv(64, n_classes)\n\n    def forward(self, x):\n        x1 = self.inc(x)\n        x2 = self.down1(x1)\n        x3 = self.down2(x2)\n        x4 = self.down3(x3)\n        x5 = self.down4(x4)\n        out = self.up1(x5, x4)\n        out = self.up2(out, x3)\n        out = self.up3(out, x2)\n        out = self.up4(out, x1)\n        out = self.output_layer(out)\n        out = torch.sigmoid(out)\n        return out\n\n    \nclass DoubleConv(nn.Module):\n    \"\"\"(convolution => [BN] => ReLU) * 2\"\"\"\n\n    def __init__(self, in_channels, out_channels, mid_channels=None):\n        super().__init__()\n        if not mid_channels:\n            mid_channels = out_channels\n        self.double_conv = nn.Sequential(\n            nn.Conv2d(in_channels, mid_channels, kernel_size=3, padding=1),\n            nn.BatchNorm2d(mid_channels),\n            nn.ReLU(inplace=True),\n            nn.Conv2d(mid_channels, out_channels, kernel_size=3, padding=1),\n            nn.BatchNorm2d(out_channels),\n            nn.ReLU(inplace=True)\n        )\n\n    def forward(self, x):\n        return self.double_conv(x)\n\n\nclass Down(nn.Module):\n    \"\"\"Downscaling with maxpool then double conv\"\"\"\n\n    def __init__(self, in_channels, out_channels):\n        super().__init__()\n        self.maxpool_conv = nn.Sequential(\n            nn.MaxPool2d(2),\n            DoubleConv(in_channels, out_channels)\n        )\n\n    def forward(self, x):\n        return self.maxpool_conv(x)\n\n\nclass Up(nn.Module):\n    \"\"\"Upscaling then double conv\"\"\"\n\n    def __init__(self, in_channels, out_channels):\n        super().__init__()\n\n        # Use the normal convolutions to reduce the number of channels\n        self.up = nn.Upsample(scale_factor=2, mode='bilinear', align_corners=True)\n        self.conv = DoubleConv(in_channels, out_channels, in_channels // 2)\n\n\n    def forward(self, x1, x2):\n        x1 = self.up(x1)\n        # input is CHW\n        diffY = x2.size()[2] - x1.size()[2]\n        diffX = x2.size()[3] - x1.size()[3]\n\n        x1 = F.pad(x1, [diffX // 2, diffX - diffX // 2,\n                        diffY // 2, diffY - diffY // 2])\n        x = torch.cat([x2, x1], dim=1)\n        return self.conv(x)\n\n\nclass OutConv(nn.Module):\n    '''\n    Simple convolution.\n    '''\n    def __init__(self, in_channels, out_channels):\n        super(OutConv, self).__init__()\n        self.conv = nn.Conv2d(in_channels, out_channels, kernel_size=1)\n\n    def forward(self, x):\n        return self.conv(x)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:16.698570Z","iopub.execute_input":"2021-05-23T16:03:16.698908Z","iopub.status.idle":"2021-05-23T16:03:16.718869Z","shell.execute_reply.started":"2021-05-23T16:03:16.698873Z","shell.execute_reply":"2021-05-23T16:03:16.718119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Settings","metadata":{}},{"cell_type":"code","source":"root_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\nlr = 1e-4\nbatch_size = 15\nnum_workers = 8\ntotal_epoch = 100","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:16.720017Z","iopub.execute_input":"2021-05-23T16:03:16.720393Z","iopub.status.idle":"2021-05-23T16:03:16.730611Z","shell.execute_reply.started":"2021-05-23T16:03:16.720358Z","shell.execute_reply":"2021-05-23T16:03:16.729893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Create Datasets and Data loaders","metadata":{}},{"cell_type":"code","source":"# Datasets\ntrain_set = RefugeDataset(root_dir, \n                          split='train')\nval_set = RefugeDataset(root_dir, \n                        split='val')\ntest_set = RefugeDataset(root_dir, \n                         split='test')\n\n# Dataloaders\ntrain_loader = DataLoader(train_set, \n                          batch_size=batch_size, \n                          shuffle=True, \n                          num_workers=num_workers,\n                          pin_memory=True,\n                         )\nval_loader = DataLoader(val_set, \n                        batch_size=batch_size, \n                        shuffle=False, \n                        num_workers=num_workers,\n                        pin_memory=True,\n                        )\ntest_loader = DataLoader(test_set, \n                        batch_size=batch_size, \n                        shuffle=False, \n                        num_workers=num_workers,\n                        pin_memory=True)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T16:03:16.731717Z","iopub.execute_input":"2021-05-23T16:03:16.732058Z","iopub.status.idle":"2021-05-23T17:03:19.013864Z","shell.execute_reply.started":"2021-05-23T16:03:16.732024Z","shell.execute_reply":"2021-05-23T17:03:19.012979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Device, model, loss and optimizer","metadata":{}},{"cell_type":"code","source":"# Device\ndevice = torch.device(\"cuda:0\")\n\n# Network\nmodel = UNet(n_channels=3, n_classes=2).to(device)\n\n# Loss\nseg_loss = torch.nn.BCELoss(reduction='mean')\n\n# Optimizer\noptimizer = optim.Adam(model.parameters(), lr=lr)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T17:03:19.015142Z","iopub.execute_input":"2021-05-23T17:03:19.015526Z","iopub.status.idle":"2021-05-23T17:03:22.860881Z","shell.execute_reply.started":"2021-05-23T17:03:19.015487Z","shell.execute_reply":"2021-05-23T17:03:22.860042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train for OC/OD segmentation","metadata":{}},{"cell_type":"code","source":"# Define parameters\nnb_train_batches = len(train_loader)\nnb_val_batches = len(val_loader)\nnb_iter = 0\nbest_val_auc = 0.\n\nwhile model.epoch < total_epoch:\n    # Accumulators\n    train_vCDRs, val_vCDRs = [], []\n    train_classif_gts, val_classif_gts = [], []\n    train_loss, val_loss = 0., 0.\n    train_dsc_od, val_dsc_od = 0., 0.\n    train_dsc_oc, val_dsc_oc = 0., 0.\n    train_vCDR_error, val_vCDR_error = 0., 0.\n    \n    ############\n    # TRAINING #\n    ############\n    model.train()\n    train_data = iter(train_loader)\n    for k in range(nb_train_batches):\n        # Loads data\n        imgs, classif_gts, seg_gts, fov_coords, names = train_data.next()\n        imgs, classif_gts, seg_gts = imgs.to(device), classif_gts.to(device), seg_gts.to(device)\n\n        # Forward pass\n        logits = model(imgs)\n        loss = seg_loss(logits, seg_gts)\n \n        # Backward pass\n        optimizer.zero_grad()\n        loss.backward()\n        optimizer.step()\n        train_loss += loss.item() / nb_train_batches\n        \n        with torch.no_grad():\n            # Compute segmentation metric\n            pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n            pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n            gt_od = seg_gts[:,0,:,:].type(torch.int8)\n            gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n            dsc_od = compute_dice_coef(pred_od, gt_od)\n            dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n            train_dsc_od += dsc_od.item()/nb_train_batches\n            train_dsc_oc += dsc_oc.item()/nb_train_batches\n\n\n            # Compute and store vCDRs\n            vCDR_error, pred_vCDR, gt_vCDR = compute_vCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n            train_vCDRs += pred_vCDR.tolist()\n            train_vCDR_error += vCDR_error / nb_train_batches\n            train_classif_gts += classif_gts.cpu().numpy().tolist()\n            \n        # Increase iterations\n        nb_iter += 1\n        \n        # Std out\n        print('Epoch {}, iter {}/{}, loss {:.6f}'.format(model.epoch+1, k+1, nb_train_batches, loss.item()) + ' '*20, \n              end='\\r')\n        \n    # Train a logistic regression on vCDRs\n    train_vCDRs = np.array(train_vCDRs).reshape(-1,1)\n    train_classif_gts = np.array(train_classif_gts)\n    clf = LogisticRegression(random_state=0, solver='lbfgs').fit(train_vCDRs, train_classif_gts)\n    train_classif_preds = clf.predict_proba(train_vCDRs)[:,1]\n    train_auc = classif_eval(train_classif_preds, train_classif_gts)\n    \n    ##############\n    # VALIDATION #\n    ##############\n    model.eval()\n    with torch.no_grad():\n        val_data = iter(val_loader)\n        for k in range(nb_val_batches):\n            # Loads data\n            imgs, classif_gts, seg_gts, fov_coords, names = val_data.next()\n            imgs, classif_gts, seg_gts = imgs.to(device), classif_gts.to(device), seg_gts.to(device)\n\n            # Forward pass\n            logits = model(imgs)\n            val_loss += seg_loss(logits, seg_gts).item() / nb_val_batches\n\n            # Std out\n            print('Validation iter {}/{}'.format(k+1, nb_val_batches) + ' '*50, \n                  end='\\r')\n            \n            # Compute segmentation metric\n            pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n            pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n            gt_od = seg_gts[:,0,:,:].type(torch.int8)\n            gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n            dsc_od = compute_dice_coef(pred_od, gt_od)\n            dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n            val_dsc_od += dsc_od.item()/nb_val_batches\n            val_dsc_oc += dsc_oc.item()/nb_val_batches\n            \n            # Compute and store vCDRs\n            vCDR_error, pred_vCDR, gt_vCDR = compute_vCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n            val_vCDRs += pred_vCDR.tolist()\n            val_vCDR_error += vCDR_error / nb_val_batches\n            val_classif_gts += classif_gts.cpu().numpy().tolist()\n            \n\n    # Glaucoma predictions from vCDRs\n    val_vCDRs = np.array(val_vCDRs).reshape(-1,1)\n    val_classif_gts = np.array(val_classif_gts)\n    val_classif_preds = clf.predict_proba(val_vCDRs)[:,1]\n    val_auc = classif_eval(val_classif_preds, val_classif_gts)\n        \n    # Validation results\n    print('VALIDATION epoch {}'.format(model.epoch+1)+' '*50)\n    print('LOSSES: {:.4f} (train), {:.4f} (val)'.format(train_loss, val_loss))\n    print('OD segmentation (Dice Score): {:.4f} (train), {:.4f} (val)'.format(train_dsc_od, val_dsc_od))\n    print('OC segmentation (Dice Score): {:.4f} (train), {:.4f} (val)'.format(train_dsc_oc, val_dsc_oc))\n    print('vCDR error: {:.4f} (train), {:.4f} (val)'.format(train_vCDR_error, val_vCDR_error))\n    print('Classification (AUC): {:.4f} (train), {:.4f} (val)'.format(train_auc, val_auc))\n    \n    # Save model if best validation AUC is reached\n    if val_auc > best_val_auc:\n        torch.save(model.state_dict(), '/kaggle/working/best_AUC_weights.pth')\n        with open('/kaggle/working/best_AUC_classifier.pkl', 'wb') as clf_file:\n            pickle.dump(clf, clf_file)\n        best_val_auc = val_auc\n        print('Best validation AUC reached. Saved model weights and classifier.')\n    print('_'*50)\n        \n    # End of epoch\n    model.epoch += 1","metadata":{"execution":{"iopub.status.busy":"2021-05-23T17:03:22.862126Z","iopub.execute_input":"2021-05-23T17:03:22.862456Z","iopub.status.idle":"2021-05-23T17:03:24.831174Z","shell.execute_reply.started":"2021-05-23T17:03:22.862422Z","shell.execute_reply":"2021-05-23T17:03:24.828507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load best model + classifier","metadata":{}},{"cell_type":"code","source":"# Load model and classifier\nmodel = UNet(n_channels=3, n_classes=2).to(device)\nmodel.load_state_dict(torch.load('/kaggle/working/best_AUC_weights.pth'))\nwith open('/kaggle/working/best_AUC_classifier.pkl', 'rb') as clf_file:\n    clf = pickle.load(clf_file)","metadata":{"execution":{"iopub.status.busy":"2021-05-23T17:03:24.832413Z","iopub.status.idle":"2021-05-23T17:03:24.832809Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check performance is maintained on validation","metadata":{}},{"cell_type":"code","source":"model.eval()\nval_vCDRs = []\nval_classif_gts = []\nval_loss = 0.\nval_dsc_od = 0.\nval_dsc_oc = 0.\nval_vCDR_error = 0.\nwith torch.no_grad():\n    val_data = iter(val_loader)\n    for k in range(nb_val_batches):\n        # Loads data\n        imgs, classif_gts, seg_gts, fov_coords, names = val_data.next()\n        imgs, classif_gts, seg_gts = imgs.to(device), classif_gts.to(device), seg_gts.to(device)\n\n        # Forward pass\n        logits = model(imgs)\n        val_loss += seg_loss(logits, seg_gts).item() / nb_val_batches\n\n        # Std out\n        print('Validation iter {}/{}'.format(k+1, nb_val_batches) + ' '*50, \n              end='\\r')\n\n        # Compute segmentation metric\n        pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n        pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n        gt_od = seg_gts[:,0,:,:].type(torch.int8)\n        gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n        dsc_od = compute_dice_coef(pred_od, gt_od)\n        dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n        val_dsc_od += dsc_od.item()/nb_val_batches\n        val_dsc_oc += dsc_oc.item()/nb_val_batches\n\n        # Compute and store vCDRs\n        vCDR_error, pred_vCDR, gt_vCDR = compute_vCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n        val_vCDRs += pred_vCDR.tolist()\n        val_vCDR_error += vCDR_error / nb_val_batches\n        val_classif_gts += classif_gts.cpu().numpy().tolist()\n\n\n# Glaucoma predictions from vCDRs\nval_vCDRs = np.array(val_vCDRs).reshape(-1,1)\nval_classif_gts = np.array(val_classif_gts)\nval_classif_preds = clf.predict_proba(val_vCDRs)[:,1]\nval_auc = classif_eval(val_classif_preds, val_classif_gts)\n\n# Validation results\nprint('VALIDATION '+' '*50)\nprint('LOSSES: {:.4f} (val)'.format(val_loss))\nprint('OD segmentation (Dice Score): {:.4f} (val)'.format(val_dsc_od))\nprint('OC segmentation (Dice Score): {:.4f} (val)'.format(val_dsc_oc))\nprint('vCDR error: {:.4f} (val)'.format(val_vCDR_error))\nprint('Classification (AUC): {:.4f} (val)'.format(val_auc))","metadata":{"execution":{"iopub.status.busy":"2021-05-23T17:03:24.833790Z","iopub.status.idle":"2021-05-23T17:03:24.834346Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Predictions on test set","metadata":{}},{"cell_type":"code","source":"nb_test_batches = len(test_loader)\nmodel.eval()\ntest_vCDRs = []\nwith torch.no_grad():\n    test_data = iter(test_loader)\n    for k in range(nb_test_batches):\n        # Loads data\n        imgs = test_data.next()\n        imgs = imgs.to(device)\n\n        # Forward pass\n        logits = model(imgs)\n\n        # Std out\n        print('Test iter {}/{}'.format(k+1, nb_test_batches) + ' '*50, \n              end='\\r')\n            \n        # Compute segmentation\n        pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n        pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n            \n        # Compute and store vCDRs\n        pred_vCDR = vertical_cup_to_disc_ratio(pred_od.cpu().numpy(), pred_oc.cpu().numpy())\n        test_vCDRs += pred_vCDR.tolist()\n            \n\n    # Glaucoma predictions from vCDRs\n    test_vCDRs = np.array(test_vCDRs).reshape(-1,1)\n    test_classif_preds = clf.predict_proba(test_vCDRs)[:,1]\n    \n# Prepare and save .csv file\ndef create_submission_csv(prediction, submission_filename='/kaggle/working/submission.csv'):\n    \"\"\"Create a sumbission file in the appropriate format for evaluation.\n\n    :param\n    prediction: list of predictions (ex: [0.12720, 0.89289, ..., 0.29829])\n    \"\"\"\n    \n    with open(submission_filename, mode='w') as csv_file:\n        fieldnames = ['Id', 'Predicted']\n        writer = csv.DictWriter(csv_file, fieldnames=fieldnames)\n        writer.writeheader()\n\n        for i, p in enumerate(prediction):\n            writer.writerow({'Id': \"T{:04d}\".format(i+1), 'Predicted': '{:f}'.format(p)})\n\ncreate_submission_csv(test_classif_preds)\n\n# The submission.csv file is under /kaggle/working/submission.csv.\n# If you want to submit it, you should download it before closing the current kernel.","metadata":{"execution":{"iopub.status.busy":"2021-05-23T17:03:24.835524Z","iopub.status.idle":"2021-05-23T17:03:24.836114Z"},"trusted":true},"execution_count":null,"outputs":[]}]}