3D Deep Leaky Noisy-or Network 论文阅读
原文:Evaluate the Malignancy of Pulmonary Nodules Using the 3D Deep Leaky Noisy-or Network
博文参考:Doublle Tree的博客中Evaluate the Malignancy of Pulmonary Nodules Using the 3D Deep Leaky Noisy-or Network 论文阅读一文。
注:本文为2017年Kaggle举办的数据科学竞赛中,第一名获奖团队的相关论文,若需查看代码可访问Github。
简介
根据CT图像的肺癌自动诊断系统包含以下步骤:
- 检测所有可疑病变;
- 评估整个肺部的恶性程度。
但目前大多数的研究主要集中于第一步,以及通过肺结节诊断肺癌存在较高的假阳性率。因此,肺癌的诊断需要对每个可疑结节进行细致分析,再联合所有结节信息进行定性诊断。针对上述问题,本文提出了一个三维深度神经网络(3D deep neural network)用于解决这些问题。该网络由两部分组成:
- 用于结节检测的3D region proposal network;
- 基于置信检测(the detection confidence)选出top-5结节并评估其癌症可能性,最后将此概率与Leaky noisy-or模型相结合评估患者患癌的可能性。
其中,上述两个模型均采用修改后的U-net模型,并使用数据增强操作避免过拟合问题。
数据集和预处理
数据集
训练集由LUNA16数据集(the Lung Nodule Analysis 2016)和NDSB3(Data Science Bowl 2017)数据集两部分组成。其中,LUNA16数据集含有888个病例,标记了1186个肺结节;在NDSB3数据集中,1397个病例用于训练,198个病例用于验证,506个病例用于测试,且人工标注了训练集中754个结节和验证集中78个结节。
对于LUNA16数据集,其存在许多较小的注释结节,且临床经验认为直径6mm以下的肺结节无危险。但在NDSB3数据集中,存在较多的大直径结节且结节多与主支气管相连。因此,针对两个数据集的差异,需去除LUNA16数据集中直径6mm的结节,同时对NDSB3数据集进行人工标注。
此处说明了Julian de Wit的解决方案中,直接设置结节直径为6mm的原因,以及为何需对NDSB3数据集进行人工标注。
上图为结节分布情况图。其中,图a为DSB(NDSB3)与LUNA(LUNA16)数据集中结节直径分布情况;图b为DSB数据集中患癌患者与健康人群的最大结节直径分布情况。
预处理
首先将所有的原始数据转变为HU值,如下图a所示,再进行如下步骤:
-
掩膜提取:在2D切片上,首先使用标准差为1的高斯滤波和阈值为-600的处理得到肺部以及周围较暗部分的掩膜,如上图b所示,然后进行连通性分析去除小于 30 m m 2 30mm^2 30mm2的connected component和离心率大于0.99的部分(some high-luminance radial imaging noise),再计算得到二值的3D矩阵中所有的3D connected component,且仅保留非边缘部分(用于去除肺部周围较暗的部分)以及体积在0.68~7.5L之间的部分,结果如上图c所示;
-
凸包与扩张:若结节与肺的外壁相连,则其将不会出现在上述提取的掩膜中。因此,对于这种情况,首先将肺部分为左右两个部分,即左肺与右肺,如上图d所示。然后分别对左右肺进行凸包处理,并向外扩张10像素,如上图f所示。但对于一些2D切片而言,肺部的底部类似与月牙形,如下图所示。若对于该类型进行凸包处理后,面积大于初始的1.5倍,则放弃凸包,从而避免引入过多的其他组织;
-
灰度标准化:将HU值([-1200, 600])线性变换至0~255内的灰度值,且掩膜以外的像素灰度值均设为170,以及扩张区域内的像素灰度值高于210则也设为170。
用于结节检测的3D卷积神经网络
该网络是基于U-net的3D版RPN(Region Proposal Network)模型。
输入数据
受限于显存,输入数据大小为 128 × 128 × 128 × 1 ( H e i g h t × L e n g t h × W i d t h × C h a n n e l ) 128×128×128×1(Height×Length×Width×Channel) 128×128×128×1(Height×Length×Width×Channel),并随机选择两种patch:一种为70%的输入数据至少包含一个结节;另一种为30%的输入数据不含结节。其中,patch超出图像部分用灰度值为170填充。
为了避免过拟合问题,数据采用数据增强方法。
此处可根据实际情况将输入数据大小调整为64或者32。
网络结构
网络由前馈路径和反馈路径组成,如下图图a所示。
前馈路径
以两层卷积核为 3 × 3 × 3 3×3×3 3×3×3的卷积(channel为24)开始,且padding为1;其后为4个残差块,其中每个残差块由3个残差单元组成(如上图图b所示),而每个残差单元由卷积、Batch Norm、ReLU激活函数、卷积和Batch Norm组成,且卷积核大小均为 3 × 3 × 3 3×3×3 3×3×3。除此之外,每个残差块均有一个最大池化层,大小为 2 × 2 × 2 2×2×2 2×2×2,步长为2。
反馈路径
反馈路径由两层反卷积(装置卷积)层和两个融合单元构成。最后,由卷积核均为 1 × 1 × 1 1×1×1 1×1×1且channel为64和15的两层卷积层将数据大小转换为 32 × 32 × 32 × 15 32×32×32×15 32×32×32×15。
反卷积层
卷积核大小为2,步长为2。
注意该部分的代码实现部分,原始U-net网络设置为不可学习。
融合单元
每个融合单元(如上图图c所示)均由一个前馈blob和反馈blob组成,其结果作为残差块的输入。
值得注意的一点,本文作者在此处引入了位置信息,作为额外的输入数据。
位置信息
proposal的位置信息可能影响是否为结节和是否为恶性的判断,因而引入位置信息。
具体方法:对于每个patch,计算其相对位置坐标,并将其大小转换为 32 × 32 × 32 × 3 32×32×32×3 32×32×32×3。
其中,位置坐标对应归一化后的X,Y和Z轴(每个轴的取值范围为-1~1,对应于肺的两端)。
输出层
输出数据为4D的tensor, 32 × 32 × 32 × 3 × 5 32×32×32×3×5 32×32×32×3×5,其中3表示anchor个数,5表示回归量( o ^ , d x ^ , d y ^ , d z ^ , d r ^ \hat{o}, \hat{d_x}, \hat{d_y}, \hat{d_z}, \hat{d_r} o^,dx^,dy^,dz^,dr^,即概率,三维坐标和bounding box直径大小)。
其中,对于 o ^ \hat{o} o^激活函数采用sigmoid函数,其余不使用任何激活函数。
损失函数
真值标签为 ( G x , G y , G z , G r ) (G_x, G_y, G_z, G_r) (Gx,Gy,Gz,Gr),每个anchor记为 ( A x , A y , A z , A r ) (A_x, A_y, A_z, A_r) (Ax,Ay,Az,Ar),IoU(Intersection over Unit)大于0.5记为Positive,小于0.02记为False,其他在训练过程中忽略。
分类损失为:
L c l s = p l o g ( p ^ ) + ( 1 − p ) l o g ( 1 − p ^ ) L_{cls}=plog(\hat{p}) + (1-p)log(1-\hat{p}) Lcls=plog(p^)+(1−p)log(1−p^)
其中, p p p为anchor box的真值标签, p ^ \hat{p} p^为预测概率。
bounding box回归标签为:
d x = ( G x − A x ) A r d y = ( G y − A y ) A r d z = ( G z − A z ) A r d r = l o g ( G r A r ) d_x = \frac{(G_x - A_x)}{A_r} \\ d_y = \frac{(G_y - A_y)}{A_r} \\ d_z = \frac{(G_z - A_z)}{A_r} \\ d_r = log(\frac{G_r}{A_r}) dx=Ar(Gx−Ax)dy=Ar(Gy−Ay)dz=Ar(Gz−Az)dr=log(ArGr)
回归总损失为:
L r e g = ∑ k ∈ { x ,   y ,   z ,   r } S ( d k ,    d k ^ ) L_{reg}=\sum\limits_{k\in\{x,\,y,\,z,\,r\}}S(d_k,\; \hat{d_k}) Lreg=k∈{x,y,z,r}∑S(dk,dk^)
其中, S S S为smoothed L1-norm函数:
S ( d , d ^ ) = { ∣ d − d ^ ∣ ,      i f    ∣ d k − d ^ ∣ > 1 ( d − d ^ ) 2 , e l s e S(d, \hat{d})=\begin{cases} |d-\hat{d}|, \;\;\quad if\; |d_k - \hat{d}|>1\\ (d-\hat{d})^2, \quad else \end{cases} S(d,d^)={∣d−d^∣,if∣dk−d^∣>1(d−d^)2,else
对于每个anchor box的损失函数为:
L = L c l s + p L r e g L = L_{cls} + pL_{reg} L=Lcls+pLreg
最后,整体的anchor box的损失函数为anchor box的损失值取平均。
正反例数据
正例数据
对于大结节而言,网络会生成较多的positive anchor box,因此为了降低训练数据之间的相关性,随机挑选其中一个。
由于结节直径大小分布不均,而NDSB3数据集多为大结节,因而对大于30mm和40mm的结节,采样频率分别是其他结节的2倍和6倍。
反例数据
对于一些易误诊为结节的反例数据,通过使用hard negative mining方法解决。
具体方法为:
- 将不同的patch输入至网络得到不同置信度的输出映射;
- 随机选择N个反例数据构成候选池;
- 侯选池中的数据以置信度值大小排序,且选出top-n的数据作为反例数据。
未选中的数据忽略且不参与损失计算。
图像分割(测试过程)
输入数据大小为 208 × 208 × 208 × 1 208×208×208×1 208×208×208×1,overlap为32像素。
输出数据为 { x i , y i , z i , r i , p i } \{x_i, y_i, z_i, r_i, p_i\} {xi,yi,zi,ri,pi},其中 x i , y i , z i x_i, y_i, z_i xi,yi,zi表示proposal中心坐标, r i r_i ri表示其半径大小, p i p_i pi表示其置信度。
输出数据且使用非极大值抑制操作来去除overlaping proposal。
肿瘤分类
由于受限于训练样本数,因而复用结节检测器阶段的N-Net网络。
输入数据为结节的proposal,大小均为 96 × 96 × 96 × 1 96×96×96×1 96×96×96×1,其仅使用了结节中心点的信息。在分类器训练阶段,随机挑选proposal,且其选中的概率与proposal的置信度成正比;在测试阶段,只挑选top-5的proposal。
经卷积核为 24 × 24 × 24 × 128 24×24×24×128 24×24×24×128的最后一个卷积层得到输出结果;随后提取每个proposal中心处 2 × 2 × 2 2×2×2 2×2×2的体素,并将其通过最大池化操作后得到128维的特征,如下图图a所示。
对比四种预测肿瘤类别的方法(Feature combining method,MaxP method,Noisy-or method和Leaky Noisy-or method),挑选出Leaky Noisy-or方法作为最终的分类方法,如下图图b所示。
Leaky Noisy-or Method
引入一个假想结节,其患癌概率为 P d P_d Pd, P d P_d Pd的值在模型训练阶段学习获得。
将特征输入至两层相同的Perceptron得到分类概率 P P P:
P = 1 − ( 1 − P d ) ∏ i ( 1 − P i ) P=1-(1-P_d)\prod\limits_i(1-P_i) P=1−(1−Pd)i∏(1−Pi)
其中, P i P_i Pi表示第i个结节癌变的概率。
训练过程
损失函数为交叉熵函数。为了避免过拟合采用了数据增强和正则化操作。
训练的步骤:
- transfer检测器训练参数后,再训练分类器;
- 采用gradient clipping方法训练分类器,随后存储BN(Batch Normalization)参数;
- 用存储的BN参数和gradient clipping方法交替训练检测器和分类器。
注:BN在训练阶段和测试阶段所计算的方法有所差异。因复用N-Net网络,分类器和检测器交替训练,因而需对BN的参数做特殊处理。