5809 hedjitter

monai/apps/pathology/transforms/stain/array.py

@@ -11,9 +11,75 @@
 
 from __future__ import annotations
 
+import numbers
+
 import numpy as np
 
 from monai.transforms.transform import Transform
+from monai.utils import optional_import
+
+color, _ = optional_import("skimage", name="color")
+
+
+class HEDJitter(Transform):
+    """Randomly perturbe the HED color space value an RGB image.
+    First, it disentangled the hematoxylin and eosin color channels by color deconvolution method using a fixed matrix.
+    Second, it perturbed the hematoxylin, eosin and DAB stains independently.
+    Third, it transformed the resulting stains into regular RGB color space.
+    Args:
+        theta (float): How much to jitter HED color space,
+         alpha is chosen from a uniform distribution [1-theta, 1+theta]
+         betti is chosen from a uniform distribution [-theta, theta]
+         the jitter formula is **s' = \alpha * s + \betti**
+
+    Note:
+        For more information refer to:
+            - the paper benchmarking different stain augmentation and normalization methods including HED:
+                Tellez et al. 2020 https://arxiv.org/abs/1902.06543
+            - the original paper on color deconvolution: Ruifrok et al. 2001
+                https://helios2.mi.parisdescartes.fr/~lomn/Cours/CV/BME/HistoPatho/Color/PythonColorDeconv/Quantification_of_histochemical_staining.pdf
+            - previous Pytorch implementation: https://github.com/gatsby2016/Augmentation-PyTorch-Transforms
+
+    """
+
+    def __init__(self, theta: float = 0.0) -> None:
+        assert isinstance(theta, numbers.Number), "theta should be a single number."
+        self.theta = theta
+        self.alpha = np.random.uniform(1 - theta, 1 + theta, (1, 3))
+        self.betti = np.random.uniform(-theta, theta, (1, 3))
+
+    def adjust_hed(self, img, alpha, betti):
+        # img = np.array(img)
+
+        s = np.reshape(color.rgb2hed(img), (-1, 3))
+        ns = alpha * s + betti  # perturbations on HED color space
+        nimg = color.hed2rgb(np.reshape(ns, img.shape))
+
+        imin = nimg.min()
+        imax = nimg.max()
+        rsimg = (255 * (nimg - imin) / (imax - imin)).astype("uint8")  # rescale to [0,255]
+        # transfer to PIL image
+        return rsimg
+
+    def __call__(self, image: np.ndarray) -> np.ndarray:
+        """
+        Args:
+            image: uint8 RGB image to perform HEDJitter on.
+
+        Returns:
+            perturbed_image: uint8 RGB image HED preturbed image.
+        """
+        if not isinstance(image, np.ndarray):
+            raise TypeError("Image must be of type numpy.ndarray.")
+        perturbed_image = self.adjust_hed(image, self.alpha, self.betti)
+        return perturbed_image
+
+    def __repr__(self):
+        format_string = self.__class__.__name__ + "("
+        format_string += f"theta={self.theta}"
+        format_string += f",alpha={self.alpha}"
+        format_string += f",betti={self.betti}"
+        return format_string
 
 
 class ExtractHEStains(Transform):

monai/apps/pathology/transforms/stain/dictionary.py

@@ -17,14 +17,41 @@
 
 from __future__ import annotations
 
+import numbers
 from collections.abc import Hashable, Mapping
 
 import numpy as np
 
 from monai.config import KeysCollection
 from monai.transforms.transform import MapTransform
 
-from .array import ExtractHEStains, NormalizeHEStains
+from .array import ExtractHEStains, HEDJitter, NormalizeHEStains
+
+
+class HEDJitterd(MapTransform):
+
+    """Dictionary-based wrapper of :py:class:`monai.apps.pathology.transforms.HEDJitter`.
+    Class to randomly perturbe HED color space values of image using stain deconvolution.
+
+    Args:
+        keys: keys of the corresponding items to be transformed.
+            See also: :py:class:`monai.transforms.compose.MapTransform`
+        theta: jitter range factor
+
+        allow_missing_keys: don't raise exception if key is missing.
+
+    """
+
+    def __init__(self, keys: KeysCollection, theta: float = 0.0, allow_missing_keys: bool = False) -> None:
+        super().__init__(keys, allow_missing_keys)
+        assert isinstance(theta, numbers.Number), "theta should be a single number."
+        self.hedjitter = HEDJitter(theta=theta)
+
+    def __call__(self, data: Mapping[Hashable, np.ndarray]) -> dict[Hashable, np.ndarray]:
+        d = dict(data)
+        for key in self.key_iterator(d):
+            d[key] = self.hedjitter(d[key])
+        return d
 
 
 class ExtractHEStainsd(MapTransform):
@@ -111,3 +138,4 @@ def __call__(self, data: Mapping[Hashable, np.ndarray]) -> dict[Hashable, np.nda
 
 ExtractHEStainsDict = ExtractHEStainsD = ExtractHEStainsd
 NormalizeHEStainsDict = NormalizeHEStainsD = NormalizeHEStainsd
+HEDJitterDict = HEDJitterD = HEDJitterd

tests/test_pathology_hedjitter.py

@@ -0,0 +1,223 @@
+# Copyright (c) MONAI Consortium
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#     http://www.apache.org/licenses/LICENSE-2.0
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import annotations
+
+import unittest
+
+import numpy as np
+from parameterized import parameterized
+
+from monai.apps.pathology.transforms import ExtractHEStainsD, NormalizeHEStainsD
+
+# None inputs
+EXTRACT_STAINS_TEST_CASE_0 = (None,)
+EXTRACT_STAINS_TEST_CASE_00 = (None, None)
+NORMALIZE_STAINS_TEST_CASE_0 = (None,)
+NORMALIZE_STAINS_TEST_CASE_00: tuple = ({}, None, None)
+
+# input pixels all transparent and below the beta absorbance threshold
+EXTRACT_STAINS_TEST_CASE_1 = [np.full((3, 2, 3), 240)]
+
+# input pixels uniformly filled, but above beta absorbance threshold
+EXTRACT_STAINS_TEST_CASE_2 = [np.full((3, 2, 3), 100)]
+
+# input pixels uniformly filled (different value), but above beta absorbance threshold
+EXTRACT_STAINS_TEST_CASE_3 = [np.full((3, 2, 3), 150)]
+
+# input pixels uniformly filled with zeros, leading to two identical stains extracted
+EXTRACT_STAINS_TEST_CASE_4 = [
+    np.zeros((3, 2, 3)),
+    np.array([[0.0, 0.0], [0.70710678, 0.70710678], [0.70710678, 0.70710678]]),
+]
+
+# input pixels not uniformly filled, leading to two different stains extracted
+EXTRACT_STAINS_TEST_CASE_5 = [
+    np.array([[[100, 0, 0], [0, 0, 0]], [[0, 0, 0], [0, 0, 0]], [[0, 0, 0], [0, 0, 0]]]),
+    np.array([[0.70710677, 0.18696113], [0.0, 0.0], [0.70710677, 0.98236734]]),
+]
+
+# input pixels all transparent and below the beta absorbance threshold
+NORMALIZE_STAINS_TEST_CASE_1 = [np.full((3, 2, 3), 240)]
+
+# input pixels uniformly filled with zeros, and target stain matrix provided
+NORMALIZE_STAINS_TEST_CASE_2 = [{"target_he": np.full((3, 2), 1)}, np.zeros((3, 2, 3)), np.full((3, 2, 3), 11)]
+
+# input pixels uniformly filled with zeros, and target stain matrix not provided
+NORMALIZE_STAINS_TEST_CASE_3 = [
+    {},
+    np.zeros((3, 2, 3)),
+    np.array([[[63, 25, 60], [63, 25, 60]], [[63, 25, 60], [63, 25, 60]], [[63, 25, 60], [63, 25, 60]]]),
+]
+
+# input pixels not uniformly filled
+NORMALIZE_STAINS_TEST_CASE_4 = [
+    {"target_he": np.full((3, 2), 1)},
+    np.array([[[100, 0, 0], [0, 0, 0]], [[0, 0, 0], [0, 0, 0]], [[0, 0, 0], [0, 0, 0]]]),
+    np.array([[[87, 87, 87], [33, 33, 33]], [[33, 33, 33], [33, 33, 33]], [[33, 33, 33], [33, 33, 33]]]),
+]
+
+
+class TestExtractHEStainsD(unittest.TestCase):
+    @parameterized.expand([EXTRACT_STAINS_TEST_CASE_0, EXTRACT_STAINS_TEST_CASE_1])
+    def test_transparent_image(self, image):
+        """
+        Test HE stain extraction on an image that comprises
+        only transparent pixels - pixels with absorbance below the
+        beta absorbance threshold. A ValueError should be raised,
+        since once the transparent pixels are removed, there are no
+        remaining pixels to compute eigenvectors.
+        """
+        key = "image"
+        if image is None:
+            with self.assertRaises(TypeError):
+                ExtractHEStainsD([key])({key: image})
+        else:
+            with self.assertRaises(ValueError):
+                ExtractHEStainsD([key])({key: image})
+
+    @parameterized.expand([EXTRACT_STAINS_TEST_CASE_0, EXTRACT_STAINS_TEST_CASE_2, EXTRACT_STAINS_TEST_CASE_3])
+    def test_identical_result_vectors(self, image):
+        """
+        Test HE stain extraction on input images that are
+        uniformly filled with pixels that have absorbance above the
+        beta absorbance threshold. Since input image is uniformly filled,
+        the two extracted stains should have the same RGB values. So,
+        we assert that the first column is equal to the second column
+        of the returned stain matrix.
+        """
+        key = "image"
+        if image is None:
+            with self.assertRaises(TypeError):
+                ExtractHEStainsD([key])({key: image})
+        else:
+            result = ExtractHEStainsD([key])({key: image})
+            np.testing.assert_array_equal(result[key][:, 0], result[key][:, 1])
+
+    @parameterized.expand([EXTRACT_STAINS_TEST_CASE_00, EXTRACT_STAINS_TEST_CASE_4, EXTRACT_STAINS_TEST_CASE_5])
+    def test_result_value(self, image, expected_data):
+        """
+        Test that an input image returns an expected stain matrix.
+
+        For test case 4:
+        - a uniformly filled input image should result in
+          eigenvectors [[1,0,0],[0,1,0],[0,0,1]]
+        - phi should be an array containing only values of
+          arctan(1) since the ratio between the eigenvectors
+          corresponding to the two largest eigenvalues is 1
+        - maximum phi and minimum phi should thus be arctan(1)
+        - thus, maximum vector and minimum vector should be
+          [[0],[0.70710677],[0.70710677]]
+        - the resulting extracted stain should be
+          [[0,0],[0.70710678,0.70710678],[0.70710678,0.70710678]]
+
+        For test case 5:
+        - the non-uniformly filled input image should result in
+          eigenvectors [[0,0,1],[1,0,0],[0,1,0]]
+        - maximum phi and minimum phi should thus be 0.785 and
+          0.188 respectively
+        - thus, maximum vector and minimum vector should be
+          [[0.18696113],[0],[0.98236734]] and
+          [[0.70710677],[0],[0.70710677]] respectively
+        - the resulting extracted stain should be
+          [[0.70710677,0.18696113],[0,0],[0.70710677,0.98236734]]
+        """
+        key = "image"
+        if image is None:
+            with self.assertRaises(TypeError):
+                ExtractHEStainsD([key])({key: image})
+        else:
+            result = ExtractHEStainsD([key])({key: image})
+            np.testing.assert_allclose(result[key], expected_data)
+
+
+class TestNormalizeHEStainsD(unittest.TestCase):
+    @parameterized.expand([NORMALIZE_STAINS_TEST_CASE_0, NORMALIZE_STAINS_TEST_CASE_1])
+    def test_transparent_image(self, image):
+        """
+        Test HE stain normalization on an image that comprises
+        only transparent pixels - pixels with absorbance below the
+        beta absorbance threshold. A ValueError should be raised,
+        since once the transparent pixels are removed, there are no
+        remaining pixels to compute eigenvectors.
+        """
+        key = "image"
+        if image is None:
+            with self.assertRaises(TypeError):
+                NormalizeHEStainsD([key])({key: image})
+        else:
+            with self.assertRaises(ValueError):
+                NormalizeHEStainsD([key])({key: image})
+
+    @parameterized.expand(
+        [
+            NORMALIZE_STAINS_TEST_CASE_00,
+            NORMALIZE_STAINS_TEST_CASE_2,
+            NORMALIZE_STAINS_TEST_CASE_3,
+            NORMALIZE_STAINS_TEST_CASE_4,
+        ]
+    )
+    def test_result_value(self, arguments, image, expected_data):
+        """
+        Test that an input image returns an expected normalized image.
+
+        For test case 2:
+        - This case tests calling the stain normalizer, after the
+          _deconvolution_extract_conc function. This is because the normalized
+          concentration returned for each pixel is the same as the reference
+          maximum stain concentrations in the case that the image is uniformly
+          filled, as in this test case. This is because the maximum concentration
+          for each stain is the same as each pixel's concentration.
+        - Thus, the normalized concentration matrix should be a (2, 6) matrix
+          with the first row having all values of 1.9705, second row all 1.0308.
+        - Taking the matrix product of the target stain matrix and the concentration
+          matrix, then using the inverse Beer-Lambert transform to obtain the RGB
+          image from the absorbance image, and finally converting to uint8,
+          we get that the stain normalized image should be a matrix of
+          dims (3, 2, 3), with all values 11.
+
+        For test case 3:
+        - This case also tests calling the stain normalizer, after the
+          _deconvolution_extract_conc function returns the image concentration
+          matrix.
+        - As in test case 2, the normalized concentration matrix should be a (2, 6) matrix
+          with the first row having all values of 1.9705, second row all 1.0308.
+        - Taking the matrix product of the target default stain matrix and the concentration
+          matrix, then using the inverse Beer-Lambert transform to obtain the RGB
+          image from the absorbance image, and finally converting to uint8,
+          we get that the stain normalized image should be [[[63, 25, 60], [63, 25, 60]],
+          [[63, 25, 60], [63, 25, 60]], [[63, 25, 60], [63, 25, 60]]]
+
+        For test case 4:
+        - For this non-uniformly filled image, the stain extracted should be
+          [[0.70710677,0.18696113],[0,0],[0.70710677,0.98236734]], as validated for the
+          ExtractHEStains class. Solving the linear least squares problem (since
+          absorbance matrix = stain matrix * concentration matrix), we obtain the concentration
+          matrix that should be [[-0.3101, 7.7508, 7.7508, 7.7508, 7.7508, 7.7508],
+          [5.8022, 0, 0, 0, 0, 0]]
+        - Normalizing the concentration matrix, taking the matrix product of the
+          target stain matrix and the concentration matrix, using the inverse
+          Beer-Lambert transform to obtain the RGB image from the absorbance
+          image, and finally converting to uint8, we get that the stain normalized
+          image should be [[[87, 87, 87], [33, 33, 33]], [[33, 33, 33], [33, 33, 33]],
+          [[33, 33, 33], [33, 33, 33]]]
+        """
+        key = "image"
+        if image is None:
+            with self.assertRaises(TypeError):
+                NormalizeHEStainsD([key])({key: image})
+        else:
+            result = NormalizeHEStainsD([key], **arguments)({key: image})
+            np.testing.assert_allclose(result[key], expected_data)
+
+
+if __name__ == "__main__":
+    unittest.main()