{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Detecting Prostate Cancer on Images of Prostate Tissue Samples using Machine Learning\nIn this notebook, we are going build a machine learning model that help detect Prostate Cancer based on images of Prostate Tissue samples and estimate severity of the disease.\n\n## 1. Problem definition\n\n> How well can we detect and estimate severity of Prostate Cancer given images of Prostate Tissue samples ?\n\nFor more description about the problem and its importance check: https://www.kaggle.com/c/prostate-cancer-grade-assessment/overview.\n\n## 2. Data\n\nThe data is downloaded from Kaggle Prostate Cancer Grade Assessment competition: \nhttps://www.kaggle.com/c/prostate-cancer-grade-assessment/data.\n\nThere are 4 main datasets:\n\n* [train/test].csv contains ID code for the images.\n\n\n* train_images: The images. Each is a large multi-level tiff file. You can expect roughly 1,000 images in the hidden test set. Note that slightly different procedures were in place for the images used in the test set than the training set. Some of the training set images have stray pen marks on them, but the test set slides are free of pen marks.\n\n* train_label_masks: Segmentation masks showing which parts of the image led to the ISUP grade. Not all training images have label masks, and there may be false positives or false negatives in the label masks for a variety of reasons. These masks are provided to assist with the development of strategies for selecting the most useful subsamples of the images.\n\n## 3. Evaluation\n\nSubmissions are scored based on the quadratic weighted kappa, which measures the agreement between two outcomes. This metric typically varies from 0 (random agreement) to 1 (complete agreement). In the event that there is less agreement than expected by chance, the metric may go below 0.\n\nFor more on the evaluation of this project check: https://www.kaggle.com/c/prostate-cancer-grade-assessment/overview/evaluation.\n\n## 4. Features\n\nKaggle provides details of all features of the data sets. Check: \nhttps://www.kaggle.com/c/prostate-cancer-grade-assessment/data"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# print out the names of the first 2 image_files (total = 4 images for train_imgaes & train_label_masks) with the train, test, submission.csv files & 5 file.hdf5\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames[:2]:\n        print(os.path.join(dirname, filename))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Naming the directories"},{"metadata":{"trusted":true},"cell_type":"code","source":"import cv2\nimport openslide\nimport skimage.io\nimport random\nimport seaborn as sns\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport PIL\nfrom IPython.display import Image, display\n\nBASE_PATH = '../input/prostate-cancer-grade-assessment'\ndata_dir = f'{BASE_PATH}/train_images'\nmask_dir = f'{BASE_PATH}/train_label_masks'\nhdf5_dir = r'/kaggle/input/radboud-database/radboud_tiles_coordinates.h5'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 1. Load database"},{"metadata":{"trusted":true},"cell_type":"code","source":"import deepdish as dd\n\ndf = dd.io.load(hdf5_dir)\nlen(df)//36, len(df[0]), df[0], len(df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 2. Create the class to load PANDA_dataset with this database"},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_data_and_mask(ID, coordinates, level = 1):\n    \"\"\"\n    Input args:\n        ID (str): img_id from the dataset\n        coordinates (list of int): list of coordinates, includes: [x_start, x_end, y_start, y_end] from h5.database\n        level (={0, 1, 2}) : level of images for loading with skimage\n    Return: 3D tiles shape 512x512 of the mask images and data images w.r.t the input_coordinates, ID and level\n    \"\"\"\n    data_img = skimage.io.MultiImage(os.path.join(data_dir, f'{ID}.tiff'))[level]\n    mask_img = skimage.io.MultiImage(os.path.join(mask_dir, f'{ID}_mask.tiff'))[level]\n    coordinates = [coordinate // 2**(2*level) for coordinate in coordinates]\n    data_tile = data_img[coordinates[0]: coordinates[1], coordinates[2]: coordinates[3], :]\n    mask_tile = mask_img[coordinates[0]: coordinates[1], coordinates[2]: coordinates[3], :]\n    data_tile = cv2.resize(data_tile, (512, 512))\n    mask_tile = cv2.resize(mask_tile, (512, 512))\n    del data_img, mask_img\n    \n    # Load and return small image\n    return data_tile, mask_tile","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### First trying with the first `3500 (img_id)` or `126000 (tiles)`"},{"metadata":{"trusted":true},"cell_type":"code","source":"from torch.utils.data import Dataset, DataLoader\nimport torch\n\nclass PANDADataset(Dataset):\n    def __init__(self, df, level = 2, transform=None):\n        self.df = df\n        self.level = level\n        self.transform = transform\n\n    def __len__(self):\n        return len(self.df)\n    def __getitem__(self, index, level = 2):\n        ID = self.df[index][0]\n        coordinate = self.df[index][1: ]\n        image, mask = load_data_and_mask(ID, coordinate, level)\n        \n        return torch.tensor(image).permute(2, 0, 1), torch.tensor(mask).permute(2, 0, 1)[0]\n    \ncls = PANDADataset(df[: 60000], 1)\n%time cls[0][0].size(), cls[0][1].size(), len(cls)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3. Build the model\n\nAdapted from https://discuss.pytorch.org/t/unet-implementation/426"},{"metadata":{"trusted":true},"cell_type":"code","source":"dataLoader = DataLoader(cls, batch_size=8, shuffle=True, num_workers=8)\ndel df, cls\n\nimport torch\nfrom torch import nn\nimport torch.nn.functional as F\n\n\nclass UNet(nn.Module):\n    def __init__(self, in_channels=1, n_classes=2, depth=5, wf=6, padding=False,\n                 batch_norm=False, up_mode='upconv'):\n        \"\"\"\n        Implementation of\n        U-Net: Convolutional Networks for Biomedical Image Segmentation\n        (Ronneberger et al., 2015)\n        https://arxiv.org/abs/1505.04597\n        Using the default arguments will yield the exact version used\n        in the original paper\n        Args:\n            in_channels (int): number of input channels\n            n_classes (int): number of output channels\n            depth (int): depth of the network\n            wf (int): number of filters in the first layer is 2**wf\n            padding (bool): if True, apply padding such that the input shape\n                            is the same as the output.\n                            This may introduce artifacts\n            batch_norm (bool): Use BatchNorm after layers with an\n                               activation function\n            up_mode (str): one of 'upconv' or 'upsample'.\n                           'upconv' will use transposed convolutions for\n                           learned upsampling.\n                           'upsample' will use bilinear upsampling.\n        \"\"\"\n        super(UNet, self).__init__()\n        assert up_mode in ('upconv', 'upsample')\n        self.padding = padding\n        self.depth = depth\n        prev_channels = in_channels\n        self.down_path = nn.ModuleList()\n        for i in range(depth):\n            self.down_path.append(UNetConvBlock(prev_channels, 2**(wf+i),\n                                                padding, batch_norm))\n            prev_channels = 2**(wf+i)\n\n        self.up_path = nn.ModuleList()\n        for i in reversed(range(depth - 1)):\n            self.up_path.append(UNetUpBlock(prev_channels, 2**(wf+i), up_mode,\n                                            padding, batch_norm))\n            prev_channels = 2**(wf+i)\n\n        self.last = nn.Conv2d(prev_channels, n_classes, kernel_size=1)\n\n    def forward(self, x):\n        blocks = []\n        for i, down in enumerate(self.down_path):\n            x = down(x)\n            if i != len(self.down_path)-1:\n                blocks.append(x)\n                x = F.avg_pool2d(x, 2)\n\n        for i, up in enumerate(self.up_path):\n            x = up(x, blocks[-i-1])\n\n        return self.last(x)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class UNetConvBlock(nn.Module):\n    def __init__(self, in_size, out_size, padding, batch_norm):\n        super(UNetConvBlock, self).__init__()\n        block = []\n\n        block.append(nn.Conv2d(in_size, out_size, kernel_size=3,\n                               padding=int(padding)))\n        block.append(nn.ReLU())\n        if batch_norm:\n            block.append(nn.BatchNorm2d(out_size))\n\n        block.append(nn.Conv2d(out_size, out_size, kernel_size=3,\n                               padding=int(padding)))\n        block.append(nn.ReLU())\n        if batch_norm:\n            block.append(nn.BatchNorm2d(out_size))\n\n        self.block = nn.Sequential(*block)\n\n    def forward(self, x):\n        out = self.block(x)\n        return out","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class UNetUpBlock(nn.Module):\n    def __init__(self, in_size, out_size, up_mode, padding, batch_norm):\n        super(UNetUpBlock, self).__init__()\n        if up_mode == 'upconv':\n            self.up = nn.ConvTranspose2d(in_size, out_size, kernel_size=2,\n                                         stride=2)\n        elif up_mode == 'upsample':\n            self.up = nn.Sequential(nn.Upsample(mode='bilinear', scale_factor=2),\n                                    nn.Conv2d(in_size, out_size, kernel_size=1))\n\n        self.conv_block = UNetConvBlock(in_size, out_size, padding, batch_norm)\n\n    def center_crop(self, layer, target_size):\n        _, _, layer_height, layer_width = layer.size()\n        diff_y = (layer_height - target_size[0]) // 2\n        diff_x = (layer_width - target_size[1]) // 2\n        return layer[:, :, diff_y:(diff_y + target_size[0]), diff_x:(diff_x + target_size[1])]\n\n    def forward(self, x, bridge):\n        up = self.up(x)\n        crop1 = self.center_crop(bridge, up.shape[2:])\n        out = torch.cat([up, crop1], 1)\n        out = self.conv_block(out)\n\n        return out","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### `Unet`-params & training params\n\nThese parameters get fed directly into the UNET class, and more description of them can be discovered there\n\nBut here, I will try with epochs = 3"},{"metadata":{"trusted":true},"cell_type":"code","source":"# --- Unet params\nn_classes= 6    # number of classes in the data mask that we'll aim to predict\n\n\nin_channels = 3  # input channel of the data, RGB = 3\npadding = True   # should levels be padded\ndepth = 5        # depth of the network \nwf = 2           # wf (int): number of filters in the first layer is 2**wf, was 6\nup_mode = 'upconv' #should we simply upsample the mask, or should we try and learn an interpolation \nbatch_norm = True #should we use batch normalization between the layers\n\n# --- training params\n\nbatch_size = 8\npatch_size = 512\nnum_epochs = 2\nedge_weight = 1.1 # edges tend to be the most poorly segmented given how little area they occupy in the training set, this paramter boosts their values along the lines of the original UNET paper\nphases = [\"train\",\"val\"] # how many phases did we create databases for?\nvalidation_phases= [\"val\"] # when should we do valiation? note that validation is time consuming, so as opposed to doing for both training and validation, we do it only for vlaidation at the end of the epoch","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 4. Decide what divice to run the model"},{"metadata":{"trusted":true},"cell_type":"code","source":"gpuid = 0\nif(torch.cuda.is_available()):\n    print(torch.cuda.get_device_properties(gpuid))\n    torch.cuda.set_device(gpuid)\n    device = torch.device(f'cuda:{gpuid}')\nelse:\n    device = torch.device(f'cpu')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 5. Fit the model according to the paramters specified above and copy it to the GPU.\n\nThen finally print out the number of trainable parameters."},{"metadata":{"trusted":true},"cell_type":"code","source":"model = UNet(n_classes = n_classes, in_channels = in_channels, \n             padding = padding, depth = depth, wf = wf, \n             up_mode = up_mode, batch_norm = batch_norm).to(device)\nprint(f\"total params: \\t{sum([np.prod(p.size()) for p in model.parameters()])}\")\n\noptim = torch.optim.Adam(model.parameters()) #adam is going to be the most robust\ncriterion = nn.CrossEntropyLoss(reduce=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Train model"},{"metadata":{"trusted":true},"cell_type":"code","source":"import time\nprint('========================================== Training started ==========================================')\nfor epoch in range(num_epochs):\n    print('======================================================================================================')\n    # model.train()  # Set model to training mode\n    running_loss = 0.0\n    total_train = 0\n    correct_train = 0\n    t0 = time.time()\n    \n    for i, data in enumerate(dataLoader, 0):\n        inputs, labels = data\n        inputs = inputs.to(device,dtype = torch.float) \n        labels = labels.type('torch.LongTensor').to(device)\n        \n        # zero the parameter gradients\n        optim.zero_grad()\n\n        # =========================== forward + backward + optimize ===========================\n        outputs = model(inputs)\n        \n        ## =========================== Loss computation ===========================\n        loss = criterion(outputs, labels)\n        loss.sum().backward()\n        optim.step()\n        \n        ## =========================== Accuracy computation ========================================\n        \n        # return the indices of max values along rows in softmax probability output\n        _, predicted = torch.max(outputs, 1)\n        \n        # number of pixel in the batch\n        total_train += labels.nelement()\n        \n        # count of the number of times the neural network has produced a correct output, and \n        # we take an accumulating sum of these correct predictions so that we can determine the accuracy of the network.\n        correct_train += (predicted == labels).sum().item()\n        \n        # =========================== print statistics ===========================\n        running_loss += loss.mean()\n        train_accuracy = correct_train / total_train\n        \n        if i % 100 == 99:    # print every 100 mini-batches\n            t1 = time.time()\n            h = (t1 - t0) // 3600\n            m = (t1 - t0 - h*3600) // 60\n            s = (t1 - t0) % 60\n            print('Eps %02d, upto %05d mnbch; after %02d (hours) %02d (minutes) and %02d (seconds);  train_loss = %.3f, train_acc = %.3f'%\n                  (epoch + 1, i + 1, h, m, s, running_loss / 100, train_accuracy))\n            running_loss = 0.0\nprint('======================================================================================================')\nprint('========================================== Finished Training =========================================')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Visualization"},{"metadata":{"trusted":true},"cell_type":"code","source":"cls_test = PANDADataset(df[: 40], 1)\na = cls_test[0][0].detach().squeeze().cpu().numpy()\nb = cls_test[0][1].detach().squeeze().cpu().numpy()\nplt.subplot(121), plt.imshow(a.reshape(512, 512, 3))\nplt.subplot(122), plt.imshow(b, cmap = cmap)\nplt.show()\nprint(a.min(), a.max(), b.min(), b.max())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"predicts = []\ndataLoader_test = DataLoader(cls_test, batch_size=8, shuffle=True, num_workers=8)\nfor i, data in enumerate(dataLoader_test, 0):\n    inputs, labels = data\n    inputs = inputs.to(device,dtype = torch.float) \n    labels = labels.type('torch.LongTensor').to(device)\n    predict = model(inputs)\n    predicts += predict\nlen(X), X[0].shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"c = X[0][0].detach().squeeze().cpu().numpy()\ncmap =  matplotlib.colors.ListedColormap(['black', 'gray', 'green', 'yellow', 'orange', 'red'])\nplt.imshow(c, cmap = cmap)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Prediction"},{"metadata":{"trusted":true},"cell_type":"markdown","source":"dataset={}\ndataLoader={}\nfor phase in phases: #now for each of the phases, we're creating the dataloader\n                     #interestingly, given the batch size, i've not seen any improvements from using a num_workers>0\n    \n    dataset[phase]=Dataset(f\"./{dataname}_{phase}.pytable\", img_transform=img_transform , mask_transform = mask_transform ,edge_weight=edge_weight)\n    dataLoader[phase]=DataLoader(dataset[phase], batch_size=batch_size, \n                                shuffle=True, num_workers=8, pin_memory=True)"},{"metadata":{},"cell_type":"markdown","source":"#### Evaluation\n\n(I will do it later in the next few days) !!!"}],"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}