Use pyramid levels instead of TIFF pages in TiffFileWSIReader

monai/data/wsi_reader.py

@@ -1319,6 +1319,27 @@ class TiffFileWSIReader(BaseWSIReader):
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
 
+    @staticmethod
+    def _get_level_page(wsi, level: int):
+        """
+        Returns the first TIFF page of the given pyramid level.
+
+        Pyramid levels are taken from the first image series, the same levels used by `wsi.asarray(level=...)`.
+        They do not always map one-to-one to TIFF pages, e.g. Aperio SVS files store a thumbnail as the second page.
+
+        Args:
+            wsi: a whole slide image object loaded from a file.
+            level: the level number.
+
+        Returns:
+            the TIFF page holding the metadata of that level.
+
+        Raises:
+            IndexError: when ``level`` is not a valid pyramid level of the image.
+
+        """
+        return wsi.series[0].levels[level].keyframe
+
     @staticmethod
     def get_level_count(wsi) -> int:
         """
@@ -1328,7 +1349,7 @@ def get_level_count(wsi) -> int:
             wsi: a whole slide image object loaded from a file.
 
         """
-        return len(wsi.pages)
+        return len(wsi.series[0].levels)
 
     def get_size(self, wsi, level: int) -> tuple[int, int]:
         """
@@ -1339,7 +1360,8 @@ def get_size(self, wsi, level: int) -> tuple[int, int]:
             level: the level number where the size is calculated.
 
         """
-        return (wsi.pages[level].imagelength, wsi.pages[level].imagewidth)
+        page = self._get_level_page(wsi, level)
+        return (page.imagelength, page.imagewidth)
 
     def get_downsample_ratio(self, wsi, level: int) -> float:
         """
@@ -1350,7 +1372,7 @@ def get_downsample_ratio(self, wsi, level: int) -> float:
             level: the level number where the downsample ratio is calculated.
 
         """
-        return float(wsi.pages[0].imagelength) / float(wsi.pages[level].imagelength)
+        return float(self._get_level_page(wsi, 0).imagelength) / float(self._get_level_page(wsi, level).imagelength)
 
     @staticmethod
     def get_file_path(wsi) -> str:
@@ -1366,13 +1388,14 @@ def get_mpp(self, wsi, level: int) -> tuple[float, float]:
             level: the level number where the mpp is calculated.
 
         """
+        page = self._get_level_page(wsi, level)
         if (
-            "XResolution" in wsi.pages[level].tags
-            and "YResolution" in wsi.pages[level].tags
-            and wsi.pages[level].tags["XResolution"].value
-            and wsi.pages[level].tags["YResolution"].value
+            "XResolution" in page.tags
+            and "YResolution" in page.tags
+            and page.tags["XResolution"].value
+            and page.tags["YResolution"].value
         ):
-            unit = wsi.pages[level].tags.get("ResolutionUnit")
+            unit = page.tags.get("ResolutionUnit")
             if unit is not None:
                 unit = str(unit.value.name)
             if unit is None or len(unit) == 0:
@@ -1382,8 +1405,8 @@ def get_mpp(self, wsi, level: int) -> tuple[float, float]:
             convert_to_micron = ConvertUnits(unit, "micrometer")
 
             # Here, x and y resolutions are rational numbers so each of them is represented by a tuple.
-            yres = wsi.pages[level].tags["YResolution"].value
-            xres = wsi.pages[level].tags["XResolution"].value
+            yres = page.tags["YResolution"].value
+            xres = page.tags["XResolution"].value
             if xres[0] and yres[0]:
                 return convert_to_micron(yres[1] / yres[0]), convert_to_micron(xres[1] / xres[0])
             else:
@@ -1417,16 +1440,16 @@ def get_wsi_at_mpp(
         """
         mpp = ensure_tuple_rep(mpp, 2)
 
-        mpp_list = [self.get_mpp(wsi, lvl) for lvl in range(len(wsi.pages))]  # Fails for some Tifffiles
+        mpp_list = [self.get_mpp(wsi, lvl) for lvl in range(self.get_level_count(wsi))]  # Fails for some Tifffiles
         closest_lvl = self._find_closest_level("mpp", mpp, mpp_list, 0, 5)
 
         within_tolerance, closest_level_is_bigger = self._compute_mpp_tolerances(closest_lvl, mpp_list, mpp, atol, rtol)
 
         if within_tolerance:
             # If the image at the desired mpp resolution is within tolerances, return the image at closest_level.
-            # TiffFile does not expose `read_region`; read the whole page instead (consistent with `_get_patch`).
+            # TiffFile does not expose `read_region`; read the whole level instead (consistent with `_get_patch`).
             pil_image, _ = optional_import("PIL", name="Image")
-            closest_lvl_wsi = pil_image.fromarray(wsi.pages[closest_lvl].asarray())
+            closest_lvl_wsi = pil_image.fromarray(wsi.asarray(level=closest_lvl))
 
         elif closest_level_is_bigger:
             # Otherwise, select the level closest to the desired mpp with a higher resolution and downsample it.
@@ -1563,7 +1586,7 @@ def _resize_to_mpp_res(self, wsi, closest_lvl, mpp_list, user_mpp: tuple):
 
         target_res_x, target_res_y = self._compute_mpp_target_res(closest_lvl, closest_lvl_dim, mpp_list, user_mpp)
 
-        closest_lvl_wsi = pil_image.fromarray(wsi.pages[closest_lvl].asarray())
+        closest_lvl_wsi = pil_image.fromarray(wsi.asarray(level=closest_lvl))
         # PIL `resize` takes (W, H); use the same row/col order as OpenSlide.
         closest_lvl_wsi = closest_lvl_wsi.resize((target_res_y, target_res_x), pil_image.BILINEAR)
 

tests/utils/enums/test_wsireader.py

@@ -40,6 +40,7 @@
 has_cucim = has_cucim and hasattr(cucim, "CuImage")
 openslide, has_osl = optional_import("openslide")
 imwrite, has_tiff = optional_import("tifffile", name="imwrite")
+TiffWriter, _ = optional_import("tifffile", name="TiffWriter")
 _, has_codec = optional_import("imagecodecs")
 has_tiff = has_tiff and has_codec
 
@@ -60,6 +61,8 @@
 
 TEST_CASE_WHOLE_0 = [WSI_GENERIC_TIFF_PATH, 2, (3, WSI_GENERIC_TIFF_HEIGHT // 4, WSI_GENERIC_TIFF_WIDTH // 4)]
 
+TEST_CASE_WHOLE_1 = [WSI_APERIO_SVS_PATH, 2, (3, 2057, 2875)]
+
 TEST_CASE_TRANSFORM_0 = [
     WSI_GENERIC_TIFF_PATH,
     4,
@@ -318,6 +321,18 @@
 
 TEST_CASE_SVS_MPP_SCALAR = [WSI_APERIO_SVS_PATH, {"mpp": 8.0}, {"openslide": (3, 2057, 2875), "cucim": (3, 2057, 2875)}]
 
+# pyramid levels of Aperio SVS (the thumbnail, label, and macro images are not levels)
+TEST_CASE_SVS_LEVELS = [WSI_APERIO_SVS_PATH, [(32914, 46000), (8228, 11500), (2057, 2875)]]
+
+# patch on tissue at the lowest resolution level of Aperio SVS
+# (the location maps to about 0.3 pixel past a level-2 pixel, and the surrounding pixels differ by at most 5,
+# so backends that round, truncate, or interpolate the level-2 location give nearly the same values)
+TEST_CASE_SVS_LEVEL_2 = [
+    WSI_APERIO_SVS_PATH,
+    {"location": (16261, 31782), "size": (2, 1), "level": 2},
+    np.array([[[106], [105]], [[64], [66]], [[115], [116]]], dtype=np.uint8),
+]
+
 TEST_CASE_DEVICE_2 = [
     WSI_GENERIC_TIFF_PATH,
     {"level": 8, "dtype": torch.float32, "device": "cuda"},
@@ -449,6 +464,33 @@ def save_gray_tiff(array: np.ndarray, filename: str):
     return filename
 
 
+def save_pyramid_tiff_with_thumbnail(array: np.ndarray, filename: str):
+    """
+    Save a two-level RGB pyramid with a thumbnail page between the levels, as in Aperio SVS files.
+    Level 0 has 1 micron per pixel and level 1 has 2 microns per pixel.
+
+    Args:
+        array: numpy ndarray with the shape of HxWx3, the full resolution image
+        filename: the filename to be used for the tiff file.
+
+    Returns:
+        the filename of the saved tiff file.
+    """
+    with TiffWriter(filename) as tif:
+        tif.write(array, photometric="rgb", tile=(64, 64), resolution=(1e4, 1e4), resolutionunit="CENTIMETER")
+        tif.write(array[::5, ::5], photometric="rgb")
+        tif.write(
+            array[::2, ::2],
+            photometric="rgb",
+            tile=(64, 64),
+            subfiletype=1,
+            resolution=(5e3, 5e3),
+            resolutionunit="CENTIMETER",
+        )
+
+    return filename
+
+
 @skipUnless(has_cucim or has_osl or has_tiff, "Requires cucim, openslide, or tifffile!")
 def setUpModule():
     download_url_or_skip_test(
@@ -490,7 +532,7 @@ def tearDown(self):
             """
             gc.collect()
 
-        @parameterized.expand([TEST_CASE_WHOLE_0])
+        @parameterized.expand([TEST_CASE_WHOLE_0, TEST_CASE_WHOLE_1])
         def test_read_whole_image(self, file_path, level, expected_shape):
             reader = WSIReader(self.backend, level=level)
             with reader.read(file_path) as img_obj:
@@ -699,6 +741,40 @@ def test_read_region_multi_thread(self, file_path, kwargs, patch_info, expected_
                     assert_allclose(meta[WSIPatchKeys.SIZE], patch_info["size"], type_test=False)
                     assert_allclose(meta[WSIPatchKeys.LOCATION], patch_info["location"], type_test=False)
 
+        @parameterized.expand([TEST_CASE_SVS_LEVELS])
+        def test_svs_levels(self, file_path, expected_sizes):
+            """Verify the level count, level sizes, and downsample ratios of an Aperio SVS pyramid.
+
+            Args:
+                file_path: Path to the Aperio SVS file.
+                expected_sizes: Expected (height, width) of each pyramid level, from level 0 down.
+            """
+            reader = WSIReader(self.backend)
+            with reader.read(file_path) as img_obj:
+                self.assertEqual(reader.get_level_count(img_obj), len(expected_sizes))
+                for level, expected_size in enumerate(expected_sizes):
+                    self.assertTupleEqual(tuple(reader.get_size(img_obj, level)), expected_size)
+                    self.assertAlmostEqual(
+                        reader.get_downsample_ratio(img_obj, level), expected_sizes[0][0] / expected_size[0], places=2
+                    )
+
+        @parameterized.expand([TEST_CASE_SVS_LEVEL_2])
+        def test_read_region_svs_level(self, file_path, patch_info, expected_img):
+            """Verify that a patch read at a reduced level of an Aperio SVS file comes from the requested location.
+
+            Args:
+                file_path: Path to the Aperio SVS file.
+                patch_info: Keyword arguments for ``get_data``: level-0 location, patch size, and level.
+                expected_img: Expected patch in channel-first layout.
+            """
+            reader = WSIReader(self.backend)
+            with reader.read(file_path) as img_obj:
+                img, meta = reader.get_data(img_obj, **patch_info)
+            self.assertTupleEqual(img.shape, expected_img.shape)
+            # backends map a level-0 location to a level-2 pixel slightly differently (see TEST_CASE_SVS_LEVEL_2)
+            assert_allclose(img, expected_img, type_test=False, atol=5)
+            self.assertEqual(meta[WSIPatchKeys.LEVEL], patch_info["level"])
+
         @parameterized.expand([TEST_CASE_MPP_0])
         def test_resolution_mpp(self, file_path, level, expected_mpp):
             reader = WSIReader(self.backend, level=level)
@@ -772,6 +848,19 @@ class TestTiffFile(WSIReaderTests.Tests):
     def setUpClass(cls):
         cls.backend = "tifffile"
 
+    def test_mpp_levels_with_thumbnail_page(self):
+        """Verify that level metadata and MPP reads skip a thumbnail page stored between pyramid levels."""
+        image = np.random.RandomState(0).randint(0, 256, (512, 512, 3), dtype=np.uint8)
+        file_path = save_pyramid_tiff_with_thumbnail(
+            image, os.path.join(TESTS_PATH, "testing_data", "temp_tiff_pyramid_thumbnail.tiff")
+        )
+        reader = WSIReader(self.backend)
+        with reader.read(file_path) as img_obj:
+            self.assertEqual(reader.get_level_count(img_obj), 2)
+            self.assertTupleEqual(tuple(reader.get_mpp(img_obj, 1)), (2.0, 2.0))
+            img = reader.get_wsi_at_mpp(img_obj, mpp=2.0)
+        assert_allclose(img, np.moveaxis(image[::2, ::2], -1, 0), type_test=False)
+
 
 if __name__ == "__main__":
     unittest.main()