文章目录

    • 数据整理
      • 整理出所有的dcm文件
      • 整理出所有的xml标注文件
      • 整理数据①——舍弃错误的标注文件
      • 整理数据②——两个标注文件指向同一个目标图片的情况
        • 封装函数，传入xml文件，显示标注效果
      • 整理数据③——将PETCT的三通道图像转成平扫CT的单通道图像格式
        • 关于PETCT三通道转单通道灰度图的一些思考
      • 建立Dataset

本文进一步对数据集进行探索，对数据进行清洗、整理、验证。
最后建立Dataset对象，以便进行后续的训练工作。
本文所用代码： 我的Github

数据整理

我们最终的目的是建立Dataset对象，所以我们的目标应该是整理出 数据+标签 的组合。
在目标检测任务中，数据即我们的CT图片（读取自dcm文件）；标签即我们的xml文件，内含 Bounding Box框 + 类别标签。
所以我们的目标是对所有可配对的 dcm文件和 xml文件 进行一一对应，剩下读取则是一件简单的工作。

整理出所有的dcm文件

dcm文件和xml文件对应的方法在上一节中已经阐释。
现在需要将数据集中所有的 dcm文件 和 对应的SOP Instance UID（后通称UID）整理出来，一共有25万条记录。

import pydicomimport matplotlib.pyplot as pltimport osfrom tqdm import tqdmimport pandas as pdimport numpy as np# 查看整个数据集下所有的dcm文件名dcm_file=[]for root, dirs, files in os.walk('manifest-1608669183333/Lung-PET-CT-Dx'):for file in files:file_path = os.path.join(root, file)if 'dcm' in file_path:dcm_file.append(file_path)print(dcm_file[0])len(dcm_file)

导入DataFrame后，获取对应的UID。
由于有25w张图片，我在固态硬盘上跑这段代码耗时约50min，跑完后保存为csv，下次就不用再跑一遍了。

df_dcm = pd.DataFrame(dcm_file)# 建立新列uid，获取dcm文件对应的 SOPInstanceUID 号，这步耗时50min左右df_dcm['uid'] = df_dcm[0].apply(lambda x: pydicom.dcmread(x).SOPInstanceUID)df_dcm['uid_str'] = df_dcm['uid'].apply(lambda x: str(x))df_dcm.to_csv('dcm_file_uid.csv')# 上面一步耗时较久，如果已经读取过一遍并已保存，这里直接读取df_dcm = pd.read_csv('dcm_file_uid.csv', index_col=0)df_dcm

简单检查一下dcm文件是否独一无二。

整理出所有的xml标注文件

我们再把所有的xml文件整理出来，附上对应的UID。

整理数据①——舍弃错误的标注文件

我们先检查一下3万多个xml文件对应的UID是否都对应在df_dcm中。

# 检查Annotation目录下所有的uid文件名是否都在df_dcm中xml_file_df = xml_file_df.loc[xml_file_df['uid_str'].isin(df_dcm['uid_str'])]xml_file_df# 我们发现 30884 ＜ 31562


结果我们发现并不是所有的xml都有对应原图片，这应该属于这个数据集的标注错误问题，有可能标注文件的对应UID写错了。
总之，我们不能将这些xml文件纳入训练，所以将这些找不到对应图片的标注文件舍弃掉。

整理数据②——两个标注文件指向同一个目标图片的情况

我们查看一下XML标注文件是否有重复的文件名。

# 查看Annotation目录下的XML文件是否有重复的文件名len(xml_file_df[0].unique()) == len(xml_file_df['uid_str'].unique())# 如有False，说明Annotation下有重复的文件


这里我们发现该数据集第二个坑，Annotation目录下存在重复的标注文件指向同一个目标图片的情况（否则上面就显示True）。
我们继续进一步探究是哪些标注文件重复了。

我们查看一下原本的XML文件数量 和 去重后的XML文件数量。

# 显示Annotation目录下去重后的XML文件数量、原本的XML文件数量print(len(xml_file_df['uid_str'].unique()), len(xml_file_df[0].unique()) )

我们获取一下那些重名（同一UID）的xml文件。

# 获取那些重复的xml文件的索引idx = xml_file_df['uid_str'][xml_file_df['uid_str'].duplicated(keep=False)].index# 对xml_file_df文件应用idx索引，获取完整文件名xml_file_name_duplicated = pd.DataFrame(xml_file_df.loc[idx,:])# 导出到csv文件(记录一下)xml_file_name_duplicated.to_csv('duplicated.csv')xml_file_name_duplicated[:6]

xml_file_name_duplicated = pd.read_csv('duplicated.csv', index_col=0)print(len(xml_file_name_duplicated))# 显示重复条目数量xml_file_name_duplicated['0'][:12].values# 显示前12行# 可以看到A0003和A0004中存在相同名称的文件。# 我们知道A0003和A0004是两个不同的患者(Clinical Data的“statistics-clinical-20201221.xlsx”文件)，理应不可能存在相同的CT影像，所以应该是xml标注文件本身出错了


以上我们可以看到A0003和A0004中存在相同名称的文件。
我们知道A0003和A0004是两个不同的患者(从Clinical Data的“statistics-clinical-20201221.xlsx”文件看出)，理应不可能存在相同的CT影像，所以应该是xml标注文件本身出错了。
那么具体是怎么出错了呢？

# 查看重复条目有多少种，如果为重复条目总数的1/2，那么均为两两重复print('重复条目有', len(xml_file_name_duplicated['uid_str'].unique()), '种，总条目为', len(xml_file_name_duplicated), '个')# 下面我们看看这些重复的标注文件的标注情况


上面我们得知重复标注的情况是两两重复。
我们再进一步查看一下重复标注的情况。

封装函数，传入xml文件，显示标注效果

我们设计一个函数，传入xml文件路径，显示图片和锚框。

# 设计一个函数，传入xml文件路径，显示图片和锚框################### 从XML文件获取bounding box信息import xml.etree.ElementTree as ETdef get_labelFromXml(xml_file): label=[]bbox_list=[]an_file = open(xml_file, encoding='utf-8')tree=ET.parse(an_file)root = tree.getroot()for object in root.findall('object'):cancer_type=object.find('name').text.upper()xmin=object.find('bndbox').find('xmin').textxmax=object.find('bndbox').find('xmax').textymin=object.find('bndbox').find('ymin').textymax=object.find('bndbox').find('ymax').textif int(xmin)==0 or int(xmax)==0 or int(ymin)==0 or int(ymax)==0:passelif int(xmin)==int(xmax) or int(ymin)==int(ymax):passelse:bbox=[int(xmin),int(ymin),int(xmax),int(ymax)]bbox_list.append(bbox)label.append(cancer_type)return bbox_list,label# 画布上添加锚框def bbox_to_rect(bbox, color):# 将边界框(左上x,左上y,右下x,右下y)格式转换成matplotlib格式：# ((左上x,左上y),宽,高)return plt.Rectangle(xy=(bbox[0], bbox[1]), width=bbox[2]-bbox[0], height=bbox[3]-bbox[1],fill=False, edgecolor=color, linewidth=2)def show_img_and_anchors(url):# 传入xml文件路径，展示标注效果，一张图片可显示多个bbox框# 获取uid = xml_file_df.loc[xml_file_df[0]==url, 'uid_str'].values[0]im = df_dcm.loc[df_dcm['uid_str']==uid].values[0][0]im = pydicom.read_file(im)# 获取像素矩阵img_arr = im.pixel_arrayfig = plt.imshow(img_arr,cmap=plt.cm.bone)plt.title("UID:{}".format(uid))bbox, label = get_labelFromXml(url)for i, j in zip(bbox, label):fig.axes.add_patch(bbox_to_rect(i, 'red'))fig.axes.text(i[0]+12, i[1]+12, j,va='center', ha='center', fontsize=12, color='red')def show_imgs_and_anchors(urls, num=5):# 传入xml文件路径组成的列表，展示多张图片发标注效果num_toshow = min(num, len(urls))print('载入', len(urls), '张图片，显示前', num, '张。')# 获取temp_df = xml_file_df.loc[xml_file_df[0].isin(urls)]urls = temp_df[0].values# 这里urls顺序被重排，需要重新赋值以此，以和uids建立一一对应uids = temp_df['uid_str'].valuesims = [df_dcm.loc[df_dcm['uid_str']==x].values[0][0] for x in uids]ims_dcm = [pydicom.read_file(x).pixel_array for x in ims]# 获取像素矩阵List# 取标注框和类别bbox_label = [get_labelFromXml(i) for i in urls]# _表示忽略不适用的变量，返回的fig用不上_, axes = plt.subplots(nrows=1, ncols=num_toshow, figsize=(48,48))# 使用zip方法，在for循环中设置axes中的各个子区域的参数并绘图for ax, uid, img, lbl in zip(axes, uids, ims_dcm, bbox_label):ax.imshow(img, cmap=plt.cm.bone)ax.set_title(uid)bbox, label = lblfor i, j in zip(bbox, label):ax.add_patch(bbox_to_rect(i, 'red'))ax.text(i[0]+12, i[1]+12, j,va='center', ha='center', fontsize=12, color='red')ax.axes.get_xaxis().set_visible(False)ax.axes.get_yaxis().set_visible(False)

展示标注效果（一个xml文件）：

展示标注效果（多个xml文件）：

我们看到上方第1行第1张图片和第2行第1张图片是标注了同一张CT图片的情况。
我们可以看到标注的质量差不多，若标注得差不多，则可以假设不同的标注框可以作为一种数据增强的形式，因此可以无需特意舍弃重复的对应同UID的xml文件。
下面抽样查看一下所有的的重复标记情况。

我们可以再随机抽20张看看

# 随机抽20张看看slice = np.random.choice(len(xml_file_name_duplicated), size=(4,5), replace=False)# 生成4×5整数矩阵for i in slice:urls = xml_file_name_duplicated['0'].values[i]show_imgs_and_anchors(urls, num=5)

效果见源码文件。

整理数据③——将PETCT的三通道图像转成平扫CT的单通道图像格式

普通CT图像是单通道（512×512），但PETCT图像是三通道（512×512×3），需转成单通道进行训练。
因为日后我们希望使用平扫的CT片子进行推理预测。

我们看看读取PETCT和普通CT片子的图片shape：

# 各取一张ct，一张petct图像的xml文件url_ct = 'Lung-PET-CT-Dx-Annotations-XML-Files-rev12222020/Annotation\\A0213\\1.3.6.1.4.1.14519.5.2.1.6655.2359.211226049141098860991228194230.xml'# 普通CT图像url_petct = 'Lung-PET-CT-Dx-Annotations-XML-Files-rev12222020/Annotation\\A0213\\1.3.6.1.4.1.14519.5.2.1.6655.2359.284461122128650673200074494931.xml'# PETCT图像# 从xml文件获取对应ct图像文件信息def get_img_pixel_array(url):uid = xml_file_df.loc[xml_file_df[0]==url, 'uid_str'].values[0]im = df_dcm.loc[df_dcm['uid_str']==uid].values[0][0]im = pydicom.read_file(im)img_arr = im.pixel_array# 获取像素矩阵return img_arrprint('ct图像矩阵的形状为，', get_img_pixel_array(url_ct).shape)print('petct图像矩阵的形状为，', get_img_pixel_array(url_petct).shape)

关于PETCT三通道转单通道灰度图的一些思考

我这里使用了opencv库的方法将PETCT的三通道转成单通道灰度图。
在opencv库中，三通道转单通道主要是两种转换方法：COLOR_BGR2GRAY or COLOR_RGB2GRAY

我们知道，直接使用cv2读取后图片的通道顺序是BGR，而pydicom读取的pixel_array的通道顺序应该是RGB（调用PhotometricInterpretation属性可以看到）。
所以对pixel_array直接使用cv2的转灰度方法时，可以使用COLOR_RGB2GRAY模式。

但是我们知道PETCT的图像信息是 正常的CT影像 和 代谢摄取增高影（高亮部分）两者重合。
我们发现使用COLOR_BGR2GRAY模式转灰度后，代谢摄取增广的高亮信息较COLOR_RGB2GRAY模式不明显。（下面的代码块可见效果）
这样似乎更适合用于训练，因为后面我们需要常用普通平扫的CT片子进行预测。

因此，转单通道我采用了cv2的COLOR_BGR2GRAY模式。

关于三通道转灰度图的一些其他参考资料：

  1. https://note.nkmk.me/python-opencv-numpy-color-to-gray/
  2. https://stackoverflow.com/questions/62855718/why-would-cv2-color-rgb2gray-and-cv2-color-bgr2gray-give-different-results

两种转单通道模式的对比：

import cv2 as cvimg_array = get_img_pixel_array(url_petct)img_array_gray = cv.cvtColor(img_array, cv.COLOR_BGR2GRAY)img_array_gray2 = cv.cvtColor(img_array, cv.COLOR_RGB2GRAY)# # img_array = np.array(img_array, dtype=np.uint16)_, axes = plt.subplots(1,3,figsize=(24,24))axes[0].imshow(img_array, cmap=plt.cm.bone)axes[0].set_title('PETCT')axes[1].imshow(img_array_gray, cmap=plt.cm.bone)axes[1].set_title('COLOR_BGR2GRAY')axes[2].imshow(img_array_gray2, cmap=plt.cm.bone)axes[2].set_title('COLOR_RGB2GRAY')plt.show()plt.close()

建立Dataset

有了前面的数据整理，我们可以建立Dataset数据集对象了。
目标是传入 [dcm图片集] 和 对应的 [xml标注集] ，即生成一个Dataset对象。

import torchfrom torch.utils.data import Dataset, DataLoaderfrom torchvision import transformsfrom PIL import Imagetransform=transforms.Compose([transforms.ToTensor(),])

但是，这里我们要注意，pydicom从CT读取的pixel_array数据类型（dtype）为uint16，从PETCT读取的pixel_array则为uint8。
这会导致Image.open（PIL库）读取并使用transforms.ToTensor()转成Tensor后数据类型不一样。前者变为torch.int16，后者变为torch.float32。（如下演示）

因此，我们需要将pixel_array统一转成float32。

创建Dataset数据集对象（Pytorch）。
注意一下Dataset输出的数据格式，img为图片Tensor，label则是一个dict，内含bbox的坐标和类别。

class_to_id={'A':1, 'B':2, 'E':3, 'G':4}# 重写get_labelFromXml，返回label的iddef get_labelidFromXml(xml_file): label=[]bbox_list=[]an_file = open(xml_file, encoding='utf-8')tree=ET.parse(an_file)root = tree.getroot()for object in root.findall('object'):cancer_type=object.find('name').text.upper()xmin=object.find('bndbox').find('xmin').textxmax=object.find('bndbox').find('xmax').textymin=object.find('bndbox').find('ymin').textymax=object.find('bndbox').find('ymax').textif int(xmin)==0 or int(xmax)==0 or int(ymin)==0 or int(ymax)==0:passelif int(xmin)==int(xmax) or int(ymin)==int(ymax):passelse:bbox=[int(xmin),int(ymin),int(xmax),int(ymax)]bbox_list.append(bbox)label.append(class_to_id[cancer_type])return bbox_list,label# 计划传入img和label，以此对应dcm文件路径集和xml文件路径集class LungDetection(Dataset):def __init__(self, img, label, transform=transform,):self.img = imgself.label = labelself.transform = transformdef __getitem__(self, index):img = self.img[index]label = self.label[index]img_open=pydicom.read_file(img)img_array=img_open.pixel_arrayif len(img_array.shape) == 3:img_array = cv.cvtColor(img_array, cv.COLOR_BGR2GRAY)img_array = np.array(img_array, dtype=np.float32)img_pic = Image.fromarray(img_array)img_tensor=self.transform(img_pic)bbox, label = get_labelidFromXml(label)bbox_tensor = torch.as_tensor(bbox,dtype=torch.float32)label_tensor = torch.as_tensor(label,dtype=torch.int64)target={}target['boxes']=bbox_tensortarget['labels']=label_tensorreturn img_tensor, targetdef __len__(self):return len(self.img)

xml_file_df 是我们之前整理出的xml标签文件，我们现在将xml文件和dcm文件对应起来。

# xml_file_df 是我们之前整理出的xml标签文件，我们现在将xml文件和dcm文件对应起来xml_file_dataset = xml_file_df.reset_index(drop=True)xml_file_dataset.columns = ['xml', 'uid_str']xml_file_dataset

根据xml的UID对应上dcm文件：
下面这段代码跑起来需要几分钟的时间，所以跑完后可以保存一下csv文件

xml_file_dataset['dcm'] = xml_file_dataset['uid_str'].apply(lambda x : df_dcm.loc[df_dcm['uid_str']==x].values[0][0])xml_file_dataset.to_csv('xml_file_dataset.csv')# 以上代码需要跑8分钟，所以保存起来xml_file_dataset = pd.read_csv('xml_file_dataset.csv', index_col=0)xml_file_dataset

传入Dataset。

img = xml_file_dataset['dcm'].valueslabel = xml_file_dataset['xml'].valueslungdetection_dataset = LungDetection(img=img, label=label, transform=transform)

顺利完成！