在前面进行了肺结节数据的预处理之后，接下来开始进入肺结节检测环节。首先附上该项目的Github链接：https://github.com/Minerva-J/DeepLung。

检测器训练

之前我们对10份数据进行划分，将分别对10折数据进行独立的9份训练和1份测试实验，以其中一个实验为例：

首先修改config_training0.py，将subset0-8作为训练集，subset9作为验证集和测试集。

代码如下：

config = {'train_data_path':['/home/zhaojie/zhaojie/Lung/data/luna16/subset_data/subset0/',
                             '/home/zhaojie/zhaojie/Lung/data/luna16/subset_data/subset1/',
                             '/home/zhaojie/zhaojie/Lung/data/luna16/subset_data/subset2/',
                             '/home/zhaojie/zhaojie/Lung/data/luna16/subset_data/subset3/',
                             '/home/zhaojie/zhaojie/Lung/data/luna16/subset_data/subset4/',
                             '/home/zhaojie/zhaojie/Lung/data/luna16/subset_data/subset5/',
                             '/home/zhaojie/zhaojie/Lung/data/luna16/subset_data/subset6/',
                             '/home/zhaojie/zhaojie/Lung/data/luna16/subset_data/subset7/',
                             '/home/zhaojie/zhaojie/Lung/data/luna16/subset_data/subset8/'],
          'val_data_path':['/home/zhaojie/zhaojie/Lung/data/luna16/subset_data/subset9/'], 
          'test_data_path':['/home/zhaojie/zhaojie/Lung/data/luna16/subset_data/subset9/'],  

          'train_preprocess_result_path':'/home/zhaojie/zhaojie/Lung/DeepLung-Minerva/Data/LUNA16PROPOCESSPATH/', # contains numpy for the data and label, which is generated by prepare.py
          'val_preprocess_result_path':'/home/zhaojie/zhaojie/Lung/DeepLung-Minerva/Data/LUNA16PROPOCESSPATH/',
          'test_preprocess_result_path':'/home/zhaojie/zhaojie/Lung/DeepLung-Minerva/Data/LUNA16PROPOCESSPATH/',

          'train_annos_path':'/home/zhaojie/zhaojie/Lung/data/luna16/CSVFILES/annotations.csv',
          'val_annos_path':'/home/zhaojie/zhaojie/Lung/data/luna16/CSVFILES/annotations.csv',
          'test_annos_path':'/home/zhaojie/zhaojie/Lung/data/luna16/CSVFILES/annotations.csv',

          'black_list':[],

          'preprocessing_backend':'python',
         } 

python main.py --model dpn3d26 -b 8 --save-dir dpn3d26/retrft96X/ --epochs 1000 --config config_training0.py

You should modify X from 0-9 to train 10 models for 10 fold cross validation

You can modify -b(batch_size) depend on your GPU memory and number.

重点介绍一下数据读取和数据增强部分

通过DataBowl3Detector类函数的crop和label_mapping.从clean.npy的体数据中截取96*96*96的立体数据和制作对应的立体标签24*24*24*3*5. 5代表4+1,4代表xyzd的精确修正值，1代表用来判断前景还是背景（是否是肺结节）

肺结节标签包含两部分：第一部分就是70%裁剪96*96*96的立方体，但是结节随机的大致处在中心，而非都在正中心。其次还要分为放大和缩小的两种。根据结节大小缩放corp_size（比如96要缩放1.2被，那么corp_size=80），再截取80*80*80(如果结节在边缘不足80，用170补全)，接着在缩放回96*96*96.另一种30%随机裁剪到96*96*96。

DataBowl3Detector代码如下：

class DataBowl3Detector(Dataset):
    def __init__(self, data_dir, split_path, config, phase='train', split_comber=None):
        assert(phase == 'train' or phase == 'val' or phase == 'test')
        self.phase = phase
        self.max_stride = config['max_stride']       
        self.stride = config['stride']#4
        sizelim = config['sizelim']/config['reso']#2.5/1
        sizelim2 = config['sizelim2']//config['reso']#10/1
        sizelim3 = config['sizelim3']//config['reso']#20/1
        self.blacklist = config['blacklist']
        self.isScale = config['aug_scale']
        self.r_rand = config['r_rand_crop']
        self.augtype = config['augtype']#{'flip':True,'swap':False,'scale':True,'rotate':False}
        self.pad_value = config['pad_value']
        self.split_comber = split_comber
        idcs = split_path # np.load(split_path)

        if phase!='test':
            idcs = [f for f in idcs if (f not in self.blacklist)]

        self.filenames = [os.path.join(data_dir, '%s_clean.npy' % idx) for idx in idcs]
        # print self.filenames
        self.kagglenames = [f for f in self.filenames]# if len(f.split('/')[-1].split('_')[0])>20]
        # self.lunanames = [f for f in self.filenames if len(f.split('/')[-1].split('_')[0])<20]
        # print('-------',idcs)
        labels = []

        # print('len(idcs)',len(idcs))
        for idx in idcs:
            # print(data_dir, idx)
            l = np.load(data_dir+idx+'_label.npy', allow_pickle = True)
            # l = np.load(data_dir+idx+'_label.npy')
            if np.all(l==0):#without nodule
                l=np.array([])
            labels.append(l)
        # print('Warning')
        self.sample_bboxes = labels#nodule
        if self.phase != 'test':
            self.bboxes = []
            for i, l in enumerate(labels):#[[81.933381 166.26051683000003 43.932055500000004 6.440878725]]
                # print('-------',self.bboxes)
                if len(l) > 0 :
                    for t in l:
                        if t[3]>sizelim:#直径>2.5
                            self.bboxes.append([np.concatenate([[i],t])])
                        if t[3]>sizelim2:#直径>10
                            self.bboxes+=[[np.concatenate([[i],t])]]*2
                        if t[3]>sizelim3:#直径>20
                            self.bboxes+=[[np.concatenate([[i],t])]]*4
            self.bboxes = np.concatenate(self.bboxes,axis = 0)# many boxes
        # print('-------',self.bboxes, self.sample_bboxes)
        self.crop = Crop(config)
        self.label_mapping = LabelMapping(config, self.phase)

    def __getitem__(self, idx,split=None):
        t = time.time()
        np.random.seed(int(str(t%1)[2:7]))#seed according to time
        isRandomImg  = False
        if self.phase !='test':
            # print('idx, len(self.bboxes)',idx, len(self.bboxes))
            if idx>=len(self.bboxes):
                isRandom = True
                idx = idx%len(self.bboxes)
                isRandomImg = np.random.randint(2)#随机产生0和1
            else:
                isRandom = False
        else:
            isRandom = False
        # print('Warning1')
        if self.phase != 'test':
            if not isRandomImg:
                bbox = self.bboxes[idx]#num,x,y,z,r,第一个参数代表的是第几个CT序列（trainlist排列好的）后面代表的是xyz和直径。当直径大于2.5写1行，大于10写1+2行，大于20写1+2+4行
                filename = self.filenames[int(bbox[0])]#第int(bbox[0])序列的路径（包含名称_clearn.npy）
                imgs = np.load(filename)

                bboxes = self.sample_bboxes[int(bbox[0])]#self.sample_bboxes就是800CT对应的结节的xyz和直径，是一个列表list（可能是空，也可能是多行）
                isScale = self.augtype['scale'] and (self.phase=='train')#True
                # print('---bbox---',bbox)
                # print('---bboxes---',isRandom)
                #imgs是_clearn.npy导入的3d矩阵(仅仅包含肺部并且分辨率为1×1*1mm)/bbox代表[556 74.7895776 249.75681577 267.72357450000004 10.57235285],bbox[1:]代表[74.7895776 249.75681577 267.72357450000004 10.57235285]
                # bboxes代表[[74.7895776 249.75681577 267.72357450000004 10.57235285],[165.8177988 172.30634478 79.42561491000001 20.484109399999998]]等于或包含bbox[1:].
                sample, target, bboxes, coord = self.crop(imgs, bbox[1:], bboxes,isScale,isRandom)
                if self.phase=='train' and not isRandom:
                     sample, target, bboxes, coord = augment(sample, target, bboxes, coord,
                        ifflip = self.augtype['flip'], ifrotate=self.augtype['rotate'], ifswap = self.augtype['swap'])
            else:
                randimid = np.random.randint(len(self.kagglenames))
                filename = self.kagglenames[randimid]
                imgs = np.load(filename)
                bboxes = self.sample_bboxes[randimid]
                isScale = self.augtype['scale'] and (self.phase=='train')
                sample, target, bboxes, coord = self.crop(imgs, [], bboxes,isScale=False,isRand=True)
            # print('Warning2')
            label = self.label_mapping(sample.shape[1:], target, bboxes, filename)
            sample = (sample.astype(np.float32)-128)/128

            # print('label',np.unique(label),label.shape,label.dtype)#array([-1.,  0.], dtype=float32), (24, 24, 24, 3, 5), dtype('float32')
            # print('coord',coord.shape, coord.dtype)#(3, 24, 24, 24)
            #if filename in self.kagglenames and self.phase=='train':
            #    label[label==-1]=0
            return torch.from_numpy(sample), torch.from_numpy(label), coord
        else:#test
            imgs = np.load(self.filenames[idx])
            # print('imgs.shape',imgs.shape)
            bboxes = self.sample_bboxes[idx]
            nz, nh, nw = imgs.shape[1:]#(197, 181, 224)
            pz = int(np.ceil(float(nz) / self.stride)) * self.stride
            ph = int(np.ceil(float(nh) / self.stride)) * self.stride
            pw = int(np.ceil(float(nw) / self.stride)) * self.stride
            imgs = np.pad(imgs, [[0,0],[0, pz - nz], [0, ph - nh], [0, pw - nw]], 'constant',constant_values = self.pad_value)#能被stride整除
            # print('1------pad',self.stride, pz,ph,pw)
            xx,yy,zz = np.meshgrid(np.linspace(-0.5,0.5,imgs.shape[1]/self.stride),
                                   np.linspace(-0.5,0.5,imgs.shape[2]/self.stride),
                                   np.linspace(-0.5,0.5,imgs.shape[3]/self.stride),indexing ='ij')
            coord = np.concatenate([xx[np.newaxis,...], yy[np.newaxis,...],zz[np.newaxis,:]],0).astype('float32')
            # print('------1',imgs.shape)
            imgs, nzhw = self.split_comber.split(imgs)
            # print('------2',imgs.shape)
            # print('------',self.split_comber.max_stride,self.stride)

            coord2, nzhw2 = self.split_comber.split(coord,
                                                   side_len = self.split_comber.side_len//self.stride,
                                                   max_stride = self.split_comber.max_stride//self.stride,
                                                   margin = self.split_comber.margin//self.stride)
            assert np.all(nzhw==nzhw2)
            imgs = (imgs.astype(np.float32)-128)/128
            # print('sample.shape',imgs.shape,coord2.shape)
            return torch.from_numpy(imgs), bboxes, torch.from_numpy(coord2), np.array(nzhw)

    def __len__(self):
        # print('Warning00')
        if self.phase == 'train':
            # print(type(len(self.bboxes)),type(1-self.r_rand))
            # print(len(self.bboxes)/(1-self.r_rand))
            # print(len(self.bboxes),(1-self.r_rand))
            # print('Warning11')
            # return len(self.bboxes)/(1-self.r_rand)
            return int(len(self.bboxes)/(1-self.r_rand))
        elif self.phase =='val':
            return len(self.bboxes)
        else:
            return len(self.sample_bboxes)

比较重要的是Crop函数

代码如下：

class Crop(object):#imgs, bbox[1:], bboxes,isScale,isRandom
    def __init__(self, config):
        self.crop_size = config['crop_size']#[96, 96, 96]
        self.bound_size = config['bound_size']#12
        self.stride = config['stride']#4
        self.pad_value = config['pad_value']#170
    def __call__(self, imgs, target, bboxes,isScale=False,isRand=False):
        if isScale:
            radiusLim = [8.,120.]
            scaleLim = [0.75,1.25]
            scaleRange = [np.min([np.max([(radiusLim[0]/target[3]),scaleLim[0]]),1]) ,np.max([np.min([(radiusLim[1]/target[3]),scaleLim[1]]),1])]#target代表bbox[1:]代表[74.7895776 249.75681577 267.72357450000004 10.57235285]
            scale = np.random.rand()*(scaleRange[1]-scaleRange[0])+scaleRange[0]#np.random.rand(d0, d1, …, dn)的随机样本位于[0, 1)中
            crop_size = (np.array(self.crop_size).astype('float')/scale).astype('int')#根据实际结节直径大小调整crop_size大小,target[3]小crop_size变大.target[3]大crop_size变小
        else:
            crop_size=self.crop_size
        bound_size = self.bound_size
        target = np.copy(target)#目标结节
        bboxes = np.copy(bboxes)#这个ct含有的所有结节
        # print('---crop---',target)
        start = []
        for i in range(3):
            if not isRand:
                r = target[3] / 2
                s = np.floor(target[i] - r)+ 1 - bound_size
                e = np.ceil (target[i] + r)+ 1 + bound_size - crop_size[i] 
            else:#
                s = np.max([imgs.shape[i+1]-crop_size[i]//2,imgs.shape[i+1]//2+bound_size])
                e = np.min([crop_size[i]//2, imgs.shape[i+1]//2-bound_size])
                target = np.array([np.nan,np.nan,np.nan,np.nan])
                # print(s,e)
            if s>e:
                start.append(np.random.randint(e,s))#!
            else:
                start.append(int(target[i])-crop_size[i]//2+np.random.randint(-bound_size//2,bound_size//2))#求取结节的3d立方矩阵最靠近原点的点


        normstart = np.array(start).astype('float32')/np.array(imgs.shape[1:])-0.5#将normstart放到3d块[-0.5:0.5,-0.5:0.5,-0.5:0.5]对应实际[-imgs.shape[1]/2:imgs.shape[1]/2,-imgs.shape[2]/2:imgs.shape[2]/2,-imgs.shape[1]/2:imgs.shape[1]/2,-imgs.shape[3]/2:imgs.shape[3]/2]
        normsize = np.array(crop_size).astype('float32')/np.array(imgs.shape[1:])
        xx,yy,zz = np.meshgrid(np.linspace(normstart[0],normstart[0]+normsize[0],self.crop_size[0]//self.stride),
                           np.linspace(normstart[1],normstart[1]+normsize[1],self.crop_size[1]//self.stride),#np.linspace创建等差数列
                           np.linspace(normstart[2],normstart[2]+normsize[2],self.crop_size[2]//self.stride),indexing ='ij')#可以这么理解，meshgrid函数用两个坐标轴上的点在平面上画网格(3d也是如此)
        coord = np.concatenate([xx[np.newaxis,...], yy[np.newaxis,...],zz[np.newaxis,:]],0).astype('float32')

        pad = []
        pad.append([0,0])
        for i in range(3):
            leftpad = max(0,-start[i])
            rightpad = max(0,start[i]+crop_size[i]-imgs.shape[i+1])
            pad.append([leftpad,rightpad])
        crop = imgs[:,#以结节位置为中心来截取
            max(start[0],0):min(start[0] + crop_size[0],imgs.shape[1]),
            max(start[1],0):min(start[1] + crop_size[1],imgs.shape[2]),
            max(start[2],0):min(start[2] + crop_size[2],imgs.shape[3])]
        # print('---crop0---',crop.shape)
        crop = np.pad(crop,pad,'constant',constant_values =self.pad_value)#越界进行填充
        # print('---crop---',max(start[0],0),min(start[0] + crop_size[0],imgs.shape[1]))
        # print('---crop1---',crop.shape)
        # print(stop)
        for i in range(3):
            target[i] = target[i] - start[i]#目标结节位置已经减去目标结节3d立方矩阵最靠近原点的开始点
        for i in range(len(bboxes)):
            for j in range(3):
                bboxes[i][j] = bboxes[i][j] - start[j]#同一ct的所有结节位置已经减去该目标结节3d立方矩阵最靠近原点的开始点

        if isScale:
            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                crop = zoom(crop,[1,scale,scale,scale],order=1)
            newpad = self.crop_size[0]-crop.shape[1:][0]
            if newpad<0:
                crop = crop[:,:-newpad,:-newpad,:-newpad]
            elif newpad>0:
                pad2 = [[0,0],[0,newpad],[0,newpad],[0,newpad]]
                crop = np.pad(crop, pad2, 'constant', constant_values=self.pad_value)
            for i in range(4):
                target[i] = target[i]*scale
            for i in range(len(bboxes)):
                for j in range(4):
                    bboxes[i][j] = bboxes[i][j]*scale
        # for i in range(crop.shape[1]):
            # name = str(i) + '.png'
            # misc.imsave(os.path.join('../Visualize/crop/', name), crop[0][i])
        # print('target',target)
        # print('bboxes',bboxes)
        return crop, target, bboxes, coord#1*96*96*96/可能是[nan,nan,nan,nan]/可能是空/3*24*24*24

模型结构

接下来介绍一下检测部分的模型，具有DPN结构的3D版本的u-net。

dual path connection如下图

主要架构是DPN，就是dense 与 residual的结合，输入特征图一部分通过residual与输出相加，另一部分与residual的结果再串联。 输入是969696的立方体，里面包含标记的结节，经过24个333的卷积核，通道数变为24，然后经过4个stage，尺寸缩减为1/16，接下来是分辨率放大阶段，采用反卷积实现，连续两个阶段都是反卷积后与低层特征串联，然后经过两个卷积操作，通道数变为15，图示中以3*5显示，是为了更清楚地表明，最后输出的proposal中，每个位置有三个，分别采用三种尺寸，设置的三个anchor尺寸是[5,10,20]，每个位置预测z，y，x，d，p分别是结节的三维坐标以及直径，置信度。

DPN函数代码如下：

class DPN(nn.Module):
    def __init__(self, cfg):
        super(DPN, self).__init__()
        in_planes, out_planes = cfg['in_planes'], cfg['out_planes']
        num_blocks, dense_depth = cfg['num_blocks'], cfg['dense_depth']

        # self.in_planes = in_planes
        # self.out_planes = out_planes
        # self.num_blocks = num_blocks
        # self.dense_depth = dense_depth

        self.conv1 = nn.Conv3d(1, 24, kernel_size=3, stride=1, padding=1, bias=False)
        self.bn1 = nn.BatchNorm3d(24)
        self.last_planes = 24
        self.layer1 = self._make_layer(in_planes[0], out_planes[0], num_blocks[0], dense_depth[0], stride=2)#stride=1)
        self.layer2 = self._make_layer(in_planes[1], out_planes[1], num_blocks[1], dense_depth[1], stride=2)
        self.layer3 = self._make_layer(in_planes[2], out_planes[2], num_blocks[2], dense_depth[2], stride=2)
        self.layer4 = self._make_layer(in_planes[3], out_planes[3], num_blocks[3], dense_depth[3], stride=2)
        self.last_planes = 216
        self.layer5 = self._make_layer(128, 128, num_blocks[2], dense_depth[2], stride=1)
        self.last_planes = 224+3
        self.layer6 = self._make_layer(224, 224, num_blocks[1], dense_depth[1], stride=1)
        self.linear = nn.Linear(out_planes[3]+(num_blocks[3]+1)*dense_depth[3], 2)#10)
        self.last_planes = 120
        self.path1 = nn.Sequential(
            nn.ConvTranspose3d(self.last_planes, self.last_planes, kernel_size = 2, stride = 2),
            nn.BatchNorm3d(self.last_planes),
            nn.ReLU(inplace = True))
        self.last_planes = 152
        self.path2 = nn.Sequential(
            nn.ConvTranspose3d(self.last_planes, self.last_planes, kernel_size = 2, stride = 2),
            nn.BatchNorm3d(self.last_planes),
            nn.ReLU(inplace = True))
        self.drop = nn.Dropout3d(p = 0.5, inplace = False)
        self.output = nn.Sequential(nn.Conv3d(248, 64, kernel_size = 1),
                                    nn.ReLU(),
                                    #nn.Dropout3d(p = 0.3),
                                   nn.Conv3d(64, 5 * len(config['anchors']), kernel_size = 1))
    def _make_layer(self, in_planes, out_planes, num_blocks, dense_depth, stride):
        strides = [stride] + [1]*(num_blocks-1)
        layers = []
        for i,stride in enumerate(strides):
            if debug: print(i, self.last_planes, in_planes, out_planes, dense_depth)
            layers.append(Bottleneck(self.last_planes, in_planes, out_planes, dense_depth, stride, i==0))
            self.last_planes = out_planes + (i+2) * dense_depth
        return nn.Sequential(*layers)

    def forward(self, x, coord):
        if debug: print('0', x.size(), 64)#, coord.size
        out0 = F.relu(self.bn1(self.conv1(x)))
        if debug: print ('1', out0.size())
        out1 = self.layer1(out0)
        if debug: print ('2', out1.size())
        out2 = self.layer2(out1)
        if debug: print ('3', out2.size())
        out3 = self.layer3(out2)
        if debug: print ('4', out3.size())
        out4 = self.layer4(out3)
        if debug: print ('5', out4.size())

        out5 = self.path1(out4)
        if debug: print ('6', out5.size(), torch.cat((out3, out5), 1).size())
        out6 = self.layer5(torch.cat((out3, out5), 1))
        if debug: print ('7', out6.size())
        out7 = self.path2(out6)
        if debug: print ('8', out7.size(), torch.cat((out2, out7), 1).size())#torch.cat((out2, out7, coord), 1).size()
        out8 = self.layer6(torch.cat((out2, out7, coord), 1))
        if debug: print ('9', out8.size())
        comb2 = self.drop(out8)
        out = self.output(comb2)
        if debug: print ('10', out.size())
        size = out.size()
        out = out.view(out.size(0), out.size(1), -1)
        if debug: print ('11', out.size())
        #out = out.transpose(1, 4).transpose(1, 2).transpose(2, 3).contiguous()
        out = out.transpose(1, 2).contiguous().view(size[0], size[2], size[3], size[4], len(config['anchors']), 5)
        if debug: print ('12', out.size())
        return out#, out_1

运行之后会在/results/dpn3d26/retrft960/路径下生成一系列.ckpt模型参数。

检测器测试

python main.py --model dpn3d26 -b 1 --resume results/dpn3d26/retrft96X/CkptFile/1000.ckpt --test 1 --save-dir dpn3d26/retrft96X/ --config config_training0.py

You should modify X from 0-9 to test 10 models for 10 fold cross validation mkdir in results/res18/retrft96X/val/#(X=0-9) mv results/res18/retrft96X/bbox/*.npy results/res18/retrft96X/val/

之后在val/路径下会有对subset9进行测试得到的肺结节。

具体流程如下：

Load np

Pad使得能被stride=4整除（便于corrd）

由side_len计算最多能取多少patch，nzhw = [nz,nh,nw]

Pad使得能被side_len整除

取patch，大小为 side_len + 2 * margin（144+2*32=208）（128+2*16=160）

返回patch的list

用缩小stride=4倍同样处理coord返回list

依次将patch输入网络，将结果拼回，

margin = 16，sidelen = 128，

Data shape ([18, 1, 160, 160, 160])

output.shape (18, 40, 40, 40, 3, 5)

You can use the ResNet or dual path net model by revising --model attribute.

在训练和测试之后, 可以使用 ./evaluationScript/frocwrtdetpepchluna16.py to validate the epoch used for test. After that, collect all the 10 folds' prediction, use ./evaluationScript/noduleCADEvaluationLUNA16.py to get the FROC for all 10 folds. You can directly run noduleCADEvaluationLUNA16.py, and get the performance in the paper.

下面是我自己的实验结果：