CN103455999A - 基于内超声图像序列的血管壁边缘自动检测方法 - Google Patents

Abstract

一种基于内超声图像序列的血管壁边缘自动提取方法，先粗检：利用血管壁处的信息特征，通过系列图像处理方法获取第一帧图像的初始边缘；再细检：改进GVF-snake算法，引入自我调节因子和自适应法向外力，在增强边缘附近数据项的梯度影响，减少扩散引起的平滑效应的同时，能依据图像边缘和当前轮廓曲线相对位置来调整力的方向，扩大活动轮廓的捕捉范围，解决了由于无法到达局部区域而不能获得所需边缘的问题，使初始边缘精确收敛于实际管壁边缘。在变形过程中加入三次B样条，以减少控制点，提高收敛速度，光顺轮廓曲线，得到第一帧图像的最终边缘。由于相邻图像之间时间、空间相关性较大，将所获最终边缘作为下一帧图像初始边缘，重复上述细检工作，以此类推，得到序列图像的血管壁边缘。

Description

基于内超声图像序列的血管壁边缘自动检测方法

技术领域

本发明涉及一种基于内超声图像序列的血管边缘自动检测方法，属于计算机医学图像分析领域，本方法实现了血管壁内外膜的自动提取，以此为基础进行图像分析，可以对血管腔的面积、血管厚度进行精确测量，有利于发现或识别早期冠状动脉粥样硬化斑块，具有重要的临床意义。

背景技术

在世界范围内，心脑血管疾病已成为人类疾病死亡的首要原因。血管内超声(Intravascular Ultrasound，IVUS)是当前诊断心脑血管疾病的最新技术。提取血管内超声超声图像中血管壁的内外边缘，并以此为基础进行图像分析，对于发现或识别早期冠状动脉粥样硬化斑块具有重要的临床意义。

近年来，基于血管内超声图像的血管边缘检测算法得到广泛研究。其中比较常用的提取算法有水平集模型，动态规划法，活动轮廓模型(Snake模型)及其各种改进型算法。这些算法需要手动选取初始轮廓，抗噪声干扰能力差，且存在边界凹陷处不易收敛的缺点。为了改善闭合边界凹陷处的曲线收敛特性，引入梯度向量流场(gradient vector flow，GVF)为新的外力来定义的GVF-Snake模型，依然存在一些不足，比如：收敛速度慢，且对初始边缘位置依赖大，当血管内超声图像有较大的噪声和伪影干扰时，提取效果不佳。

发明内容

本发明的目的在于，通过提供基于内超声图像的血管边缘自动检测方法，以提高边缘提取的准确率和运行效率。

本发明是采用以下技术手段实现的：

一种基于内超声图像序列的血管壁边缘自动检测方法，利用图像序列的时间特性和经验知识来抑制噪声，并检测冠状动脉血管壁的初始边缘，然后，利用三次B样条来光顺轮廓；再结合灰度图像的梯度、均值、方差信息，结合自适应算法，在梯度矢量流轮廓提取算法中引入新的调节因子和自我调整法向外力，用于对初始轮廓进行变形；最后获得冠状动脉血管内超声图像血管壁的最终边缘；包括下述步骤：

1.1.根据血管内超声图像的特点，综合利用图像处理方法，分别对血管壁内外边缘进行预提取，获得血管壁的初始边缘；

1.2.通过三次B样条曲线的局部性，构造矩阵如下：

$P(s-n)=\frac{1}{6}\left[\begin{array}{cccc}{s}^3& s^2& s& 1\end{array}\right]\left[\begin{array}{cccc}-1& 3& -3& 1\\ 3& -6& 3& 0\\ -3& 0& 3& 0\\ 1& 4& 1& 0\end{array}\right]\left[\begin{array}{c}{C}_{n-1}\\ C_n\\ C_{n+1}\\ C_{n+2}\end{array}\right]---\left(1\right)$

其中0≤s≤1，在初始边缘上通过径向搜索得到均匀采样点P_i(s)＝(x_i，y_i)，经改进的GVF-snake变形后到达的点的位置坐标为P_j(s)＝(x_j，y_j)，采样点位置坐标的变化与所在改进后的矢量场的关系满足下式：

P_j(s)-P_i(s)＝μV_GVF(x_i，y_i)           (2)

μ为比例系数。通过B样条的光顺特性用来完成整个闭合轮廓曲线平滑程度的改善；

1.3.设置参数λ₁和λ₂来调节整个梯度矢量流场V_GVF(x，y)＝[u(x，y)，v(x，y)]，该矢量场的外力方程式为：

其中，λ₁和λ₂的大小可以通过迭代过程中的循环变量i来控制，

$\left\{\begin{array}{c}{l}_1=m*i/k\\ l_2=1-l_1=1-m*i/k\end{array}\right.---\left(4\right)$

k是灰度图像I(x，y)中x和y的最大值，将矢量场外力方程中的u和v作为时间函数进行求解：

对上式进行中心差分，得到离散方程：

其中，n为迭代次数，Δt为时间间隔；

1.4.在GVF基础上添加活动轮廓上控制点处的单位法向量，用来产生一个垂直于活动轮廓曲线的自适应外力：

$F_N=\mathrm{sign}\left(\right[u,v]\×\stackrel{\→}{n})\frac{1}{1+k{\left(\right[u,v]\×\stackrel{\→}{n})}^2}---\left(7\right)$

其中，为活动轮廓上控制点处所取外力的符号，

是单位法向量，k为梯度矢量和法向矢量内积的加权系数；通过不断迭代变形，获得所需的最终边缘。

前述的对血管壁外边缘进行预提取，获得血管壁的初始边缘的步骤如下：

2.1.选取连续10帧图像，计算其时间方差图，利用帧的连续性来移除随时间不变性而遗留下来的一些噪声(如环晕伪像等)和连接边缘处的一些断点(如声影区域)：

$V(x,y)=\frac{1}{n-1}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{[I_m(x,y)-\stackrel{\‾}{I}(x,y)]}^2---\left(8\right)$

其中，I_m(x，y)是第m帧灰度图像，

是连续n帧的平均灰度图像；

2.2.对灰度图像V(x，y)进行高斯滤波来减少血液散射引起的噪声，并对降噪后的图像进行阈值分割，获得二值图像，

V_g(x，y)＝G_σ(x，y)V(x，y)          (9)

式中，G_σ(x，y)是一个标准差为σ的二维高斯函数；

2.3.利用连通分量标记法，找出二值图像区域中所有的闭合轮廓，计算轮廓周长，设定周长阈值为L，去除一些来自环晕伪像和导丝伪像等产生的小面积区域干扰；

2.4.设所需外边缘点数为N_out，图像中心点为(x₀，y₀)，以Δθ＝2π/N_out间隔对二值图像做径向搜索，那么在径向上值为1的点即为边缘点(r_i，θ_i)，i＝1～N_out，通过边缘点即可得到所需的初始外边缘。

步骤1.1所述的对血管壁内边缘进行预提取，获得血管壁的初始边缘的步骤如下：

3.1.对血管内超声图像做直方图均衡处理，为了进一步拉伸图像的灰度级，引入基于灰度图像的指数变换：

g(x，y)＝b^{c|f(x，y)-a|}-1              (10)

f(x，y)和g(x，y)分别为指数变换前后的图像，参数a改变曲线的起始位置，参数c改变曲线的变化速率；

3.2.通过数学形态学处理将图像中存在的空洞填满以及将毛刺滤掉，引用开启操作消除与结构元素相比尺寸较小的亮细节而保持图像整体灰度值和大的亮的区域基本不受影响，进行阈值分割，获得二值图像；

3.3.设所需内边缘点数为N_int，图像中心点为(x₀，y₀)，以Δθ＝2π/N_int间隔对二值图像做径向搜索，那么在径向上值为1的点即为边缘点(r_i，θ_i)，i＝1～N_int，通过边缘点即可得到所需的初始内边缘。

前述的外边缘点数N_out为36，图像帧数n取为10。

前述的内边缘点数N_int为36。

本发明采用改进的GVF-Snake方法对血管内超声进行边缘自动检测，与现有方法相比，其优点在于：

依据血管壁内外膜处的图像特征，通过不同的图像处理方法来获取初始边缘；采用三次B样条来减少控制点，提高收敛速度，光顺轮廓曲线。在此基础上，利用改进的GVF-snake算法，使初始边缘精确收敛于实际管壁边缘。首先在GVF-snake的基础上引入自我调整调节因子，增强边缘附近数据项的梯度影响，减少由于扩散引起的平滑效应；然后给活动轮廓添加一个自适应外力，依据图像边缘和当前轮廓曲线的相对位置来调整力的方向。该外力的加入增大了梯度向量场的驱动力，在强噪声下也能精确地进行弱边缘图像的轮廓逼近，解决了由于无法到达局部区域而不能获得所需边缘的问题。

附图说明

图1是连续10帧血管内超声图像；

图2(a)是连通分量标记技术示意图；

图2(b)是径向搜索示意图；

图2(c)是初始外边缘获取结果示意图；

图3(a)是指数拉伸示意图；

图3(b)是形态学处理结果示意图；

图3(c)是初始内边缘获取结果示意图；

图4(a)是血管内超声图像血管壁边缘预提取结果示意图；

图4(b)是血管内超声图像血管壁最终边缘提取结果示意图；

图5是图像序列边缘提取结果显示；

其中，a为第二帧图像，b为第三帧图像，c为第四帧图像，d为第五帧图像，e为第六帧图像，f为第七帧图像，g为第八帧图像，h为第九帧图像，i为第十帧图像；

图6是本发明的流程示意图。

具体实施方式

一种基于血管内超声边缘自动检测方法，具体步骤如下：

步骤1、根据血管内超声图像的特点，综合利用图像处理方法，对血管壁外边缘进行预提取，获得血管壁的初始边缘。

步骤1.1、如图1所示，选取连续10帧图像，计算其时间方差图，利用帧的连续性来移除随时间不变性而遗留下来的一些噪声(如环晕伪像等)和连接边缘处的一些断点(如声影区域)：

$V(x,y)=\frac{1}{n-1}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{[I_m(x,y)-\stackrel{\‾}{I}(x,y)]}^2---\left(6\right)$

其中，I_m(x，y)是第m帧灰度图像，

是连续n帧的平均灰度图像。

步骤1.2、对灰度图像V(x，y)进行高斯滤波来减少血液散射引起的噪声，并对降噪后的图像进行阈值分割，获得二值图像。

V_g(x，y)＝G_σ(x，y)V(x，y)       (7)

式中，G_σ(x，y)是一个标准差为σ的二维高斯函数。

步骤1.3、利用连通分量标记法，找出二值图像区域中所有的闭合轮廓，计算轮廓周长，设定周长阈值为L，去除如图2(a)所示的一些来自环晕伪像和导丝伪像等产生的小面积区域干扰。

步骤1.4、如图2(b)所示，设所需外边缘点数为N_out，图像中心点为(x₀，y₀)，以Δθ＝2π/N_out间隔对二值图像做径向搜索，那么在径向上值为1的点即为边缘点(r_i，θ_i)，i＝1～N_out，通过边缘点即可得到所需的初始外边缘，如图2(c)所示。

步骤2、根据血管内超声图像的特点，综合利用图像处理方法，对血管壁内边缘进行预提取，获得血管壁的初始边缘。

步骤2.1、对血管内超声图像做直方图均衡处理，为了进一步拉伸图像的灰度级，引入基于灰度图像的指数变换：

g(x，y)＝b^c|f(x，y)|-1             (8)

f(x，y)和g(x，y)分别为指数变换前后的图像，参数a，b和c是曲线形状位置的控制参数，其中，a为曲线的左右偏移量，b控制曲线弯曲的速率，c控制曲线的弯曲程度，所得结果如图3(a)所示。

步骤2.2、通过数学形态学处理将图像中存在的空洞填满以及将毛刺滤掉，引用开启操作消除与结构元素相比尺寸较小的亮细节而保持图像整体灰度值和大的亮的区域基本不受影响，如图3(b)所示，再进行阈值分割，获得二值图像。

步骤2.3、设所需内边缘点数为N_int，图像中心点为(x₀，y₀)，以Δθ＝2π/N_int间隔对二值图像做径向搜索，那么在径向上值为1的点即为边缘点(r_i，θ_i)，i＝1～N_int，通过边缘点即可得到所需的初始内边缘，如图3(c)所示。由此可得血管内超声图像血管壁内外膜边缘预提取结果如图4(a)所示。

步骤3、通过三次B样条曲线的局部性，构造矩阵如下：

$P(s-n)=\frac{1}{6}\left[\begin{array}{cccc}{s}^3& s^2& s& 1\end{array}\right]\left[\begin{array}{cccc}-1& 3& -3& 1\\ 3& -6& 3& 0\\ -3& 0& 3& 0\\ 1& 4& 1& 0\end{array}\right]\left[\begin{array}{c}{C}_{n-1}\\ C_n\\ C_{n+1}\\ C_{n+2}\end{array}\right]---\left(9\right)$

其中0≤s≤1，在初始边缘上通过径向搜索得到均匀采样点P_i(s)＝(x_i，y_i)，经改进的GVF-snake变形后到达的点的位置坐标为P_j(s)＝(x_j，y_j)，采样点位置坐标的变化与所在改进后的矢量场的关系满足下式：

P_j(s)-P_i(s)＝μV_GVF(x_i，y_i)       (10)

μ为比例系数。通过B样条的光顺特性用来完成整个闭合轮廓曲线平滑程度的改善。

步骤4、设置参数λ₁和λ₂来调节整个梯度矢量流场V_GVF(x，y)＝[u(x，y)，v(x，y)]。该矢量场的外力方程式为：

其中，λ₁和λ₂的大小可以通过迭代过程中的循环变量i来控制。

$\left\{\begin{array}{c}{l}_1=m*i/k\\ l_2=1-l_1=1-m*i/k\end{array}\right.---\left(12\right)$

k是灰度图像I(x，y)中x和y的最大值。将矢量场外力方程中的u和v作为时间函数进行求解：

对上式进行中心差分，得到离散方程：

其中，n为迭代次数，Δt为时间间隔。

步骤5、在GVF基础上添加活动轮廓上控制点处的单位法向量，用来产生一个垂直于活动轮廓曲线的自适应外力：

$F_N=\mathrm{sign}\left(\right[u,v]\×\stackrel{\→}{n})\frac{1}{1+k{\left(\right[u,v]\×\stackrel{\→}{n})}^2}---\left(15\right)$

其中，

为活动轮廓上控制点处所取外力的符号，

是单位法向量，k为梯度矢量和法向矢量内积的加权系数。通过不断迭代变形，获得所需的最终边缘，如图4(b)所示。

以该帧图像所得最终边缘作为下一帧图像的初始边缘，运用本发明中的形变方法逐步逼近得到最终边缘，以此类推，得到其余九帧图像的边缘检测结果如图5所示。

最后应说明的是：以上实施例仅用以说明本发明而并非限制本发明所描述的技术方案；因此，尽管本说明书参照上述的各个实施例对本发明已进行了详细的说明，但是，本领域的普通技术人员应当理解，仍然可以对本发明进行修改或等同替换；而一切不脱离发明的精神和范围的技术方案及其改进，其均应涵盖在本发明的权利要求范围当中。

Claims (5)

1.一种基于内超声图像序列的血管壁边缘自动检测方法，利用图像序列的时间特性和经验知识来抑制噪声，并检测冠状动脉血管壁的初始边缘，然后，利用三次B样条来光顺轮廓；再结合灰度图像的梯度、均值、方差信息，结合自适应算法，在梯度矢量流轮廓提取算法中引入新的调节因子和自我调整法向外力，用于对初始轮廓进行变形；最后获得冠状动脉血管内超声图像血管壁的最终边缘；其特征在于：包括下述步骤：

1.1.根据血管内超声图像的特点，综合利用图像处理方法，分别对血管壁内外边缘进行预提取，获得血管壁的初始边缘；

1.2.通过三次B样条曲线的局部性，构造矩阵如下：

$P(s-n)=\frac{1}{6}\left[\begin{array}{cccc}{s}^3& s^2& s& 1\end{array}\right]\left[\begin{array}{cccc}-1& 3& -3& 1\\ 3& -6& 3& 0\\ -3& 0& 3& 0\\ 1& 4& 1& 0\end{array}\right]\left[\begin{array}{c}{C}_{n-1}\\ C_n\\ C_{n+1}\\ C_{n+2}\end{array}\right]---\left(1\right)$

其中0≤s≤1，在初始边缘上通过径向搜索得到均匀采样点P_i(s)＝(x_i，y_i)，经改进的GVF-snake变形后到达的点的位置坐标为P_j(s)＝(x_j，y_j)，采样点位置坐标的变化与所在改进后的矢量场的关系满足下式：

P_j(s)-P_i(s)＝μV_GVF(x_i，y_i)      (2)

μ为比例系数。通过B样条的光顺特性用来完成整个闭合轮廓曲线平滑程度的改善；

1.3.设置参数λ₁和λ₂来调节整个梯度矢量流场V_GVF(x，y)＝[u(x，y)，v(x，y)]，该矢量场的外力方程式为：

(3)

其中，λ₁和λ₂的大小可以通过迭代过程中的循环变量i来控制，

$\left\{\begin{array}{c}{l}_1=m*i/k\\ l_2=1-l_1=1-m*i/k\end{array}\right.---\left(4\right)$

k是灰度图像I(x，y)中x和y的最大值，将矢量场外力方程中的u和v作为时间函数进行求解：

对上式进行中心差分，得到离散方程：

其中，n为迭代次数，Δt为时间间隔；

1.4.在GVF基础上添加活动轮廓上控制点处的单位法向量，用来产生一个垂直于活动轮廓曲线的自适应外力：

$F_N=\mathrm{sign}\left(\right[u,v]\×\stackrel{\→}{n})\frac{1}{1+k{\left(\right[u,v]\×\stackrel{\→}{n})}^2}---\left(7\right)$

其中，

为活动轮廓上控制点处所取外力的符号，

是单位法向量，k为梯度矢量和法向矢量内积的加权系数；通过不断迭代变形，获得所需的最终边缘。

2.根据权利要求1所述基于内超声图像序列的血管壁边缘自动检测方法，其特征在于：步骤1.1所述的对血管壁外边缘进行预提取，获得血管壁的初始边缘的步骤如下：

2.1.选取连续10帧图像，计算其时间方差图，利用帧的连续性来移除随时间不变性而遗留下来的一些噪声(如环晕伪像等)和连接边缘处的一些断点(如声影区域)：

$V(x,y)=\frac{1}{n-1}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{[I_m(x,y)-\stackrel{\‾}{I}(x,y)]}^2---\left(8\right)$

其中，I_m(x，y)是第m帧灰度图像，

是连续n帧的平均灰度图像；

2.2.对灰度图像V(x，y)进行高斯滤波来减少血液散射引起的噪声，并对降噪后的图像进行阈值分割，获得二值图像，

V_g(x，y)＝G_σ(x，y)V(x，y)                (9)

式中，G_σ(x，y)是一个标准差为σ的二维高斯函数；

2.3.利用连通分量标记法，找出二值图像区域中所有的闭合轮廓，计算轮廓周长，设定周长阈值为L，去除一些来自环晕伪像和导丝伪像等产生的小面积区域干扰；

2.4.设所需外边缘点数为N_out，图像中心点为(x₀，y₀)，以Δθ＝2π/N_out间隔对二值图像做径向搜索，那么在径向上值为1的点即为边缘点(r_i，θ_i)，i＝1～N_out，通过边缘点即可得到所需的初始外边缘。

3.根据权利要求1所述基于内超声图像序列的血管壁边缘自动检测方法，其特征在于：步骤1.1所述的对血管壁内边缘进行预提取，获得血管壁的初始边缘的步骤如下：

3.1.对血管内超声图像做直方图均衡处理，为了进一步拉伸图像的灰度级，引入基于灰度图像的指数变换：

g(x，y)＝b^{c|f(x，y)-a|}-1              (10)

f(x，y)和g(x，y)分别为指数变换前后的图像，参数a改变曲线的起始位置，参数c改变曲线的变化速率；

3.2.通过数学形态学处理将图像中存在的空洞填满以及将毛刺滤掉，引用开启操作消除与结构元素相比尺寸较小的亮细节而保持图像整体灰度值和大的亮的区域基本不受影响，进行阈值分割，获得二值图像；

3.3.设所需内边缘点数为N_int，图像中心点为(x₀，y₀)，以Δθ＝2π/N_int间隔对二值图像做径向搜索，那么在径向上值为1的点即为边缘点(r_i，θ_i)，i＝1～N_int，通过边缘点即可得到所需的初始内边缘。

4.根据权利要求1或2所述基于内超声图像序列的血管壁边缘自动检测方法，其特征在于：所述外边缘点数N_out为36，图像帧数n取为10。

5.根据权利要求1或3所述基于内超声图像序列的血管壁边缘自动检测方法，其特征在于：所述内边缘点数N_int为36。

Images