在前面进行了肺结节数据的预处理之后,接下来开始进入肺结节检测环节。首先附上该项目的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.
下面是我自己的实验结果:
本文分享自微信公众号 - Python编程和深度学习(Python_Deeplearning),作者:Python编程和深度学习
原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。
原始发表时间:2020-01-01
本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。
我来说两句