脑血肿CT图像分割的改进快速行进方法

杨泽富，战荫伟，杨荣骞

(1.广东工业大学 计算机学院，广东 广州 510006；2.华南理工大学 材料科学与工程学院，广东 广州 510640)

0 引 言

脑血肿CT图像的分割是规划脑出血手术的基础，由于手动分割及半自动分割均非常耗时，所以应采用效率更高的全自动分割方法。目前医学图像分割领域已有大量的研究工作，其中多数算法用于分割二维切片。如文献[1]利用极限树设计分类器识别病灶，解决了手动分割重复性差的问题且稳定性较强。文献[2]使用随机森林识别血肿区域，但对小面积区域识别较差。文献[3]利用离散曲率波变换对图像进行多尺度几何分析，进而区分正常脑组织和病变组织。文献[4]虽然改进并减少水平集的演化时间，但仍需初始化零水平集，分割效率较低。

临床上对病灶进行三维重建可以辅助放射科医师更加直观地了解病症，减小误诊率。因为体数据是三维重建的基础，所以针对体数据分割技术的研究逐渐成为热点。如文献[5]采用有监督学习方式分析病灶的邻域数据从而实现自动分割，分割精度较高。文献[6]结合社会群体优化算法和模糊熵理论提出一种半自动图像处理方法，在ISLES 2015竞赛中有较好的表现。文献[7,8]使用卷积神经网络识别病灶，分割结果较为准确，但算法对数据集的依赖性较强，稳定性较差。以上算法虽然都能取得较好的分割结果，但分割速度较慢，为此文献[9]提出一种特殊的三维水平集算法，通过增加多个零水平集来加快分割速度。文献[10]在FM方法的基础上提出一种Eikonal方程快速解法，以此来加快FM方法的分割速度。

文献[9,10]虽然分割速度较快，但距离临床时间要求仍有一定的距离，由于手术计划的制定过程是串行的，下一步的执行必须使用上一步的结果，因此如果分割速度较慢，将增加整个手术计划的制定时间，降低医生的工作效率，为此本文提出一种分割精度较高且速度更快的分割方法，辅助医生更加快捷地制定手术计划。

1 算法设计基础

1.1 体数据空间结构表示

图1左侧展示了一组连续的CT图像切片序列，右侧为该序列按扫描顺序排列后得到的扫描部位体数据。体数据由体素组成，体素是基本体积元素，可以理解为三维空间内排列的点或者一小块区域。设Nx、Ny和Nz均为大于1的正整数，体数据中每一张CT图像切片可视作一个大小为Nx×Ny的二维网格。若将切片记为Ik，k表示切片的序号，则可定义

Ik={V(i,j,k);1≤i≤Nx,1≤j≤Ny}

其中，V(i,j,k) 表示空间坐标为 (i,j,k) 的体素，(i,j,k)∈Z3,Z表示整数集。同理，体数据可视为一个规模大小等于Nx×Ny×Nz的三维空间区域，若用V表示体数据，参照二维切片的集合定义，则有

V={V(i,j,k);1≤i≤Nx,1≤j≤Ny,1≤k≤Nz}

其中，Nx和Ny表示单张切片的长和宽，Nz表示体数据中切片的张数。

图1 CT图像序列和体数据

1.2 FM方法介绍

FM方法的基本思想是根据输入的速度图像进行轮廓曲线演化，并在时间图像中提取演化后的轮廓曲线，其输出结果为一幅二值掩模图像。其中速度图像上每点表示轮廓曲线在该点的演化速度，时间图像上每点表示轮廓曲线到达该点所需的时间。如果用F(x,y) 表示速度图像上点 (x,y) 的值，用T(x,y) 表示时间图像上点 (x,y) 的值，则可满足Eikonal方程

(1)

其中

因此将式(1)展开可得

结合上述方程式，图2描述了FM方法中轮廓曲线的演化过程。

图2 FM方法中的轮廓曲线演化

在FM方法中轮廓曲线从初始种子点开始进行演化，图2中使用黑色网格点表示已知点，灰色网格点表示实验点，白色网格点表示未知点。在轮廓曲线演化过程中，首先设置一点为种子点，并将该种子点作为已知点，同时将其余点设为未知点，然后通过下述步骤进行曲线演化：

(1)设置已知点的初始值，一般将初始值设为0，表示曲线演化到该点所需的时间等于0；

(2)从已知点出发向外访问其四邻域点，如果邻域点是未知点则将其状态改为实验点；

(3)根据式(1)计算并更新所有实验点的值；

(4)在所有实验点中选择值最小的一点，将其状态改为已知点；

(5)返回步骤(2)继续循环操作。

2 改进的FM方法

图3为改进的FM方法流程图，本节将根据该流程图具体描述改进的FM方法。

图3 算法流程

2.1 自动设置初始种子点集

FM方法需要手动设置种子点，属于半自动分割模式，本文设计一种自动设置初始种子点集的方式，将FM方法改进为全自动分割模式。因为CT图像中脑血肿区域与相邻正常组织之间存在差异，所以可以使用阈值处理提取种子点，但考虑到头皮组织对阈值处理的影响，因此在阈值处理前应该先提取脑实质，图4描述了具体的流程。

图4 自动设置初始种子点集流程

在自动设置种子点的过程中，体数据的每张切片大小均为512×512，设采样因子为4对体数据进行下采样操作，即在切片的横向和纵向上每隔4个体素采样一次，得到大小为128×128的较小切片。然后对下采样结果进行阈值处理提取颅骨，由于颅骨的HU最低约为150[11]，与其它组织差异较为明显，因此可将阈值函数设置为

(2)

随后根据提取的颅骨计算脑实质区域的二值模板，即统计颅骨内边界中包含的所有体素，并将它们作为二值模板的前景，剩下的作为模板的背景。

体数据乘法滤波器对输入的两份大小相等的体数据，进行体素级别的乘法操作，即将两份体数据中坐标相同的体素对应相乘，并将结果作为一份新的体数据输出。由于上述二值模板的前景对应脑实质区域且值为1，背景值为0，因此体数据乘法滤波器的输出结果为脑实质区域。下一步对脑实质区域进行阈值处理，将[40，100]设为阈值区间[12]，则有阈值函数

(3)

然后对阈值处理的结果进行数学形态学二值腐蚀操作，消除离散的噪声点同时收缩种子点区域，得到待筛选的种子点集合，即图4(e)白色区域。在确认种子点时，逐一判断每张切片上的种子点个数，如果数量大于设定值20，说明该切片上包含血肿块，保留种子点；反之则舍去该切片上的所有种子点。

2.2 确定感兴趣区域包围盒

分割过程中应该先提取体数据内的脑血肿区域作为感兴趣区域(region of interest，ROI)来减少计算量。ROI一般通过包围盒的形式表示。为了便于理解三维包围盒的计算过程，先介绍二维包围盒的计算过程。

二维包围盒的计算即在图像上对ROI加框。若用R1表示二维包围盒，则可定义

R1={p(x,y);x0≤x≤x0+l,y0≤y≤y0+w}

其中，p(x,y) 表示图像坐标等于 (x,y) 的像素，(x0,y0) 为包围盒矩形的左上角起点坐标，l和w分别为包围盒的长和宽，如图5所示。

图5 二维包围盒计算过程

二维包围盒的计算过程为：首先对原始图像进行中值滤波去除脉冲噪声，然后计算图像的梯度得到梯度图像，再通过Sigmoid函数对梯度图像进行非线性映射求得速度图像。其中Sigmoid函数公式为

(4)

X表示输入像素大小，Y表示经变换后输出的像素大小，P和Q分别表示输出像素的最大值和最小值，α和β为大于零的常数。若令P=1，Q=0，式(4)可化简为

(5)

式中：当X→-∞时，Y→1，当X→+∞时，Y→0，即对于定义域X∈(-∞,+∞)，值域Y∈(0,1)。 此处采用式(5)对梯度图像做归一化处理，经多次实验验证，当α=0.3，β=2.0时处理结果较为稳定。随后根据速度图像使用FM方法分割，遍历输出的二值图像，从而确定包围盒起点坐标、长度和宽度，求得包围盒。

在二维平面上，包围盒是矩形，而在三维空间中，包围盒是长方体。若用R2表示三维包围盒，则

R2={V(x,y,z);x0≤x≤x0+a,y0≤y≤y0+b;z0≤z≤z0+c}

体数据下采样前后同一相对位置坐标之间存在比例关系。假设空间坐标系x轴和y轴方向采样因子等于n，z轴方向采样因子等于1，则包围盒R2的起点坐标 (x0,y0,z0) 在保持相对位置不变的情况下，经下采样后将变为(x0/n,y0/n,z0)，长宽高也分别变为a/n，b/n和c，从而得到下采样后的包围盒R′2，即

R′2={V(x,y,z);(x0/n)≤x≤(x0/n)+(a/n),(y0/n)≤y≤(y0/n)+(b/n),z0≤z≤z0+c}

因此可以先求体数据下采样之后的包围盒，再根据采样因子逆推下采样之前的包围盒，从而节省计算时间。逆推包围盒的求解过程如图6所示。

图6 根据采样因子逆推包围盒

图7 感兴趣区域金字塔迭代分割

2.3 ROI金字塔迭代分割

上文介绍了快速求取ROI的方法，在该方法的基础上参考图像金字塔原理可以建立一种特殊的金字塔，ROI金字塔。ROI金字塔是一系列以金字塔形状排列的体数据集合，金字塔的每一层是三维的体而不是二维的面。ROI金字塔的底部为原始体数据ROI，往上每一层都通过对相邻层进行下采样获得，层数越高，体数据尺寸越小，图7中虚线框内为具体示例。

ROI金字塔迭代分割算法：

输入：脑血肿CT图像序列，即原始体数据；初始种子点集。

输出：分割后的二值图像序列。

(1)根据前文介绍的方法提取体数据ROI，并通过下采样建立ROI金字塔；

(2)从金字塔顶层的体数据开始，输入前文求得的初始种子点集，使用三维FM方法分割；

(3)采用最近邻插值法对分割结果上采样，上采样结果的大小和金字塔下一层相等，然后使用半径为一的探针做二值腐蚀操作，提取腐蚀结果中的白色体素作为金字塔下一层的种子点集；

(4)输入新的种子点集，使用三维FM方法对金字塔的下一层进行分割；

(5)返回步骤(3)，迭代分割直至金字塔底层为止。

分析FM方法可以发现，其总是从种子点出发向外寻找边界，且如果出发位置距离边界越近，演化计算的次数越少，反之则越多。因此如果通过某种方式使种子点尽量贴近边界内侧，就可大量减少计算量，缩短FM方法的分割时间。本文根据该点提出上述ROI金字塔迭代分割算法，即除初始种子点集之外，ROI金字塔每一层的种子点集，都是通过对上一层分割结果先上采样后腐蚀获得的，这样处理得到的种子点都将包含于分割目标内部，且种子点集合的最外侧点将贴近于目标边界，又因为种子点集合的内部区域均为已知点，FM方法会忽略内部区域直接从最外侧种子点开始计算，所以只需要少量的计算即可完成分割。总而言之，在ROI金字塔迭代分割中，通过使用上一层的分割结果来估计下一层的边界，即上一层的分割可视为粗分割，当前层在粗分割的基础上使用FM方法再做一次细分割，如此层层迭代分割下去，从而避免大量多余的计算，达到提高分割速度的目的。

3 实验结果与分析

实验系统配置为Intel(R) Core(TM) i5-8400 CPU @ 2.80 GHz，16 GB RAM，操作系统为Windows 10专业版，软件环境为Microsoft Visual Studio 2015。实验数据来源于中山大学附属第一医院神经外科的8套临床CT数据集，共计2536张，数据集层厚均为0.5 mm，单张CT图像大小为512×512。实验中使用专家的手工分割结果作为分割的金标准，同时将本文方法与原FM方法，以及文献[9,10]进行对比，以此验证本文方法的有效性。

3.1 评估指标

本文对分割结果的评估指标主要包括：Dice相似系数(Dice coefficient)，灵敏度(Sensitivity)，特异性(Specifi-city)和准确率(Accuracy)。Dice相似系数用来衡量真实结果与预测结果之间的相似性，预测结果越接近真实结果，Dice值越大。若用U1表示Dice的大小，则

(6)

灵敏度又称召回率(Recall)，或称真阳性率(true positive rate，TPR)，用来衡量算法对正样本的识别能力，若用U2表示TPR的值，则

(7)

特异性又称为真阴性率(true negative rate，TNR)，用来衡量算法对负样本的识别能力，若用U3表示TNR的值，则

(8)

准确率表示分割正确的像素个数占总像素个数的比例，若用U4表示准确率，则

(9)

3.2 实验对比

为了直观地观察算法分割效果，图8给出了包括本文方法在内的3种不同方法，分别提取脑血肿轮廓的结果对比。

图8 脑血肿图像分割结果对比

图8(a)为原始CT图像；图8(b)为专家手工分割结果，并将其作为实验分割的金标准；图8(c)为文献[10]的分割结果，文献[10]是目前关于FM方法最新的改进方法；图8(d)为本文方法的分割结果。通过对比可以发现，本文方法的结果与专家手工分割结果较为接近，分割较为精确。

为了更加客观地反映本文分割算法的有效性，将本文方法分割结果的各项评估指标与原FM方法，文献[9]和文献[10]的各项评估指标进行对比，其中文献[9]为最新的水平集分割脑血肿方法，对比结果见表1。

表1 本文方法与其它方法的评估指标对比

由表1可知，在本文方法分割结果的各项评估指标中，除TNR一项之外，其它评估指标均为最优，说明本文方法是有效的。

在临床应用中，分割的时间效率也是需要考虑的重要因素。因为每套CT数据集的张数不同，所以本文选用单张CT图像的平均处理时间来评估分割的时间效率。为了进一步说明本文方法在时间效率上的优势，将其与其它方法进行对比，对比结果见表2。从表2可以看出，本文方法分割单张CT图像平均用时0.043 s，分割效率明显优于其它3种方法，特别是在处理的图像数量较大时，分割效率的优势将更加明显。

表2 单张CT图像的平均分割时间/s

4 结束语

本文对FM方法进行改进，提出一种三维全自动脑血肿CT图像序列分割方法，实现了全自动、迅速、准确的脑血肿区域轮廓提取。通过简单的阈值处理提取种子点，解决了FM方法需要手动设置种子点的问题，使算法得以全自动执行。建立一种特殊的感兴趣区域金字塔进行迭代分割，并采用先做粗分割，再做细分割的设计思想，从而在保证分割精度的同时，极大地减少了算法的计算量。实验结果表明，本文方法对比其它分割方法，不仅在分割精度上占优，而且在分割效率上也有着明显的优势，说明其具有高效性，实用性和有效性，能够为临床手术计划的制定提供有利条件。