Session 1 -- Intro to QuPath and OpenSlide.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "17dc3ea1",
   "metadata": {},
   "source": [
    "# Understanding whole-slide images:\n",
    "\n",
    "This first hands-on session will guide you to learn how to visualize and manipulate giga-pixel whole-slide images (WSIs).\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "735431cc",
   "metadata": {},
   "source": [
    "## Qupath  \n",
    "\n",
    "![](https://github.com/guillaumejaume/mil-tutorial/blob/main/support/wsi_example.png?raw=true)\n",
    "\n",
    "We will be using QuPath, an open-source software that allows visualizing, annotating bounding boxes, and extracting information from WSIs.\n",
    "\n",
    "When building a computational pathology pipeline, you will typically need to:\n",
    "- Assess the image quality: ensure staining quality, images are not blurred, or have ink stains\n",
    "- Validate the quality of the annotations, either at slide-level or regions (boxes or dense pixel annotations)\n",
    "- understand the type of task/prior that needs to be considered when learning on your images, e.g.,\n",
    "   - if classifying different tumor subtypes, what is the percentage of the slide that occupies the tumor \n",
    "   - if detecting small objects, such as mitosis, what is the appropriate resolution\n",
    "- Validate the quality of the Deep Learning system on fined-grained annotations\n",
    "\n",
    "QuPath is the best tool for all these tasks! \n",
    "\n",
    "Instructions:\n",
    "\n",
    "- Download QuPath (https://qupath.readthedocs.io/en/0.4/docs/intro/installation.html#download-install)\n",
    "   - Support Windows, Mac and Ubuntu \n",
    "   - Requires 4GB of RAM (runs smoothly on most modern computers)\n",
    "\n",
    "Download WSIs from here:\n",
    "\n",
    "- https://drive.google.com/file/d/1GWtuOLRf6-7C8zCG6b3JdYDPk_C5IS3_/view?usp=sharing\n",
    "- https://drive.google.com/file/d/1o9Es1qdY15y8q_gsjX3KVkn7HS8amlSF/view?usp=sharing"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "aa3d9b41",
   "metadata": {},
   "source": [
    "Using QuPath, play around with the two WSIs:\n",
    "- a necropsy of a rat liver with presence of a necrotic region\n",
    "- a TCGA sample of invasive breast cancer (BRCA) "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "89911a45",
   "metadata": {},
   "source": [
    "## OpenSlide \n",
    "\n",
    "QuPath is helpful for a quick exploration of a new dataset, make small annotations etc, but it cannot be used in a robust deep learning pipeline. Instead, we can use OpenSlide, an open-source Python package that can open, read, and save WSIs (or typically a small region of the slide). \n",
    "\n",
    "OpenSlide is a very small package with a fairly small and neat API (see https://openslide.org/api/python/). It also supports most pyramidal formats, such as tiff, svs, ndpi, mrxs, etc. These formats tyically depends on the scanner used to digitize the slides. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cdcd4742",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: openslide-python in /Users/V/miniconda3/envs/mil/lib/python3.10/site-packages (1.2.0)\n",
      "Requirement already satisfied: Pillow in /Users/V/miniconda3/envs/mil/lib/python3.10/site-packages (from openslide-python) (10.0.0)\n"
     ]
    }
   ],
   "source": [
    "# Install openslide \n",
    "!pip install openslide-python"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "899212b1",
   "metadata": {},
   "source": [
    "### Accessing slide properties and reading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "51712c96",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image size is:  61751  x  35769\n",
      "Number of levels:  3\n",
      "All levels:\n",
      "    - ( 61751  x  35769 )\n",
      "    - ( 15437  x  8942 )\n",
      "    - ( 3859  x  2235 )\n",
      "file size: 671.22 MB\n"
     ]
    }
   ],
   "source": [
    "# Understand slide properties \n",
    "\n",
    "import os\n",
    "from openslide import OpenSlide\n",
    "\n",
    "slide_path = os.path.join('data', 'slides', 'necrosis.tiff')\n",
    "slide = OpenSlide(slide_path)\n",
    "\n",
    "# Check the image dimensions at the highest resolution (20x here): \n",
    "print('Image size is: ', slide.dimensions[0], ' x ', slide.dimensions[1])\n",
    "\n",
    "# Check the number of pre-extracted levels:\n",
    "print('Number of levels: ', slide.level_count)\n",
    "\n",
    "# Check the dimensions of each level:\n",
    "print('All levels:')\n",
    "for dims in slide.level_dimensions:\n",
    "    print('    - (', dims[0], ' x ', dims[1], ')')\n",
    "\n",
    "# Check slide size \n",
    "print(f\"file size: {round(os.path.getsize(slide_path) / 1024 ** 2, 2)} MB\")\n",
    "\n",
    "# Optional: Check all slide properties \n",
    "# pprint(dict(slide.properties))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "64776287",
   "metadata": {},
   "source": [
    "![](https://github.com/guillaumejaume/mil-tutorial/blob/main/support/context_resolution.png?raw=true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d5080f52",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'plt' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[3], line 10\u001b[0m\n\u001b[1;32m      8\u001b[0m region \u001b[39m=\u001b[39m slide\u001b[39m.\u001b[39mread_region(topleft, \u001b[39m0\u001b[39m, size)\n\u001b[1;32m      9\u001b[0m region \u001b[39m=\u001b[39m region\u001b[39m.\u001b[39mresize((\u001b[39m256\u001b[39m, \u001b[39m256\u001b[39m))\n\u001b[0;32m---> 10\u001b[0m plt\u001b[39m.\u001b[39mfigure(figsize\u001b[39m=\u001b[39m(\u001b[39m10\u001b[39m, \u001b[39m10\u001b[39m))\n\u001b[1;32m     11\u001b[0m plt\u001b[39m.\u001b[39mimshow(region)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'plt' is not defined"
     ]
    }
   ],
   "source": [
    "# Extracting and visualizing regions\n",
    "\n",
    "topleft = (25000, 25000)\n",
    "size = (256, 256)\n",
    "\n",
    "for l in range(4):\n",
    "    size = (size[0] * 2, size[1] * 2)\n",
    "    region = slide.read_region(topleft, 0, size)\n",
    "    region = region.resize((256, 256))\n",
    "    plt.figure(figsize=(10, 10))\n",
    "    plt.imshow(region)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "4bbfe6e7",
   "metadata": {},
   "source": [
    "### Tissue segmentation\n",
    "\n",
    "![](https://github.com/guillaumejaume/mil-tutorial/blob/main/support/segmentation.png?raw=true)\n",
    "\n",
    "#### Importance of tissue segmentation: \n",
    "\n",
    "- **Removal of Non-Biological Background**: In histopathology images, the tissue of interest is embedded in a background that often contains artifacts, stains, penmarks or other non-biological elements. These elements do not carry any meaningful biological information and can interfere with the accurate analysis of tissue structures. By segmenting and removing the background, computational pathology algorithms can focus solely on the relevant tissue regions, leading to more precise and reliable results.\n",
    "\n",
    "\n",
    "- **Elimination of Spurious Correlations**: Background regions can introduce spurious correlations in the downstream analysis. If the background is not properly removed, image analysis algorithms may inadvertently assign importance to irrelevant features or patterns present in the background. This can lead to false associations or misleading results, affecting the accuracy and reliability of computational pathology findings. Tissue segmentation helps avoid such spurious correlations by isolating and analyzing only the biologically relevant components.\n",
    "\n",
    "\n",
    "- **Reduction of Computational Complexity**: Background regions often contribute to the overall complexity of image analysis algorithms without providing any meaningful information. Including these regions in computational processes unnecessarily increases the computational burden and time required for analysis. By segmenting and removing the background, the computational complexity can be significantly reduced, allowing for more efficient and streamlined analysis workflows. This optimization is particularly valuable in large-scale studies or applications that require processing a vast number of histopathology images.\n",
    "\n",
    "\n",
    "#### Main techniques for tissue segmentation:\n",
    "\n",
    "- **Thresholding**: This is a straightforward technique where a specific intensity threshold is set to distinguish tissue from the background. Pixels with intensity values above the threshold are classified as tissue, while those below it are considered background. This method assumes a clear contrast between the tissue and background, making it suitable for images with well-defined boundaries.\n",
    "\n",
    "\n",
    "- **Region-based Methods**: These methods utilize image properties such as color, texture, or intensity variations to identify tissue regions. Techniques like region growing, watershed segmentation, or graph cuts analyze neighboring pixel properties to iteratively expand or partition regions until tissue boundaries are determined. Region-based methods are useful when there are variations in tissue appearance or when the tissue and background exhibit different characteristics.\n",
    "\n",
    "\n",
    "- **Machine Learning-based Methods**: With advancements in machine learning, techniques like supervised or unsupervised learning can be used for tissue segmentation. Supervised learning involves training a model on a set of labeled images, where the model learns to distinguish tissue and background regions based on the provided annotations. Unsupervised learning methods, such as clustering or autoencoders, can identify patterns or clusters within the image data, enabling the separation of tissue from the background.\n",
    "\n",
    "\n",
    "- **Deep Learning-based Methods**: Deep learning approaches, particularly convolutional neural networks (CNNs), have shown remarkable success in tissue segmentation tasks. These models can learn hierarchical features directly from the image data and make pixel-level predictions. Fully Convolutional Networks (FCNs), U-Net, or Mask R-CNN are commonly used architectures for tissue segmentation. \n",
    "\n",
    "\n",
    "#### Existing libraries that include tissue segmentations:\n",
    "\n",
    "- CLAM (see https://github.com/mahmoodlab/CLAM)\n",
    "- HistoCartography (see https://github.com/BiomedSciAI/histocartography)\n",
    "- TIA ToolBox (see https://github.com/TissueImageAnalytics/tiatoolbox)\n",
    "\n",
    "An overview of existing tools can be found [here](https://proceedings.mlr.press/v156/jaume21a/jaume21a.pdf). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8f4ee5b7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of contours detected: 2\n",
      "Number of holes detected: 4\n"
     ]
    }
   ],
   "source": [
    "# Extract contours where tissue is present (also removing holes!)\n",
    "\n",
    "import numpy as np\n",
    "import cv2\n",
    "\n",
    "def segment_tissue(\n",
    "    wsi,\n",
    "    seg_level=2,\n",
    "    sthresh=20,\n",
    "    sthresh_up = 255,\n",
    "    mthresh=7,\n",
    "    close=0,\n",
    "    use_otsu=False, \n",
    "    filter_params={'a_t':100, 'a_h': 4, 'max_n_holes': 5},\n",
    "    ref_patch_size=512,\n",
    "    exclude_ids=[],\n",
    "    keep_ids=[]):\n",
    "        \"\"\"\n",
    "            Segment the tissue via HSV -> Median thresholding -> Binary threshold\n",
    "        \"\"\"\n",
    "        \n",
    "        def _filter_contours(contours, hierarchy, filter_params):\n",
    "            \"\"\"\n",
    "                Filter contours by: area.\n",
    "            \"\"\"\n",
    "            filtered = []\n",
    "\n",
    "            # find indices of foreground contours (parent == -1)\n",
    "            hierarchy_1 = np.flatnonzero(hierarchy[:,1] == -1)\n",
    "            all_holes = []\n",
    "            \n",
    "            # loop through foreground contour indices\n",
    "            for cont_idx in hierarchy_1:\n",
    "                # actual contour\n",
    "                cont = contours[cont_idx]\n",
    "                # indices of holes contained in this contour (children of parent contour)\n",
    "                holes = np.flatnonzero(hierarchy[:, 1] == cont_idx)\n",
    "                # take contour area (includes holes)\n",
    "                a = cv2.contourArea(cont)\n",
    "                # calculate the contour area of each hole\n",
    "                hole_areas = [cv2.contourArea(contours[hole_idx]) for hole_idx in holes]\n",
    "                # actual area of foreground contour region\n",
    "                a = a - np.array(hole_areas).sum()\n",
    "                if a == 0: continue\n",
    "                if tuple((filter_params['a_t'],)) < tuple((a,)): \n",
    "                    filtered.append(cont_idx)\n",
    "                    all_holes.append(holes)\n",
    "\n",
    "\n",
    "            foreground_contours = [contours[cont_idx] for cont_idx in filtered]\n",
    "            \n",
    "            hole_contours = []\n",
    "\n",
    "            for hole_ids in all_holes:\n",
    "                unfiltered_holes = [contours[idx] for idx in hole_ids ]\n",
    "                unfilered_holes = sorted(unfiltered_holes, key=cv2.contourArea, reverse=True)\n",
    "                # take max_n_holes largest holes by area\n",
    "                unfilered_holes = unfilered_holes[:filter_params['max_n_holes']]\n",
    "                filtered_holes = []\n",
    "                \n",
    "                # filter these holes\n",
    "                for hole in unfilered_holes:\n",
    "                    if cv2.contourArea(hole) > filter_params['a_h']:\n",
    "                        filtered_holes.append(hole)\n",
    "\n",
    "                hole_contours.append(filtered_holes)\n",
    "\n",
    "            return foreground_contours, hole_contours\n",
    "        \n",
    "        img = np.array(wsi.read_region((0,0), seg_level, wsi.level_dimensions[seg_level]))\n",
    "        img_hsv = cv2.cvtColor(img, cv2.COLOR_RGB2HSV)     # Convert to HSV space\n",
    "        img_med = cv2.medianBlur(img_hsv[:,:,1], mthresh)  # Apply median blurring\n",
    "                \n",
    "        # Thresholding\n",
    "        if use_otsu:\n",
    "            _, img_otsu = cv2.threshold(img_med, 0, sthresh_up, cv2.THRESH_OTSU+cv2.THRESH_BINARY)\n",
    "        else:\n",
    "            _, img_otsu = cv2.threshold(img_med, sthresh, sthresh_up, cv2.THRESH_BINARY)\n",
    "\n",
    "        # Morphological closing\n",
    "        if close > 0:\n",
    "            kernel = np.ones((close, close), np.uint8)\n",
    "            img_otsu = cv2.morphologyEx(img_otsu, cv2.MORPH_CLOSE, kernel)                 \n",
    "\n",
    "        scale = int(wsi.level_dimensions[0][0] / wsi.level_dimensions[seg_level][0])\n",
    "        scaled_ref_patch_area = int(ref_patch_size**2 / (scale * scale))\n",
    "        filter_params = filter_params.copy()\n",
    "        filter_params['a_t'] = filter_params['a_t'] * scaled_ref_patch_area\n",
    "        filter_params['a_h'] = filter_params['a_h'] * scaled_ref_patch_area\n",
    "        \n",
    "        # Find and filter contours\n",
    "        contours, hierarchy = cv2.findContours(img_otsu, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_NONE) # Find contours \n",
    "        hierarchy = np.squeeze(hierarchy, axis=(0,))[:, 2:]\n",
    "        if filter_params:\n",
    "            foreground_contours, hole_contours = _filter_contours(contours, hierarchy, filter_params)  # Necessary for filtering out artifacts\n",
    "\n",
    "        contours_tissue = scale_contour_dim(foreground_contours, scale)\n",
    "        holes_tissue = scale_holes_dim(hole_contours, scale)\n",
    "\n",
    "        if len(keep_ids) > 0:\n",
    "            contour_ids = set(keep_ids) - set(exclude_ids)\n",
    "        else:\n",
    "            contour_ids = set(np.arange(len(contours_tissue))) - set(exclude_ids)\n",
    "\n",
    "        contours_tissue = [contours_tissue[i] for i in contour_ids]\n",
    "        holes_tissue = [holes_tissue[i] for i in contour_ids]\n",
    "        \n",
    "        return contours_tissue, holes_tissue\n",
    "\n",
    "def scale_contour_dim(contours, scale):\n",
    "    return [np.array(cont * scale, dtype='int32') for cont in contours]\n",
    "\n",
    "def scale_holes_dim(contours, scale):\n",
    "    return [[np.array(hole * scale, dtype = 'int32') for hole in holes] for holes in contours]\n",
    "\n",
    "\n",
    "contours_tissue, holes_tissue = segment_tissue(slide)\n",
    "\n",
    "print('Number of contours detected:', len(contours_tissue))\n",
    "print('Number of holes detected:', sum([len(entry) for entry in holes_tissue]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "896a1e7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualizing segmentation quality on thumbnails\n",
    "import math \n",
    "from PIL import Image\n",
    "\n",
    "def vis_wsi(wsi,\n",
    "            contours_tissue,\n",
    "            holes_tissue,\n",
    "            coords=None, \n",
    "            seg_level=2,\n",
    "            color=(0,255,0),\n",
    "            hole_color=(0,0,255),\n",
    "            line_thickness=250,\n",
    "            custom_downsample=1,\n",
    "            number_contours=False,\n",
    "    ):\n",
    "\n",
    "    downsample = int(wsi.level_dimensions[0][0] / wsi.level_dimensions[seg_level][0])\n",
    "    scale = [1/downsample, 1/downsample]\n",
    "    top_left = (0,0)\n",
    "    region_size = wsi.level_dimensions[seg_level]\n",
    "\n",
    "    img = np.array(wsi.read_region(top_left, seg_level, region_size).convert(\"RGB\"))\n",
    "\n",
    "    offset = tuple(-(np.array(top_left) * scale).astype(int))\n",
    "    line_thickness = int(line_thickness * math.sqrt(scale[0] * scale[1]))\n",
    "    \n",
    "    # show contours \n",
    "    if not number_contours:\n",
    "        cv2.drawContours(img, scale_contour_dim(contours_tissue, scale), \n",
    "                         -1, color, line_thickness, lineType=cv2.LINE_8, offset=offset)\n",
    "    else: # add numbering to each contour\n",
    "        for idx, cont in enumerate(contours_tissue):\n",
    "            contour = np.array(scale_contour_dim(cont, scale))\n",
    "            M = cv2.moments(contour)\n",
    "            cX = int(M[\"m10\"] / (M[\"m00\"] + 1e-9))\n",
    "            cY = int(M[\"m01\"] / (M[\"m00\"] + 1e-9))\n",
    "            # draw the contour and put text next to center\n",
    "            cv2.drawContours(img,  [contour], -1, color, line_thickness, lineType=cv2.LINE_8, offset=offset)\n",
    "            cv2.putText(img, \"{}\".format(idx), (cX, cY),\n",
    "                    cv2.FONT_HERSHEY_SIMPLEX, 2, (255, 0, 0), 10)\n",
    "\n",
    "    for holes in holes_tissue:\n",
    "        cv2.drawContours(img, scale_contour_dim(holes, scale), \n",
    "                         -1, hole_color, line_thickness, lineType=cv2.LINE_8)\n",
    "        \n",
    "    if coords is not None:\n",
    "        for coord in coords:\n",
    "            downsample = 16 \n",
    "            patch_size = int(256 / downsample)\n",
    "            x1 = int(coord[0] / downsample)\n",
    "            y1 = int(coord[1] / downsample)\n",
    "            x2 = int(x1 + patch_size)\n",
    "            y2 = int(y1 + patch_size)\n",
    "            cv2.rectangle(img, (x1, y1), (x2, y2), (0,0,0), 2)\n",
    "\n",
    "    img = Image.fromarray(img)\n",
    "\n",
    "    w, h = img.size\n",
    "    if custom_downsample > 1:\n",
    "        img = img.resize((int(w/custom_downsample), int(h/custom_downsample)))\n",
    "\n",
    "    return img\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6baf4ccc",
   "metadata": {},
   "outputs": [],
   "source": [
    "img = vis_wsi(slide, contours_tissue, holes_tissue)\n",
    "img.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "5044450f",
   "metadata": {},
   "source": [
    "### Slide patching\n",
    "\n",
    "Given the tissue segmentation, we will now extract a series of patches of a fixed (pre-determined) size. These patches will form the basis for any downstream application based on deep learning. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f6fe615e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of contours to process:  2\n",
      "Bounding Box: 4896 20160 50945 13953\n",
      "Contour Area: 388059008.0\n",
      "Extracted 5854 coordinates\n",
      "Bounding Box: 320 3664 57105 23633\n",
      "Contour Area: 573257856.0\n",
      "Extracted 8910 coordinates\n"
     ]
    }
   ],
   "source": [
    "# extract coords of all valid patches \n",
    "\n",
    "class IsInContour():\n",
    "    def __init__(self, contour, patch_size, center_shift=0.5):\n",
    "        self.cont = contour\n",
    "        self.patch_size = patch_size\n",
    "        self.shift = int(patch_size//2*center_shift)\n",
    "        \n",
    "    def __call__(self, pt): \n",
    "        center = (pt[0]+self.patch_size//2, pt[1]+self.patch_size//2)\n",
    "        if self.shift > 0:\n",
    "            all_points = [(center[0]-self.shift, center[1]-self.shift),\n",
    "                          (center[0]+self.shift, center[1]+self.shift),\n",
    "                          (center[0]+self.shift, center[1]-self.shift),\n",
    "                          (center[0]-self.shift, center[1]+self.shift)\n",
    "                          ]\n",
    "        else:\n",
    "            all_points = [center]\n",
    "\n",
    "        for points in all_points:\n",
    "            if cv2.pointPolygonTest(self.cont, tuple(np.array(points).astype(float)), False) >= 0:\n",
    "                return 1\n",
    "        return 0\n",
    "\n",
    "def is_in_holes(holes, pt, patch_size):\n",
    "    for hole in holes:\n",
    "        if cv2.pointPolygonTest(hole, (pt[0]+patch_size/2, pt[1]+patch_size/2), False) > 0:\n",
    "            return 1\n",
    "    return 0\n",
    "\n",
    "\n",
    "def is_in_contours(cont_check_fn, pt, holes=None, patch_size=256):\n",
    "    if cont_check_fn(pt):\n",
    "        if holes is not None:\n",
    "            return not is_in_holes(holes, pt, patch_size)\n",
    "        else:\n",
    "            return 1\n",
    "    return 0\n",
    "\n",
    "\n",
    "def process_contours(wsi, contours_tissue, holes_tissue, patch_level=0, patch_size=256, step_size=256):\n",
    "    n_contours = len(contours_tissue)\n",
    "    print(\"Total number of contours to process: \", n_contours)\n",
    "    fp_chunk_size = math.ceil(n_contours * 0.05)\n",
    "    init = True\n",
    "    all_coords = []\n",
    "    for idx, cont in enumerate(contours_tissue):\n",
    "        if (idx + 1) % fp_chunk_size == fp_chunk_size:\n",
    "            print('Processing contour {}/{}'.format(idx, n_contours))\n",
    "        coords = process_contour(wsi, cont, holes_tissue[idx], patch_level, patch_size, step_size)\n",
    "        all_coords.append(coords)\n",
    "    flatten_coords = []\n",
    "    for entry in all_coords:\n",
    "        for coord in entry:\n",
    "            flatten_coords.append(coord)\n",
    "    return flatten_coords\n",
    "\n",
    "\n",
    "def process_contour(\n",
    "        wsi, \n",
    "        cont,\n",
    "        contour_holes,\n",
    "        patch_level,\n",
    "        patch_size=256,\n",
    "        step_size=256,\n",
    "        use_padding=True,\n",
    "    ):\n",
    "    \n",
    "    start_x, start_y, w, h = cv2.boundingRect(cont) \n",
    "    ref_patch_size = (patch_size, patch_size)\n",
    "    \n",
    "    img_w, img_h = wsi.level_dimensions[0]\n",
    "    if use_padding:\n",
    "        stop_y = start_y+h\n",
    "        stop_x = start_x+w\n",
    "    else:\n",
    "        stop_y = min(start_y+h, img_h-ref_patch_size[1]+1)\n",
    "        stop_x = min(start_x+w, img_w-ref_patch_size[0]+1)\n",
    "\n",
    "    print(\"Bounding Box:\", start_x, start_y, w, h)\n",
    "    print(\"Contour Area:\", cv2.contourArea(cont))\n",
    "\n",
    "    cont_check_fn = IsInContour(contour=cont, patch_size=ref_patch_size[0], center_shift=0.5)\n",
    "\n",
    "    step_size_x = step_size \n",
    "    step_size_y = step_size \n",
    "\n",
    "    x_range = np.arange(start_x, stop_x, step=step_size_x)\n",
    "    y_range = np.arange(start_y, stop_y, step=step_size_y)\n",
    "    x_coords, y_coords = np.meshgrid(x_range, y_range, indexing='ij')\n",
    "    coord_candidates = np.array([x_coords.flatten(), y_coords.flatten()]).transpose()\n",
    "\n",
    "    results = []\n",
    "    for coord in coord_candidates:\n",
    "        if is_in_contours(cont_check_fn, coord, contour_holes, ref_patch_size[0]):\n",
    "            results.append(coord) \n",
    "    results = np.array(results)\n",
    "    print('Extracted {} coordinates'.format(len(results)))\n",
    "    return results\n",
    "\n",
    "coords = process_contours(slide, contours_tissue, holes_tissue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "360f661a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualize the patches \n",
    "img = vis_wsi(slide, contours_tissue, holes_tissue, coords)\n",
    "img.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}