使用多模板匹配进行DICOM图像的自动勾画

目的

• 放射治疗是治疗恶性肿瘤的主要手段之一。其基本原则是最大限度地保证射线剂量全部集中在放疗靶区内（包括GTV、CTV和PTV）,并且使危及器官（Organs at Risk, OARs）少受辐射剂量。
• 目前医生需要对图像的切片进行逐层勾画。过程耗时费力，影响了医生的工作效率，并且勾画极度依赖医生的知识和个人经验。
• 因而需一种辅助医生勾画的自动勾画系统

什么是模板匹配

• 定义
  模板匹配是一种查找与模板图像(patch)匹配(相似)的图像区域的技术。
• 工作过程

1. 源图像（I）：我们期望找到与模板图像匹配的图像
2. 模板图像（T）：将与源图像进行比较的图像块
3. 要识别匹配区域，我们必须通过滑动来比较模板图像与源图像。通过滑动，我们将图像块一次移动一个像素（从左到右，从上到下）。在每个位置，计算度量，以便它表示在该位置处的匹配的“好”还是“坏”（或者与图像的特定区域相似）。
   
4. 对于 T 覆盖在 I 上的每个位置,你把度量值 保存 到 结果图像矩阵 ® 中. 在 R 中的每个位置 (x,y) 都包含匹配度量值
   
5. 上图就是 TM_CCORR_NORMED 方法处理后的结果图像 R . 最白的位置代表最高的匹配. 正如您所见, 红色椭圆框住的位置很可能是结果图像矩阵中的最大数值, 所以这个区域 (以这个点为顶点,长宽和模板图像一样大小的矩阵) 被认为是匹配的.
6. 实际上, 我们使用函数 minMaxLoc 来定位在矩阵 R 中的最大值点 (或者最小值, 根据函数输入的匹配参数) .

• 匹配算法：

8. 标准平方差匹配 method=CV_TM_SQDIFF_NORMED
9. 标准平方差匹配 method=CV_TM_SQDIFF_NORMED
10. 标准相关匹配 method=CV_TM_CCORR_NORMED

什么是自动勾画

• 目前临床上主要应用的是基于Atlas的自动勾画，常用软件有瓦里安的Velocity、MIMvista、基于图谱库ABAS自动勾画软件等
• 以卷积神经网络为代表的深度学习方法已经广泛应于自动分割。
• 本文探讨的是传统方法

如何进行基于多模板匹配的自动勾画

• 将已勾画好的DICOM图像裁剪后作为模板图像集中起来，用于多模板匹配。
• 新（源）图像进来时，通过与模板库中图像比较，找出最匹配的候选者。
• 将候选者对应的mask label复制给源图像即可。
• 模板匹配得到的自动勾画可能与真实勾画差异很大，其效果并不如atlas方法和深度学习方法。

template_path = './template_path/'#存放模板的位置
label_path = './label_path/'# 存放勾画的位置
image_path = './image_path/'#欲匹配的图像位置

//读取templates图像
def read_templates(path):
    file_list = os.listdir(path)
    listTemplates = []
    for file in file_list:
        with gzip.open(template_path + file, 'rb') as unzipped_data:
            info_item = []
            info_item.append(file)
            template_data = pickle.load(unzipped_data)
            info_item.append(template_data)
            listTemplates.append(info_item)
    return listTemplates

//给出匹配的模板名字和匹配得分
def findMatches(listTemplates,target,score_threshold):
    listHit = []
    for templateName, templateValue in listTemplates:
        if templateValue.dtype == "float64": 
            raise ValueError("64-bit not supported, max 32-bit")
        else:
            corrMap = cv2.matchTemplate(target,templateValue, method=cv2.TM_CCOEFF_NORMED)
            minVal, maxVal, minLoc, maxLoc = cv2.minMaxLoc(corrMap)
            coeff = maxVal
            newHit = {'TemplateName':templateName,'Score':coeff}
            listHit.append(newHit)
    return pd.DataFrame(listHit)

//使用pydicom读取dicom图片
dcms = glob.glob(os.path.join(image_path,'*.dcm'))
slices = [pydicom.read_file(dcm) for dcm in dcms]
slices.sort(key = lambda x : float(x.ImagePositionPatient[2]))
images = np.stack([s.pixel_array for s in slices],axis = -1).astype(np.float32)
targets = images * slices[0].RescaleSlope + slices[0].RescaleIntercept
Hits = findMatches(listTemplates,targets[:,:,0],score_threshold=0.6)
print(Hits)
#找出得分最高的那张slice名称
#HH = Hits.sort_values('Score',ascending=False)
Hits_group = Hits.groupby(by='Score',as_index=False).max()
print(Hits_group)
num = len(Hits_group.index)-1
print(num)
rightName = Hits_group.iloc[[num],[1]].values[0][0]#iloc是基于position进行索引的
print(rightName)