{"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":"# Challenge 2 Team 9\n\nIn this challenge, to attempted to diagnose glaucoma from iris images.\nWe first analyse the data to see with what we are working.\nWe noticed that there is only a region from the images that is inerreting, thus we decided to create a crop algorithm that extract that regio of interest.\nWe did some work on the model already implemented by chaging some parameters ( number of epochs, changing proportion of train/val dataset, using the images from the crop algorithm).\n\nThen we wanted to do some Data Augmentation and Transfer Learning with the Resnet50 model.\n\nWe finally obtained a AUC Classification of 0.9815 thanks to a high pass filter.\n","metadata":{}},{"cell_type":"code","source":"import os\nimport json\nimport csv\nimport random\nimport pickle\nimport cv2\nimport numpy as np\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torch.optim as optim\nimport torchvision.transforms as transforms\n\n\nfrom PIL import Image, ImageOps\nfrom tensorflow.keras.preprocessing.image import ImageDataGenerator\nfrom torch.utils.data import Dataset, DataLoader\nfrom scipy.ndimage.measurements import label\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.neural_network import MLPClassifier\nfrom sklearn import svm\nfrom sklearn.metrics import roc_auc_score, roc_curve\n\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.datasets import make_classification\nfrom sklearn.tree import DecisionTreeClassifier\nfrom xgboost import XGBClassifier\nimport xgboost as xgb\nfrom sklearn.ensemble import ExtraTreesClassifier\n\nimport matplotlib.pyplot as plt\nimport tensorflow as tf\n\nfrom tensorflow.keras.applications import VGG19, Xception, ResNet50, MobileNet\nfrom keras.models import Sequential\nfrom keras.models import Model\nfrom keras.layers import Input \nimport keras\nfrom keras import layers\nfrom keras import models\nfrom keras import optimizers\nfrom keras.models import Sequential","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Image cropping\nThe idea is to crop the area of interest with this cropping algorithm","metadata":{}},{"cell_type":"code","source":"def Norme(p1,p2,p3,p4):\n        n = np.sqrt((p1[0]-p3[0])*(p1[0]-p3[0]) + (p2[0]-p4[0])*(p2[0]-p4[0]))\n        return n\n    \ndef contouring(img):\n    img1 = img.copy()\n    seuil = 25\n    colonne, ligne = img.size\n    for i in range(1,ligne-1):\n        for j in range(1,colonne-1):\n            p1 = img.getpixel((j-1,i))\n            p2 = img.getpixel((j,i-1))\n            p3 = img.getpixel((j+1,i))\n            p4 = img.getpixel((j,i+1))\n            n = Norme(p1,p2,p3,p4)\n            if n < seuil:\n                p = (255,255,255)\n            else:\n                p = (0,0,0)\n            img1.putpixel((j-1,i-1),p)\n    return(img1)\n\ndef crop_fovea(img,seg=None ) :\n    gray = ImageOps.grayscale(img)\n    width, height = img.size\n    arr = np.asarray(gray)\n    arg = arr.argmax()\n    x = arg%width\n    y = arg//width \n    step = 256\n\n\n    left = x-step\n    top = y-step\n    right = x+step\n    bottom = y+step\n\n    # Cropped image of above dimension\n    # (It will not change orginal image)\n    im1 = img.crop((left, top, right, bottom))\n    if seg : \n        seg1=seg.crop((left, top, right, bottom))\n        return(im1,seg1)\n    else : return(im1)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data analysis","metadata":{}},{"cell_type":"markdown","source":"First let's look at the dataset properties.\n\nHere we use the json file to better know the data, specially the proportion of glaucome/healthy eyes.","metadata":{}},{"cell_type":"code","source":"import pandas as pd\ndf_train = pd.read_json (r'/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/train/index.json')\ndf_train=df_train.T\ndf_train.head()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The Glaucoma diagnosis: a 1 for a glaucomatous eye, 0 for a healthy eye.\n\nLet analyses the difference between a glaucomatous eye and a healthy one.","metadata":{}},{"cell_type":"code","source":"df_train['Label'].value_counts().plot.pie()\nprint(\"number of healthy eyes :  \" , len(df_train[df_train[\"Label\"]==0]))\nprint(\"number of  glaucomatous eyes :  \" , len(df_train[df_train[\"Label\"]==1]))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It is worth to notice that there are much more healthy eyes and glaucomatous eyes. This means that the dataset is unbalanced and that this could lead to biaised models.\n\nWe may add also that the dataset is quite small ( only 400 datas).\n\nThe training, the validation and the test sets are all builts the same way( 400 images each)\n\nNow let's have a look at the images.","metadata":{}},{"cell_type":"code","source":"print(\"Glaucomatous eye: we see that the OD/OC area is big\")\nseg=\"/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/train/gts/g0001.bmp\"\nim=Image.open(r\"/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/train/images/g0001.jpg\")\n#im=im.convert('L')\nim=im.resize((256,256))\nim\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The segmentation files tells us where is the OD/OC area","metadata":{}},{"cell_type":"code","source":"seg=Image.open(seg)\nseg=seg.resize((256,256))\nim1arr = np.asarray(im)\nim2arr = np.asarray(seg)\n\n#addition = im1arr + im2arr\n#resultImage = Image.fromarray(addition)\n\n#resultImage\nseg.paste(im, (0,0), seg)\nseg\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"healthy eye: we see that the OD/OC area is smaller\")\nseg=\"/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/train/gts/n0089.bmp\"\nim=Image.open(r\"/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/train/images/n0089.jpg\")\n#im=im.convert('L')\nim=im.resize((256,256))\nim\nseg=Image.open(seg)\nseg=seg.resize((256,256))\nim1arr = np.asarray(im)\nim2arr = np.asarray(seg)\n\n#addition = (im1arr + im2arr)//2\n#resultImage = Image.fromarray(addition)\n\n#resultImage\nseg.paste(im, (0, 0), seg)\nseg","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here are the images from the crop algorithm","metadata":{}},{"cell_type":"code","source":"#images from the cropped algorithm\nseg=\"/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/train/gts/g0001.bmp\"\nim=\"/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/train/images/g0001.jpg\"\nseg=Image.open(seg)\nim=Image.open(im)\nim,seg=crop_fovea(im,seg)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"im=im.resize((128,128))\nim","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"seg=seg.resize((128,128))\nseg","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dataset class","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        \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            \n        self.images = []\n        for k in range(len(self.index)):\n            print('Loading {} image {}/{}...'.format(split, k, len(self.index)), 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            #img = Image.open(img_name).convert('RGB')\n            #img = img.filter(ImageFilter.Kernel((3, 3),\n                (-1, -2, -1, 0, 0, 0, 1, 2, 1), 0.1, 100)\n            #img = np.array(img)\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            \n        # Load ground truth for 'train' and 'val' sets\n        if split != 'test':\n            self.segs = []\n            for k in range(len(self.index)):\n                print('Loading {} segmentation {}/{}...'.format(split, k, len(self.index)), 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                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                \n        print('Succesfully loaded {} dataset.'.format(split) + ' '*50)\n            \n            \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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#algorithm to have the cropped images directly built in the dataset\nclass RefugeDatasetCrop(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        \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        if split == 'test':    \n            self.images = []\n            for k in range(len(self.index)):\n                print('Loading {} image {}/{}...'.format(split, k, len(self.index)), end='\\r')\n                img_name = os.path.join(self.root_dir, self.split, 'images', self.index[str(k)]['ImgName'])\n                img = Image.open(img_name)\n                img=np.array(crop_fovea(img).convert('RGB'))\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            \n        # Load ground truth for 'train' and 'val' sets\n        if split != 'test':\n            self.segs = []\n            self.images = []\n            for k in range(len(self.index)):\n                print('Loading {} image {}/{}...'.format(split, k, len(self.index)), end='\\r')\n                img_name = os.path.join(self.root_dir, self.split, 'images', self.index[str(k)]['ImgName'])\n                print('Loading {} segmentation {}/{}...'.format(split, k, len(self.index)), end='\\r')\n                seg_name = os.path.join(self.root_dir, self.split, 'gts', self.index[str(k)]['ImgName'].split('.')[0]+'.bmp')\n                img = Image.open(img_name)\n                seg=Image.open(seg_name)\n                img,seg=crop_fovea(img,seg)\n                img=np.array((img).convert('RGB'))\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                seg = np.array(seg).copy()\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                \n        print('Succesfully loaded {} dataset.'.format(split) + ' '*50)\n            \n            \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":{"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":{"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":{"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":{"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-6\nbatch_size = 16\nnum_workers = 8\ntotal_epoch = 100","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create datasets and data loaders\nAll image files are loaded in RAM in order to speed up the pipeline. Therefore, each dataset creation should take a few minutes.","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\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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We rearrange the train/val size proportion","metadata":{}},{"cell_type":"code","source":"#We split the dataset at 720/80 (90%)\n#be aware the first 40th are labelled 1 .\ndataset=train_set+val_set\ntrain_set1=[]\nval_set1=[]\nfor i in range(432) :\n    train_set1.append(dataset[i])\nfor i in range(432,512) :\n    val_set1.append(dataset[i])\nfor i in range(512,800) :\n    train_set1.append(dataset[i])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_loader1 = DataLoader(train_set1, \n                          batch_size=batch_size, \n                          shuffle=True, \n                          num_workers=num_workers,\n                          pin_memory=True,\n                         )\nval_loader1 = DataLoader(val_set1, \n                        batch_size=batch_size, \n                        shuffle=True, \n                        num_workers=num_workers,\n                        pin_memory=True,\n                        )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We create the cropped datasets","metadata":{}},{"cell_type":"code","source":"# Datasets\ntrain_set_crop = RefugeDatasetCrop(root_dir, \n                          split='train')\nval_set_crop = RefugeDatasetCrop(root_dir, \n                        split='val')\ntest_set_crop = RefugeDatasetCrop(root_dir, \n                         split='test')\n\n\n# Dataloaders\ntrain_loader_crop = DataLoader(train_set_crop, \n                          batch_size=batch_size, \n                          shuffle=True, \n                          num_workers=num_workers,\n                          pin_memory=True,\n                         )\nval_loader_crop = DataLoader(val_set_crop, \n                        batch_size=batch_size, \n                        shuffle=False, \n                        num_workers=num_workers,\n                        pin_memory=True,\n                        )\ntest_loader_crop = DataLoader(test_set_crop, \n                        batch_size=batch_size, \n                        shuffle=False, \n                        num_workers=num_workers,\n                        pin_memory=True)","metadata":{"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\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":{"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_crop)\nnb_val_batches = len(val_loader_crop)\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_crop)\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        # Forward pass\n        logits = model(imgs)\n        loss = seg_loss(logits, seg_gts)\n\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           \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().fit(train_vCDRs, train_classif_gts)\n    #clf = RandomForestClassifier(max_depth=2, random_state=0).fit(train_vCDRs, train_classif_gts)\n    #clf = XGBClassifier().fit(train_vCDRs, train_classif_gts)\n    #clf = ExtraTreesClassifier(n_estimators=100, random_state=0).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_crop)\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\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            \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\n        \n","metadata":{"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":{"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_crop)\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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Data Augmentation and Transfer Learning\n\nWe add some trubble using the datasets thus we had to reaarnage it first.","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport cv2\n\ntrain_classif_gts, val_classif_gts = [], []\n\nnb_train_batches = len(train_loader_crop)\nnb_val_batches = len(val_loader_crop)\n\nim_crop_train = []\nlabelled_train = []\n\n\ntrain_data = iter(train_loader_crop)\nfor k in range(nb_train_batches):\n    with torch.no_grad():\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        imgs = imgs.cpu()\n        train_classif_gts += classif_gts.cpu().numpy().tolist()\n        \n        for p in range (len(imgs[:,0,:,:])):\n            img = imgs[p,:,:,:].numpy()\n            labelled_train.append(classif_gts.cpu().numpy()[p])   \n            ima = np.zeros((256,256,3))\n            for i in range(3):\n                ima[:,:,i] = img[i,:,:]\n            im_crop_train.append(np.array(ima))\ntrain_classif_gts = np.array(train_classif_gts)\n\nim_crop_val= []\nlabelled_val = []\n\nwith torch.no_grad():\n    val_data = iter(val_loader_crop)\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        imgs = imgs.cpu()\n        val_classif_gts += classif_gts.cpu().numpy().tolist()\n        for p in range (len(imgs[:,0,:,:])):\n            img = imgs[p,:,:,:].numpy()\n            labelled_val.append(classif_gts.cpu().numpy()[p])  \n            ima = np.zeros((256,256,3))\n            for i in range(3):\n                ima[:,:,i] = img[i,:,:]\n            im_crop_val.append(np.array(ima))\n            \nval_classif_gts = np.array(val_classif_gts)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"labelled_train = np.array(labelled_train).reshape(-1,1)\nim_crop_train=np.array(im_crop_train)\nim_crop_train= im_crop_train.reshape(im_crop_train.shape[0], 256, 256, 3)\n\nlabelled_val = np.array(labelled_val).reshape(-1,1)\nim_crop_val=np.array(im_crop_val)\nim_crop_val = im_crop_val.reshape(im_crop_val.shape[0], 256, 256, 3)\n\ntrain_datagen = ImageDataGenerator(rotation_range = 40,width_shift_range=0.2,height_shift_range=0.2,\n                                  shear_range=0.2,zoom_range=0.2, brightness_range=[0.2,1.0])\n\ntrain_generator = train_datagen.flow(im_crop_train,labelled_train,batch_size=batch_size)\n\nval_datagen = ImageDataGenerator()\n\nval_generator =  val_datagen.flow(im_crop_val,labelled_val,batch_size=batch_size)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Resnet Model\n","metadata":{}},{"cell_type":"code","source":"from keras.applications.resnet50 import ResNet50","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"BATCH_SIZE = 8\nFC_LAYERS = [1024, 512, 256]\ndropout = 0.5\nNUM_EPOCHS = 100\n\n\ndef 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 = layers.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 = Model(inputs = base_model.input, outputs = preditions)\n    return finetune_model\n\nbase_model_1 = ResNet50(weights = 'imagenet',\n                       include_top = False,\n                       input_shape = (256, 256, 3))\n\n\nresnet50_model = build_model(base_model_1,\n                                      dropout = dropout,\n                                      fc_layers = FC_LAYERS,\n                                      num_classes = 2)\n\nadam = optimizers.Adam(lr = 1e-6)\nresnet50_model.compile(adam, loss=\"binary_crossentropy\", metrics=[\"accuracy\"])\n\n#the next commented line code allows to see the entire architecture of the model\n#keras.utils.plot_model(resnet50_model, \"multi_input_and_output_model.png\", show_shapes=True)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history = resnet50_model.fit(train_generator, validation_data = val_generator, epochs=15, verbose=1,batch_size=8,shuffle = True)\n                ","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(history.history['accuracy'])\nplt.plot(history.history['val_accuracy'])\nplt.title('model accuracy')\nplt.ylabel('accuracy')\nplt.xlabel('epoch')\nplt.legend(['train', 'test'], loc='upper left')\nplt.show()\n\nval_classif_preds = resnet50_model.predict(val_generator)\nval_preds=[]\nfor i in range(len(val_classif_preds)) :\n    val_preds.append(np.argmax(val_classif_preds[i]))\n\nval_auc = classif_eval(val_preds, labelled_val)\nprint(val_auc)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### VGG19 model\nhere is another model using transfer learning.","metadata":{}},{"cell_type":"code","source":"vgg19 = VGG19(weights='imagenet', include_top=False, input_shape=(60, 60, 3))\nfor layer in vgg19.layers[:]:\n    layer.trainable = False\nmodelVGG = Sequential()\nmodelVGG.add(vgg19)\nmodelVGG.add(layers.Flatten())\nmodelVGG.add(layers.Dense(1024, activation='relu'))\nmodelVGG.add(layers.Dense(1, activation='softmax'))\n\n\nmodelVGG.compile(loss='binary_crossentropy', optimizer=optimizers.Adam(lr=1e-4),metrics=[\"accuracy\"])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"historyVGG = modelVGG.fit(train_generator, validation_data = val_generator, epochs=15, verbose=1,batch_size=64)\nplt.plot(historyVGG.history['accuracy'])\nplt.plot(historyVGG.history['val_accuracy'])\nplt.title('model accuracy')\nplt.ylabel('accuracy')\nplt.xlabel('epoch')\nplt.legend(['train', 'test'], loc='upper left')\nplt.show()\n\nval_classif_predsVGG = modelVGG.predict(train_generator)\n\nval_aucVGG = classif_eval(val_classif_predsVGG, train_classif_gts)\nprint(val_aucVGG)","metadata":{"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        # Forward pass\n        logits = model(npimg)\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/submission700.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":{"trusted":true},"execution_count":null,"outputs":[]}]}