ENH: Swap pygeos to shapely (#202)

- Replace pygeos with shapely to resolve the GEO binary installation issue.

Co-authored-by: @vqdang

Author: vqdang

docs/installation.rst

@@ -29,19 +29,19 @@ Windows
 1. Download OpenSlide binaries from `this page <https://openslide.org/download/>`_. Extract the folder and add `bin` and `lib` subdirectories to
 Windows `system path <https://docs.microsoft.com/en-us/previous-versions/office/developer/sharepoint-2010/ee537574(v=office.14)>`_.
 
-2. Install
-TIAToolbox.
+2. Install OpenJPEG. The easiest way is to install OpenJpeg is through conda
+using
 
 .. code-block:: console
 
-    $ pip install tiatoolbox
+    C:\> conda install -c conda-forge openjpeg
 
-3. Install OpenJPEG. The easiest way is to install OpenJpeg is through conda
-using
+3. Install
+TIAToolbox.
 
 .. code-block:: console
 
-    C:\> conda install -c conda-forge openjpeg pygeos
+    C:\> pip install tiatoolbox
 
 Linux (Ubuntu)
 ^^^^^^^^^^^^^^

requirements.conda.yml

@@ -23,7 +23,6 @@ dependencies:
   - pip
   - pixman<0.38.0
   - python=3.7
-  - pygeos
   - python=3.7
   - pytorch
   - pyyaml>=5.1

requirements.dev.conda.yml

@@ -26,7 +26,6 @@ dependencies:
   - pip
   - pixman<0.38.0
   - pre-commit
-  - pygeos
   - pytest-cov
   - pytest-runner==5.2
   - pytest==6.2.5

requirements.txt

@@ -10,7 +10,6 @@ opencv-python>=4.0
 openslide-python==1.1.2
 pandas
 pillow
-pygeos
 pyyaml>=5.1
 requests
 scikit-image

requirements.win64.conda.yml

@@ -21,7 +21,6 @@ dependencies:
   - pillow
   - pip
   - pixman<0.38.0
-  - pygeos
   - python=3.7
   - pytorch
   - pyyaml>=5.1

requirements_dev.txt

@@ -16,7 +16,6 @@ pandas
 pillow
 pip>=20.0.2
 pre-commit
-pygeos
 pytest-cov==2.9.0
 pytest-runner==5.2
 pytest==6.2.5

setup.py

@@ -23,6 +23,7 @@
     "glymur",
     "scikit-learn>=0.23.2",
     "scikit-image>=0.17",
+    "shapely",
     "torchvision==0.10.1",
     "torch==1.9.1",
     "tqdm==4.60.0",

tests/models/test_nucleus_instance_segmentor.py

@@ -20,6 +20,9 @@
 """Tests for Nucleus Instance Segmentor."""
 
 import copy
+
+# ! The garbage collector
+import gc
 import pathlib
 import shutil
 
@@ -298,6 +301,7 @@ def test_crash_segmentor(remote_sample, tmp_path):
 
 def test_functionality_travis(remote_sample, tmp_path):
     """Functionality test for nuclei instance segmentor."""
+    gc.collect()
     root_save_dir = pathlib.Path(tmp_path)
     mini_wsi_svs = pathlib.Path(remote_sample("wsi4_512_512_svs"))
 
@@ -319,7 +323,7 @@ def test_functionality_travis(remote_sample, tmp_path):
             {"units": "mpp", "resolution": resolution},
         ],
         margin=128,
-        tile_shape=[512, 512],
+        tile_shape=[1024, 1024],
         patch_input_shape=[256, 256],
         patch_output_shape=[164, 164],
         stride_shape=[164, 164],

tiatoolbox/models/engine/nucleus_instance_segmentor.py

@@ -26,9 +26,10 @@
 # replace with the sql database once the PR in place
 import joblib
 import numpy as np
-import pygeos
 import torch
 import tqdm
+from shapely.geometry import box as shapely_box
+from shapely.strtree import STRtree
 
 from tiatoolbox.models.engine.semantic_segmentor import (
     IOSegmentorConfig,
@@ -134,32 +135,31 @@ def _process_tile_predictions(
     if len(inst_dict) == 0:
         return {}, []
 
-    # ! DEPRECATION:
-    # !     will be deprecated upon finalization of SQL annotation store
+    # !
     m = ioconfig.margin
     w, h = tile_shape
     inst_boxes = [v["box"] for v in inst_dict.values()]
     inst_boxes = np.array(inst_boxes)
-    tile_rtree = pygeos.STRtree(
-        pygeos.box(
-            inst_boxes[:, 0], inst_boxes[:, 1], inst_boxes[:, 2], inst_boxes[:, 3]
-        )
-    )
+
+    geometries = [shapely_box(*bounds) for bounds in inst_boxes]
+    # An auxiliary dictionary to actually query the index within the source list
+    index_by_id = {id(geo): idx for idx, geo in enumerate(geometries)}
+    tile_rtree = STRtree(geometries)
     # !
 
     # create margin bounding box, ordering should match with
     # created tile info flag (top, bottom, left, right)
     boundary_lines = [
-        pygeos.box(0, 0, w, 1),  # noqa top egde
-        pygeos.box(0, h - 1, w, h),  # noqa bottom edge
-        pygeos.box(0, 0, 1, h),  # noqa left
-        pygeos.box(w - 1, 0, w, h),  # noqa right
+        shapely_box(0, 0, w, 1),  # noqa top egde
+        shapely_box(0, h - 1, w, h),  # noqa bottom edge
+        shapely_box(0, 0, 1, h),  # noqa left
+        shapely_box(w - 1, 0, w, h),  # noqa right
     ]
     margin_boxes = [
-        pygeos.box(0, 0, w, m),  # noqa top egde
-        pygeos.box(0, h - m, w, h),  # noqa bottom edge
-        pygeos.box(0, 0, m, h),  # noqa left
-        pygeos.box(w - m, 0, w, h),  # noqa right
+        shapely_box(0, 0, w, m),  # noqa top egde
+        shapely_box(0, h - m, w, h),  # noqa bottom edge
+        shapely_box(0, 0, m, h),  # noqa left
+        shapely_box(w - m, 0, w, h),  # noqa right
     ]
     # ! this is wrt to WSI coord space, not tile
     margin_lines = [
@@ -169,7 +169,7 @@ def _process_tile_predictions(
         [[w - m, m], [w - m, h - m]],  # noqa right
     ]
     margin_lines = np.array(margin_lines) + tile_tl[None, None]
-    margin_lines = [pygeos.box(*v.flatten().tolist()) for v in margin_lines]
+    margin_lines = [shapely_box(*v.flatten().tolist()) for v in margin_lines]
 
     # the ids within this match with those within `inst_map`, not UUID
     sel_indices = []
@@ -182,8 +182,12 @@ def _process_tile_predictions(
             for idx, box in enumerate(margin_boxes)
             if tile_flag[idx] or tile_mode == 3
         ]
+
         sel_indices = [
-            tile_rtree.query(bounds, predicate="contains") for bounds in sel_boxes
+            index_by_id[id(geo)]
+            for bounds in sel_boxes
+            for geo in tile_rtree.query(bounds)
+            if bounds.contains(geo)
         ]
     elif tile_mode in [1, 2]:
         # for `horizontal/vertical strip` tiles
@@ -196,8 +200,11 @@ def _process_tile_predictions(
             margin_boxes[idx] if flag else boundary_lines[idx]
             for idx, flag in enumerate(tile_flag)
         ]
+
         sel_indices = [
-            tile_rtree.query(bounds, predicate="intersects") for bounds in sel_boxes
+            index_by_id[id(geo)]
+            for bounds in sel_boxes
+            for geo in tile_rtree.query(bounds)
         ]
     else:
         raise ValueError(f"Unknown tile mode {tile_mode}.")
@@ -208,7 +215,6 @@ def retrieve_sel_uids(sel_indices, inst_dict):
         if len(sel_indices) > 0:
             # not sure how costly this is in large dict
             inst_uids = list(inst_dict.keys())
-            sel_indices = [idx for sub_sel in sel_indices for idx in sub_sel]
             sel_uids = [inst_uids[idx] for idx in sel_indices]
         return sel_uids
 
@@ -218,26 +224,19 @@ def retrieve_sel_uids(sel_indices, inst_dict):
     # this one should contain UUID with the reference database
     remove_insts_in_orig = []
     if tile_mode == 3:
-        # ! DEPRECATION:
-        # !     will be deprecated upon finalization of SQL annotation store
         inst_boxes = [v["box"] for v in ref_inst_dict.values()]
         inst_boxes = np.array(inst_boxes)
-        ref_inst_rtree = pygeos.STRtree(
-            pygeos.box(
-                inst_boxes[:, 0],
-                inst_boxes[:, 1],
-                inst_boxes[:, 2],
-                inst_boxes[:, 3],
-            )
-        )
-        # !
 
-        # remove existing instances in old prediction which intersect
-        # with the margin lines
+        geometries = [shapely_box(*bounds) for bounds in inst_boxes]
+        # An auxiliary dictionary to actually query the index within the source list
+        index_by_id = {id(geo): idx for idx, geo in enumerate(geometries)}
+        ref_inst_rtree = STRtree(geometries)
         sel_indices = [
-            ref_inst_rtree.query(bounds, predicate="intersects")
+            index_by_id[id(geo)]
             for bounds in margin_lines
+            for geo in ref_inst_rtree.query(bounds)
         ]
+
         remove_insts_in_orig = retrieve_sel_uids(sel_indices, ref_inst_dict)
 
     # move inst position from tile space back to WSI space
@@ -390,24 +389,30 @@ def _get_tile_info(
         # * === be removed in postproc callback
         boxes = tile_outputs
 
+        # This saves computation time if the image is smaller than the expected tile
+        if np.all(image_shape <= tile_shape):
+            flag = np.zeros([boxes.shape[0], 4], dtype=np.int32)
+            return [[boxes, flag]]
+
         # * remove all sides for boxes
         # unset for those lie within the selection
         def unset_removal_flag(boxes, removal_flag):
             """Unset removal flags for tiles intersecting image boundaries."""
-            # ! DEPRECATION:
-            # !     will be deprecated upon finalization of SQL annotation store
             sel_boxes = [
-                pygeos.box(0, 0, w, 0),  # top edge
-                pygeos.box(0, h, w, h),  # bottom edge
-                pygeos.box(0, 0, 0, h),  # left
-                pygeos.box(w, 0, w, h),  # right
+                shapely_box(0, 0, w, 0),  # top edge
+                shapely_box(0, h, w, h),  # bottom edge
+                shapely_box(0, 0, 0, h),  # left
+                shapely_box(w, 0, w, h),  # right
             ]
-            spatial_indexer = pygeos.STRtree(
-                pygeos.box(boxes[:, 0], boxes[:, 1], boxes[:, 2], boxes[:, 3])
-            )
-            # !
+            geometries = [shapely_box(*bounds) for bounds in boxes]
+            # An auxiliary dictionary to actually query the index within the source list
+            index_by_id = {id(geo): idx for idx, geo in enumerate(geometries)}
+            spatial_indexer = STRtree(geometries)
+
             for idx, sel_box in enumerate(sel_boxes):
-                sel_indices = spatial_indexer.query(sel_box)
+                sel_indices = [
+                    index_by_id[id(geo)] for geo in spatial_indexer.query(sel_box)
+                ]
                 removal_flag[sel_indices, idx] = 0
             return removal_flag
 
@@ -621,17 +626,10 @@ def _predict_one_wsi(
             patch_inputs = patch_inputs[sel]
 
         # assume to be in [top_left_x, top_left_y, bot_right_x, bot_right_y]
-        # ! DEPRECATION:
-        # !     will be deprecated upon finalization of SQL annotation store
-        spatial_indexer = pygeos.STRtree(
-            pygeos.box(
-                patch_outputs[:, 0],
-                patch_outputs[:, 1],
-                patch_outputs[:, 2],
-                patch_outputs[:, 3],
-            )
-        )
-        # !
+        geometries = [shapely_box(*bounds) for bounds in patch_outputs]
+        # An auxiliary dictionary to actually query the index within the source list
+        index_by_id = {id(geo): idx for idx, geo in enumerate(geometries)}
+        spatial_indexer = STRtree(geometries)
 
         # * retrieve tile placement and tile info flag
         # tile shape will always be corrected to be multiple of output
@@ -651,7 +649,11 @@ def _predict_one_wsi(
 
                 # select any patches that have their output
                 # within the current tile
-                sel_indices = spatial_indexer.query(pygeos.box(*tile_bounds))
+                sel_box = shapely_box(*tile_bounds)
+                sel_indices = [
+                    index_by_id[id(geo)] for geo in spatial_indexer.query(sel_box)
+                ]
+
                 tile_patch_inputs = patch_inputs[sel_indices]
                 tile_patch_outputs = patch_outputs[sel_indices]
                 self._to_shared_space(wsi_idx, tile_patch_inputs, tile_patch_outputs)