{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Histopathologic Cancer Detection\n\n#### Ömer Faruk Yaşar - 21527577\n#### Furkan Kaya - 21527161\n\nWe will establish an algorithm that tries to predict whether metastatic cancer is present in small image patches (96x96px) from high resolution pathology scans. Our data set will be a modified data set obtained by subtracting duplicates from the PatchCamelyon (PCam) dataset. Using Convolutional Neural Network, we aim to achieve a successful binary classification. To measure our success, we will use the area under the ROC curve metric.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Table of Content\n\n[Problem](#problem)   \n[Data Understanding](#data_understanding)   \n[Data Preparation](#data_preparation)   \n[Modeling](#modeling)   \n[Evaluation](#evaluation)   \n[References](#references)","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Problem <a class=\"anchor\" id=\"problem\"></a>[](http://)","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"For pathologists, the diagnostic procedure is a very time-consuming and effortful process. They need to make a very detailed examination on **high resolution digital pathological scans** (like 100.000x100.000px). And they are likely to overlook small metastases in this large image. Due to this small overlook, the doctor can diagnose the patient is not cancer, even though the person has cancer. Although we will not completely solve this problem, we want to detect small metastases (**binary classification problem**)  with machine learning techniques and reduce the workload of pathologists, using these high-resolution digital scans fragmented as 96x96px. In this way, pathologists will know which regions to focus on in these large digital pathological scans and will spend less time.\n\nIn short, we have **96x96px** images taken from large pathological scans, and we want to find out whether we can classify whether these images have cancerous tissue at high accuracy using **deep learning techniques**.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Data Understanding<a class=\"anchor\" id=\"data_understanding\"></a>","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"The dataset that we will use is a subset of the original [PCam dataset](https://github.com/basveeling/pcam) which in the end is derived from the [Camelyon16 Challenge dataset](https://camelyon16.grand-challenge.org/Data/), which contains 400 H&E stained whole slide images of sentinel lymph node sections that were acquired and digitized at 2 different centers using a 40x objective. The PCam's dataset including this one uses 10x undersampling to increase the field of view, which gives the resultant pixel resolution of 2.43 microns.\n\nIt consists of **277.493** color images (96x96px) extracted from **histopathologic scans of lymph node sections**. Each image is annotated with a binary label indicating presence of metastatic tissue. The train data we have here contains **220,025** images and the test set contains **57,468** images.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\n\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nimport seaborn as sns\nsns.set(style=\"darkgrid\")\nfrom PIL import Image\n\nimport torch\nfrom torch.utils.data import TensorDataset, DataLoader,Dataset\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.utils.data import DataLoader\nfrom torch.utils.data.sampler import SubsetRandomSampler\nfrom torch.optim.lr_scheduler import StepLR, ReduceLROnPlateau, CosineAnnealingLR\nimport torch.optim as optim\nfrom torch.optim import lr_scheduler\n\nimport torchvision\nimport torchvision.transforms as transforms\n\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import roc_auc_score","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_images_folder_path = \"../input/histopathologic-cancer-detection/train\"\ntest_images_folder_path = \"../input/histopathologic-cancer-detection/test\"\nlabels_file_path = '../input/histopathologic-cancer-detection/train_labels.csv'","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"labels = pd.read_csv(labels_file_path)\nlabels.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here **id**s are used to match the filenames of the images. **Label 0** means negative that is an image without cancer tissue. **Label 1** means positive that is an image with a cancerous tissue.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(5,5))\nax.pie([labels.label.value_counts()[0], labels.label.value_counts()[1]], labels=['Negative', 'Positive'], autopct='%1.1f%%');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have a training data set with a distribution of approximately **60:40** (60% negative, 40% positive).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.countplot(x='label', data=labels);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> *A positive label indicates that the center 32x32px region of a patch contains at least one pixel of tumor tissue. Tumor tissue in the outer region of the patch does not influence the label. This outer region is provided to enable the design of fully-convolutional models that do not use any zero-padding, to ensure consistent behavior when applied to a whole-slide image.* (Source: https://github.com/basveeling/pcam)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(2,5, figsize=(20,8))\nfig.suptitle('Histopathologic scans of lymph node sections')\n\n# Negative Samples\nfor i, image_id in enumerate(labels[labels.label == 0]['id'][:5]):\n    path = os.path.join(train_images_folder_path, image_id)\n    ax[0,i].imshow(Image.open(path + '.tif'))\n    box = patches.Rectangle((32,32), 32, 32, linewidth=5, edgecolor='g', facecolor='none')\n    ax[0,i].add_patch(box)\nax0 = ax[0,0].set_ylabel('Negative samples')\n\n# Positive Samples\nfor i, image_id in enumerate(labels[labels.label == 1]['id'][:5]):\n    path = os.path.join(train_images_folder_path, image_id)\n    ax[1,i].imshow(Image.open(path + '.tif'))\n    box = patches.Rectangle((32,32), 32, 32, linewidth=5, edgecolor='r', facecolor='none')\n    ax[1,i].add_patch(box)\nax1 = ax[1,0].set_ylabel('Positive samples')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see from randomly received image files, it is impossible for us to look at the picture and make a label estimate, it is a very challenging process even for the professionals of this domain. Since it is decided to classify this by taking many features together such as cell density, cell color distribution etc., finding an estimate using all these features using deep learning methods will make things much easier.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Data Preparation<a class=\"anchor\" id=\"data_preparation\"></a>","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"In machine learning overfitting and genaralization are big problems. To solve this problems you should have diverse and huge dataset. To have that kind of dataset naturally is not possible most of the time because collecting and labeling the data is challenging and expensive process. So in order to increase diversity in our dataset we apply some **data augmentation** techniques to our dataset for both avoiding overfitting and obtain better model. Random horizontal flip and vertical flip, small degrees of random rotation and normalization were applied on train images. All this operations do not increase dataset size each of them applied to images during training with 0.5 probability by that way we obtain different versions of same images in different epochs. This process increase our dataset's diversity. To do data augmentation we use pytorch's data loader and transformation classes which is very handy and easy to use. First we describe which transforms going to perform and than apply while creating train and test sets. We also split **70% to 30%** of our train data for validation.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"train_indices, validation_indices = train_test_split(labels.label, stratify=labels.label, test_size=0.3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data_transformations_train = transforms.Compose([\n    transforms.Pad(64, padding_mode='reflect'),\n    transforms.RandomHorizontalFlip(),\n    transforms.RandomVerticalFlip(),\n    transforms.RandomRotation(15),\n    transforms.ToTensor(),\n    transforms.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225])\n])\n\ndata_transformations_test = transforms.Compose([\n    transforms.Pad(64, padding_mode='reflect'),\n    transforms.ToTensor(),\n    transforms.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225])\n])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,6, figsize=(24,4))\n\nimage_path = os.path.join(train_images_folder_path, 'c18f2d887b7ae4f6742ee445113fa1aef383ed77.tif')\nimage_original = Image.open(image_path)\n\ncenter_crop = transforms.CenterCrop(64)\n\ncomposition = transforms.Compose([\n    transforms.Pad(64, padding_mode='reflect'),\n    transforms.RandomHorizontalFlip(),\n    transforms.RandomVerticalFlip(),\n    transforms.RandomRotation(15)\n])\n\nax[0].imshow(image_original)\nax[1].imshow(center_crop(image_original))\nax[2].imshow(transforms.functional.hflip(image_original))\nax[3].imshow(transforms.functional.vflip(image_original))\nax[4].imshow(transforms.functional.rotate(image_original, 15))\nax[5].imshow(composition(image_original))\n\nax0 = ax[0].set_xlabel('Original image')\nax1 = ax[1].set_xlabel('Pad')\nax2 = ax[2].set_xlabel('Horizontal Flip')\nax3 = ax[3].set_xlabel('Vertical Flip')\nax4 = ax[4].set_xlabel('Rotation')\nax5 = ax[5].set_xlabel('Composition of these transformations')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since some transformation occurs randomly in the composition section, the result may vary each time.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"class PCamDataset(Dataset):\n    def __init__(self, data_folder, data_type, transform, labels_dict={}):\n        self.data_folder = data_folder\n        self.data_type = data_type\n        self.image_files_list = [image_file_name for image_file_name in os.listdir(data_folder)]\n        self.transform = transform\n        self.labels_dict = labels_dict\n        if self.data_type == 'train':\n            self.labels = [labels_dict[image_file_name.split('.')[0]] for image_file_name in self.image_files_list]\n        else:\n            self.labels = [0 for _ in range(len(self.image_files_list))]\n\n    def __len__(self):\n        return len(self.image_files_list)\n\n    def __getitem__(self, index):\n        image_file_name = os.path.join(self.data_folder, self.image_files_list[index])\n        image = Image.open(image_file_name)\n        image = self.transform(image)\n        image_id = self.image_files_list[index].split('.')[0]\n        if self.data_type == 'train':\n            label = self.labels_dict[image_id]\n        else:\n            label = 0\n        return image, label","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image_labels_dict = { image:label for image, label in zip(labels.id, labels.label) }\n\ndataset = PCamDataset(data_folder=train_images_folder_path, data_type='train', transform=data_transformations_train, labels_dict=image_labels_dict)\ntest_set = PCamDataset(data_folder=test_images_folder_path, data_type='test', transform=data_transformations_test)\n\ntrain_sampler = SubsetRandomSampler(list(train_indices.index))\nvalid_sampler = SubsetRandomSampler(list(validation_indices.index))\n\nbatch_size = 64\n\ntrain_loader = DataLoader(dataset, batch_size=batch_size, sampler=train_sampler)\nvalid_loader = DataLoader(dataset, batch_size=batch_size, sampler=valid_sampler)\ntest_loader = DataLoader(test_set, batch_size=batch_size)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Modeling<a class=\"anchor\" id=\"modeling\"></a>","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"As we mentioned earlier parts of this notebook, we try to solve this binary classification problem with deep learning approach to be more specific we are going to use Convolutional Neural Networks which is widely used in visual classification, detection problems lately. There are lots of popular CNN architectures with pretrained models which are VGG, ResNet etc. And we are going to choose and adopt from one of them to our problem. By doing that, we are avoiding from trainin a network from scratch and also increase our accuracy. So it's efficent in terms of both time and accuracy. This process called as transfer learning in literature. It's widely used in deep learning solutions. There are two versions of transfer learning; you can train the weights that you tuned to your network or freeze it. We are going to try both approach in our experiments. As a choice of network we choose ResNet (Residual Networks) because of several reasons; ResNet are widely used in many visual problem in machine learning because in deep learning while the layers increased, we start to have a problems while trying to optimize and train our networks due to vanishing/exploiding gradients. ResNet offer a solution for this problem with residual connections between layers and avoid to lose earlier features that we learned in earlier layers. We decided on using ResNet with pretrained weigths on imagenet dataset but ResNet has several versions according to layer size. We train and evaluate on validation set with ResNet18,ResNet50 and ResNet101 and observe that all of them perform almost same on our problem so we decided to continue our project with ResNet18 because it's the shallowest network among them so it's easier to train. Than after we try to pick best hyperparameters for our network. Hyperparameter is a parameter that defined before training process and can not updated during this process so it is something that we can not learn during our training.Batch size, learning rate , loss function, epoch number are hyperparameters that we have in our project so we are trying to pick optimal one by doing experiments with different values for each of them.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"model = torchvision.models.resnet18(pretrained=True)\nfor i, param in model.named_parameters():\n    param.requires_grad = False","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We add a 2-feature linear layer so that the model can do binary classification into the last fully connected layer.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"num_features = model.fc.in_features\nmodel.fc = nn.Linear(num_features, 2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.cuda()\ncriterion = nn.BCEWithLogitsLoss()\noptimizer = optim.SGD(model.fc.parameters(), lr=0.001, momentum=0.9)\nexp_lr_scheduler = lr_scheduler.StepLR(optimizer, step_size=5, gamma=0.1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n_epochs = 4\npatience = 5\np = 0\nstop = False\nvalid_loss_min = np.Inf\n\ntrain_loss_epoch = []\nval_loss_epoch = []\nval_auc_epoch = []\n\n\nfor epoch in range(1, n_epochs+1):\n\n    train_loss = []\n    train_auc = []\n\n    for batch_i, (data, target) in enumerate(train_loader):\n\n        data, target = data.cuda(), target.cuda()\n\n        optimizer.zero_grad()\n        output = model(data)\n        loss = criterion(output[:,1], target.float())\n        train_loss.append(loss.item())\n        \n        a = target.data.cpu().numpy()\n        b = output[:,-1].detach().cpu().numpy()\n        train_auc.append(roc_auc_score(a, b))\n\n        loss.backward()\n        optimizer.step()\n    \n    exp_lr_scheduler.step()\n    \n    train_loss_epoch.append(np.mean(train_loss))\n    \n    model.eval()\n    \n    val_loss = []\n    val_auc = []\n    \n    for batch_i, (data, target) in enumerate(valid_loader):\n        data, target = data.cuda(), target.cuda()\n        output = model(data)\n\n        loss = criterion(output[:,1], target.float())\n\n        val_loss.append(loss.item()) \n        a = target.data.cpu().numpy()\n        b = output[:,-1].detach().cpu().numpy()\n        val_auc.append(roc_auc_score(a, b))\n\n    val_loss_epoch.append(np.mean(val_loss))\n    val_auc_epoch.append(np.mean(val_auc))\n    \n    print(f'Epoch {epoch}, train loss: {np.mean(train_loss):.4f}, valid loss: {np.mean(val_loss):.4f}, train auc: {np.mean(train_auc):.4f}, valid auc: {np.mean(val_auc):.4f}')\n    \n    valid_loss = np.mean(val_loss)\n    if valid_loss <= valid_loss_min:\n        print('Validation loss decreased ({:.6f} --> {:.6f}).  Saving model ...'.format(\n        valid_loss_min,\n        valid_loss))\n        torch.save(model.state_dict(), 'model.pt')\n        valid_loss_min = valid_loss\n        p = 0\n\n    # check if validation loss didn't improve\n    if valid_loss > valid_loss_min:\n        p += 1\n        print(f'{p} epochs of increasing val loss')\n        if p > patience:\n            print('Stopping training')\n            stop = True\n            break        \n    if stop:\n        break","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Evaluation<a class=\"anchor\" id=\"evaluation\"></a>","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Our evaluation metric is the **Area Under the Curve (AUC) score**. It is one of the most important evaluation metrics for checking any classification model’s performance. It tells how much model is capable of distinguishing between classes. Higher AUC, better the model is at distinguishing between healthy tissue and cancerous tissue.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"model = torchvision.models.resnet18(pretrained=True)\n\nnum_features = model.fc.in_features\nmodel.fc = nn.Linear(num_features, 2)\n\nmodel.cuda()\nmodel.load_state_dict(torch.load('model.pt'))\nmodel.eval();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(train_loss_epoch, label='Train Loss')\nplt.plot(val_loss_epoch, label='Validation Loss')\nplt.xlabel(\"Epochs\")\nplt.ylabel(\"Loss\")\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(val_auc_epoch, label='Validation AUC')\nplt.xlabel(\"Epochs\")\nplt.ylabel(\"Area Under the Curve\")\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"preds = []\nfor batch_i, (data, target) in enumerate(test_loader):\n    data = data.cuda()\n    output = model(data)\n    pr = output[:,1].detach().cpu().numpy()\n    for i in pr:\n        preds.append(i)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_preds = pd.DataFrame({'imgs': test_set.image_files_list, 'preds': preds})\ntest_preds['imgs'] = test_preds['imgs'].apply(lambda x: x.split('.')[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = pd.read_csv('../input/histopathologic-cancer-detection/sample_submission.csv')\nsub = pd.merge(sub, test_preds, left_on='id', right_on='imgs')\nsub = sub[['id', 'preds']]\nsub.columns = ['id', 'label']\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## References<a class=\"anchor\" id=\"references\"></a>","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"[1] B. S. Veeling, J. Linmans, J. Winkens, T. Cohen, M. Welling. \"Rotation Equivariant CNNs for Digital Pathology\". arXiv:1806.03962\n\n[2] Ehteshami Bejnordi et al. Diagnostic Assessment of Deep Learning Algorithms for Detection of Lymph Node Metastases in Women With Breast Cancer. JAMA: The Journal of the American Medical Association, 318(22), 2199–2210. doi:jama.2017.14585\n\n[3] He, Kaiming, Xiangyu Zhang, Shaoqing Ren and Jian Sun. “Deep Residual Learning for Image Recognition.” 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2016): 770-778.\n\n[4] Encyclopedia of Science Education, 2015. Transfer of Learning. pp.1079-1079.\n\nhttps://www.kaggle.com/qitvision/a-complete-ml-pipeline-fast-ai\n\nhttps://www.kaggle.com/artgor/simple-eda-and-model-in-pytorch\n\nhttps://www.kaggle.com/ashishpatel26/hc-detection-using-pytorch-resnet-101\n\nhttps://www.kaggle.com/abhinand05/histopathologic-cancer-detection-using-cnns","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"**Disclaimer!** <font color='grey'>This notebook was prepared by <student name(s)> as a term project for the *BBM469 - Data Intensive Applications Laboratory* class. The notebook is available for educational purposes only. There is no guarantee on the correctness of the content provided as it is a student work.\n\nIf you think there is any copyright violation, please let us [know](https://forms.gle/BNNRB2kR8ZHVEREq8). \n</font>","execution_count":null}],"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":4,"nbformat_minor":4}