CN102243759B - 一种基于几何形变模型的三维肺血管图像分割方法 - Google Patents

Abstract

本发明提供一种基于几何形变模型的三维肺血管图像分割方法，步骤为：(1)根据人体生理结构特点确定血管分割计算区域，区域的选择不仅将待分目标完全包含，而且其形状特征稳定，避免了对全局区域进行计算，提高了分割速度；(2)血管区域均值计算，定位目标内外同质区域；(3)血管边缘能量计算，使曲面沿着在图像梯度方向的二阶导数进行演化，以使其精确收敛到目标边缘；(4)在上述基础上，建立三维血管分割曲面演化模型，有效结合肺血管内外区域均值和边缘能量；(5)根据所构建的形变模型，采用优化水平集演化求解，根据水平集函数曲面演化隐含地求解曲面运动。通过大量肺CT图像实验证明，该方法对肺血管分割快速，准确，鲁棒性强。

Description

一种基于几何形变模型的三维肺血管图像分割方法

技术领域

本发明属于数字图像处理与医学成像技术的交叉领域，特别涉及一种基于几何形变模型的三维肺血管图像分割方法。

背景技术

由于CT设备能够提供高清晰度的图像，并且为图像中各人体组织提供很高的对比度，通常应用于肺部疾病的诊断。肺血管是人体最为重要且拓扑结构最为复杂的组织器官之一。肺血管的精确分割是肺癌病灶检测的关键，这是由于在CT图像中肺血管多分布于肺癌病灶附近，而且二者灰度值相近，是肺癌病灶检测的重要干扰。此外，肺血管分割对于肺动脉栓塞、肺心病等肺部疾病的计算机辅助诊断研究也具有重要意义。

目前血管分割方法主要包含基于形变模型的方法、人工智能方法、神经网络方法及管状目标检测方法等。其中，形变模型是一类被广泛研究的方法，非常适用于血管分割。从总体上看，这些形变模型的能量方程一般包括两种类型的能量项：  一种为内力能量项，它描述了形变曲线(或曲面)自身的几何性质在变形过程中的影响；另一种为外力能量项，它描述了包括图像特征在内的与形变曲线(或曲面)本身几何性质无关的外在势力对曲线形变的影响。形变过程就是这两种力量彼此消长的过程，最后达到两种力量的平衡。

形变模型可以分为参数形变模型与几何形变模型。参数形变模型也称为Snake模型，该模型在变形过程中用参数形式显式地表示曲线和曲面，可以控制轮廓的光滑性和连续性，可以在一定程度上克服图像噪声和边缘狭缝的影响。 但是该模型仍存在几个缺点：(1)分割结果依赖初始轮廓；(2)从低维向高维扩展困难；(3)难以分割拓扑结构复杂的物体，即对分割目标的融合与分裂处理困难。

几何形变模型很好的解决了参数形变模型所带来的问题。该模型对初始轮廓的位置没有严格的要求，从而增加了该方法的鲁棒性；无弹性系数，从而避免了计算轮廓切线和顶点法线的过程；可以扩展到高维，并且从高维面函数较容易得到低一维曲线或曲面的内在几何特征，如法向量、曲率等。  现有几何形变模型主要可以分为基于目标边缘与目标区域两类。典型的基于目标边缘的几何形变模型为GAC模型，基本思想为利用目标边缘梯度、曲率等信息，使活动轮廓的运动可以表现为如何在一个势场上，找到两点之间的最小加权距离；典型的基于目标区域的几何形变模型为Chan-Vese模型，基本思想为假设图像被闭合曲线划分为目标和背景两个同质区域，可以适用于梯度有意义和无意义的轮廓检测，也就是对边缘平滑或离散的图像分割同样适用，提高了图像的分割精度和鲁棒性。

但是，现有这两类几何形变模型只包含了图像的区域与目标边缘属性中的一类信息。而且，大都只适用于三维图像分割。由于肺血管拓扑结构复杂，呈“树”状结构，从肺主动脉至分支将近17级；而且直径较细，范围从几十毫米到几毫米不等。现有几何形变模型并不能灵活适应三维血管拓扑结构的改变，从而使血管分割的准确性受到了限制。

发明内容

针对现有技术中采用的肺血管分割方法存在速度、精度不能满足应用要求等问题，本发明要解决的技术问题是提供一种满足准确性、快速性与鲁棒性的三维肺血管分割方法。

本发明提供的一种基于几何形变模型的三维肺血管图像分割方法，包括以下步骤：

(1)初始分割区域确定

读入肺CT图像数据后，由于每层图像尺寸很大，而且通常一个患者的肺序列图像有几十层至上百层。为了加快算法的运行速度，减少不必要的计算，需要在待分割的肺血管周围确定一个区域，使曲面演化只在这个区域内进行，而且要求该区域形状特征稳定，即不受不同患者及肺实质内病变影响。该步骤又包含以下步骤：

a.肋骨分割

根据解剖学知识，人体肋骨与肺壁外侧是紧密相邻的，且形状等特征不受肺实质内病变影响，故可作为初始区域的可靠参照物。由于肋骨CT值较高，与周围组织对比明显，可以采用阈值分割方法获得，阈值T可以依据肋骨CT值范围选取。

b.初始分割区域计算

在肋骨分割基础上，肋骨区域靠近图像中心近的像素点被作为肋骨关键点(P_i，i＝1，2..N)保留，然后采用三次B样条插值计算连续轮廓，插值公式如下：

$P_{i,n}\left(t\right)=\underset{k=0}{\overset{n}{\mathrm{\Σ}}}{P}_{i+k}\·F_{k,n}\left(t\right)$

$F_{k,n}\left(t\right)=\frac{1}{n!}\underset{j=0}{\overset{n-k}{\mathrm{\Σ}}}{(-1)}^j\·C_{n+1}^j\·{(t+n-k-j)}^n$

其中，0≤t≤1，P_i+k代表已有离散点，P_i，n代表插值后的连续点，n代表次数。

(2)血管区域均值计算

设原图像I(x，y)被活动轮廓线C划分为目标和背景两个区域，各个区域的平均灰度值为C₁和C₂，血管区域均值函数定义为：

E_REG＝α[∫_inside(c)|I(x，y)c₁|²dxdy+∫_outside(c)|I(x，y)c₂|²dxdy]+r∫_c ds

其中第一项为图像区域能量拟合项，第二项为几何约束项，要求曲线C尽可能短，α，r为正常数(依据经验值设定)，ds为弧长微元。

(3)血管边缘能量计算

设I_εε为在图像梯度方向的二阶导数，I_ηη为其垂直矢量。由于

血管边缘能量函数定义为：

$E_{\mathrm{EDGE}}={\&Integral;\&Integral;}_s<\&dtri;I,\stackrel{\&RightArrow;}{n}>\mathrm{da}-{\&Integral;\&Integral;\&Integral;}_{\mathrm{\&Omega;}}{H}_I|\&dtri;I|\mathrm{dxdydz}$

其中H_I为曲面的曲率，da为曲面表面区域，

为曲面单位方向矢量。

(4)建立血管分割形变模型

通过目标边缘能量的计算，可以定位目标边缘，但是当初始曲面远离实际目标时，很难精确定位所有目标边缘。本发明通过与图像目标内外区域信息的集合，建立血管分割形变模型，可以较好的解决这一问题。因此，总能量函数定义为：

$E=\mathrm{\&alpha;}[{\&Integral;}_{\mathrm{inside}\left(c\right)}{|I(x,y)-c_1|}^2\mathrm{dxdy}+{\&Integral;}_{\mathrm{outside}\left(c\right)}{|I(x,y)-c_2|}^2\mathrm{dxdy}]-({\&Integral;\&Integral;}_s<\&dtri;I,\stackrel{\&RightArrow;}{n}>\mathrm{da}-\&Integral;{\&Integral;\&Integral;}_{{\mathrm{\&Omega;}}_s}k|\&dtri;I\left|\mathrm{dxdydz}\right)+r{\&Integral;}_c\mathrm{ds}$

该公式由三项构成，第一项为图像区域均值计算；第二项为图像边缘能量计算；第三项为曲面几何约束项。其中，α为图像均值项系数，r为几何约束项系数，均为正常数。对于α与r系数值的确定，根据图像特点的不同，依据经验值选取。

(5)水平集演化

水平集方法是求解几何形变模型的有效计算工具，可以计算封闭运动曲面随时间演化过程中的几何拓扑变化。本发明采用水平集方法求解能量最小的优化问题，根据已建立的曲面演化模型，优化水平集函数定义为：

$S_t=\left\{\mathrm{\&alpha;}\right[{(I-c_1)}^2-{(I-c_2)}^2]-I_{\mathrm{\&epsiv;\&epsiv;}}-r{H}_I\}\stackrel{\&RightArrow;}{n}$

其中 $I_{\mathrm{\&epsiv;\&epsiv;}}=\&dtri;I-H_I|\&dtri;I|$

对于三维血管图像，把S_t作为高一维函数φ(x，y，z)的零水平集，满足S＝{{x，y，z}|φ(x，y，z)＝0}，相应的水平集函数定义为：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\left\{\mathrm{\&alpha;}\right[{(I_0-c_1)}^2-{(I_0-c_2)}^2]-I_{\mathrm{\&epsiv;\&epsiv;}}-r{H}_I\}\left|\mathrm{\&Delta;\&phi;}\right|$

曲面停止演化后，所得到的零水平即为最终的血管分割结果。

综上所述，本发明针对现有肺血管分割方法存在的问题，提出一种基于几何形变模型的三维肺血管图像分割方法。首先，本发明根据人体生理结构特点限定血管分割计算区域，区域的选择不仅将待分目标完全包含，而且其形状特征稳定，不受不同患者及肺实质内病变影响，避免了对全局区域进行计算，大幅度提高了分割速度。第二，本发明有效结合CT图像中肺血管内外区域均值和边缘能量，在精确定位目标内外同质区域的基础上，建立三维血管分割曲面演化模型，通过目标边缘能量的计算可以使曲面沿着在图像梯度方向的二阶导数进行演化，以使其精确收敛到目标边缘。第三，根据所构建形变模型，本发明采用优化水平集演化求解。该方法具体用有限差分方法，将闭合曲线隐含地表达为高一维的曲面函数的水平集，即具有相同函数值的点集，根据水平集函数曲面演化隐含地求解曲面运动。通过大量肺CT图像实验证明，该方法不依赖于初始位置，参数选择固定，肺血管分割快速，准确，鲁棒性强。

附图说明

图1为基于几何形变模型的三维肺血管图像分割方法流程图；

图2为血管分割初始区域计算图，其中a为CT原始图像，b为基于阈值的肋骨分割图，c为采用三次B样条插值后的初始区域图。

图3为血管局部分割演化结果。a-d分别为在不同迭代次数的分割结果，显示了形变模型从初始曲面至最终形态的分割过程。

图4为三维肺血管分割图。a为CT原始图像，b为基于几何形变模型的三维肺血管图像分割图。

图5为另一组CT图像三维肺血管分割图。a为CT原始图像，b为基于几何形变模型的三维肺血管图像分割图。

具体实施方式

本发明的实施采用多排螺旋CT(HRCT)图像数据。由于CT图像能够提供高清晰度的图像，并且为图像中各组织提供很高的对比度，通常应用于肺部疾病的诊断。随着多排螺旋CT的发展，医生可以获取更高分辨率的图像，最大限度减少局部容积效应，通过一次检测能够获得更多的患者信息，进一步拓展了CT图像的应用。

结合附图，基于几何形变模型的三维肺血管图像分割流程图如图1所示，本发明详细的分割方法，包括以下五个步骤：

(1)初始分割区域确定；

(2)血管区域均值计算；

(3)血管边缘能量计算；

(4)建立血管分割形变模型；

(5)水平集演化。

本发明方法具体实现过程如下：

步骤(1)：

读入肺CT图像数据后，由于每层图像尺寸很大，而且通常一个患者的肺序列图像有几十层至上百层。为了加快算法的运行速度，减少不必要的计算，需要在待分割的肺血管周围确定一个区域，使曲面演化只在这个区域内进行，而且要求该区域形状特征稳定，即不受不同患者及肺实质内病变影响。该步骤又包含以下步骤：

a.肋骨分割

根据解剖学知识，人体肋骨与肺壁外侧是紧密相邻的，且形状等特征不受肺实质内病变影响，故可作为初始区域的可靠参照物。由于肋骨CT值较高，与周围组织对比明显，可以采用阈值分割方法获得，阈值T可以依据肋骨CT值范围人工选取。在具体实施过程中，骨骼CT值范围选定在200-800HU之间。肋骨分割结果如图2b。

b.初始分割区域计算

在肋骨分割基础上，肋骨区域靠近图像中心近的像素点被作为肋骨关键点(P_i，i＝1，2..N)保留，然后采用三次B样条插值计算连续轮廓，插值公式如下：

$P_{i,n}\left(t\right)=\underset{k=0}{\overset{n}{\mathrm{\&Sigma;}}}{P}_{i+k}\&CenterDot;F_{k,n}\left(t\right)$

$F_{k,n}\left(t\right)=\frac{1}{n!}\underset{j=0}{\overset{n-k}{\mathrm{\&Sigma;}}}{(-1)}^j\&CenterDot;C_{n+1}^j\&CenterDot;{(t+n-k-j)}^n$

其中，0≤t≤1，P_i+k代表已有离散点，P_i，n代表插值后的连续点，n代表次数。初始区域限定结果如图2c。

步骤(2)：

设原图像I(x，y)被活动轮廓线C划分为目标和背景两个区域，各个区域的平均灰度值为C₁和C₂，血管区域均值函数定义为：

E_REG＝α[∫_inside(c)|I(x，y)c₁|²dxdy+∫_outside(c)|I(x，y)c₂|²dxdy]+r∫_cds

其中第一项为图像区域能量拟合项，第二项为几何约束项，要求曲线C尽可能短，α，r为正常数(依据经验值设定，γ比α小，当图像噪声较多时，α应当比较大，否则应该减小)，ds为弧长微元。

步骤(3)：

设图像由I(x，y)表示，其中I_x和I_y是图像在水平和垂直方向的一阶导数，图像梯度方向矢量定义为：

与此相垂直的矢量定义为：

$\stackrel{\&RightArrow;}{\mathrm{\&eta;}}(x,y)=\frac{\&dtri;I}{|\&dtri;I|}=\frac{\{-I_y,I_x\}}{\sqrt{I_x^2+I_y^2}},$ 因此， $<\stackrel{\&RightArrow;}{\mathrm{\&xi;}},\stackrel{\&RightArrow;}{\mathrm{\&eta;}}>=0.$

设I_εε为在图像梯度方向的二阶导数，I_ηη为其垂直矢量。由于血管边缘能量函数定义为：

$E_{\mathrm{EDGE}}={\&Integral;\&Integral;}_s<\&dtri;I,\stackrel{\&RightArrow;}{n}>\mathrm{da}-{\&Integral;\&Integral;\&Integral;}_{{\mathrm{\&Omega;}}_s}{H}_I|\&dtri;I|\mathrm{dxdydz}$

其中H_I为曲面的曲率，da为曲面表面区域，为曲面单位方向矢量。

步骤(4)：

通过目标边缘能量的计算，可以定位目标边缘，但是当初始曲面远离实际目标时，很难精确定位所有目标边缘。本发明通过与图像目标内外区域信息的集合，建立血管分割形变模型，可以较好的解决这一问题。因此，总能量函数定义为：

$E=\mathrm{\&alpha;}[{\&Integral;}_{\mathrm{inside}\left(c\right)}{|I(x,y)-c_1|}^2\mathrm{dxdy}+{\&Integral;}_{\mathrm{outside}\left(c\right)}{|I(x,y)-c_2|}^2\mathrm{dxdy}]-({\&Integral;\&Integral;}_s<\&dtri;I,\stackrel{\&RightArrow;}{n}>\mathrm{da}-{\&Integral;\&Integral;\&Integral;}_{{\mathrm{\&Omega;}}_s}k|\&dtri;I\left|\mathrm{dxdydz}\right)+r{\&Integral;}_c\mathrm{ds}$

该公式由三项构成，第一项为图像区域均值计算；第二项为图像边缘能量计算；第三项为曲面几何约束项。其中，α为图像均值项系数，r为几何约束项系数，均为正常数。对于α与r系数值的确定，根据图像特点的不同，依据经验值选取。

步骤(5)：

本发明采用水平集方法求解能量最小的优化问题，根据已建立的曲面演化模型，优化水平集函数定义为：

$S_t=\left\{\mathrm{\&alpha;}\right[{(I-c_1)}^2-{(I-c_2)}^2]-I_{\mathrm{\&epsiv;\&epsiv;}}-r{H}_I\}\stackrel{\&RightArrow;}{n}$

其中 $I_{\mathrm{\&epsiv;\&epsiv;}}=\&dtri;I-H_I|\&dtri;I|$

对于三维血管图像，把S_t作为高一维函数φ(x，y，z)的零水平集，满足S＝{{x，y，z}|φ(x，y，z)＝0}，相应的水平集函数定义为：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\left\{\mathrm{\&alpha;}\right[{(I_0-c_1)}^2-{(I_0-c_2)}^2]-I_{\mathrm{\&epsiv;\&epsiv;}}-r{H}_I\}\left|\mathrm{\&Delta;\&phi;}\right|$

在具体实施过程中，曲面演化是对水平集函数φ的演化过程。在演化之前，水平集函数φ要进行初始化。本文模型对φ的初始化值并没有特殊的要求，为了提高分割精度，应该将φ的初始化尽量接近最终的血管分割结果。因此，利用计算区域中三维区域增长结果作为初始曲而，并把初始曲面内所包含的像素对应的φ值设为-c，计算区域内的其它点设为c，c为常数。确定初始曲面后，进行演化。其中α，γ都为常数，C₁表示血管区内部的CT值平均值，C₂表示背景区的CT值平均值，I₀表示每个点的CT值。实际计算中，参数α，γ通过经验确定，I_εε采用upwind数值方法进行迭代求解，C₁通过种子点像素周围的均值来代替，C₂可以通过求计算区域内的背景区域的像素均值来确定。

曲面停止演化后，所得到的零水平集即为最终的血管分割结果，迭代过程如图3，血管分割结果如图4与图5。

Claims (5)

1.一种基于几何形变模型的三维肺血管图像分割方法，其特征在于本发明方法由以下五个步骤组成：

(1)肋骨区域分割，并选取区域内关键点，采用三次B样条插值计算初始轮廓；

(2)肺血管区域均值计算；

(3)肺血管边缘能量计算；

(4)基于肺血管区域均值、边缘能量与曲面几何约束信息，计算形变模型总能量；

形变模型总能量函数定义为：

$E=\mathrm{\&alpha;}[{\&Integral;}_{\mathrm{inside}\left(c\right)}{|I(x,y)-c_1|}^2\mathrm{dxdy}+{\&Integral;}_{\mathrm{outside}\left(c\right)}{|I(x,y)-c_2|}^2\mathrm{dxdy}]-({\&Integral;\&Integral;}_s<\&dtri;I,\stackrel{\&RightArrow;}{n}>\mathrm{da}-{\&Integral;\&Integral;\&Integral;}_{{\mathrm{\&Omega;}}_s}-H_I|\&dtri;I\left|\mathrm{dxdydz}\right)+r{\&Integral;}_c\mathrm{ds}$

该公式由三项构成，第一项为肺血管图像区域均值计算；第二项为肺血管图像边缘能量计算；第三项为曲面几何约束项：其中，I(x，y)为原始图像，C为活动轮廓线，C₁和C₂分别表示轮廓线内外区域的平均灰度值，α为图像均值项系数，r为几何约束项系数，s为演化曲面，Ω_s表示曲面演化区域，da为曲面区域，

为曲面单位方向矢量，H_I为曲面曲率；

(5)采用基于肺血管影像的水平集演化方法求解形变模型能量最小化。

2.根据权利要求1所述的一种基于几何形变模型的三维肺血管图像分割方法，其特征在于步骤(1)由以下步骤组成：

(1)肋骨分割

由于肋骨CT值较高，与周围组织对比明显，可以采用阈值分割方法获得，阈值T可以依据肋骨CT值范围选取；

(2)初始分割区域计算

在肋骨分割基础上，肋骨区域靠近图像中心近的像素点被作为肋骨关键点P_i，i＝1，2…N保留，然后采用三次B样条插值计算连续轮廓，插值公式如下：

$P_{i,n}\left(t\right)=\underset{k=0}{\overset{n}{\mathrm{\&Sigma;}}}{P}_{i+k}\&CenterDot;F_{k,n}\left(t\right)$

$F_{k,n}\left(t\right)=\frac{1}{n!}\underset{j=0}{\overset{n-k}{\mathrm{\&Sigma;}}}{(-1)}^j\&CenterDot;C_{n+1}^j\&CenterDot;{(t+n-k-j)}^n$

其中，0≤t≤1，P_i+k代表已有离散点，P_i，n代表插值后的连续点，n代表次数。

3.根据权利要求1所述的一种基于几何形变模型的三维肺血管图像分割方法，其特征在于步骤(2)由以下步骤组成：

设原图像I(x，y)被活动轮廓线C划分为目标和背景两个区域，各个区域的平均灰度值为C₁和C₂，肺血管区域均值函数定义为：

$E_{\mathrm{REG}}=\mathrm{\&alpha;}[{\&Integral;}_{\mathrm{inside}\left(c\right)}{|I(x,y)-c_1|}^2\mathrm{dxdy}+{\&Integral;}_{\mathrm{outside}\left(c\right)}{|I(x,y)-c_2|}^2\mathrm{dxdy}]+r{\&Integral;}_c\mathrm{ds}$

其中第一项为图像区域能量拟合项，第二项为几何约束项，C表示活动轮廓线，要求曲线C尽可能短，α为图像均值项系数，r为几何约束项系数，它们依据经验值设定，ds为弧长微元。

4.根据权利要求1所述的一种基于几何形变模型的三维肺血管图像分割方法，其特征在于步骤(3)由以下步骤组成：

设I_εε为在图像梯度方向的二阶导数，I_ηη为其垂直矢量，肺血管边缘能量函数定义为：

$E_{\mathrm{EDGE}}={\&Integral;\&Integral;}_s<\&dtri;I,\stackrel{\&RightArrow;}{n}>\mathrm{da}-{\&Integral;\&Integral;\&Integral;}_{{\mathrm{\&Omega;}}_s}{H}_I|\&dtri;I|\mathrm{dxdydz}$

其中H_I为曲面的曲率，s为演化曲面，Ω_s表示曲面演化区域，da为曲面表面区域，

为曲面单位方向矢量。

5.根据权利要求1所述的一种基于几何形变模型的三维肺血管图像分割方法，其特征在于步骤(5)由以下步骤组成：

本发明采用水平集方法求解能量最小的优化问题，根据已建立的曲面演化模型，优化水平集函数S_t定义为：

$S_t=\left\{\mathrm{\&alpha;}\right[{(I-c_1)}^2{-(I-c_2)}^2]-I_{\mathrm{\&epsiv;\&epsiv;}}-r{H}_I\}\stackrel{\&RightArrow;}{n}$

其中 $I_{\mathrm{\&epsiv;\&epsiv;}}=\&dtri;I-H_I|\&dtri;I|$

对于三维肺血管图像，把S_t作为高一维函数φ(x，y，z)的零水平集，满足S＝{{x，y，z}|φ(x，y，z)＝0}，相应的水平集函数φ定义为：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\left\{\mathrm{\&alpha;}\right[{(I_0-c_1)}^2-{(I_0-c_2)}^2]-I_{\mathrm{\&epsiv;\&epsiv;}}-r{H}_I\}\left|\mathrm{\&Delta;\&phi;}\right|$

其中，I表示原始图像，I_εε表示在图像梯度方向的二阶导数，C₁和C₂分别表示轮廓线内外区域的平均灰度值，α为图像均值项系数，r为几何约束项系数，H_I为曲面曲率，

为曲面单位方向矢量，I₀为像素CT值。

Images