基于LevelSet的超声乳腺肿瘤图像的轮廓提取

黄韫栀 HUANG Yunzhi

刘 奇 LIU Qi

四川大学电气信息学院医学信息工程系四川成都 610065

乳腺癌发病率居女性癌症首位，近年呈不断上升趋势，因此开展乳腺癌的诊断和防治研究具有重要意义。超声诊断以其无损伤、非侵入性、灵敏度高、适用于鉴别软组织等特点，越来越多地应用于乳腺肿瘤的辅助检测，而对乳腺肿瘤超声图像进行分割以获取肿瘤的边界是判断乳腺肿瘤良恶性的关键。由于超声图像存在大量的斑点噪声，图像质量明显下降，造成医学图像的复杂性和多样性，从而大大增加了乳腺肿瘤分割的难度。本研究旨在采用曲率各向异性扩散的方法，在保留主要边界信息的前提下对超声乳腺图像加以平滑，并利用水平集的方法自动提取乳腺肿瘤轮廓，为后续研究打下基础。

1 资料与方法

1.1 肿瘤粗定位 锁定目标，提取肿瘤所在的方形区域（图1），以减小后期进行自动轮廓提取的迭代时间，并增加迭代的成功率。

1.2 保留边缘平滑滤波 传统的图像平滑算法如均值滤波、中值滤波和高斯滤波等，属于各向同性扩散，这些平滑算子在去除噪声的同时会模糊甚至破坏图像的边缘[1]。为了在去除噪声的同时保持图像的边缘，Perona等[2]提出在不同的方向设置不同分布系数的各向异性扩散滤波（anisotropic diffusion），建立了各向异性扩散方程：

图1 肿瘤粗定位

其中原始图像f (x, y)作为其初始值，在不同时刻可以得到在此时刻下扩散滤波后的图像。若设置：

为扩散程度的主控因子，其中X＝（x,y）代表图像中像素点的位置，g（X）表示每个像素点的梯度强度[3]。显然，C（X）是梯度幅值的减函数。采用该方法对乳腺超声图像的滤波效果见图2。

图2 A.包含肿瘤的区域；B.基于梯度各向异性扩散滤波（步长=0.25，迭代次数=5，传导系数=5）

本次实验中，以公式（3）作为主控因子表达式的曲率各向异性扩散（MCDE）[4]对图像进行滤波。MCDE方程见公式（4）。

对比图2B和图3B可知，与经典的Perona-Malik扩散方式相比，MCDE具有更低的敏感度，在去除噪声和保留超声图像的边界信息方面优于梯度各向异性扩散，并能够更好地保留图像中的结构细节。此外，对比图2B和图3C可知，MCDE使用更少的迭代即可以得到一个可接受的结果。因此，采用MCDE对超声医学图像进行预处理，不仅能够减少图像中的斑点和噪声，同时能够保留乳腺肿瘤的对比度和形状。1.3 边缘的潜在图像 图像的梯度强度广泛应用于图像分析，帮助检测物体的轮廓以及分离均匀区域。为了能够得到更加准确的轮廓，先将经过MCDE滤波后的图像再次进行滤波——高斯滤波（选择合适的标准差，控制图像边缘的影响范围），然后计算图像的梯度。滤波后的梯度图像见图4，其中σ＝1.0。

图3 A.包含肿瘤的区域；B、C.基于曲率各向异性扩散滤波（B：步长=0.25，迭代次数=5，传导系数=5；C：步长=0.25，迭代次数=3，传导系数=5）

图4 梯度图像

在曲面演化的过程中，需要设置演化终止时的条件，本实验中终止的条件就是曲面到达了图像边缘。为了产生边缘的潜在图像，本次实验中选用Sigmoid函数[5,6]：

其中g是输入像素的梯度，g'是输出像素的梯度，gmax、gmin是经过MCDE滤波以后所得图像的最大梯度值和最小梯度值，α定义了输入梯度值范围的宽度，β定义了围绕在范围中心的梯度。所得潜在边界见图5。

图5 潜在的肿瘤边界

Osher等[7,8]于1982年提出了基于几何变形模型的水平集方法。利用水平集实现图像分割的典型用法为：由用户先定义出待分割对象的初始轮廓，再利用水平集的进化功能，使水平集一直进化到符合待分割对象的解剖结构为止。

在寻找肿瘤边界的过程中，不能保证速度函数F恒正[9-12]，因而不能利用FastMarching[13]提取乳腺肿瘤轮廓，需采用动态水平集：

其中初始条件为距离函数：

2.1 简单的初始轮廓 按照以下方式生成一个初始水平集曲面：①待分割区域内用户提供一个用于轮廓扩张的种子点（a，b），一个好的种子点集将增大分割一个复杂对象而不丢失数据的可能性。在肿瘤轮廓提取实验中，由于肿瘤的特殊形状，一般设置轴线上的中心位置为种子点。②以种子点为圆心，建立一个半径为r的圆，该圆作为初始的轮廓曲线，初始水平集曲面的方程表达式为：

水平集的初始轮廓见图6。

图6 水平集的初始轮廓

2.2 乳腺轮廓的提取 曲面随着时间而变化，嵌入其中的零水平集曲线的拓扑结构也会随之发生改变。在曲面演化过程中，可利用图像信息，如灰度均值、梯度和边缘等来控制曲面演化的过程。通用的曲面运动偏微分方程（PDE）为：

其中A（X）是水平对流系数，将水平线集引向对象的边界；P（X）是传播（膨胀）系数，用于初始轮廓的扩散；其中Z（X）是一个曲率均值的空间调节器系数：高曲率区域的平滑系数，进一步将噪声平滑，曲率越大，分割结果越平滑，然而曲率缩放比例太大会把轮廓移出形状边界[15]。常数α、β、γ是每个系数在界面运作上相关影响的权值。因而曲面演化的方程表达式即为：

2.2.1 测地活动轮廓（GAC） Caselles等[16]阐述了公式（10）的具体实现，即GAC的求解。演化过程需要定义达到收敛的标准或达到迭代器的最大数量，和边界的潜在图像一起构成水平集的停止条件。本实验的收敛标准是以水平集函数的均方根变换的形式来定义的。如果均方根变换在用户指定的一个门限之下，运动将会聚合到一点。具体参数设置：收敛标准：RMS=0.02，迭代器的最大数量为1200。

图7显示了利用GAC采用表1中的参数所得的图像轮廓。实验过程中测地线方法中的参数和传播参数很难控制，并且从测试的结果也发现最后得到的速度图像收敛到的边缘处并不清晰，与原图结构有很大差异。图7B中孔洞即是分割的泄露。

图7 测地活动轮廓水平集分割肿瘤

表1 GAC水平集分割的参数设置（曲率缩放=1.0，水平缩放=1.0）

对于超声乳腺肿瘤图像的分割，需要一个相对较大的传播缩放比例：①超声图像的边界对比较低；②乳腺肿瘤生理结构的边界形状比较复杂。但是这些缩放比例参数的最优值只能通过不断地实验得到。

2.2.2 阈值水平集分割 由于速度函数的定义直接关系到水平集方法的分割质量，针对GAC方法在连接成员方案中出现的“泄漏”，因此对GAC水平集的分割算法中的速度函数进行修改，并且定义一个灰度值范围来对相关的组织类型继续分类，然后求出对那个亮度范围基于水平集等式上的传播系数。任意点的速度函数F：

缺省了公式（10）中的水平对流系数。速度函数F中，需满足图8所示的要求传播项：

其中g（X）为输入的待分割图像中各点的灰度值，U为阈值上限，L为阈值下限[17]。

对于公式（12）中阈值L和U的选取，可以参照以下2点：①乳腺肿瘤超声中，目标处于低灰度值范围，因此L值一般采用ROI区域中的最小灰度值。下门限设置成负数，以便确保被分割对象的内部出现在二值区域内。②U值采用大津阈值法（Ostu）[18]，即类间方差最大化或类中差异最小化进行估计。

图8 传播系数函数

实验结果见图9，其中各参数设置见表2。下门限设置成一个大的负数以便确保被分割对象的内部将出现在二值区域内。阈值水平集的方法改进了GAC和FastMarching分割方法在弱边界处出现“泄漏”，并且提取到与图10手工勾边较为相似的轮廓，较好地包裹了肿瘤区域。

图9 阈值水平集分割结果。图10 医师手动勾画的轮廓

表2 阈值水平集分割的参数设置（曲率缩放=1.0）

2.2.3 阈值水平集分割定量分析 为了评价分割结果的精度，以医师手工勾画的结果作为“金标准”，图11显示了自动提取病灶区域与“金标准”之间的关系。本文采用Udupa[19]提出的3种不同的测量误差对99幅乳腺肿瘤（61幅良性肿瘤，38幅恶性肿瘤）图像的分割结果进行定量分析：假阳性（FP）、假阴性（FN）、真阳性（TP）：

图12和表3显示了99例肿瘤利用阈值水平集进行自动分割后得到的FP、FN、TP的测量结果。由于获取的乳腺超声图像会出现边界极不清晰的情况，使得自动分割未获得满意的结果（<5%，如图12A中第42例和图12B中第26例），此时不能完全依靠自动分割，还需要医师进一步调整肿瘤轮廓以获得比较满意的结果。

图11 自动分割区域和手动分割区域真阳性、假阳性、假阴性比较

图12 阈值水平集获取的超声乳腺肿瘤轮廓（图像均采自四川大学华西医院超声科）

表3 99例肿瘤分割的误差测度（99例图像均采自四川大学华西医院超声科）

3 结束语

本文首先定位到超声图像中的肿瘤区域，然后针对乳腺肿瘤超声图像自身的一些缺点如图像对比度低、斑点噪声大、部分边缘缺失、肿瘤内部微细结构分布复杂（如血管、钙化灶等），特别是恶性肿瘤还具有复杂形状等，提出了针对性的MCDE的方法进行滤波处理，平滑了噪声，但边缘信息丢失甚少，对于后期图像处理至关重要。乳腺肿瘤边界的提取采用LevelSet方法，其中最关键的是速度函数的确定，不同速度函数导致不同的分割效果，不同图像的目标分割需要的速度函数也不尽相同。本实验首先设定肿瘤的初始轮廓，然后测试了GAC以及阈值水平集方法在肿瘤自动分割中的效果。然而GAC方法的实验参数较难控制，容易出现分割泄露，分割效果不佳；相对而言，阈值水平集分割只需要通过Ostu方法找到阈值，通过小范围的偏移就可以获得与手动勾边相似的轮廓，大幅度改善了GAC的分割效果。自动提取肿瘤边界为后续肿瘤良恶性的判断打下了重要基础，也正朝着解放医师读片工作、自动诊断肿瘤的方向上前进。

[1] Gonzalez RC, Woods RE. Digital image processing. 2nd ed.Beijing: Publishing House of Electronics Industry, 2008: 183-215.

[2] Perona P, Malik J. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on Pattern Analysis Machine Intelligence, 1990, 12(7): 629-639.

[3] Black MJ. Robust anisotropic diffusion. IEEE Trans Image Process, 1998, 7(3): 421-432.

[4] Huang YL, Jiang YR, Chen DR, et al. Level set contouring for breast tumor in sonography. J Digit Imaging, 2007, 20(3): 238-247.

[5] 周振环, 伍云智, 赵明. 医学图像编程技术. 北京:电子工业出版社, 2010: 9-11.

[6] 周振环, 王安明,王京阳. 医学图像分割与配准Ⅰ. 成都:电子科技大学出版社, 2007: 111-113.

[7] Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on hamilton-jacobi formulation. J Computer Physics, 1988, 79(1): 12-49.

[8] Sethian JA. Level set methods and fast marching methods. 2nd ed. Cambridge: Cambridge University Press, 1999.

[9] 徐倩. 基于Level Set模型的脑血管图像骨架提取研究. 苏州:苏州大学硕士学位论文, 2009.

[10] 陈波,赖剑煌.用于图像分割的活动轮廓模型综述.中国图象图形学报, 2007, 12(1): 11-20.

[11] 师为礼,苗语,颜雁. 基于ITK的水平集医学图像分割算法.长春工业大学学报(自然科学版), 2008, 29(5): 516-519.

[12] 王峥, 杨新, 李俊, 等. 基于水平集人机交互模型的医学图像分割. 模式识别与人工智能, 2002, 15(4): 392-396.

[13] Sethian JA. A fast marching level set method for monotonically advancing fronts. PNAS, 1996, 93(4): 1591-1595.

[14] 郑国贤. LevelSet方程的求解与应用. 杭州: 浙江大学硕士学位论文, 2006.

[15] Sethian JA. Level set methods: evolving interfaces in geometry,fluid mechanics, computer vision and materials science.Cambridge: Cambridge University Press, 1996.

[16] Caselles V, Kimmel R, Sapiro G. Geodesic active contours.International Journal on Computer Vision, 1997, 22(1): 61-97.

[17] Ibanez L, Schroeder W, Ng L, et al. ITK software guide. 2nd ed. New York: Kitware Inc, 2005: 531-558.

[18] 夏良正. 数字图像处理(修订版). 南京:东南大学出版社,1999: 222-223.

[19] Udupa JK. A methodology for evaluating image segmentation algorithms. Proc SPIE, 2002, 4684(2): 266-277.