🆕 Improve Support for QuPath GeoJson and Multipoint Geometry (#841)

This PR makes a collection of changes related to reviewer comments in the TIAViz paper.

Main changes are:

1. Layer names in the interface for static image overlays now are named after filename stem, instead of generic 'layer1', 'layer2' etc to make it easier to distinguish the layers.
2. Adds multipoint geometry rendering as qupath exports collections of points as multipoints not individual point annotations
3. Robustifies qupath import from geojson to handle qupath special case of `nucleusGeometry`.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: measty

tests/test_annotation_tilerendering.py

@@ -270,15 +270,37 @@ def sub_tile_level_count(self: MockTileGenerator) -> int:
     assert tile.size == (112, 112)
 
 
+def test_multi_point() -> None:
+    """Test multi-point rendering."""
+    renderer = AnnotationRenderer(max_scale=8, edge_thickness=0)
+    tile = np.zeros((256, 256, 3), dtype=np.uint8)
+    renderer.render_by_type(
+        tile=tile,
+        annotation=Annotation(MultiPoint([(5.0, 5.0), (10.0, 10.0)])),
+        top_left=(0, 0),
+        scale=1,
+    )
+    # check point locations are now non-zero
+    assert np.any(tile[5, 5, :] > 0)
+    assert np.any(tile[10, 10, :] > 0)
+    # check a non-point region still zero
+    assert np.all(tile[100:150, 100:150, :] == 0)
+
+
 def test_unknown_geometry(
-    fill_store: Callable,  # noqa: ARG001
     caplog: pytest.LogCaptureFixture,
 ) -> None:
     """Test warning when unknown geometries cannot be rendered."""
     renderer = AnnotationRenderer(max_scale=8, edge_thickness=0)
+    unknown_ann = Annotation(
+        MultiPolygon(
+            [Polygon.from_bounds(0, 0, 1, 1), Polygon.from_bounds(2, 2, 3, 3)]
+        ),
+        {},
+    )
     renderer.render_by_type(
         tile=np.zeros((256, 256, 3), dtype=np.uint8),
-        annotation=Annotation(MultiPoint([(5.0, 5.0), (10.0, 10.0)])),
+        annotation=unknown_ann,
         top_left=(0, 0),
         scale=1,
     )

tests/test_app_bokeh.py

@@ -285,7 +285,7 @@ def test_add_annotation_layer(doc: Document, data_path: pytest.TempPathFactory)
     # trigger an event to select the geojson file
     click = MenuItemClick(layer_drop, str(data_path["geojson_anns"]))
     layer_drop._trigger_event(click)
-    assert main.UI["vstate"].types == ["annotation"]
+    assert set(main.UI["vstate"].types) == {"nucleus", "cell", "annotation"}
 
     # test the name2type function.
     assert main.name2type("annotation") == '"annotation"'
@@ -434,7 +434,8 @@ def test_load_img_overlay(doc: Document, data_path: pytest.TempPathFactory) -> N
     # trigger an event to select the image overlay
     click = MenuItemClick(layer_drop, str(data_path["img_overlay"]))
     layer_drop._trigger_event(click)
-    layer_slider = doc.get_model_by_name("layer2_slider")
+    l_name = data_path["img_overlay"].stem
+    layer_slider = doc.get_model_by_name(f"{l_name}_slider")
     assert layer_slider is not None
 
     # check alpha controls
@@ -443,14 +444,22 @@ def test_load_img_overlay(doc: Document, data_path: pytest.TempPathFactory) -> N
     assert type_column_list[-1].active
     # toggle off and check alpha is 0
     type_column_list[-1].active = False
-    assert main.UI["p"].renderers[main.UI["vstate"].layer_dict["layer2"]].alpha == 0
+    assert main.UI["p"].renderers[main.UI["vstate"].layer_dict[l_name]].alpha == 0
     # toggle back on and check alpha is back to default 0.75
     type_column_list[-1].active = True
-    assert main.UI["p"].renderers[main.UI["vstate"].layer_dict["layer2"]].alpha == 0.75
+    assert main.UI["p"].renderers[main.UI["vstate"].layer_dict[l_name]].alpha == 0.75
     # set alpha to 0.4
     layer_slider.value = 0.4
     # check that the alpha values have been set correctly
-    assert main.UI["p"].renderers[main.UI["vstate"].layer_dict["layer2"]].alpha == 0.4
+    assert main.UI["p"].renderers[main.UI["vstate"].layer_dict[l_name]].alpha == 0.4
+
+    # check loading a new layer with same stem uses full file name to disambiguate
+    click = MenuItemClick(layer_drop, str(data_path["annotations"]))
+    layer_drop._trigger_event(click)
+    click = MenuItemClick(layer_drop, str(data_path["img_overlay"]))
+    layer_drop._trigger_event(click)
+    full_name = data_path["img_overlay"].name
+    assert full_name in main.UI["vstate"].layer_dict
 
 
 def test_hovernet_on_box(doc: Document, data_path: pytest.TempPathFactory) -> None:

tests/test_tileserver.py

@@ -100,11 +100,7 @@ def app(remote_sample: Callable, tmp_path: Path) -> TileServer:
     imwrite(thumb_path, thumb)
 
     sample_store = Path(remote_sample("annotation_store_svs_1"))
-    store = SQLiteStore(sample_store)
-    geo_path = tmp_path / "test.geojson"
-    store.to_geojson(geo_path)
-    store.commit()
-    store.close()
+    geo_path = Path(remote_sample("geojson_cmu_1"))
 
     # make tileserver with layers representing all the types
     # of things it should be able to handle
@@ -115,7 +111,7 @@ def app(remote_sample: Callable, tmp_path: Path) -> TileServer:
             "tile": str(thumb_path),
             "im_array": np.zeros(wsi.slide_dimensions(1.25, "power"), dtype=np.uint8).T,
             "overlay": str(sample_store),
-            "store_geojson": tmp_path / "test.geojson",
+            "store_geojson": str(geo_path),
         },
     )
     app.config.from_mapping({"TESTING": True})
@@ -446,14 +442,22 @@ def test_change_overlay(  # noqa: PLR0915
         assert response.status_code == 200
         response = client.put(
             "/tileserver/overlay",
-            data={"overlay_path": safe_str(geo_path)},
+            data={"overlay_path": safe_str(sample_store)},
         )
         assert response.status_code == 200
         assert response.content_type == "text/html; charset=utf-8"
         assert set(json.loads(response.data)) == {0, 1, 2, 3, 4}
         # check that the annotations have been correctly loaded
         assert len(empty_app.pyramids[session_id]["overlay"].store) == num_annotations
 
+        # check loading a geojson from qupath with cells/nuclei
+        qpath_geo_path = Path(remote_sample("geojson_cmu_1"))
+        response = client.put(
+            "/tileserver/overlay",
+            data={"overlay_path": safe_str(qpath_geo_path)},
+        )
+        assert response.status_code == 200
+
         # add another image layer
         response = client.put(
             "/tileserver/overlay",
@@ -462,10 +466,18 @@ def test_change_overlay(  # noqa: PLR0915
         assert response.status_code == 200
         assert response.content_type == "text/html; charset=utf-8"
         # check that the overlay has been correctly added
-        lname = f"layer{len(empty_app.pyramids[session_id]) - 1}"
+        lname = Path(overlay_path).stem
         layer = empty_app.pyramids[session_id][lname]
         assert layer.wsi.info.file_path == overlay_path
 
+        # add same image again to check if duplicate names triggers disambiguation
+        response = client.put(
+            "/tileserver/overlay",
+            data={"overlay_path": safe_str(overlay_path)},
+        )
+        assert response.status_code == 200
+        assert Path(overlay_path).name in empty_app.pyramids[session_id]
+
         # replace existing store overlay
         response = client.put(
             "/tileserver/overlay",
@@ -494,7 +506,7 @@ def test_change_overlay(  # noqa: PLR0915
             data={"overlay_path": safe_str(jpg_path)},
         )
         # check that the overlay has been correctly added
-        lname = f"layer{len(empty_app.pyramids[session_id]) - 1}"
+        lname = Path(jpg_path).stem
         layer = empty_app.pyramids[session_id][lname]
         assert np.all(layer.wsi.img == imread(jpg_path))
 
@@ -518,7 +530,7 @@ def test_change_overlay(  # noqa: PLR0915
             data={"overlay_path": safe_str(tiff_path)},
         )
         # check that the overlay has been correctly added
-        lname = f"layer{len(empty_app.pyramids[session_id]) - 1}"
+        lname = Path(tiff_path).stem
         layer = empty_app.pyramids[session_id][lname]
         assert layer.wsi.info.file_path == tiff_path
 

tiatoolbox/annotation/storage.py

@@ -1995,6 +1995,20 @@ def transform_geometry(geom: Geometry) -> Geometry:
             )
             for feature in geojson["features"]
         ]
+        # check for presence of 'nucleusGeometry' key in features
+        # if present, add them (support qupath export format)
+        annotations += [
+            transform(
+                Annotation(
+                    transform_geometry(
+                        feature2geometry(feature["nucleusGeometry"]),
+                    ),
+                    {"type": "nucleus"},
+                ),
+            )
+            for feature in geojson["features"]
+            if "nucleusGeometry" in feature
+        ]
 
         logger.info("Adding %d annotations.", len(annotations))
         self.append_many(annotations)

tiatoolbox/utils/visualization.py

@@ -963,11 +963,41 @@ def render_pt(
                 top_left,
                 scale,
             )[0][0],
-            np.maximum(self.edge_thickness, 1),
+            np.maximum(int(16 / scale**0.5), 1),
             col,
-            thickness=self.thickness,
+            thickness=-1,
         )
 
+    def render_pts(
+        self: AnnotationRenderer,
+        tile: np.ndarray,
+        annotation: Annotation,
+        top_left: tuple[float, float],
+        scale: float,
+    ) -> None:
+        """Render a multipoint annotation onto a tile using cv2.
+
+        Args:
+            tile (ndarray):
+                The rgb(a) tile image to render onto.
+            annotation (Annotation):
+                The annotation to render.
+            top_left (tuple):
+                The top left corner of the tile in wsi.
+            scale (float):
+                The zoom scale at which we are rendering.
+
+        """
+        col = self.get_color(annotation, edge=False)
+        for pt in annotation.coords:
+            cv2.circle(
+                tile,
+                self.to_tile_coords([pt], top_left, scale)[0][0],
+                np.maximum(int(16 / scale**0.5), 1),
+                col,
+                thickness=-1,
+            )
+
     def render_line(
         self: AnnotationRenderer,
         tile: np.ndarray,
@@ -1172,5 +1202,7 @@ def render_by_type(
             self.render_poly(tile, annotation, top_left, scale)
         elif geom_type == GeometryType.LINE_STRING:
             self.render_line(tile, annotation, top_left, scale)
+        elif geom_type == GeometryType.MULTI_POINT:
+            self.render_pts(tile, annotation, top_left, scale)
         else:
             logger.warning("Unknown geometry: %s", geom_type, stacklevel=3)

tiatoolbox/visualization/bokeh_app/main.py

@@ -510,7 +510,6 @@ def add_layer(lname: str) -> None:
             end=1,
             value=0.75,
             step=0.01,
-            title=lname,
             height=40,
             width=100,
             max_width=90,
@@ -1069,7 +1068,9 @@ def layer_slider_cb(
             UI["vstate"].layer_dict[obj.name.split("_")[0]]
         ].glyph.line_alpha = new
     else:
-        UI["p"].renderers[UI["vstate"].layer_dict[obj.name.split("_")[0]]].alpha = new
+        UI["p"].renderers[
+            UI["vstate"].layer_dict["_".join(obj.name.split("_")[0:-1])]
+        ].alpha = new
 
 
 def color_input_cb(

tiatoolbox/visualization/tileserver.py

@@ -35,6 +35,7 @@
 if TYPE_CHECKING:  # pragma: no cover
     from matplotlib.colors import Colormap
 
+    from tiatoolbox.annotation.storage import Annotation
     from tiatoolbox.wsicore import WSIMeta
 
 
@@ -538,7 +539,7 @@ def _handle_registration_overlay(
         self,
         session_id: str,
         overlay_path: Path,
-        other_session_id: str,
+        other_session_id: str | None,
     ) -> str:
         def _apply_transform(source_fp: str, target_fp: str) -> None:
             # loading a registration transformation
@@ -559,10 +560,11 @@ def _apply_transform(source_fp: str, target_fp: str) -> None:
             )
             return json.dumps("slide")
 
-        layer_keys = [k for k in self.layers[session_id] if "layer" in k]
-        for key in sorted(
-            layer_keys, key=lambda x: int(x.replace("layer", "")), reverse=True
-        ):
+        layer_keys = [
+            k for k in self.layers[session_id] if k not in ["slide", "overlay"]
+        ]
+        layer_keys.reverse()  # try newest first
+        for key in layer_keys:
             target_fp = self.layers[session_id][key].info.file_path
             if Path(target_fp).suffix in [".tiff", ".svs", ".ndpi", ".mrxs"]:
                 logger.warning(
@@ -585,7 +587,10 @@ def _apply_transform(source_fp: str, target_fp: str) -> None:
         return json.dumps("slide")
 
     def _add_image_overlay(self, session_id: str, overlay_path: Path) -> str:
-        layer = f"layer{len(self.pyramids[session_id])}"
+        layer = overlay_path.stem
+        if layer in self.layers[session_id]:
+            # use full file name to disambiguate
+            layer = overlay_path.name
         if overlay_path.suffix == ".tiff":
             self.layers[session_id][layer] = OpenSlideWSIReader(
                 overlay_path,
@@ -608,7 +613,19 @@ def _add_image_overlay(self, session_id: str, overlay_path: Path) -> str:
 
     def _add_annotation_overlay(self, session_id: str, overlay_path: Path) -> str:
         if overlay_path.suffix == ".geojson":
-            sq = SQLiteStore.from_geojson(overlay_path)
+
+            def unpack_qupath(ann: Annotation) -> Annotation:
+                # Helper function to unpack QuPath measurements if present.
+                props = ann.properties
+                if "measurements" in props:
+                    measurements = props.pop("measurements")
+                    for k, v in measurements.items():
+                        props[k] = v
+                if "objectType" in props:
+                    props["type"] = props.pop("objectType")
+                return ann
+
+            sq = SQLiteStore.from_geojson(overlay_path, transform=unpack_qupath)
         elif overlay_path.suffix == ".dat":
             sq = store_from_dat(overlay_path)
         if overlay_path.suffix == ".db":