{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.9","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":22990,"datasetId":1136396,"databundleVersionId":2048213},{"sourceType":"kernelVersion","sourceId":56694565}],"dockerImageVersionId":30066,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# HuBMAP masks and probability density function\n\n> Idea: Knowing that the glomeruli are mainly found in the cortex, we should focus on this region during training.","metadata":{}},{"cell_type":"markdown","source":"## Background\n\nA glomerulus is a network of small blood vessels located at the beginning of a nephron in the kidney ([Wikipedia](https://en.wikipedia.org/wiki/Glomerulus_(kidney))). Glomeruli are mainly found in the renal **cortex**, while the renal **medulla** contains mainly the renal tubule. Since we are dealing with biological structures, the separation is not not absolute and the transitions are not always perfectly sharp.\n\n![Diagram of a nephron](http://s3-us-west-2.amazonaws.com/courses-images/wp-content/uploads/sites/1842/2017/05/26234530/m9skcbftjqzrokkkopam.png)\n[Diagram of a nephron from libretexts.org, Introductory and General Biology](https://bio.libretexts.org/Bookshelves/Introductory_and_General_Biology/Book%3A_General_Biology_(Boundless)/41%3A_Osmotic_Regulation_and_the_Excretory_System/41.4%3A_Human_Osmoregulatory_and_Excretory_Systems/41.4B%3A_Nephron%3A_The_Functional_Unit_of_the_Kidney)\n\n## Key Idea\n\nA common approach to deal with the very large (>500MB - 5GB) TIFF files in the dataset is to decompose the images in smaller patches/tiles, for instance by using a sliding window apporach.\n> **Knowing that the glomeruli are mainly found in the cortex, we should focus on this region during training**. \n\nInstead of preprocessing the images by saving them into fixed tiles, we prepare the data to sample random tiles from the entire images with a higher probability to sample tiles that contain glomeruli and cortex.\n\n**We combine two sampling apporaches**\n\n> 1. Sampling tiles that contain all glomeruli - this ensures that each glomerulus is seen during one epoch of training at least once. \n> 2. Sampling random tiles based on region probabilites (e.g., medulla; cortex)\n\n*Note: We startet out by only using the second approach (see our public kernel https://www.kaggle.com/matjes/hubmap-labels-pdf-0-5-0-25-0-01). However, we realized that due to the randomness of the sampling sometimes not all glomeruli were drawn during training. Thus, we added the first step.*\n\n## Advantages of this approaches\n\nIn combination with [deepflash2](https://github.com/matjesg/deepflash2/tree/master/) and the deepflash2 [pytorch datasets](https://matjesg.github.io/deepflash2/data.html#Datasets) in particular, this approach has several advantages:\n- no preprocessing of the data (only saving them to .zarr files for memory efficient loading)\n    - flexible tile shapes (input shapes, e.g. 1024, 512, 256) at runtime\n    - flexible scaling (e.g., by facors of 2,3,4)\n- faster convergence during traing\n    - focusing on the relevant regions (e.g., tiles that contain glomeruli and cortex)\n    - \"additional\" data augmentation from random sampling (compared to fixed windows)","metadata":{}},{"cell_type":"markdown","source":"## Implementation (masks and probability density function)\n\nThe data preprocessing steps in this notebook comprise\n- Creating segmentations masks from RLE\n- Creating a data frame containing position information on each glomerulus in the segmentation mask\n- Creating probability density function (PDF) for efficient sampling from mask and anatmical structure\n\n**Inputs**\n- RLE glomeruli segmentations ([challenge dataset](https://www.kaggle.com/c/hubmap-kidney-segmentation/data))\n- Anatomical region masks ([.zarr converted](https://www.kaggle.com/matjes/hubmap-pdf-zarr))\n   \n**Settings**\n\n- Probability for cortex regions (`cortex_p`): 0.7 \n- Probability for medulla regions (`medulla_p`): 0.2\n- Probability for other regions (`bg_p` (background)): 0.1    \n\nWe chose the region sampling probabilites based on expert knowledge and some experiements. However, these settings are hyperparamteres that can be tuned or adjusted if the unterlying data changes.\n","metadata":{}},{"cell_type":"code","source":"# Install zarr and load packages\n!pip install -qq zarr\nimport cv2, zarr, gc\nimport matplotlib.pyplot as plt, numpy as np, pandas as pd\nfrom pathlib import Path\ngc.enable()","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from https://www.kaggle.com/paulorzp/rle-functions-run-lenght-encode-decode\ndef rle2mask(mask_rle, shape):\n    '''\n    mask_rle: run-length as string formated (start length)\n    shape: (width,height) of array to return \n    Returns numpy array, 1 - mask, 0 - background\n\n    '''\n    s = mask_rle.split()\n    starts, lengths = [\n        np.asarray(x, dtype=int) for x in (s[0:][::2], s[1:][::2])\n    ]\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    return img.reshape(shape).T","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Settings","metadata":{}},{"cell_type":"code","source":"class CONFIG:\n    path = Path('/kaggle/input/hubmap-kidney-segmentation')\n    anatomy = '../input/hubmap-anatomy-zarr/anatomy'\n    \n    scale = 2        # Downscale final mask by factor 2\n    cdf_size = 512   # Downscale CDF for memory efficient loading during training\n    bg_p = 0.1       # Background Probability\n    cortex_p = 0.7   # Cortex Probability\n    medulla_p = 0.2  # Medulla Probability\n    \ncfg = CONFIG()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load inputs and define outputs**","metadata":{}},{"cell_type":"code","source":"# Input \ndf_train = pd.read_csv(cfg.path/\"train.csv\",index_col='id')\ndf_info = pd.read_csv(cfg.path/\"HuBMAP-20-dataset_information.csv\")\ngrp_pdf = zarr.open_group(cfg.anatomy) # Anatomy labels sever as basis for PDF\n\n# Output\nroot = zarr.group(f'/kaggle/working/masks_scale{cfg.scale}')\n# Saving cdf in 'pdfs' due to naming conventions for sampling during training in deepflash2\ng_msk, g_pdf, g_cdf = root.create_groups('labels', 'pdfs', 'cdfs', overwrite=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Loop over files to create...**\n1. The segmentation mask (.zarr)\n1. A list of glomeruli for *glomeruli sampling*\n1. The probability density function for *region sampling*","metadata":{}},{"cell_type":"code","source":"df_list = []\nfor idx, row in df_train.iterrows():\n    \n    # Get image info\n    img_info = df_info[df_info.image_file==f'{idx}.tiff']\n    shape = (img_info.height_pixels.values[0], img_info.width_pixels.values[0])\n    print(idx, 'with shape', shape)\n    \n    msk = rle2mask(row.encoding, (shape[1], shape[0])).astype('uint8')\n    \n    # Plot\n    fig, ax = plt.subplots(ncols=2, figsize=(15,15))\n    resize_w = int((msk.shape[1]/msk.shape[0])*cfg.cdf_size)\n    ax[0].imshow(cv2.resize(msk, dsize=(resize_w, cfg.cdf_size)))\n    ax[0].set_title('Mask')\n    ax[0].set_axis_off()\n    \n    pdf = grp_pdf[idx][:]\n    \n    if cfg.scale!=1:\n        new_size = (msk.shape[1] // cfg.scale, msk.shape[0] // cfg.scale)\n        print('Scaling to', new_size)\n        msk = cv2.resize(msk, new_size)\n        pdf = cv2.resize(pdf, new_size)\n        \n    pdf = pdf.astype('float32')          \n    pdf[pdf==0] = cfg.bg_p/np.sum(pdf==0)\n    pdf[msk>0] = 0\n    pdf[pdf==1] = cfg.cortex_p/np.sum(pdf==1)\n    pdf[pdf==2] = cfg.medulla_p/np.sum(pdf==2)    \n      \n    print('Getting glomeruli stats')\n    nb_components, output, stats, centroids = cv2.connectedComponentsWithStats(msk, connectivity=4)\n    print(f'Found {nb_components} glomeruli')\n    df_centroids = pd.DataFrame(centroids[1:], columns=['cy', 'cx'])\n    df_centroids = df_centroids.join(pd.DataFrame(stats[1:], columns=['left', 'top', 'width', 'height', 'area']))\n    df_centroids['idx'] = idx \n    df_centroids.reset_index(inplace=True)\n    df_centroids.set_index(['idx', 'index'], inplace=True)\n    df_list.append(df_centroids)\n    \n    # Saving \n    g_msk[idx] = msk\n    g_pdf[idx] = pdf\n    \n    # Saving cdf\n    pdf = cv2.resize(pdf, dsize=(resize_w, cfg.cdf_size))      \n    g_cdf[idx] = np.cumsum(pdf/np.sum(pdf)) \n            \n    ax[1].imshow(pdf)\n    ax[1].set_title('Probability density function for sampling')\n    ax[1].set_axis_off() \n    plt.show()\n\n\ndf_stats = pd.concat(df_list)\ndf_stats.to_csv(f'/kaggle/working/masks_scale{cfg.scale}/roi_stats.csv')\ndf_stats","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}