{"cells":[{"metadata":{},"cell_type":"markdown","source":"This kernel is inspired from [this kernel](http://https://www.kaggle.com/dannellyz/tissue-detect-scaling-bounding-boxes-4xfaster) . Thanks Zac Dannelly for your amazing kernel. Please upvote his kernel. Code for making dataset is taken from [iafosss's kernel](http://https://www.kaggle.com/iafoss/panda-16x128x128-tiles). Note that I haven't used tiles.\n\nChanges:\n1. As mentioned in the comments of Zac's kernel, resizing images can lead to deformations which are undesired. \n2. So, here I dont resize the images and save the images in the same shape as they are.\n3. I'll try to use Spatial Pyramid Pooling layer to account for different shape inputs during training. Thanks to Nanashi for telling about SPP. Here is the SPP paper: https://arxiv.org/pdf/1406.4729.pdf","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport cv2\nimport skimage.io\nfrom tqdm.notebook import tqdm\nimport zipfile\nimport numpy as np\n# All imports\nimport time\n\nimport matplotlib.pyplot as plt\nimport openslide\nimport pandas as pd\nfrom skimage import morphology\nimport pprint\nimport json","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"TRAIN = '../input/prostate-cancer-grade-assessment/train_images/'\nOUT_TRAIN = 'train.zip'\n\n# Set up example slide and run pipeline on low resolution\nslide_dir = \"../input/prostate-cancer-grade-assessment/train_images/\"\nannotation_dir = \"../input/prostate-cancer-grade-assessment/train_label_masks/\"\nexample_id = \"0032bfa835ce0f43a92ae0bbab6871cb\"\nexample_slide = f\"{slide_dir}{example_id}.tiff\"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"WonderFul code by Zac Dannelly. I removed the resize part and adjusted the code to work.","execution_count":null},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def get_disk_size(numpy_image):\n    \"\"\"Return disk size of a numpy array\"\"\"\n    return (numpy_image.size * numpy_image.itemsize) / 1000000\n\n\ndef detect_tissue_external(input_slide, sensitivity=3000):\n    \n    \"\"\"\n    Description\n    ----------\n    Find RoIs containing tissue in WSI and only return the external most.\n    Generate mask locating tissue in an WSI. Inspired by method used by\n    Wang et al. [1]_.\n    .. [1] Dayong Wang, Aditya Khosla, Rishab Gargeya, Humayun Irshad, Andrew\n    H. Beck, \"Deep Learning for Identifying Metastatic Breast Cancer\",\n    arXiv:1606.05718\n    Credit: Github-wsipre\n    \n    Parameters\n    ----------\n    input_slide: numpy array\n        Slide to detect tissue on.\n    sensitivity: int\n        The desired sensitivty of the model to detect tissue. The baseline is set\n        at 3000 and should be adjusted down to capture more potential issue and\n        adjusted up to be more agressive with trimming the slide.\n        \n    Returns (3)\n    -------\n    -Tissue binary mask as numpy 2D array, \n    -Tiers investigated,\n    -Time Stamps from running tissue detection pipeline\n    \"\"\"\n    \n    # For timing\n    time_stamps = {}\n    time_stamps[\"start\"] = time.time()\n\n    # Convert from RGB to HSV color space\n    slide_hsv = cv2.cvtColor(input_slide, cv2.COLOR_BGR2HSV)\n    time_stamps[\"re-color\"] = time.time()\n    # Compute optimal threshold values in each channel using Otsu algorithm\n    _, saturation, _ = np.split(slide_hsv, 3, axis=2)\n\n    mask = otsu_filter(saturation, gaussian_blur=True)\n    time_stamps[\"filter\"] = time.time()\n    # Make mask boolean\n    mask = mask != 0\n\n    mask = morphology.remove_small_holes(mask, area_threshold=sensitivity)\n    mask = morphology.remove_small_objects(mask, min_size=sensitivity)\n    time_stamps[\"morph\"] = time.time()\n    mask = mask.astype(np.uint8)\n    mask_contours, tier = cv2.findContours(mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n    time_stamps[\"contour\"] = time.time()\n    time_stamps = {\n        key: (value - time_stamps[\"start\"]) * 1000 for key, value in time_stamps.items()\n    }\n    \n    return mask_contours, tier, time_stamps\n\n\ndef new_detect_and_crop(image_location=\"\",sensitivity: int = 3000, downsample_lvl = -1,\n                        show_plots= \"simple\", out_lvl=-2):\n    \"\"\"\n    Description\n    ----------\n    This method performs the pipeline as described in the notebook:\n    https://www.kaggle.com/dannellyz/panda-tissue-detect-scaling-bounding-boxes-fast\n    \n    Parameters\n    ----------\n    image_location:str\n        Location of the slide image to process\n    sensitivity:int\n        The desired sensitivty of the model to detect tissue. The baseline is set\n        at 3000 and should be adjusted down to capture more potential issue and\n        adjusted up to be more agressive with trimming the slide.\n    downsample_lvl: int\n        The level at which to downsample the slide. This can be referenced in\n        reverse order to access the lowest resoltuion items first.\n        [-1] = lowest resolution\n        [0] = highest resolution\n    show_plots: str (verbose|simple|none)\n        The types of plots to display:\n            - verbose - show all steps of process\n            - simple - show only last step\n            - none - show none of the plots\n    out_lvl: int\n        The level at which the final slide should sample at. This can be referenced in\n        reverse order to access the lowest resoltuion items first.\n        [-1] = lowest resolution\n        [0] = highest resolution\n    shape: touple\n        (height, width) of the desired produciton(prod) image\n        \n    Returns (4)\n    -------\n    - Numpy array of final produciton(prod) slide\n    - Percent memory reduciton from original slide\n    - Time stamps from stages of the pipeline\n    - Time stamps from the Tissue Detect pipeline\n    \"\"\"\n    # For timing\n    time_stamps = {}\n    time_stamps[\"start\"] = time.time()\n\n    # Open Small Slide\n    wsi_small = skimage.io.MultiImage(image_location)[downsample_lvl]\n    time_stamps[\"open_small\"] = time.time()\n\n    # Get returns from detect_tissue() ons mall image\n    (   tissue_contours,\n        tier,\n        time_stamps_detect,\n    ) = detect_tissue_external(wsi_small, sensitivity)\n    \n    base_slide_mask = np.zeros(wsi_small.shape[:2])\n    # Get minimal bounding rectangle for all tissue contours\n    if len(tissue_contours) == 0:\n        img_id = image_location.split(\"/\")[-1]\n        print(f\"No Tissue Contours - ID: {img_id}\")\n        return None, 0, None, None\n    \n    # Open Big Slide\n    wsi_big = skimage.io.MultiImage(image_location)[out_lvl]\n    time_stamps[\"open_big\"] = time.time()\n    \n    #Get small boudning rect and scale\n    bounding_rect_small = cv2.minAreaRect(np.concatenate(tissue_contours))\n\n    # Scale Rectagle to larger image\n    scale = int(wsi_big.shape[0] / wsi_small.shape[0])\n    scaled_rect = (\n        (bounding_rect_small[0][0] * scale, bounding_rect_small[0][1] * scale),\n        (bounding_rect_small[1][0] * scale, bounding_rect_small[1][1] * scale),\n        bounding_rect_small[2],\n    )\n    # Crop bigger image with getSubImage()\n    scaled_crop = getSubImage(wsi_big, scaled_rect)\n    time_stamps[\"scale_bounding\"] = time.time()\n    \n    #Cut out white\n    white_cut = color_cut(scaled_crop)\n    time_stamps[\"white_cut_big\"] = time.time()\n    \n\n    # Get returns from detect_tissue() on small image\n    (   tissue_contours_big,\n        tier_big,\n        time_stamps_detect,\n    ) = detect_tissue_external(white_cut, sensitivity)\n    prod_slide = tissue_cutout(white_cut, tissue_contours_big)\n    time_stamps[\"remove_tissue\"] = time.time()\n\n    # Get size change\n    base_size_high = get_disk_size(wsi_big)\n    final_size = get_disk_size(prod_slide)\n    pct_change = final_size / base_size_high\n    \n    if show_plots == \"simple\":\n        print(f\"Percent Reduced from Base Slide to Final: {(1- pct_change)*100:.2f}\")\n        plt.imshow(smart_bounding_crop)\n        plt.show()\n    elif show_plots == \"verbose\":\n        # Set-up dictionary for plotting\n        verbose_plots = {}\n        # Add Base Slide to verbose print\n        verbose_plots[f\"Smaller Slide\\n{get_disk_size(wsi_small):.2f}MB\"] = wsi_small\n        # Add Tissue Only to verbose print\n        verbose_plots[f\"Tissue Detect Low\\nNo Change\"] = wsi_big\n        # Add Larger Plot cut with bounding boxes\n        verbose_plots[f\"Larger scaled\\n{get_disk_size(scaled_crop):.2f}MB\"] = scaled_crop\n        # Add Bounding Boxes to verbose print\n        verbose_plots[\n            f\"Final Produciton\\n{get_disk_size(prod_slide):.2f}MB\"\n        ] = prod_slide\n        print(f\"Percent Reduced from Base Slide to Final: {(1- pct_change)*100:.2f}\")\n        plt = plot_figures(verbose_plots, 2, 2)\n    elif show_plots == \"none\":\n        pass\n    else:\n        pass\n    time_stamps = {\n        key: (value - time_stamps[\"start\"]) * 1000 for key, value in time_stamps.items()\n    }\n    return prod_slide, (1 - pct_change), time_stamps, time_stamps_detect\n\ndef otsu_filter(channel, gaussian_blur=True):\n    \n    \"\"\"Otsu filter.\"\"\"\n    \n    if gaussian_blur:\n        channel = cv2.GaussianBlur(channel, (5, 5), 0)\n    channel = channel.reshape((channel.shape[0], channel.shape[1]))\n\n    return cv2.threshold(channel, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)[1]\n\ndef getSubImage(input_slide, rect):\n    \n    \"\"\"\n    Description\n    ----------\n    Take a cv2 rectagle object and remove its contents from\n    a source image.\n    Credit: https://stackoverflow.com/a/48553593\n    \n    Parameters\n    ----------\n    input_slide: numpy array \n            Slide to pull subimage off \n    rect: cv2 rect\n        cv2 rectagle object with a shape of-\n            ((center_x,center_y), (hight,width), angle)\n    \n    Returns (1)\n    -------\n    - Numpy array of rectalge data cut from input slide\n    \"\"\"\n    \n    width = int(rect[1][0])\n    height = int(rect[1][1])\n    box = cv2.boxPoints(rect)\n\n    src_pts = box.astype(\"float32\")\n    dst_pts = np.array(\n        [[0, height - 1], [0, 0], [width - 1, 0], [width - 1, height - 1]],\n        dtype=\"float32\",\n    )\n    M = cv2.getPerspectiveTransform(src_pts, dst_pts)\n    output_slide = cv2.warpPerspective(input_slide, M, (width, height))\n    return output_slide\n\ndef color_cut(in_slide, color = [255,255,255]):\n    \n    \"\"\"\n    Description\n    ----------\n    Take a input image and remove all rows or columns that\n    are only made of the input color [R,G,B]. The default color\n    to cut from image is white.\n    \n    Parameters\n    ----------\n    input_slide: numpy array \n        Slide to cut white cols/rows \n    color: list\n        List of [R,G,B] pixels to cut from the input slide\n    \n    Returns (1)\n    -------\n    - Numpy array of input_slide with white removed\n    \"\"\"\n    #Remove by row\n    row_not_blank = [row.all() for row in ~np.all(in_slide == color, axis=1)]\n    output_slide = in_slide[row_not_blank, :]\n    \n    #Remove by col\n    col_not_blank = [col.all() for col in ~np.all(output_slide == color, axis=0)]\n    output_slide = output_slide[:, col_not_blank]\n    return output_slide\n\ndef tissue_cutout(input_slide, tissue_contours):\n    \n    \"\"\"\n    Description\n    ----------\n    Set all parts of the in_slide to black except for those\n    within the provided tissue contours\n    Credit: https://stackoverflow.com/a/28759496\n    \n    Parameters\n    ----------\n    input_slide: numpy array\n            Slide to cut non-tissue backgound out\n    tissue_contours: numpy array \n            These are the identified tissue regions as cv2 contours\n            \n    Returns (1)\n    -------\n    - Numpy array of slide with non-tissue set to black\n    \"\"\"\n    \n    # Get intermediate slide\n    base_slide_mask = np.zeros(input_slide.shape[:2])\n    # Create mask where white is what we want, black otherwise\n    crop_mask = np.zeros_like(base_slide_mask) \n    \n    # Draw filled contour in mask\n    cv2.drawContours(crop_mask, tissue_contours, -1, 255, -1) \n    \n    # Extract out the object and place into output image\n    tissue_only_slide = np.zeros_like(input_slide)  \n    tissue_only_slide[crop_mask == 255] = input_slide[crop_mask == 255]\n    \n    return tissue_only_slide\n\ndef plot_figures(figures, nrows=1, ncols=1):\n    \n    \"\"\"\n    Description\n    ----------\n    Plot a dictionary of figures.\n    Credit: https://stackoverflow.com/a/11172032\n\n    Parameters\n    ----------\n    figures: dict \n        <title, figure> for those to plot\n    ncols: int \n        number of columns of subplots wanted in the display\n    nrows: int \n        number of rows of subplots wanted in the figure\n    \n    Returns(0)\n    ----------\n    \"\"\"\n\n    fig, axeslist = plt.subplots(ncols=ncols, nrows=nrows)\n    for ind, title in enumerate(figures):\n        axeslist.ravel()[ind].imshow(figures[title], aspect=\"auto\")\n        axeslist.ravel()[ind].set_title(title)\n    plt.tight_layout()\n    plt.show()\n    return ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Processed on level 1 image (Intermediate Resolution)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"(processed_slide_med, pct_change_med, \n time_stamps_pipeline_med, detect_time_med) = new_detect_and_crop(\n                                        image_location=example_slide, out_lvl=-2,\n                                        show_plots=\"verbose\"\n                                        )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Processed on level 2 image (Lowest Resolution)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"(processed_slide_med, pct_change_med, \n time_stamps_pipeline_med, detect_time_med) = new_detect_and_crop(\n                                        image_location=example_slide, out_lvl=-1,\n                                        show_plots=\"verbose\"\n                                        )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Making Dataset","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"x_tot,x2_tot = [],[]\nnames = [name[:-5] for name in os.listdir(TRAIN)]\nwith zipfile.ZipFile(OUT_TRAIN, 'w') as img_out:\n    for name in tqdm(names):\n        img_path = f\"{slide_dir}{name}.tiff\"\n        img,_,_,_ = new_detect_and_crop(image_location=img_path, show_plots=None, downsample_lvl=-1, out_lvl=-1)\n        if img is None:\n            continue\n        x_tot.append((img/255.0).reshape(-1,3).mean(0))\n        x2_tot.append(((img/255.0)**2).reshape(-1,3).mean(0)) \n         #if read with PIL RGB turns into BGR\n        img = cv2.imencode('.png',cv2.cvtColor(img, cv2.COLOR_RGB2BGR))[1]\n        img_out.writestr(f'{name}.png', img)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Statistics","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"img_avr =  np.array(x_tot).mean(0)\nimg_std =  np.sqrt(np.array(x2_tot).mean(0) - img_avr**2)\nprint('mean:',img_avr, ', std:', np.sqrt(img_std))","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}