自动初始化的局域活动轮廓模型心脑血管分割方法

摘要

本发明提供了一种自动初始化的局域活动轮廓模型心脑血管分割方法，包括下列步骤：根据医学影像数据中血管的管状特征对三维血管进行各向异性滤波得到初始轮廓线；构造局域能量函数；通过梯度下降法最小化局域能量函数，获得对应的梯度下降流；然后，依据医学影像数据中的血管形状度量构造血管矢量场；最后联合梯度下降流和血管矢量场，得到血管演化方程。本发明根据血管形状自动设置初始轮廓线；分割精确，效率高，鲁棒性强。

说明书

技术领域

本发明涉及计算机应用技术领域，特别涉及一种利用局域活动轮廓模型和血管矢量场进 行血管分割的方法，主要用于医学教学和计算法辅助诊疗等领域。

背景技术

当前，心脑血管疾病已成为人类健康的“第一杀手”。全世界每年死于心脑血管疾病的人 数高达1500万人，中国心脑血管疾病患者已超过2.7亿，每年死于心脑血管疾病近300万人， 占总死亡数量的51%。而幸存下来的患者75%不同程度丧失劳动能力，40%重残。

心脑血管疾病主要是由血管狭窄、动脉瘤和动静脉畸形引起。获取血管的三维形态信息 对于该类疾病的辅助诊疗具有重要作用，如组织的定量分析、计算机辅助诊断、病灶定位、 解剖结构研究及计算机引导手术等。为此，通常需要对医学影像数据集进行分割，以获得精 确的血管组织信息。心脑血管相比其它血管，有其自身特点，如拓扑结构的复杂性，灰度的 动态变化性及部分血管结构的细小微弱性等，这些都增加了心脑血管分割的难度。基于血管 造影图像的心脑血管分割已成为国际上计算机科学、信息学、生物医学工程等相关领域的研 究热点。

近年来，基于几何活动轮廓模型在图像分析中获得了广泛的应用，由于同一类的目标形 状上具有一定的相似性，通过在一定条件下控制变形得到另一个目标形状，从而进行分类、 形似度比较等分析。如果把形变面放入图像中，由图像的各种特征，如灰度、形状、梯度等， 控制目标的形变，使之形状向图像中的目标变化，从而提取出图像的特定目标。Farag等人提 出了基于水平集方法的血管树三维分割方法，实现了从MRA数据集中抽取3D脑血管，该方 法属于有监督的分类方法，需要事先知道分类的类数及每类的概率分布。通常，不同器官、 不同数据的概率分布模型建立比较困难。Gao等人提出一种全自动的TOF-MRA图像脑血管 分割方法。首先构造有限混合模型，估计多个分量的分布参数，然后在原始数据空间中对脑 血管进行粗分割；接下来，借助统计信息通过对能量函数中的速度项和区域项进行修改，从 而实现脑血管的精细分割。但上述算法都没有结合血管自身的几何特征。Hernandez和Frangi 提出了非参数测地线活动区域模型的脑动脉分割方法。在图像的高阶特征空间中，应用K-N 近邻非参数估计方法对血管几何特征的联合分布进行估计，然后将分布信息嵌入能量函数中， 为了对参数进行估计，K-N近邻法需要一定数量的样本进行学习。针对血管组织灰度分布不 均的问题，尤其是细小血管，难以准确分割。当前方法在轮廓线演化过程中，通常图像域中 所有像素参与运算，运算量大。此外，传统的基于几何活动轮廓模型的血管分割方法容易陷 入局部极小值，对初始轮廓线的设置要求较为严格。

我们首先将几何活动轮廓模型局域化，然后借助构造的血管矢量场对血管进行精确分割。 该方法的特点在于：对血管造影体数集进行血管滤波增强，获得血管的大致轮廓，以此作为 初始轮廓线，在血管矢量场的辅助下，通过局域能量的最小化实现血管的精确分割；无需手 工设置轮廓线，自动化程度高；分割精确，效率高，鲁棒性强。

发明内容

本发明的目的就是克服现有技术中的不足，提供一种自动设置初始轮廓线的，高效、精 确、鲁棒性强的局域活动轮廓模型心脑血管分割方法。

为实现上述目的，本发明的技术方案是利用血管的形状特征设置初始轮廓线，局域化几 何活动轮廓模型，借助血管矢量场实现血管组织的精确分割。

一种自动初始化的局域活动轮廓模型心脑血管分割方法，包括下列步骤：

（1）血管形状滤波及初始轮廓线设置：根据医学影像数据中血管的管状特征，借助 Hessian矩阵，获取血管中心线方向和横切面方向，对三维血管进行各向异性滤波；通过阈值 法和边缘检测算子提取血管轮廓作为初始轮廓线；

（2）构造局域能量函数：定义活动轮廓线Γ:{x|φ(x)=0}上各点的小邻域O_k；通过局域 加权函数对活动轮廓线的点进行加权拟合，同时结合加权因子构造局域能量函数：

$E={\int }_{{\Omega }_x}\delta \left(\phi \left(x\right)\right){\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·{(H\left(\phi \left(y\right)\right)I\left(y\right)-f_{\mathrm{in}}\left(x\right))}^2$

$\left(4\right)$

$+(1-H\left(\phi \left(y\right)\right)){(I\left(y\right)-f_{\mathrm{out}}\left(x\right))}^2)\mathrm{dydx}$

其中，x为活动轮廓线Γ:{x |φ(x)=0}上一点，以点x为中心定义一小邻域O_k，δ(φ(x)) 为函数H(x)的一阶导数，H(x)用于指定小邻域O_k与活动轮廓线相交的内部区域；

（3）获得梯度下降流：通过梯度下降法最小化局域能量函数，获得对应的梯度下降流

（4）构造血管矢量场：依据医学影像数据中的血管形状度量构造血管矢量场： $\overrightarrow{V}\left(x\right)=\left\{{\overrightarrow{V}}_{\sigma }\left(x\right)\right|{\alpha }_{\sigma }{R}_{\sigma }\left(x\right)=R\left(x\right)\},$其中：$R\left(x\right)=\underset{{\sigma }_{\min }\le \sigma \le {\sigma }_{\max }}{\max }\left\{{\alpha }_{\sigma }{R}_{\sigma }\left(x\right)\right\},$α_σ为尺度权重，

所述血管矢量场的方向平行于血管方向，大小是和血管相关的函数；

（5）得到血管演化方程：联合步骤（3）的梯度下降流和步骤（4）构造的血管矢量场， 得到血管演化方程：$\frac{\partial \phi }{\partial t}\left(x\right)=f_g^{\mathrm{LRF}}+f\left(R\left(x\right)\right)|\overrightarrow{V}\left(x\right)·▿\phi \left(x\right)|,$其中表示梯度下降流，f(R(x)) 是一表示血管矢量场大小的、关于R(x)的函数。

优选地，步骤（1）中所述的各向异性滤波是基于血管形状特征的。

优选地，所述步骤（2）的具体步骤为：

①设表示图像域，I:Ω→R表示给定的图像，x为活动轮廓线上一点，以点x为 中心定义一小邻域O_k，活动轮廓线Γ与小邻域O_k相交的内部区域均值μ_in(x)由内部区域点的 灰度I(y)及与x距离的加权获得，该加权值由引入的高斯核函数K_σ(x,y)得到；外部区域均 值μ_out(x)也通过高斯核函数K_σ(x,y)的方法得到：

${\mu }_{\mathrm{in}}\left(x\right)=\frac{{\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·H\left(\phi \left(y\right)\right)·I\left(y\right)\mathrm{dy}}{{\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·H\left(\phi \left(y\right)\right)\mathrm{dy}}---\left(1\right)$

${\mu }_{\mathrm{out}}\left(x\right)=\frac{{\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·(1-H\left(\phi \left(y\right)\right))·I\left(y\right)\mathrm{dy}}{{\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·(1-H\left(\phi \left(y\right)\right))\mathrm{dy}}---\left(2\right)$

H(x)用于指定O_k与Γ相交的内部区域，φ(y)是I(y)对应的水平集函数；

②将点x∈Γ的能量函数定义为由小邻域内的所有点的灰度值与 μ_in(x)和μ_out(x)的差通过加权平均得到：

$E_x^{\mathrm{LRF}}(\Gamma ,f_{\mathrm{in}}\left(x\right),f_{\mathrm{out}}\left(x\right))={\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·(H{\left(\phi \left(y\right)\right)(I\left(y\right)-f_{\mathrm{in}}\left(x\right))}^2$

$\left(3\right)$

$+(1-H\left(\phi \left(y\right)\right)){(I\left(y\right)-f_{\mathrm{out}}\left(x\right))}^2)\mathrm{dy}$

进而得到局部能量函数。

优选地，所述步骤（4）的具体步骤为：

a．通过对医学影像获得的3D灰度图像的Hessian矩阵H_I的特征分析得到局部血管的方 向和管状结构度量R：Hessian矩阵H_I为：

$H_I=\left(\begin{array}{ccc}{I}_{\mathrm{xx}}& I_{\mathrm{xy}}& I_{\mathrm{xz}}\\ I_{\mathrm{yx}}& I_{\mathrm{yy}}& I_{\mathrm{yz}}\\ I_{\mathrm{zx}}& I_{\mathrm{zy}}& I_{\mathrm{zz}}\end{array}\right)---\left(5\right)$

λ₁、λ₂和λ₃是H_I的特征值，设|λ₁|≤|λ₂|≤|λ₃|，定义：和三个相应的特征向量 相互垂直形成局部笛卡尔坐标，表示灰度变化量最小的方向为血管方向，和形成垂直 于血管的纵切面，血管矢量场为在血管区域的特征向量的集合：

τ为一设定阈值0.05，R(x)为管状结构度量，用于滤掉非血管区域的向量，

式（7）中，A=|λ₂|/|λ₃|，c是一个小的常数0.7； 对于每一像素点，存在两个相反的血管方向，即和为了驱使活动轮廓演化到血管内部， 矢量场应该和活动轮廓线的法向量方向一致，于是有为此，血管矢量场定 义为：

$\overrightarrow{V}\left(x\right)=\left(\begin{array}{ccc}\overrightarrow{V}\left(x\right)& \mathrm{if}& \overrightarrow{V}\left(x\right)·▿\phi \left(x\right)\ge 0\\ -\overrightarrow{V}\left(x\right)& \mathrm{if}& \overrightarrow{V}\left(x\right)·▿\phi \left(x\right)<0\end{array}\right)---\left(8\right)$

b．管状结构度量取加权的所有尺度下最大的度量：

$R\left(x\right)=\underset{{\sigma }_{\min }\le \sigma \le {\sigma }_{\max }}{\max }\left\{{\alpha }_{\sigma }{R}_{\sigma }\left(x\right)\right\},$α_σ＝1.6；

在多尺度框架下计算血管矢量场：$\overrightarrow{V}\left(x\right)=\left\{{\overrightarrow{V}}_{\sigma }\left(x\right)\right|{\alpha }_{\sigma }{R}_{\sigma }\left(x\right)=R\left(x\right)\}.$

本发明基于血管形状的局域活动轮廓模型心脑血管分割方法，将血管造影图像进行各向 异性滤波，在有效增强血管的对比度的同时也保留了血管边缘细节，进而得到血管的大致轮 廓作为初始轮廓线，采用局域化Chan-Vese活动轮廓模型，结合加权因子函数，仅由轮廓线 上的点构造局域能量函数，通过对能量函数最小化，得到其梯度下降流；定义血管矢量场， 引导活动轮廓线沿着血管方向演化。本方法自动设置初始轮廓线；分割精确，效率高，鲁棒 性强。

附图说明

图1为本发明所述基于血管形状的局域活动轮廓模型心脑血管分割方法的流程图。

具体实施方式

下面结合附图对本发明作进一步详细的说明，下面通过附图将本发明基于血管形状的局 域活动轮廓模型心脑血管分割方法进行详尽的描述。

由于医学影像数据的血管与周围组织的灰度对比度较小，高斯滤波等传统方法难以将血 管有效增强，本发明结合心脑血管的形状特征有针对性地增强对比度；能量函数的构造仅通 过活动轮廓线及周围小邻域的点实现，提高演化效率；依据血管中心线位置和活动轮廓线的 方向定义血管矢量场，辅助活动轮廓线沿着血管方向演化。本发明所述基于血管形状的局域 活动轮廓模型心脑血管分割方法的总体流程如图1所示，具体包括以下步骤：

步骤S1011：血管滤波增强及初始轮廓线设置：

借助体数据场的Hessian矩阵，从心脑血管的医学影像数据中获取血管中心线方向和横 切面方向，在不同方向对三维血管进行各项异性滤波，在保留血管边缘细节的同时，又能有 效增强血管的对比度；然后，通过阈值法和边缘检测算子得到大致的血管轮廓作为活动轮廓 模型的初始轮廓线。所述心脑血管的医学影像数据为现有技术中常见的医学图像，例如序列 CTA或MRA断层切片。

步骤S1012：构造局域能量函数：

Chan和Vese提出Chan-Vese模型是基于图像全局统计信息的，尽管能够达到全局最优， 但难以分割灰度分布不均匀的血管图像。为此，本发明对Chan-Vese模型局域化，仅考虑活 动轮廓线周围的点。在步骤S1011得到的初始轮廓线的保证下，能精确演化至血管的真实边 缘。活动轮廓线从初始轮廓线开始，在能量函数能量最小化的过程中，不断进行演化，当能 量极小时，演化结束，此时的轮廓线对应的就是血管的边缘。其步骤为：

①设表示图像域，I:Ω→R表示给定的图像，选取活动轮廓线Γ:{x |φ(x)=0}上 一点x，以点x为中心定义一小邻域O_k。活动轮廓线Γ与O_k相交的内部区域均值μ_in(x)由内 部区域点的灰度I(y)及与x距离的加权获得，该加权值由引入的高斯核函数K_σ(x,y)得到； 外部区域均值μ_out(x)也通过高斯核函数K_σ(x,y)的方法得到：

${\mu }_{\mathrm{in}}\left(x\right)=\frac{{\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·H\left(\phi \left(y\right)\right)·I\left(y\right)\mathrm{dy}}{{\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·H\left(\phi \left(y\right)\right)\mathrm{dy}}---\left(1\right)$

${\mu }_{\mathrm{out}}\left(x\right)=\frac{{\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·(1-H\left(\phi \left(y\right)\right))·I\left(y\right)\mathrm{dy}}{{\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·(1-H\left(\phi \left(y\right)\right))\mathrm{dy}}---\left(2\right)$

H(x)用于指定O_k与Γ相交的内部区域，φ(y)是I(y)对应的水平集函数；

②对于活动轮廓线Γ上的一点x∈Γ，将该点x∈Γ的能量函数定义为 由小邻域O_k内所有点的灰度值与μ_in(x)和μ_out(x)的差通过加权平均得 到；不仅考虑了小邻域内点I(y)的灰度，而且考虑距离点x的距离K_σ(x,y)，即：

$E_x^{\mathrm{LRF}}(\Gamma ,f_{\mathrm{in}}\left(x\right),f_{\mathrm{out}}\left(x\right))={\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·(H{\left(\phi \left(y\right)\right)(I\left(y\right)-f_{\mathrm{in}}\left(x\right))}^2$

$\left(3\right)$

$+(1-H\left(\phi \left(y\right)\right)){(I\left(y\right)-f_{\mathrm{out}}\left(x\right))}^2)\mathrm{dy}$

进而得到活动轮廓线Γ上所有点的局部能量函数：

$E={\int }_{{\Omega }_x}\delta \left(\phi \left(x\right)\right){\int }_{{\Omega }_y}{K}_{\sigma }(x,y)·{(H\left(\phi \left(y\right)\right)I\left(y\right)-f_{\mathrm{in}}\left(x\right))}^2$

$\left(4\right)$

$+(1-H\left(\phi \left(y\right)\right)){(I\left(y\right)-f_{\mathrm{out}}\left(x\right))}^2)\mathrm{dydx}$

式（4）中，δ(φ(x))为函数H(x)的一阶导数；

步骤S1013：获得梯度下降流：

借助梯度下降法，最小化活动轮廓线Γ上所有点的局部能量函数，即最小化式(4)所述的 能量函数，获得对应的梯度下降流

步骤S102：构造血管矢量场：

为了对细小血管精确分割，构造一血管矢量场，驱使活动轮廓线在血管内部以较高的速 度沿着血管方向运动。矢量场的方向平行于血管方向，大小是和血管相关的函数。分别通过 如下步骤实现：

a．通过对医学影像获得的3D灰度图像的Hessian矩阵H_I的特征分析得到局部血管的方 向和管状结构度量R。Hessian矩阵H_I为：

$H_I=\left(\begin{array}{ccc}{I}_{\mathrm{xx}}& I_{\mathrm{xy}}& I_{\mathrm{xz}}\\ I_{\mathrm{yx}}& I_{\mathrm{yy}}& I_{\mathrm{yz}}\\ I_{\mathrm{zx}}& I_{\mathrm{zy}}& I_{\mathrm{zz}}\end{array}\right)---\left(5\right)$

λ₁、λ₂和λ₃是H_I的特征值，设|λ₁|≤|λ₂|≤|λ₃|，和是相应的特征向量，和 三个特征向量相互垂直形成局部笛卡尔坐标。对于血管结构，表示灰度变化量最小的方 向，也就是血管方向；和形成垂直于血管的纵切面。血管矢量场就是在血管区域的特征 向量的集合，即：

R(x)为管状结构度量，用于滤掉非血管区域的向量；τ为一设定阈值，取τ=0.05。本 发明采用Manniesing提出的管状结构度量：

式（7）中，A＝|λ₂|/|λ₃|，c是一个小的常数，取c=0.7。

b．由于管状结构是一个与尺度相关的概念，不同直径的管状结构，只有在特定的尺度下 表现出管状的特征，因此，需要在多尺度框架下计算血管矢量场，最终的管状结构度量取加 权的所有尺度下最大的度量。即，血管矢量场为：

$\overrightarrow{V}\left(x\right)=\left\{{\overrightarrow{V}}_{\sigma }\left(x\right)\right|{\alpha }_{\sigma }{R}_{\sigma }\left(x\right)=R\left(x\right)\},$

其中：$R\left(x\right)=\underset{{\sigma }_{\min }\le \sigma \le {\sigma }_{\max }}{\max }\left\{{\alpha }_{\sigma }{R}_{\sigma }\left(x\right)\right\},$α_σ为尺度权重，α_σ＝1.6。

对于每个像素，存在两个相反的血管方向，即和为了驱使活动轮廓演化到血管内 部，矢量场应该和活动轮廓线的法向量方向一致，于是有最终血管矢量场 定义为：

$\overrightarrow{V}\left(x\right)=\left(\begin{array}{ccc}\overrightarrow{V}\left(x\right)& \mathrm{if}& \overrightarrow{V}\left(x\right)·▿\phi \left(x\right)\ge 0\\ -\overrightarrow{V}\left(x\right)& \mathrm{if}& \overrightarrow{V}\left(x\right)·▿\phi \left(x\right)<0\end{array}\right)---\left(8\right)$

步骤S103：构造血管演化方程：

根据步骤S1013获得的梯度下降流和步骤S102构建的血管矢量场，得到血管演化 方程：$\frac{\partial \phi }{\partial t}\left(x\right)=f_g^{\mathrm{LRF}}+f\left(R\left(x\right)\right)|\overrightarrow{V}\left(x\right)·▿\phi \left(x\right)|$

其中，f(R(x))是一表示血管矢量场大小的、关于R(x)的函数。

总之，本发明的实施例公布的是其较佳的实施方式，但并不限于此。本领域的普通技术 人员极易根据上述实施例，领会本发明的精神，并做出不同的引申和变化，但只要不脱离本 发明的精神，都在本发明的保护范围之内。