RGB to HED conversion after version 0.17

[nicolebussola]

Description

After upgrading our project environment to version 0.18.x we have noticed that the RGB to HED conversion behaves unexpectedly according to the biological meaning of the HED color space. We have done some research and carefully read the PR changing the stain deconvolution method, however the HED conversion in version 0.17.x seemed to be more appropriate; in particular the H channel highlights the nuclei while the DAB channel highlights the brown stain:

version 0.17.x

We are not able to find the same biological meaning in version 0.18.x:

version 0.18.x

Regards,
Nicole, Alessia, Ernesto (histolab team)

[emmanuelle]

Hi @nicolebussola thanks for reporting this. It's true that the doc example https://scikit-image.org/docs/dev/auto_examples/color_exposure/plot_ihc_color_separation.html does not look great. @alexdesiqueira since you worked on this PR could you please comment?

[alexdesiqueira]

Hi @nicolebussola, thank you for opening this issue. Yeah @emmanuelle, the doc example doesn't "feel right" anymore.
I will need some time to check that, Nicole. As you saw, there was a lot of research — and people! — involved in #4725, and according to the papers, the solution we have now seem appropriate. However, as seen, something doesn't feel right...
My course of action here would be:

• Check out @crisluengo's implementation in DIPLib. I'll run a couple of tests using their library; let's see their results using the same images. Maybe there's something missing in our latest implementation.
• Check if there's something we need to improve that example. I'd bet that a better colormap is needed to present the results. Not sure, though.

That said, how does that sound, Nicole? Is there something I'm missing, or would you like to add something else?
Thank you again!

[crisluengo]

In skimage.color.separate_stains() we're missing a clip operation. Adding ihc_hed[ihc_hed<0] = 0 before display corrects things. I missed this in our rewrite of these functions, my apologies.

Current code for separate_stains:

rgb = _prepare_colorarray(rgb, force_copy=True)
np.maximum(rgb, 1E-6, out=rgb)  # avoiding log artifacts
log_adjust = np.log(1E-6)  # used to compensate the sum above
stains = (np.log(rgb) / log_adjust) @ conv_matrix

In the last line, np.log(rgb) / log_adjust produces a correct image in the range [0,1]. But the rotation (@ conv_matrix) can lead to negative values. These negative values spoil the fun (actually spoil the display, which has a min-max stretch built in). Adding

np.minimum(stains, 0, out= stains)

should fix the function.

If we create the input artificially by combining positive values and using the colors as assumed in conv_matrix, then the function will be guaranteed to produce only positive values, but the real world is not ideal: the stain colors are not the same as those given in conv_matrix, and the image has noise also.

Yes, this clipping makes the round-trip not idempotent. But if you start with a synthetic image, the clip will do nothing and hence the tests should continue to work. On the other hand,

ihc_hed = rgb2hed(ihc_rgb)
reconstructed = hed2rgb(ihc_hed)

shows in reconstructed which colors for each stain are chosen in conv_matrix (actually hed_from_rgb). These colors don't exactly match the data: the DAB is more red, less yellow, and the hematoxylin is somewhat more purple. This also shows that the unmixing is inexact: the only way to get the exact colors is to measure them specifically for the microscope and chemistry sets used. Every lab produces different colors for the same stains...

[crisluengo]

Also, the example at https://scikit-image.org/docs/dev/auto_examples/color_exposure/plot_ihc_color_separation.html has a bug:

cmap_eosin = LinearSegmentedColormap.from_list('mycmap', ['darkviolet', 'white'])

The two colors are reversed, this should be:

cmap_eosin = LinearSegmentedColormap.from_list('mycmap', ['white', 'darkviolet'])

Nobody noticed this before because the image has no Eosin, and so this channel has only residue (= negative values). Reversing the colormap actually looks good in this case, but the image should be all zeros (white), which we do see after clipping negative values.

Finally, the DAB channel (brown) looks rather dim. This is because there is one pixel that happens to be really dense (maybe noise?) and the image is scaled to that. But this is just a display issue. If you want to avoid this, you can do like in the ImageJ plugin, where they don't use 1e-6 but 1/255 in the initial clip. This causes the darkest stain to be less dark. This is a question about what you prefer: better display or better data.

[crisluengo]

OK, one last comment from me and then I'll shut up. :)

I'd change the demo as follows for proper display of the unmixed stains. What we're doing here is creating an image in "hed color space" that has only one stain, then using hed2rgb to produce a correct representation of that stain. There is now no worry about scaling of the images. The result of hed2rgb should be comparable in intensities to the original input.

import matplotlib.pyplot as plt
import numpy as np

from skimage import data
from skimage.color import rgb2hed, hed2rgb


ihc_rgb = data.immunohistochemistry()
ihc_hed = rgb2hed(ihc_rgb)
ihc_hed[ihc_hed<0] = 0  # TODO: remove this after fixing `rgb2hed()`

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))

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_h)
ax[1].set_title("Hematoxylin")

ax[2].imshow(ihc_e)
ax[2].set_title("Eosin")

ax[3].imshow(ihc_d)
ax[3].set_title("DAB")

for a in ax.ravel():
    a.axis('off')

fig.tight_layout()
plt.show()

It is still fine to use ihc_hed[:,:,0] etc. as in the second part of the example.

[grlee77]

Thanks for sorting this out @crisluengo! This is looking much better and the fix is pretty straightforward. Let us know if you want to make a PR, otherwise we can take care of it.

[alexdesiqueira]

Hi @nicolebussola, please let us know if the current code solves these issues. Thank you again @crisluengo!