项目笔记《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)

这段代码属实是长，慢慢看。　　

主要分为以下几步。

1. 加载原始数据和掩码，用的是load_itk_image函数
2. 求取掩码的边界，即非零部分的边缘，求出一个box，然后对其应用新的分辨率，也就是重采样，将分辨率统一，采用的函数是resample
3. 将数据clip至-1200~600，此范围外的数据置为-1200或600，然后再将数据归一化至0~255，采用的是lum_trans函数
4. 对掩码进行一下膨胀操作，去除肺部的小空洞，采用的函数是process_mask，然后对原始数据应用新掩码，并将掩码外的数据值为170（水的HU值经过归一化后的新数值）
5. 将原始数据重采样，再截取box内的数据即可。
6. 读取标签，将其转换为体素坐标，采用的函数是worldToVoxelCoord，再对其应用新的分辨率，最后注意，数据是box内的数据，所以坐标是相对box的坐标。
7. 将预处理后的数据和标签以.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