Merge pull request #301 from Dana-Farber-AIOS/update-quantifyMIF

Update quantify mif

pathml/preprocessing/transforms.py

@@ -12,11 +12,15 @@
 import pandas as pd
 import pathml.core
 import pathml.core.slide_data
-from pathml.utils import (RGB_to_GREY, RGB_to_HSI, RGB_to_HSV, RGB_to_OD,
-                          normalize_matrix_cols)
+from pathml.utils import (
+    RGB_to_GREY,
+    RGB_to_HSI,
+    RGB_to_HSV,
+    RGB_to_OD,
+    normalize_matrix_cols,
+)
 from skimage import restoration
-from skimage.exposure import (equalize_adapthist, equalize_hist,
-                              rescale_intensity)
+from skimage.exposure import equalize_adapthist, equalize_hist, rescale_intensity
 from skimage.measure import regionprops_table
 
 
@@ -1418,7 +1422,7 @@ class QuantifyMIF(Transform):
     """
     Convert segmented image into anndata.AnnData counts object `AnnData <https://anndata.readthedocs.io/en/latest/>`_.
     Counts objects are used to interface with the Python single cell analysis ecosystem `Scanpy <https://scanpy.readthedocs.io/en/stable/>`_.
-    The counts object contains a summary of protein expression statistics in each cell along with its coordinate.
+    The counts object contains a summary of channel statistics in each cell along with its coordinate.
 
     Args:
         segmentation_mask (str): key indicating which mask to use as label image
@@ -1430,14 +1434,30 @@ def __init__(self, segmentation_mask):
     def __repr__(self):
         return f"QuantifyMIF(segmentation_mask={self.segmentation_mask})"
 
-    def F(self, tile):
-        # pass (x, y, channel) image and (x, y) segmentation
-        img = tile.image.copy()
-        segmentation = tile.masks[self.segmentation_mask][:, :, 0]
+    def F(self, img, segmentation, coords_offset=(0, 0)):
+        """
+        Functional implementation
+
+        Args:
+            img (np.ndarray): Input image of shape (i, j, n_channels)
+            segmentation (np.ndarray): Segmentation map of shape (i, j) or (i, j, 1). Zeros are background. Regions should be
+                labelled with unique integers.
+            coords_offset (tuple, optional): Coordinates (i, j) used to convert tile-level coordinates to slide-level.
+                Defaults to (0, 0) for no offset.
+
+        Returns:
+            Counts matrix
+        """
+        if segmentation.ndim != 2:
+            assert (
+                segmentation.shape[2] == 1
+            ), f"input segmentation is of shape {segmentation.shape}. must be (x, y) or (x, y, 1)"
+            segmentation = segmentation.squeeze(2)
         countsdataframe = regionprops_table(
             label_image=segmentation,
             intensity_image=img,
             properties=[
+                "label",
                 "coords",
                 "max_intensity",
                 "mean_intensity",
@@ -1453,15 +1473,17 @@ def F(self, tile):
         for i in range(img.shape[-1]):
             X[i] = countsdataframe[f"mean_intensity-{i}"]
         # populate anndata object
+        # i,j are relative to the input image (0 to img.shape). Adding offset converts to slide-level coordinates
         counts = anndata.AnnData(
             X=X,
             obs=[
-                tuple([i + tile.coords[0], j + tile.coords[1]])
+                tuple([i + coords_offset[0], j + coords_offset[1]])
                 for i, j in zip(
                     countsdataframe["centroid-0"], countsdataframe["centroid-1"]
                 )
             ],
         )
+        counts.obs["label"] = countsdataframe["label"]
         counts.obs = counts.obs.rename(columns={0: "y", 1: "x"})
         counts.obs["filled_area"] = countsdataframe["filled_area"]
         counts.obs["euler_number"] = countsdataframe["euler_number"]
@@ -1489,7 +1511,11 @@ def apply(self, tile):
         assert (
             tile.slide_type.stain == "Fluor"
         ), f"Tile has slide_type.stain='{tile.slide_type.stain}', but must be 'Fluor'"
-        tile.counts = self.F(tile)
+        tile.counts = self.F(
+            img=tile.image,
+            segmentation=tile.masks[self.segmentation_mask],
+            coords_offset=tile.coords,
+        )
 
 
 class CollapseRunsVectra(Transform):

tests/preprocessing_tests/test_pipeline.py

@@ -100,7 +100,9 @@ def test_pipeline_mif(tileVectra):
     ).F(collapsed_im)
     orig_tile.image = collapsed_im
     orig_tile.masks["segmentation_mask"] = cell_segmentation
-    adata = QuantifyMIF(segmentation_mask="segmentation_mask").F(orig_tile)
+    adata = QuantifyMIF(segmentation_mask="segmentation_mask").F(
+        orig_tile.image, orig_tile.masks["segmentation_mask"], orig_tile.coords
+    )
 
     assert np.array_equal(tileVectra.masks["cell_segmentation"], cell_segmentation)
     pd.testing.assert_frame_equal(adata.obs, tileVectra.counts.obs)