{"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:: Medical image segmentation\n\n## Team 23\n\n**Paul-Friedrich Dugué**\n\n**Pierre Piovesan**","metadata":{}},{"cell_type":"markdown","source":"## I) Overview of the process","metadata":{}},{"cell_type":"markdown","source":"General imports","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\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","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The inputs are retinal fundus images from normal and glaucomatous eyes.","metadata":{}},{"cell_type":"markdown","source":"Visualization of the input.","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\nimage = plt.imread(\"/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/train/images/g0001.jpg\")\nplt.title(\"Retinal fundus image 2124x2056\")\nplt.imshow(image)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As the resolution of the image is to high, it is recommended to resize, and more generally pre-process, the images.","metadata":{}},{"cell_type":"code","source":"resized = cv2.resize(image, dsize=(256, 256), interpolation=cv2.INTER_LINEAR)\nplt.title(\"Retinal fundus image 256x256\")\nplt.imshow(resized)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As the ratio between the heights of the optical disc and the optical cup is symptomatic for glocaumas, the first step is to segment these two eye areas.","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport cv2\n\nimage = plt.imread(\"/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/train/gts/g0001.bmp\")\nresized = cv2.resize(image, dsize=(256, 256), interpolation=cv2.INTER_LINEAR)\nplt.title(\"Optical disc and optical cup segmented\")\nplt.imshow(resized)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, depending of this ratio, we have to classify if the eye has a glaucoma or not.","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport json    \n                                                                                               \ndf = pd.read_json(\"/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data/train/index.json\")\n\nl = [\"normal\", \"glaucomous\"]\nprint(\"In this case, the eye is a \" + str(l[df.iloc[5][0]]) + \" eye.\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## II) Data implementation","metadata":{}},{"cell_type":"markdown","source":"Useful metrics","metadata":{}},{"cell_type":"code","source":"EPS = 1e-7\n\ndef compute_dice_coef(input, target):\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    Input:\n        binary_segmentation: a boolean 2D numpy array representing a region of interest.\n    Output:\n        diameter: the vertical diameter of the structure, defined as the largest diameter between the upper and the lower interfaces\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    The vertical cup to disc ratio is defined as here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1722393/pdf/v082p01118.pdf\n    Input:\n        segmentation: binary 2D numpy array representing a segmentation, with 0: optic cup, 128: optic disc, 255: elsewhere.\n    Output:\n        cdr: vertical cup to disc ratio\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    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 get_sensitivity_at_given_specificity(sensitivity, specificity, specificity_reference=0.85):\n    '''\n    Get the sensitivity for a given specificity reference\n    Input:\n        sensitivity: sensitivity values\n        specificity: specificity values\n        [specificity_reference]: reference value for evaluation\n    Output:\n        sensitivity_value: sensitivity value at the specificity reference\n    '''\n\n    # interpolate a continue curve based on sensitivity and specificity\n    sensitivity_interp = interp1d(specificity, sensitivity)\n    # get the sensitivity value for the specificity reference\n    sensitivity_value = sensitivity_interp(specificity_reference)\n\n    return sensitivity_value\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\ndef fov_error(seg_fov, fov_coord):\n    mass_centers = [scipy.ndimage.measurements.center_of_mass(seg_fov[k,:,:]) for k in range(seg_fov.shape[0])]\n    mass_centers = np.array([[elt[1], elt[0]] for elt in mass_centers])\n    err = np.sqrt(np.sum((fov_coord-mass_centers)**2, axis=1)).mean()\n    return err","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Test to add horizontal CDR ###\n\ndef horizontal_diameter(binary_segmentation):\n    '''\n    Get the horizontal diameter from a binary segmentation.\n    The horizontal diameter is defined as the \"fattest\" area of the binary_segmentation parameter.\n    Input:\n        binary_segmentation: a boolean 2D numpy array representing a region of interest.\n    Output:\n        diameter: the horizontal diameter of the structure, defined as the largest diameter between the upper and the lower interfaces\n    '''\n\n    # get the sum of the pixels in the horizontal axis\n    horizontal_axis_diameter = np.sum(binary_segmentation, axis=2)\n\n # pick the maximum value\n    diameter = np.max(horizontal_axis_diameter, axis=1)\n\n    # return it\n    return diameter\n\n\n\ndef horizontal_cup_to_disc_ratio(od, oc):\n    '''\n    Compute the horizontal cup-to-disc ratio from a given labelling map.\n    The horizontal cup to disc ratio is defined as here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1722393/pdf/v082p01118.pdf\n    Input:\n        segmentation: binary 2D numpy array representing a segmentation, with 0: optic cup, 128: optic disc, 255: elsewhere.\n    Output:\n        cdr: horizontal cup to disc ratio\n    '''\n    # compute the cup diameter\n    cup_diameter = horizontal_diameter(oc)\n    # compute the disc diameter\n    disc_diameter = horizontal_diameter(od)\n\n    return cup_diameter / (disc_diameter + EPS)\n\ndef compute_hCDR_error(pred_od, pred_oc, gt_od, gt_oc):\n    pred_hCDR = horizontal_cup_to_disc_ratio(pred_od, pred_oc)\n    gt_hCDR = horizontal_cup_to_disc_ratio(gt_od, gt_oc)\n    hCDR_err = np.mean(np.abs(gt_hCDR - pred_hCDR))\n    return hCDR_err, pred_hCDR, gt_hCDR\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Create the dataset class.","metadata":{}},{"cell_type":"code","source":"## taken from the notebook\n\nclass 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":"Useful values for the network","metadata":{}},{"cell_type":"code","source":"root_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\nlr = 1e-4\nbatch_size = 40\nnum_workers = 40\ntotal_epoch = 150 # 100","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Create the dataset","metadata":{}},{"cell_type":"code","source":"train_set = RefugeDataset(root_dir, split='train')\nval_set = RefugeDataset(root_dir, split='val')\ntest_set = RefugeDataset(root_dir, split='test')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Check what is inside the training set ","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(1,3, figsize=(20,20))\n\n(img, lab, seg, fov, img_name)= train_set.__getitem__(0)\ntemp_img =img.numpy()\ntemp_img = temp_img.transpose((1, 2, 0)) # the channels in torch are inverted \nax[0].imshow(temp_img)\nax[0].set_title(str(img_name))\n\ntemp_seg =seg.numpy()\n\nod_seg = temp_seg[0,:,:]\noc_seg = temp_seg[1,:,:]\n\n\n\nax[1].imshow(od_seg)\nax[1].set_title(\"Optical disc\")\n\nax[2].imshow(oc_seg)\nax[2].set_title(\"Optical cup\")\n\nprint(\"Label = \" + str(lab.item()) )\n\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Load the dataset","metadata":{}},{"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                         )\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_data = iter(train_loader)\nnb_train_batches = len(train_loader)\nimgs, 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\ntrain_vCDRs, val_vCDRs = [], []\ntrain_vCDR_error, val_vCDR_error = 0., 0.\n\ntrain_hCDRs, val_hCDRs = [], []\ntrain_hCDR_error, val_hCDR_error = 0., 0.\n\n\ngt_od = seg_gts[:,0,:,:].type(torch.int8)\ngt_oc = seg_gts[:,1,:,:].type(torch.int8)\n\n\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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":"****","metadata":{}},{"cell_type":"markdown","source":"# Test new network","metadata":{}},{"cell_type":"markdown","source":"Device, model, loss and optimizer","metadata":{}},{"cell_type":"code","source":"!pip install segmentation_models_pytorch","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import segmentation_models_pytorch as smp\nCUDA_LAUNCH_BLOCKING=\"1\"\n\n# Device\ndevice = torch.device(\"cuda:0\")\n\n# Network\nmodel = smp.DeepLabV3Plus(\n    encoder_name=\"resnet34\",        # choose encoder, e.g. resnet\n    in_channels=3,                  # model input channels (1 for grayscale images, 3 for RGB, etc.)\n    classes=2,                      # model output channels (number of classes in your dataset)\n    activation = \"sigmoid\"\n).to(device)\n\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":"First, we can feed the logistic regression model only with the vCDR.","metadata":{}},{"cell_type":"code","source":"# Define parameters\nnb_train_batches = len(train_loader)\nnb_val_batches = len(val_loader)\nnb_iter = 0\nbest_val_auc = 0.\n\nepo = 0\n\nwhile epo < 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    ############\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).to(device)\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            \n            \n            \n        # Increase iterations\n        nb_iter += 1\n        \n        # Std out\n        print('Epoch {}, iter {}/{}, loss {:.6f}'.format(epo+1, k+1, nb_train_batches, loss.item()) + ' '*20, \n              end='\\r')\n        \n    # Train a logistic regression on vCDRs\n    train_vCDRs = np.array(train_vCDRs).reshape(-1,1)\n    train_classif_gts = np.array(train_classif_gts)\n    clf = LogisticRegression(random_state=0, solver='lbfgs').fit(train_vCDRs, train_classif_gts)\n    train_classif_preds = clf.predict_proba(train_vCDRs)[:,1]\n    train_auc = classif_eval(train_classif_preds, train_classif_gts)\n    ##############\n    # 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            # 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(epo+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    # 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    epo += 1","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then, we can also feed the logistic regression model with the vCDR and the hCDR.","metadata":{}},{"cell_type":"code","source":"# Define parameters\nnb_train_batches = len(train_loader)\nnb_val_batches = len(val_loader)\nnb_iter = 0\nbest_val_auc = 0.\n\nepo = 0\n\nwhile epo < 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    ############\n    # for horizontal CDR #\n    train_hCDRs, val_hCDRs = [], []\n    train_hCDR_error, val_hCDR_error = 0., 0.\n    ############\n    \n    \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).to(device)\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            \n            train_classif_gts += classif_gts.cpu().numpy().tolist()\n            \n            \n            # Compute and store hCDRs\n            hCDR_error, pred_hCDR, gt_hCDR = compute_hCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n            train_hCDRs += pred_hCDR.tolist()\n            train_hCDR_error += hCDR_error / nb_train_batches\n\n            \n            \n            \n            \n        # Increase iterations\n        nb_iter += 1\n        \n        # Std out\n        print('Epoch {}, iter {}/{}, loss {:.6f}'.format(epo+1, k+1, nb_train_batches, loss.item()) + ' '*20, \n              end='\\r')\n        \n    # Train a logistic regression on vCDRs AND hCDRs \n    train_vCDRs = np.array(train_vCDRs).reshape(-1,1)\n    train_hCDRs = np.array(train_hCDRs).reshape(-1,1)\n\n    train_vhCDRs = np.concatenate((train_vCDRs,train_hCDRs),axis=1)\n    \n    train_classif_gts = np.array(train_classif_gts)\n    clf = LogisticRegression(random_state=0, solver='lbfgs').fit(train_vhCDRs, train_classif_gts)\n    train_classif_preds = clf.predict_proba(train_vhCDRs)[:,1]\n    train_auc = classif_eval(train_classif_preds, train_classif_gts)\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            # 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 hCDRs\n            hCDR_error, pred_hCDR, gt_hCDR = compute_hCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n            val_hCDRs += pred_hCDR.tolist()\n            val_hCDR_error += hCDR_error / nb_val_batches\n\n    # Glaucoma predictions from vCDRs AND hCDRs\n    val_vCDRs = np.array(val_vCDRs).reshape(-1,1)\n    val_hCDRs = np.array(val_hCDRs).reshape(-1,1)\n    \n    val_vhCDRs = np.concatenate((val_vCDRs,val_hCDRs),axis=1)\n    \n    val_classif_gts = np.array(val_classif_gts)\n    val_classif_preds = clf.predict_proba(val_vhCDRs)[:,1]\n    val_auc = classif_eval(val_classif_preds, val_classif_gts)\n        \n    # Validation results\n    print('VALIDATION epoch {}'.format(epo+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('hCDR error: {:.4f} (train), {:.4f} (val)'.format(train_hCDR_error, val_hCDR_error))\n    print('Classification (AUC): {:.4f} (train), {:.4f} (val)'.format(train_auc, val_auc))\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    epo += 1","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Re-check result on validation set","metadata":{}},{"cell_type":"code","source":"model.eval()\nval_vCDRs = []\nval_hCDRs = []\nval_classif_gts = []\nval_loss = 0.\nval_dsc_od = 0.\nval_dsc_oc = 0.\nval_vCDR_error = 0.\nval_hCDR_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 hCDRs\n        hCDR_error, pred_hCDR, gt_hCDR = compute_hCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n        val_hCDRs += pred_hCDR.tolist()\n        val_hCDR_error += hCDR_error / nb_val_batches\n        \n        \n        \n\n# Glaucoma predictions from vCDRs AND hCDRs\nval_vCDRs = np.array(val_vCDRs).reshape(-1,1)\nval_hCDRs = np.array(val_hCDRs).reshape(-1,1)\n\nval_vhCDRs = np.concatenate((val_vCDRs,val_hCDRs),axis=1)\n\nval_classif_gts = np.array(val_classif_gts)\nval_classif_preds = clf.predict_proba(val_vhCDRs)[:,1]\nval_auc = classif_eval(val_classif_preds, val_classif_gts)\n        \n    \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('hCDR error: {:.4f} (val)'.format( val_hCDR_error))\nprint('Classification (AUC): {:.4f} (val)'.format(val_auc))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Prediction on test set","metadata":{}},{"cell_type":"code","source":"nb_test_batches = len(test_loader)\nmodel.eval()\ntest_vCDRs = []\nwith torch.no_grad():\n    test_data = iter(test_loader)\n    for k in range(nb_test_batches):\n        # Loads data\n        imgs = test_data.next()\n        imgs = imgs.to(device)\n\n        # Forward pass\n        logits = model(imgs)\n\n        # Std out\n        print('Test iter {}/{}'.format(k+1, nb_test_batches) + ' '*50, \n              end='\\r')\n            \n        # Compute segmentation\n        pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n        pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n            \n        # Compute and store vCDRs\n        pred_vCDR = vertical_cup_to_disc_ratio(pred_od.cpu().numpy(), pred_oc.cpu().numpy())\n        test_vCDRs += pred_vCDR.tolist()\n            \n\n    # Glaucoma predictions from vCDRs\n    test_vCDRs = np.array(test_vCDRs).reshape(-1,1)\n    test_classif_preds = clf.predict_proba(test_vCDRs)[:,1]\n\n       \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.val_hCDRs ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nb_test_batches = len(test_loader)\nmodel.eval()\ntest_vCDRs = []\ntest_hCDRs = []\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 hCDRs\n        pred_hCDR = horizontal_cup_to_disc_ratio(pred_od.cpu().numpy(), pred_oc.cpu().numpy())\n        test_hCDRs += pred_hCDR.tolist()\n            \n\n    # Glaucoma predictions from vCDRs\n    test_vCDRs = np.array(test_vCDRs).reshape(-1,1)\n    test_hCDRs = np.array(test_hCDRs).reshape(-1,1)\n    test_vhCDRs = np.concatenate((test_vCDRs, test_hCDRs), axis = 1)\n    \n    test_classif_preds = clf.predict_proba(test_vhCDRs)[:,1]\n\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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}