项目笔记《DeepLung:Deep 3D Dual Path Nets for Automated Pulmonary Nodule Detection and Classification》（一）预处理
最近一个月都在做肺结节的检测，学到了不少东西，运行的项目主要是基于这篇论文，在github上可以查到项目代码。
我个人总结的肺结节检测可以分为三个阶段，数据预处理，网络搭建及训练，结果评估。
这篇博客主要分析一下项目预处理部分的代码实现。
预处理的全部代码都在prepare.py中，对原始数据进行处理，输出预处理后的数据。
首先是主函数
def preprocess_luna():
luna_segment = config['luna_segment']#存放CT掩码的路径
savepath = config['preprocess_result_path']#存放预处理后数据的路径
luna_data = config['luna_data']#LUNA16的原始数据
luna_label = config['luna_label']#存放所有病例标签的文件
finished_flag = '.flag_preprocessluna'#是否已经预处理过的标志
print('starting preprocessing luna')
if not os.path.exists(finished_flag):
annos = np.array(pandas.read_csv(luna_label))
pool = Pool()#开启线程池
if not os.path.exists(savepath):
os.mkdir(savepath)
for setidx in xrange(10):#十份数据
print 'process subset', setidx
filelist = [f.split('.mhd')[0] for f in os.listdir(luna_data+'subset'+str(setidx)) if f.endswith('.mhd') ]#原始数据为.mhd文件，只保留文件名，去掉.mhd后缀
if not os.path.exists(savepath+'subset'+str(setidx)):
os.mkdir(savepath+'subset'+str(setidx))#为每份数据创建存放预处理结果的文件夹
partial_savenpy_luna = partial(savenpy_luna, annos=annos, filelist=filelist,#函数修饰器，将一些参数预先设定，后面调用更简洁
luna_segment=luna_segment, luna_data=luna_data+'subset'+str(setidx)+'/',
savepath=savepath+'subset'+str(setidx)+'/')
N = len(filelist)
#savenpy(1)
_=pool.map(partial_savenpy_luna,range(N))#将函数调用在序列的每个元素上，返回一个含有所有返回值的列表
pool.close()#关闭线程池
pool.join()
print('end preprocessing luna')
f= open(finished_flag,"w+")#预处理结束，写入结束标志
上面的代码就是预处理的全部代码，当然里面还调用了其它函数，主要的流程就是针对十份数据中的每一份，针对每一份中的每一个case(CT图像，以.mhd格式存储），分别处理，为加快速度，在每份数据中，针对每个文件分别开启一个线程，我试过不用线程，采取循环处理，速度确实慢了很多。
上面的代码中最关键的就是数据预处理函数savenpy_luna，在主函数中已经设定好一部分参数，如文件名列表，标签，掩码，原始数据，预处理存储路径，万事俱备，只欠实现，接下来就看一下savenpy_luna的代码。
def savenpy_luna(id, annos, filelist, luna_segment, luna_data,savepath):
islabel = True
isClean = True
resolution = np.array([1,1,1])
# resolution = np.array([2,2,2])
name = filelist[id]
sliceim,origin,spacing,isflip = load_itk_image(os.path.join(luna_data,name+'.mhd'))#加载原始数据
Mask,origin,spacing,isflip = load_itk_image(os.path.join(luna_segment,name+'.mhd'))#加载相应的掩码
if isflip: #这一步没看懂
Mask = Mask[:,::-1,::-1]
newshape = np.round(np.array(Mask.shape)*spacing/resolution).astype('int')#获取mask在新分辨率下的尺寸
m1 = Mask==3 #LUNA16的掩码有两种值，3和4
m2 = Mask==4
Mask = m1+m2 #将两种掩码合并
xx,yy,zz= np.where(Mask) #确定掩码的边界
box = np.array([[np.min(xx),np.max(xx)],[np.min(yy),np.max(yy)],[np.min(zz),np.max(zz)]])
box = box*np.expand_dims(spacing,1)/np.expand_dims(resolution,1) #对边界即掩码的最小外部长方体应用新分辨率
box = np.floor(box).astype('int')
margin = 5
extendbox = np.vstack([np.max([[0,0,0],box[:,0]-margin],0),np.min([newshape,box[:,1]+2*margin],axis=0).T]).T #对box留置一定空白
this_annos = np.copy(annos[annos[:,0]==(name)]) #读取该病例对应标签
if isClean:
convex_mask = m1
dm1 = process_mask(m1) #对掩码采取膨胀操作，去除肺部黑洞
dm2 = process_mask(m2)
dilatedMask = dm1+dm2
Mask = m1+m2
extramask = dilatedMask ^ Mask
bone_thresh = 210
pad_value = 170
if isflip:
sliceim = sliceim[:,::-1,::-1]
print('flip!')
sliceim = lumTrans(sliceim) #对原始数据阈值化，并归一化
sliceim = sliceim*dilatedMask+pad_value*(1-dilatedMask).astype('uint8') #170对应归一化话后的水，掩码外的区域补充为水
bones = (sliceim*extramask)>bone_thresh #210对应归一化后的骨头，凡是大于骨头的区域都填充为水
sliceim[bones] = pad_value
sliceim1,_ = resample(sliceim,spacing,resolution,order=1) #对原始数据重采样，即采用新分辨率
sliceim2 = sliceim1[extendbox[0,0]:extendbox[0,1], #将extendbox内数据取出作为最后结果
extendbox[1,0]:extendbox[1,1],
extendbox[2,0]:extendbox[2,1]]
sliceim = sliceim2[np.newaxis,...]
np.save(os.path.join(savepath, name+'_clean.npy'), sliceim)
np.save(os.path.join(savepath, name+'_spacing.npy'), spacing)
np.save(os.path.join(savepath, name+'_extendbox.npy'), extendbox)
np.save(os.path.join(savepath, name+'_origin.npy'), origin)
np.save(os.path.join(savepath, name+'_mask.npy'), Mask)
if islabel:
this_annos = np.copy(annos[annos[:,0]==(name)]) #一行代表一个结节，所以一个病例可能对应多行标签
label = []
if len(this_annos)>0:
for c in this_annos:
pos = worldToVoxelCoord(c[1:4][::-1],origin=origin,spacing=spacing) #将世界坐标转换为体素坐标
if isflip:
pos[1:] = Mask.shape[1:3]-pos[1:]
label.append(np.concatenate([pos,[c[4]/spacing[1]]]))
label = np.array(label)
if len(label)==0:
label2 = np.array([[0,0,0,0]]) #若没有结节则设为全0
else:
label2 = np.copy(label).T
label2[:3] = label2[:3]*np.expand_dims(spacing,1)/np.expand_dims(resolution,1) #对标签应用新的分辨率
label2[3] = label2[3]*spacing[1]/resolution[1] #对直径应用新的分辨率
label2[:3] = label2[:3]-np.expand_dims(extendbox[:,0],1) #将box外的长度砍掉，也就是相对于box的坐标
label2 = label2[:4].T
np.save(os.path.join(savepath,name+'_label.npy'), label2)
print(name)
这段代码属实是长，慢慢看。
主要分为以下几步。
- 加载原始数据和掩码，用的是load_itk_image函数
- 求取掩码的边界，即非零部分的边缘，求出一个box，然后对其应用新的分辨率，也就是重采样，将分辨率统一，采用的函数是resample
- 将数据clip至-1200~600，此范围外的数据置为-1200或600，然后再将数据归一化至0~255，采用的是lum_trans函数
- 对掩码进行一下膨胀操作，去除肺部的小空洞，采用的函数是process_mask，然后对原始数据应用新掩码，并将掩码外的数据值为170（水的HU值经过归一化后的新数值）
- 将原始数据重采样，再截取box内的数据即可。
- 读取标签，将其转换为体素坐标，采用的函数是worldToVoxelCoord，再对其应用新的分辨率，最后注意，数据是box内的数据，所以坐标是相对box的坐标。
- 将预处理后的数据和标签以.npy格式存储
针对上面用到的各种工具函数，分别解析下。
加载图像：主要用到sitk模块，读取原始数据的numpy表示，以及origin和space，origin就是真实坐标的原点，space就是分辨率
def load_itk_image(filename):
with open(filename) as f:
contents = f.readlines()
line = [k for k in contents if k.startswith('TransformMatrix')][0]
transformM = np.array(line.split(' = ')[1].split(' ')).astype('float')
transformM = np.round(transformM)
if np.any( transformM!=np.array([1,0,0, 0, 1, 0, 0, 0, 1])):
isflip = True
else:
isflip = False
itkimage = sitk.ReadImage(filename)
numpyImage = sitk.GetArrayFromImage(itkimage)
numpyOrigin = np.array(list(reversed(itkimage.GetOrigin())))
numpySpacing = np.array(list(reversed(itkimage.GetSpacing())))
return numpyImage, numpyOrigin, numpySpacing,isflip
重采样：原始CT分辨率往往不一致，为便于应用网络，需要统一分辨率
def resample(imgs, spacing, new_spacing,order=2):
if len(imgs.shape)==3:
new_shape = np.round(imgs.shape * spacing / new_spacing)
true_spacing = spacing * imgs.shape / new_shape
resize_factor = new_shape / imgs.shape
imgs = zoom(imgs, resize_factor, mode = 'nearest',order=order)
return imgs, true_spacing
elif len(imgs.shape)==4:
n = imgs.shape[-1]
newimg = []
for i in range(n):
slice = imgs[:,:,:,i]
newslice,true_spacing = resample(slice,spacing,new_spacing)
newimg.append(newslice)
newimg=np.transpose(np.array(newimg),[1,2,3,0])
return newimg,true_spacing
else:
raise ValueError('wrong shape')
坐标转换：给定的标签是世界坐标，单位是mm，需要转换为体素坐标，也就是在像素体内的坐标
def worldToVoxelCoord(worldCoord, origin, spacing):
stretchedVoxelCoord = np.absolute(worldCoord - origin)
voxelCoord = stretchedVoxelCoord / spacing
return voxelCoord
掩码处理：这里对掩码进行膨胀处理
def process_mask(mask):
convex_mask = np.copy(mask)
for i_layer in range(convex_mask.shape[0]):
mask1 = np.ascontiguousarray(mask[i_layer])
if np.sum(mask1)>0:
mask2 = convex_hull_image(mask1)
if np.sum(mask2)>1.5*np.sum(mask1):
mask2 = mask1
else:
mask2 = mask1
convex_mask[i_layer] = mask2
struct = generate_binary_structure(3,1)
dilatedMask = binary_dilation(convex_mask,structure=struct,iterations=10)
return dilatedMask
归一化：将数据归一化至0~255
def lumTrans(img):
lungwin = np.array([-1200.,600.])
newimg = (img-lungwin[0])/(lungwin[1]-lungwin[0])
newimg[newimg<0]=0
newimg[newimg>1]=1
newimg = (newimg*255).astype('uint8')
return newimg
