基于三维曲面重建的CT图像肺边界缺陷修复

张欣，徐永潇，王兵

(1.河北大学 电子信息工程学院，河北 保定 071002；2. 河北大学 数学与信息科学学院，河北 保定 071002)

肺癌目前已经是致死率非常高的恶性肿瘤，而且肺癌的发病率也在不断上升[1].计算机断层扫描技术(computed tomography,CT)和计算机辅助诊断(computer aided diagnosis,CAD)为肺癌的诊断提供了更多的方法[2].为了提高医生对肿瘤诊断的精确性和可靠性，需要对肺及肺部病变部位进行准确的分析.有些CT图像中肿瘤粘连着肺边界，其密度与胸廓或心脏等部位的密度相差很小,并且肿瘤的形状、大小和粘连的位置都存在很大的差异[3]，导致对肺实质与肿瘤的分割没有一个全面有效的方法.

针对肺实质的分割与边界修复问题，多采用阈值法[4]、聚类法[5]、边缘检测法和区域生长[6]等方法对CT肺图像进行粗分割，然后对肺边界粘连型肿瘤造成的肺边界缺陷部分进行修复.目前的修复方法有滚球法、链码法[7]、凸包算法[8]等.文献[5]提出了一种基于灰度积分投影与模糊C均值聚类并结合滚球法修复缺陷部分的分割方法，但滚球半径的选择对肺边界修复效果影响较大.文献[7]提出了一种基于改进链码和Bresenham算法结合的缺陷肺实质修复方法，但该算法容易出现过分割现象.凸包算法的修复效果受到缺陷处形状的影响较大，当需要修复的肺实质边界是曲率变化较小的平滑曲线时，凸包算法可以作为首要的修复方法.代双凤等[8]提出一种基于3D区域增长和凸包算法相结合的分割方法，但是针对缺陷位于二维肺实质边界曲率变化比较大处的情况，例如当缺陷位于肺叶尖端时，二维肺部CT图像中缺少指导肺边界缺陷修复的信息，效果不理想.因此，本文提出一种基于三维曲面重建[9]的方法，实验证明该方法用于修复肺实质边界曲率变化较大处的缺陷效果良好，修复后的肺实质有助于准确分割CT图像中的肺肿瘤.

1 算法描述与分析

基于三维曲面重建算法的肺边界缺陷修复方法主要由3部分组成：CT图像提取肺实质、肺实质的三维重建、修复具有缺陷的三维肺实质，具体流程如图1所示.

图1 肺实质修复流程Fig.1 Flow chart of pulmonary parenchyma restore

1.1 CT图像提取肺实质

从CT图像中提取出肺实质区域：首先将肺部CT图像进行格式和灰度转换，再进行高斯滤波去除噪声，得到较为平滑的人体肺部CT图像；然后为了更加快速地提取序列CT图像中的肺实质区域，使用质心灰度法改进区域生长算法；最后通过形态学操作填补肺结节和气管等造成的空洞，得到肺实质区域.

由于传统的区域生长算法对每一张二维图像都需要人为选取种子点，操作起来较为繁琐.本文使用相邻层生长区域的质心位置结合邻域灰度阈值筛选的方法自动确定当前层切片图像的生长种子点，适用于从肺部序列CT图像中提取肺实质.具体步骤如下.

步骤1：在某层肺实质区域内选取1个或多个种子点，标记为已生长区域R[t]，令t为迭代次数.

步骤2：设平均灰度值即为生长区域的相似性质，计算已生长区域R[t]的平均灰度值

(1)

其中amount(R[t])表示R[t]内像素点的数量，I(y)表示像素y点的灰度值.

步骤5：在相邻切片上选择qc+1的8邻域像素点中满足min|I(Qc+1)-μc|条件的像素点qo,其中Qc+1为qc+1的8邻域像素点集，I(Qc+1)表示点集Qc+1中像素点的灰度值.将qo作为该层切片的初始种子点，重复步骤2～5，直至到达最底层切片中肺部消失.此时该切片上没有肺实质，中心点的灰度值与相邻切片的平均灰度值相差较大，设置合理的终止阈值φ，当无法满足条件I(qo)<φ时停止生长，即得到完整序列的二维肺实质图像.

1.2 肺实质的三维重建

因为后续缺陷修复需要三维曲面信息，故先采用面绘制中最具有代表性的MC(marching cubes)算法[10]进行肺实质的三维重建，即通过阈值判断提取三维模型表面的灰度等值面实现重建.再使用拉普拉斯平滑算法对生成的三维模型进行平滑处理，最终得到表面光滑的三维肺实质.具体步骤如下.

步骤1：读取经过预处理的二维肺实质图像.在三维数据场中扫描所有空间像素点并构建体素(体素是由相邻的2个切片垂直对应的8个像素点组成的立方体).由于灰度等值面与不同体素相交后，体素内部会形成不同的剖分面，将立方体中所有可能产生的剖分面进行分类，创建立方体构型的索引表用于响应体素.

步骤2：提取表面体素.与等值面相交的体素为表面体素，即体素棱边与等值面的交点大于等于1.通过线性插值法求出上述交点即为等值面片的顶点，计算出等值面片顶点的法向量用以提供三维模型的表面光照变化信息，使得视觉效果更好.

步骤3：将等值面片进行连接，选定光照模型并进行平滑处理，完成三维肺实质的重建.

1.3 修复肺实质缺陷

根据三维曲面信息修复肺实质边界缺陷的方法分为：点云数据集构建和Possion三维曲面重建.

1.3.1 构建点云数据集

在肺实质缺陷处周围的表面进行采样，将所得点集记为P.具体步骤如下.

步骤1：采用凸包理论中常用的Graham扫描法和阈值判断法相结合的方法[8]找出序列CT图像中肺边界的缺陷部位.然后将缺陷部位扩大τ个像素的距离并采集二维缺陷的边界信息，将所得信息映射到重建的三维肺实质上，即可得到肺实质缺陷处的三维轮廓点集.对三维点集进行插值和连接，最终得到肺实质缺陷处的三维轮廓曲线.

步骤2：计算三维轮廓点集中质心与最远点的距离r.构建一个以质心为球心、2r为半径的球体，对球体区域与缺陷三维肺实质区域进行布尔运算中的“交运算”，获得感兴趣区域ROI(region of interest).针对ROI区域进行操作可以提高运算效率.

步骤3：采集ROI区域的点云数据，删除位于ROI区域的缺陷轮廓曲线内点云，得到点云数据集P.

1.3.2 Possion三维曲面重建算法修复肺实质

Possion曲面重建算法是一种隐式曲面重建算法[11]，具有局部和全局拟合的特点，鲁棒性强，可以使重建出来的肺实质三维轮廓更精准平滑.该算法的核心就是构建空间元素的向量场来近似模拟三维肺实质表面的指示函数[9]，将向量场的求解转化为求解Possion等式的问题，通过求解的Possion等式得到曲面重建所需要的等值面.具体步骤如下.

步骤1：节点混合函数的定义.首先对P建立一个节点为o的八叉树ϑ.然后令基函数F∶R3→R表示函数空间，该基函数具有节点积分和尺度缩放性质.所得函数空间具有多尺度结构[9].设q为函数空间中的变量，为函数空间中的ϑ的每个节点o添加一个对应于单位积分的混合函数Fo，将混合函数以节点o的中心和宽度进行展开可得

(2)

其中o.c和o.w为节点的中心和宽度.所有节点的混合函数Fo线性求和可以表示向量场V.

步骤2：基函数的定义.采用盒式滤波器的n维卷积近似基函数F，该方法能够使节点函数Fo(q)的线性求和更加精确的表示向量场V，基函数

(3)

步骤3：向量场的构造.先针对当前节点临近的8个节点采用三次线性插值的方法进行分配，然后再构造梯度域函数的向量场

(4)

对点云数据P进行平均采样，所得样本数据s∈P，其中，样本数据中的点s.p∈s，s.N为s的表面内法向量，NgbrD(s)表示深度为D的八叉树中最靠近s.p的节点区域，αo,s为权重系数.因为三维模型表面都存在一个几乎不会发生变化的指示函数，可以看做一个常量.因此指示函数的梯度在任何除了表面的节点上都为0，即可以用上面构造的向量场近似估计三维肺实质ROI区域表面.

(5)

步骤5：等值面的提取.通过均值法计算等值，再从指示函数中提取等值面，使用MC算法对等值面进行曲面重建即可将三维肺实质ROI区域的表面修复完整.对修复面与带有缺陷的肺实质进行“并运算”即可将肺实质修复完好.最后将切片轮廓从修复完好的肺实质中提取出来，进行填充处理后制成掩模，该掩膜用于从原始CT图像中分割肺实质.

2 实验分析

实验数据共选取了10组肺部带有粘连型肿瘤的三维CT图像数据.一部分数据来源于河北大学附属医院CT室，数据格式为DICOM，CT图像扫描厚度为3 mm，分辨率为512×512像素. 另一部分来源于美国国家癌症研究中心公开的DSB2017肺癌预测数据集(data science bowl)，CT图像的格式为DICOM，扫描厚度为2 mm，分辨率为512×512像素.实验硬件平台基于Windows10操作系统，Intel(R) Core(TM) i5处理器，8 G内存， 2.50 GHz主频，软件平台基于Matlab R2016a，VTK，Visual Studio 2017.

2.1 应用本文方法修复肺实质缺陷

基于Possion三维曲面重建算法的肺实质修复方法实验过程效果如图2所示.首先选取带有边界粘连型肿瘤的肺部CT图像，缺陷位于左肺上半部分的尖端处标注的地方，如图2a所示.使用改进的区域生长法提取肺部序列图像中的肺实质图像，经过多次实验对比后设置生长阈值θ=25,φ=180.再使用形态学操作填补肺实质中由肺结节和血管等造成的孔洞，得到肺实质的掩模图像如图2b所示.通过掩模得到肺实质的初步分割结果图如2c所示，从该结果可以明显看到缺陷部位.然后使用MC算法对所有切片初分割后的肺实质进行三维重建，如图2d所示.提取二维肺实质缺陷边界并映射到三维空间，经过多次实验对比后设置τ=3.构建球体后提取ROI区域，同时标记缺陷轮廓边界，如图2e所示.再对肺实质的缺陷部位进行重采样获得点云，如图2f所示.使用本文的Possion曲面重建算法对缺陷部位进行三维曲面重建，依据文献[9]中的经验值设置八叉树的最大深度D=8，如图2g所示.将重建后的缺陷部位与原始缺陷肺实质进行结合，结果如图2h所示.最后从修复完整的肺实质中提取二维切片，如图2i所示.对二维切面进行填充获得掩模后再对CT图像中肺实质进行分割，得到的最终的肺实质分割结果如图2j所示.从分割结果可以看出，采用本文的修复方法可以修复标注的肺实质缺陷，且修复的边缘更加平滑，同时将肿瘤完整的包含在了肺实质中.

a.原始肺部CT图像；b.区域生长得到掩摸；c.通过掩摸提取肺实质；d.三维重建肺实质；e.提取缺陷部位；f.重采样获得点云；g.三维重建缺陷部位；h.修复完整的肺实质；i.提取切片处轮廓；j.最终分割结果.图2 本文方法修复肺实质的实验效果Fig.2 Experimental effect of restoring lung parenchyma with the method used in this paper

2.2 本文的修复方法与经典的“滚球法”和“凸包算法”相比较

在肺实质的修复方法中，最经典的2种方法分别是“滚球法”和“凸包算法”.它们具有计算速度快、原理简单易操作的优点，但是在修复效果上都存在着一定的缺陷.将本文所使用的修复方法与滚球法和凸包算法应用到同一张切片上比较，修复效果如图3所示.选用的切片原图为图3a，从图3b和图3c可以看出使用滚球法修复肺实质的效果很依赖滚球半径.当滚球半径选的较小时，如r=5，会出现肺实质的“欠分割”，如图3b中圈起来的部分所示，缺陷并没有得到修复.当滚球半径选的较大时，如r=15，从图3c中上圆标注的位置可以看出缺陷处得到一定程度的修复，但是在肺门处却出现了“过分割”.从图3d中可以看出使用凸包算法修复缺陷肺实质时有一定的局限性.对肺实质边界曲率变化较大处的缺陷，由于二维图像中可用于修复缺陷的参考特征信息较少，所以凸包算法的修复效果差.此外，从标记的地方可以看出在其他正常边界处容易出现“过分割”的现象.而本文采取的修复方法可以将粘连型肿瘤导致的肺实质缺失部分补充出来，修复后的肺实质可以完整的包含肿瘤，且边缘连续性和平滑性都较好，如图3e所示.从修复方法的复杂程度来看，滚球法和凸包算法的过程都较为简单，而且仅在二维图像上即可完成，算法的运行速度也较快.本文所使用的方法运算过程较为复杂，需要三维建模等较为耗时的运算，但可以直接从三维角度观察肺实质，使得修复后的效果更加直观清晰.

a.原图;b.r=5，滚球法;c.r=15，滚球法;d.凸包算法;e.本文方法.图3 本文的修复方法与经典的“滚球法”和“凸包算法”相比较Fig.3 Patching algorithm in this paper is compared with the classical “ball rolling method” and “convex hull algorithm”

2.3 分割正确率

使用本文方法对选取的10组肺部CT数据进行全肺分割，计算分割正确率

.

(6)

其中正确分割的图像由放射科医生进行判断.经过实验统计计算，表1为应用本文方法和应用滚球法与凸包算法对全肺进行分割结果的正确率统计表.从统计结果可以看出，应用本文的方法对肺实质进行分割的正确率优于滚球法和凸包算法.对失败的几个案例进行分析可知其主要原因是粘连的肿瘤过大，导致肺实质的边界缺陷过大，无法进行三维重建.这种情况下继续使用本文算法则会出现欠分割的现象，所以对于有大肿瘤粘连的肺实质的修复方法需要进一步研究.

表1 不同分割方法的正确率比较

3 结束语

在实际的计算机辅助临床诊断中，准确分割肺实质对后续肺肿瘤的分割以及肺肿瘤的检测和分析等都具有十分重要的意义.本文针对肺实质边界曲率变化较大处因肿瘤粘连导致的肺实质边界缺陷，由于该位置处的缺陷在二维图像上缺少用于指导修复的特征信息，所以经典的凸包算法或滚球法的修复效果较差，故提出一种基于三维曲面重建的修复方法.从实验结果看出，该方法能够通过提取肺实质缺陷处的三维表面信息对修复曲面进行指导，修复的效果更理想，只是相应的计算时间较长.因此，若缺陷位于肺实质边界较平滑的位置，使用凸包算法进行修复即可得到良好的修复效果，并且凸包算法需要的人工干预较少，计算速度更快.若缺陷位于肺实质边界曲率变化较大的位置，采用本文的曲面重建修复方法可以得到更好的修复效果.临床医生对使用本文方法进行肺实质修复的结果进行分析后认为，该方法能够有效地修复此类肺实质缺陷，对临床病理学的分析研究有一定的辅助作用.通过本文的修复方法得到的肺实质可以完整地包含肿瘤，对肿瘤的分割也具有参考价值，后续将从肺肿瘤的三维重建方向进行深入研究，以实现肺肿瘤精确和有效地分割.