目录
介绍
颜色反卷积算法的设计针对RGB摄像机获取的颜色信息，基于免疫组化技术使用的染色剂RGB分量光的特异性吸收，分别计算每种染色剂对图像的作用效果。免疫组织化学图像处理通常用的染色包括DAB、H&E。
颜色反卷积可应用于H-DAB图像和组织病理学中常用的染色剂（H-E,H AEC,FAST Red，DAB），广泛应用于免疫组化染色图像颜色分离。
环境
#### Language de programmation : Python 2.7
#### Libs : numpy, matplotlib, sikit-learn, PySide, OpenCV2
运行步骤
1.运行ColorDeconvolution.py或者hsd.py，可以生成反卷积后的图像（一共有6个图），依据情况选择结果。
2.得到的HSI_Teinte_t1.png图像，灰度值比较低，在PS里面先用：图像-调整-曲线-自动，增加灰度值，再图像-调整-反相
- 一共三个文件，放在一个文件夹即可
OpticalDensity.py
import numpy as np
from matplotlib import pyplot as plt
def rgb_2_od(img):
s = img.shape
od = np.zeros(s)
if (len(img)==3):
for i in range(s[0]):
for j in range(s[1]):
for k in range(s[2]):
if img[i, j, k] != 0:
od[i, j, k] = np.log10(img[i, j, k])
else:
for i in range(s[0]):
for j in range(s[1]):
od[i, j] = np.log10(img[i, j])
plt.title("Espace OD")
plt.imshow(od)
plt.show()
return od
def rgb_2_od2(img):
s = img.shape
od = np.zeros(s)
if (len(img)==3):
for i in range(s[0]):
for j in range(s[1]):
for k in range(s[2]):
if img[i, j, k] != 0:
od[i, j, k] = np.log10(img[i, j, k])
else:
for i in range(s[0]):
for j in range(s[1]):
od[i, j] = np.log10(img[i, j])
return od
if __name__ == "__main__":
# reading the image
# path="Tumor_CD31_LoRes.png"
path = "../DataSet/BreastCancerCell_dataset/ytma10_010704_benign1_ccd.tif"
img = plt.imread(path)
img = img[:, :, 0:3].astype(np.double)
rgb_2_od(img)hsd.py
#coding: utf-8
import numpy as np
import matplotlib.pyplot as plt
import OpticalDensity as od
from sklearn.cluster import KMeans
from matplotlib import cm
import cv2
from PIL import Image
class HSD:
def __init__(self, img,path,type):
"""
:param img: image en entrer pour changer d'espace de couleur
chroma : matrice de chromaticité. de dimension
[
l: nbr de ligne ,
c: nbr de colone,
d: dimension 2 qui represente les coordonée chromatique
]
:fonctions :
- rgb_2_gray : transfomé en gris
- chromaticite : calculer les coordonnées chromatiques
- calcule_HSI : calculer le nouvel espace Hue Saturation Intensity
"""
self.img_0 = img
self.img_hsi = None
self.chroma = None
self.imgReconstruite=None
self.intensity=None
self.path=path
self.path2="../Color-Deconvolution/Resultat"
self.type=type
self.imageBinariser=None
def getRed(self):
return self.img_0[:,:,0]
def getGreen(self):
return self.img_0[:,:,1]
def getBlue(self):
return self.img_0[:,:,2]
def GlobalIntensity(self):
self.intensity = np.zeros([self.img_0.shape[0], self.img_0.shape[1]])
for i in range(self.img_0.shape[0]):
for j in range(self.img_0.shape[1]):
self.intensity[i, j] = (self.getRed()[i,j]+ self.getGreen()[i,j] + self.getBlue()[i,j]) / 3
def RGB_2_GRAY(self):
[l, c, d] = self.img_0.shape
gray = np.zeros([l, c])
for i in range(l):
for j in range(c):
gray[i, j] = sum(self.img_0[i, j])/3
return gray
def chromaticite(self):
[l, c, d] = self.img_0.shape
gray = self.RGB_2_GRAY()
self.chroma = np.zeros([l, c, 2])
for i in range(l):
for j in range(c):
self.chroma[i, j ,0] = (self.img_0[i, j, 0] / gray[i, j]) - 1 if (gray[i, j] - 1) != 0 else 0
self.chroma[i, j, 1] = (self.img_0[i, j, 1] - self.img_0[i, j, 2]) / (gray[i, j] * np.sqrt(3)) if (gray[i, j] * np.sqrt(3)) !=0 else 0
def calcule_HSI(self):
[l, c, d] = self.img_0.shape
gray = self.RGB_2_GRAY()
self.img_hsi = np.zeros([l, c, 3])
self.img_hsi[:, :, 2] = gray
for i in range(l):
for j in range(c):
self.img_hsi[i, j, 0] = self.getHue2(self.chroma[i, j, 0], self.chroma[i, j ,1])
self.img_hsi[i, j, 1] = self.getSaturation2(self.chroma[i, j ,0], self.chroma[i, j ,1])
def getSaturation2(self, cx, cy):
return np.sqrt(np.square(cx)+np.square(cy))
def getHue2(self, cx, cy):
return 0 if cx==0 else np.arctan(cy/cx)
def getCX(self,hue,saturation):
return saturation*np.cos(hue)
def getCY(self,hue,saturation):
return saturation*np.sin(hue)
# plotHSD permet de afficher les images de chaque canal de la hsd en gray
def plotHSD(self):
plt.subplot(1,3,1)
plt.imshow(self.img_hsi[:,:,0], cmap="gray")
plt.subplot(1,3,2)
plt.imshow(self.img_hsi[:, :, 1], cmap="gray")
plt.subplot(1, 3, 3)
plt.imshow(self.img_hsi[:, :, 2], cmap="gray")
plt.show()
# Afficher le plot manuellement selon les valeurs de chroma
def plotHS(self,Z):
for i in range(Z.shape[0]):
for j in range(Z.shape[1]):
if Z[i][j]==0:
c1 = plt.scatter(self.chroma[i,j,0], self.chroma[i,j, 1], c='r', marker='+')
elif Z[i][j]==1:
c2 = plt.scatter(self.chroma[i, j, 0], self.chroma[i, j, 1], c='b', marker='o')
elif Z[i][j]==2:
c3 = plt.scatter(self.chroma[i, j, 0], self.chroma[i, j, 1], c='g', marker='*')
elif Z[i][j]==3:
c4 = plt.scatter(self.chroma[i, j, 0], self.chroma[i, j, 1], c='black', marker='l')
plt.legend([c1, c2, c3], ['Cluster 0', 'Cluster 1','Cluster 2'])
plt.title('K-means clusters into 3 clusters')
plt.show()
def plotHS2(self,Z):
Z=np.reshape(Z,(self.chroma.shape[0]*self.chroma.shape[1],1))
print (Z.shape)
c1=0
vec = np.reshape(self.chroma, (self.chroma.shape[0] * self.chroma.shape[1], 2))
vec1 = []
vec2 = []
vec3 = []
for i in range(Z.shape[0]):
if Z[i]==0:
vec1.append(vec[i])
elif Z[i]==1:
vec2.append(vec[i])
elif Z[i]==2:
vec3.append(vec[i])
# print len(vec1)
c1 = plt.scatter(vec1[1:,0], vec1[1:,1], c='r', marker='+')
c2 = plt.scatter(vec2[1:,0], vec2[1:,1], c='g', marker='o')
c3 = plt.scatter(vec3[1:,0], vec3[1:,1], c='b', marker='*')
plt.legend([c1, c2, c3], ['Cluster 0', 'Cluster 1', 'Cluster 2'])
plt.title('K-means clusters into 3 clusters')
plt.show()
print ("finish")
def Kmeans(self,cluster):
np.random.seed(42)
vec=np.reshape(self.chroma, (self.chroma.shape[0]*self.chroma.shape[1],2))
kmeans = KMeans(n_clusters=cluster, random_state=0)
kmeans.fit(vec)
kmeans.labels_=np.reshape(kmeans.labels_,(self.chroma.shape[0],self.chroma.shape[1]))
Z=np.reshape(kmeans.labels_,self.chroma.shape[0]*self.chroma.shape[1],2)
Z=np.reshape(Z[:100],(10,10))
return Z
def Kmeans2(self,cluster):
np.random.seed(42)
vec=np.reshape(self.chroma, (self.chroma.shape[0]*self.chroma.shape[1],2))
kmeans = KMeans(n_clusters=cluster, random_state=0)
kmeans.fit(vec)
plt.scatter(vec[:1000, 0], vec[:1000, 1], c=kmeans.labels_[:1000])
plt.show()
def binarisation(self):
# Otsu's thresholding after Gaussian filtering
# il faut lui donner image self.imgReconstruite et binariser chaque canal hue saturation et density
blur1 = cv2.GaussianBlur((self.imgReconstruite[:, :, 0]*255).astype(np.uint8), (5, 5), 0)
ret1, th1 = cv2.threshold(blur1, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
blur2 = cv2.GaussianBlur((self.imgReconstruite[:, :, 1]*255).astype(np.uint8), (5, 5), 0)
ret2, th2 = cv2.threshold(blur2, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
blur3 = cv2.GaussianBlur((self.imgReconstruite[:, :, 2]*255).astype(np.uint8), (5, 5), 0)
ret3, th3 = cv2.threshold(blur3, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
plt.subplot(1,3,1)
plt.title("hue")
plt.imshow(th1,cmap="gray")
plt.subplot(1, 3, 2)
plt.title("Saturation")
plt.imshow(th2,cmap="gray")
plt.subplot(1, 3, 3)
plt.title(self.type)
plt.imshow(th3,cmap="gray")
plt.show()
self.saveOpencv(th1,"Hue_Binariser.png")
self.saveOpencv(th2, "Saturation_Binariser.png")
self.saveOpencv(th3, "Intensity_Binariser.png")
return th3
def recontructionToRGB(self):
#Calcul de intensity global de chaque pixel on en a besoin pour la reconstruction
self.GlobalIntensity()
[l,c,d]=self.img_hsi.shape
self.imgReconstruite=np.zeros([l,c,3])
for i in range(l):
for j in range(c):
cx = self.getCX(self.img_hsi[i, j, 0],self.img_hsi[i, j, 1])
cy = self.getCY(self.img_hsi[i, j, 0],self.img_hsi[i, j, 1])
self.imgReconstruite[i, j, 0] = self.intensity[i, j] * (cx + 1)
self.imgReconstruite[i, j, 1] = 0.5 * self.intensity[i, j] * (2 - cx - np.sqrt(3) * cy)
self.imgReconstruite[i, j, 2] = 0.5 * self.intensity[i, j] * (2 - cx - np.sqrt(3) * cy)
plt.subplot(2,3,1)
plt.title("hue")
plt.imshow(self.imgReconstruite[:, :, 0], cmap="gray")
# self.saveOpencv((self.imgReconstruite[:, :, 0]*255).astype(np.uint8), "HSI_Teinte.png")
cv2.imwrite("HSI_Teinte.png",(self.imgReconstruite[:,:,0]*255).astype(np.uint8))
#plt.colorbar(ticks=[0, 60, 120, 179], orientation='horizontal', cmap=cm.hsv)
plt.subplot(2,3,2)
plt.title("saturation")
plt.imshow(self.imgReconstruite[:, :, 1], cmap="gray")
# save HSI
cv2.imwrite("HSI_Saturation.png",(self.imgReconstruite[:,:,1]*255).astype(np.uint8))
# self.saveOpencv((self.imgReconstruite[:, :, 1]*255).astype(np.uint8), "HSI_Saturation.png")
plt.subplot(2,3,3)
plt.title(self.type)
plt.imshow(self.imgReconstruite[:, :, 2], cmap="gray")
# save HSD
#self.savePillow((self.imgReconstruite[:,:,2]*255).astype(np.uint8), self.path+"HSD_Density.tif")
cv2.imwrite("HSI_Density.png",(self.imgReconstruite[:,:,2]*255).astype(np.uint8))
# self.saveOpencv((self.imgReconstruite[:,:,2]*255).astype(np.uint8), "HSI_Density.png")
#img = Image.fromarray(self.img_0[:, :, 0])
#img.save(self.path + "img.tif")
plt.show()
def saveOpencv(self,img,path):
cv2.imwrite(self.path2+"/"+self.path[19:]+"/"+path, img)
def savePillow(self,img,path):
img_to_save = Image.fromarray(img)
img_to_save.save(self.path2+"/"+self.path[19:]+"/"+path)
if __name__ == "__main__":
# reading the image
# path="Tumor_CD31_LoRes.png"
#path = "../DataSet/BreastCancerCell_dataset/ytma10_010704_benign1_ccd.tif"
# path = "../DataSet_Lomenie/tumor.png"
path="t2.png"
img = plt.imread(path)
img = img[:, :, 0:3].astype(np.double)
od=od.rgb_2_od2(img)
# hsi
h = HSD(od, path, "Density")
h.chromaticite()
h.calcule_HSI()
# h.Kmeans2(3)
# h.plotHS(Z)
h.recontructionToRGB()
h.binarisation()
ColorDeconvolution.py
#!/usr/bin/python
import matplotlib.pyplot as plt
import numpy as np
import cv2
from PIL import Image
import hsd
from matplotlib import cm
class ColorDeconvolution:
def __init__(self, img,path):
self.img_0 = img
self.stains = None
self.od = None
self.path=path
self.path2 = "autodl-tmp/color-deconv/Color-Deconvolution/Resultat"
def getStains(self):
return self.stains
def setImage(self, img):
if img != None:
self.img_0 = img
def RGB_2_OD(self):
[l, c, d] = self.img_0.shape
self.od = np.zeros([l, c, d])
for i in range(l):
for j in range(c):
for k in range(d):
if self.img_0[i,j,k] != 0 :
self.od[i,j,k] = np.log(self.img_0[i, j, k])
return self.od
def norm(self, vector):
n = 0
for i in vector:
n = n + np.square(i)
return np.sqrt(n)
def separateStain(self):
# He = [0.65; 0.70; 0.29];
# Eo = [0.07; 0.99; 0.11];
# DAB = [0.27; 0.57; 0.78];
He = np.array([0.6500286, 0.704031, 0.2860126]) # Hematoxylin
Eo = np.array([0.07, 0.99, 0.11]) # Eosine
DAB = np.array([0.26814753, 0.57031375, 0.77642715]) # DAB
Res = np.array([0.7110272, 0.42318153, 0.5615672]) # residual
HDABtoRGB = np.array([He / self.norm(He), Eo / self.norm(Eo), DAB / self.norm(DAB)])
RGBtoHDAB = np.linalg.inv(HDABtoRGB)
[l,c,d] = self.img_0.shape
self.stains = np.zeros([l, c, d])
for i in range(l):
for j in range(c):
a = np.dot(self.od[i, j], RGBtoHDAB)
b=self.od[i,j]
self.stains[i, j, 0] = a[0]
self.stains[i, j, 1] = a[1]
self.stains[i, j, 2] = a[2]
return self.stains
# this methods shows the differents stains of each canal for the method of colorDeconvolution
def binarisation(self,image):
# Otsu's thresholding after Gaussian filtering
blur = cv2.GaussianBlur(image, (5, 5), 0)
ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
return th3
def showStains(self):
plt.subplot(1, 4, 1)
plt.title("original")
plt.imshow(self.img_0)
plt.subplot(1, 4, 2)
plt.title('Hematoxylin')
plt.imshow(self.stains[:, :, 0], cmap="gray")
cv2.imwrite("colorDeconvolution_Hematoxylin.png",(self.stains[:, :, 0]*255).astype(np.uint8))
# self.saveOpencv((self.stains[:, :, 0]*255).astype(np.uint8), " colorDeconvolution_Hematoxylin.tif")
plt.subplot(1, 4, 3)
plt.title('Eosine')
plt.imshow(self.stains[:, :, 1], cmap="gray")
cv2.imwrite("colorDeconvolution_Eosin.png",(self.stains[:, :, 1]*255).astype(np.uint8))
# self.saveOpencv((self.stains[:, :, 1]*255).astype(np.uint8), " colorDeconvolution_Eosin.tif")
plt.subplot(1, 4, 4)
plt.title('DAB')
plt.imshow(self.stains[:, :, 2], cmap="gray")
cv2.imwrite(" colorDeconvolution_DAB.png",(self.stains[:, :, 2]*255).astype(np.uint8))
# self.saveOpencv((self.stains[:, :, 2]*255).astype(np.uint8)," colorDeconvolution_DAB.tif")
plt.show()
self.binarisation()
def saveOpencv(self,img,path):
cv2.imwrite(self.path2 + "/" + self.path[19:] + "/" + path, img)
def savePillow(self,img,path):
img_to_save = Image.fromarray(img)
img_to_save.save(path)
def binarisation(self):
# Otsu's thresholding after Gaussian filtering
# il faut lui donner image self.imgReconstruite et binariser chaque canal hue saturation et density
blur1 = cv2.GaussianBlur((self.stains[:, :, 0] * 255).astype(np.uint8), (5, 5), 0)
ret1, th1 = cv2.threshold(blur1, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
blur2 = cv2.GaussianBlur((self.stains[:, :, 1] * 255).astype(np.uint8), (5, 5), 0)
ret2, th2 = cv2.threshold(blur2, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
blur3 = cv2.GaussianBlur((self.stains[:, :, 2] * 255).astype(np.uint8), (5, 5), 0)
ret3, th3 = cv2.threshold(blur3, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
plt.subplot(1, 3, 1)
plt.title("Hematoxyline")
plt.imshow(th1,cmap="gray")
plt.subplot(1, 3, 2)
plt.title("Eosine")
plt.imshow(th2,cmap="gray")
plt.subplot(1, 3, 3)
plt.title("DAB")
plt.imshow(th3,cmap="gray")
plt.show()
self.saveOpencv(th1, "Hematoxyline.png")
self.saveOpencv(th2, "Eosine.png")
self.saveOpencv(th3, "DAB.png")
return th3
def saveOpencv(self,img, path):
cv2.imwrite(self.path2 + "/" + self.path[19:] + "/" + path, img)
if __name__=="__main__":
# path="../DataSet_Lomenie/arn2.tif"
path = "t2.png"
img = plt.imread(path)
img = img[:, :, 0:3]
# deconvolution de couleur
satin = ColorDeconvolution(img,path)
satin.RGB_2_OD()
image=satin.separateStain()
satin.showStains()
一共会有6张图的结果，最后一张是原图：
最后
看起来是倒数第二张图结果比较好，最后在PS里面先用：图像-调整-曲线-自动，增加灰度值，再图像-调整-反相，得到我们需要的结果：
1001
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
