一种融合边缘和区域特征的ct图像中主动脉的分割方法
【专利摘要】本发明一种融入边缘和区域特征的CT图像中主动脉的分割方法，融入边缘和区域特征，借助水平集思想，能得到清晰准确的目标轮廓，实现CT图像中主动脉的分割，运行速度非常快且分割结果准确，较好地解决了CT图像中主动脉的提取问题，为医学诊断提供一个参考。
【专利说明】
一种融合边缘和区域特征的CT图像中主动脉的分割方法
技术领域
[0001] 本发明属于医学图像处理技术领域，尤其涉及一种基于水平集方法融合图像边缘 和区域特征的CT图像中主动脉的分割方法。
【背景技术】
[0002] 图像分割是图像处理中的一项关键技术。CT(Computed Tomography即计算机断层 扫描)能够提供高分辨率高质量的CT图像。CT图像对医学组织器官的解剖结构和功能作了 准确的描述。将CT图像中医生感兴趣的区域提取出来可以方便医生的临床诊断和手术治 疗。CT图像不同于一般网络图像，其成像机理复杂，数据特殊，传统的一些分割方法已无法 满足其图像分割的需要。
[0003] 尽管医学图像分割的研究众多，终究没找到一种具有普适性的分割方法。对于CT 图像中主动脉的分割，孙浩、陆雪松、刘海华等人在论文《CT图像中的主动脉分割方法》中提 出联合CV模型边界的区域分割算法，该算法以基于种子点的区域生长作为区域分割的生长 方式，以CV模型得到的最终演化轮廓曲线为约束条件，两者结合来达到分割主动脉区域的 目的。该方法对清晰度高噪音少的完整的主动脉可以较精确的分割出来，但是一旦图像有 噪声，或者出现主动脉夹层现象，效果便不太理想。沈阳市第四人民医院在其申请的专利 "一种基于延展性主动轮廓模型的眼底图像血管分割方法"（专利申请号CN201510274619， 公开号CN104867151A)中公开了一种利用眼底血管呈线状这一特征构建方向场，采用延展 性主动轮廓分割模型对眼底血管进行分割的方法，该方法利用了血管的特性，却不适合用 于主动脉的分割。因此，我们需要找到一种精确度高，适应性强的适用于主动脉的分割方 法。
【发明内容】
[0004] 本发明的目的在于针对上述现有的方法不适于主动脉分割的应用，提出一种融入 边缘和区域特征的CT图像中主动脉的分割方法。
[0005] 本发明一种融合边缘和区域特征的CT图像中主动脉的分割方法，包括如下步骤：
[0006] 步骤1、构建传统CV模型的能量函数E:
[0007] 设图像定义域为Ω的分割图像I(x，y)被闭合轮廓线C划分为两个同质区域，即闭 合轮廓线C的内部和外部，各区域的平均灰度值分别为&和(： 2，能量泛函E(C)可表示如下：
[0008]
⑴
[0009] 其中，x，y为像素点在图像中的横纵坐标；
[0010] 若闭合轮廓线C位于分害_像1&，7)的实际边界(：〇的外部时』怔1(〇>0而(〇~ 〇,若闭合轮廓线C位于分割图像I(x，y)的实际边界Co的内部时，则E 1(C)~0，E2(C)>0;若闭 合轮廓线C跨越分割图像的实际边界Co的内外部，E 1(C)XhE2(C)X);若闭合轮廓线C位于 分割图像I(x，y)的实际边界Co上时，E 1(C)~0，E2(C)~0,当闭合轮廓线C与实际边界Co重合 时，E(C)取得最小值；
[0011]在上式能量泛函E(C)中加入长度惩罚项μ，得到CV水平集模型的能量函数EW1, C2)为：
[0012]
[0013]其中，Length(C)表示闭合轮廓线C的长度，Pj1J2是各个能量项权重的预置系数， μ ^ 0，λι、入2 > 0;
[0014] 为了求得能量函数EW1, C2)的最小值，使用水平集的思想，将未知的演化曲线C， 即上述提到的闭合轮廓曲线C，用水平集函数Φ(χ， 7)来代替，并且设定如果点(x，y)在演化 曲线C的内部，则水平集函数Φ (X，y )>0;如果点（X，y)在演化曲线C的外部时，则水平集函 数Mx，y)<〇;而如果点（x，y)恰好在演化曲线C上面时，则Φ(χ，γ) = 〇,通过引入水平集函 数Φ (x，y)，能量泛函E(C)可被重写为如下的形式：
[0015]
[0016] 上式中，V为梯度符号，即表示分别对函数Φ(χ，γ)中的x，y求导，Ηε(ζ_δε( ζ)分 别是海氏函数H(Z)和狄拉克函数δ(ζ)的正则化形式，海氏函数H(Z)的正则化形式Η ε(ζ)为：
[0017]
[0018]
[0019]
[0020] 上式中ε为正则化参数；
[0021] 步骤2、在传统CV模型中加入能量惩罚项，使得水平集函数Φ(Χ，7)在演化过程中 保持为近似的符号距离函数，避免重新初始化水平集函数，此能量惩罚项表达式为：
[0022]
(6)
[0023] 其中V表示梯度符号
[0024] 步骤3、在传统CV模型中加入圆形先验知识约束形变模型的收敛形状，使曲线的演 化结果与先验形状特征相一致，使用符号距离函数Φ来隐含表示目标的形状，构造先验形 状约束项：
[0025]设图像定义域Ω，符号距离函数φ表示演化曲线，φρ为先验形状的符号距离函 数，则形状先验能量项Ερ( Φ，ΦΡ)表示为：
[0026]
[0027] 定义圆形特祉公式为：
[0028]
[0029] 则形状约束能量泛函Ecircle表示为：
[0030]
(9)
[0031 ]上式中，XQ, yo为圆的圆心坐标，r为圆的半径，他们随着曲线的演化也在不断地变 化；
[0032] 步骤4、综合上述式(3) (6) (9)构造最终的能量函数E:
[0033]
[0034] 通过固定水平集函数Φ，对图像的灰度均值函数&和&进行求导，通过变分法可以 得到(Μ Φ )和&( Φ )的如下表达形式：
[0037] 利用欧拉-拉格朗日法推导式（10)得到水平集函数Φ的欧拉-拉格朗日方程，引入 时间变量t后，可以得到如下的水平集演化方程：
[0035] (12)
[0036] (13)
[0038]
[0039] 公式(11)便为最终得到的新的水平集演化方程；
[0040] 步骤5、利用步骤4得到的新的能量函数E对CT图像中的主动脉进行分割：
[00411步骤5.1、将输入的CT图像序列文件统一转换为.jpg格式，CT图像属于中分辨率图 像(512*512*12)，所以格式转换后的图像大小为512*512;
[0042] 步骤5.2、读入步骤5.1中统一格式的CT图像序列文件，并对第一张图像进行剪裁 和滤波预处理，得到初始感兴趣区域；
[0043] 步骤5.3、将步骤5.2得到的初始感兴趣区域图像的中心设置为中心点，设置初始 轮廓线为圆形；
[0044] 步骤5.4、设置水平集方法的参数：
[0045] 步骤5.5、迭代求解式(11)的水平集演化方程；
[0046] 步骤5.6、从水平集函数中提取零水平集，即演化曲线；
[0047] 步骤5.7、设置水平集自适应停止函数，在水平集演化过程中，演化曲线不断向目 标轮廓逼近，当演化曲线最终演化到真实轮廓的时候，不管是否到达设置的迭代次数，演化 曲线都应该停止演化，即如果自适应停止函数得到满足，算法结束，则第一张图像分割完 毕，将分割结果保存后转步骤5.8，否则返回步骤5.6;
[0048] 步骤5.8、下一张图像以前一张图像的分割结果为依据，以前一张图像分割区域中 心点为中心点，进行感兴趣区域的初始裁剪，并将此点作为初始轮廓线的中心点；
[0049] 步骤5.9、重复步骤5.3至5.8，其中步骤5.4的参数设置在第一次设置好之后可省 略，待所有CT图像分割完毕，得到最终的分割结果。
[0050] 所述步骤5.2得到初始感兴趣区域的具体做法为：使用Matlab中图像裁剪函数 imcrop将主动脉在CT图像中的区域作为感兴趣区域裁剪出来，然后对感兴趣区域进行中值 滤波，去除椒盐噪声，得到初始感兴趣区域。
[0051] 所述步骤5.4设置的水平集方法的参数为:迭代次数n = 100,时间步长At = O. 1， 正则化参数ε = 1，长度惩罚项μ = 0，内部能量项权重参数λι = 0.5，λ2 = 1，形状约束项控制 参数入3=1。
[0052] 采用本发明的技术方案，融入边缘和区域特征，借助水平集思想，能得到清晰准确 的目标轮廓，实现CT图像中主动脉的分割，运行速度非常快且分割结果准确，较好地解决了 CT图像中主动脉的提取问题，为医学诊断提供一个参考。
【具体实施方式】
[0053]本发明一种融合边缘和区域特征的CT图像中主动脉的分割方法，在传统CV水平集 模型中融入分割图像的边缘特征，加入圆形作为先验能量项，构建出新的能量函数Ε，其具 体实现步骤如下：
[0054]步骤1、构建传统CV模型的能量函数E:
[0055] CV模型是一个典型的基于区域的水平集模型，设图像定义域为Ω的分割图像Ι(χ， y)被闭合轮廓线C划分为两个同质区域，即闭合轮廓线C的内部和外部，各区域的平均灰度 值分别为&和&，能量泛函E(C)可表示如下：
[0056]
⑴
[0057]在上式中，x，y分别表示像素点在图像中的横纵坐标；
[0058] 若闭合轮廓线C位于分割图像I(x，y)的实际边界Co的外部时，则E1(C)X), E2(C)~ 〇,若闭合轮廓线C位于分割图像I(x，y)的实际边界Co的内部时，则E1(C)~0，E 2(C)>0;若闭 合轮廓线C跨越分割图像的实际边界Co的内外部，E1(C)XhE 2(C)X);若闭合轮廓线C位于 分割图像I(x，y)的实际边界Co上时，E1(C)~0，E 2(C)~0,因此当闭合轮廓线C与实际边界 Co重合时，E(C)取得最小值；
[00611其中，Length(C)表示闭合轮廓线C的长度，μ彡0山八2>0是各个能量项权重的系 数，这些系数都是预先设置好的，一般情况下可固定取P = ^1 = A2 = I;
[0059] 在上式能量泛函E(C)中加入长度惩罚项μ，得到CV水平集模型的能量函数EW1, C2)为 r
[0060]
[0062] 为了求得能量函数EW1, C2)的最小值，使用水平集的思想，将未知的演化曲线C， 即上述提到的闭合轮廓曲线C，用水平集函数Φ(χ， 7)来代替，并且设定如果点(x，y)在演化 曲线C的内部，则水平集函数Φ (X，y )>0;如果点（X，y)在演化曲线C的外部时，则水平集函 数Mx，y)<〇;而如果点（x，y)恰好在演化曲线C上面时，则Φ(χ，γ) = 〇,通过引入水平集函 数Φ(χ.ν)，能量泛函E(C)可被重写为如下的形式：
[0063]
[0064] 上式中，▽为梯度符号，即表示分别对函数Φ (x，y)中的x，y求导，Ηε(ζ)和δε(ζ)分 别是海氏(Heaviside)函数H(Z)和狄拉克(Dirac)函数δ (ζ)的正则化形式，海氏函数的作用 是划分演化区域，狄拉克函数的作用是用来限定演化在零水平集函数周围取值，在实际问 题中，用数值方法直接计算海氏函数较为困难，因此选择其非紧支撑、光滑和严格单调的正 则化形式ft (Z):
[0065]
[0066] 上式中ε为正则化参数；
[0067] 而狄拉克函数δ(ζ)正则化形式δε(ζ)可用下式来计算：
[0068]
[0069] 步骤2、在传统CV模型中加入能量惩罚项避免重新初始化水平集函数：
[0070] 在很多情况下，水平集函数Φ(χ，7)在演化了一段时间之后会发生震荡，并渐渐失 去光滑性和距离函数特性，随着演化的继续，误差会不断累积并造成最终的曲线偏离真实 的目标边缘。在一般水平集方法中，一般会重新初始化水平集函数。本发明不采用重新初始 化方式，而是引入一个能量惩罚项，使得水平集函数Φ(χ，7)在演化过程中保持为近似的符 号距离函数，此能暈惩罚项由李纯明教授在2005年提出，其表达式为：
[0071]
(6)
[0072] 其中▽表示梯度符号 [0073]步骤3、构造形状能量约束项
[0074]在传统CV模型中加入形状先验知识可以约束形变模型的收敛形状，使曲线的演化 结果与先验形状特征相一致，可以解决非目标区域信息的干扰，一般使用符号距离函数Φ 来隐含表示目标的形状，构造先验形状约束项：
[0075]设图像定义域Ω，符号距离函数φ表示演化曲线，φρ为先验形状的符号距离函 数，则形状先验能量项Ερ( Φ，ΦΡ)表示为：
[0076] Ερ( Φ , Φρ)=/ω(Η( Φ (x,y))-H( ΦP(x,y)))2dxdy (7)
[0077] 本发明中，定义圆形特征公式为：
[0078]
[0079]
[0080] ⑶
[0081 ]上式中，xo，yo为圆的圆心坐标，r为圆的半径，他们随着曲线的演化也在不断地变 化；
[0082] 步骤4、综合式(3) (6) (9)构造最终的能量函数E:
[0083]
[0084] 通过固定水平集函数Φ，对图像的灰度均值函数Ci和C2进行求导，通过变分法可以 得到(Μ Φ )和&( Φ )的如下表达形式：
[0085] (12)
[0086] (⑶
[0087] 利用欧拉-拉格朗日法推导式（10)得到水平集函数Φ的欧拉-拉格朗日方程，引入 时间变量t后，可以得到如下的水平集演化方程：
[0088]
[0089] 公11)使为最终得到的新的水干集演化万桎；
[0090] 步骤5、利用步骤4得到的新的能量函数E对CT图像中的主动脉进行分割：
[0091] 步骤5.1、将输入的CT图像序列文件统一转换为.jPg格式，CT图像属于中分辨率图 像(512*512*12)，所以格式转换后的图像大小为512*512;
[0092]步骤5.2、读入步骤5.1中统一格式的CT图像序列文件，并对第一张图像进行剪裁 和滤波预处理，得到初始感兴趣区域，具体为：
[0093]使用Matlab中图像裁剪函数imcrop将主动脉在CT图像中的区域作为感兴趣区域 裁剪出来，然后对感兴趣区域进行中值滤波，去除椒盐噪声，得到初始感兴趣区域；
[0094]步骤5.3、对步骤5.2得到的初始感兴趣区域设置初始轮廓线：
[0095]本发明主要依据基于区域特征的水平集方法，其基本思想就是利用迭代的方法使 初始轮廓线不断接近目标边缘，本发明中将初始轮廓线设为圆形，其中心点为初始感兴趣 区域图像的中心，半径为10;
[0096I步骤5.4、水平集方法的参数设置：
[0097I 设置本发明的迭代次数n = 100，时间步长Λ t = 0.1，正则化参数ε = 1，长度惩罚项 μ=0，内部能量项参数λι=0. 5，λ2 = 1，形状约束项控制参数λ3 = 1;
[0098]步骤5.5、迭代求解式(11)的水平集演化方程；
[0099]步骤5.6、从水平集函数中提取零水平集，即演化曲线；
[0100] 步骤5.7、设置水平集自适应停止函数，判断自适应停止函数是否得到满足，如果 满足，算法结束，则第一张图像分割完毕，将分割结果保存后转步骤5.8，否则返回步骤5.6:
[0101] 在水平集演化过程中，演化曲线不断向目标轮廓逼近，当演化曲线最终到达真实 轮廓的时候，曲线应该停止演化;本发明设计一种自适应停止函数使得演化曲线到达边缘 处，不管是否到达设置的迭代次数都要停止演化，其设计原理为当演化曲线接近边缘区域 的时候，演化速度会变慢，因此当曲线长度变化的绝对值在一个给定的迭代次数里一直小 于某个事先设置的阈值，则水平集演化应终止；
[0102] 步骤5.8、下一张图像以前一张图像的分割结果为依据，以前一张图像分割区域中 心点为中心点，进行感兴趣区域的初始裁剪，并将此点作为初始轮廓线的中心点；
[0103] 步骤5.9、重复步骤5.3至5.8,其中步骤5.4的参数设置在第一次设置好之后可省 略，待所有CT图像分割完毕，得到最终的分割结果。
[0104]以上所述，并非对本发明的技术范围作任何限制，故凡是依据本发明的技术实质 对以上实施例所作的任何细微修改、等同变化与修饰，均仍属于本发明技术方案的范围内。
【主权项】
1. 一种融合边缘和区域特征的CT图像中主动脉的分割方法，其特征在于包括如下步 骤： 步骤1、构建传统CV模型的能量函数E: 设图像定义域为Ω的分割图像I(x，y)被闭合轮廓线C划分为两个同质区域，即闭合轮 廓线C的内部和外部，各区域的平均灰度值分别为&和(：2，能量泛函E(C):其中，x，y为像素点在图像中的横纵坐标； 若闭合轮廓线C位于分割图像I(x，y)的实际边界Co的外部时，则EKOX^EKC)~0,若 闭合轮廓线C位于分割图像I(x，y)的实际边界Co的内部时，则ΕΚΟΟ,ΕΚΟΧ);若闭合轮 廓线C跨越分割图像的实际边界Co的内外部，EKOXhEdOX);若闭合轮廓线C位于分割 图像I (X，y)的实际边界Co上时，EKC)~0，E2(C)~0，当闭合轮廓线C与实际边界Co重合时，E (C)取得最小值； 在上式能量泛函E(C)中加入长度惩罚项μ，得到CV水平集模型的能量函数 为：其中，Length(C)表示闭合轮廓线C的长度，是各个能量项权重的预置系数，μ^： 〇，入1、入2>0 ; 为了求得能量函数的最小值，使用水平集的思想，将未知的演化曲线C，即上 述提到的闭合轮廓曲线C，用水平集函数Φ (x，y)来代替，并且设定如果点（x，y)在演化曲 线C的内部，则水平集函数Φ (X，y) >0;如果点(X，y)在演化曲线C的外部时，则水平集函数 小（1，7)<〇;而如果点(义，7)恰好在演化曲线(：上面时，则小（1，7)=〇，通过引入水平集函数 Φ (x，y)，能量泛函E(C)可被重写为如下的形式：上式中，▽为梯度符号，即表示分别对函数Φ(χ，7)中的x，y求导，Ηε(ζΜΡδε( ζ)分别是 海氏函数Η(ζ)和狄拉克函数δ(ζ)的正则化形式，海氏函数Η(ζ)的正则化形式Ηε(ζ)为：狄拉克函数W ζ)正则化形式δε (ζ)为：上式中ε为正则化参数； 步骤2、在传统CV模型中加入能量惩罚项，使得水平集函数Φ(Χ，7)在演化过程中保持 为近似的符号距离函数，避免重新初始化水平集函数，此能量惩罚项表达式为： 其中▽表示梯度符号步骤3、在传统CV模型中加入圆形先验知识约束形变模型的收敛形状，使曲线的演化结 果与先验形状特征相一致，使用符号距离函数Φ来隐含表示目标的形状，构造先验形状约 束项： 设图像定义域Ω，符号距离函数Φ表示演化曲线，ΦΡ为先验形状的符号距离函数，则形 状先验能量项Ερ( Φ，ΦΡ)表示为： Ερ( Φ , Φρ) = /ω(Η( Φ (x,y))-H( ΦΡ(χ,γ)))2?χ?γ (7) 定义圆形特征公式为：则形状约束能量泛函Ecircle表示为：上式中，xq , yo为圆的圆心坐标，:r为圆的半径，他们随着曲线的演化也在不断地变化； 步骤4、综合上述式(3) (6) (9)构造最终的能量函数E:通过固定水平集函数Φ，对图像的灰度均值函数&和(：2进行求导，通过变分法可以得到 Cl( Φ )和C2( Φ )的如下表达形式：利用欧拉-拉格朗日法推导式（10)得到水平集函数Φ的欧拉-拉格朗日方程，引入时 间变量t后，可以得到如下的水平集演化方程：公式(11)便为最终得到的新的水平集演化方程； 步骤5、利用步骤4得到的新的能量函数E对CT图像中的主动脉进行分割： 步骤5.1、将输入的CT图像序列文件统一转换为.jpg格式，CT图像属于中分辨率图像， 分辨率为512*512*12，则格式转换后的图像大小为512*512; 步骤5.2、读入步骤5.1中统一格式的CT图像序列文件，并对第一张图像进行剪裁和滤 波预处理，得到初始感兴趣区域； 步骤5.3、将步骤5.2得到的初始感兴趣区域图像的中心设置为中心点，设置初始轮廓 线为圆形； 步骤5.4、设置水平集方法的参数： 步骤5.5、迭代求解式(11)的水平集演化方程； 步骤5.6、从水平集函数中提取零水平集，即演化曲线； 步骤5.7、设置水平集自适应停止函数，在水平集演化过程中，演化曲线不断向目标轮 廓逼近，当演化曲线最终到达真实轮廓的时候，不管是否到达设置的迭代次数，演化曲线应 该停止演化，即自适应停止函数如果得到满足，算法结束，则第一张图像分割完毕，将分割 结果保存后转步骤5.8，否则返回步骤5.6; 步骤5.8、下一张图像以前一张图像的分割结果为依据，以前一张图像分割区域中心 点为中心点，进行感兴趣区域的初始裁剪，并将此点作为初始轮廓线的中心点； 步骤5.9、重复步骤5.3至5.8,其中步骤5.4的参数设置在第一次设置好之后可省略，待 所有CT图像分割完毕，得到最终的分割结果。2. 根据权利要求1所述一种融入边缘和区域特征的CT图像中主动脉的分割方法，其特 征在于，所述步骤5.2得到初始感兴趣区域的具体做法为：使用Matlab中图像裁剪函数 imcrop将主动脉在CT图像中的区域作为感兴趣区域裁剪出来，然后对感兴趣区域进行中值 滤波，去除椒盐噪声，得到初始感兴趣区域。3. 根据权利要求1所述一种融入边缘和区域特征的CT图像中主动脉的分割方法，其特 征在于，所述步骤5.4设置的水平集方法的参数为:迭代次数11=100，时间步长八1 = 0.1，正 则化参数ε = 1，长度惩罚项μ = 0，内部能量项权重参数λι = 0.5，λ2 = 1，形状约束项控制参数 入3 = 1 〇
【文档编号】G06T7/00GK106056576SQ201610319535
【公开日】2016年10月26日
【申请日】2016年5月13日
【发明人】骆炎民, 杨珺, 柳培忠
【申请人】华侨大学, 泉州市中仿宏业信息科技有限公司