dye_unmixing_algorithm.py

# -*- coding: utf-8 -*-
"""Colin_images_analysis.ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1wYV-FoqlWBQKdKniLevRz5PBweXaWqir
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from skimage import io
import scipy.io as sio
from scipy.optimize import minimize
import os
from PIL import Image
import torch
from torchvision import transforms
from torchvision.utils import save_image
import torch.nn.functional as F
import sys

from sklearn.preprocessing import normalize

import cv2

# pip install histomicstk --find-links https://girder.github.io/large_image_wheels --quiet

# import histomicstk as htk

image = io.imread('sample_image.jpg')
image.shape

# convert RGB to OD https://github.com/Peter554/StainTools/blob/master/staintools/utils/optical_density_conversion.py
# np.maximum(-1 * np.log(I / 255), 1e-6)

od_image = np.maximum(-1 * np.log(image / 255), 1e-6)
od_image.shape

image[:2,:2,0]

od_image[:2,:2,0]

# stain matrix 4 stains
# Hx, Eo, DAB, MG
mat_4stain = np.array([[0.65,0.70,0.28],[0.07,0.99,0.10],[0.26,0.57,0.77],[0.98,0.14,0.13]])

# stain matrix stain_row_wise x RGB_col_wise
mat_4stain

# normalize RGB OD to 1
normalized_mat_4stain = normalize(mat_4stain, axis=1, norm='l1')

normalized_mat_4stain

# for linear algebra calculations transpose to RGB_row_wise X stain_col_wise
normalized_transposed_mat_4stain = np.transpose(normalized_mat_4stain)
normalized_transposed_mat_4stain.shape

normalized_transposed_mat_4stain

# calculate norm of normalized stain matrix
np.linalg.norm(normalized_mat_4stain), np.linalg.norm(normalized_transposed_mat_4stain)

# reshape the od image
reshaped_od_image = np.reshape(od_image,(-1,3))

image.shape, od_image.shape, reshaped_od_image.shape

# optimization algorithm
def unmix_stain4(data, H):
  b = data
  # calculate the minimum distance between OD values and chosen x values
  fun = lambda x: np.linalg.norm(np.dot(H, x) - b)
  #fun = lambda x: x[0]
  sol = minimize(fun, np.asarray([0.5, 0.5, 0.5,0.5]), method='SLSQP')
  # print(f"sol shape is {sol['x'].shape}")
  #return np.array([sol['x'][0], sol['x'][1], sol['x'][2]],sol['x'][3])
  return sol['x']

np.asarray([0.5, 0.5, 0.5,0.5]).shape

def lu_4stain(data, H,chann):
    data = np.reshape(np.array(data), (-1, chann))
    print(f'reshaped image data {data.shape}')
    def perform_spectral_unmixing(_data):
      # zeros matrix with same size as image and 4 columns
      res = np.zeros((_data.shape[0],4))
      print(f'res matrix shape {res.shape}')
      for idx in range(len(_data)):
        # print(f'input data shape {data[idx].shape}')
        res[idx,:] = unmix_stain4(data[idx], H) # fill the stain concentration values in the zeros matrix
        # print(f"unmix output {res.shape}")
        #break
      return res
    lu_result = perform_spectral_unmixing(data)
    print(f"unmixed output data {lu_result.shape}")
    print('computation finished')
    return np.reshape(lu_result, (img_size, img_size, 4))
    #return lu_result

img_size = od_image.shape[0]

lu_result4 = lu_4stain(od_image, normalized_transposed_mat_4stain,3)

lu_result4.shape

fig,ax = plt.subplots(1,5,figsize = (10,10))
ax = ax.ravel()
ax[0].imshow(image)
ax[0].set_title('image')
ax[1].imshow(lu_result4[:,:,0],cmap="gray")
ax[1].set_title('Hx')
ax[2].imshow(lu_result4[:,:,1],cmap="gray")
ax[2].set_title('Eo')
ax[3].imshow(lu_result4[:,:,2],cmap="gray")
ax[3].set_title('DAB')
ax[4].imshow(lu_result4[:,:,3],cmap="gray")
ax[4].set_title('MG')
plt.tight_layout()
plt.show()

Hx_conc = lu_result4[:,:,0]
Eo_conc = lu_result4[:,:,1]
DAB_conc = lu_result4[:,:,2]
MG_conc = lu_result4[:,:,3]

dye_concentration_per_pixel_dict = {'Hx_conc':Hx_conc,'Eo_conc':Eo_conc,'DAB_conc':DAB_conc,'MG_conc':MG_conc}

Hx_conc = Hx_conc.reshape(-1)
Hx_conc.shape

Hx_absorption_factor = normalized_mat_4stain[0,:]

Hx_abs_fac_columns = Hx_absorption_factor.reshape(-1,1).T
Hx_abs_fac_columns.shape

Hx_conc = Hx_conc.reshape(-1,1)
Hx_conc.shape

Hx_RGB_OD = Hx_conc*Hx_abs_fac_columns
Hx_RGB_OD.shape

# reshape it in the size of the image
Hx_RGB_OD = np.reshape(Hx_RGB_OD, (img_size, img_size, 3))
Hx_RGB_OD.shape

plt.figure(figsize = (2,2))
plt.imshow(Hx_RGB_OD)

# convert the dye concentration* absorption factor to ratio of intensites
Hx_RGB_intensities = np.power(10,-(Hx_conc*Hx_abs_fac_columns))

# convert ratio of intensities to RGB
Hx_RGB= (Hx_RGB_intensities*255).astype('uint8')

# reshape it in the size of the image
Hx_RGB = np.reshape(Hx_RGB, (img_size, img_size, 3))
Hx_RGB.shape

plt.figure(figsize = (2,2))
plt.imshow(Hx_RGB)

"""### generalizing to all stains"""

Hx_absorption_factor = normalized_mat_4stain[0,:]
Eo_absorption_factor = normalized_mat_4stain[1,:]
DAB_absorption_factor = normalized_mat_4stain[2,:]
MG_absorption_factor = normalized_mat_4stain[3,:]

def convert_dye_concentration_to_RGB(dye_concentration_dict, stain_matrix,image_size):
  normalized_stain_matrix = normalize(stain_matrix, axis=1, norm='l1') # normalize across RGB
  Hx_absorption_factor = normalized_stain_matrix[0,:] # hematoxylin specific absorption factor
  Eo_absorption_factor = normalized_stain_matrix[1,:] # eosin specific absorption factor
  DAB_absorption_factor = normalized_stain_matrix[2,:] # DAB specific absorption factor
  MG_absorption_factor = normalized_stain_matrix[3,:]  # methyl green specific absorption factor
  abs_factor_colums = 0
  for pixel_stain in dye_concentration_dict:
    pixel_dye_concentration = dye_concentration_dict[pixel_stain].reshape(-1)
    pixel_dye_concentration = pixel_dye_concentration.reshape(-1,1) # create a matrix of shape num_pixel_as_rows X 1
    #print(pixel_stain)
    #print(pixel_dye_concentration.shape)
    # OD = dye_concentration X absorption factor
    if pixel_stain == "Hx_conc":
      abs_factor_columns = Hx_absorption_factor.reshape(-1,1).T # reshape asborption factor as 1 X RGB
      Hx_RGB_OD = pixel_dye_concentration*abs_factor_columns    # calculate RGB OD based on dye concentration
      # convert the dye concentration* absorption factor to ratio of intensites
      Hx_RGB_intensities = np.power(10,-Hx_RGB_OD)
      # convert ratio of intensities to RGB
      Hx_RGB = (Hx_RGB_intensities*255).astype('uint8')
      # reshape it in the size of the image
      Hx_RGB = np.reshape(Hx_RGB, (img_size, img_size, 3))
    elif pixel_stain == "Eo_conc":
      abs_factor_columns = Eo_absorption_factor.reshape(-1,1).T # reshape asborption factor as 1 X RGB
      Eo_RGB_OD = pixel_dye_concentration*abs_factor_columns
      # convert the dye concentration* absorption factor to ratio of intensites
      Eo_RGB_intensities = np.power(10,-Eo_RGB_OD)
      # convert ratio of intensities to RGB
      Eo_RGB = (Eo_RGB_intensities*255).astype('uint8')
      # reshape it in the size of the image
      Eo_RGB = np.reshape(Eo_RGB, (img_size, img_size, 3))
    elif pixel_stain == "DAB_conc":
      abs_factor_columns = DAB_absorption_factor.reshape(-1,1).T # reshape asborption factor as 1 X RGB
      DAB_RGB_OD = pixel_dye_concentration*abs_factor_columns
      # convert the dye concentration* absorption factor to ratio of intensites
      DAB_RGB_intensities = np.power(10,-DAB_RGB_OD)
      # convert ratio of intensities to RGB
      DAB_RGB = (DAB_RGB_intensities*255).astype('uint8')
      # reshape it in the size of the image
      DAB_RGB = np.reshape(DAB_RGB, (img_size, img_size, 3))
    elif pixel_stain == "MG_conc":
      abs_factor_columns = MG_absorption_factor.reshape(-1,1).T # reshape asborption factor as 1 X RGB
      MG_RGB_OD = pixel_dye_concentration*abs_factor_columns
      # convert the dye concentration* absorption factor to ratio of intensites
      MG_RGB_intensities = np.power(10,-MG_RGB_OD)
      # convert ratio of intensities to RGB
      MG_RGB = (MG_RGB_intensities*255).astype('uint8')
      # reshape it in the size of the image
      MG_RGB = np.reshape(MG_RGB, (img_size, img_size, 3))
    else:
      print("desired stain concentration not found")
      break
  return (Hx_RGB, Eo_RGB ,DAB_RGB ,MG_RGB )
  #return concentration*absorption_factor

Hx_RGB, Eo_RGB ,DAB_RGB ,MG_RGB = convert_dye_concentration_to_RGB(dye_concentration_per_pixel_dict,mat_4stain,img_size)

fig,ax = plt.subplots(1,5,figsize = (10,10))
ax = ax.ravel()
ax[0].imshow(image)
ax[0].set_title('image')
ax[1].imshow(Hx_RGB)
ax[1].set_title('Hx')
ax[2].imshow(Eo_RGB)
ax[2].set_title('Eo')
ax[3].imshow(DAB_RGB)
ax[3].set_title('DAB')
ax[4].imshow(MG_RGB)
ax[4].set_title('MG')
plt.tight_layout()
plt.show()

print(f"Hx conc {np.linalg.norm(Hx_conc)}")
print(f"Eo conc {np.linalg.norm(Eo_conc)}")
print(f"DAB conc {np.linalg.norm(DAB_conc)}")
print(f"MG conc {np.linalg.norm(MG_conc)}")