{"cells":[{"metadata":{},"cell_type":"markdown","source":"### References"},{"metadata":{},"cell_type":"markdown","source":"[Pytorch FCN-Resnet50 in 20 minute](https://www.kaggle.com/finlay/pytorch-fcn-resnet50-in-20-minute)"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"scrolled":false},"cell_type":"code","source":"import numpy as np\nimport pathlib\nimport os\nimport time\nimport cv2\nimport gc\nfrom tqdm.notebook import tqdm\nimport warnings\n\n# data structure\nimport pandas as pd\n\n# graphics\nimport matplotlib.pyplot as plt\n\n# tiff file\nimport tifffile\nimport rasterio\nfrom rasterio.windows import Window\nfrom PIL import Image\n\n# transforms\nimport albumentations as A\nimport torch.utils.data as D\n\nwarnings.filterwarnings('ignore')\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"BASE_DIR = '../input/hubmap-kidney-segmentation'\nTRAIN_DIR = BASE_DIR + \"/train\"\nTEST_DIR = BASE_DIR + \"/test\"\nTEST_PRED_DIR = \"../input/submission-predictions\"\nWINDOW = 1024\nNEW_SIZE = 512\nOVERLAP = 32\nTHRESHOLD = 50","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Tools"},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def rle2mask(mask_rle, shape=(NEW_SIZE, NEW_SIZE)):\n    '''\n    mask_rle: run-length as string formated (start length)\n    shape: (height,width) of array to return\n    Returns numpy array, 1 - mask, 0 - background\n\n    '''\n    s = mask_rle.split()\n    starts, lengths = [np.asarray(x, dtype=int)\n                       for x in (s[0:][::2], s[1:][::2])]\n    starts -= 1\n    ends = starts + lengths\n    img = np.zeros(shape[0]*shape[1], dtype=np.uint8)\n    for lo, hi in zip(starts, ends):\n        img[lo:hi] = 1\n\n    del s, starts, lengths, ends\n    gc.collect()\n    return img.reshape(shape).T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def make_grid(shape, window=WINDOW, min_overlap=OVERLAP):\n    \"\"\"\n        Return Array of size (N,4), where N - number of tiles,\n        2nd axis represente slices: x1,x2,y1,y2\n    \"\"\"\n    x, y = shape\n\n    nx = x // (window - min_overlap) + 1\n    x1 = np.linspace(0, x, num=nx, endpoint=False, dtype=np.int64)\n    x1[-1] = x - window\n    x2 = (x1 + window).clip(0, x)\n\n    ny = y // (window - min_overlap) + 1\n    y1 = np.linspace(0, y, num=ny, endpoint=False, dtype=np.int64)\n    y1[-1] = y - window\n    y2 = (y1 + window).clip(0, y)\n\n    slices = np.zeros((nx, ny, 4), dtype=np.int64)\n    for i in range(nx):\n        for j in range(ny):\n            slices[i, j] = x1[i], x2[i], y1[j], y2[j]\n\n    del x, y, x1, x2, y1, y2\n    gc.collect()\n    return slices.reshape(nx*ny, 4)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def read_image(item, path, scale=None, verbose=1, with_mask=True):\n    # image\n    image_id = item.get('id')\n    image = tifffile.imread(os.path.join(path, f\"{image_id}.tiff\"))\n    if len(image.shape) == 5:\n        image = image.squeeze().transpose(1, 2, 0)\n\n    # mask\n    if with_mask:\n        column_name = 'encoding' if path == TRAIN_DIR else 'predicted'\n        mask = rle2mask(item.get(column_name), (image.shape[1], image.shape[0]))\n\n    if verbose:\n        print(f\"[{image_id}] Image shape: {image.shape}\")\n        if with_mask:\n            print(f\"[{image_id}] Mask shape: {mask.shape}\")\n\n    # scale\n    if scale:\n        new_size = (image.shape[1] // scale, image.shape[0] // scale)\n        image = cv2.resize(image, new_size)\n        if with_mask:\n            mask = cv2.resize(mask, new_size)\n\n        if verbose:\n            print(f\"[{image_id}] Resized Image shape: {image.shape}\")\n            if with_mask:\n                print(f\"[{image_id}] Resized Mask shape: {mask.shape}\")\n                \n    if with_mask:\n        return image, mask\n    \n    return image, None","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def read_mask_only(item, path, verbose=1):\n    # special function with memory optimization\n    image = tifffile.imread(os.path.join(path, f\"{item.get('id')}.tiff\"))\n    if len(image.shape) == 5:\n        image = image.squeeze().transpose(1, 2, 0)\n    image_shape = image.shape\n    del image\n    gc.collect()\n\n    column_name = 'encoding' if path == TRAIN_DIR else 'predicted'\n    mask = rle2mask(item.get(column_name), (image_shape[1], image_shape[0]))\n    if verbose:\n        print(f\"[{image_id}] Mask shape: {mask.shape}\")\n\n    return mask","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def plot_image(item, dir_path, scale=None, verbose=1):\n    image, mask = read_image(item, dir_path, scale, verbose)\n\n    plt.figure(figsize=(25, 10))\n\n    plt.subplot(131)\n    plt.imshow(image)\n    plt.title(f\"Image {item.get('id')}\")\n    plt.axis(\"off\")\n\n    plt.subplot(132)\n    plt.imshow(mask, cmap=\"hot\")\n    plt.title(\"Mask\")\n    plt.axis(\"off\")\n\n    plt.subplot(133)\n    plt.imshow(image)\n    plt.imshow(mask, cmap=\"hot\", alpha=0.5)\n    plt.title(\"Image + mask\")\n    plt.axis(\"off\")\n\n    plt.savefig(f\"full_image_mask_{item.get('id')}.png\", bbox_inches='tight')\n    plt.show()\n    del image, mask","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def to_tensor(x, **kwargs):\n    return x.transpose(2, 0, 1).astype('float32')\n\n\ndef get_preprocessing():\n    return A.Compose([A.Lambda(image=to_tensor, mask=to_tensor)])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"identity = rasterio.Affine(1, 0, 0, 0, 1, 0)\n\n\nclass HubDataset(D.Dataset):\n    def __init__(self, path, mode):\n        self.csv = pd.read_csv(path, index_col=[0])\n        self.mode = mode\n        self.x, self.y = [], []\n        self.build_slices()\n\n    def build_slices(self):\n        self.masks = []\n        self.files = []\n        self.slices = []\n        column_name = 'encoding' if self.mode == 'train' else 'predicted'\n        for i, file_id in enumerate(self.csv.index.values):\n            filepath = f'{BASE_DIR}/{self.mode}/{file_id}.tiff'\n            self.files.append(filepath)\n            print('Transform', filepath)\n            with rasterio.open(filepath,\n                               transform=identity, driver='GTiff') as dataset:\n                self.masks.append(\n                    rle2mask(\n                        self.csv.loc[file_id, column_name],\n                        (dataset.shape[1], dataset.shape[0])\n                    )\n                )\n                slices = make_grid(dataset.shape,\n                                   window=WINDOW,\n                                   min_overlap=OVERLAP)\n                for (x1, x2, y1, y2) in tqdm(slices):\n                    use_slice = self.masks[-1][x1:x2, y1:y2].sum() > THRESHOLD\n                    if use_slice:\n                        self.slices.append([i, x1, x2, y1, y2])\n                        image = dataset.read(\n                            [1, 2, 3],  # (r, g, b), h, w\n                            window=Window.from_slices((x1, x2), (y1, y2)))\n                        image = np.moveaxis(image, 0, -1)  # h, w, (r, g, b)\n                        self.x.append(image)\n\n                        mask = self.masks[-1][x1:x2, y1:y2]\n                        mask = np.expand_dims(mask, axis=2)\n                        self.y.append(mask)\n\n                # visualise few patchs only\n                del slices, image, mask, dataset, column_name, use_slice\n                gc.collect()\n                break\n\n    # get data operation\n    def __getitem__(self, index):\n        preproc = get_preprocessing()(\n            image=self.x[index],\n            mask=self.y[index])  # (r, g, b), h, w\n        image, mask = preproc['image'], preproc['mask']\n        return image, mask\n\n    def __len__(self):\n        return len(self.x)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"def stats_blobs(mask, connectivity):\n    output = cv2.connectedComponentsWithStats(mask, connectivity)\n    num_glomerulus = output[0] - 1\n    print(\"Number of glomerulus : {}\".format(num_glomerulus))\n    stats = output[2]\n\n    # areas\n    largest_area = np.max(stats[1:, cv2.CC_STAT_AREA])\n    print(\"Largest area of a glomerulus: {}\".format(largest_area))\n    smallest_area = np.min(stats[1:, cv2.CC_STAT_AREA])\n    print(\"Smallest area of a glomerulus: {}\".format(smallest_area))\n    mean_area = np.mean(stats[1:, cv2.CC_STAT_AREA])\n    print(\"Mean area of a glomerulus: {}\".format(mean_area))\n    covered_area = 100 * len(\n        mask[np.nonzero(mask)]) / (mask.shape[0]*mask.shape[1])\n    print(\"Covered area by glomerulus: {}%\\n\".format(covered_area))\n\n    del output\n    gc.collect()\n\n    return stats","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# A. Train dataset "},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"df_train = pd.read_csv(os.path.join(BASE_DIR, 'train.csv'))\nprint(\"The train dataset has {} images\".format(df_train.shape[0]))\ndf_train = df_train.sort_values(by=['id'])\ndf_train","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## A.1. Visualisation"},{"metadata":{},"cell_type":"markdown","source":"Let's visualize full images of the train dataset on the first time"},{"metadata":{"trusted":true,"scrolled":true},"cell_type":"code","source":"for idx in range(df_train.shape[0]):\n    item = df_train.iloc[idx]\n    plot_image(item, TRAIN_DIR, scale=20, verbose=1)\n    del item\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## A.2. Blob analysis"},{"metadata":{},"cell_type":"markdown","source":"Let's now analyze glomerulus on the train dataset to identify some features : \n- size\n- area\n- distribution\n- numbers"},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"plt.figure(figsize=(10, 7))\nbboxes = []\nfor idx in range(df_train.shape[0]):\n    item = df_train.iloc[idx]\n    print(\"Image id: {}\".format(item.get('id')))\n    tmp_mask = read_mask_only(item, TRAIN_DIR, verbose=0)\n    stats = stats_blobs(tmp_mask, 4)\n\n    # plot histogram\n    plt.hist(stats[1:, cv2.CC_STAT_AREA],\n             histtype='step',\n             label=item.get('id'))\n\n    # get bboxes\n    if idx == 0:\n        bboxes.append([\n            stats[1:30, cv2.CC_STAT_LEFT],\n            stats[1:30, cv2.CC_STAT_TOP],\n            stats[1:30, cv2.CC_STAT_LEFT] + stats[1:30, cv2.CC_STAT_WIDTH],\n            stats[1:30, cv2.CC_STAT_TOP] + stats[1:30, cv2.CC_STAT_HEIGHT]])\n\n    del tmp_mask, item, stats\n    gc.collect()\n\nplt.xlabel('Area (pixels²)')\nplt.ylabel('Frequency')\nplt.legend()\nplt.savefig('train_glomerulus_area_hist.png', bbox_inches='tight')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"These feature allow us to define a realistic range of area for glomerulus and therefore, delete too small glomerulus."},{"metadata":{},"cell_type":"markdown","source":"### Zoom on glomerulus"},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"image, mask = read_image(df_train.iloc[0], TRAIN_DIR, verbose=0)\nbbox = bboxes[0]\n\nfig, axes = plt.subplots(2, 5, figsize=(20, 8))\nval = 10\nfor i in range(2):\n    for j in range(5):\n        axes[i, j].imshow(\n            image[bbox[1][val]:bbox[3][val],\n                  bbox[0][val]:bbox[2][val]])\n        axes[i, j].imshow(\n            mask[bbox[1][val]:bbox[3][val],\n                 bbox[0][val]:bbox[2][val]],\n            alpha=0.3)\n        axes[i, j].axis('off')\n        val += 1\nplt.savefig('train_glomerulus_zoom.png', bbox_inches='tight')\nplt.show()\n\ndel image, mask, bbox, bboxes, df_train\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Glomerulus in the train patches "},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"ds = HubDataset('../input/hubmap-kidney-segmentation/train.csv', 'train')\n\nplt.figure(figsize=(30, 10))\nidx = 1\nfor i in range(30, 42):\n    image, mask = ds[i]\n\n    plt.subplot(3, 8, idx)\n    plt.imshow(image[0])\n    idx += 1\n\n    plt.subplot(3, 8, idx)\n    plt.imshow(image[0])\n    plt.imshow(mask[0], cmap='hot', alpha=0.5)\n    idx += 1\n\n    del image, mask\ndel ds\ngc.collect()\nplt.savefig('train_patches.png', bbox_inches='tight')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# B. Test dataset "},{"metadata":{},"cell_type":"markdown","source":"We can do the same identification on the test dataset and compare our results"},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"df_test = pd.read_csv(os.path.join(TEST_PRED_DIR, 'submission.csv'))\nprint(\"The test dataset has {} images\".format(df_test.shape[0]))\ndf_test = df_test.sort_values(by=['id'])\ndf_test","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## B.1. Visualisation "},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"for idx in range(df_test.shape[0]):\n    item = df_test.iloc[idx]\n    plot_image(item, TEST_DIR, scale=20, verbose=0)\n    del item\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## B.2 Blob analysis"},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"bboxes = []\nfor idx in range(df_test.shape[0]):\n    item = df_test.iloc[idx]\n    print(\"Image id: {}\".format(item.get('id')))\n    tmp_mask = read_mask_only(item, TEST_DIR, verbose=0)\n    stats = stats_blobs(tmp_mask, 4)\n\n    # plot histogram\n    plt.hist(stats[1:, cv2.CC_STAT_AREA],\n             histtype='step',\n             label=item.get('id'))\n\n    # get bboxes\n    if idx == 0:\n        bboxes.append([\n            stats[1:30, cv2.CC_STAT_LEFT],\n            stats[1:30, cv2.CC_STAT_TOP],\n            stats[1:30, cv2.CC_STAT_LEFT] + stats[1:30, cv2.CC_STAT_WIDTH],\n            stats[1:30, cv2.CC_STAT_TOP] + stats[1:30, cv2.CC_STAT_HEIGHT]])\n\n    del tmp_mask, item, stats\n    gc.collect()\nplt.xlabel('Area (pixels²)')\nplt.ylabel('Frequency')\nplt.legend()\nplt.savefig('test_glomerulus_area_hist.png', bbox_inches='tight')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that we have too smal glomerulus in our predictions, we should do some erosion/dilatation to remove these glomerulus and have a best score. I didn't done that because of a lack of memory"},{"metadata":{},"cell_type":"markdown","source":"### Zoom on glomerulus "},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"image, mask = read_image(df_test.iloc[0], TEST_DIR, verbose=0)\nbbox = bboxes[0]\n\nfig, axes = plt.subplots(2, 5, figsize=(20, 8))\nval = 10\nfor i in range(2):\n    for j in range(5):\n        axes[i, j].imshow(\n            image[bbox[1][val]:bbox[3][val],\n                  bbox[0][val]:bbox[2][val]])\n        axes[i, j].imshow(\n            mask[bbox[1][val]:bbox[3][val],\n                 bbox[0][val]:bbox[2][val]],\n            alpha=0.3)\n        axes[i, j].axis('off')\n        val += 1\n\ndel image, mask\ngc.collect()\nplt.savefig('test_glomerulus_zoom.png', bbox_inches='tight')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Glomerulus on test patches "},{"metadata":{"trusted":true},"cell_type":"code","source":"ds = HubDataset('../input/submission-predictions/submission.csv', 'test')\n\nplt.figure(figsize=(30, 10))\nidx = 1\nfor i in range(30, 42):\n    image, mask = ds[i]\n\n    plt.subplot(3, 8, idx)\n    plt.imshow(image[0])\n    idx += 1\n\n    plt.subplot(3, 8, idx)\n    plt.imshow(image[0])\n    plt.imshow(mask[0], cmap='hot', alpha=0.5)\n    idx += 1\n\n    del image, mask\ndel ds\ngc.collect()\nplt.savefig('test_patches.png', bbox_inches='tight')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## B.3. Influence of probability threshold on predictions"},{"metadata":{},"cell_type":"markdown","source":"During my research of the best model for this competition, I've tried several threshold of binarization on my probability mask. You can visualisze below the effect of these thresholds on the final mask.\nMy best results were with threshold=0.5 as on the pictures."},{"metadata":{"trusted":true,"scrolled":false},"cell_type":"code","source":"p_files = pathlib.Path('../input/tresholds/')\nthresholds = ['0.2', '0.5', '0.7']\ncols = 4\nrows = 5\nimage, mask = read_image(df_test.iloc[0],\n                         TEST_DIR,\n                         verbose=0,\n                         with_mask=False)\n\nfig, axes = plt.subplots(rows, cols, figsize=(20, 15))\nfor j, filename in enumerate(p_files.glob('*.csv')):\n    print(filename)\n    val = 10\n    df_temp = pd.read_csv(filename)\n    df_temp = df_temp.sort_values(by=['id'])\n    mask = read_mask_only(df_temp.iloc[0], TEST_DIR, verbose=0)\n    for i in range(5):\n        if j == 0:\n            axes[i, j].imshow(\n                image[bbox[1][val] - 50:bbox[3][val] + 50,\n                      bbox[0][val] - 50:bbox[2][val] + 50])\n            axes[i, j].axis('off')\n            axes[i, j].set_title(f'original image')\n\n        axes[i, j+1].imshow(\n            image[bbox[1][val] - 50:bbox[3][val] + 50,\n                  bbox[0][val] - 50:bbox[2][val] + 50])\n        axes[i, j+1].imshow(\n            mask[bbox[1][val] - 50:bbox[3][val] + 50,\n                 bbox[0][val] - 50:bbox[2][val] + 50],\n            alpha=0.3)\n        axes[i, j+1].axis('off')\n        axes[i, j+1].set_title(\n            f'image + mask with th={thresholds[j]}')\n        val += 1\n\n    del mask, df_temp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del bboxes, image, bbox\ngc.collect()","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}