在前面进行了肺结节数据的预处理之后,接下来开始进入肺结节检测环节。首先附上该项目的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编程和深度学习 微信公众号,前往查看
如有侵权,请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!