• 病理图像反卷积


    目录

    介绍

    环境

    运行步骤

    OpticalDensity.py

     hsd.py

     ColorDeconvolution.py

     一共会有6张图的结果,最后一张是原图:

     最后

    源代码: 


    介绍

    颜色反卷积算法的设计针对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

    1. import numpy as np
    2. from matplotlib import pyplot as plt
    3. def rgb_2_od(img):
    4. s = img.shape
    5. od = np.zeros(s)
    6. if (len(img)==3):
    7. for i in range(s[0]):
    8. for j in range(s[1]):
    9. for k in range(s[2]):
    10. if img[i, j, k] != 0:
    11. od[i, j, k] = np.log10(img[i, j, k])
    12. else:
    13. for i in range(s[0]):
    14. for j in range(s[1]):
    15. od[i, j] = np.log10(img[i, j])
    16. plt.title("Espace OD")
    17. plt.imshow(od)
    18. plt.show()
    19. return od
    20. def rgb_2_od2(img):
    21. s = img.shape
    22. od = np.zeros(s)
    23. if (len(img)==3):
    24. for i in range(s[0]):
    25. for j in range(s[1]):
    26. for k in range(s[2]):
    27. if img[i, j, k] != 0:
    28. od[i, j, k] = np.log10(img[i, j, k])
    29. else:
    30. for i in range(s[0]):
    31. for j in range(s[1]):
    32. od[i, j] = np.log10(img[i, j])
    33. return od
    34. if __name__ == "__main__":
    35. # reading the image
    36. # path="Tumor_CD31_LoRes.png"
    37. path = "../DataSet/BreastCancerCell_dataset/ytma10_010704_benign1_ccd.tif"
    38. img = plt.imread(path)
    39. img = img[:, :, 0:3].astype(np.double)
    40. rgb_2_od(img)

     hsd.py

    1. #coding: utf-8
    2. import numpy as np
    3. import matplotlib.pyplot as plt
    4. import OpticalDensity as od
    5. from sklearn.cluster import KMeans
    6. from matplotlib import cm
    7. import cv2
    8. from PIL import Image
    9. class HSD:
    10. def __init__(self, img,path,type):
    11. """
    12. :param img: image en entrer pour changer d'espace de couleur
    13. chroma : matrice de chromaticité. de dimension
    14. [
    15. l: nbr de ligne ,
    16. c: nbr de colone,
    17. d: dimension 2 qui represente les coordonée chromatique
    18. ]
    19. :fonctions :
    20. - rgb_2_gray : transfomé en gris
    21. - chromaticite : calculer les coordonnées chromatiques
    22. - calcule_HSI : calculer le nouvel espace Hue Saturation Intensity
    23. """
    24. self.img_0 = img
    25. self.img_hsi = None
    26. self.chroma = None
    27. self.imgReconstruite=None
    28. self.intensity=None
    29. self.path=path
    30. self.path2="../Color-Deconvolution/Resultat"
    31. self.type=type
    32. self.imageBinariser=None
    33. def getRed(self):
    34. return self.img_0[:,:,0]
    35. def getGreen(self):
    36. return self.img_0[:,:,1]
    37. def getBlue(self):
    38. return self.img_0[:,:,2]
    39. def GlobalIntensity(self):
    40. self.intensity = np.zeros([self.img_0.shape[0], self.img_0.shape[1]])
    41. for i in range(self.img_0.shape[0]):
    42. for j in range(self.img_0.shape[1]):
    43. self.intensity[i, j] = (self.getRed()[i,j]+ self.getGreen()[i,j] + self.getBlue()[i,j]) / 3
    44. def RGB_2_GRAY(self):
    45. [l, c, d] = self.img_0.shape
    46. gray = np.zeros([l, c])
    47. for i in range(l):
    48. for j in range(c):
    49. gray[i, j] = sum(self.img_0[i, j])/3
    50. return gray
    51. def chromaticite(self):
    52. [l, c, d] = self.img_0.shape
    53. gray = self.RGB_2_GRAY()
    54. self.chroma = np.zeros([l, c, 2])
    55. for i in range(l):
    56. for j in range(c):
    57. self.chroma[i, j ,0] = (self.img_0[i, j, 0] / gray[i, j]) - 1 if (gray[i, j] - 1) != 0 else 0
    58. 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
    59. def calcule_HSI(self):
    60. [l, c, d] = self.img_0.shape
    61. gray = self.RGB_2_GRAY()
    62. self.img_hsi = np.zeros([l, c, 3])
    63. self.img_hsi[:, :, 2] = gray
    64. for i in range(l):
    65. for j in range(c):
    66. self.img_hsi[i, j, 0] = self.getHue2(self.chroma[i, j, 0], self.chroma[i, j ,1])
    67. self.img_hsi[i, j, 1] = self.getSaturation2(self.chroma[i, j ,0], self.chroma[i, j ,1])
    68. def getSaturation2(self, cx, cy):
    69. return np.sqrt(np.square(cx)+np.square(cy))
    70. def getHue2(self, cx, cy):
    71. return 0 if cx==0 else np.arctan(cy/cx)
    72. def getCX(self,hue,saturation):
    73. return saturation*np.cos(hue)
    74. def getCY(self,hue,saturation):
    75. return saturation*np.sin(hue)
    76. # plotHSD permet de afficher les images de chaque canal de la hsd en gray
    77. def plotHSD(self):
    78. plt.subplot(1,3,1)
    79. plt.imshow(self.img_hsi[:,:,0], cmap="gray")
    80. plt.subplot(1,3,2)
    81. plt.imshow(self.img_hsi[:, :, 1], cmap="gray")
    82. plt.subplot(1, 3, 3)
    83. plt.imshow(self.img_hsi[:, :, 2], cmap="gray")
    84. plt.show()
    85. # Afficher le plot manuellement selon les valeurs de chroma
    86. def plotHS(self,Z):
    87. for i in range(Z.shape[0]):
    88. for j in range(Z.shape[1]):
    89. if Z[i][j]==0:
    90. c1 = plt.scatter(self.chroma[i,j,0], self.chroma[i,j, 1], c='r', marker='+')
    91. elif Z[i][j]==1:
    92. c2 = plt.scatter(self.chroma[i, j, 0], self.chroma[i, j, 1], c='b', marker='o')
    93. elif Z[i][j]==2:
    94. c3 = plt.scatter(self.chroma[i, j, 0], self.chroma[i, j, 1], c='g', marker='*')
    95. elif Z[i][j]==3:
    96. c4 = plt.scatter(self.chroma[i, j, 0], self.chroma[i, j, 1], c='black', marker='l')
    97. plt.legend([c1, c2, c3], ['Cluster 0', 'Cluster 1','Cluster 2'])
    98. plt.title('K-means clusters into 3 clusters')
    99. plt.show()
    100. def plotHS2(self,Z):
    101. Z=np.reshape(Z,(self.chroma.shape[0]*self.chroma.shape[1],1))
    102. print (Z.shape)
    103. c1=0
    104. vec = np.reshape(self.chroma, (self.chroma.shape[0] * self.chroma.shape[1], 2))
    105. vec1 = []
    106. vec2 = []
    107. vec3 = []
    108. for i in range(Z.shape[0]):
    109. if Z[i]==0:
    110. vec1.append(vec[i])
    111. elif Z[i]==1:
    112. vec2.append(vec[i])
    113. elif Z[i]==2:
    114. vec3.append(vec[i])
    115. # print len(vec1)
    116. c1 = plt.scatter(vec1[1:,0], vec1[1:,1], c='r', marker='+')
    117. c2 = plt.scatter(vec2[1:,0], vec2[1:,1], c='g', marker='o')
    118. c3 = plt.scatter(vec3[1:,0], vec3[1:,1], c='b', marker='*')
    119. plt.legend([c1, c2, c3], ['Cluster 0', 'Cluster 1', 'Cluster 2'])
    120. plt.title('K-means clusters into 3 clusters')
    121. plt.show()
    122. print ("finish")
    123. def Kmeans(self,cluster):
    124. np.random.seed(42)
    125. vec=np.reshape(self.chroma, (self.chroma.shape[0]*self.chroma.shape[1],2))
    126. kmeans = KMeans(n_clusters=cluster, random_state=0)
    127. kmeans.fit(vec)
    128. kmeans.labels_=np.reshape(kmeans.labels_,(self.chroma.shape[0],self.chroma.shape[1]))
    129. Z=np.reshape(kmeans.labels_,self.chroma.shape[0]*self.chroma.shape[1],2)
    130. Z=np.reshape(Z[:100],(10,10))
    131. return Z
    132. def Kmeans2(self,cluster):
    133. np.random.seed(42)
    134. vec=np.reshape(self.chroma, (self.chroma.shape[0]*self.chroma.shape[1],2))
    135. kmeans = KMeans(n_clusters=cluster, random_state=0)
    136. kmeans.fit(vec)
    137. plt.scatter(vec[:1000, 0], vec[:1000, 1], c=kmeans.labels_[:1000])
    138. plt.show()
    139. def binarisation(self):
    140. # Otsu's thresholding after Gaussian filtering
    141. # il faut lui donner image self.imgReconstruite et binariser chaque canal hue saturation et density
    142. blur1 = cv2.GaussianBlur((self.imgReconstruite[:, :, 0]*255).astype(np.uint8), (5, 5), 0)
    143. ret1, th1 = cv2.threshold(blur1, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    144. blur2 = cv2.GaussianBlur((self.imgReconstruite[:, :, 1]*255).astype(np.uint8), (5, 5), 0)
    145. ret2, th2 = cv2.threshold(blur2, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    146. blur3 = cv2.GaussianBlur((self.imgReconstruite[:, :, 2]*255).astype(np.uint8), (5, 5), 0)
    147. ret3, th3 = cv2.threshold(blur3, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    148. plt.subplot(1,3,1)
    149. plt.title("hue")
    150. plt.imshow(th1,cmap="gray")
    151. plt.subplot(1, 3, 2)
    152. plt.title("Saturation")
    153. plt.imshow(th2,cmap="gray")
    154. plt.subplot(1, 3, 3)
    155. plt.title(self.type)
    156. plt.imshow(th3,cmap="gray")
    157. plt.show()
    158. self.saveOpencv(th1,"Hue_Binariser.png")
    159. self.saveOpencv(th2, "Saturation_Binariser.png")
    160. self.saveOpencv(th3, "Intensity_Binariser.png")
    161. return th3
    162. def recontructionToRGB(self):
    163. #Calcul de intensity global de chaque pixel on en a besoin pour la reconstruction
    164. self.GlobalIntensity()
    165. [l,c,d]=self.img_hsi.shape
    166. self.imgReconstruite=np.zeros([l,c,3])
    167. for i in range(l):
    168. for j in range(c):
    169. cx = self.getCX(self.img_hsi[i, j, 0],self.img_hsi[i, j, 1])
    170. cy = self.getCY(self.img_hsi[i, j, 0],self.img_hsi[i, j, 1])
    171. self.imgReconstruite[i, j, 0] = self.intensity[i, j] * (cx + 1)
    172. self.imgReconstruite[i, j, 1] = 0.5 * self.intensity[i, j] * (2 - cx - np.sqrt(3) * cy)
    173. self.imgReconstruite[i, j, 2] = 0.5 * self.intensity[i, j] * (2 - cx - np.sqrt(3) * cy)
    174. plt.subplot(2,3,1)
    175. plt.title("hue")
    176. plt.imshow(self.imgReconstruite[:, :, 0], cmap="gray")
    177. # self.saveOpencv((self.imgReconstruite[:, :, 0]*255).astype(np.uint8), "HSI_Teinte.png")
    178. cv2.imwrite("HSI_Teinte.png",(self.imgReconstruite[:,:,0]*255).astype(np.uint8))
    179. #plt.colorbar(ticks=[0, 60, 120, 179], orientation='horizontal', cmap=cm.hsv)
    180. plt.subplot(2,3,2)
    181. plt.title("saturation")
    182. plt.imshow(self.imgReconstruite[:, :, 1], cmap="gray")
    183. # save HSI
    184. cv2.imwrite("HSI_Saturation.png",(self.imgReconstruite[:,:,1]*255).astype(np.uint8))
    185. # self.saveOpencv((self.imgReconstruite[:, :, 1]*255).astype(np.uint8), "HSI_Saturation.png")
    186. plt.subplot(2,3,3)
    187. plt.title(self.type)
    188. plt.imshow(self.imgReconstruite[:, :, 2], cmap="gray")
    189. # save HSD
    190. #self.savePillow((self.imgReconstruite[:,:,2]*255).astype(np.uint8), self.path+"HSD_Density.tif")
    191. cv2.imwrite("HSI_Density.png",(self.imgReconstruite[:,:,2]*255).astype(np.uint8))
    192. # self.saveOpencv((self.imgReconstruite[:,:,2]*255).astype(np.uint8), "HSI_Density.png")
    193. #img = Image.fromarray(self.img_0[:, :, 0])
    194. #img.save(self.path + "img.tif")
    195. plt.show()
    196. def saveOpencv(self,img,path):
    197. cv2.imwrite(self.path2+"/"+self.path[19:]+"/"+path, img)
    198. def savePillow(self,img,path):
    199. img_to_save = Image.fromarray(img)
    200. img_to_save.save(self.path2+"/"+self.path[19:]+"/"+path)
    201. if __name__ == "__main__":
    202. # reading the image
    203. # path="Tumor_CD31_LoRes.png"
    204. #path = "../DataSet/BreastCancerCell_dataset/ytma10_010704_benign1_ccd.tif"
    205. # path = "../DataSet_Lomenie/tumor.png"
    206. path="t2.png"
    207. img = plt.imread(path)
    208. img = img[:, :, 0:3].astype(np.double)
    209. od=od.rgb_2_od2(img)
    210. # hsi
    211. h = HSD(od, path, "Density")
    212. h.chromaticite()
    213. h.calcule_HSI()
    214. # h.Kmeans2(3)
    215. # h.plotHS(Z)
    216. h.recontructionToRGB()
    217. h.binarisation()

     ColorDeconvolution.py

    1. #!/usr/bin/python
    2. import matplotlib.pyplot as plt
    3. import numpy as np
    4. import cv2
    5. from PIL import Image
    6. import hsd
    7. from matplotlib import cm
    8. class ColorDeconvolution:
    9. def __init__(self, img,path):
    10. self.img_0 = img
    11. self.stains = None
    12. self.od = None
    13. self.path=path
    14. self.path2 = "autodl-tmp/color-deconv/Color-Deconvolution/Resultat"
    15. def getStains(self):
    16. return self.stains
    17. def setImage(self, img):
    18. if img != None:
    19. self.img_0 = img
    20. def RGB_2_OD(self):
    21. [l, c, d] = self.img_0.shape
    22. self.od = np.zeros([l, c, d])
    23. for i in range(l):
    24. for j in range(c):
    25. for k in range(d):
    26. if self.img_0[i,j,k] != 0 :
    27. self.od[i,j,k] = np.log(self.img_0[i, j, k])
    28. return self.od
    29. def norm(self, vector):
    30. n = 0
    31. for i in vector:
    32. n = n + np.square(i)
    33. return np.sqrt(n)
    34. def separateStain(self):
    35. # He = [0.65; 0.70; 0.29];
    36. # Eo = [0.07; 0.99; 0.11];
    37. # DAB = [0.27; 0.57; 0.78];
    38. He = np.array([0.6500286, 0.704031, 0.2860126]) # Hematoxylin
    39. Eo = np.array([0.07, 0.99, 0.11]) # Eosine
    40. DAB = np.array([0.26814753, 0.57031375, 0.77642715]) # DAB
    41. Res = np.array([0.7110272, 0.42318153, 0.5615672]) # residual
    42. HDABtoRGB = np.array([He / self.norm(He), Eo / self.norm(Eo), DAB / self.norm(DAB)])
    43. RGBtoHDAB = np.linalg.inv(HDABtoRGB)
    44. [l,c,d] = self.img_0.shape
    45. self.stains = np.zeros([l, c, d])
    46. for i in range(l):
    47. for j in range(c):
    48. a = np.dot(self.od[i, j], RGBtoHDAB)
    49. b=self.od[i,j]
    50. self.stains[i, j, 0] = a[0]
    51. self.stains[i, j, 1] = a[1]
    52. self.stains[i, j, 2] = a[2]
    53. return self.stains
    54. # this methods shows the differents stains of each canal for the method of colorDeconvolution
    55. def binarisation(self,image):
    56. # Otsu's thresholding after Gaussian filtering
    57. blur = cv2.GaussianBlur(image, (5, 5), 0)
    58. ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    59. return th3
    60. def showStains(self):
    61. plt.subplot(1, 4, 1)
    62. plt.title("original")
    63. plt.imshow(self.img_0)
    64. plt.subplot(1, 4, 2)
    65. plt.title('Hematoxylin')
    66. plt.imshow(self.stains[:, :, 0], cmap="gray")
    67. cv2.imwrite("colorDeconvolution_Hematoxylin.png",(self.stains[:, :, 0]*255).astype(np.uint8))
    68. # self.saveOpencv((self.stains[:, :, 0]*255).astype(np.uint8), " colorDeconvolution_Hematoxylin.tif")
    69. plt.subplot(1, 4, 3)
    70. plt.title('Eosine')
    71. plt.imshow(self.stains[:, :, 1], cmap="gray")
    72. cv2.imwrite("colorDeconvolution_Eosin.png",(self.stains[:, :, 1]*255).astype(np.uint8))
    73. # self.saveOpencv((self.stains[:, :, 1]*255).astype(np.uint8), " colorDeconvolution_Eosin.tif")
    74. plt.subplot(1, 4, 4)
    75. plt.title('DAB')
    76. plt.imshow(self.stains[:, :, 2], cmap="gray")
    77. cv2.imwrite(" colorDeconvolution_DAB.png",(self.stains[:, :, 2]*255).astype(np.uint8))
    78. # self.saveOpencv((self.stains[:, :, 2]*255).astype(np.uint8)," colorDeconvolution_DAB.tif")
    79. plt.show()
    80. self.binarisation()
    81. def saveOpencv(self,img,path):
    82. cv2.imwrite(self.path2 + "/" + self.path[19:] + "/" + path, img)
    83. def savePillow(self,img,path):
    84. img_to_save = Image.fromarray(img)
    85. img_to_save.save(path)
    86. def binarisation(self):
    87. # Otsu's thresholding after Gaussian filtering
    88. # il faut lui donner image self.imgReconstruite et binariser chaque canal hue saturation et density
    89. blur1 = cv2.GaussianBlur((self.stains[:, :, 0] * 255).astype(np.uint8), (5, 5), 0)
    90. ret1, th1 = cv2.threshold(blur1, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    91. blur2 = cv2.GaussianBlur((self.stains[:, :, 1] * 255).astype(np.uint8), (5, 5), 0)
    92. ret2, th2 = cv2.threshold(blur2, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    93. blur3 = cv2.GaussianBlur((self.stains[:, :, 2] * 255).astype(np.uint8), (5, 5), 0)
    94. ret3, th3 = cv2.threshold(blur3, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    95. plt.subplot(1, 3, 1)
    96. plt.title("Hematoxyline")
    97. plt.imshow(th1,cmap="gray")
    98. plt.subplot(1, 3, 2)
    99. plt.title("Eosine")
    100. plt.imshow(th2,cmap="gray")
    101. plt.subplot(1, 3, 3)
    102. plt.title("DAB")
    103. plt.imshow(th3,cmap="gray")
    104. plt.show()
    105. self.saveOpencv(th1, "Hematoxyline.png")
    106. self.saveOpencv(th2, "Eosine.png")
    107. self.saveOpencv(th3, "DAB.png")
    108. return th3
    109. def saveOpencv(self,img, path):
    110. cv2.imwrite(self.path2 + "/" + self.path[19:] + "/" + path, img)
    111. if __name__=="__main__":
    112. # path="../DataSet_Lomenie/arn2.tif"
    113. path = "t2.png"
    114. img = plt.imread(path)
    115. img = img[:, :, 0:3]
    116. # deconvolution de couleur
    117. satin = ColorDeconvolution(img,path)
    118. satin.RGB_2_OD()
    119. image=satin.separateStain()
    120. satin.showStains()

     一共会有6张图的结果,最后一张是原图:

     

     

     最后

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

    源代码: 

    病理图像颜色反卷积,颜色反卷积可应用于H-DAB图像和组织病理学中常用的染色剂(H-E,HAEC,FASTRed,DAB)-深度学习文档类资源-CSDN下载介绍颜色反卷积算法的设计针对RGB摄像机获取的颜色信息,基于免疫组化技术使用的染色剂RGB分量光的更多下载资源、学习资料请访问CSDN下载频道.https://download.csdn.net/download/weixin_42934729/86438199

  • 相关阅读:
    [MT8766][Android12] 使用谷歌LPA实现ESIM功能的流程
    php解析html类库simple_html_dom(3)
    2.4GHz、DA14530-00000FX2射频收发器/LSM6DSOTR 6 轴运动传感器/SKY66423-11射频前端 860MHz 至 930MHz
    node-sass安装不上的问题
    pytest测试框架使用基础07 fixture—parametrize获取参数的几种常用形式
    一个Python中优雅的数据分块方法
    如何搭建测试环境?(详解)
    算法训练营第十三天 | 239. 滑动窗口最大值、347.前 K 个高频元素
    银行互联网类业务基于分布式存储的架构设计与实施运维分享
    Plant Simulation 与Web交互 V3.0工具
  • 原文地址:https://blog.csdn.net/weixin_42934729/article/details/126486336