2. import os
    3. import shutil
    4. import lxml.etree as ET
    6. from PIL import ImageDraw, Image
    7. import numpy as np
    8. import copy
    9. import cv2
    11. current_path = os.path.dirname(__file__)
    13. def dist(a, b):
    14.     return round(abs(a[0]-b[0]) + abs(a[1]-b[1]))
    16. def del_contours(contours, delete_list):
    17.     delta = 0
    18.     for i in range(len(delete_list)):
    19.         del contours[delete_list[i] - delta]
    20.         delta = delta + 1
    21.     return contours
    23. def xml_to_region(xml_file):
    24.     """
    25.     parse XML label file and get the points
    26.     :param xml_file: xml file
    27.     :return: region list,region_class
    28.     """
    30.     tree = ET.parse(xml_file)
    31.     region_list = []
    32.     region_class = []
    33.     for color in tree.findall('.//Annotation'):
    34.         if color.attrib['LineColor'] in ['65280','255','65535']:
    35.             # '65280'是绿色,'255'是红色,可以根据自己的实际情况更改这个判断条件(或者直接if True)
    36.             for region in color.findall('./Regions/Region'):
    37.                 vertex_list = []
    38.                 #region.attrib.get('Type')=='0':
    39.                 region_class.append(region.attrib.get('Type'))
    40.                 for vertex in region.findall('.//Vertices/Vertex'):
    41.                     # parse the 'X' and 'Y' for the vertex
    42.                     vertex_list.append(vertex.attrib)
    43.                 region_list.append(vertex_list)
    45.     return region_list,region_class
    48. def region_handler(im, region_list,region_class, level_downsample):
    49.     """
    50.     handle region label point to discrete point, and draw the region point to line
    51.     :param im: the image painted in region line
    52.     :param region_list: region list, region point,
    53.                     eg : [[{'X': '27381.168113', 'Y': '37358.653791'}], [{'X': '27381.168113', 'Y': '37358.653791'}]]
    54.     :param region_class : list,keep the value of region.attrib.get('Type') in elements of region list
    55.                     eg : [0,0,0,1,2,3]
    56.     :param level_downsample: slide level down sample
    57.     :return: image painted in region line of numpy array format
    58.     """
    60.     dr = ImageDraw.Draw(im)
    61.     for r_class, region in enumerate(region_list):
    62.         point_list = []
    63.         if region_class[r_class] == '0' or region_class[r_class] == '3':
    64.             for __, point in enumerate(region):
    65.                 X, Y = int(float(point['X'])/level_downsample), int(float(point['Y'])/level_downsample)
    66.                 point_list.append((X, Y))
    69.             #        points_length = len(point_list)
    70.             x_max = max(point_list, key=lambda point: point[0])[0]
    71.             x_min = min(point_list, key=lambda point: point[0])[0]
    72.             y_max = max(point_list, key=lambda point: point[1])[1]
    73.             y_min = min(point_list, key=lambda point: point[1])[1]
    74.             # mislabeled, here checked by x and y coordinate max and min difference
    75.             if (x_max - x_min < 50) or (y_max - y_min < 50): continue
    76.             if region_class[r_class] == '3':
    77.                 dr.arc(point_list, 0, 360, fill='#000000', width=12)
    78.             else:
    79.                 dr.line(point_list, fill="#000000", width=12)
    81.     return im
    83. def region_binary_image(tile, region_list,region_class, level_downsample):
    84.     """
    85.     convert the region labeled or not by doctor to binary image
    86.     :param tile: a return image based on the method of Slide class object in 'utils.openslide_utils'
    87.     :param region_list: region list, region point,
    88.                     eg : [[{'X': '27381.168113', 'Y': '37358.653791'}], [{'X': '27381.168113', 'Y': '37358.653791'}]]
    89.     :param region_class : list,keep the value of region.attrib.get('Type') in elements of region list
    90.                     eg : [0,0,0,1,2,3]
    91.     :param level_downsample: slide level down sample
    92.     :return: image painted in region line of numpy array format
    93.     """
    94.     im = Image.new(mode="1", size=tile.size)
    95.     dr = ImageDraw.Draw(im)
    96.     regions_list = []
    97.     for r_class, region in enumerate(region_list):
    98.         point_list = []
    99.         if region_class[r_class] == '0':
    100.             for __, point in enumerate(region):
    101.                 X, Y = int(float(point['X'])/level_downsample), int(float(point['Y'])/level_downsample)
    102.                 point_list.append((X, Y))
    104.             regions_list.append(point_list)
    106.     #由于医生的标注会出现不连续(非闭合)的情况，导致提取出来的标注坐标列表，会分成多段，比如：
    107.     #            0 (1979, 798) (2144, 1479)
    108.     #            1 (2139, 1483) (2319, 2162)
    109.     #            2 (2308, 2160) (3003, 1646)
    110.     # 正常情况下，标注坐标列表应该收尾闭合(前后坐标一致)，上下列表之间差异应该较大，如：
    111.     #            12 (1177, 2986) (1177, 2986)
    112.     #            13 (1507, 2942) (1507, 2940)
    113.     # 针对上述第一种情况，需要对提出出来的标注坐标列表进行循环判断，对收尾不一而且和其他标注列表的首坐标相对较近的话，进行合并处理
    115.     pin_jie_flag = [] #存储已经被拼接过的标注坐标列表序号
    116.     single_list = [] #存储新标注坐标列表的列表
    117.     for j,p_list in enumerate(regions_list):
    118.         if dist(p_list[0], p_list[-1]) < 50 and j not in pin_jie_flag:
    119.             #如果收尾坐标距离相差在150范围内(曼哈顿距离)，且未成被拼接过，直接认为这个组坐标无须拼接，存储起来
    120.             single_list.append(p_list)
    121.         elif dist(p_list[0], p_list[-1]) > 50 and j not in pin_jie_flag:
    122.             #如果收尾坐标距离相差在150范围外(曼哈顿距离)，且未成被拼接过，说明这组坐标是残缺非闭合的，需要对其余标注坐标进行新一轮的循环判断
    123.             for j_2,p_list_2 in enumerate(regions_list):
    124.                 if j_2 != j and j_2 not in pin_jie_flag:
    126.                     if dist(p_list[-1],p_list_2[0]) < 50 :
    127.                         p_list = p_list + p_list_2.copy()
    128.                         pin_jie_flag.append(j_2)
    129.                     elif dist(p_list[0],p_list_2[-1]) < 50 :
    130.                         p_list = p_list_2.copy() + p_list
    131.                         pin_jie_flag.append(j_2)
    132.                     elif dist(p_list[-1],p_list_2[-1]) < 50 :
    133.                         p_list_2_new = copy.deepcopy(p_list_2)
    134.                         p_list_2_new.reverse()
    135.                         p_list = p_list + p_list_2_new
    136.                         pin_jie_flag.append(j_2)
    137.                     elif dist(p_list[0],p_list_2[0]) < 50 :
    138.                         p_list_2_new = copy.deepcopy(p_list_2)
    139.                         p_list_2_new.reverse()
    140.                         p_list = p_list_2_new + p_list
    141.                         pin_jie_flag.append(j_2)
    142.                     # 当这组非闭合的尾坐标和其他组坐标的首坐标接近到一定范围时(距离是150内),就让当前的非闭合的坐标列表和该组坐标列表相加
    143.                     # 处理完毕之后，将该组坐标的序号增加到已拼接坐标的列表中，确保后续循环不会再判断这个列表
    144.             single_list.append(p_list)
    147.     for points in single_list:
    148.         dr.polygon(points, fill="#ffffff")
    150.     #由于医生的标注除了出现不连续(非闭合)的情况外，还存在多余勾画的情况，对这种情况暂时没有完整的思路予以接近，先用
    151.     # opencv中的开闭操作组合来进行修补
    152.     kernel = np.ones((20,20),np.uint8)
    153.     filter_matrix = np.array(im).astype(np.uint8)
    154.     filter_matrix = cv2.morphologyEx(filter_matrix, cv2.MORPH_OPEN, kernel)
    155.     filter_matrix = cv2.morphologyEx(filter_matrix, cv2.MORPH_CLOSE, kernel)
    156.     #    plt.imshow(filter_matrix)
    157.     return filter_matrix
    160. class Region:
    161.     """"
    162.     handle the template xml format file to insert label svs region
    163.     """
    164.     def __init__(self, xml_file):
    165.         parser = ET.XMLParser(remove_blank_text=True)
    166.         if not os.path.isfile(xml_file):
    167.             template = os.path.join(current_path, "template.xml")
    168.             shutil.copy(template, xml_file)
    169.         self._xml_file = xml_file
    170.         self._tree = ET.parse(xml_file, parser)
    172.     def get_region(self, region_id):
    173.         """
    174.         get region by region id
    175.         :param region_id: region id, 0: green, 1: yellow, 2: red, see the template.xml
    176.         :return: the region
    177.         """
    179.         return self._tree.findall(".//Annotation/Regions")[region_id]
    181.     def add(self, region_id, points):
    182.         """
    183.         add one region to the specified region by region id, the added region is ellipse
    184.         and the parameter points is a rectangle bounded by an ellipse
    185.         :param points: list with two element(upper-left, bottom-right), is the rectangle bounded by an ellipse
    186.         :return:
    187.         """
    189.         region = self.get_region(region_id)
    190.         region_num = len(region.findall(".//Region"))
    191.         region_attr = {
    192.             "Id": str(region_num+1),
    193.             "Type": "2",
    194.             "Zoom": "1",
    195.             "Selected": "1",
    196.             "ImageLocation": "",
    197.             "ImageFocus": "0",
    198.             "Length": "80",
    199.             "Area": "400",
    200.             "LengthMicrons": "20",
    201.             "AreaMicrons": "30",
    202.             "Text": "",
    203.             "NegativeROA": "0",
    204.             "InputRegionId": "0",
    205.             "Analyze": "0",
    206.             "DisplayId": "1"
    207.         }
    209.         region_tag = ET.Element("Region", region_attr)
    210.         region.append(region_tag)
    212.         attributes = ET.SubElement(region_tag, "Attributes")
    213.         vertices = ET.Element("Vertices")
    214.         region_tag.append(vertices)
    216.         for point in points:
    217.             # insert point
    218.             ET.SubElement(vertices, "Vertex", attrib=point)
    220.     def save(self):
    221.         """
    222.         save the xml file
    223.         :return:
    224.         """
    226.         self._tree.write(self._xml_file, pretty_print=True)
    229. def color_int_to_str(color_int):
    230.     color_str = hex(color_int)[2:]
    231.     assert len(color_str) <= 6, 'Found unknow color!'
    232.     pad_count = 6 - len(color_str)
    233.     color_str = ''.join(['0'] * pad_count) + color_str
    234.     b, g, r = color_str[0:2], color_str[2:4], color_str[4:6]
    235.     return r+g+b
    238. def color_str_to_int(color_str):
    239.     assert len(color_str) == 6, 'Found unknow color!'
    240.     r, g, b = color_str[0:2], color_str[2:4], color_str[4:6]
    241.     color_int = (int(r, 16)) + (int(g, 16) << 8) + (int(b, 16) << 16)
    242.     return color_int
    245. def contours_to_xml(savepath,contours,mode,if_add = True,level_downsample = 16,mpp= "0.252100",linecolor ="16711680",contour_area_threshold=0):
    246.     """
    247.     based on a mask of svs file(mask sure the size of the mask equals the size of the svs file 's level_dimensions in level 2) to make a
    248.     xml format lable file for this svs file
    249.     :param savepath :  the xml format lable file save file path
    250.     :param contours :  contours list return from cv2.findContours of the mask
    251.     :param if_add : Added niew Annotation to an exits xml format label file or not ,defaut False
    252.     :param level_downsample : the value of slide.level_downsamples[2]
    253.     :param mpp : the value of MicronsPerPixel in slide.properties['openslide.mpp-x']
    254.     :param linecolor : the value of decimal color code to draw contours in xml format lable file ,default color is blue
    255.     :param contour_area_threshold : the threshold to drop small contours base on cv2.contourArea,which helps to keep the big area contours in xml format lable  file
    256.     :return:
    257.     """
    259.     ann_begin_tag = 1
    260.     Annotations = ET.Element('Annotations', {'MicronsPerPixel': mpp})
    261.     origin_color_list = []
    262.     if if_add and os.path.exists(savepath):
    263.         origin = ET.parse(savepath)
    264.         ann_begin_tag = len(origin.findall('.//Annotation')) + 1
    265.         for ann in origin.findall('.//Annotation'):
    266.             origin_color_list.append(ann.attrib['LineColor'])
    267.             Annotations.append(ann)
    268.     if mode == 'fp':
    269.         linecolor = '65535'
    270.     elif mode == 'fn':
    271.         linecolor = "65280"
    272.     elif mode == 'all':
    273.         linecolor = '65280'
    274.     # if linecolor in origin_color_list: linecolor = "13382297"
    275.     Annotation = ET.SubElement(Annotations, 'Annotation',
    276.                                {'Id': str(ann_begin_tag), 'Name': '', 'ReadOnly': '0', 'NameReadOnly': '0',
    277.                                 'LineColorReadOnly': '0', 'Incremental': '0', 'Type': '4',
    278.                                 'LineColor': linecolor, 'Visible': '1', 'Selected': '1',
    279.                                 'MarkupImagePath': '', 'MacroName': ''})
    280.     Attributes = ET.SubElement(Annotation, 'Attributes')
    281.     ET.SubElement(Attributes, 'Attribute', {'Name': '', 'Id': '0', 'Value': ''})
    282.     Regions = ET.SubElement(Annotation, 'Regions')
    283.     RegionAttributeHeaders = ET.SubElement(Regions, 'RegionAttributeHeaders')
    284.     ET.SubElement(RegionAttributeHeaders, 'AttributeHeader',
    285.                   {'Id': "9999", 'Name': 'Region', 'ColumnWidth': '-1'})
    286.     ET.SubElement(RegionAttributeHeaders, 'AttributeHeader',
    287.                   {'Id': "9997", 'Name': 'Length', 'ColumnWidth': '-1'})
    288.     ET.SubElement(RegionAttributeHeaders, 'AttributeHeader',
    289.                   {'Id': "9996", 'Name': 'Area', 'ColumnWidth': '-1'})
    290.     ET.SubElement(RegionAttributeHeaders, 'AttributeHeader',
    291.                   {'Id': "9998", 'Name': 'Text', 'ColumnWidth': '-1'})
    292.     ET.SubElement(RegionAttributeHeaders, 'AttributeHeader',
    293.                   {'Id': "1", 'Name': 'Description', 'ColumnWidth': '-1'})
    294.     i = 1
    295.     contour_area_threshold = 10000 / (level_downsample ** 2)
    296.     # contour_area_threshold = 0
    297.     for cnt in contours:
    298.         contour_area = cv2.contourArea(cnt)
    299.         if contour_area > contour_area_threshold:
    300.             Region = ET.SubElement(Regions, 'Region',
    301.                                    {'Id': str(i), 'Type': '0', 'Zoom': '0.011', 'Selected': '0',
    302.                                     'ImageLocation': '', 'ImageFocus': '-1', 'Length': str(cnt.shape[0]), 'Area': str(level_downsample**2*contour_area),
    303.                                     'LengthMicrons': '0', 'AreaMicrons': '0', 'Text': '', 'NegativeROA': '0',
    304.                                     'InputRegionId': '0', 'Analyze': '1', 'DisplayId': str(i)})
    305.             ET.SubElement(Region, 'Attributes')
    306.             Vertices = ET.SubElement(Region, 'Vertices')
    307.             cnt = np.squeeze(np.asarray(cnt))
    308.             for j in range(cnt.shape[0]):
    309.                 ET.SubElement(Vertices, 'Vertex', {'X': str(int(cnt[j,0]*level_downsample)), 'Y': str(int(cnt[j,1]*level_downsample))})
    310.             i = i + 1
    311.     ET.SubElement(Annotation, 'Plots')
    312.     doc = ET.ElementTree(Annotations)
    313.     doc.write(open(savepath, "wb"), pretty_print=True)
    315. if __name__ == '__main__':
    316.     import sys,glob,os
    317.     sys.path.append('../')
    318.     from utils.openslide_utils import Slide
    319.     import matplotlib.pyplot as plt
    320.     plt.rcParams['figure.figsize'] = 15, 15
    321.     # all_svs_path = sorted(glob.glob('/media/totem_disk/totem/Data_20201014/LABIAL_SALIVARY_GLAND_BIOPSY/*.svs'))
    322.     all_svs_path = sorted(glob.glob('/media/totem_disk/totem/data_20200903/second_batch/*.svs'))
    323.     for svs_path in all_svs_path:
    324.         slide = Slide(svs_path)
    325.         img_name = os.path.splitext(svs_path)[0].split('/')[-1]
    326.         # draw_mask = cv2.imread('/media/totem_disk/totem/guozunhu/Project/lp_3_2_u2net/predict/12.18_efficient_b0_s1:d0:bce+lo_so + d6:bce_s2:d0:bce+lo_so/new_dataset_batch2/hotpic/{}_4.png'.format(img_name),0)
    327.         draw_mask = cv2.imread('/media/totem_disk/totem/guozunhu/Project/lp_4_hooknet/xml_out/new_dataset_batch2/{}_hm_m.png'.format(img_name),0)
    328.         draw_mask = np.where(draw_mask<128,0,draw_mask)
    329.         draw_mask = np.where(draw_mask>=128,255,draw_mask)
    331.         # kernel = cv2.getStructuringElement(cv2.MORPH_RECT,(3,3))
    332.         # draw_mask = cv2.morphologyEx(draw_mask,cv2.MORPH_CLOSE,kernel)
    334.         region_process_mask = draw_mask
    335.         if cv2.__version__[0]=='4':
    336.             cnts, _ = cv2.findContours(region_process_mask,mode=cv2.RETR_EXTERNAL,method=cv2.CHAIN_APPROX_SIMPLE)
    337.         else:
    338.             _, cnts , _ = cv2.findContours(region_process_mask,mode=cv2.RETR_EXTERNAL,method=cv2.CHAIN_APPROX_SIMPLE)
    339.         # delete_list = []
    340.         # delete_show_list = []
    341.         # for i in range(len(cnts)):
    342.         #     if cv2.contourArea(cnts[i]) < 1000:
    343.         #     # if cv2.contourArea(cnts[i]) < draw_mask.shape[0]*draw_mask.shape[1]*0.01:
    344.         #         delete_list.append(i)
    345.         #         delete_show_list.append(cnts[i])
    346.         # cnts = del_contours(cnts,delete_list)
    347.         # draw_mask_bgr = draw_mask.copy()
    348.         # draw_mask_bgr = cv2.cvtColor(draw_mask_bgr,cv2.COLOR_GRAY2BGR)
    349.         # cv2.drawContours(draw_mask_bgr,delete_show_list,-1,(0,0,0),-1)
    350.         # save_filter_mask = cv2.cvtColor(draw_mask_bgr,cv2.COLOR_BGR2GRAY)
    351.         # cv2.imwrite('/media/totem_disk/totem/guozunhu/Project/lp_2/cache_ndpi_he_cc/hot_pic/c_xml/filter/{}.png'.format(img_name),save_filter_mask)
    352.         savepath = '/media/totem_disk/totem/guozunhu/Project/lp_4_hooknet/xml_out/new_dataset_batch2/xml/{}.xml'.format(img_name)
    353.         mpp = str(slide.get_mpp()*1000)
    354.         # level_dimension = slide.get_level_dimension(level=0)
    355.         level_downsample = slide.get_level_downsample(level=1)
    356.         level_downsample = 8
    357.         mode = 'all'
    358.         # need_cnts_list = []
    359.         # for i in range(len(cnts)):
    360.         #     if cv2.contourArea(cnts[i]) > 10000 / (level_downsample ** 2):
    361.         #         need_cnts_list.append(cnts[i])
    362.         # if img_name in ['S04','S10','S11','S13']:
    363.         #     level_downsample = level_downsample // 2
    364.         contours_to_xml(savepath,cnts,mode,if_add = False,level_downsample=level_downsample)
    365.         print(img_name)