{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"#  Medical Whole Slide Processing: Speed Up!","metadata":{}},{"cell_type":"markdown","source":"The first challenge with **medical whole slide images** is probably their huge size and pre-processing required to get them prepared for model training and further analysis. Major concerns are:\n* Generating Tiles\n* Thresholding (discard tiles with no structural data, e.g. mostly background tiles)\n* Resource management (RAM and Storage)\n* Time Management (Heavy Computational Processes)\n\nHere I share my experience with some tricky slides from [this kaggle competition](https://www.kaggle.com/competitions/mayo-clinic-strip-ai) and will discuss some solutions.\n\nTo view the code in detail you will need to open the notebook in kaggle or view it from google colab.","metadata":{}},{"cell_type":"code","source":"import warnings\nwarnings.filterwarnings('ignore')\n\n\n# install fast kaggle\ntry: import fastkaggle\nexcept ModuleNotFoundError:\n    !pip install -Uqq fastkaggle\n    \nfrom fastkaggle import *","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:13:43.622673Z","iopub.execute_input":"2022-09-05T03:13:43.624099Z","iopub.status.idle":"2022-09-05T03:13:43.70823Z","shell.execute_reply.started":"2022-09-05T03:13:43.62399Z","shell.execute_reply":"2022-09-05T03:13:43.706932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# install pyvips\n!conda install -y -qq --channel conda-forge pyvips > quiet.txt","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:13:43.719318Z","iopub.execute_input":"2022-09-05T03:13:43.719743Z","iopub.status.idle":"2022-09-05T03:14:52.488928Z","shell.execute_reply.started":"2022-09-05T03:13:43.719705Z","shell.execute_reply":"2022-09-05T03:14:52.487302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# set up the competition\ncomp = 'mayo-clinic-strip-ai'\npath = setup_comp(comp, install='\"fastcore>=1.4.5\" \"fastai>=2.7.1\" \"timm>=0.6.2.dev0\"')\npath.ls()","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:14:52.493565Z","iopub.execute_input":"2022-09-05T03:14:52.493986Z","iopub.status.idle":"2022-09-05T03:15:04.591947Z","shell.execute_reply.started":"2022-09-05T03:14:52.493948Z","shell.execute_reply":"2022-09-05T03:15:04.590645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from fastai.vision.all import *\nfrom fastcore.parallel import *\nimport pyvips","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:15:04.596167Z","iopub.execute_input":"2022-09-05T03:15:04.597562Z","iopub.status.idle":"2022-09-05T03:15:08.345644Z","shell.execute_reply.started":"2022-09-05T03:15:04.597482Z","shell.execute_reply":"2022-09-05T03:15:08.344061Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trn_path = path/'train'\ntrn_slides = get_image_files(trn_path)\n\ntst_path = path/'test'\ntst_slides = get_image_files(tst_path)","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:15:08.347933Z","iopub.execute_input":"2022-09-05T03:15:08.349239Z","iopub.status.idle":"2022-09-05T03:15:08.529259Z","shell.execute_reply.started":"2022-09-05T03:15:08.349192Z","shell.execute_reply":"2022-09-05T03:15:08.527802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## What do these slides look like anyway?","metadata":{}},{"cell_type":"markdown","source":"Let’s get on to it and see what we are dealing with.","metadata":{}},{"cell_type":"code","source":"slides = tst_slides[3], tst_slides[1]\nfig = plt.figure(figsize=(12,8))\n\nax = fig.add_subplot(1,2,1)\nwst = pyvips.Image.thumbnail(slides[0], 400)\nwsi = pyvips.Image.new_from_file(slides[0])\nax.set_title(f'{wsi.width} x {wsi.height} pixel')\nimplot = plt.imshow(wst)\n\nax = fig.add_subplot(1,2,2)\nwst = pyvips.Image.thumbnail(slides[1], 400)\nwsi = pyvips.Image.new_from_file(slides[1])\nax.set_title(f'{wsi.width} x {wsi.height} pixel')\nimplot = plt.imshow(wst)","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:15:08.532719Z","iopub.execute_input":"2022-09-05T03:15:08.533412Z","iopub.status.idle":"2022-09-05T03:16:05.845532Z","shell.execute_reply.started":"2022-09-05T03:15:08.533372Z","shell.execute_reply":"2022-09-05T03:16:05.844216Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Well, if nothing they are huge in size (pixel) and colour tone could be very different for each slide apparently because of the microscope, lab lighting, and some other factors which is not my concern here.","metadata":{}},{"cell_type":"markdown","source":"## What is a Tile about?","metadata":{}},{"cell_type":"markdown","source":"To make these whole slides manageable, we need to break them into smaller pieces. Well it would be worth of a discussion to ask:\n* What tile size is OK?\n* How many of them? Is sampling OK?\n* What should each tile be labelled?\n\nVery important questions I suppose, but not for this blog post.\n\nLet’s have a closer look at some of those tiles from one slide.","metadata":{}},{"cell_type":"code","source":"slide = tst_slides[1]\nwsi = pyvips.Image.new_from_file(slide)\ntiles = [wsi.crop(1500, 3000, 2000, 2000),\n         wsi.crop(3500, 3000, 2000, 2000),\n         wsi.crop(5000, 3000, 2000, 2000)]\n_,axs = subplots(1, 3)\nfor i, ax in enumerate(axs):\n    show_image(tiles[i], ctx=ax, title=f'tile-0{i+1}')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:05.847024Z","iopub.execute_input":"2022-09-05T03:16:05.847458Z","iopub.status.idle":"2022-09-05T03:16:07.959772Z","shell.execute_reply.started":"2022-09-05T03:16:05.847414Z","shell.execute_reply":"2022-09-05T03:16:07.958291Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Obviously no one wants to waste time and compute resources on background tiles like **tile-03** which do not contain any structural data.\n\nWhat about **tile-02**? Is it important to keep that one?\n\nWhat if you have hundreds of tiles like **tile-01** from one slide? Is it still worth it to keep **tile-02**?\n\nI would say there is no definite answer unless through experiments!\n\nFor now let’s see how thresholding is going to help.","metadata":{}},{"cell_type":"markdown","source":"## Thresholding: When Things Do Not Work as Expected!","metadata":{}},{"cell_type":"markdown","source":"Thresholding  is basically separating foreground from background. So it is used to discard the tiles or part of the image that does not contain any useful data. [skimage](https://scikit-image.org/docs/stable/auto_examples/applications/plot_thresholding.html) is probably  a good place to dive into more details.\n\nAnyway, there are different algorithms for thresholding and in practice it becomes to some extent clear why so many! Because they do not perform consistently depending on the image they are processing.\n\nThresholding is usually performed at grayscale (single band) level. One fun fact to know is that medical images do not necessarily look grey at grayscale :)","metadata":{}},{"cell_type":"code","source":"_,axs = subplots(1, 3)\nfor i, ax in enumerate(axs):\n    show_image(tiles[i].colourspace('b-w'), ctx=ax, title=f'tile-0{i+1}-grey')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:07.961581Z","iopub.execute_input":"2022-09-05T03:16:07.963847Z","iopub.status.idle":"2022-09-05T03:16:09.76944Z","shell.execute_reply.started":"2022-09-05T03:16:07.963793Z","shell.execute_reply":"2022-09-05T03:16:09.768271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let’s start with a fairly simple method known as [threshold minimum](https://scikit-image.org/docs/stable/api/skimage.filters.html#skimage.filters.threshold_minimum).","metadata":{}},{"cell_type":"code","source":"from skimage.filters import threshold_minimum\n\n_,axs = subplots(1, 3)\nfor i, ax in enumerate(axs):\n    tile = tiles[i].colourspace('b-w')\n    tile = tile > threshold_minimum(tile.numpy())\n    show_image( tile, ctx=ax, title=f'tile-0{i+1} ({tile.avg()/255:.0%})')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:09.77091Z","iopub.execute_input":"2022-09-05T03:16:09.771251Z","iopub.status.idle":"2022-09-05T03:16:12.738863Z","shell.execute_reply.started":"2022-09-05T03:16:09.771219Z","shell.execute_reply":"2022-09-05T03:16:12.737554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Well, the first thing to point out is that I am apparently measuring the background rather than foreground. By definition, `image > threshold` should grab the foreground. Which does not seem to be the case here.\n\nApart from that, the results seem promising! Don't they!?\n","metadata":{}},{"cell_type":"markdown","source":"Let’s grab another set of tiles from a different slide and see how this method works out!\n\nThese are the new tiles:","metadata":{}},{"cell_type":"code","source":"slide = tst_slides[3]\nwsi = pyvips.Image.new_from_file(slide)\ntiles = [wsi.crop(4000, 14000, 2000, 2000),\n         wsi.crop(10000, 12000, 2000, 2000),\n         wsi.crop(10000, 16000, 2000, 2000)]\n\n_,axs = subplots(1, 3)\nfor i, ax in enumerate(axs):\n    show_image(tiles[i], ctx=ax, title=f'tile-0{i+4}')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:12.740435Z","iopub.execute_input":"2022-09-05T03:16:12.741241Z","iopub.status.idle":"2022-09-05T03:16:15.016368Z","shell.execute_reply.started":"2022-09-05T03:16:12.741187Z","shell.execute_reply":"2022-09-05T03:16:15.014975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And these are the thresholding results:","metadata":{}},{"cell_type":"code","source":"_,axs = subplots(1, 3)\nfor i, ax in enumerate(axs):\n    tile = tiles[i].colourspace('b-w')\n    tile = tile > threshold_minimum(tile.numpy())\n    show_image( tile, ctx=ax, title=f'tile-0{i+4} ({tile.avg()/255:.0%})')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:15.018015Z","iopub.execute_input":"2022-09-05T03:16:15.018412Z","iopub.status.idle":"2022-09-05T03:16:17.814198Z","shell.execute_reply.started":"2022-09-05T03:16:15.018373Z","shell.execute_reply":"2022-09-05T03:16:17.812786Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I expected both **tile-05** and **tile-06** would be discarded as garbage, but the method begs to differ. Now what? ","metadata":{}},{"cell_type":"markdown","source":"## There Should be Better Methods to Do This!?","metadata":{}},{"cell_type":"markdown","source":"There are various thresholding methods implemented by `skimage` that you can give a try. I am going to try one or two more methods here. Let's try a well known method called [**otsu thresholding**](https://scikit-image.org/docs/stable/api/skimage.filters.html#skimage.filters.threshold_otsu).\n\nLet’s see how `otsu` will perform on the latest tiles.","metadata":{}},{"cell_type":"code","source":"from skimage.filters import threshold_otsu\n\n_,axs = subplots(1, 3)\nfor i, ax in enumerate(axs):\n    tile = tiles[i].colourspace('b-w')\n    tile = tile > threshold_otsu(tile.numpy())\n    show_image( tile, ctx=ax, title=f'tile-0{i+4} (otsu: {tile.avg()/255:.0%})')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:17.816233Z","iopub.execute_input":"2022-09-05T03:16:17.816742Z","iopub.status.idle":"2022-09-05T03:16:20.598812Z","shell.execute_reply.started":"2022-09-05T03:16:17.816694Z","shell.execute_reply":"2022-09-05T03:16:20.597283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Does not sound good! Does it!? Gorgeous output though! :)\n\nThere is another method I like to try which is called [yen thresholding](https://scikit-image.org/docs/stable/api/skimage.filters.html#skimage.filters.threshold_yen) method.\n\nThe output looks like below:","metadata":{}},{"cell_type":"code","source":"from skimage.filters import threshold_yen\n\n_,axs = subplots(1, 3)\nfor i, ax in enumerate(axs):\n    tile = tiles[i].colourspace('b-w')\n    tile = tile > threshold_yen(tile.numpy())\n    show_image( tile, ctx=ax, title=f'tile-0{i+4} (yen: {tile.avg()/255:.0%})')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:20.603948Z","iopub.execute_input":"2022-09-05T03:16:20.604579Z","iopub.status.idle":"2022-09-05T03:16:23.391543Z","shell.execute_reply.started":"2022-09-05T03:16:20.604521Z","shell.execute_reply":"2022-09-05T03:16:23.390348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"On one hand it has successfully recognised both background tiles, but on the other hand it tends to evaluate low the **tile-04** which is mostly of valuable data. I could set the discard margin like 20% and live with it. But is this a consistent result really? You have to keep trying!\n\nLet's give it another try on a third slide.","metadata":{}},{"cell_type":"code","source":"from skimage.filters import threshold_yen\n\nslide = trn_slides[157]\nwsi = pyvips.Image.new_from_file(slide)\ntiles = [wsi.crop(14000, 20000, 2000, 2000),\n         wsi.crop(14000, 25000, 2000, 2000),\n         wsi.crop(14000, 27000, 2000, 2000),\n         wsi.crop(13000, 28000, 2000, 2000)]\n\n_,axs = subplots(1, 4)\nfor i, ax in enumerate(axs):\n    show_image(tiles[i], ctx=ax, title=f'tile-0{i+7}')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:23.392986Z","iopub.execute_input":"2022-09-05T03:16:23.393777Z","iopub.status.idle":"2022-09-05T03:16:26.071713Z","shell.execute_reply.started":"2022-09-05T03:16:23.393739Z","shell.execute_reply":"2022-09-05T03:16:26.070519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Funny right! But they are honestly all from a single slide! :)\n\nLet check out the `yen` method:\n","metadata":{}},{"cell_type":"code","source":"slide = trn_slides[157]\nwsi = pyvips.Image.new_from_file(slide)\ntiles = [wsi.crop(14000, 20000, 2000, 2000),\n         wsi.crop(14000, 25000, 2000, 2000),\n         wsi.crop(14000, 27000, 2000, 2000),\n         wsi.crop(13000, 28000, 2000, 2000)]\n\n_,axs = subplots(1, 4)\nfor i, ax in enumerate(axs):\n    tile = tiles[i].colourspace('b-w')\n    tile = tile > threshold_yen(tile.numpy())\n    show_image( tile, ctx=ax, title=f'tile-0{i+7} (yen: {tile.avg()/255:.0%})')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:26.073463Z","iopub.execute_input":"2022-09-05T03:16:26.074032Z","iopub.status.idle":"2022-09-05T03:16:29.451405Z","shell.execute_reply.started":"2022-09-05T03:16:26.073993Z","shell.execute_reply":"2022-09-05T03:16:29.450048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Mission accomplished! The `yen thresholding method` is broken successfully!\n\nBut it is not the end of the world yet.","metadata":{}},{"cell_type":"markdown","source":"## Morphology: Robust but Expensive!!","metadata":{}},{"cell_type":"markdown","source":"There is a totally separate chapter on the image processing as called **Morphology**. It is a collection of `non-linear` operations to manipulate an image. I suppose [skimage documentation](https://scikit-image.org/docs/stable/auto_examples/applications/plot_morphology.html) would be a good start.\n\nI did a bit of research and apparently a very common pipeline is dilation or erosion over an image with edges marked/improved/whatever!\n\nI will implement a `canny()` edge detector plus a `binary_dilation()` operation with a `mask` of size `11`.","metadata":{}},{"cell_type":"code","source":"from skimage.feature import canny\nfrom skimage.morphology import binary_dilation, disk\n\ndef skimage_morph(tile):\n    tile = tile.colourspace('b-w').numpy()\n    tile = canny(tile)\n    tile = binary_dilation(tile, disk(6))\n    return tile","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:29.453209Z","iopub.execute_input":"2022-09-05T03:16:29.454425Z","iopub.status.idle":"2022-09-05T03:16:29.568612Z","shell.execute_reply.started":"2022-09-05T03:16:29.454374Z","shell.execute_reply":"2022-09-05T03:16:29.567388Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see how it will perform on all above tiles!","metadata":{}},{"cell_type":"code","source":"slide = tst_slides[1]\nwsi = pyvips.Image.new_from_file(slide)\ntiles = [wsi.crop(1500, 3000, 2000, 2000),\n         wsi.crop(3500, 3000, 2000, 2000),\n         wsi.crop(5000, 3000, 2000, 2000)]\n_,axs = subplots(1, 3)\nfor i, ax in enumerate(axs):\n    tile = skimage_morph(tiles[i])\n    show_image(tile, ctx=ax, title=f'tile-0{i+1} ({tile.mean():.0%})')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:29.56995Z","iopub.execute_input":"2022-09-05T03:16:29.570325Z","iopub.status.idle":"2022-09-05T03:16:35.996217Z","shell.execute_reply.started":"2022-09-05T03:16:29.57029Z","shell.execute_reply":"2022-09-05T03:16:35.99449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slide = tst_slides[3]\nwsi = pyvips.Image.new_from_file(slide)\ntiles = [wsi.crop(4000, 14000, 2000, 2000),\n         wsi.crop(10000, 12000, 2000, 2000),\n         wsi.crop(10000, 16000, 2000, 2000)]\n_,axs = subplots(1, 3)\nfor i, ax in enumerate(axs):\n    tile = skimage_morph(tiles[i])\n    show_image(tile, ctx=ax, title=f'tile-0{i+4} ({tile.mean():.0%})')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:35.998113Z","iopub.execute_input":"2022-09-05T03:16:35.99883Z","iopub.status.idle":"2022-09-05T03:16:43.043949Z","shell.execute_reply.started":"2022-09-05T03:16:35.998776Z","shell.execute_reply":"2022-09-05T03:16:43.042701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"slide = trn_slides[157]\nwsi = pyvips.Image.new_from_file(slide)\ntiles = [wsi.crop(14000, 20000, 2000, 2000),\n         wsi.crop(14000, 25000, 2000, 2000),\n         wsi.crop(14000, 27000, 2000, 2000),\n         wsi.crop(13000, 28000, 2000, 2000)]\n_,axs = subplots(1, 4)\nfor i, ax in enumerate(axs):\n    tile = skimage_morph(tiles[i])\n    show_image(tile, ctx=ax, title=f'tile-0{i+7} ({tile.mean():.0%})')","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:43.045519Z","iopub.execute_input":"2022-09-05T03:16:43.045926Z","iopub.status.idle":"2022-09-05T03:16:52.225315Z","shell.execute_reply.started":"2022-09-05T03:16:43.045888Z","shell.execute_reply":"2022-09-05T03:16:52.223899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Voila! Can not be better! Right?\n\nWell not so fast!\n\nLet’s put it in a more practical test and see how it will perform in terms of both the number of tiles generated and required time!","metadata":{}},{"cell_type":"markdown","source":"## “skimage” Time Based Test","metadata":{}},{"cell_type":"markdown","source":"I am going to run the pipeline on the slide below which is quite massive and messy! Kudos to the lab technician! :)","metadata":{}},{"cell_type":"code","source":"slide = trn_slides[157]\nwst = pyvips.Image.thumbnail(slide, 600)\nwsi = pyvips.Image.new_from_file(slide)\ntitle = f'{wsi.width} x {wsi.height} pixel'\nshow_image(wst, title=title, figsize=(8,10));","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:16:52.227052Z","iopub.execute_input":"2022-09-05T03:16:52.227408Z","iopub.status.idle":"2022-09-05T03:19:10.002582Z","shell.execute_reply.started":"2022-09-05T03:16:52.227376Z","shell.execute_reply":"2022-09-05T03:19:10.001269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I am running the test on a **Kaggle** notebook instance. At the time of writing, it is powered up by a Linux machine with 4 `CPU core`s (allowing two threads per core) and  16GB of `RAM`. I also run the pipeline through a parallel process of `ThreadPoolExcuter` to utilise the resources as much as possible.\n\nWhat I am doing here is basically cropping the slide in 2k x 2k pixels, downsizing to scale of 0.25 if contains like 40% blood clot, then save the tile.\n\n**One more thing!**\n\nDuring my research I found out that apparently with medical slides, the green band (of `RGB ` spectrum image) contains the most structural (/useful) data. So I will use the green layer here instead of converting it to `grayscale`","metadata":{}},{"cell_type":"code","source":"slide = trn_slides[176]\nwsi = pyvips.Image.new_from_file(slide)\ntile = wsi.crop(2000, 1000, 2000, 2000)\ntiles = [tile, tile.colourspace('b-w'), tile[1]]\ntitles = ['Original', 'Greyscale', 'G band of RGB']\n_,axs = subplots(1, 3)\nfor i, ax in enumerate(axs):\n    show_image(tiles[i], ctx=ax, title=titles[i])","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:19:10.004285Z","iopub.execute_input":"2022-09-05T03:19:10.004704Z","iopub.status.idle":"2022-09-05T03:19:12.075412Z","shell.execute_reply.started":"2022-09-05T03:19:10.004667Z","shell.execute_reply":"2022-09-05T03:19:12.073911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!rm -rf /tmp/tiles/\n\ntpath = Path('/tmp/tiles/')\nfolders = ['skimage', 'pyvips']\n\nfor f in folders:\n    if not (tpath/f).exists():\n        (tpath/f).mkdir(exist_ok=True, parents=True)\n\ndef skimage_pct(tile, thresh=0.4):\n    tile = tile[1].numpy()\n    tile = canny(tile)\n    tile = binary_dilation(tile, disk(6))\n    return tile.mean() > thresh\n    \nmask = pyvips.Image.mask_ideal(11, 11, 1, optical=True, reject=True)\nmask = (mask * 128 + 128).cast(\"uchar\").copy_memory()\n\ndef pyvips_pct(tile, thresh=0.4, mask=mask):\n    tile = tile[1].canny(precision='integer').dilate(mask)\n    return tile.avg()/255 > thresh\n\ndef get_tiles(slide, t_size=2_000):\n    slide = pyvips.Image.new_from_file(slide)\n    return (\n        slide.crop(x, y, min(t_size, slide.width - x), min(t_size, slide.height - y)\n                   ).gravity('south-east', t_size, t_size, extend='repeat')\n        for y in range(0, slide.height, t_size)\n        for x in range(0, slide.width, t_size))\n\n\ndef save_jpeg(tiles, tname, clot_pct, folder, tcount=1000):\n    for tile in tiles:\n        if clot_pct(tile):\n            tile = tile.resize(0.25, kernel='linear')\n            fname = f'{tname}_{tcount}.jpg'\n            tile.write_to_file(tpath/folder/fname)\n            tcount += 1\n            \ndef make_jpeg_skimage(slide, tname='skimage', clot_pct=skimage_pct, folder='skimage'):\n    tiles = get_tiles(slide)\n    save_jpeg(tiles=tiles, tname=tname, clot_pct=clot_pct, folder=folder)\n    \ndef make_jpeg_pyvips(slide, tname='pyvips', clot_pct=pyvips_pct, folder='pyvips'):\n    tiles = get_tiles(slide)\n    save_jpeg(tiles=tiles, tname=tname, clot_pct=clot_pct, folder=folder)\n    \ndef walltime(s, e):\n    t = e - s\n    m = int(t//60)\n    s = t - (m * 60)\n    if m==0: return f\"{s:.0f}s\"\n    else: return f\"{m}min {s:.0f}s\"\n    \nelapsed=list()\nt_num=list()","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:19:12.077786Z","iopub.execute_input":"2022-09-05T03:19:12.078303Z","iopub.status.idle":"2022-09-05T03:19:13.218751Z","shell.execute_reply.started":"2022-09-05T03:19:12.078254Z","shell.execute_reply":"2022-09-05T03:19:13.217146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"start = time.time()\nparallel(make_jpeg_skimage, trn_slides[157:158], n_workers=8, progress=False, threadpool=True)\nend = time.time()\nelapsed.append(walltime(s=start, e=end))\nf = get_image_files(tpath/'skimage')\nt_num.append(len(f))","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:19:13.220819Z","iopub.execute_input":"2022-09-05T03:19:13.22191Z","iopub.status.idle":"2022-09-05T03:49:49.932249Z","shell.execute_reply.started":"2022-09-05T03:19:13.221815Z","shell.execute_reply":"2022-09-05T03:49:49.930729Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"OK, let’s see how the performance is!","metadata":{}},{"cell_type":"code","source":"df = pd.DataFrame({\n    'library': ['skimage'],\n    'time elapsed': elapsed,\n    'tiles generated': t_num\n})\ndf","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:49:49.934746Z","iopub.execute_input":"2022-09-05T03:49:49.936244Z","iopub.status.idle":"2022-09-05T03:49:49.96547Z","shell.execute_reply.started":"2022-09-05T03:49:49.936183Z","shell.execute_reply":"2022-09-05T03:49:49.964465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To my expectation, that is terrible!\n\nConsidering a collection of 1,000 slides and 10 mins average on each slide (being optimistic), we are talking about days of pre-processing. With the trial and experimental nature of deep learning, this does not sound practical at all.\n","metadata":{}},{"cell_type":"markdown","source":"The tiles generated look good though!","metadata":{}},{"cell_type":"code","source":"def label_func(f): return f.stem\n\ndbl = DataBlock(\n    blocks = (ImageBlock, CategoryBlock),\n    get_items = get_image_files,\n    get_y = label_func)\n\n\ndls = dbl.dataloaders(tpath/'skimage', bs=64)\n\ndls.show_batch(max_n=12)","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:49:49.96942Z","iopub.execute_input":"2022-09-05T03:49:49.970585Z","iopub.status.idle":"2022-09-05T03:49:52.512077Z","shell.execute_reply.started":"2022-09-05T03:49:49.97054Z","shell.execute_reply":"2022-09-05T03:49:52.510598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## “pyvips” is Out There to Astonish  You!","metadata":{}},{"cell_type":"markdown","source":"Before I run the test for [pyvips](https://libvips.github.io/pyvips/vimage.html), I must make it clear that even with [skimage](https://scikit-image.org/docs/stable/), the backbone processes like loading those massive slides in memory, cropping, and grabbing tiles all have been handled  by `pyvips`.","metadata":{}},{"cell_type":"code","source":"start = time.time()\nparallel(make_jpeg_pyvips, trn_slides[157:158], n_workers=8, progress=False, threadpool=True)\nend = time.time()\nelapsed.append(walltime(s=start, e=end))\nf = get_image_files(tpath/'skimage')\nt_num.append(len(f))","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:49:52.514037Z","iopub.execute_input":"2022-09-05T03:49:52.514814Z","iopub.status.idle":"2022-09-05T03:56:54.75615Z","shell.execute_reply.started":"2022-09-05T03:49:52.514774Z","shell.execute_reply":"2022-09-05T03:56:54.754586Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Having said that, let’s see how `pyvips` performs on same condition:","metadata":{}},{"cell_type":"code","source":"df = pd.DataFrame({\n    'library': folders,\n    'time elapsed': elapsed,\n    'tiles generated': t_num\n})\ndf","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:56:54.759Z","iopub.execute_input":"2022-09-05T03:56:54.759576Z","iopub.status.idle":"2022-09-05T03:56:54.776851Z","shell.execute_reply.started":"2022-09-05T03:56:54.759487Z","shell.execute_reply":"2022-09-05T03:56:54.775558Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That is something! Isn’t it? That is about five times quicker. \n\n`pyvips` has generated the same number of tiles as well. They are for sure up to expectation:\n","metadata":{}},{"cell_type":"code","source":"dbl = DataBlock(\n    blocks = (ImageBlock, CategoryBlock),\n    get_items = get_image_files,\n    get_y = label_func)\n\ndls = dbl.dataloaders(tpath/'pyvips', bs=64)\n\ndls.show_batch(max_n=12)","metadata":{"execution":{"iopub.status.busy":"2022-09-05T03:56:54.778941Z","iopub.execute_input":"2022-09-05T03:56:54.779758Z","iopub.status.idle":"2022-09-05T03:56:57.509881Z","shell.execute_reply.started":"2022-09-05T03:56:54.779706Z","shell.execute_reply":"2022-09-05T03:56:57.508615Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Satisfied!?","metadata":{}},{"cell_type":"markdown","source":"Honestly, not quite there yet!\n\nAlthough it has amazingly improved thanks to `pyvips`, it will still take about a day or so to run it through all the slides. Which I do not feel comfortable to call it practical.\n\nWhat we can do better!?\n\nWell, one solution could be using both `histogram` based operations (e.g. `otsu`) and morphology back to back. You may lose some useful data, but it could shorten the process to the scale of *a couple of hours*.\n\n**I would say it is a trade-off worth trying!**\n","metadata":{}}]}