\n
If you find this notebook useful, do give me an upvote, it motivates me a lot.
This notebook is still a work in progress. Keep checking for further developments!😊
","metadata":{}},{"cell_type":"markdown","source":"
\n\n
0. MOTIVATION ⤒
\n\n# *Problems*\n* **Whole Slide Images (WSIs)**: large (multiple GB in size) -> need to break up into tiles/patches to further analysis. (preprocessing)\n\n\n\n* **Artefacts due to slide preparation, large portions of background**: need to be removed to speed up the training process and potentially improve performance-> deep tissue detector (artifacts, background, tissue)\n\n\n\n# *Data preprocessing*\n**1. Image Scaling:** scale down whole-slide images (WSIs) in the training data to 1024x1024 resolution\n\n**2. Tissue Segmentation:** identify tissue regions in scaled-down images using the following foreground/background filtering approaches:\n* [Otsu's Binarization](https://docs.opencv.org/4.x/d7/d4d/tutorial_py_thresholding.html)\n* [Triangle Binarization](https://subscription.packtpub.com/book/data/9781789344912/9/ch09lvl1sec80/the-triangle-binarization-algorithm)\n \n**3. Fast Tiling + Tile Selection:** select tiles based on thresholds calculated from step 2. Note that tile size should be large enough that feature relevant to the task are visible to the model to learn\n\n**4. Tile Retrieval + Saving:** retrive and save selected tiles for further analysis (future work)\n\n[**Reference:** Apply filters for tissue segmentation](https://developer.ibm.com/articles/an-automatic-method-to-identify-tissues-from-big-whole-slide-images-pt2/)\n","metadata":{}},{"cell_type":"markdown","source":"
\n\n
1. WSI ANALYSIS ⤒
","metadata":{}},{"cell_type":"markdown","source":"# *Imports*\n\n[OpenSlide](https://openslide.org/) can read virtual slides in the following formats:\n\n* **Aperio (.svs, .tif)**: our image data is in the .tif format\n* Philips (.tiff)\n* Sakura (.svslide)\n* Generic tiled TIFF (.tif)","metadata":{}},{"cell_type":"code","source":"# the OpenSlide package needs C library dependencies \n!apt-get install openslide-tools\n# install openslide-python package in google colab \n!pip install openslide-python","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-08-10T11:01:08.271598Z","iopub.execute_input":"2022-08-10T11:01:08.272275Z","iopub.status.idle":"2022-08-10T11:01:24.616215Z","shell.execute_reply.started":"2022-08-10T11:01:08.272174Z","shell.execute_reply":"2022-08-10T11:01:24.614432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os \nimport sys\nimport logging \nimport time \nimport numpy as np \nimport matplotlib.pyplot as plt \nfrom PIL import Image\nimport tifffile as tiff\nimport pandas as pd\nimport gc\nimport math\nimport cv2\nimport seaborn as sns\n\nimport openslide # see installation guide which requires a 3rd party library\nfrom openslide import open_slide\nfrom openslide.deepzoom import DeepZoomGenerator","metadata":{"execution":{"iopub.status.busy":"2022-08-10T11:01:24.618644Z","iopub.execute_input":"2022-08-10T11:01:24.619012Z","iopub.status.idle":"2022-08-10T11:01:25.615078Z","shell.execute_reply.started":"2022-08-10T11:01:24.618981Z","shell.execute_reply":"2022-08-10T11:01:25.613926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *Helper functions*","metadata":{}},{"cell_type":"code","source":"def rescale_wsi(slide, scale_factor):\n dims = slide.level_dimensions[0] #level 0 is the highest resolution\n large_w, large_h = slide.dimensions\n new_w = math.floor(large_w/scale_factor)\n new_h = math.floor(large_h/scale_factor)\n wsi = slide.read_region((0,0), 0, dims) #read region at level 0\n wsi = wsi.convert(\"RGB\") #since img is in RGBA, need to convert to RGB for further preprocessing\n rescale_img = wsi.resize((new_w, new_h), Image.Resampling.BILINEAR) # PIL image\n return rescale_img\n\ndef create_property_df(slide):\n slide_properties = slide.properties\n descriptions = ['The number of levels in the slide. Levels are numbered from 0 (highest resolution) to level_count - 1 (lowest resolution)',\n 'A list of downsample factors for each level of the slide. level_downsamples[k] is the downsample factor of level k',\n 'Height at level k',\n 'Tile height at level k',\n 'Tile width at level k',\n 'Width at level k',\n 'The name of the property containing an identification of the vendor',\n 'Resolution Unit of the slide',\n 'X Resolution',\n 'Y Resolution',\n 'A (width, height) tuple for level 0 of the slide',\n 'Image mode/attribute'\n ]\n properties = {}\n smaller_region = slide.read_region((slide.dimensions[0]//2,slide.dimensions[1]//2), 0, (1024,1024))\n for key, value in slide_properties.items():\n properties[key] = value\n properties['slide-dimension'] = slide.dimensions # widthxheight\n properties['image-mode'] = smaller_region.mode\n df_properties = pd.DataFrame.from_dict(properties, orient='index').reset_index() #orient index is to create the df using dict keys as rows\n df_properties.columns = ['Slide Property', 'Value']\n df_properties['Description'] = descriptions\n df_properties = df_properties[['Slide Property','Description','Value']]\n return df_properties\n\ndef visualize_slide(df, loc, img):\n plt.figure(figsize=(30,10))\n plt.imshow(img)\n plt.title(f'Label: {df.loc[loc].label}\\nImage ID: {df.loc[loc].image_id}\\nCenter ID: {df.loc[loc].center_id}\\nImage shape: {img.size}')\n plt.show()\n \ndef highlight(row):\n df = lambda x: ['background: #8DE7E3' if x.name in row\n else '' for i in x]\n return df","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2022-08-10T11:01:35.17803Z","iopub.execute_input":"2022-08-10T11:01:35.178474Z","iopub.status.idle":"2022-08-10T11:01:35.192646Z","shell.execute_reply.started":"2022-08-10T11:01:35.178446Z","shell.execute_reply":"2022-08-10T11:01:35.191412Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *Slide properties*","metadata":{}},{"cell_type":"code","source":"#Load the slide file (svs, tif) into an object\nslide = open_slide(\"../input/mayo-clinic-strip-ai/train/00c058_0.tif\")\nrescale_img = rescale_wsi(slide, 32)\ndf = create_property_df(slide)\ndf.style.hide_index().apply(highlight([0,2,5,7,11]), axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-08-10T11:01:37.66426Z","iopub.execute_input":"2022-08-10T11:01:37.665499Z","iopub.status.idle":"2022-08-10T11:04:05.551163Z","shell.execute_reply.started":"2022-08-10T11:01:37.66542Z","shell.execute_reply":"2022-08-10T11:04:05.546334Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *Slide visualization*","metadata":{}},{"cell_type":"code","source":"df = pd.read_csv('../input/mayo-clinic-strip-ai/train.csv')\ndf.head().style.hide_index().apply(highlight([2]), axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-08-10T11:04:05.554315Z","iopub.execute_input":"2022-08-10T11:04:05.554757Z","iopub.status.idle":"2022-08-10T11:04:05.595084Z","shell.execute_reply.started":"2022-08-10T11:04:05.554723Z","shell.execute_reply":"2022-08-10T11:04:05.594017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"visualize_slide(df, 2, rescale_img)","metadata":{"execution":{"iopub.status.busy":"2022-08-10T11:04:05.597025Z","iopub.execute_input":"2022-08-10T11:04:05.597496Z","iopub.status.idle":"2022-08-10T11:04:05.97572Z","shell.execute_reply.started":"2022-08-10T11:04:05.59747Z","shell.execute_reply":"2022-08-10T11:04:05.973643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n\n
2. DATA PREPROCESSING ⤒
","metadata":{}},{"cell_type":"markdown","source":"# *Tissue Segmentation*\n\n## Goal:¶\nMask out non-tissue by setting non-tissue pixels to 0 for their red, green, and blue channels.\n\nConceptually and mathematically, it is often useful to have background values close to or equal to 0 (complement image)\n\nComplement image: simply subtract each pixel value from the maximum pixel value supported by the class (class uint8, the max value of pixel can be 255) and store in the output image array. In the output image, dark areas become lighter and light areas become darker\n\n## Foreground/Background Thresholding Techniques:\n* [Otsu](https://docs.opencv.org/4.x/d7/d4d/tutorial_py_thresholding.html): global image thresholding algorithm, applied to images which are bimodal, which means the image is having 2 peaks in the histogram\n* [Triangle binarization](https://subscription.packtpub.com/book/data/9781789344912/9/ch09lvl1sec80/the-triangle-binarization-algorithm): a line is drawn from the highest bar to the end of the histogram. Then for choosing the optimal threshold, distance of each bar is calculated from the line, whichever is least becomes the threshold value\n\n**Side Note:** [Understanding Histograms in Image Processing](https://towardsdatascience.com/histograms-in-image-processing-with-skimage-python-be5938962935)","metadata":{}},{"cell_type":"markdown","source":"# *Helper functions*","metadata":{}},{"cell_type":"code","source":"# refer to my notebook: https://www.kaggle.com/code/nghihuynh/wsi-preprocessing-tiling-tissue-segmentation/notebook\ndef thresholding(img, method='otsu'):\n # convert to grayscale complement image\n grayscale_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n img_c = 255 - grayscale_img\n thres, thres_img = 0, img_c.copy()\n if method == 'otsu':\n thres, thres_img = cv2.threshold(img_c, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)\n elif method == 'triangle':\n thres, thres_img = cv2.threshold(img_c, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_TRIANGLE)\n return thres, thres_img, img_c\n\ndef histogram(img, thres_img, img_c, thres):\n \"\"\"\n style: ['color', 'grayscale']\n \"\"\" \n plt.figure(figsize=(15,15))\n \n plt.subplot(3,2,1)\n plt.imshow(img)\n plt.title('Scaled-down image')\n \n plt.subplot(3,2,2)\n sns.histplot(img.ravel(), bins=np.arange(0,256), color='orange', alpha=0.5)\n sns.histplot(img[:,:,0].ravel(), bins=np.arange(0,256), color='red', alpha=0.5)\n sns.histplot(img[:,:,1].ravel(), bins=np.arange(0,256), color='Green', alpha=0.5)\n sns.histplot(img[:,:,2].ravel(), bins=np.arange(0,256), color='Blue', alpha=0.5)\n plt.legend(['Total', 'Red_Channel', 'Green_Channel', 'Blue_Channel'])\n plt.ylim(0,0.05e6)\n plt.xlabel('Intensity value')\n plt.title('Color histogram')\n \n plt.subplot(3,2,3)\n plt.imshow(img_c, cmap='gist_gray')\n plt.title('Complement grayscale image')\n \n plt.subplot(3,2,4)\n sns.histplot(img_c.ravel(), bins=np.arange(0,256))\n plt.axvline(thres, c='red', linestyle=\"--\")\n plt.ylim(0,0.05e6)\n plt.xlabel('Intensity value')\n plt.title('Grayscale complement histogram')\n \n plt.subplot(3,2,5)\n plt.imshow(thres_img, cmap='gist_gray')\n plt.title('Thresholded image')\n \n plt.subplot(3,2,6)\n sns.histplot(thres_img.ravel(), bins=np.arange(0,256))\n plt.axvline(thres, c='red', linestyle=\"--\")\n plt.ylim(0,0.05e6)\n plt.xlabel('Intensity value')\n plt.title('Thresholded histogram')\n \n plt.tight_layout()\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-08-10T11:05:04.748827Z","iopub.execute_input":"2022-08-10T11:05:04.749852Z","iopub.status.idle":"2022-08-10T11:05:04.842059Z","shell.execute_reply.started":"2022-08-10T11:05:04.74981Z","shell.execute_reply":"2022-08-10T11:05:04.839185Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *Otsu's binarization*","metadata":{}},{"cell_type":"code","source":"img_np = np.array(rescale_wsi(slide, 16))","metadata":{"execution":{"iopub.status.busy":"2022-08-10T11:05:07.78618Z","iopub.execute_input":"2022-08-10T11:05:07.786511Z","iopub.status.idle":"2022-08-10T11:06:54.889167Z","shell.execute_reply.started":"2022-08-10T11:05:07.786484Z","shell.execute_reply":"2022-08-10T11:06:54.887999Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"thres_otsu, thres_img, img_c = thresholding(img_np, method='otsu')\nhistogram(img_np, thres_img, img_c, thres_otsu)","metadata":{"execution":{"iopub.status.busy":"2022-08-10T11:06:54.890925Z","iopub.execute_input":"2022-08-10T11:06:54.891248Z","iopub.status.idle":"2022-08-10T11:07:15.899479Z","shell.execute_reply.started":"2022-08-10T11:06:54.891188Z","shell.execute_reply":"2022-08-10T11:07:15.898419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *Triangle binarization*","metadata":{}},{"cell_type":"code","source":"thres_triangle, thres_img, img_c = thresholding(img_np, method='triangle')\nhistogram(img_np, thres_img, img_c, thres_triangle)","metadata":{"execution":{"iopub.status.busy":"2022-08-10T11:07:15.900346Z","iopub.execute_input":"2022-08-10T11:07:15.900593Z","iopub.status.idle":"2022-08-10T11:07:37.850824Z","shell.execute_reply.started":"2022-08-10T11:07:15.900565Z","shell.execute_reply":"2022-08-10T11:07:37.849622Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n\n","metadata":{}},{"cell_type":"markdown","source":"# *Helper functions*","metadata":{}},{"cell_type":"code","source":"def tile_properties(tiles):\n tile_properties = {'Level count': ['The number of Deep Zoom levels in the image',tiles.level_count], \n 'Tile count': ['The total number of Deep Zoom tiles in the image',tiles.tile_count],\n 'Level tiles': ['A list of (tiles_x, tiles_y) tuples for each Deep Zoom level. level_tiles[k] are the tile counts of level k',tiles.level_tiles],\n 'Level dimensions': ['A list of (pixels_x, pixels_y) tuples for each Deep Zoom level. level_dimensions[k] are the dimensions of level k',tiles.level_dimensions]\n }\n tile_df = pd.DataFrame.from_dict(tile_properties, orient='index').reset_index()\n tile_df.columns = ['Tile property', 'Description', 'Value']\n return tile_df\n\ndef visualize_tiles(tiles, level, thres, tile_size):\n rows, cols = tiles.level_tiles[level]\n plt.figure(figsize=(10,10))\n count = 0\n img_size = 0\n for row in range(rows):\n for col in range(cols):\n count += 1\n plt.subplot(4,4,count)\n plt.xticks([])\n plt.yticks([])\n single_tile = tiles.get_tile(level, (row,col))\n single_tile_RGB = single_tile.convert('RGB')\n img_test = np.array(single_tile_RGB.resize((tile_size, tile_size), Image.Resampling.BILINEAR))\n img_size = img_test.shape\n plt.imshow(img_test)\n plt.tight_layout()\n plt.suptitle(f'Tile level: {level}\\nNum tiles: {count}\\nTile size: {img_size}')\n plt.tight_layout()\n plt.show()\n \n plt.figure(figsize=(10,10))\n count = 0\n selected_tiles = 0\n img_size = 0\n for row in range(rows):\n for col in range(cols):\n count += 1\n plt.subplot(4,4,count)\n plt.xticks([])\n plt.yticks([])\n single_tile = tiles.get_tile(level, (row,col))\n single_tile_RGB = single_tile.convert('RGB')\n img_test = np.array(single_tile_RGB.resize((tile_size, tile_size), Image.Resampling.BILINEAR))\n img_size = img_test.shape\n img_c = 255 - img_test\n if img_c.mean() > thres:\n selected_tiles += 1\n plt.imshow(img_test)\n plt.tight_layout()\n plt.suptitle(f'Tile level: {level}\\nNum selected tiles: {selected_tiles}\\nTile size: {img_size}')\n plt.tight_layout()\n plt.show()\n ","metadata":{"execution":{"iopub.status.busy":"2022-08-10T11:07:37.853493Z","iopub.execute_input":"2022-08-10T11:07:37.853842Z","iopub.status.idle":"2022-08-10T11:07:37.868033Z","shell.execute_reply.started":"2022-08-10T11:07:37.853815Z","shell.execute_reply":"2022-08-10T11:07:37.867184Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# *Tile properties*","metadata":{}},{"cell_type":"code","source":"tiles = DeepZoomGenerator(slide, tile_size=1024, overlap=0, limit_bounds=False)\ntile_df = tile_properties(tiles)\ntile_df.style.hide_index().apply(highlight([1,3]), axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-08-10T11:07:37.869413Z","iopub.execute_input":"2022-08-10T11:07:37.870354Z","iopub.status.idle":"2022-08-10T11:07:37.899656Z","shell.execute_reply.started":"2022-08-10T11:07:37.87029Z","shell.execute_reply":"2022-08-10T11:07:37.898333Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> #### *Working with lower level dimensions might be beneficial for quick prototyping. In contrast, we should be cautious that higher level dimensions yield more tiles with a given tile size. If we worry about OOR, we should start with lower dimensions and build up from that.*","metadata":{}},{"cell_type":"markdown","source":"# *Tile selection*","metadata":{}},{"cell_type":"code","source":"visualize_tiles(tiles, 13, thres_triangle, 1024)","metadata":{"execution":{"iopub.status.busy":"2022-08-10T11:07:37.90102Z","iopub.execute_input":"2022-08-10T11:07:37.901491Z","iopub.status.idle":"2022-08-10T11:10:49.252017Z","shell.execute_reply.started":"2022-08-10T11:07:37.901463Z","shell.execute_reply":"2022-08-10T11:10:49.250578Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> #### *We can discard many uneccessary tiles such as background or artefacts after these thresholding techniques. However, there's still a large portion of background in the selected tiles. We might need to choose smaller tile sizes to increase the quality of the training data.* \n> #### *WSI are often labeled at the slide level, not at the tile level. Thus, we might need to annotate them at tile level to derive the ground truth labels of individual tiles from slide-wide label data*","metadata":{}},{"cell_type":"markdown","source":"
\n\n
4. CONCLUSION ⤒
\n\n## Recap: WSI preprocessing framework:\n### Steps:\n1. **Scale** down images\n\n2. **Filter** foreground/background (tissue/non-tissue)\n\n3. **Break** down images into **tiles**, and **select** tiles containing large portion of **tissue regions**\n\n4. **Retrieve** and **save** selected tiles for further analysis (future work)\n\n### Remarks:\n\n* Working with **lower level dimensions** might be beneficial for quick prototyping. In contrast, we should be cautious that **higher level dimensions** yield more tiles with a given tile size. If we worry about OOR, we should start with lower dimensions and build up from that.\n* We can discard many **uneccessary tiles** such as background or artefacts after these thresholding techniques. However, there's still a large portion of background in the selected tiles. We might need to choose smaller **tile sizes** to increase the quality of the training data.\n* WSI are often labeled at the **slide level**, not at the tile level. Thus, we might need to **annotate** them at **tile level** to derive the ground truth labels of individual tiles from slide-wide label data","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}