{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport 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\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\n\nimport matplotlib.pyplot as plt\n\nimport torchvision.transforms.functional as TF\nimport torchvision","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Original Dataloader","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 = 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":"markdown","source":"# Metrics","metadata":{}},{"cell_type":"code","source":"def 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 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    if np.ndim(binary_segmentation) == 2:\n        # get the sum of the pixels in the vertical axis\n        vertical_axis_diameter = np.sum(binary_segmentation, axis=1)\n        # pick the maximum value\n        diameter = np.max(vertical_axis_diameter)\n    else:\n        # get the sum of the pixels in the vertical axis\n        vertical_axis_diameter = np.sum(binary_segmentation, axis=1)\n        # pick the maximum value\n        diameter = np.max(vertical_axis_diameter, axis=1)\n\n    # return it\n    return diameter\n\ndef area(binary_segmentation):\n    '''\n    Get the area of the segmentation mask\n    '''\n    if np.ndim(binary_segmentation) == 2:\n        return np.sum(binary_segmentation)\n    else:\n        return np.sum(binary_segmentation, axis=(1,2))\n\n\ndef center(binary_segmentation):\n    '''\n    Get the coordinates of the center\n    '''\n    \n    center = np.nonzero(binary_segmentation)\n    c_center = int(np.mean(center[1]))  #row of the center\n    r_center = int(np.mean(center[0]))  #column of the center\n    \n    return np.array([c_center, r_center])\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\n\ndef area_cup_to_disc_ratio(od, oc):\n    '''\n    Compute the area cup-to-disc ratio from a given labelling map.\n    '''\n    ad = area(od)\n    ac = area(oc)\n    return ac / (ad + EPS)\n\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\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\ndef compute_aCDR_error(pred_od, pred_oc, gt_od, gt_oc):\n    '''\n    Compute aCDR prediction error, along with predicted aCDR and ground truth aCDR.\n    '''\n    pred_aCDR = area_cup_to_disc_ratio(pred_od, pred_oc)\n    gt_aCDR = area_cup_to_disc_ratio(gt_od, gt_oc)\n    aCDR_err = np.mean(np.abs(gt_aCDR - pred_aCDR))\n    return aCDR_err, pred_aCDR, gt_aCDR\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dataset Analysis\n\nIn this section we are going to explore the dataset, how the classes are distributed and the distribution of the feature we are interested in. ","metadata":{}},{"cell_type":"code","source":"# Load data index\nroot_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\n\nsplit = 'train'\nwith open(os.path.join(root_dir, split, 'index.json')) as f:\n    index_train = json.load(f)\n    \n#Class distribution in Training Dataset\nwith_glaucoma = 0\nfor k in range(len(index_train)):\n    if( index_train[str(k)]['Label'] == 1):\n        with_glaucoma = with_glaucoma + 1\nprint(f\"Train with Glaucoma:{with_glaucoma}/{len(index_train)} {with_glaucoma/len(index_train)*100}%\")\n\n\n#Class distribution in Validation Dataset\n\nsplit = 'val'\nwith open(os.path.join(root_dir, split, 'index.json')) as f:\n    index_val = json.load(f)\n\nwith_glaucoma = 0\nfor k in range(len(index_val)):\n    if( index_val[str(k)]['Label'] == 1):\n        with_glaucoma = with_glaucoma + 1\nprint(f\"Val with Glaucoma:{with_glaucoma}/{len(index_train)} {with_glaucoma/len(index_train)*100}%\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can notice that the class distribution is highly unbalanced, hence we will need to expand the number of minority samples.","metadata":{}},{"cell_type":"markdown","source":"Since as first alternative we are interested in the classification of Glaucoma based of the vCDR, we are going to use the groud truths masks to see how this feature is distributed based on the class.","metadata":{}},{"cell_type":"code","source":"from sklearn import preprocessing\nfrom sklearn.metrics import roc_auc_score, roc_curve","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_workers = 8\nbatch_size = 8\n\n# Datasets\ntrain_set = RefugeDataset(root_dir, \n                          split='train')\nval_set = RefugeDataset(root_dir, \n                        split='val')\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                        )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Store as array\ntrain_od = np.zeros((len(train_set.segs)))\ntrain_oc = np.zeros((len(train_set.segs)))\ntrain_vCDR = np.zeros((len(train_set.segs)))\n\ntrain_ad = np.zeros((len(train_set.segs))) #area disc\ntrain_ac = np.zeros((len(train_set.segs))) #area cup\ntrain_aCDR = np.zeros((len(train_set.segs)))\n\ntrain_label = np.zeros((len(train_set.segs)))\n\nfor cnt, element in enumerate(train_set.segs):\n    train_od[cnt] = vertical_diameter(element.numpy()[0, :, :])\n    train_oc[cnt] = vertical_diameter(element.numpy()[1, :, :])\n    train_vCDR[cnt] = vertical_cup_to_disc_ratio(element.numpy()[0, :, :], element.numpy()[1, :, :])\n    train_label[cnt] = train_set.index[str(cnt)]['Label'] \n    \n    train_ad[cnt] = area(element.numpy()[0,:,:])\n    train_ac[cnt] = area(element.numpy()[1,:,:])\n    train_aCDR[cnt] = area_cup_to_disc_ratio(element.numpy()[0, :, :],element.numpy()[1, :, :])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looking at the distribution of the training set, the two distributions are mostly divisible with a good degree of precision.","metadata":{}},{"cell_type":"code","source":"glaucoma_train = (train_label == 1)\nno_glaucoma_train = (train_label == 0)\n\nbins = 20\nplt.hist(train_vCDR[glaucoma_train], bins, alpha=0.8, density=True, stacked=True,  label='glaucoma')\nplt.hist(train_vCDR[no_glaucoma_train], bins, alpha=0.5, density=True, stacked=True, label='no glaucoma')\nplt.xlabel(\"Vertical Cup-Disk Ratio\")\nplt.legend(loc='upper right')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.scatter(train_aCDR[glaucoma_train],train_vCDR[glaucoma_train])\nplt.scatter(train_aCDR[no_glaucoma_train], train_vCDR[no_glaucoma_train])\nplt.xlabel(\"Area Cup-Disk Ratio\")\nplt.ylabel(\"Vertical Cup-Disk Ratio\")\n\n# Train a logistic regression on vCDRs\ntrain_vCDRs = np.array(train_vCDR).reshape(-1,1)\ntrain_aCDRs = np.array(train_aCDR).reshape(-1,1)\nX = np.concatenate((train_vCDRs, train_aCDRs),axis=1)\ntrain_classif_gts = np.array(train_label)\nclf = LogisticRegression(random_state=0, solver='lbfgs').fit(X, train_classif_gts)\ntrain_classif_preds = clf.predict_proba(X)[:,1]\ntrain_auc = classif_eval(train_classif_preds, train_classif_gts)\nprint(f\"Train_AUC (Logistic Regression): {train_auc}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Store as array\nval_od = np.zeros((len(val_set.segs)))\nval_oc = np.zeros((len(val_set.segs)))\nval_vCDR = np.zeros((len(val_set.segs)))\nval_label = np.zeros((len(val_set.segs)))\n\nval_ad = np.zeros((len(val_set.segs)))\nval_ac = np.zeros((len(val_set.segs)))\nval_aCDR = np.zeros((len(val_set.segs)))\n\nval_dist_cdc = np.zeros((len(val_set.segs))) #centers distance\n\nfor cnt, element in enumerate(val_set.segs):\n    val_od[cnt] = vertical_diameter(element.numpy()[0, :, :])\n    val_oc[cnt] = vertical_diameter(element.numpy()[1, :, :])\n    val_vCDR[cnt] = vertical_cup_to_disc_ratio(element.numpy()[0, :, :], element.numpy()[1, :, :])\n    val_label[cnt] = val_set.index[str(cnt)]['Label'] \n    \n    val_ad[cnt] = area(element.numpy()[0,:,:])\n    val_ac[cnt] = area(element.numpy()[1,:,:])\n    val_aCDR[cnt] = area_cup_to_disc_ratio(element.numpy()[0, :, :],element.numpy()[1, :, :])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"val_label.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, the validation set does not have the same distribution of the training set. The distribution of the two classes are mostly overlapping when using as feature only the vCDR. That's why we thought to introduce a second feature, the Area Cup-Disk Ratio (aCDR) in order divide the two distributions. The problem that comes with the use of such feature is that, due to the segmentation process, the area estimation is much more sensitive to errors.","metadata":{}},{"cell_type":"code","source":"glaucoma_val = (val_label == 1)\nno_glaucoma_val = (val_label == 0)\n\nbins = 20\nplt.hist(val_vCDR[glaucoma_val], bins, alpha=0.8, density=True, stacked=True,  label='glaucoma')\nplt.hist(val_vCDR[no_glaucoma_val], bins, alpha=0.5, density=True, stacked=True, label='no glaucoma')\nplt.xlabel(\"Vertical Cup-Disk Ratio\")\nplt.legend(loc='upper right')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.scatter(val_aCDR[glaucoma_val], val_vCDR[glaucoma_val])\nplt.scatter(val_aCDR[no_glaucoma_val], val_vCDR[no_glaucoma_val])\n\n# Train a logistic regression on vCDRs\nval_vCDRs = np.array(val_vCDR).reshape(-1,1)\nval_aCDRs = np.array(val_aCDR).reshape(-1,1)\nX = np.concatenate((val_vCDRs, val_aCDRs),axis=1)\nval_classif_gts = np.array(val_label)\nval_classif_preds = clf.predict_proba(X)[:,1]\nval_auc = classif_eval(val_classif_preds, val_classif_gts)\nprint(f\"Validation_AUC (Logistic Regression): {val_auc}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Images for report","metadata":{}},{"cell_type":"code","source":"bins = 20\n\nfig, axs = plt.subplots(2,1,figsize=(10,10))\n\naxs[0].hist(train_vCDR[no_glaucoma_train], bins, alpha=0.5, density=True, stacked=True, label='no glaucoma')\naxs[0].hist(train_vCDR[glaucoma_train], bins, alpha=0.8, density=True, stacked=True,  label='glaucoma')\naxs[0].set_xlabel(\"Vertical Cup-Disk Ratio\")\naxs[0].legend(loc='upper right')\naxs[0].set_title(\"Training set\")\n\naxs[1].hist(val_vCDR[no_glaucoma_val], bins, alpha=0.5, density=True, stacked=True, label='no glaucoma')\naxs[1].hist(val_vCDR[glaucoma_val], bins, alpha=0.8, density=True, stacked=True,  label='glaucoma')\naxs[1].set_xlabel(\"Vertical Cup-Disk Ratio\")\naxs[1].legend(loc='upper right')\naxs[1].set_title(\"Validation set\")\n\nfig.savefig('/kaggle/working/distribution_vDCR.pdf', format='pdf')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(2,1,figsize=(10,10))\n\naxs[0].scatter(train_aCDR[no_glaucoma_train], train_vCDR[no_glaucoma_train], label = \"no glaucoma\")\naxs[0].scatter(train_aCDR[glaucoma_train],train_vCDR[glaucoma_train], label='glaucoma')\naxs[0].set_xlabel(\"Area Cup-Disk Ratio\")\naxs[0].set_ylabel(\"Vertical Cup-Disk Ratio\")\naxs[0].set_title(\"Training set\")\naxs[0].legend(loc='lower right')\n\naxs[1].scatter(val_aCDR[no_glaucoma_val], val_vCDR[no_glaucoma_val], label='no glaucoma')\naxs[1].scatter(val_aCDR[glaucoma_val], val_vCDR[glaucoma_val], label='glaucoma')\naxs[1].set_xlabel(\"Area Cup-Disk Ratio\")\naxs[1].set_ylabel(\"Vertical Cup-Disk Ratio\")\naxs[1].set_title(\"Validation set\")\naxs[1].legend(loc='lower right')\n\nfig.savefig('/kaggle/working/distribution_aDCR_vCDR.pdf', format='pdf')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dataset Manipulation and Minority Augmentation\n\nIn this section we are going to manipulate the dataset in order to balance the minority class and have better images of the disk area that we are interested in.","metadata":{}},{"cell_type":"markdown","source":"## Image Cropping","metadata":{}},{"cell_type":"code","source":"root_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\nsplit = 'train'\nout_size = (256,256)\n\nwith open(os.path.join(root_dir, split, 'index.json')) as f:\n    index = json.load(f)\n    \nk = 103 # Select a random Image\n\n_, axs = plt.subplots(1,3,figsize=(10, 10))\n\n#Load Image\nprint('Loading {} image {}/{}...'.format(split, k, len(index)), end='\\r')\nimg_name = os.path.join(root_dir, split, 'images', index[str(k)]['ImgName'])\nimg = np.array(Image.open(img_name).convert('RGB'))\n\nchannels = ['Red','Green','Blue']\nfor i,channel in enumerate(channels):\n    axs[i].imshow(img[:,:,i])\n    axs[i].axis('off')\n    axs[i].set_xlabel(channel)\n\n#We assume that the cup is on the left side of the image and crop the upper\n#part since there could be strange reflection due to photography\nimg_g = img[:,:,1]\nh, w = img_g.shape\nimg_g = img_g[int(h/4):,:int(w/2)]\n\n#Thresholding, since the cup is the most luminous point of the image\nmask = (img_g > np.max(img_g)*0.9).astype(np.float)\ncenter = np.nonzero(mask)\nx_center = np.mean(center[1]).astype(int)\ny_center = np.mean(center[0]).astype(int)+int(h/4)\n\nfig, axs = plt.subplots(1,2, figsize=(10, 5))\naxs[0].imshow(img)\naxs[0].axis('off')\naxs[0].scatter(x_center, y_center, marker=\".\", s=500)\n\n#Crop on the cup\nif x_center > 300:\n    img_crop = img[y_center-300:y_center+300, x_center-300:x_center+300]\nelse:\n    img_crop = img[y_center-300:y_center+300, 0:x_center+300]\naxs[1].imshow(img_crop)\naxs[1].axis('off')\n\nfig.savefig('/kaggle/working/cropping.pdf', format='pdf')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The following pipeline will calculate the cup locations of all disks in the dataset. The entire cropping pipeline has been implemented on in our Dataloader.","metadata":{}},{"cell_type":"code","source":"def cup_center(img):\n    \n    #We can assume that the cup is on the left side of the image\n    #We will use the green component of the image since it has the\n    #best contrast between all channels\n    img_g = img[:,:,1]\n    h, w = img_g.shape\n    img_g = img_g[int(h/4):,:int(w/2)]\n\n    #smoothing to remove reflection noise\n    kernel = np.ones((31,31),np.float32)/31**2\n    img_g = cv2.filter2D(img_g,-1,kernel)\n\n    #Thresholding to keep only brightest part \n    mask = (img_g > np.max(img_g)*0.9).astype(np.float)\n    center = np.nonzero(mask)\n    c_center = int(np.mean(center[1]))  #row of the center\n    r_center = int(np.mean(center[0])) + int(h/4)  #column of the center\n\n    return r_center, c_center\n\ndef cup_crop(img, center, max_img_side):\n    \n    x_center = center[1]\n    y_center = center[0]\n\n    #Crop on the cup\n    l = int(max_img_side/2)\n    if x_center > l:\n        img_crop = img[y_center-l:y_center+l, x_center-l:x_center+l]\n    else:\n        img_crop = img[y_center-l:y_center+l, 0:x_center+l]\n\n    return img_crop\n\ndef sharpening(img):\n    # Create kernel\n    kernel = np.array([[-1, -1, -1], \n                   [-1,  9, -1], \n                   [-1, -1, -1]])\n\n    # Sharpen image\n    return cv2.filter2D(img, -1, kernel)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pickle\n\nroot_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\nsplits = ['train','val','test']\n\ncenters = {}\nfor split in splits:\n    \n    centers[split] = {}\n    with open(os.path.join(root_dir, split, 'index.json')) as f:\n        index = json.load(f)\n\n    for k in range(len(index)):\n        #Load Image\n        print('Loading {} image {}/{}...'.format(split, k, len(index)), end='\\r')\n        img_name = os.path.join(root_dir, split, 'images', index[str(k)]['ImgName'])\n        img = np.array(Image.open(img_name).convert('RGB'))\n        r_center, c_center = cup_center(img)\n        \n        centers[split][k] = {}\n        centers[split][k][\"Cup_X\"] = c_center\n        centers[split][k][\"Cup_Y\"] = r_center","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dirName = f\"/kaggle/working/centers/\"\nif not os.path.exists(dirName):\n    os.makedirs(dirName)\n    print(\"Directory \" , dirName ,  \" Created \")\nelse:    \n    print(\"Directory \" , dirName ,  \" already exists\") \n        \n        \nfor split in splits:\n    json_name = f\"/kaggle/working/centers/{split}.json\"\n    with open(json_name, 'w') as fp:\n        json.dump(centers[split], fp)\n        print(f\"{split}.json saved.\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Minority Class Augmentation\n\nGiven the highly unbalanced class distribution, we thought of introducing an image manipulation pipeline in the Dataloader in order to create replicas of the minority samples introducing rotations, scaling and Random Colour Jitter.","metadata":{}},{"cell_type":"code","source":"import torchvision.transforms.functional as TF\nfrom torchvision.transforms import ColorJitter\nfrom torchvision.transforms import RandomAffine\n\nroot_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\ncenters_dir =  '/kaggle/input/centers'\nsplit = 'train'\nout_size = (256,256)\n\nwith open(os.path.join(root_dir, split, 'index.json')) as f:\n    index = json.load(f)\n\nwith open(os.path.join(centers_dir, f\"{split}.json\")) as f:\n    centers = json.load(f)\n            \nk = 103 # Select a random Image\n\n#Load Image\nimg_name = os.path.join(root_dir, split, 'images', index[str(k)]['ImgName'])\nimg = np.array(Image.open(img_name).convert('RGB'))\n\ncup_x = centers[str(k)][\"Cup_X\"]\ncup_y = centers[str(k)][\"Cup_Y\"]\nimg = cup_crop(img,(cup_y,cup_x),500)\n\nimg = transforms.functional.to_tensor(img)\n\n\n_, axs = plt.subplots(1,4,figsize=(30, 10))\ntmp = (img.numpy()*255).astype(np.uint8)\ntmp = np.moveaxis(tmp,[0,1,2],[-1,-3,-2])\naxs[0].imshow(tmp)\n\n\nra =  RandomAffine(10, scale=(1,1.3))\ncj = ColorJitter(brightness = (0.8,1.2),contrast = (0.8,1.2) , saturation = (0.8,1.2), hue = (-0.01,0.01))\ntmp2 = ra(img)\ntmp2 = cj(tmp2)\ntmp2 = (tmp2.numpy()*255).astype(np.uint8)\ntmp2 = np.moveaxis(tmp2,[0,1,2],[-1,-3,-2])\naxs[1].imshow(tmp2)\n\ntmp3 = TF.rotate(img, 6)\ntmp3 = (tmp3.numpy()*255).astype(np.uint8)\ntmp3 = np.moveaxis(tmp3,[0,1,2],[-1,-3,-2])\naxs[2].imshow(tmp3)\n\ntmp3 = TF.adjust_contrast(img, 1.5)\ntmp3 = (tmp3.numpy()*255).astype(np.uint8)\ntmp3 = np.moveaxis(tmp3,[0,1,2],[-1,-3,-2])\naxs[2].imshow(tmp3)\n\ntmp4 = TF.adjust_contrast(img, 1.5)\ntmp4 = TF.adjust_gamma(tmp4, 1.5)\ntmp4 = (tmp4.numpy()*255).astype(np.uint8)\ntmp4 = np.moveaxis(tmp4,[0,1,2],[-1,-3,-2])\naxs[3].imshow(tmp4)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Polar Transformation","metadata":{}},{"cell_type":"code","source":"_, axs = plt.subplots(1,2,figsize=(30, 10))\ntmp = (img.numpy()*255).astype(np.uint8)\ntmp = np.moveaxis(tmp,[0,1,2],[-1,-3,-2])\naxs[0].imshow(tmp)\n\nvalue = np.sqrt(((img.shape[-1]/2.0)**2.0)+((img.shape[-2]/2.0)**2.0))\ntmp4 = (img.numpy()*255).astype(np.uint8)\ntmp4 = np.moveaxis(tmp4,[0,1,2],[-1,-3,-2])\ntmp = cv2.linearPolar(tmp4,(tmp4.shape[0]/2, tmp4.shape[1]/2), value*0.7, cv2.WARP_FILL_OUTLIERS)\naxs[1].imshow(tmp)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Custom Dataloader\nWe have modified the original dataloder to accomodate all the discussed transformations in order to have an on-the-flight loading of the modified images.","metadata":{}},{"cell_type":"code","source":"from torchvision.transforms import ColorJitter\nimport torchvision.transforms.functional as TF\n\nclass RefugeDataset(Dataset):\n\n    def __init__(self, root_dir, working_dir, split='train', output_size=(256,256), minority_augmentation = False, polar_transform = False):\n        # Define attributes\n        self.output_size = output_size\n        self.root_dir = root_dir\n        self.working_dir = working_dir\n        self.split = split\n        \n        self.rotations = [-10, -5, 5, 10]\n        self.zoom = [1.2, 1.3]\n        brightness = (0.8,1.2)\n        contrast = (0.8,1.2) \n        saturation = (0.8,1.2)\n        hue = (-0.05,0.05)\n        self.cj = ColorJitter(brightness = brightness,contrast = contrast , saturation = saturation, hue = hue)\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        with open(os.path.join(self.working_dir, f\"{self.split}.json\")) as f:\n            self.centers = json.load(f)\n            \n        self.images = []\n        self.new_index = {}\n        j = 0\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            \n            #We need to insert here the function to modify the images\n            \n            #Crop Image on Cup\n            cup_x = self.centers[str(k)][\"Cup_X\"]\n            cup_y = self.centers[str(k)][\"Cup_Y\"]\n            img = cup_crop(img,(cup_y,cup_x),500)\n            \n            #Sharpening of the Image\n            img = sharpening(img)\n            \n            if polar_transform:\n                value = np.sqrt(((img.shape[0]/2.0)**2.0)+((img.shape[1]/2.0)**2.0))\n                img = cv2.linearPolar(img,(img.shape[0]/2, img.shape[1]/2), value*0.7, cv2.WARP_FILL_OUTLIERS)\n            \n            img = transforms.functional.to_tensor(img)\n            img = transforms.functional.resize(img, self.output_size, interpolation=Image.BILINEAR)\n            original = img\n            \n            #Adjust Contrast\n            img = TF.adjust_contrast(img, 1.5)\n            #img = TF.adjust_gamma(img, 1.5)\n            \n            self.images.append(img)\n            self.new_index[str(j)] = self.index[str(k)]\n            j= j + 1\n            \n            if minority_augmentation and (self.index[str(k)]['Label'] == 1):\n                for rotation in self.rotations:\n                    for zoom in self.zoom:\n                        img = original\n\n                        img = TF.rotate(img,rotation)\n                        img = TF.affine(img,angle = 0, translate = (0,0), shear = 0, scale = zoom)\n                        img = self.cj(img)\n                        \n                        #Adjust Contrast\n                        img = TF.adjust_contrast(img, 1.5)\n                        #img = TF.adjust_gamma(img, 1.5)\n                                                \n                        self.images.append(img)\n                        self.new_index[str(j)] = self.index[str(k)]\n                        j= j + 1\n                        \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                \n                #We need to insert here the function to modify the images\n                cup_x = self.centers[str(k)][\"Cup_X\"]\n                cup_y = self.centers[str(k)][\"Cup_Y\"]\n                seg = cup_crop(seg,(cup_y,cup_x),500)\n                \n                if polar_transform:\n                    value = np.sqrt(((seg.shape[0]/2.0)**2.0)+((seg.shape[1]/2.0)**2.0))\n                    seg = cv2.linearPolar(seg,(seg.shape[0]/2, seg.shape[1]/2), value*0.7, cv2.WARP_FILL_OUTLIERS)\n                \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                \n                    \n                original_seg = seg\n                self.segs.append(seg)\n                \n                if minority_augmentation and (self.index[str(k)]['Label'] == 1):\n                    for rotation in self.rotations:\n                        for zoom in self.zoom:\n                            seg = original_seg\n\n                            seg = TF.rotate(seg,rotation)\n                            seg = TF.affine(seg,angle = 0, translate = (0,0), shear = 0, scale = zoom)\n\n                            self.segs.append(seg)\n                \n                \n        self.index = self.new_index        \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":"# Datasets\n\nroot_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\nworking_dir = '/kaggle/input/centers'\ntrain_set = RefugeDataset(root_dir, working_dir, split='train', minority_augmentation = True, polar_transform = False)\nprint(f\"Dataset: {len(train_set.index)}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"batch_size = 8\nnum_workers = 8\n\ntrain_loader = DataLoader(train_set, \n                          batch_size=batch_size, \n                          shuffle=True, \n                          num_workers=num_workers,\n                          pin_memory=True,\n                         )\ntrain_data = iter(train_loader)\nimgs, classif_gts, seg_gts, fov_coords, names = train_data.next()\nprint(classif_gts)\n\n_, axs = plt.subplots(2,len(imgs),figsize=(30, 10))\nfor cnt,(img,seg) in enumerate(zip(imgs,seg_gts)):\n    tmp = (img.numpy()*255).astype(np.uint8)\n    tmp = np.moveaxis(tmp,[0,1,2],[-1,-3,-2])\n    axs[0,cnt].imshow(tmp)\n    axs[0,cnt].axis('off')\n    axs[1,cnt].imshow(seg[1] + seg[0]*0.5)\n    axs[1,cnt].axis('off')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(2, 2,figsize=(8, 8))\naxs = axs.flatten()\nfor cnt,(img,seg) in enumerate(zip(imgs,seg_gts)):\n    if cnt == 4: break\n    tmp = (img.numpy()*255).astype(np.uint8)\n    tmp = np.moveaxis(tmp,[0,1,2],[-1,-3,-2])\n    axs[cnt].imshow(tmp)\n    axs[cnt].axis('off')\n    \nfig.savefig('/kaggle/working/sample_augmentation.pdf', format='pdf')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Datasets\n\nroot_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\nworking_dir = '/kaggle/input/centers'\ntrain_set = RefugeDataset(root_dir, working_dir, split='train', minority_augmentation = True, polar_transform = True)\nprint(f\"Dataset: {len(train_set.index)}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"batch_size = 8\nnum_workers = 8\n\ntrain_loader = DataLoader(train_set, \n                          batch_size=batch_size, \n                          shuffle=True, \n                          num_workers=num_workers,\n                          pin_memory=True,\n                         )\ntrain_data = iter(train_loader)\nimgs, classif_gts, seg_gts, fov_coords, names = train_data.next()\nprint(classif_gts)\n\n_, axs = plt.subplots(2,len(imgs),figsize=(30, 10))\nfor cnt,(img,seg) in enumerate(zip(imgs,seg_gts)):\n    tmp = (img.numpy()*255).astype(np.uint8)\n    tmp = np.moveaxis(tmp,[0,1,2],[-1,-3,-2])\n    axs[0,cnt].imshow(tmp)\n    axs[0,cnt].axis('off')\n    axs[1,cnt].imshow(seg[1] + seg[0]*0.5)\n    axs[1,cnt].axis('off')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Classification Using Segmentation\nIn this section we are going to explore the classification of Glaucoma using features crafted through segmentation.","metadata":{}},{"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: UNet","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'\nworking_dir = '/kaggle/input/centers'\nlr = 1e-4\nbatch_size = 20\nnum_workers = 8\ntotal_epoch = 100","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Dataset and Dataloaders","metadata":{}},{"cell_type":"code","source":"# Datasets\ntrain_set = RefugeDataset(root_dir, working_dir,\n                          split='train', minority_augmentation = True )\nval_set = RefugeDataset(root_dir, working_dir,\n                        split='val', minority_augmentation = False)\ntest_set = RefugeDataset(root_dir, working_dir,\n                         split='test', minority_augmentation = False)\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                         )\n\nval_loader = DataLoader(val_set, \n                        batch_size=batch_size, \n                        shuffle=False, \n                        num_workers=num_workers,\n                        pin_memory=True,\n                        )\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":"### Device, Model, Loss and Optimizer","metadata":{}},{"cell_type":"code","source":"import torch\nimport torch.nn as nn\n\ndef dice_loss(pred, target, smooth = 1.):\n    pred = pred.contiguous()\n    target = target.contiguous()    \n\n    intersection = (pred * target).sum(dim=2).sum(dim=2)\n    \n    loss = (1 - ((2. * intersection + smooth) / (pred.sum(dim=2).sum(dim=2) + target.sum(dim=2).sum(dim=2) + smooth)))\n    \n    return loss.mean()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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# From the results, the use of dice loss creates more irregular contours than BCE\n#seg_loss = dice_loss\n\n# Optimizer\noptimizer = optim.Adam(model.parameters(), lr=lr)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Training Loop","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"use_area = False\n\n# 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 = [], [] #accumulators vertical ratios\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    if use_area:\n        train_aCDR_error, val_aCDR_error = 0., 0.\n        train_aCDRs, val_aCDRs = [], [] #accumulators area ratios\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            if use_area:\n                #Compute and store aCDRs\n                aCDR_error, pred_aCDR, gt_aCDR = compute_aCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n                train_aCDRs += pred_aCDR.tolist()\n                train_aCDR_error += aCDR_error / nb_train_batches\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    if use_area:\n        train_aCDRs = np.array(train_aCDRs).reshape(-1,1)\n        X = np.concatenate((train_vCDRs, train_aCDRs),axis=1)\n    else:\n        X = train_vCDRs\n        \n    train_classif_gts = np.array(train_classif_gts)\n    clf = LogisticRegression(random_state=0, solver='lbfgs').fit(X, train_classif_gts)\n    train_classif_preds = clf.predict_proba(X)[:,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            if use_area:\n                #Compute and store aCDRs\n                aCDR_error, pred_aCDR, gt_aCDR = compute_aCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n                val_aCDRs += pred_aCDR.tolist()\n                val_aCDR_error += aCDR_error / nb_val_batches\n            \n\n    # Glaucoma predictions from vCDRs\n    val_vCDRs = np.array(val_vCDRs).reshape(-1,1)\n    if use_area:\n        val_aCDRs = np.array(val_aCDRs).reshape(-1,1)\n        X = np.concatenate((val_vCDRs, val_aCDRs),axis=1)\n    else:\n        X = val_vCDRs\n        \n    val_classif_gts = np.array(val_classif_gts)\n    val_classif_preds = clf.predict_proba(X)[:,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    if use_area:\n        print('aCDR error: {:.4f} (train), {:.4f} (val)'.format(train_aCDR_error, val_aCDR_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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.max(train_vCDRs)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Segmentation Results, Validation and Submission","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":"### Segmentation Result","metadata":{}},{"cell_type":"code","source":"val_data = iter(val_loader)\nwith torch.no_grad():\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    logits = model(imgs)\n    pred_ods = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu())\n    pred_ocs = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu())\n\n    _, axs = plt.subplots(len(imgs),3, figsize = (10,100))\n    for cnt,(img, seg_gt, pred_od, pred_oc) in enumerate(zip(imgs.cpu(),seg_gts.cpu(),pred_ods, pred_ocs)):\n        tmp = (img.numpy()*255).astype(np.uint8)\n        tmp = np.moveaxis(tmp,[0,1,2],[-1,-3,-2])\n        axs[cnt,0].imshow(tmp)\n        axs[cnt,0].axis('off')\n        axs[cnt,0].set_title(\"Sample\")\n\n        axs[cnt,1].imshow(seg_gt[1] + seg_gt[0]*0.5)\n        axs[cnt,1].axis('off')\n        axs[cnt,1].set_title(\"Sample\")\n        \n        axs[cnt,2].imshow(pred_od + pred_oc*0.5)\n        axs[cnt,2].axis('off')\n        axs[cnt,2].set_title(\"Estimated\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(1,3, figsize = (10,5))\nfor cnt,(img, seg_gt, pred_od, pred_oc) in enumerate(zip(imgs.cpu(),seg_gts.cpu(),pred_ods, pred_ocs)):\n    tmp = (img.numpy()*255).astype(np.uint8)\n    tmp = np.moveaxis(tmp,[0,1,2],[-1,-3,-2])\n    axs[0].imshow(tmp)\n    axs[0].axis('off')\n    axs[0].set_title(\"Sample\")\n\n    axs[1].imshow(seg_gt[1] + seg_gt[0]*0.5)\n    axs[1].axis('off')\n    axs[1].set_title(\"Sample\")\n\n    axs[2].imshow(pred_od + pred_oc*0.5)\n    axs[2].axis('off')\n    axs[2].set_title(\"Estimated\")Lo\n    break\n    \nfig.savefig('/kaggle/working/segmentation_results.pdf', format='pdf')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Checking Performances on Validation","metadata":{}},{"cell_type":"code","source":"model.eval()\nval_vCDRs, val_aCDRs = [], []\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        #Compute and store aCDRs\n        if use_area:\n            val_aCDRs += area_cup_to_disc_ratio(ad, ac).tolist()\n\n\n# Glaucoma predictions from vCDRs\nval_vCDRs = np.array(val_vCDRs).reshape(-1,1)\nif use_area:\n    val_aCDRs = np.array(val_aCDRs).reshape(-1,1)\n    X = np.concatenate((val_vCDRs, val_aCDRs),axis=1)\nelse:\n    X = val_vCDRs\nval_classif_gts = np.array(val_classif_gts)\nval_classif_preds = clf.predict_proba(X)[:,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":"### Results on Test: Submission","metadata":{}},{"cell_type":"code","source":"nb_test_batches = len(test_loader)\nmodel.eval()\ntest_vCDRs, test_aCDRs = [], []\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        #Compute and store aCDRs\n        ad = area(pred_od.cpu().numpy())\n        ac = area(pred_oc.cpu().numpy())\n        test_aCDRs += area_cup_to_disc_ratio(ad, ac).tolist()\n\n    # Glaucoma predictions from vCDRs\n    test_vCDRs = np.array(test_vCDRs).reshape(-1,1)\n    val_aCDRs = np.array(val_aCDRs).reshape(-1,1)\n    X = np.concatenate((val_vCDRs, val_aCDRs),axis=1)\n    test_classif_preds = clf.predict_proba(X)[:,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_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Classification with ResNet","metadata":{}},{"cell_type":"markdown","source":"### Dataset and Settings","metadata":{}},{"cell_type":"code","source":"root_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\nlr = 1e-4\nbatch_size = 8\nnum_workers = 8\ntotal_epochs = 100\nworking_dir = '/kaggle/input/centers'\noutput_size = (64, 64) # since we are working with pretrained resnet","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Datasets\ntrain_set = RefugeDataset(root_dir, working_dir,\n                          split='train', minority_augmentation = True, polar_transform = False )\nval_set = RefugeDataset(root_dir, working_dir,\n                        split='val', minority_augmentation = False, polar_transform = False )\ntest_set = RefugeDataset(root_dir, working_dir,\n                         split='test', minority_augmentation = False, polar_transform = False )\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":"### Model Selector","metadata":{}},{"cell_type":"code","source":"from torch.optim import lr_scheduler","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def modelSelector(model_sel):\n    if (model_sel == \"resnet\"):\n        model = torchvision.models.resnet50(pretrained=True).to(device)\n        for param in model.parameters():\n            param.requires_grad = False\n        model.fc = torch.nn.Linear(model.fc.in_features, 1).to(device)\n        model = torch.nn.Sequential(\n        model,\n        torch.nn.Sigmoid()\n        )\n    elif (model_sel == \"squeezenet\"):\n        model = torchvision.models.squeezenet1_1(pretrained=True).to(device)\n        for param in model.features.parameters():\n            param.requires_grad = False\n        model.classifier[1] = torch.nn.Conv2d(512, 1, kernel_size=(1, 1), stride=(1, 1)).to(device)\n        model = torch.nn.Sequential(\n        model,\n        torch.nn.Sigmoid()\n        )\n    elif (model_sel == \"vgg16\"):\n        model = torchvision.models.vgg16_bn(pretrained=True).to(device)\n        for param in model.features.parameters():\n            param.requires_grad = False\n        model.classifier[6] = torch.nn.Linear(4096, 1).to(device)\n        model = torch.nn.Sequential(\n        model,\n        torch.nn.Sigmoid()\n        )\n    elif (model_sel == \"vgg11\"):\n        model = torchvision.models.vgg11_bn(pretrained=True).to(device)\n        for param in model.features.parameters():\n            param.requires_grad = False\n        model.classifier[6] = torch.nn.Linear(4096, 1).to(device)\n        model = torch.nn.Sequential(\n        model,\n        torch.nn.Sigmoid()\n        )\n    return model\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport json\nimport csv\nimport random\nimport pickle\nimport cv2\nimport numpy as np\nimport copy\n\nimport torch, torchvision\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.svm import SVC\nfrom sklearn.metrics import roc_auc_score, roc_curve","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"device = torch.device(\"cuda:0\")\nmodel = modelSelector(\"vgg11\")\noptimizer = optim.Adam(model.parameters(), lr=lr)\ncriterion = torch.nn.BCELoss()\nscheduler = lr_scheduler.StepLR(optimizer, step_size=7, gamma=0.1)\ndataloaders = {'train': train_loader, 'val': val_loader}\ndataset_sizes = {'train': len(train_set), 'val': len(val_set)}","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(model[0].classifier[0].in_features)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"best_model_wts = copy.deepcopy(model.state_dict())\nbest_acc = 0.0\n\n\nfor epoch in range(total_epochs):\n    print('Epoch {}/{}'.format(epoch, total_epochs - 1))\n    print('-' * 10)\n\n    # Each epoch has a training and validation phase\n    for phase in ['train', 'val']:\n        if phase == 'train':\n            model.train()  # Set model to training mode\n        else:\n            model.eval()   # Set model to evaluate mode\n\n        running_loss = 0.0\n        running_corrects = 0\n\n        # Iterate over data.\n        first_iter = True\n        for inputs, labels, seg, fov, img_name in dataloaders[phase]:\n            inputs = inputs.to(device)\n            labels = labels.flatten().to(device)\n            \n            # zero the parameter gradients\n            optimizer.zero_grad()\n\n            # forward\n            # track history if only in train\n            with torch.set_grad_enabled(phase == 'train'):\n                outputs = model(inputs)\n                #preds = torch.sigmoid(outputs)\n                preds=outputs\n                loss = criterion(outputs.flatten(), labels)\n                \n                output_cpu = preds.cpu().detach().numpy()\n                labels_epoch_cpu = labels.cpu().detach().numpy()\n                if (first_iter):\n                    epoch_outputs = output_cpu\n                    labels_epoch = labels_epoch_cpu\n                    first_iter = False\n                else:\n                    epoch_outputs = np.concatenate((epoch_outputs, output_cpu), axis=0)\n                    labels_epoch = np.concatenate((labels_epoch, labels_epoch_cpu), axis=0)\n                \n                # backward + optimize only if in training phase\n                if phase == 'train':\n                    loss.backward()\n                    optimizer.step()\n\n            # statistics\n            running_loss += loss.item() * inputs.size(0)\n            running_corrects += torch.sum(torch.round(preds) == labels.data)\n        if phase == 'train':\n            scheduler.step()\n        if phase == \"val\":\n            pass\n            #print(np.unique(epoch_outputs))\n        epoch_loss = running_loss / dataset_sizes[phase]\n        epoch_acc = roc_auc_score(labels_epoch, epoch_outputs)\n\n        print('{} Loss: {:.4f} ROC: {:.4f}'.format(\n            phase, epoch_loss, epoch_acc))\n\n        # deep copy the model\n        if phase == 'val' and epoch_acc > best_acc:\n            best_acc = epoch_acc\n            best_model_wts = copy.deepcopy(model.state_dict())\n\n    print()\n\nprint('Best val Acc: {:4f}'.format(best_acc))\n\n# load best model weights\nmodel.load_state_dict(best_model_wts)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.load_state_dict(best_model_wts)\n\nnb_test_batches = len(test_loader)\nmodel.eval()\ntest_vCDRs, test_aCDRs = [], []\nwith torch.no_grad():\n    test_data = iter(test_loader)\n    first_iter = True\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        preds=logits   \n        output_cpu = preds.cpu().detach().numpy()\n        if (first_iter):\n            epoch_outputs = output_cpu\n            first_iter = False\n        else:\n            epoch_outputs = np.concatenate((epoch_outputs, output_cpu), axis=0)\n\n\n    test_classif_preds = epoch_outputs\n\nprint(\"\\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(float(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_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ensamble Method: Ensamble of Deep Neural Network","metadata":{}},{"cell_type":"markdown","source":"## Polar and Cartesian coordinate Ensemble","metadata":{}},{"cell_type":"code","source":"from torchvision.transforms import ColorJitter\nimport torchvision.transforms.functional as TF\n\nclass RefugeDataset_Polar_Regular(Dataset):\n\n    def __init__(self, root_dir, working_dir, split='train', output_size=(256,256), minority_augmentation = False):\n        # Define attributes\n        self.output_size = output_size\n        self.root_dir = root_dir\n        self.working_dir = working_dir\n        self.split = split\n        \n        self.rotations = [-10, -5, 5, 10]\n        self.zoom = [1.2, 1.3]\n        brightness = (0.8,1.2)\n        contrast = (0.8,1.2) \n        saturation = (0.8,1.2)\n        hue = (-0.05,0.05)\n        self.cj = ColorJitter(brightness = brightness,contrast = contrast , saturation = saturation, hue = hue)\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        with open(os.path.join(self.working_dir, f\"{self.split}.json\")) as f:\n            self.centers = json.load(f)\n            \n        self.images_regular = []\n        self.images_polar = []\n        self.new_index = {}\n        j = 0\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            \n            #We need to insert here the function to modify the images\n            \n            #Crop Image on Cup\n            cup_x = self.centers[str(k)][\"Cup_X\"]\n            cup_y = self.centers[str(k)][\"Cup_Y\"]\n            img = cup_crop(img,(cup_y,cup_x),500)\n            \n            #Sharpening of the Image\n            img = sharpening(img)\n            \n            \n            value = np.sqrt(((img.shape[0]/2.0)**2.0)+((img.shape[1]/2.0)**2.0))\n            img_polar = cv2.linearPolar(img,(img.shape[0]/2, img.shape[1]/2), value*0.7, cv2.WARP_FILL_OUTLIERS)\n            img_polar = transforms.functional.to_tensor(img_polar)\n            img_polar = transforms.functional.resize(img_polar, self.output_size, interpolation=Image.BILINEAR)\n            original_polar = img_polar\n            \n            img = transforms.functional.to_tensor(img)\n            img = transforms.functional.resize(img, self.output_size, interpolation=Image.BILINEAR)\n            original = img\n            \n            #Adjust Contrast\n            img = TF.adjust_contrast(img, 1.5)\n            img_polar = TF.adjust_contrast(img_polar, 1.5)\n            \n            self.images_regular.append(img)\n            self.images_polar.append(img_polar)\n            self.new_index[str(j)] = self.index[str(k)]\n            j= j + 1\n            \n            if minority_augmentation and (self.index[str(k)]['Label'] == 1):\n                for rotation in self.rotations:\n                    for zoom in self.zoom:\n                        \n                        img = original\n                        img = TF.rotate(img,rotation)\n                        img = TF.affine(img,angle = 0, translate = (0,0), shear = 0, scale = zoom)\n                        img = self.cj(img)\n                        \n                        #Adjust Contrast\n                        img = TF.adjust_contrast(img, 1.5)\n                                                \n                        self.images_regular.append(img)\n                        \n                        img_polar = original_polar\n                        img_polar = TF.rotate(img_polar,rotation)\n                        img_polar = TF.affine(img_polar,angle = 0, translate = (0,0), shear = 0, scale = zoom)\n                        img_polar = self.cj(img_polar)\n                        \n                        #Adjust Contrast\n                        img_polar = TF.adjust_contrast(img_polar, 1.5)\n                                                \n                        self.images_polar.append(img_polar)\n                                                \n                        \n                        self.new_index[str(j)] = self.index[str(k)]\n                        j= j + 1\n                                \n                \n        self.index = self.new_index        \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_regular = self.images_regular[idx]\n        img_polar = self.images_polar[idx]\n    \n        # Return only images for 'test' set\n        if self.split == 'test':\n            return img_regular, img_polar\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            # 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_regular, img_polar, lab, fov, self.index[str(idx)]['ImgName']","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import torch.nn as nn\n\nclass Ensamble(nn.Module):\n    def __init__(self, device):\n        \n        super(Ensamble,self).__init__()\n        self.core_regular = torchvision.models.vgg11_bn(pretrained=True).to(device)\n        self.core_polar =  torchvision.models.vgg11_bn(pretrained=True).to(device)\n        \n        in_features_classifier = self.core_regular.classifier[0].in_features * 2    \n        \n        for param in self.core_regular.features.parameters():\n            param.requires_grad = False\n        self.core_regular.classifier = torch.nn.Identity()\n        \n        \n        \n        for param in self.core_polar.features.parameters():\n            param.requires_grad = False  \n        self.core_polar.classifier = torch.nn.Identity()\n            \n        self.classifier = torch.nn.Sequential(\n            torch.nn.Linear(in_features_classifier, 4096, bias = True),\n            torch.nn.ReLU(inplace=True),\n            torch.nn.Dropout(p = 0.5, inplace = False),\n            torch.nn.Linear(4096, 4096, bias = True),\n            torch.nn.ReLU(inplace=True),\n            torch.nn.Dropout(p = 0.5, inplace = False),\n            torch.nn.Linear(4096, 1, bias=True),\n            torch.nn.Sigmoid()\n        ).to(device)\n        \n    def forward(self,img_regular, img_polar):\n        features_regular = self.core_regular(img_regular)\n        features_polar = self.core_polar(img_polar)\n        \n        features = torch.cat((features_regular, features_polar), -1)\n        \n        return self.classifier(features)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"root_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\nlr = 1e-4\nbatch_size = 8\nnum_workers = 8\ntotal_epochs = 100\nworking_dir = '/kaggle/input/centers'\noutput_size = (64, 64) # since we are working with pretrained resnet","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Datasets\ntrain_set = RefugeDataset_Polar_Regular(root_dir, working_dir, \n                            split='train', minority_augmentation = True)\nval_set = RefugeDataset_Polar_Regular(root_dir, working_dir,\n                        split='val', minority_augmentation = False)\ntest_set = RefugeDataset_Polar_Regular(root_dir, working_dir,\n                         split='test', minority_augmentation = False)\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":"code","source":"train_loader = DataLoader(train_set, \n                          batch_size=batch_size, \n                          shuffle=True, \n                          num_workers=num_workers,\n                          pin_memory=True,\n                         )\nimgs_reg, imgs_pol, classif_gts, fov_coords, names = iter(train_loader).next()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img_pol.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"_, axs = plt.subplots(2,len(imgs_reg),figsize=(30, 10))\nfor cnt,(img_reg,img_pol) in enumerate(zip(imgs_reg,imgs_pol)):\n    tmp = (img_reg.numpy()*255).astype(np.uint8)\n    tmp = np.moveaxis(tmp,[0,1,2],[-1,-3,-2])\n    axs[0,cnt].imshow(tmp)\n    axs[0,cnt].axis('off')\n    \n    tmp = (img_pol.numpy()*255).astype(np.uint8)\n    tmp = np.moveaxis(tmp,[0,1,2],[-1,-3,-2])\n    axs[1,cnt].imshow(tmp)\n    axs[1,cnt].axis('off')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"device = torch.device(\"cuda:0\")\nmodel = Ensamble(device)\noptimizer = optim.Adam(model.parameters(), lr=lr)\ncriterion = torch.nn.BCELoss()\nscheduler = lr_scheduler.StepLR(optimizer, step_size=7, gamma=0.1)\ndataloaders = {'train': train_loader, 'val': val_loader}\ndataset_sizes = {'train': len(train_set), 'val': len(val_set)}","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"best_model_wts = copy.deepcopy(model.state_dict())\nbest_acc = 0.0\n\n\nfor epoch in range(total_epochs):\n    print('Epoch {}/{}'.format(epoch, total_epochs - 1))\n    print('-' * 10)\n\n    # Each epoch has a training and validation phase\n    for phase in ['train', 'val']:\n        if phase == 'train':\n            model.train()  # Set model to training mode\n        else:\n            model.eval()   # Set model to evaluate mode\n\n        running_loss = 0.0\n        running_corrects = 0\n\n        # Iterate over data.\n        first_iter = True\n        for inputs_reg, inputs_pol, labels, fov, img_name in dataloaders[phase]:\n            inputs_reg = inputs_reg.to(device)\n            inputs_pol = inputs_pol.to(device)\n            labels = labels.flatten().to(device)\n            \n            # zero the parameter gradients\n            optimizer.zero_grad()\n\n            # forward\n            # track history if only in train\n            with torch.set_grad_enabled(phase == 'train'):\n                outputs = model(inputs_reg, inputs_pol)\n                #preds = torch.sigmoid(outputs)\n                preds=outputs\n                loss = criterion(outputs.flatten(), labels)\n                \n                output_cpu = preds.cpu().detach().numpy()\n                labels_epoch_cpu = labels.cpu().detach().numpy()\n                if (first_iter):\n                    epoch_outputs = output_cpu\n                    labels_epoch = labels_epoch_cpu\n                    first_iter = False\n                else:\n                    epoch_outputs = np.concatenate((epoch_outputs, output_cpu), axis=0)\n                    labels_epoch = np.concatenate((labels_epoch, labels_epoch_cpu), axis=0)\n                \n                # backward + optimize only if in training phase\n                if phase == 'train':\n                    loss.backward()\n                    optimizer.step()\n\n            # statistics\n            running_loss += loss.item() * inputs_reg.size(0)\n            running_corrects += torch.sum(torch.round(preds) == labels.data)\n        if phase == 'train':\n            scheduler.step()\n        if phase == \"val\":\n            pass\n            #print(np.unique(epoch_outputs))\n        epoch_loss = running_loss / dataset_sizes[phase]\n        epoch_acc = roc_auc_score(labels_epoch, epoch_outputs)\n\n        print('{} Loss: {:.4f} ROC: {:.4f}'.format(\n            phase, epoch_loss, epoch_acc))\n\n        # deep copy the model\n        if phase == 'val' and epoch_acc > best_acc:\n            best_acc = epoch_acc\n            best_model_wts = copy.deepcopy(model.state_dict())\n\n    print()\n\nprint('Best val Acc: {:4f}'.format(best_acc))\n\n# load best model weights\nmodel.load_state_dict(best_model_wts)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Crop and no Crop Ensemble","metadata":{}},{"cell_type":"code","source":"class Ensamble(nn.Module):\n    def __init__(self, device):\n        \n        super(Ensamble,self).__init__()\n        self.core_crop = torchvision.models.vgg11_bn(pretrained=True).to(device)\n        self.core_nocrop =  torchvision.models.vgg11_bn(pretrained=True).to(device)\n        \n        in_features_classifier = self.core_crop.classifier[0].in_features * 2    \n        \n        for param in self.core_crop.features.parameters():\n            param.requires_grad = False\n        self.core_crop.classifier = torch.nn.Identity()\n        \n        \n        \n        for param in self.core_nocrop.features.parameters():\n            param.requires_grad = False  \n        self.core_nocrop.classifier = torch.nn.Identity()\n            \n        self.classifier = torch.nn.Sequential(\n            torch.nn.Linear(in_features_classifier, 4096, bias = True),\n            torch.nn.ReLU(inplace=True),\n            torch.nn.Dropout(p = 0.5, inplace = False),\n            torch.nn.Linear(4096, 4096, bias = True),\n            torch.nn.ReLU(inplace=True),\n            torch.nn.Dropout(p = 0.5, inplace = False),\n            torch.nn.Linear(4096, 1, bias=True),\n            torch.nn.Sigmoid()\n        ).to(device)\n        \n    def forward(self,img_crop, img_nocrop):\n        features_crop = self.core_crop(img_crop)\n        features_nocrop = self.core_nocrop(img_nocrop)\n        \n        features = torch.cat((features_crop, features_nocrop), -1)\n        \n        return self.classifier(features)\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"root_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\nlr = 1e-4\nbatch_size = 8\nnum_workers = 8\ntotal_epochs = 50\nworking_dir = '/kaggle/working/centers'\noutput_size = (256, 256) # since we are working with pretrained resnet\naugmented = True\n\n# Datasets\ntrain_set_nocrop = RefugeDataset_Original_Cropped(root_dir, working_dir, \n                            split='train', minority_augmentation = True)\nval_set_nocrop = RefugeDataset_Original_Cropped(root_dir, working_dir,\n                        split='val', minority_augmentation = False)\ntest_set_nocrop = RefugeDataset_Original_Cropped(root_dir, working_dir,\n                         split='test', minority_augmentation = False)\n\n# Dataloaders\ntrain_loader_nocrop = DataLoader(train_set_nocrop, \n                          batch_size=batch_size, \n                          shuffle=True, \n                          num_workers=num_workers,\n                          pin_memory=True,\n                         )\nval_loader_nocrop = DataLoader(val_set_nocrop, \n                        batch_size=batch_size, \n                        shuffle=False, \n                        num_workers=num_workers,\n                        pin_memory=True,\n                        )\ntest_loader_nocrop = DataLoader(test_set_nocrop, \n                        batch_size=batch_size, \n                        shuffle=False, \n                        num_workers=num_workers,\n                        pin_memory=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def modelSelector(model_sel):\n    if (model_sel == \"resnet\"):\n        model = torchvision.models.resnet18(pretrained=True).to(device)\n    elif (model_sel == \"squeezenet\"):\n        model = torchvision.models.squeezenet1_1(pretrained=True).to(device)\n    elif (model_sel == \"vgg16\"):\n        model = torchvision.models.vgg16(pretrained=True).to(device)\n    elif (model_sel == \"vgg16_bn\"):\n        model = torchvision.models.vgg16_bn(pretrained=True).to(device)\n    elif (model_sel == \"vgg11\"):\n        model = torchvision.models.vgg11(pretrained=True).to(device)\n    elif (model_sel == \"vgg11_bn\"):\n        model = torchvision.models.vgg11_bn(pretrained=True).to(device)\n\n    return model","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for selected_model in [\"vgg11\", \"vgg11_bn\", \"vgg16\", \"vgg16_bn\"]:\n    device = torch.device(\"cuda:0\")\n    model = Ensamble(device, modelSelector(selected_model))\n    optimizer = optim.Adam(model.parameters(), lr=lr)\n    criterion = torch.nn.BCELoss()\n    scheduler = lr_scheduler.StepLR(optimizer, step_size=7, gamma=0.1)\n    dataloaders = {'train': train_loader_nocrop, 'val': val_loader_nocrop}\n    dataset_sizes = {'train': len(train_set_nocrop), 'val': len(val_set_nocrop)}\n    best_model_wts = copy.deepcopy(model.state_dict())\n    best_acc = 0.0\n\n\n    for epoch in range(total_epochs):\n        print('Epoch {}/{}'.format(epoch, total_epochs - 1))\n        print('-' * 10)\n\n        # Each epoch has a training and validation phase\n        for phase in ['train', 'val']:\n            if phase == 'train':\n                model.train()  # Set model to training mode\n            else:\n                model.eval()   # Set model to evaluate mode\n\n            running_loss = 0.0\n            running_corrects = 0\n\n            # Iterate over data.\n            first_iter = True\n            for inputs_reg, inputs_pol, labels, fov, img_name in dataloaders[phase]:\n                inputs_reg = inputs_reg.to(device)\n                inputs_pol = inputs_pol.to(device)\n                labels = labels.flatten().to(device)\n\n                # zero the parameter gradients\n                optimizer.zero_grad()\n\n                # forward\n                # track history if only in train\n                with torch.set_grad_enabled(phase == 'train'):\n                    outputs = model(inputs_reg, inputs_pol)\n                    #preds = torch.sigmoid(outputs)\n                    preds=outputs\n                    loss = criterion(outputs.flatten(), labels)\n\n                    output_cpu = preds.cpu().detach().numpy()\n                    labels_epoch_cpu = labels.cpu().detach().numpy()\n                    if (first_iter):\n                        epoch_outputs = output_cpu\n                        labels_epoch = labels_epoch_cpu\n                        first_iter = False\n                    else:\n                        epoch_outputs = np.concatenate((epoch_outputs, output_cpu), axis=0)\n                        labels_epoch = np.concatenate((labels_epoch, labels_epoch_cpu), axis=0)\n\n                    # backward + optimize only if in training phase\n                    if phase == 'train':\n                        loss.backward()\n                        optimizer.step()\n\n                # statistics\n                running_loss += loss.item() * inputs_reg.size(0)\n                running_corrects += torch.sum(torch.round(preds) == labels.data)\n            if phase == 'train':\n                scheduler.step()\n            if phase == \"val\":\n                pass\n                #print(np.unique(epoch_outputs))\n            epoch_loss = running_loss / dataset_sizes[phase]\n            epoch_acc = roc_auc_score(labels_epoch, epoch_outputs)\n\n            print('{} Loss: {:.4f} ROC: {:.4f}'.format(\n                phase, epoch_loss, epoch_acc))\n\n            # deep copy the model\n            if phase == 'val' and epoch_acc > best_acc:\n                best_acc = epoch_acc\n                best_model_wts = copy.deepcopy(model.state_dict())\n\n        print()\n\n    print(selected_model + \" finished training\")\n    print('Best val Acc: {:4f}'.format(best_acc))\n\n    \n    \n\n    # load best model weights\n    model.load_state_dict(best_model_wts)\n    \n    nb_test_batches = len(test_loader_nocrop)\n    model.eval()\n    test_vCDRs, test_aCDRs = [], []\n    with torch.no_grad():\n        test_data = iter(test_loader_nocrop)\n        first_iter = True\n        for k in range(nb_test_batches):\n            # Loads data\n            img_orig, img_crop = test_data.next()\n            img_orig = img_orig.to(device)\n            img_crop = img_crop.to(device)\n\n            # Forward pass\n            logits = model(img_orig, img_crop)\n\n            # Std out\n            print('Test iter {}/{}'.format(k+1, nb_test_batches) + ' '*50, \n                  end='\\r')\n            preds=logits   \n            output_cpu = preds.cpu().detach().numpy()\n            if (first_iter):\n                epoch_outputs = output_cpu\n                first_iter = False\n            else:\n                epoch_outputs = np.concatenate((epoch_outputs, output_cpu), axis=0)\n\n\n        test_classif_preds = epoch_outputs\n\n    print(\"\\n\")\n    # Prepare and save .csv file\n    def 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(float(p))})\n    filename = 'submission_' + selected_model + '.csv'\n    create_submission_csv(test_classif_preds, submission_filename=filename)","metadata":{},"execution_count":null,"outputs":[]}]}