Stain unmixing fix

doc/examples/color_exposure/plot_ihc_color_separation.py

@@ -19,35 +19,38 @@
        American Society of Cytology, vol. 23, no. 4, pp. 291-9, Aug. 2001.
 
 """
+import numpy as np
 import matplotlib.pyplot as plt
 
 from skimage import data
-from skimage.color import rgb2hed
-from matplotlib.colors import LinearSegmentedColormap
-
-# Create an artificial color close to the original one
-cmap_hema = LinearSegmentedColormap.from_list('mycmap', ['white', 'navy'])
-cmap_dab = LinearSegmentedColormap.from_list('mycmap', ['white',
-                                             'saddlebrown'])
-cmap_eosin = LinearSegmentedColormap.from_list('mycmap', ['darkviolet',
-                                               'white'])
+from skimage.color import rgb2hed, hed2rgb
 
+# Example IHC image
 ihc_rgb = data.immunohistochemistry()
+
+# Separate the stains from the IHC image
 ihc_hed = rgb2hed(ihc_rgb)
 
+# Create an RGB image for each of the stains
+null = np.zeros_like(ihc_hed[:, :, 0])
+ihc_h = hed2rgb(np.stack((ihc_hed[:, :, 0], null, null), axis=-1))
+ihc_e = hed2rgb(np.stack((null, ihc_hed[:, :, 1], null), axis=-1))
+ihc_d = hed2rgb(np.stack((null, null, ihc_hed[:, :, 2]), axis=-1))
+
+# Display
 fig, axes = plt.subplots(2, 2, figsize=(7, 6), sharex=True, sharey=True)
 ax = axes.ravel()
 
 ax[0].imshow(ihc_rgb)
 ax[0].set_title("Original image")
 
-ax[1].imshow(ihc_hed[:, :, 0], cmap=cmap_hema)
+ax[1].imshow(ihc_h)
 ax[1].set_title("Hematoxylin")
 
-ax[2].imshow(ihc_hed[:, :, 1], cmap=cmap_eosin)
-ax[2].set_title("Eosin")
+ax[2].imshow(ihc_e)
+ax[2].set_title("Eosin")  # Note that there is no Eosin stain in this image
 
-ax[3].imshow(ihc_hed[:, :, 2], cmap=cmap_dab)
+ax[3].imshow(ihc_d)
 ax[3].set_title("DAB")
 
 for a in ax.ravel():
@@ -59,13 +62,16 @@
 ######################################################################
 # Now we can easily manipulate the hematoxylin and DAB "channels":
 
-import numpy as np
 from skimage.exposure import rescale_intensity
 
 # Rescale hematoxylin and DAB signals and give them a fluorescence look
-h = rescale_intensity(ihc_hed[:, :, 0], out_range=(0, 1))
-d = rescale_intensity(ihc_hed[:, :, 2], out_range=(0, 1))
-zdh = np.dstack((np.zeros_like(h), d, h))
+h = rescale_intensity(ihc_hed[:, :, 0], out_range=(0, 1),
+                      in_range=(0, np.percentile(ihc_hed[:, :, 0], 99)))
+d = rescale_intensity(ihc_hed[:, :, 2], out_range=(0, 1),
+                      in_range=(0, np.percentile(ihc_hed[:, :, 2], 99)))
+
+# Put the two channels into an RGB image as green and blue channels
+zdh = np.dstack((null, d, h))
 
 fig = plt.figure()
 axis = plt.subplot(1, 1, 1, sharex=ax[0], sharey=ax[0])

skimage/color/colorconv.py

@@ -1451,6 +1451,8 @@ def separate_stains(rgb, conv_matrix):
 
     stains = (np.log(rgb) / log_adjust) @ conv_matrix
 
+    np.maximum(stains, 0, out=stains)
+
     return stains
 
 

skimage/color/tests/test_colorconv.py

@@ -47,6 +47,7 @@ class TestColorconv(TestCase):
     img_rgba = np.array([[[0, 0.5, 1, 0],
                           [0, 0.5, 1, 1],
                           [0, 0.5, 1, 0.5]]]).astype(float)
+    img_stains = img_as_float(img_rgb) * 0.3
 
     colbars = np.array([[1, 1, 0, 0, 1, 1, 0, 0],
                         [1, 1, 1, 1, 0, 0, 0, 0],
@@ -183,32 +184,33 @@ def test_xyz_rgb_roundtrip(self):
         img_rgb = img_as_float(self.img_rgb)
         assert_array_almost_equal(xyz2rgb(rgb2xyz(img_rgb)), img_rgb)
 
-    # RGB<->HED roundtrip with ubyte image
+    # HED<->RGB roundtrip with ubyte image
     def test_hed_rgb_roundtrip(self):
-        img_rgb = img_as_ubyte(self.img_rgb)
-        new = img_as_ubyte(hed2rgb(rgb2hed(img_rgb)))
-        assert_equal(new, img_rgb)
+        img_in = img_as_ubyte(self.img_stains)
+        img_out = rgb2hed(hed2rgb(img_in))
+        assert_equal(img_as_ubyte(img_out), img_in)
 
-    # RGB<->HED roundtrip with float image
+    # HED<->RGB roundtrip with float image
     def test_hed_rgb_float_roundtrip(self):
-        img_rgb = img_as_float(self.img_rgb)
-        assert_array_almost_equal(hed2rgb(rgb2hed(img_rgb)), img_rgb)
-
-    # RGB<->HDX roundtrip with ubyte image
-    def test_hdx_rgb_roundtrip(self):
-        from skimage.color.colorconv import hdx_from_rgb, rgb_from_hdx
-        img_rgb = self.img_rgb
-        conv = combine_stains(separate_stains(img_rgb, hdx_from_rgb),
-                              rgb_from_hdx)
-        assert_equal(img_as_ubyte(conv), img_rgb)
-
-    # RGB<->HDX roundtrip with float image
-    def test_hdx_rgb_roundtrip_float(self):
-        from skimage.color.colorconv import hdx_from_rgb, rgb_from_hdx
-        img_rgb = img_as_float(self.img_rgb)
-        conv = combine_stains(separate_stains(img_rgb, hdx_from_rgb),
-                              rgb_from_hdx)
-        assert_array_almost_equal(conv, img_rgb)
+        img_in = self.img_stains
+        img_out = rgb2hed(hed2rgb(img_in))
+        assert_array_almost_equal(img_out, img_in)
+
+    # BRO<->RGB roundtrip with ubyte image
+    def test_bro_rgb_roundtrip(self):
+        from skimage.color.colorconv import bro_from_rgb, rgb_from_bro
+        img_in = img_as_ubyte(self.img_stains)
+        img_out = combine_stains(img_in, rgb_from_bro)
+        img_out = separate_stains(img_out, bro_from_rgb)
+        assert_equal(img_as_ubyte(img_out), img_in)
+
+    # BRO<->RGB roundtrip with float image
+    def test_bro_rgb_roundtrip_float(self):
+        from skimage.color.colorconv import bro_from_rgb, rgb_from_bro
+        img_in = self.img_stains
+        img_out = combine_stains(img_in, rgb_from_bro)
+        img_out = separate_stains(img_out, bro_from_rgb)
+        assert_array_almost_equal(img_out, img_in)
 
     # RGB to RGB CIE
     def test_rgb2rgbcie_conversion(self):