🐛 Fix DICOM Reading (Issue #933)

tests/conftest.py

@@ -162,6 +162,17 @@ def sample_jp2(remote_sample: Callable) -> Path:
     return remote_sample("jp2-omnyx-small")
 
 
+@pytest.fixture(scope="session")
+def sample_dicom(remote_sample: Callable) -> Path:
+    """Sample pytest fixture for DICOM images.
+
+    This fixture downloads a sample DICOM file in a standard format for testing.
+    The file represents a single DICOM image and is stored in a temporary directory.
+
+    """
+    return remote_sample("dicom-1")
+
+
 @pytest.fixture(scope="session")
 def sample_all_wsis(
     sample_ndpi: Path,

tests/test_wsireader.py

@@ -85,6 +85,41 @@
 # -------------------------------------------------------------------------------------
 # Utility Test Functions
 # -------------------------------------------------------------------------------------
+def get_tissue_com_tile(reader: WSIReader, size: int) -> IntBounds:
+    """Returns bounds of a tile located approximately at COM of the tissue.
+
+    Uses reader.tissue_mask() to find the center of mass of the tissue
+    and returns a tile centered at that point, at requested size. Used
+    to ensure we are looking at a tissue region when doing level consistency
+     tests etc.
+
+    Args:
+        reader (WSIReader): WSIReader instance.
+        size (int): Size at baseline of the tile to return.
+
+    Returns:
+        IntBounds: Baseline bounds of the tile centered at COM of the tissue.
+
+    """
+    mask = reader.tissue_mask(resolution=8.0, units="mpp").img
+
+    # Find the center of mass of the tissue
+    ys, xs = np.nonzero(mask)
+    com_y = int(ys.mean())
+    com_x = int(xs.mean())
+    # convert to baseline coordinates
+    com_x = int(com_x * (8.0 / reader.info.mpp[0]))
+    com_y = int(com_y * (8.0 / reader.info.mpp[1]))
+
+    # Calculate bounds for the tile centered at COM
+    half_size = size // 2
+    bounds = (
+        max(0, com_x - half_size),
+        max(0, com_y - half_size),
+        min(reader.info.slide_dimensions[0], com_x + half_size),
+        min(reader.info.slide_dimensions[1], com_y + half_size),
+    )
+    return np.array(bounds)
 
 
 def strictly_increasing(sequence: Iterable) -> bool:
@@ -2690,7 +2725,8 @@ def test_read_rect_level_consistency(wsi: WSIReader) -> None:
     they are aligned.
 
     """
-    location = (0, 0)
+    bounds = get_tissue_com_tile(wsi, 1024)
+    location = bounds[:2]
     size = np.array([1024, 1024])
     # Avoid testing very small levels (e.g. as in Omnyx JP2) because
     # MSE for very small levels is noisy.
@@ -2725,7 +2761,7 @@ def test_read_bounds_level_consistency(wsi: WSIReader) -> None:
     they are aligned.
 
     """
-    bounds = (0, 0, 1024, 1024)
+    bounds = get_tissue_com_tile(wsi, 1024)
     # This logic can be moved from the helper to here when other
     # reader classes have been parameterised into scenarios also.
     read_bounds_level_consistency(wsi, bounds)
@@ -2770,15 +2806,17 @@ def test_read_rect_coord_space_consistency(wsi: WSIReader) -> None:
     will not be of the same size, but the field of view will match.
 
     """
+    bounds = get_tissue_com_tile(wsi, 2000)
+    location = (bounds[:2] // 2) * 2  # ensure even coordinates
     roi1 = wsi.read_rect(
-        np.array([500, 500]),
+        location,
         np.array([2000, 2000]),
         coord_space="baseline",
         resolution=1.00,
         units="baseline",
     )
     roi2 = wsi.read_rect(
-        np.array([250, 250]),
+        location // 2,
         np.array([1000, 1000]),
         coord_space="resolution",
         resolution=0.5,
@@ -2979,3 +3017,22 @@ def test_fsspec_reader_open_pass_empty_json(tmp_path: Path) -> None:
     json_path.write_text("{}")
 
     assert not FsspecJsonWSIReader.is_valid_zarr_fsspec(str(json_path))
+
+
+def test_oob_read_dicom(sample_dicom: Path) -> None:
+    """Test that out of bounds returns background value.
+
+    For consistency with openslide, our readers should return a
+    background tile when reading out of bounds.
+
+    """
+    wsi = DICOMWSIReader(sample_dicom)
+    # Read a region that is out of bounds
+    region = wsi.read_rect(
+        location=(200000, 200),
+        size=(100, 100),
+    )
+    # Check that the region is the same size as the requested size
+    assert region.shape == (100, 100, 3)
+    # Check that the region is white (255)
+    assert np.all(region == 255)

tiatoolbox/visualization/bokeh_app/main.py

@@ -725,7 +725,16 @@ def populate_slide_list(slide_folder: Path, search_txt: str | None = None) -> No
     """Populate the slide list with the available slides."""
     file_list = []
     len_slidepath = len(slide_folder.parts)
-    for ext in ["*.svs", "*ndpi", "*.tiff", "*.mrxs", "*.jpg", "*.png", "*.tif"]:
+    for ext in [
+        "*.svs",
+        "*ndpi",
+        "*.tiff",
+        "*.mrxs",
+        "*.jpg",
+        "*.png",
+        "*.tif",
+        "*.dcm",
+    ]:
         file_list.extend(list(Path(slide_folder).glob(str(Path("*") / ext))))
         file_list.extend(list(Path(slide_folder).glob(ext)))
     if search_txt is None:
@@ -2086,7 +2095,16 @@ def setup_doc(self: DocConfig, base_doc: Document) -> tuple[Row, Tabs]:
 
         # Set initial slide to first one in base folder
         slide_list = []
-        for ext in ["*.svs", "*ndpi", "*.tiff", "*.tif", "*.mrxs", "*.png", "*.jpg"]:
+        for ext in [
+            "*.svs",
+            "*ndpi",
+            "*.tiff",
+            "*.tif",
+            "*.mrxs",
+            "*.png",
+            "*.jpg",
+            "*.dcm",
+        ]:
             slide_list.extend(list(doc_config["slide_folder"].glob(ext)))
             slide_list.extend(
                 list(doc_config["slide_folder"].glob(str(Path("*") / ext))),

tiatoolbox/wsicore/wsireader.py

@@ -4595,6 +4595,7 @@ def _info(self: DICOMWSIReader) -> WSIMeta:
             mpp=mpp,
             level_count=len(level_dimensions),
             vendor=dataset.Manufacturer,
+            file_path=self.input_path,
         )
 
     def read_rect(
@@ -4826,8 +4827,19 @@ def read_rect(
         _, constrained_read_size = utils.transforms.bounds2locsize(
             constrained_read_bounds,
         )
+
+        # if out of bounds, return empty image consistent with openslide
+        if np.any(np.array(constrained_read_size) <= 0):
+            return (
+                np.ones(
+                    shape=(int(size[1]), int(size[0]), 3),
+                    dtype=np.uint8,
+                )
+                * 255
+            )
+
         dicom_level = wsi.levels[read_level].level
-        im_region = wsi.read_region(location, dicom_level, constrained_read_size)
+        im_region = wsi.read_region(level_location, dicom_level, constrained_read_size)
         im_region = np.array(im_region)
 
         # Apply padding outside the slide area
@@ -5003,7 +5015,6 @@ class docstrings for more information.
         wsi = self.wsi
 
         # Read at optimal level and corrected read size
-        location_at_baseline = bounds_at_baseline[:2]
         level_location, size_at_read_level = utils.transforms.bounds2locsize(
             bounds_at_read_level,
         )
@@ -5016,7 +5027,7 @@ class docstrings for more information.
         _, read_size = utils.transforms.bounds2locsize(read_bounds)
         dicom_level = wsi.levels[read_level].level
         im_region = wsi.read_region(
-            location=location_at_baseline,
+            location=level_location,
             level=dicom_level,
             size=read_size,
         )