【论文】医学影像 | Elastic Boundary Projection for 3D Medical Image Segmentation 论文汇报 (CVPR 2019)

医学影像 | Elastic Boundary Projection for 3D Medical Image Segmentation 论文汇报 (CVPR 2019)

原创

联系邮箱：xiangchenzhang@jmu.edu.cn

论文地址：https://arxiv.org/abs/1812.00518

开源代码：https://github.com/twni2016/Elastic-Boundary-Projection

论文作者及研究团队：

Tianwei Ni, Lingxi Xie, Huangjie Zheng, Elliot K. Fishman, Alan L. Yuille

Peking University, Johns Hopkins University, Noah’s Ark Lab, Huawei Inc., Shanghai Jiao Tong University, Johns Hopkins Medical Institute

介绍：

医学成像分析（MedIA），特别是3D器官分割，是计算机辅助诊断（CAD）的重要先决条件。近年来，随着深度学习的蓬勃发展，卷积神经网络已被广泛应用于该领域，这大大提高了基于手工特征的传统分割方法的性能，甚至在许多器官和病灶的分割效果上超过了人类的精确度。

自然图像和医学图像之间最显着的差异之一在于数据维度：自然图像是平面的（2D），而诸如CT和MRI扫描的医学数据是立体的（3D）。

主流方法：

用于医学成像分割的现有深度神经网络可以分类为两种类型，在处理对象的维度上彼此不同。第一种类型将3D体积切割成2D/2.5D切片，并训练2D网络以单独处理每个切片或顺序处理。第二个是直接训练3D网络来处理体积数据。虽然后者可以更好的处理空间信息，但它有两个缺点：

缺乏预训练的模型使得训练过程不稳定，不易收敛。
大内存消耗使得难以接收整个数据作为输入，但将数据切割成patch会留下融合所有patch输出的最终阶段。

AH-Net是由周少华(曾在西门子负责研发与医学图像相关的创新产品，现任于中科院计算所研究员)于2018年发表在MICCAI上的论文《3D Anisotropic Hybrid Network: Transferring Convolutional Features from 2D Images to 3D Anisotropic Volumes》中提出的网络，AH-NEt试图利用迁移学习，把二维特征有效地转到三维。

数据准备

基于原始注释V ，我们需要定义一个新的矩阵C，其中的每一个C(x, y, z)储存着每个整数坐标(x, y, z)和边界∂V的有符号距离。

C的符号表示(x, y, z)是否落在边界内。（正:内，负:外，零:边界上）
∣ C ( x , y , z ) ∣ = min ⁡ ( x ′ , y ′ , z ′ ) ∈ ∂ V ( D i s t [ ( x , y , z ) , ( x ′ , y ′ , z ′ ) ] ) |C(x, y, z)|=\min_{(x',y',z')\in\partial V}(Dist[(x, y, z), (x', y', z')]) ∣C(x,y,z)∣=(x′,y′,z′)∈∂Vmin(Dist[(x,y,z),(x′,y′,z′)])
距离使用的为欧氏距离
l p − d i s t a n c e : D i s t [ ( x , y , z ) , ( x ′ , y ′ , z ′ ) ] = ( ∣ x − x ′ ∣ 2 + ∣ y − y ′ ∣ 2 + ∣ z − z ′ ∣ 2 ) − 1 2 l_p-distance: Dist[(x, y, z), (x', y', z')] = (|x-x'|^2+|y-y'|^2+|z-z'|^2)^{-\frac12} lp−distance:Dist[(x,y,z),(x′,y′,z′)]=(∣x−x′∣2+∣y−y′∣2+∣z−z′∣2)−21
应用KD树来处理数据，可以极大的加速整体计算速度。

整体结构

论文提出了一种弥合2D网络和3D分割之间差距的新方法。它使用球面坐标系将不规则边界投影到矩形中，通过找到2D表面的边界来完成分割。这个方法称为弹性边界投影（EBP），可以将其应用于2D网络。

EBP先初始化一些位于目标器官内外的枢轴点和围绕它的外壳。外壳由半径沿不同方向被初始化为一个完美的球体训练过程中不断调整壳体，使其最终收敛到目标器官的边界，为此我们训练2D网络以预测每个终点是否位于器官内部或外部，并相应地增加或减少该方向的半径。

B = ∂ V B=\partial V B=∂V
在ROI区域中随机抽样pivot，并建立球坐标系，d是以p为起点的固定方向。

首先确定D的方向：
D = { d 1 , d 2 , . . . , d N } D=\{d_1, d_2,..., d_N\} D={d1,d2,...,dN}
我们将D表示为Ma 元素的方位角集和直角集Mp 的笛卡尔积，其中:
M a × M p = M M^a\times M^p=M Ma×Mp=M
Ma方位角是均匀分布的，即：
α m 1 = 2 m 1 π / M a α m 1 ∈ [ 0 , 2 π ) \alpha_{m_1}=2m_1\pi/M^a\\ \alpha_{m_1}\in [0, 2\pi) αm1=2m1π/Maαm1∈[0,2π)
Mp极角在赤道附近具有更密集的分布，即：
φ m 2 = cos ⁡ − 1 ( 2 m 2 M p + 1 − 1 ) − π 2 φ m 2 ∈ [ − π 2 , π 2 ] \varphi_{m_2}=\cos^{-1}(\frac{2m_2}{M^p+1}-1)-\frac\pi2\\ \varphi_{m_2}\in[-\frac\pi2, \frac\pi2] φm2=cos−1(Mp+12m2−1)−2πφm2∈[−2π,2π]
M个单位矢量近似均匀地分布在球体上，单位矢量(x, y, z)满足：
x ^ = cos ⁡ ( α m 1 ⋅ cos ⁡ ( φ 2 ) ) y ^ = sin ⁡ ( α m 1 ⋅ cos ⁡ ( φ 2 ) ) z ^ = sin ⁡ ( φ 2 ) \hat{x}=\cos(\alpha_{m_1}\cdot\cos(\varphi_2))\\ \hat{y}=\sin(\alpha_{m_1}\cdot\cos(\varphi_2))\\ \hat{z}=\sin(\varphi_2) x^=cos(αm1⋅cos(φ2))y^=sin(αm1⋅cos(φ2))z^=sin(φ2)
网络的输入与输出：

输入：（256，120，120）
I = U ( p n , { d m , r n , m } ) = U ( e n , m ) I=U(p_n, \{d_m, r_{n,m}\})=U(e_{n,m}) I=U(pn,{dm,rn,m})=U(en,m)
训练：（256，120，120）
M : O = f ( I ; θ ) M:O=f(I;\theta ) M:O=f(I;θ)

训练：数据生成和优化

然后，我们在ROI区域内采样获取一组枢轴P：
P = p 1 , p 2 , . . . , p N P={p_1, p_2, ..., p_N} P=p1,p2,...,pN
在每个P_n中，我们构造一个半径为R_0的单位球体，t次迭代后，每个终点的坐标通过以下公式计算。
e n , m ( t ) = p n + r n , m ( t ) ⋅ d m e^{(t)}_{n,m}=p_n+r^{(t)}_{n,m}\cdot d_m en,m(t)=pn+rn,m(t)⋅dm
在数据标注中找到所有的边界坐标对应的输入-输出数据对：
I n , m ( t ) = U ( { e n , m ( t ) } ) O n , m ( t ) = C ( { e n , m ( t ) } ) I^{(t)}_{n,m}=U(\{e^{(t)}_{n,m}\})\\ O^{(t)}_{n,m}=C(\{e^{(t)}_{n,m}\})\\ In,m(t)=U({en,m(t)})On,m(t)=C({en,m(t)})
迭代直至收敛，收敛后所有终点都落在边界上或坍缩到p自身。
r n , m ( t + 1 ) = max ⁡ { r n , m ( t ) + C ( e n , m ( t ) ) , 0 } r^{(t+1)}_{n,m}=\max\{r^{(t)}_{n,m}+C(e^{(t)}_{n,m}), 0\} rn,m(t+1)=max{rn,m(t)+C(en,m(t)),0}
有时候射线可能与边界具有多个交叉点，在这种情况下，算法会收敛到最接近初始球体的模型。作者没有在训练和测试中特别处理这个问题，情况建立在以下两个假设之上：
1. 大多数的终点都接近边界
2. 枢轴足够密集

构建完成训练集之后，采用常规方法优化，例如，使用随机梯度下降。

测试：迭代和推理

从一组随机放置的枢轴开始，围绕每个枢轴建立一个单位球体，对参数θ迭代直到收敛或大导最大轮数。（测试过程中，迭代不一定会收敛。）收集除了坍缩到枢轴上所有枢轴的所有终点，进行三维重建。

论文中提到了使用了四种方法，来提高测试精度。

第一，输入图像，只包含了一层边界信息，这样的信息太过稀疏单一，通常不足以精确预测，因此，采用增加更多的通道数来进行补充。有以下La,Lb两种通道:
s n , l A , m ( t ) = r n , m ( A ) + l A − ( L A + 1 ) / 2 s n , l B , m ( t ) = l B L B + 1 [ r n , m ( t ) − ( L A + 1 ) / 2 ] s^{(t)}_{n,l^A,m}=r^{(A)}_{n,m}+l^A-(L^A+1)/2\\ s^{(t)}_{n,l^B,m}=\frac{l^B}{L^B+1}[r^{(t)}_{n,m}-(L^A+1)/2] sn,lA,m(t)=rn,m(A)+lA−(LA+1)/2sn,lB,m(t)=LB+1lB[rn,m(t)−(LA+1)/2]
当LA=1，LB=0,生成的是单层边界。论文使用相对较大的LA=LB=5, 这样可以获得更多的上下文信息，但是网络仍然使用2D网络。无论LA，LB 如何，得到的O中的通道数均保持为1。

其次，我们利用分布的空间的一致性(spatial consistency)来提高准确性，当提供迭代rt的半径值时，我们可以随机抽样M个ε~N(0, σ²)，σ的值很小，将他们作为噪音加入到r。

通过空间一致性，每个方向m总是满足以下近似值：
C ( e n , m ( t ) + ϵ m ⋅ d m ) = C ( e n , m ( t ) ) + ϵ m ⋅ cos ⁡ β ( d m , e n , m ( t ) ) β ( d m , e n , m ( t ) ) 是 d m 和法向量之间的角度 C(e^{(t)}_{n,m}+\epsilon_m\cdot d_m)=C(e^{(t)}_{n,m})+\epsilon_m\cdot \cos\beta(d_m,e^{(t)}_{n,m})\\ \beta(d_m,e^{(t)}_{n,m})是d_m和法向量之间的角度 C(en,m(t)+ϵm⋅dm)=C(en,m(t))+ϵm⋅cosβ(dm,en,m(t))β(dm,en,m(t))是dm和法向量之间的角度
尽管该角度通常难以i算，但是我们可以将等式6的左侧作为εm的线性函数，并使用εm的多个样本将其值估计为0。（注：这句我不是很理解。）

第三，缩小了培训和测试数据分布之间的差距。

请注意！

在训练阶段，所有r_(n, m)^((t) )的值都是使用标注生成的。

而在测试阶段，它们是通过网络预测累积的。

因此，如果模型从未在这样的无标注数据上训练过，他们可能会随着迭代积累误差。

最后，发现测试阶段的大多数假阳性都是由外部枢轴产生的，尤其是具有相似物理特性的另一个器官内的枢轴。在这种情况下，边界可能汇聚到其他的边界上，影响分割效果。

为了解决这个问题，作者还引入了一个额外的阶段来确定那些枢轴位于目标器官内。随机抽取感兴趣区域（ROI）中的几千个点并计算相邻两个壳体的IOU。

然后可以将整个图形分成两部分，内部部分定义为内部枢轴的集合，在3D重建中仅考虑由内部枢轴产生的边界点。

三维重建

最后一步是收集所有终点并重建三维体积的表面，由于模型并不完美，因此有很多预测的终点不会落在实际边界上，采用核密度估计（KDE）来消除他们.

使用带宽为1的Epanechnikov内核，并保留所有整数坐标，其对数释然性不小于-14

最后，应用一个基本的图形框架来实现这个目标，其工作原理如下。

首先使用Delaunay三角剖分在剩下的终点上构建网格结构，然后去除周长大于α的不正确四面体。在我们获得alpha形状之后，我们使用细分算法将其体素化。这保证了最终获得的分割有封闭的边界。

绿色和蓝色的点分别表示ground-truth和预测的内部枢轴。

分析

创新点：

EBP的核心贡献是为3D分割提供了基于2D的方法。目前，这个想法还是第一次在深度学习文献中进行过研究。类似的传统的分割方法，如，GraphCut 和GrabCut 转换了2D分割，一维最小化目标函数的轮廓。
局限性：

有时候，即使是专业的医生也无法找到准确的边界位置，这种不精确性会在很大程度上影响边界体素采样的效果。这也解释了之前提到的该公式提供更多内部信息公式的有效性。
s n , l A , m ( t ) = r n , m ( A ) + l A − ( L A + 1 ) / 2 s n , l B , m ( t ) = l B L B + 1 [ r n , m ( t ) − ( L A + 1 ) / 2 ] s^{(t)}_{n,l^A,m}=r^{(A)}_{n,m}+l^A-(L^A+1)/2\\ s^{(t)}_{n,l^B,m}=\frac{l^B}{L^B+1}[r^{(t)}_{n,m}-(L^A+1)/2] sn,lA,m(t)=rn,m(A)+lA−(LA+1)/2sn,lB,m(t)=LB+1lB[rn,m(t)−(LA+1)/2]
其次，由于其的工作机制，EBP难以处理具有严重不规则几何形状的3D目标，例如，微小的血管。

实验

VNet是一个基于3D的网络，从原始补丁中随机裁减成128×128×64的patch进行训练，以便不超过GPU内存限制（12GB）。在测试阶段，相同大小的滑动窗口沿三个轴有规律地移动，并且每个体素在其获得的所有预测上取平均值。

RSTN是一个基于2D的网络，它使用3个相邻切边，并将三个（冠状，矢状，和轴向）视图融合到最终预测中。

分析表明，是RSTN在几个连续的2D切片中未能检测到该器官的一部分，但是通过检测边界，EBP成功的恢复了这种情况。这表明，在许多情况下，EBP和RSTN可以为器官分割提供补充信息。

MSD脾脏数据集具有相对较低的空间分辨率，（沿长轴31-168个体素），这给VNet 带来了额外的困难。

当人工标注出现错误的情况：
可以明显看出来标注没有贴合边界

医学影像中的标注通常是不完美的。
DSC不能准确的反应分割的绝对质量，更高的DSC不能保证分割效果就更好

EBP如何找到边界

收敛性:

图中的15条细曲线表示十五个枢轴点，粗曲线为它们的平均值。

从图中可以看到大部分枢轴，都会在第25次迭代时汇聚到边界。在实践中，为了提升速度，只进行10次迭代。虽然没有实现收敛，并且许多终点不落在边界。但是3D重建算法可以过滤这些异常值，而不影响重建效果。

一致性：

横坐标为枢轴沿着固定方向前进直到它离开边界时记录的轴间DSC。

可以观察到所有的曲线都会在某个地方突然下降，这通常表示枢轴从内部到外部的时刻。

内-内枢轴平均 DSC(73.46%) 内-外枢轴平均 DSC(51.37%)

实验表明，如果两个相邻的枢轴都位于目标内，则两个枢轴的相似度会很高。

结论

本文介绍了EBP，这是一种训练2D深度网络进行三维物体分割的新方法。核心思想是建立一个弹性壳并对其进行调整，直到它收敛到目标的实际边界。

可以应用2D网络（低计算开销，更少参数，预训练模型等）来处理三维信息（更丰富的上下文信息）。在腹部CT扫描中对几个器官进行实验，并且EBP实现可以与2D网络和3D网络相竞争的性能。

从这项工作中了解到，高维数据通常包含冗余（例如，并非每个体素在3D体积中都有用），并且挖掘判别部分虽然具有挑战性，但通常会导致更高效的模型。

应用

目标内，则两个枢轴的相似度会很高。

结论

本文介绍了EBP，这是一种训练2D深度网络进行三维物体分割的新方法。核心思想是建立一个弹性壳并对其进行调整，直到它收敛到目标的实际边界。

应用

该方法或许可以应用到修复还原不完整的三维点云数据上，应用前景广泛。但是计算时间偏长，数据处理占用内存大。（论文使用的数据集将1.5G的原始CT数据信息制作成了900多G的标注信息。）