目录
颜色反卷积算法的设计针对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里面先用:图像-调整-曲线-自动,增加灰度值,再图像-调整-反相
![]()
![]()
![]()
- 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)
- #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()
-
- #!/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()
-
-






看起来是倒数第二张图结果比较好,最后在PS里面先用:图像-调整-曲线-自动,增加灰度值,再图像-调整-反相,得到我们需要的结果:
