V-net实例-数据处理：
参考代码链接：https://github.com/menjingru/vnet
b站链接： https://b23.tv/dgFY3TY
数据描述
记录肺结节类3D数据处理
数据集：LUNA16 + 对应的一个模板数据
观察数据：LUNA16中有很多的subset文件夹，每个文件夹中包含很多对<.mhd,.raw>文件，其中.mhd文件中包含的源图片的信息：
①对应的.raw地址 ②结节坐标<x,y,z> ③结节半径 ④每个坐标系的体素间距等等
annotations.csv文件的表头是： 图片名称、结节位置<x,y,z>、半径
模板数据：存储的就是每个图片肺结节的每个坐标点，用来画标签的。
数据处理目标
- 制作标签
- 图与肺膜相乘
- 对图和肺膜平均像素间距（重采样）
- 将图归一化、去均值
制作标签
- 得到有结节的图的信息，存储到Excel文件中
bbox_annos_():产生bbox_annos文件，处理所有图的坐标转换
def bbox_annos_():
c = np.array(pd.read_csv(annos_csv))#这个就是读取那个annotations.csv
d = [] #准备放坐标转换后的["name",[[结点1],[结点2]...]]
for i in range(10): #默认10个subset都下完了
file_list = os.listdir(luna_path + fengefu+"subset%d"%i) #打开subset文件夹
for ii in file_list: #遍历 如 subset0内所有文件
if len(ii.split(".m"))==2:#如果文件名是mhd的话
name = ii.split(".m")[0] #取出文件名，去掉后缀，得到图名
ct_image_path = find_mhd_path(name) #把文件名拿去找对应的mhd文件的绝对地址
#spacing是像素间隔 fanzhuan是是否需要翻转
numpyImage,origin,spacing,fanzhuan = read_data(ct_image_path) #读取这个mhd文件
#get_raw_label(图片名字，numpyImage，读取的annotations.csv，原点，fanzhuan？？)
one_annos = get_raw_label(name,numpyImage,c,origin,fanzhuan) #进行坐标转换
d.append([name,one_annos])
bbox_annos = pd.DataFrame(d)#把d转换为Excel文件
bbox_annos.to_excel(new_bbox_annos_path)#保存到new_bbox_annos_path
find_mhd_path(name):在bbox_annos_中实际上直接用ii 就好了。这里封装了一个函数，用于后面的根据名字寻找对应的mhd文件路径
def find_mhd_path(name1):
for file_list in os.listdir(luna_path):#遍历luna16文件夹下所有文件
if file_list.find("subset")!=-1: #在有subset的文件夹下查找，这一句是为了避免找到别的例如seg-lungs-LUNA16文件夹里边去
for ii in os.listdir(luna_path+fengfu+file_list):#打开luna16文件夹下的文件夹比如subset0，遍历里面的文件
if len(ii.split(".m"))>1:#如果文件中有".m"字符，len为2，分为两部分，len>1
if ii.split(".m")[0] == name1:#找到图片所在的mhd文件 拿着图片名字对应着找
path = luna_path + fengefu +file_list +fengefu +ii
print(path)
return path
read_data(mhdPath):读取mhd文件的图像数据
def read_data(mhd_file): #从mhd读取图像数据
with open(mhd_file) as f:
mhd_data = f.readlines()
for i in mhd_data: #判断是否反转，其中transformerMatrix = 1 0 0 0 1 0 0 0 1\n 代表为正true #因为100代表x,010代表y，001代表z
if i.startwith('TransformMatrix'):#取得以transformerMatrix开头的这一行
zheng = i.split('=')[1]
if zheng == '1 0 0 0 1 0 0 0 1\n':
zheng = True #代表是正的，不需要反转
if zheng!= True:
zheng = False
itkimage = sitk.ReadImage(mhd_file)#读取mhd文件
numpyImage = sitk.GetArrayFromImage(itkimage)#从mhd读取raw也就是图 DEPTH WIDTH HEIGHT
origin = itkimage.GetOrigin()#从mhd读取origin，也就是原点坐标
spacing = itkimage.GetSpacing() #像素间隔
return numpyImage,origin,spacing,zheng
get_raw_label(img_name,img,annos,origin,zheng):坐标变换
把结节坐标从世界坐标转换到图像坐标([1,1,1]像距)
def get_raw_label(img_name,img,annos,origin,zheng):
#一个图好几个结节
annos_these_list = [] #准备装["name",[[结点1],[结点2]...]]装所有结节的图的名字和坐标点（原始）
for i in range(len(annos)):#遍历annotations所有数据
if annos[i][0] == img_name: #如果名字相符
annos_these_list.append(list(annos[i]))#["name",x,y,z,diam],["name",x,y,z,diam]最后一个是直径
return_list = []
for one_annos_list in annos_these_list: #打开第一个结节数据
w_center = [one_annos_list[1],one_annos_list[2],one_annos_list[3]]#世界坐标xyz
v_center = list(abs(w_center-np.array(origin)))#/np.array(spacing)像素间隔为[1,1,1]因此不用再除 abs是绝对值
if zheng id False:
#所有和图相关的都是zyx的 -1是因为从0开始 -v_center[0]减去x的坐标就是翻转 z轴不需要反转
v_center = [img.shape[2]-1-v_center[0],img.shape[1]-1-v_center[1],v_center[2]]
diam = one_annos_list[4]#直径
one_annos = []
one_annos.append(v_center[0])#x
one_annos.append(v_center[1])#y
one_annos.append(v_center[2])#z
one_annos.append(diam/2)#半径
return_list.append(one_annos)#收集这个结节到return_list
return return_list
结节的位置：原始结节位置是<x,y,z>并且是基于3D图像的坐标值，每张3D图像是基于世界坐标原点需要变换到图像坐标中，所以：
v_center = list(abs(w_center-np.array(origin)))代表结节坐标减去世界坐标原点
img中是<z,y,x>坐标系，还需要进行坐标系的反转
实质上就是进行一个坐标系的变换
小结：遍历mhd文件里面的信息，把annotations.csv里面的结节按照图片名称分类，组成一个新的csv文件(相当于将mhd和annotations.csv进行联表查询)
这个新的csv文件包含了结节在图像中的位置，以及半径大小
- 根据带结节的图名获取对应<.xml,.mhd>文件
find_xml_path(imgName):寻找图名对应的xml文件
def find_xml_path(name1):
list1 = []
for file_list in os.listdir(xml_file_path): #找到xml_file_path文件夹下的所有文件
for ii in os.listdir(xml_file_path + fengefu + file_list): #遍历
aim_path = xml_file_path + fengefu + file_list + ii
with open(aim_path) as f:
if name(f) == name1:
path = aim_path #绝对地址
list1.append(path)
if list1 != []:
return list1
name(xmlPath):从xml文件中拿到本文件中代表的图名
def name(xml_path):#从xml文件中读取name
child = "ResponseHeader" #响应 头
child_child = "SeriesInstanceUid" #案例名
child_child1 = "CTSeriesInstanceUid" #案例名 之所以有两个是因为dataset标注不规范
dom = xml.dom.minidom.parse(xml_path) #读取xml文件
root = dom.documentElement #获得树根
a = root.getElementsByTagName(child) #获取树根下的ResponseHeader
child_node = a[0].getElementsByTagName(child_child)
if child_node == []:
child_node = a[0].getElementsByTagName(child_child1)
child_value = child_node[0].childNodes[0].nodeValue
return child_value
小结：这个部分图片对应的xml文件，简单的遍历实现
- 根据<.xml,.mhd>文件画mask制作标签
利用mhd文件中的图片，得到图片的大小作为mask的画布（3D的，就是3维的），拿到确定大小的画布后，根据xml文件中结节坐标画mask。一张一张的画mask：数据是3D图像，怎么画？
把z看做层，每一层是2D的，在xml文件中把结节位置拥有同一个z的拿出来，每层每层的画
#一个处理每张图像的函数
#输入名字+一个空列表 输出上色好的mask + mhd的绝对地址用于后续存mask用 + 出错的图的名字
def for_one_(name,wrong):
xml_path_list = find_xml_path(name)
ct_image_path = find_mhd_path(name)
ct_image,origin,spacing,fanzhuan = read_data(ct_image_path)#从mhd文件读取信息
s = ct_image.shape
mm = np.zeros((s[0],s[1],s[2]),dtype=np.int32) #mm为全0的mask 注意 图.shape是zyx的 所以顺序不用变
#取截面 描点
for i in xml_path_list: #取得xml文件的绝对地址
list1 = point(i,origin[2]) #在这个绝对地址内获取所有[[z层，点],[z层，点],...]
for ii in list1: #遍历所有层
ceng = ii[0]
pts = ii[1]#该层所有的点位
color = 1 #(0,255,0)
#开始着色 层/像素间距 从0开始 -1 thickness不往外漫延
mm[int(ceng/spacing[2]-1),:,:] = cv.drawContours(mm[int(ceng/spacing[2]-1),:,:],[pts],-1,color=color,thickness = -1)
#弥补空洞 就是这个结节内部空洞部分进行填充
mm[int(ceng/spacing[2]-1),:,:] = scipy.ndimage.binary_fill_holes( mm[int(ceng/spacing[2]-1),:,:],structure=None,output=None,origin=0)
#.all()代表完全一致 这个wrong是总函数的位置 它的意思就是这个wrong能够存储所有图片的错误，不只是自己这张图片的错误
if (mm==np.zeros((s[0],s[1],s[2]),dtype=np.int32)).all():#如果没染上色 即仍是全0数组
wrong.append(name)
return mm,ct_image_path,wrong
point(xml_path,origin2):返回每层的结节位置<x,y> 就是对xml文件找对应标签里面的数据
def point(xml_path,origin2):#需要xml_path和图像的原点坐标origin的z轴坐标
a = [] #该图案的所有z轴和z轴上点位的列表 [ [z1,[[x1,y1],[x2,y2],...]],[z2,[[x1,y1],[x2,y2],...]],...]
dom = xml.dom.minidom.parse(xml_path) #读取xml文件
root = dom.documentElement #获得树根
nodeid = root.getElementsByTagName("readingSession")
for u in nodeid:
child = u.getElementsByTagName("unblindedReadNodule")
for i in child:
id = i.getElementsByTagName("noduleID")
id1 = id[0].childNodes[0].nodeValue
for id1: #如果noduleID值不为0
one_all_iou = i.getElementsByTagName("roi")
for i in one_all_iou: #遍历这些roi点
z = r.getElementsByTagName("imageZposition")
z1 = float(z[0].childNodes[0].nodeValue)-origin2 #得到图像坐标
ioux = r.getElementByTagName("xCoord") #取得该z轴切片上的xCoord
iouy = r.getElementByTagName("yCoord") #取得该z轴切片上的yCoord
ioux1 = np.array([int(k.childNodes[0].nodeValue) for k in ioux])#变成数组
iouy2 = np.array([int(k.childNodes[0].nodeValue) for k in iouy])#变成数组
iou = np.array([ioux1,iouy1])#数组合并 得到[[x1,x2,...],[y1,y2,...]]
point1 = np.transpose(iou) #转置 得到[[x1,y1],[x2,y2],...]
a.append([z1,point1]) #在这一层的z的坐标上的结节的轮廓
return a
小结：理解数据是3D的，分层上色
重采样（平均体素间距）
def resample(imgs,spacing,new_spacing=[1,1,1]):#重采样 把原图像的像素间隔统一
###重采样 坐标原点位置为0
if len(img.shape)==3:#如果是三维的话
new_shape=[]
for i in range(3):#对每个维度0 1 2 -》 z y x
new_zyx = np.round(imgs.shape[i]*spacing[-i-1]/new_spacing[-i-1])#原顺序为xyz spacing[-i-1]顺序为zyx 不除也行new_spacing=[1,1,1]
new_shape.append(new_zyx)
resize_factor=[] #新图尺寸/原图尺寸 即缩放比例 如原像素间隔为2.5 新像素间隔为1 缩放比例为1/2.5
for i in range(3): #z,y,x
resize_zyx = new_shape[i]/imgs.shape[i]
resize_factor.append(resize_zyx)
imgs = zoom(imgs,resize_factor,mode='_nearest')#放缩 边缘使用最近邻 插值默认使用三线性插值
return imgs
else:
raise ValueError('wrong shape') #本代码只能处理三维数据
这部分有个概念：体素间距
在肺结节这种3D图像中，数据.mhd里面本身就有一个spacing，这个就代表体素间距
每个体素代表的大小不是1，他要描述整个器官的大小，那每个体素之间的距离就拿来描述。但是把数据拿来训练，每个3D图有自己的体素间距显然不合适，所以把这个体素间距统一，可以看成是对每个维度按照体素间距进行比例缩放。
上面的代码resize_factor可以直接用spacing来表示，但是这里绕了一圈（未知原因）
主函数
anno_name_list = annos()#有结节的图的名字
print(len(anno_name_list))
wrony = [] #染色失败的标签
for name in anno_name_list: #遍历所有结节图的名字
mask,ct_image_path,wrong = for_one_(name,wrong)#染色
path = ct_image_path.split("LUNA16")[1].split(".m")[0] #取LUNA16后 .mhd前的字符串
if xitong = "linux":
path = path.replace(r"\s","/s")
path = path.replace(r"\1","/1")
else:
pass
print(path)
#path形式 ： \subset1\imagename没有mhd
#图
ct_image,origin,spacing,isflip = read_data(ct_image_path)
#mask 肺部掩膜
ct_image1,origin1,spacing1,isflip1 = read_data(mask_path + fengefu + name + ".mhd")
ct_image1[ct_image1>1]=1 #LUNA16提供的肺部掩膜分左右膜 左肺为右肺为4 需要统一1
ct_image = ct_image*ct_image1 #将图与掩膜相乘，肺外区域归0
image = resample(ct_image,spacing) #图 重采样
mask = resample(mask,spacing)# 标签 重采样
#LUNA16竞赛中常用来做归一化处理的阈值时-1000和400
max_num = 400
min_num = -1000
image = (image - min_num)/(max_num-min_num)#归一化公式
image[image > 1] = 1. #高于1 的归为1 float格式
image[image <0 ] = 0. #低于0 的归为0 float格式
#LUNA16竞赛中的均值大约是0.25
img = image - 0.25 #去均值 图的操作完成
#构建文件
if not os.path.exists(output_path):
os.mkdir(output_path)
if not os.path.exists(output_path+fengefu+"bbox_image"):
os.mkdir(output_path+fengefu+"bbox_image")
if not os.path.exists(output_path+fengefu+"bbox_mask"):
os.mkdir(output_path+fengefu+"bbox_mask")
sub_path = ct_image_path.split("LUNA16")[1].split("1.")[0] #\subset1\
if not os.path.exists(output_path+fengefu+"bbox_image"+sub_path):
os.mkdir(output_path+fengefu+"bbox_image"+sub_path)
if not os.path.exists(output_path+fengefu+"bbox_mask"+sub_path):
os.mkdir(output_path+fengefu+"bbox_mask"+sub_path)
np.save(output_path+fengefu+"bbox_image"+path,img)
np.save(output_path+fengefu+"bbox_mask"+path,mask)
wrong_img = pandas.DataFrame(wrong)#保存未染色图
wrong_img.to_excel(wrong_img_path)#保存到wrong_img_path
其实不太理解为什么mask在point中z轴除了spacing，然后在重采样的时候还除了一次
总结
- 首先获取所有有结节的图片列表
- 按照图名，获取标签
- 图像进行归一化和去均值
不好理解的地方就是体素间距和坐标系转换。
1332
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
