Merge pull request #441 from spotter/color_deconv

ENH: Adding color deconvolution for immunohistochemical images.

CONTRIBUTORS.txt

@@ -136,3 +136,6 @@
 
 - Thouis Jones
   Vectorized operators for arrays of 16-bit ints.
+
+- Xavier Moles Lopez
+  Color separation (color deconvolution) for several stainings.

doc/examples/plot_ihc_color_separation.py

@@ -0,0 +1,70 @@
+"""
+==============================================
+Immunohistochemical staining colors separation
+==============================================
+
+In this example we separate the immunohistochemical (IHC) staining
+from the hematoxylin counterstaining. The separation is achieved with the
+method described in [1]_, known as "color deconvolution".
+
+The IHC staining expression of the FHL2 protein is here revealed with
+Diaminobenzidine (DAB) which gives a brown color.
+
+
+.. [1] A. C. Ruifrok and D. A. Johnston, "Quantification of histochemical
+       staining by color deconvolution.," Analytical and quantitative
+       cytology and histology / the International Academy of Cytology [and]
+       American Society of Cytology, vol. 23, no. 4, pp. 291-9, Aug. 2001.
+"""
+import matplotlib.pyplot as plt
+
+from skimage import data
+from skimage.color import rgb2hed
+
+ihc_rgb = data.immunohistochemistry()
+ihc_hed = rgb2hed(ihc_rgb)
+
+fig, axes = plt.subplots(2, 2, figsize=(7, 6))
+ax0, ax1, ax2, ax3 = axes.ravel()
+
+ax0.imshow(ihc_rgb)
+ax0.set_title("Original image")
+
+ax1.imshow(ihc_hed[:, :, 0], cmap=plt.cm.gray)
+ax1.set_title("Hematoxylin")
+
+ax2.imshow(ihc_hed[:, :, 1], cmap=plt.cm.gray)
+ax2.set_title("Eosin")
+
+ax3.imshow(ihc_hed[:, :, 2], cmap=plt.cm.gray)
+ax3.set_title("DAB")
+
+for ax in axes.ravel():
+    ax.axis('off')
+
+fig.subplots_adjust(hspace=0.3)
+
+
+"""
+.. image:: PLOT2RST.current_figure
+
+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))
+
+plt.figure()
+plt.imshow(zdh)
+plt.title("Stain separated image (rescaled)")
+plt.axis('off')
+plt.show()
+
+"""
+.. image:: PLOT2RST.current_figure
+"""

skimage/color/colorconv.py

@@ -45,7 +45,14 @@ from __future__ import division
 
 __all__ = ['convert_colorspace', 'rgb2hsv', 'hsv2rgb', 'rgb2xyz', 'xyz2rgb',
            'rgb2rgbcie', 'rgbcie2rgb', 'rgb2grey', 'rgb2gray', 'gray2rgb',
-           'xyz2lab', 'lab2xyz', 'lab2rgb', 'rgb2lab', 'is_rgb', 'is_gray'
+           'xyz2lab', 'lab2xyz', 'lab2rgb', 'rgb2lab', 'rgb2hed', 'hed2rgb',
+           'separate_stains', 'combine_stains', 'rgb_from_hed', 'hed_from_rgb',
+           'rgb_from_hdx', 'hdx_from_rgb', 'rgb_from_fgx', 'fgx_from_rgb',
+           'rgb_from_bex', 'bex_from_rgb', 'rgb_from_rbd', 'rbd_from_rgb',
+           'rgb_from_gdx', 'gdx_from_rgb', 'rgb_from_hax', 'hax_from_rgb',
+           'rgb_from_bro', 'bro_from_rgb', 'rgb_from_bpx', 'bpx_from_rgb',
+           'rgb_from_ahx', 'ahx_from_rgb', 'rgb_from_hpx', 'hpx_from_rgb',
+           'is_rgb', 'is_gray'
            ]
 
 __docformat__ = "restructuredtext en"
@@ -312,6 +319,90 @@ gray_from_rgb = np.array([[0.2125, 0.7154, 0.0721],
 # CIE LAB constants for Observer= 2A, Illuminant= D65
 lab_ref_white = np.array([0.95047, 1., 1.08883])
 
+
+# Haematoxylin-Eosin-DAB colorspace
+# From original Ruifrok's paper: A. C. Ruifrok and D. A. Johnston,
+# "Quantification of histochemical staining by color deconvolution.,"
+# Analytical and quantitative cytology and histology / the International
+# Academy of Cytology [and] American Society of Cytology, vol. 23, no. 4,
+# pp. 291-9, Aug. 2001.
+rgb_from_hed = np.array([[0.65, 0.70, 0.29],
+                         [0.07, 0.99, 0.11],
+                         [0.27, 0.57, 0.78]])
+hed_from_rgb = linalg.inv(rgb_from_hed)
+
+# Following matrices are adapted form the Java code written by G.Landini.
+# The original code is available at:
+# http://www.dentistry.bham.ac.uk/landinig/software/cdeconv/cdeconv.html
+
+# Hematoxylin + DAB
+rgb_from_hdx = np.array([[0.650, 0.704, 0.286],
+                         [0.268, 0.570, 0.776],
+                         [0.0, 0.0, 0.0]])
+rgb_from_hdx[2, :] = np.cross(rgb_from_hdx[0, :], rgb_from_hdx[1, :])
+hdx_from_rgb = linalg.inv(rgb_from_hdx)
+
+# Feulgen + Light Green
+rgb_from_fgx = np.array([[0.46420921, 0.83008335, 0.30827187],
+                         [0.94705542, 0.25373821, 0.19650764],
+                         [0.0, 0.0, 0.0]])
+rgb_from_fgx[2, :] = np.cross(rgb_from_fgx[0, :], rgb_from_fgx[1, :])
+fgx_from_rgb = linalg.inv(rgb_from_fgx)
+
+# Giemsa: Methyl Blue + Eosin
+rgb_from_bex = np.array([[0.834750233, 0.513556283, 0.196330403],
+                         [0.092789, 0.954111, 0.283111],
+                         [0.0, 0.0, 0.0]])
+rgb_from_bex[2, :] = np.cross(rgb_from_bex[0, :], rgb_from_bex[1, :])
+bex_from_rgb = linalg.inv(rgb_from_bex)
+
+# FastRed + FastBlue +  DAB
+rgb_from_rbd = np.array([[0.21393921, 0.85112669, 0.47794022],
+                         [0.74890292, 0.60624161, 0.26731082],
+                         [0.268, 0.570, 0.776]])
+rbd_from_rgb = linalg.inv(rgb_from_rbd)
+
+# Methyl Green + DAB
+rgb_from_gdx = np.array([[0.98003, 0.144316, 0.133146],
+                         [0.268, 0.570, 0.776],
+                         [0.0, 0.0, 0.0]])
+rgb_from_gdx[2, :] = np.cross(rgb_from_gdx[0, :], rgb_from_gdx[1, :])
+gdx_from_rgb = linalg.inv(rgb_from_gdx)
+
+# Hematoxylin + AEC
+rgb_from_hax = np.array([[0.650, 0.704, 0.286],
+                         [0.2743, 0.6796, 0.6803],
+                         [0.0, 0.0, 0.0]])
+rgb_from_hax[2, :] = np.cross(rgb_from_hax[0, :], rgb_from_hax[1, :])
+hax_from_rgb = linalg.inv(rgb_from_hax)
+
+# Blue matrix Anilline Blue + Red matrix Azocarmine + Orange matrix Orange-G
+rgb_from_bro = np.array([[0.853033, 0.508733, 0.112656],
+                         [0.09289875, 0.8662008, 0.49098468],
+                         [0.10732849, 0.36765403, 0.9237484]])
+bro_from_rgb = linalg.inv(rgb_from_bro)
+
+# Methyl Blue + Ponceau Fuchsin
+rgb_from_bpx = np.array([[0.7995107, 0.5913521, 0.10528667],
+                         [0.09997159, 0.73738605, 0.6680326],
+                         [0.0, 0.0, 0.0]])
+rgb_from_bpx[2, :] = np.cross(rgb_from_bpx[0, :], rgb_from_bpx[1, :])
+bpx_from_rgb = linalg.inv(rgb_from_bpx)
+
+# Alcian Blue + Hematoxylin
+rgb_from_ahx = np.array([[0.874622, 0.457711, 0.158256],
+                         [0.552556, 0.7544, 0.353744],
+                         [0.0, 0.0, 0.0]])
+rgb_from_ahx[2, :] = np.cross(rgb_from_ahx[0, :], rgb_from_ahx[1, :])
+ahx_from_rgb = linalg.inv(rgb_from_ahx)
+
+# Hematoxylin + PAS
+rgb_from_hpx = np.array([[0.644211, 0.716556, 0.266844],
+                         [0.175411, 0.972178, 0.154589],
+                         [0.0, 0.0, 0.0]])
+rgb_from_hpx[2, :] = np.cross(rgb_from_hpx[0, :], rgb_from_hpx[1, :])
+hpx_from_rgb = linalg.inv(rgb_from_hpx)
+
 #-------------------------------------------------------------
 # The conversion functions that make use of the matrices above
 #-------------------------------------------------------------
@@ -721,3 +812,189 @@ def lab2rgb(lab):
     This function uses lab2xyz and xyz2rgb.
     """
     return xyz2rgb(lab2xyz(lab))
+
+
+def rgb2hed(rgb):
+    """RGB to Haematoxylin-Eosin-DAB (HED) color space conversion.
+
+    Parameters
+    ----------
+    rgb : array_like
+        The image in RGB format, in a 3-D array of shape (.., .., 3).
+
+    Returns
+    -------
+    out : ndarray
+        The image in HED format, in a 3-D array of shape (.., .., 3).
+
+    Raises
+    ------
+    ValueError
+        If `rgb` is not a 3-D array of shape (.., .., 3).
+
+
+    References
+    ----------
+    .. [1] A. C. Ruifrok and D. A. Johnston, "Quantification of histochemical
+           staining by color deconvolution.," Analytical and quantitative
+           cytology and histology / the International Academy of Cytology [and]
+           American Society of Cytology, vol. 23, no. 4, pp. 291-9, Aug. 2001.
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.color import rgb2hed
+    >>> ihc = data.immunohistochemistry()
+    >>> ihc_hed = rgb2hed(ihc)
+    """
+    return separate_stains(rgb, hed_from_rgb)
+
+
+def hed2rgb(hed):
+    """Haematoxylin-Eosin-DAB (HED) to RGB color space conversion.
+
+    Parameters
+    ----------
+    hed : array_like
+        The image in the HED color space, in a 3-D array of shape (.., .., 3).
+
+    Returns
+    -------
+    out : ndarray
+        The image in RGB, in a 3-D array of shape (.., .., 3).
+
+    Raises
+    ------
+    ValueError
+        If `hed` is not a 3-D array of shape (.., .., 3).
+
+    References
+    ----------
+    .. [1] A. C. Ruifrok and D. A. Johnston, "Quantification of histochemical
+           staining by color deconvolution.," Analytical and quantitative
+           cytology and histology / the International Academy of Cytology [and]
+           American Society of Cytology, vol. 23, no. 4, pp. 291-9, Aug. 2001.
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.color import rgb2hed, hed2rgb
+    >>> ihc = data.immunohistochemistry()
+    >>> ihc_hed = rgb2hed(ihc)
+    >>> ihc_rgb = hed2rgb(ihc_hed)
+    """
+    return combine_stains(hed, rgb_from_hed)
+
+
+def separate_stains(rgb, conv_matrix):
+    """RGB to stain color space conversion.
+
+    Parameters
+    ----------
+    rgb : array_like
+        The image in RGB format, in a 3-D array of shape (.., .., 3).
+    conv_matrix: ndarray
+        The stain separation matrix as described by G. Landini [1]_.
+
+    Returns
+    -------
+    out : ndarray
+        The image in stain color space, in a 3-D array of shape (.., .., 3).
+
+    Raises
+    ------
+    ValueError
+        If `rgb` is not a 3-D array of shape (.., .., 3).
+
+    Notes
+    -----
+    Stain separation matrices available in the ``color`` module and their
+    respective colorspace:
+
+    * ``hed_from_rgb``: Hematoxylin + Eosin + DAB
+    * ``hdx_from_rgb``: Hematoxylin + DAB
+    * ``fgx_from_rgb``: Feulgen + Light Green
+    * ``bex_from_rgb``: Giemsa stain : Methyl Blue + Eosin
+    * ``rbd_from_rgb``: FastRed + FastBlue +  DAB
+    * ``gdx_from_rgb``: Methyl Green + DAB
+    * ``hax_from_rgb``: Hematoxylin + AEC
+    * ``bro_from_rgb``: Blue matrix Anilline Blue + Red matrix Azocarmine\
+                        + Orange matrix Orange-G
+    * ``bpx_from_rgb``: Methyl Blue + Ponceau Fuchsin
+    * ``ahx_from_rgb``: Alcian Blue + Hematoxylin
+    * ``hpx_from_rgb``: Hematoxylin + PAS
+
+    References
+    ----------
+    .. [1] http://www.dentistry.bham.ac.uk/landinig/software/cdeconv/cdeconv.html
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.color import separate_stains, hdx_from_rgb
+    >>> ihc = data.immunohistochemistry()
+    >>> ihc_hdx = separate_stains(ihc, hdx_from_rgb)
+    """
+    rgb = dtype.img_as_float(rgb) + 2
+    stains = np.dot(np.reshape(-np.log(rgb), (-1, 3)), conv_matrix)
+    return np.reshape(stains, rgb.shape)
+
+
+def combine_stains(stains, conv_matrix):
+    """Stain to RGB color space conversion.
+
+    Parameters
+    ----------
+    stains : array_like
+        The image in stain color space, in a 3-D array of shape (.., .., 3).
+    conv_matrix: ndarray
+        The stain separation matrix as described by G. Landini [1]_.
+
+    Returns
+    -------
+    out : ndarray
+        The image in RGB format, in a 3-D array of shape (.., .., 3).
+
+    Raises
+    ------
+    ValueError
+        If `stains` is not a 3-D array of shape (.., .., 3).
+
+    Notes
+    -----
+    Stain combination matrices available in the ``color`` module and their
+    respective colorspace:
+
+    * ``rgb_from_hed``: Hematoxylin + Eosin + DAB
+    * ``rgb_from_hdx``: Hematoxylin + DAB
+    * ``rgb_from_fgx``: Feulgen + Light Green
+    * ``rgb_from_bex``: Giemsa stain : Methyl Blue + Eosin
+    * ``rgb_from_rbd``: FastRed + FastBlue +  DAB
+    * ``rgb_from_gdx``: Methyl Green + DAB
+    * ``rgb_from_hax``: Hematoxylin + AEC
+    * ``rgb_from_bro``: Blue matrix Anilline Blue + Red matrix Azocarmine\
+                        + Orange matrix Orange-G
+    * ``rgb_from_bpx``: Methyl Blue + Ponceau Fuchsin
+    * ``rgb_from_ahx``: Alcian Blue + Hematoxylin
+    * ``rgb_from_hpx``: Hematoxylin + PAS
+
+    References
+    ----------
+    .. [1] http://www.dentistry.bham.ac.uk/landinig/software/cdeconv/cdeconv.html
+
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.color import (separate_stains, combine_stains,
+    ...                            hdx_from_rgb, rgb_from_hdx)
+    >>> ihc = data.immunohistochemistry()
+    >>> ihc_hdx = separate_stains(ihc, hdx_from_rgb)
+    >>> ihc_rgb = combine_stains(ihc_hdx, rgb_from_hdx)
+    """
+    from ..exposure import rescale_intensity
+
+    stains = dtype.img_as_float(stains)
+    logrgb2 = np.dot(-np.reshape(stains, (-1, 3)), conv_matrix)
+    rgb2 = np.exp(logrgb2)
+    return rescale_intensity(np.reshape(rgb2 - 2, stains.shape), in_range=(-1, 1))

skimage/color/tests/test_colorconv.py

@@ -16,11 +16,14 @@ import os.path
 import numpy as np
 from numpy.testing import *
 
-from skimage import img_as_float
+from skimage import img_as_float, img_as_ubyte
 from skimage.io import imread
 from skimage.color import (
     rgb2hsv, hsv2rgb,
     rgb2xyz, xyz2rgb,
+    rgb2hed, hed2rgb,
+    separate_stains,
+    combine_stains,
     rgb2rgbcie, rgbcie2rgb,
     convert_colorspace,
     rgb2grey, gray2rgb,
@@ -121,6 +124,32 @@ class TestColorconv(TestCase):
         img_rgb = img_as_float(self.img_rgb)
         assert_array_almost_equal(xyz2rgb(rgb2xyz(img_rgb)), img_rgb)
 
+    # RGB<->HED roundtrip with ubyte image
+    def test_hed_rgb_roundtrip(self):
+        img_rgb = self.img_rgb
+        assert_equal(img_as_ubyte(hed2rgb(rgb2hed(img_rgb))), img_rgb)
+
+    # RGB<->HED 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 ubyte image
+    def test_hdx_rgb_roundtrip(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)
+
     # RGB to RGB CIE
     def test_rgb2rgbcie_conversion(self):
         gt = np.array([[[ 0.1488856 ,  0.18288098,  0.19277574],

skimage/data/__init__.py

@@ -127,3 +127,19 @@ def clock():
 
     """
     return load("clock_motion.png")
+
+
+def immunohistochemistry():
+    """Immunohistochemical (IHC) staining with hematoxylin counterstaining.
+
+    This picture shows colonic glands where the IHC expression of FHL2 protein
+    is revealed with DAB. Hematoxylin counterstaining is applied to enhance the
+    negative parts of the tissue.
+
+    This image was acquired at the Center for Microscopy And Molecular Imaging
+    (CMMI).
+
+    No known copyright restrictions.
+
+    """
+    return load("ihc.jpg")