{"cells":[{"metadata":{},"cell_type":"markdown","source":"## 0. List files in input_folder","execution_count":null},{"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","execution_count":null},{"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","execution_count":null},{"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","execution_count":null},{"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)`","execution_count":null},{"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[:,:,0])\n    \ncls = PANDADataset(df, 1)\n%time cls[0][0].size(), cls[0][1].size(), len(cls)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(cls[300][1])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3. Build the model\n\nAdapted from https://discuss.pytorch.org/t/unet-implementation/426","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"dataLoader = DataLoader(cls, batch_size=8, shuffle=True, num_workers=8)\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","execution_count":null},{"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 = 1\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","execution_count":null},{"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":{"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] = PANDADataset(df)\n    dataLoader[phase] = DataLoader(dataset[phase], batch_size=batch_size, \n                                shuffle=True, num_workers=8, pin_memory=True)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"markdown","source":"#visualize a single example to verify that it is correct\nimg, patch_mask_weight = dataset[\"train\"][7]\nfig, ax = plt.subplots(1, 2, figsize=(10,4))  # 1 row, 2 columns\nprint(img.shape, patch_mask_weight.shape)\n#build output showing original patch  (after augmentation), class = 1 mask, weighting mask, overall mask (to see any ignored classes)\nax[0].imshow(np.moveaxis(img.numpy(),0,-1))\nax[1].imshow(patch_mask_weight)\nplt.show()","execution_count":null},{"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.","execution_count":null},{"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()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Train model","execution_count":null},{"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    train_accuracy = 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.to(device,dtype = torch.int64)\n        \n        # zero the parameter gradients\n        optim.zero_grad()\n        \n        # =========================== forward + backward + optimize ===========================\n        outputs = model(inputs)\n        #_, outputs = torch.max(model(inputs), axis = 1)\n        #outputs = torch.argmax(model(inputs), axis = 1)\n        \n        ## =========================== Loss computation ===========================\n        loss = criterion(outputs, labels)\n        loss.sum().backward()\n        optim.step()       \n        \n        ## =========================== Accuracy computation ========================================        \n        # return the indices of max values along rows in softmax probability output\n        predicted = torch.argmax(outputs, axis = 1)\n        \n        # number of pixel in the batch\n        total_train += labels.nelement()        \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        #print(labels == predicted)\n        correct_train += (labels == predicted).sum().item()        \n        \n        # =========================== print statistics ===========================\n        running_loss += loss.mean()\n        \n        train_accuracy = correct_train / total_train\n        \n        if i % 300 == 299:    # print every 2000 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 / 300, train_accuracy))\n            running_loss = 0.0\nprint('==============================================================================================================')\nprint('============================================== Finished Training =============================================')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Prediction","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"df = dd.io.load(hdf5_dir)\ncls_test = PANDADataset(df[ : 800], 1)\n\nplt.figure(figsize = (20, 10))\nfor k in range(5):\n    idx = np.random.randint(0, 800)\n    a = cls_test[idx][0].permute((1, 2, 0)).detach().squeeze().cpu().numpy()\n    b = cls_test[idx][1].detach().squeeze().cpu().numpy()\n    cmap =  matplotlib.colors.ListedColormap(['black', 'gray', 'green', 'yellow', 'orange', 'red'])\n    plt.subplot(2, 5, k+1), plt.imshow(a)\n    plt.subplot(2, 5, k+6), 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":"data_inp = []\npredicts = []\ntrue_mask = []\ncls_test = PANDADataset(df[ : 80], 1)\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.to(device,dtype = torch.int64) ## type('torch.LongTensor').to(device)\n    \n    predict = torch.argmax(model(inputs), axis = 1)\n    \n    ## append\n    predicts += predict\n    data_inp += inputs\n    true_mask += labels\n    ## freeze\n    del inputs, labels, predict, data\n    \nprint(len(predicts), predicts[0].shape)\nprint(len(data_inp), data_inp[0].shape)\nprint(len(true_mask), true_mask[0].shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1, 4, figsize=(25, 5.5))\nfor k in range(4):\n    c = predicts[k].detach().squeeze().cpu().numpy()\n    ax[k].imshow(c, cmap = cmap), ax[k].set_title('predict %s'%(k+1))\nplt.show()\n\nfig, ax = plt.subplots(1, 4, figsize=(25, 5.5))\nfor k in range(4):\n    d = true_mask[k].detach().squeeze().cpu().numpy()\n    ax[k].imshow(d, cmap = cmap), ax[k].set_title('true_mask %s'%(k+1))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plt.imshow(torch.argmax(predicts[0], axis = 0).detach().squeeze().cpu().numpy())\ntorch.argmax(predicts[0], axis = 0)[0]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Evaluation\n\n(I will do it later in the next few days) !!!","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Grading","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"ID_list = os.listdir(mask_dir)\nID_list = [u.replace('_mask.tiff', '') for u in ID_list]\n\nID = ID_list[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def split_get_tiles(img_id, crit = 0.0005, size=512, n_tiles=36):    \n    \"\"\"\n    ==================================================================================================\n    Input:  img_id (str): image_id from the train dataset, such as '004dd32d9cd167d9cc31c13b704498af'  \n            crit (float) in (0, 1): the proportion of the dark_region over whole image (size 256 x 256)\n            size (int) : image size\n            n_tiles : number of tiles\n    return: \n            list of (img_id, x_start, x_end, y_start, y_end) images size 512x512    \n            ==========================================================================================\n    writen by Nhan\n    ==================================================================================================\n    \"\"\"\n    img = skimage.io.MultiImage(os.path.join(data_dir, f'{img_id}.tiff'))[0]\n    tile_size = 512\n    h, w = img.shape[: 2]\n    nc = int(w / 512)\n    nr = int(h / 512)\n    img_ls = []\n    coord_ls = []\n    S_img_tile = 512*512*3\n    \n    for i in range(nr):\n        for j in range(nc):\n            x_start, y_start = int(i*512), int(j*512)\n            image_dt = img[ x_start : x_start + 512, y_start : y_start + 512 , :]\n            if (image_dt.min() < 185):\n                count = len(image_dt[image_dt <= 121])\n                if count/(S_img_tile) >= crit:\n                    image_dt = cv2.resize(image_dt, (size, size), interpolation = cv2.INTER_AREA)\n                    img_ls.append(image_dt)\n                    del image_dt, x_start, y_start\n\n    ## choose n_tiles image has a best-view_range \n    img3_dt_ = np.array(img_ls)\n    idxs_dt_ = np.argsort(img3_dt_.reshape(img3_dt_.shape[0],-1).sum(-1))[:n_tiles]\n    \n    ## attach\n    list_image = []\n    for final_index in idxs_dt_:\n        list_image.append(img_ls[final_index])\n    for i in range(8): \n        yield list_image[i: i+8]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Function to make tiles from one image\n# ...\n\n# Funtion to calculate the isup from the list of outputs.\ndef ISUP(result):\n    # result: a list of masks\n    # Translation matrix of gleason scores to isup\n    import numpy as np\n    isup_mat = np.array([[1, 2, 4],[3, 4, 5],[4, 5, 5]])\n    # calculate the most dominant gleason score\n    p = np.zeros(6)\n    for mask in result:\n        for i in range(3, 6):\n            p[i] += len(np.nonzeros((mask==i)*1)) \n    gscore1 = max(argmax(p), 3)\n    p[gscore1] = 0\n    if argmax(p) ==0:\n        gscore2 = gscore1\n    else:\n        gscore2 = max(argmax(p), 3)\n        \n    return isup_mat[gscore1-3, gscore2-3]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Configuration\nimport openslide\nimg_id = '6d1a11077fe4183a4109d649cf319923'\n# Load image\n#osh = openslide.OpenSlide(file)\n\n# Load model\nmodel = 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)\ncheckpoint = torch.load(PATH)\nmodel.load_state_dict(checkpoint['model_state_dict'])\n\n# Define and add device\ngpuid = 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')\n    \nmodel.to(device)\n\n# Evaluation mode\nmodel.eval()\n\n# Cut image into tiles\ntiles = split_get_tiles(img_id, crit = 0.0005, size=512, n_tiles=36)\n\nresult = []\n# Make masks\nfor inputs in tiles: \n    predicts = torch.argmax(model(inputs), axis = 1)\n    result = result.append(predicts)\n\n# Make ISUP grade\nisup = ISUP(result)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}