{"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 pandas as pd\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\n\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.svm import SVC as svc\nfrom sklearn.metrics import roc_auc_score, roc_curve\nfrom sklearn import feature_selection\nfrom sklearn import preprocessing\nfrom sklearn.ensemble import RandomForestClassifier as RFC\n\nfrom tensorflow.keras.preprocessing.image import ImageDataGenerator\nfrom tensorflow.keras.applications import VGG19, Xception, ResNet50, MobileNet\nfrom keras.models import Model\nfrom keras.layers import Input \nimport keras\nfrom keras import layers\nfrom keras import models\nfrom keras import optimizers\nfrom keras.models import Sequential\n\nfrom sklearn.preprocessing import StandardScaler\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We are next going to import the necessary tools to use ResNet18","metadata":{}},{"cell_type":"code","source":"pip install git+https://github.com/qubvel/classification_models.git","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from classification_models.keras import Classifiers","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Dataset class","metadata":{}},{"cell_type":"code","source":"class RefugeDataset(Dataset):\n\n    def __init__(self, root_dir, split='train', output_size=(256,256)):\n        # Define attributes\n        self.output_size = output_size\n        self.root_dir = root_dir\n        self.split = split\n        \n        # Load data index\n        with open(os.path.join(self.root_dir, self.split, 'index.json')) as f:\n            self.index = json.load(f)\n            \n        self.images = []\n        for k in range(len(self.index)):\n            print('Loading {} image {}/{}...'.format(split, k, len(self.index)), end='\\r')\n            img_name = os.path.join(self.root_dir, self.split, 'images', self.index[str(k)]['ImgName'])\n            img = np.array(Image.open(img_name).convert('RGB'))\n            img = 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":"EPS = 1e-7\n\ndef compute_dice_coef(input, target):\n    '''\n    Compute dice score metric.\n    '''\n    batch_size = input.shape[0]\n    return sum([dice_coef_sample(input[k,:,:], target[k,:,:]) for k in range(batch_size)])/batch_size\n\ndef dice_coef_sample(input, target):\n    iflat = input.contiguous().view(-1)\n    tflat = target.contiguous().view(-1)\n    intersection = (iflat * tflat).sum()\n    return (2. * intersection) / (iflat.sum() + tflat.sum())\n\n\ndef vertical_diameter(binary_segmentation):\n    '''\n    Get the vertical diameter from a binary segmentation.\n    The vertical diameter is defined as the \"fattest\" area of the binary_segmentation parameter.\n    '''\n\n    # get the sum of the pixels in the vertical axis\n    vertical_axis_diameter = np.sum(binary_segmentation, axis=1)\n\n    # pick the maximum value\n    diameter = np.max(vertical_axis_diameter, axis=1)\n\n    # return it\n    return diameter\n\n\n\ndef vertical_cup_to_disc_ratio(od, oc):\n    '''\n    Compute the vertical cup-to-disc ratio from a given labelling map.\n    '''\n    # compute the cup diameter\n    cup_diameter = vertical_diameter(oc)\n    # compute the disc diameter\n    disc_diameter = vertical_diameter(od)\n\n    return cup_diameter / (disc_diameter + EPS)\n\ndef compute_vCDR_error(pred_od, pred_oc, gt_od, gt_oc):\n    '''\n    Compute vCDR prediction error, along with predicted vCDR and ground truth vCDR.\n    '''\n    pred_vCDR = vertical_cup_to_disc_ratio(pred_od, pred_oc)\n    gt_vCDR = vertical_cup_to_disc_ratio(gt_od, gt_oc)\n    vCDR_err = np.mean(np.abs(gt_vCDR - pred_vCDR))\n    return vCDR_err, pred_vCDR, gt_vCDR\n\n\ndef classif_eval(classif_preds, classif_gts):\n    '''\n    Compute AUC classification score.\n    '''\n    auc = roc_auc_score(classif_gts, classif_preds)\n    return auc\n\n\ndef fov_error(pred_fov, gt_fov):\n    '''\n    Fovea localization error metric (mean root squared error).\n    '''\n    err = np.sqrt(np.sum((gt_fov-pred_fov)**2, axis=1)).mean()\n    return err","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Post-processing functions","metadata":{}},{"cell_type":"code","source":"def refine_seg(pred):\n    '''\n    Only retain the biggest connected component of a segmentation map.\n    '''\n    np_pred = pred.numpy()\n        \n    largest_ccs = []\n    for i in range(np_pred.shape[0]):\n        labeled, ncomponents = label(np_pred[i,:,:])\n        bincounts = np.bincount(labeled.flat)[1:]\n        if len(bincounts) == 0:\n            largest_cc = labeled == 0\n        else:\n            largest_cc = labeled == np.argmax(bincounts)+1\n        largest_cc = torch.tensor(largest_cc, dtype=torch.float32)\n        largest_ccs.append(largest_cc)\n    largest_ccs = torch.stack(largest_ccs)\n    \n    return largest_ccs","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Network","metadata":{}},{"cell_type":"code","source":"class UNet(nn.Module):\n    def __init__(self, n_channels=3, n_classes=2):\n        super(UNet, self).__init__()\n        self.n_channels = n_channels\n        self.n_classes = n_classes\n        self.epoch = 0\n\n        self.inc = DoubleConv(n_channels, 64)\n        self.down1 = Down(64, 128)\n        self.down2 = Down(128, 256)\n        self.down3 = Down(256, 512)\n        factor = 2 \n        self.down4 = Down(512, 1024 // factor)\n        self.up1 = Up(1024, 512 // factor)\n        self.up2 = Up(512, 256 // factor)\n        self.up3 = Up(256, 128 // factor)\n        self.up4 = Up(128, 64)\n        self.output_layer = OutConv(64, n_classes)\n\n    def forward(self, x):\n        x1 = self.inc(x)\n        x2 = self.down1(x1)\n        x3 = self.down2(x2)\n        x4 = self.down3(x3)\n        x5 = self.down4(x4)\n        out = self.up1(x5, x4)\n        out = self.up2(out, x3)\n        out = self.up3(out, x2)\n        out = self.up4(out, x1)\n        out = self.output_layer(out)\n        out = torch.sigmoid(out)\n        return out\n\n    \nclass DoubleConv(nn.Module):\n    \"\"\"(convolution => [BN] => ReLU) * 2\"\"\"\n\n    def __init__(self, in_channels, out_channels, mid_channels=None):\n        super().__init__()\n        if not mid_channels:\n            mid_channels = out_channels\n        self.double_conv = nn.Sequential(\n            nn.Conv2d(in_channels, mid_channels, kernel_size=3, padding=1),\n            nn.BatchNorm2d(mid_channels),\n            nn.ReLU(inplace=True),\n            nn.Conv2d(mid_channels, out_channels, kernel_size=3, padding=1),\n            nn.BatchNorm2d(out_channels),\n            nn.ReLU(inplace=True)\n        )\n\n    def forward(self, x):\n        return self.double_conv(x)\n\n\nclass Down(nn.Module):\n    \"\"\"Downscaling with maxpool then double conv\"\"\"\n\n    def __init__(self, in_channels, out_channels):\n        super().__init__()\n        self.maxpool_conv = nn.Sequential(\n            nn.MaxPool2d(2),\n            DoubleConv(in_channels, out_channels)\n        )\n\n    def forward(self, x):\n        return self.maxpool_conv(x)\n\n\nclass Up(nn.Module):\n    \"\"\"Upscaling then double conv\"\"\"\n\n    def __init__(self, in_channels, out_channels):\n        super().__init__()\n\n        # Use the normal convolutions to reduce the number of channels\n        self.up = nn.Upsample(scale_factor=2, mode='bilinear', align_corners=True)\n        self.conv = DoubleConv(in_channels, out_channels, in_channels // 2)\n\n\n    def forward(self, x1, x2):\n        x1 = self.up(x1)\n        # input is CHW\n        diffY = x2.size()[2] - x1.size()[2]\n        diffX = x2.size()[3] - x1.size()[3]\n\n        x1 = F.pad(x1, [diffX // 2, diffX - diffX // 2,\n                        diffY // 2, diffY - diffY // 2])\n        x = torch.cat([x2, x1], dim=1)\n        return self.conv(x)\n\n\nclass OutConv(nn.Module):\n    '''\n    Simple convolution.\n    '''\n    def __init__(self, in_channels, out_channels):\n        super(OutConv, self).__init__()\n        self.conv = nn.Conv2d(in_channels, out_channels, kernel_size=1)\n\n    def forward(self, x):\n        return self.conv(x)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Settings","metadata":{}},{"cell_type":"code","source":"root_dir = '/kaggle/input/eurecom-aml-2021-challenge-2/refuge_data/refuge_data'\nlr = 1e-4\nbatch_size = 8\nnum_workers = 8\ntotal_epoch = 100","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create datasets and data loaders\nAll image files are loaded in RAM in order to speed up the pipeline. Therefore, each dataset creation should take a few minutes.","metadata":{}},{"cell_type":"code","source":"# Datasets\ntrain_set = RefugeDataset(root_dir, \n                          split='train')\nval_set = RefugeDataset(root_dir, \n                        split='val')\ntest_set = RefugeDataset(root_dir, \n                         split='test')\n\n# Dataloaders\ntrain_loader = DataLoader(train_set, \n                          batch_size=batch_size, \n                          shuffle=False, \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":"nb_train_batches = len(train_loader)\nnb_val_batches = len(val_loader)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Device, model, loss and optimizer","metadata":{}},{"cell_type":"code","source":"# Device\ndevice = torch.device(\"cuda:0\")\n\n# Network\nmodel = UNet(n_channels=3, n_classes=2).to(device)\n\n# Loss\nseg_loss = torch.nn.BCELoss(reduction='mean')\n\n\n# Optimizer\noptimizer = optim.Adam(model.parameters(), lr=lr)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train for OC/OD segmentation","metadata":{}},{"cell_type":"markdown","source":"We decided to not shuffle the train data so that we can build other features for the moment outside the loop of the U-Net","metadata":{}},{"cell_type":"code","source":"# Define parameters\nnb_train_batches = len(train_loader)\nnb_val_batches = len(val_loader)\nnb_iter = 0\nbest_val_auc = 0.\n\nwhile model.epoch < total_epoch:\n    # Accumulators\n    train_vCDRs, val_vCDRs = [], [] \n    train_classif_gts, val_classif_gts = [], []\n    train_loss, val_loss = 0., 0.\n    train_dsc_od, val_dsc_od = 0., 0.\n    train_dsc_oc, val_dsc_oc = 0., 0.\n    train_vCDR_error, val_vCDR_error = 0., 0.\n    \n    ############\n    # TRAINING #\n    ############\n    model.train()\n    train_data = iter(train_loader)\n    for k in range(nb_train_batches):\n        # Loads data\n        imgs, classif_gts, seg_gts, fov_coords, names = train_data.next()\n        imgs, classif_gts, seg_gts = imgs.to(device), classif_gts.to(device), seg_gts.to(device)\n\n        # Forward pass\n        logits = model(imgs)\n        loss = seg_loss(logits, seg_gts)\n \n        # Backward pass\n        optimizer.zero_grad()\n        loss.backward()\n        optimizer.step()\n        train_loss += loss.item() / nb_train_batches\n        \n        with torch.no_grad():\n            # Compute segmentation metric\n            pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n            pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n            gt_od = seg_gts[:,0,:,:].type(torch.int8)\n            gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n            dsc_od = compute_dice_coef(pred_od, gt_od)\n            dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n            train_dsc_od += dsc_od.item()/nb_train_batches\n            train_dsc_oc += dsc_oc.item()/nb_train_batches\n\n\n            # Compute and store vCDRs\n            vCDR_error, pred_vCDR, gt_vCDR = compute_vCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n            train_vCDRs += pred_vCDR.tolist()\n            train_vCDR_error += vCDR_error / nb_train_batches\n            train_classif_gts += classif_gts.cpu().numpy().tolist()\n            \n        # Increase iterations\n        nb_iter += 1\n        \n        # Std out\n        print('Epoch {}, iter {}/{}, loss {:.6f}'.format(model.epoch+1, k+1, nb_train_batches, loss.item()) + ' '*20, \n              end='\\r')\n        \n    # Train a logistic regression on vCDRs\n    train_vCDRs = np.array(train_vCDRs).reshape(-1,1)\n    train_classif_gts = np.array(train_classif_gts)\n    clf = LogisticRegression(random_state=0, solver='lbfgs').fit(train_vCDRs, train_classif_gts)\n    train_classif_preds = clf.predict_proba(train_vCDRs)[:,1]\n    train_auc = classif_eval(train_classif_preds, train_classif_gts)\n    \n    ##############\n    # VALIDATION #\n    ##############\n    model.eval()\n    with torch.no_grad():\n        val_data = iter(val_loader)\n        for k in range(nb_val_batches):\n            # Loads data\n            imgs, classif_gts, seg_gts, fov_coords, names = val_data.next()\n            imgs, classif_gts, seg_gts = imgs.to(device), classif_gts.to(device), seg_gts.to(device)\n\n            # Forward pass\n            logits = model(imgs)\n            val_loss += seg_loss(logits, seg_gts).item() / nb_val_batches\n\n            # Std out\n            print('Validation iter {}/{}'.format(k+1, nb_val_batches) + ' '*50, \n                  end='\\r')\n            \n            # Compute segmentation metric\n            pred_od = refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n            pred_oc = refine_seg((logits[:,1,:,:]>=0.5).type(torch.int8).cpu()).to(device)\n            gt_od = seg_gts[:,0,:,:].type(torch.int8)\n            gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n            dsc_od = compute_dice_coef(pred_od, gt_od)\n            dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n            val_dsc_od += dsc_od.item()/nb_val_batches\n            val_dsc_oc += dsc_oc.item()/nb_val_batches\n            \n            # Compute and store vCDRs\n            vCDR_error, pred_vCDR, gt_vCDR = compute_vCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n            val_vCDRs += pred_vCDR.tolist()\n            val_vCDR_error += vCDR_error / nb_val_batches\n            val_classif_gts += classif_gts.cpu().numpy().tolist()\n            \n            \n\n    # Glaucoma predictions from vCDRs\n    val_vCDRs = np.array(val_vCDRs).reshape(-1,1)\n    val_classif_gts = np.array(val_classif_gts)\n    val_classif_preds = clf.predict_proba(val_vCDRs)[:,1]\n    val_auc = classif_eval(val_classif_preds, val_classif_gts)\n        \n    # Validation results\n    print('VALIDATION epoch {}'.format(model.epoch+1)+' '*50)\n    print('LOSSES: {:.4f} (train), {:.4f} (val)'.format(train_loss, val_loss))\n    print('OD segmentation (Dice Score): {:.4f} (train), {:.4f} (val)'.format(train_dsc_od, val_dsc_od))\n    print('OC segmentation (Dice Score): {:.4f} (train), {:.4f} (val)'.format(train_dsc_oc, val_dsc_oc))\n    print('vCDR error: {:.4f} (train), {:.4f} (val)'.format(train_vCDR_error, val_vCDR_error))\n    print('Classification (AUC): {:.4f} (train), {:.4f} (val)'.format(train_auc, val_auc))\n    \n    # Save model if best validation AUC is reached\n    if val_auc > best_val_auc:\n        torch.save(model.state_dict(), '/kaggle/working/best_AUC_weights.pth')\n        with open('/kaggle/working/best_AUC_classifier.pkl', 'wb') as clf_file:\n            pickle.dump(clf, clf_file)\n        best_val_auc = val_auc\n        print('Best validation AUC reached. Saved model weights and classifier.')\n    print('_'*50)\n        \n    # End of epoch\n    model.epoch += 1\n        ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Getting weights of Model\n","metadata":{}},{"cell_type":"markdown","source":"We decided to save our weights fromerly computed so that we don't have to execute the U-net anymore.","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('../input/best-auc-weightspth/best_AUC_weights.pth'))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Adding Features","metadata":{}},{"cell_type":"markdown","source":"The first method that came to our mind to improve the accuracy was to create some features for each image and take them into account to train and use our classifier. We therefore created 26 features (some are very similar) and then tried to select the most pertinent one.","metadata":{}},{"cell_type":"markdown","source":"**Feature Engineering**","metadata":{}},{"cell_type":"code","source":"from skimage.measure import regionprops\nfrom skimage.measure import find_contours\nimport skimage.morphology as morpho\n\n\ndef average_contrast(ori, region):\n    \n    \"\"\"This function gives the features of the image which corresponds to the average contrast of the region\n    input : ori (original image), region (segmentation)\n    ouput : contrast_R, contrast_G, contrast_B the set of average contrast on every channel\"\"\"\n    \n    ori = np.array(ori)    #computing the gradient norm of the image\n    se=morpho.selem.disk(1)\n        \n    contrast = []\n    \n    \n    for i in range(len(ori)) :\n        im = ori[i]\n        grad = morpho.dilation(im[:,:], se)- morpho.erosion(im[:,:], se)\n        \n        seg = region[i]\n        contrast.append(np.mean(grad[seg>0],axis=0))\n       \n        \n    return contrast\ndef norm_grad(im):\n    sobelx = cv2.Sobel(im, cv2.CV_64F, 1, 0, ksize=5)\n    sobely = cv2.Sobel(im, cv2.CV_64F, 0, 1, ksize=5)\n    return(sobelx**2 + sobely**2)**0.5\n\ndef featuresComputation(im, segCup, segDisc):\n    sumN = np.sum(segDisc)\n    features = []\n    #  ratio and areas\n    features.append(np.sum(segCup))\n    features.append(np.sum(segDisc))\n    features.append(np.sum(segDisc)/np.sum(segCup))\n    \n    # relative position\n    cCup = regionprops(segCup)[0].centroid\n    if sumN!=0:\n        cDisc = regionprops(segDisc)[0].centroid\n    else:\n        cDisc = (0,0)\n    features.append(cCup[0]-cDisc[0])  # decalage en x\n    features.append(cCup[1]-cDisc[1])  # decalage en y\n    \n    # eccentricity\n    features.append(regionprops(segCup)[0].eccentricity)\n    if sumN !=0:\n        features.append(regionprops(segDisc)[0].eccentricity)\n    else :\n        features.append(0)\n    \n    # major_axis_length\n    features.append(regionprops(segCup)[0].major_axis_length)\n    if sumN !=0:\n        features.append(regionprops(segDisc)[0].major_axis_length)\n    else :\n        features.append(0)\n    \n    # minor_axis_length\n    features.append(regionprops(segCup)[0].minor_axis_length)\n    if sumN !=0:\n        features.append(regionprops(segDisc)[0].minor_axis_length)\n    else :\n        features.append(0)\n    \n    # perimeter\n    features.append(regionprops(segCup)[0].perimeter)\n    if sumN !=0:\n        features.append(regionprops(segDisc)[0].perimeter)\n    else :\n        features.append(0)\n    \n    # orientation\n    features.append(regionprops(segCup)[0].orientation)\n    if sumN !=0:\n        features.append(regionprops(segDisc)[0].orientation)\n    else :\n        features.append(0)\n    \n    if sumN !=0:\n        contours_Cup = find_contours(segCup,0.5)\n        contours_Disc = find_contours(segDisc,0.5)\n        tension1 = np.gradient(contours_Cup[0],axis=0)\n        tension2 = np.gradient(contours_Disc[0],axis=0)\n        features.append(np.abs(tension1).mean())\n        features.append(np.abs(tension2).mean())\n    \n        courbure1 = np.gradient(tension1,axis=0)\n        courbure2 = np.gradient(tension2,axis=0)\n        features.append(np.abs(courbure1).mean())\n        features.append(np.abs(courbure2).mean())\n    else : \n        features.append(0)\n        features.append(0)\n        features.append(0)\n        features.append(0)\n    \n    \n    \n    # average intensity\n    av_Cup = (im*segCup).mean()\n    av_Disc = (im*segDisc).mean()\n    features.append(av_Cup)\n    features.append(av_Disc)\n        \n    # standard deviation \n    std_Cup = (im*segCup).std()\n    std_Disc = (im*segDisc).std()\n    features.append(std_Cup)\n    features.append(std_Disc)\n        \n    # gradient\n    grad_Cup = norm_grad(im*segCup).mean()\n    grad_Disc = norm_grad(im*segDisc).mean()\n    features.append(grad_Cup)\n    features.append(grad_Disc)        \n        \n    return features\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now that we have functions that can compute the wanted features we apply it to the dataset in order to train the classifier with them.\n","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n#Computing the right VCDRs (with the good weight)\ntrain_vCDRs = []\ntrain_classif_gts = []\ntrain_loss = 0.\ntrain_dsc_od = 0.\ntrain_dsc_oc = 0.\ntrain_vCDR_error = 0.\n\nContrast = []\n\nFeaturesList = []\n\nwith torch.no_grad():\n    train_data = iter(train_loader)\n    for k in range(nb_train_batches):\n    \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            logits = model(imgs)\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            train_dsc_od += dsc_od.item()/nb_val_batches\n            train_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            train_vCDRs += pred_vCDR.tolist()\n            train_vCDR_error += vCDR_error / nb_val_batches\n            train_classif_gts += classif_gts.cpu().numpy().tolist()\n            \n            \n            # Computing the contrast Feature\n            \n            imgs = imgs.cpu()\n            logits = logits.cpu()\n\n            ctrs = average_contrast(imgs[:,0,:,:],logits[:,0,:,:])    \n            Contrast += ctrs\n            # Computing the other features\n            features = []\n            for i in range (8):\n                features.append(featuresComputation(imgs[i,0,:,:].numpy(), (logits[i,1,:,:].numpy()>0.1).astype(int), (logits[i,0,:,:].numpy()>0.1).astype(int)))\n            FeaturesList += features","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then we combine every features in the same array and make a DataFrame for readability and analysis. We need to select the most pertinent feature. (too much feature is no good for our accuracy)\n","metadata":{}},{"cell_type":"code","source":"FeaturesList = np.array(FeaturesList).reshape(-1,25)\ntrain_vCDRs = np.array(train_vCDRs).reshape(-1,1)\nContrast = np.array(Contrast).reshape(-1,1)\nfeatures2 = np.concatenate((train_vCDRs,Contrast),axis=1)\nfeaturesAll = np.concatenate((train_vCDRs,Contrast,FeaturesList),axis=1)\n\nFeatureNames = ['VCDR','Contrast','SumCup','SumDisc','RationSum','X_Gap','y_Gap','Cup_Ecc','Disc_Ecc','Major_Axis_Cup','Major_axis_Disc','Minor_axis_Cup','Minor_axis_disk','Perimeter_cup','Perimeter_Disk','Orientation_Cup','Orientation_Disk','Tenstion_cup','Tension_disk','Courbure_cup','Courbure_disk','av_Cup','av_Disk','std_cup','std_Disk','gradient_cup','gradient_Disk','Label']\n\nTable = pd.DataFrame(np.concatenate((featuresAll,np.array(train_classif_gts).reshape(-1,1)),axis=1),columns=FeatureNames)\n\nprint(Table.head())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Feature Selection**\n\nIn this case we have a numerical input, categorical output. This is a classification predictive modeling problem with numerical input variables. Again, the most common techniques are correlation based, although in this case, they must take the categorical target into account.\n\nBecause of this we decided to use ANOVA correlation coefficient (linear) and the mutual information.\n\nLet's use Anova -f from sklearn:","metadata":{}},{"cell_type":"code","source":"# configure to select all features\nfs = feature_selection.SelectKBest(k='all')\n\nX_train = Table[['VCDR','Contrast','SumCup','SumDisc','RationSum','X_Gap','y_Gap','Cup_Ecc','Disc_Ecc','Major_Axis_Cup','Major_axis_Disc','Minor_axis_Cup','Minor_axis_disk','Perimeter_cup','Perimeter_Disk','Orientation_Cup','Orientation_Disk','Tenstion_cup','Tension_disk','Courbure_cup','Courbure_disk','av_Cup','av_Disk','std_cup','std_Disk','gradient_cup','gradient_Disk']]\ny_train = np.array(Table['Label']).reshape(-1,1)\n# learn relationship from training data\nfs.fit(X_train,y_train)\n\n# plot the scores\nplt.bar(['VCDR','Contrast','SumCup','SumDisc','RationSum','X_Gap','y_Gap','Cup_Ecc','Disc_Ecc','Major_Axis_Cup','Major_axis_Disc','Minor_axis_Cup','Minor_axis_disk','Perimeter_cup','Perimeter_Disk','Orientation_Cup','Orientation_Disk','Tenstion_cup','Tension_disk','Courbure_cup','Courbure_disk','av_Cup','av_Disk','std_cup','std_Disk','gradient_cup','gradient_Disk'], fs.scores_)\nplt.xticks(rotation=90,fontsize=8, fontname='monospace')\nplt.yticks(fontsize=20, fontname='monospace')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This clearly shows that feature VCDRs / SumCup / Major_Axis_Cup / Minor_axis_Cup / Perimeter_Cup /av_cup / std_cup / gradient_cup might be the most relevant (according to test) and that perhaps six of the eight input features are the more relevant.","metadata":{}},{"cell_type":"markdown","source":"Let's compute **the mutual information** to confirm our result from adova. The scikit-learn machine learning library provides an implementation of mutual information for feature selection with numeric input and categorical output variables via the mutual_info_classif() function.","metadata":{}},{"cell_type":"code","source":" # configure to select all features\nfs = feature_selection.SelectKBest(score_func=feature_selection.mutual_info_classif, k='all')\n# learn relationship from training data\nfs.fit(X_train, y_train)\n\n# plot the scores\nplt.bar(['VCDR','Contrast','SumCup','SumDisc','RationSum','X_Gap','y_Gap','Cup_Ecc','Disc_Ecc','Major_Axis_Cup','Major_axis_Disc','Minor_axis_Cup','Minor_axis_disk','Perimeter_cup','Perimeter_Disk','Orientation_Cup','Orientation_Disk','Tenstion_cup','Tension_disk','Courbure_cup','Courbure_disk','av_Cup','av_Disk','std_cup','std_Disk','gradient_cup','gradient_Disk'], fs.scores_)\nplt.xticks(rotation=90,fontsize=8, fontname='monospace')\nplt.yticks(fontsize=20, fontname='monospace')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This confirms the selection of feature that ADOVA gave. ","metadata":{}},{"cell_type":"markdown","source":"We now need to select a combinaison of those that gives us the best results. We didn't find any other way than the very simple empiric way of testing with several combinaison. The results on the validation set are at the end of this notebook.","metadata":{}},{"cell_type":"code","source":"# We test different combinaison of features (from the ones with some correlation)\n\n#a simple pre-precessing step\nX_scaler = StandardScaler().fit(featuresAll)\nfeaturesAll = X_scaler.transform(featuresAll)\n\nselectedFeatures1 = featuresAll[:,(0,2,4,9,11,13)]\nselectedFeatures2 = featuresAll[:,(0,2,4,9,11,13,21,23,25)]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Classifier\n\nWe decided to separate the feature extraction and the classification so that we have a clearer way of introducing new \nfeatures to the process. We are training several classifier to compare the results.\n\n","metadata":{}},{"cell_type":"code","source":"# Train a logistic regression on vCDRs with the good weight\ntrain_classif_gts = np.array(train_classif_gts).reshape(-1,1)\ntrain_vCDRs =  np.array(train_vCDRs).reshape(-1,1)\nclf = LogisticRegression(random_state=0, solver='lbfgs').fit(train_vCDRs, train_classif_gts)\ntrain_classif_preds = clf.predict_proba(train_vCDRs)[:,1]\n\n        ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Train a logistic regression with the selected features\n\nclf2 = LogisticRegression(random_state=0, solver='lbfgs').fit(selectedFeatures1, train_classif_gts)\n\nclf3 = LogisticRegression(random_state=0, solver='lbfgs').fit(selectedFeatures2, train_classif_gts)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The result aren't that encouraging when using logistic regression on a selection of features. Let's try to use random forest, this model is very transparent an would give us information on where to improve our features.","metadata":{}},{"cell_type":"code","source":"rfc1 = RFC(random_state=0)\nrfc1.fit(selectedFeatures1,train_classif_gts)\n\nrfc2 = RFC(random_state=0)\nrfc2.fit(selectedFeatures2,train_classif_gts)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Let's try a CNN directly on the picture \n\nGlaucoma displays its main clinical symptoms in the optic disc region. Based on this mechanism, with optic disc segmentation predictions, we first crop the regions of interest and resize the patches to size 60*60 to exclude irrelevant\nbackground contexts. \n\nIn order to improve our results we use Histogram equalization, CLAHE, for contrast enhancement\nand normalization so that even though the images are not taken in the exact same condition it should be processed correclty as a whole. \n\nThe image are therefore enhanced and way smaller so they can be processed by a CNN directly.\n","metadata":{}},{"cell_type":"code","source":"pip install progressbar","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport cv2\nfrom progressbar import ProgressBar\npbar = ProgressBar()\n\n# create a CLAHE object .\nclahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8,8))\n\nimagesTrain = []\nsick_vector_train = []\nthreshold = 0.5\n\nwith torch.no_grad():\n    train_data = iter(train_loader)\n    for k in pbar(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        logits = model(imgs)\n\n        imgs = imgs.cpu()\n        logits = logits.cpu()\n        masks = np.array(refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8)))\n        for p in range (len(imgs[:,0,:,:])):\n            x = []\n            y = []\n            mask = masks[p]\n            sick_vector_train.append(classif_gts[p].item())\n            for i in range (mask.shape[0]):\n                for j in range (mask.shape[1]):\n                    if mask[i][j]>threshold:\n                        x.append(i)\n                        y.append(j)\n                       \n            imga = imgs[p,:,:,:].numpy()[:,max(round(np.mean(x))-30,0):min(round(np.mean(x))+30,256),max(round(np.mean(y))-30,0):min(round(np.mean(y))+30,256)]\n            \n            \n             #Let's resize those that weren't well done \n            if(imga.shape) != (3,60, 60):\n                image = np.zeros((3,60,60))\n                for idx in range(len(imga)):\n                    img = imga[idx, :, :]\n                    img_sm = cv2.resize(img, (60, 60), interpolation=cv2.INTER_CUBIC)\n                    image[idx, :, :] = img_sm\n                imga = image\n            \n            #that allows to have the dimension in the right order !\n            imae = np.zeros((60,60,3))\n            for i in range(3):\n                imae[:,:,i] = clahe.apply((imga[i,:,:]*255).astype(np.uint8))\n                \n            imagesTrain.append(np.array(imae))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Let's do the same with the validation data \npbar = ProgressBar()\n\nimagesVal = []\nsick_vector_Val = []\nwith torch.no_grad():\n    val_data = iter(val_loader)\n    for k in pbar(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        logits = model(imgs)\n\n        imgs = imgs.cpu()\n        logits = logits.cpu()\n        masks = np.array(refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8)))\n        for p in range (len(imgs[:,0,:,:])):\n            x = []\n            y = []\n            mask = masks[p]\n            sick_vector_Val.append(classif_gts[p].item())\n            for i in range (mask.shape[0]):\n                for j in range (mask.shape[1]):\n                    if mask[i][j]>threshold:\n                        x.append(i)\n                        y.append(j)\n                       \n          \n            imga = imgs[p,:,:,:].numpy()[:,max(round(np.mean(x))-30,0):min(round(np.mean(x))+30,256),max(round(np.mean(y))-30,0):min(round(np.mean(y))+30,256)]\n            \n            \n            \n            \n            #Let's resize those that weren't well done \n            if(imga.shape) != (3,60, 60):\n                image = np.zeros((3,60,60))\n                for idx in range(len(imga)):\n                    img = imga[idx, :, :]\n                    img_sm = cv2.resize(img, (60, 60), interpolation=cv2.INTER_CUBIC)\n                    image[idx, :, :] = img_sm\n                imga = image\n            \n\n            #that allows to have the dimension in the right order !\n            imae = np.zeros((60,60,3))\n            for i in range(3):\n                #We apply the histogram equalization CLAHE\n                imae[:,:,i] = clahe.apply((imga[i,:,:]*255).astype(np.uint8))\n\n                    \n            imagesVal.append(np.array(imae))\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Note : We wanted to use the high resolution image and then apply to it the segmentation mask but getting the real image means having another dimensions and therefore the mask won't work on the high resolution image.","metadata":{}},{"cell_type":"markdown","source":"Now that we have the right images (focused on the usefull zone) we are going to pre-process it so that it can bu use in several neural-networks.","metadata":{}},{"cell_type":"code","source":"sick_vector_Train = np.array(sick_vector_train).reshape(-1,1)\nimagesTrain=np.array(imagesTrain)\nX_train_cnn = imagesTrain.reshape(imagesTrain.shape[0], 60, 60, 3)\n\nsick_vector_Val = np.array(sick_vector_Val).reshape(-1,1)\nimagesVal=np.array(imagesVal)\nX_val_cnn = imagesVal.reshape(imagesVal.shape[0], 60, 60, 3)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To improve the performance of the models we decided to use a part of the validation set to train the models.\nSo now we have the training and validation set to train and we took out a 100 images out of the validation set \nto validate.","metadata":{}},{"cell_type":"code","source":"#Let's create the new training set\nnew_sick_vector_Train = np.concatenate((sick_vector_Train,sick_vector_Val[-250:]),axis=0)\nnew_Xtrain_cnn = np.concatenate((X_train_cnn,X_val_cnn[-250:]),axis=0)\n\n\n\n#Let's create the new validation set \n\nnew_sick_vector_Val = sick_vector_Val[:150]\nnew_Xval_cnn = X_val_cnn[:150]\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We decided to use Data Augmentation : We are creating some images from the ones we have in order to improve the accuracy of our models. Those tools also allows to have the right format input objects for our Keras nets.\n","metadata":{}},{"cell_type":"code","source":"# All images will be rescaled by 1./255. , the rest should prevent from overfitting, since we are using a small data set it can\n#be a problem\ntrain_datagen = ImageDataGenerator( rescale = 1.0/255.,rotation_range = 40,width_shift_range=0.2,height_shift_range=0.2,\n                                  shear_range=0.2,zoom_range=0.2)\n# --------------------\n# Flow training images in batches of batch_size using train_datagen generator\n# --------------------\ntrain_generator = train_datagen.flow(new_Xtrain_cnn,new_sick_vector_Train,batch_size=batch_size)\n\n# All images will be rescaled by 1./255.\nval_datagen = ImageDataGenerator( rescale = 1.0/255.)\n# --------------------\n# Flow training images in batches of batch_size using train_datagen generator\n# --------------------\nval_generator =  val_datagen.flow(new_Xval_cnn,new_sick_vector_Val,batch_size=batch_size)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Now let's try several neural net**\n\nOur goal is to try a few nets and compare the results to select the most suited one.\n\n**The first one is a simple homemade CNN**","metadata":{}},{"cell_type":"code","source":"model_cnn2 = Sequential()\nmodel_cnn2.add(layers.Conv2D(input_shape=(60,60,3),filters=64,kernel_size=(3,3),padding=\"same\", activation=\"relu\"))\nmodel_cnn2.add(layers.BatchNormalization())\nmodel_cnn2.add(layers.MaxPooling2D(pool_size=(2,2),strides=(2,2)))\n\nmodel_cnn2.add(layers.Conv2D(filters=128,kernel_size=(3,3),padding=\"same\", activation=\"relu\"))\nmodel_cnn2.add(layers.BatchNormalization())\nmodel_cnn2.add(layers.MaxPooling2D(pool_size=(2,2),strides=(2,2)))\n\nmodel_cnn2.add(layers.Conv2D(filters=256,kernel_size=(3,3),padding=\"same\", activation=\"relu\"))\nmodel_cnn2.add(layers.BatchNormalization())\nmodel_cnn2.add(layers.MaxPooling2D(pool_size=(2,2),strides=(2,2)))\n\nmodel_cnn2.add(layers.Conv2D(filters=512,kernel_size=(3,3),padding=\"same\", activation=\"relu\"))\nmodel_cnn2.add(layers.BatchNormalization())\nmodel_cnn2.add(layers.MaxPooling2D(pool_size=(2,2),strides=(2,2)))\n\nmodel_cnn2.add(layers.Flatten())\nmodel_cnn2.add(layers.Dense(units=1024,activation=\"relu\"))\nmodel_cnn2.add(layers.Dense(units=1024,activation=\"relu\"))\nmodel_cnn2.add(layers.Dense(units=1, activation=\"sigmoid\"))\n\n# Compiling the CNN\nmodel_cnn2.compile(loss = 'binary_crossentropy',\n              optimizer=optimizers.Adam(lr=1e-4),\n              metrics = ['acc'])\n\n\n#Saving the best model\n\ncheckpoint_path = \"./modelCnn.h5\"\n\n\ncp_callback = keras.callbacks.ModelCheckpoint(filepath=checkpoint_path,\n                                                 save_weights_only=True,\n                                                 verbose=1,\n                                                 save_best_only = True)\nmodel_cnn2.summary()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history = model_cnn2.fit(train_generator, validation_data = val_generator, epochs=30, verbose=1,callbacks = [cp_callback])\n\n,\nplt.plot(history.history['acc'])\nplt.plot(history.history['val_acc'])\nplt.title('model accuracy')\nplt.ylabel('accuracy')\nplt.xlabel('epoch')\nplt.legend(['train', 'test'], loc='upper left')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_cnn2.load_weights('modelCnn.h5')\nmodel_cnn2.evaluate(new_Xtest_cnn, new_sick_vector_Test)\ntest_classif_preds = model_cnn2.predict(new_X_val_cnn)\nval_aucCnn = classif_eval(test_classif_preds, new_sick_vector_Test)\nprint('The AUC score value is :',val_aucCnn)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Comment**\n\nThis model gives us pretty a very good AUC. What we notice is that we had to find a balance in the complexity of the model. Too much layers created a very important overfitting phenomenon. (That's actually what's happening for the resNet).\n\nWe also notice that the model with the best val_accuracy isn't the one that generalize the best, maybe that's because we use a large part of the validation set to train and therefor the result are better for some weights. (to explore)","metadata":{}},{"cell_type":"markdown","source":"**We also tried a pre-trained one, ResNet18 :** \n\nResNet18 is a 18 layers deep model, so the process time is short, and according to what we read is very suited for image classification like the one we do here.","metadata":{}},{"cell_type":"code","source":"ResNet18, preprocess_input = Classifiers.get('resnet18')\nmodelResNet18 = ResNet18((60, 60, 3), weights='imagenet')\n\nfor layer in modelResNet18.layers[:]:\n    layer.trainable = False\n    \nmodelRes = Sequential()\nmodelRes.add(modelResNet18)\nmodelRes.add(layers.BatchNormalization())\nmodelRes.add(layers.Flatten())\nmodelRes.add(layers.Dense(units=1024,activation=\"relu\"))\nmodelRes.add(layers.Dense(units=1024,activation=\"relu\"))\nmodelRes.add(layers.Dense(units=1, activation=\"sigmoid\"))\n\nmodelRes.compile(loss='binary_crossentropy', optimizer=optimizers.Adam(lr=1e-6),metrics=[\"accuracy\"])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"historyRes = modelRes.fit(train_generator, validation_data = val_generator, epochs=15, verbose=1,batch_size=8)\n\nplt.plot(historyRes.history['accuracy'])\nplt.plot(historyRes.history['val_accuracy'])\nplt.title('model accuracy')\nplt.ylabel('accuracy')\nplt.xlabel('epoch')\nplt.legend(['train', 'test'], loc='upper left')\nplt.show()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Performance Comments**\n\n","metadata":{}},{"cell_type":"markdown","source":"What we think is happening here is an important overfitting, that seems to explain why the val_accuracy doesn't change.\nTo prevent that we tried adding batchNormalization et use the imageGenerator tool but we haven't been able to avoid this problem. We would need more data to train the resNet. \n","metadata":{}},{"cell_type":"markdown","source":"# Check performance is maintained on validation","metadata":{}},{"cell_type":"code","source":"model.eval()\nval_vCDRs = []\nval_classif_gts = []\nval_loss = 0.\nval_dsc_od = 0.\nval_dsc_oc = 0.\nval_vCDR_error = 0.\nfeature2 = []\nFeaturesList = []\nimagesCroppedVal = []\n\n# create a CLAHE object .\nclahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8,8))\n\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        \n        gt_od = seg_gts[:,0,:,:].type(torch.int8)\n        gt_oc = seg_gts[:,1,:,:].type(torch.int8)\n        dsc_od = compute_dice_coef(pred_od, gt_od)\n        dsc_oc = compute_dice_coef(pred_oc, gt_oc)\n        val_dsc_od += dsc_od.item()/nb_val_batches\n        val_dsc_oc += dsc_oc.item()/nb_val_batches\n\n        # Compute and store vCDRs\n        vCDR_error, pred_vCDR, gt_vCDR = compute_vCDR_error(pred_od.cpu().numpy(), pred_oc.cpu().numpy(), gt_od.cpu().numpy(), gt_oc.cpu().numpy())\n        val_vCDRs += pred_vCDR.tolist()\n        val_vCDR_error += vCDR_error / nb_val_batches\n        val_classif_gts += classif_gts.cpu().numpy().tolist()\n        \n        \n        logits = model(imgs)\n        imgs = imgs.cpu()\n        logits = logits.cpu()\n        \n        contrast = average_contrast(imgs[:,0,:,:],logits[:,0,:,:])    \n        feature2 += contrast\n        \n        features = []\n    \n        for i in range (8):\n            features.append(featuresComputation(imgs[i,0,:,:].numpy(), (logits[i,0,:,:].numpy()>0.1).astype(int), (logits[i,1,:,:].numpy()>0.1).astype(int)))\n        FeaturesList += features\n        \n        \n        #Now the CNN method\n        \n        masks = np.array(refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8)))\n        for p in range (len(imgs[:,0,:,:])):\n            x = []\n            y = []\n            mask = masks[p]\n            \n            for i in range (mask.shape[0]):\n                for j in range (mask.shape[1]):\n                    if mask[i][j]>threshold:\n                        x.append(i)\n                        y.append(j)\n                       \n            imga = imgs[p,:,:,:].numpy()[:,max(round(np.mean(x))-30,0):min(round(np.mean(x))+30,256),max(round(np.mean(y))-30,0):min(round(np.mean(y))+30,256)]\n            \n            \n             #Let's resize those that weren't well done \n            if(imga.shape) != (3,60, 60):\n                image = np.zeros((3,60,60))\n                for idx in range(len(imga)):\n                    img = imga[idx, :, :]\n                    img_sm = cv2.resize(img, (60, 60), interpolation=cv2.INTER_CUBIC)\n                    image[idx, :, :] = img_sm\n                imga = image\n            \n            #that allows to have the dimension in the right order !\n            imae = np.zeros((60,60,3))\n            for i in range(3):\n                imae[:,:,i] = clahe.apply((imga[i,:,:]*255).astype(np.uint8))\n                \n            imagesCroppedVal.append(np.array(imae))\n            \n            \nimagesCroppedVal=np.array(imagesCroppedVal)\nX_cnn = imagesCroppedVal.reshape(imagesCroppedVal.shape[0], 60, 60, 3)\n\nmodel_cnn2.load_weights('modelCnn.h5')\npreds_cnn2 = model_cnn2.predict(X_cnn)\nval_aucCnn  = classif_eval(preds_cnn2.reshape(1,-1)[0], val_classif_gts)\n        \npreds_Res = modelRes.predict(X_cnn)\nval_aucRes  = classif_eval(preds_Res.reshape(1,-1)[0], val_classif_gts)\n        \n        \n        \nfeature2 = np.array(feature2).reshape(-1,1)\nval_vCDRs = np.array(val_vCDRs).reshape(-1,1)\nFeaturesList = np.array(FeaturesList).reshape(-1,25)\n\n#Then we concatenate them to have one table of features\nfeatures2 = np.concatenate((val_vCDRs,feature2),axis=1)\nfeaturesAll = np.concatenate((features2,FeaturesList),axis=1)\nfeaturesAll = X_scaler.transform(featuresAll)\n\n#And we select the ones usefull\nselectedFeatures1 = featuresAll[:,(0,2,4,9,11,13)]\nselectedFeatures2 = featuresAll[:,(0,2,4,9,11,13,21,23,25)]\n        \n# Glaucoma predictions from vCDRs\n\nval_classif_gts = np.array(val_classif_gts)\n\nval_classif_preds = clf.predict_proba(val_vCDRs)[:,1]\nval_auc = classif_eval(val_classif_preds, val_classif_gts)\n\nval_classif_preds2 = clf2.predict_proba(selectedFeatures1)[:,1]\nval_auc2 = classif_eval(val_classif_preds2, val_classif_gts)\n\nval_classif_preds3 = clf3.predict_proba(selectedFeatures2)[:,1]\nval_auc3 = classif_eval(val_classif_preds3, val_classif_gts)\n\n\n# Let's try the RFC\n\nval_classif_preds_RFC_1 = rfc1.predict_proba(selectedFeatures1)[:,1]\nval_auc_rfc_1 = classif_eval(val_classif_preds_RFC_1, val_classif_gts)\n\nval_classif_preds_RFC_2 = rfc2.predict_proba(selectedFeatures2)[:,1]\nval_auc_rfc_2  = classif_eval(val_classif_preds_RFC_2, val_classif_gts)\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('Classification (AUC) with vCDRs: {:.4f} (val)'.format(val_auc))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Let's plot the results of our different methods\nResults = [val_auc,val_auc2,val_auc3,val_auc_rfc_1,val_auc_rfc_2,val_aucCnn,val_aucRes]\nlabels = [\"LR_vCDRS\", \"LR_select1\", \"LR_select2\", \"RF_select1\", \"RF_select2\", \"CNN\", \"ResNet\"]\n\nwidth = 0.35  # the width of the bars\nrects1 = plt.bar(labels, np.round(Results,2), width)\nplt.xticks(rotation=90,fontsize=8, fontname='monospace')\nplt.bar_label(rects1, padding=3)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Results Comment : Methode 1 - using features :**\n\nWhat we see here is that no combination of features beats the use of vCDRs alone. We have tried a lot of combinaison from 25 features but every time the AUC goes down.\nWe tried using the random forst classifier but it also fails to push the accuracy up. \nWhat we concluded of this results is that there is few data to train correctly those models (400 inputs) so the models don't generalize correctly.\n\nWe have another method using CNN direclty on the zone of interest and we believe that this second method have much more potential than the classic use of features.\n\n**Results Comment : Methode 2 - using CNN on the zone of interest :**\n","metadata":{}},{"cell_type":"markdown","source":"# Predictions on test set","metadata":{}},{"cell_type":"code","source":"nb_test_batches = len(test_loader)\nmodel.eval()\ntest_vCDRs = []\nfeature = []\nimagesCroppedTest = []\n\n# create a CLAHE object .\nclahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8,8))\n\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\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        imgs = imgs.cpu()\n        logits = logits.cpu()\n        \n        masks = np.array(refine_seg((logits[:,0,:,:]>=0.5).type(torch.int8)))\n        for p in range (len(imgs[:,0,:,:])):\n            x = []\n            y = []\n            mask = masks[p]\n            \n            for i in range (mask.shape[0]):\n                for j in range (mask.shape[1]):\n                    if mask[i][j]>threshold:\n                        x.append(i)\n                        y.append(j)\n                       \n            imga = imgs[p,:,:,:].numpy()[:,max(round(np.mean(x))-30,0):min(round(np.mean(x))+30,256),max(round(np.mean(y))-30,0):min(round(np.mean(y))+30,256)]\n            \n            \n             #Let's resize those that weren't well done \n            if(imga.shape) != (3,60, 60):\n                image = np.zeros((3,60,60))\n                for idx in range(len(imga)):\n                    img = imga[idx, :, :]\n                    img_sm = cv2.resize(img, (60, 60), interpolation=cv2.INTER_CUBIC)\n                    image[idx, :, :] = img_sm\n                imga = image\n            \n            #that allows to have the dimension in the right order !\n            imae = np.zeros((60,60,3))\n            for i in range(3):\n                imae[:,:,i] = clahe.apply((imga[i,:,:]*255).astype(np.uint8))\n                \n            imagesCroppedTest.append(np.array(imae))\n            \n            \n    imagesCroppedTest=np.array(imagesCroppedTest)\n    X_cnn = imagesCroppedTest.reshape(imagesCroppedTest.shape[0], 60, 60, 3)\n            \n    model_cnn2.load_weights('modelCnn.h5')\n    test_classif_preds = model_cnn2.predict(X_cnn)\n    \n    #Glaucoma predictions from vCDRs\n    #test_vCDRs = np.array(test_vCDRs).reshape(-1,1)\n    #features = np.concatenate((test_vCDRs,feature2),axis=1)\n    #test_classif_preds = clf2.predict_proba(features)[:,1]\n    \n# Prepare and save .csv file\ndef create_submission_csv(prediction, submission_filename='/kaggle/working/submission3.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            \n            writer.writerow({'Id': \"T{:04d}\".format(i+1), 'Predicted': '{:f}'.format(p)})\n\ncreate_submission_csv(test_classif_preds.reshape(1,-1)[0])\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":[]}]}