CN103886599A - 一种基于血管内超声图像的血管roi分割方法 - Google Patents

Abstract

本发明涉及一种基于血管内超声图像的血管ROI分割方法。所述方法首先分割血管的管腔区域以及管腔膜轮廓，通过定位管腔区域的中心以获得参数主动轮廓模型的初始轮廓，继而通过收敛得到血管的中外膜轮廓曲线，中外膜的提取充分利用了管腔区域信息的先验知识。最后，将中外膜轮廓曲线以内区域作为ROI，通过全局最小化活动轮廓模型实现对血管斑块的分割。本发明实现了血管ROI的管腔膜中外膜以及斑块的轮廓信息可视化，与基于统计学的IVUS图像分割方法相比，摒弃了其复杂的统计建模过程且分割结果不受IVUS图像伪影和斑块特征的影响；省掉了对IVUS图像中外膜边缘进行初始轮廓的预分割步骤，提高了分割效率。

Description

一种基于血管内超声图像的血管ROI分割方法

技术领域

本发明涉及计算机医学图像分析领域，特别涉及一种基于血管内超声IVUS（IntravascularUltrasound）图像的血管感兴趣区域ROI（region of interest）分割方法。

背景技术

目前在大多数国家，心脑血管疾病已逐步成为人类死亡的最主要因素之一。冠状动脉粥样硬化病变是引起心肌梗塞和脑梗塞的主要病因，如果在早期就能够做到对硬化病变部分形态的识别以及诊断，这将对冠状动脉疾病的诊疗有着重大的意义。IVUS便是这样一种心血管疾病的超声诊断方法，血管内超声图像可以为医生显示实时的血管壁形态斑块形态进而对心血管疾病的临床诊断提供了指导意义。

然而，对于采集来的IVUS图像，尽管其成像可以显示血管壁的结构以及斑块的形态，但在实际诊疗中，IVUS图像中医生所感兴趣区域（ROI）为血管管腔区域，中外膜区域，以及管腔与中外膜之间的斑块区域。在IVUS图像的感兴趣区域中，医生往往需要目测或者根据经验确定管腔膜中外膜以及斑块轮廓的边缘，这样便不可避免地导致诊断结果不能非常客观的反映实际情况，同时也会给医生带来一定的工作困难。

因此，通过计算机对血管内超声图像进行快速准确自动地分割就显得很有必要，目前，血管内超声图像的计算机分割算法主要有三种：第一种为基于统计学的图像分割方法，文献1（G.Mendizabal-Ruiz,M.Rivera,et al.,“A probabilistic segmentation method for the identificationof luminal borders in intrvascular ultrasound images”,IEEE Conference on Computer Vision andpattern Recognition,pp.1-8,2008.）通过对IVUS图像的灰度分布进行统计学建模实现血管内超声图像的分割，但是由于血管内超声图像中的伪影以及斑块等较为复杂的图像特征的存在会大大降低统计建模的准确性。第二种方法主要是通过机器学习的手段实现血管内超声图像分割，但是这种方法模型过于复杂，在实际的运用中往往会受到较多的限制，第三种是基于活动轮廓模型的算法，文献2(宫延新孙丰荣等，“基于血液斑点噪声抑制和T-Snake模型的血管内超声图像边缘提取”，中国图像图形学报，Vol.12,No.4,pp.655-660,2007)提出的方法虽然能够得到良好的分割结果，但依赖于初始轮廓线的选取，尤其对于边缘比较模糊、纹理自相似性较高的IVUS图像，其初始轮廓更不易确定，从而也在一定程度上影响了这种分割方法的准确性。

发明内容

针对现有技术中存在的上述问题，本发明提供一种基于血管内超声图像的血管ROI分割方法，采用改进的水平集模型算法结合窄带法首先分割血管的管腔区域以及管腔膜轮廓，通过定位管腔区域的中心以获得参数主动轮廓模型（Snake）的初始轮廓，继而通过收敛得到血管的中外膜轮廓曲线，中外膜的提取充分利用了管腔区域信息的先验知识。最后，将中外膜轮廓曲线以内区域作为ROI，通过全局最小化活动轮廓模型实现对血管斑块的分割。

基于血管内超声图像的血管ROI分割方法，是以血管内超声图像为基础，根据IVUS图像的特点，对血管内超声图像依次进行血管管腔膜分割，血管中外膜分割，血管斑块分割。该方法能得到管腔区域轮廓信息，中外膜边缘信息以及斑块形状轮廓信息等，在不增加附加设备的前提下，充分利用血管内超声设备提供的超声图像本身的灰度信息以及区域特性并结合各种分割算法，实现血管ROI的分割。

本发明的特征如下：

步骤1，利用血管内超声仪，匀速回拉导管，获得人体冠状动脉的血管内超声视频影像。

步骤2，将步骤1得到的血管内超声视频影像导入计算机，从视频中截取连续的血管内超声图像序列作为实验图像，图像分辨率为2x₀*2y₀，以下简称为超声图像；

步骤3，对上述超声图像进行超声图像管腔膜的分割。

步骤3.1，创建与待分割IVUS图像大小相同的初始管腔轮廓模板矩阵。

设待分割图像为I(x,y)，创建大小同样为2x₀*2y₀的空矩阵M_init，确立中心坐标为O(x₀,y₀)，令M_init(x₀-a:x₀+b,y₀-c:y₀+d)＝1，即创建了长为c+d宽为a+b的矩形初始化区域模板。

步骤3.2，初始化符号距离映射矩阵。

设M_init为步骤3.1中设置的初始轮廓模板矩阵，计算M_init中当前像素点与最近的非零像素点的欧氏距离，单位像素距离为1，设符号距离函数映射矩阵为：

φ＝bwdist(M_init)-bwdist(1-M_init)-0.5

其中，bwdist为距离变换函数，默认计算二值图中当前像素点与最近的非0像素点的距离，并返回与原二值图相同大小的结果矩阵。

步骤3.3，设置窄带的宽度。

得到符号距离函数映射矩阵φ之后，设置窄带W_narrowband的宽度为2k，窄带W_narrowband为满足-k≤φ≤k条件的像素点组成的集合。

步骤3.4，计算窄带区域内的图像数据力F_image，公式如下：

F_image＝[I_Wnarrowband(x,y)-c_a]²-[I_Wnarrowband(x,y)-c_b]²

其中c_a、c_b分别是待分割图像I(x,y)内外两个区域上的灰度均值。

步骤3.5，计算曲率速度项V_curve。

对于窄带区域的点φ(x,y)(-k≤φ≤k)，通过八邻域的方法分别构造φ_x、φ_y、φ_xx、φ_yy和φ_xy：

φ_x＝φ(x+1,y)-φ(x-1,y)

φ_y＝φ(x,y+1)-φ(x,y-1)

φ_xx＝φ(x-1,y)+φ(x+1,y)-2φ(x,y)

φ_yy＝φ(x,y-¹)+φ(x,y+1)-2φ(x,y)

φx_y＝-0.25*φ(x-1,y-1)-0.25*φ(x+1,y+1)

+0.25*φ(x+1,y-1)+0.25*φ(x-1,y+1)

曲率速度项V_curve为：

$V_{\mathrm{curve}}=\frac{{\mathrm{\φ}}_{\mathrm{xx}}{\mathrm{\φ}}_y^2-2{\mathrm{\φ}}_y{\mathrm{\φ}}_x{\mathrm{\φ}}_{\mathrm{xy}}+{\mathrm{\φ}}_{\mathrm{yy}}{\mathrm{\φ}}_x^2}{{({\mathrm{\φ}}_x^2+{\mathrm{\φ}}_y^2)}^{\frac{3}{2}}}$

步骤3.6，计算速度函数，公式如下：

$F=\frac{F_{\mathrm{image}}}{\max \left(\right|F_{\mathrm{image}}\left|\right)}+a\·\frac{{\mathrm{\φ}}_{\mathrm{xx}}{\mathrm{\φ}}_y^2-{2\mathrm{\φ}}_y{\mathrm{\φ}}_x{\mathrm{\φ}}_{\mathrm{xy}}+{\mathrm{\φ}}_{\mathrm{yy}}{\mathrm{\φ}}_x^2}{{({\mathrm{\φ}}_x^2+{\mathrm{\φ}}_y^2)}^{\frac{3}{2}}}$

其中，α为光顺权重系数。

步骤3.7，对窄带W_narrowband内符号距离函数φ进行迭代，迭代公式为：

${\mathrm{\φ}}_W^{n+1}={\mathrm{\φ}}_W^n+{\mathrm{\Δt}}_1\·F$

式中，N为迭代次数，Δt₁＝0.45/max(F)为时间步长。

步骤3.8，经过N次迭代，最终迭代结果φ即为血管管腔区域的二值矩阵，记做管腔区域轮廓图像为f_init，提取管腔膜边缘即得到管腔膜分割结果。

步骤4，实现超声图像的中外膜分割。

步骤4.1，确定管腔轮廓区域图像的中心。

计算由步骤3中得到的血管管腔区域轮廓图像f_init的p+q阶矩m_pq，公式如下：

$m_{\mathrm{pq}}=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{i}^p{j}^p{f}_{\mathrm{init}}(i,j)$

其中，p和q都是非负整数。

管腔轮廓区域图像的中心C_centre(x_c,y_c)的坐标为：

$\left\{\begin{array}{c}{x}_c=\frac{m_{10}}{m_{00}}\\ y_c=\frac{m_{01}}{m_{00}}\end{array}\right.$

其中，m₀₀表示图像f_init的0阶矩，m₁₀和m₀₁表示图像f_init的一阶矩。

步骤4.2，初始化Snake初始轮廓。

设Snake初始轮廓曲线为C₀(s)，其坐标点为(x(s),y(s))，令 $\left\{\begin{array}{c}x\left(s\right)=x_c+r\cos \mathrm{\θ}\\ y\left(s\right)=y_c+r\sin \mathrm{\θ}\end{array}\right.,s=\frac{\mathrm{\θ}}{2\mathrm{\π}},$ s∈[0,1]，θ∈[0,2π]。

步骤4.3，求Snake轮廓曲线满足的方程。

定义C(s)＝(x(s),y(s))的能量：

$E_{\mathrm{snake}}=\underset{0}{\overset{1}{\∫}}\left\{\frac{1}{2}\right[\mathrm{\α}\left|C^{\′}\left(s\right)\right|{]}^2+\mathrm{\β}\left|C^{\′\′}\left(s\right){|}^2\right]+E_{\mathrm{ext}}\left(C\left(s\right)\right)\}\mathrm{ds}$

其中，E_ext＝▽[G_σ*I(x,y)]²为图像能量，I(x,y)为待分割图像。

对s求导得：

αC''(s)-βC''''(s)-▽E_ext＝0

设空间步长为h，轮廓线上点的坐标为X_i＝(x_i,y_i)，将导数用差分近似：

$X_i^{\′\′}=\frac{X_{i-1}-{2X}_i+X_{i+1}}{h^2},X_i^{\′\′\′\′}=\frac{[X_i-{2X}_{i-1}+X_i]-2[X_{i-1}-{2X}_i+X_{i+1}]+[X_i-{2X}_{i+1}+X_{i+2}]}{h^4}$

将上式化为X_i的矢量形式的线性方程组：

$\left\{\begin{array}{c}\mathrm{Ax}+f_x(x,y)=0\\ \mathrm{Ay}+f_y(x,y)=0\end{array}\right.$

其中，A为由α,β,h构成的对称五对角循环矩阵，(x,y)为轮廓线上点的横坐标和纵坐标组成的矢量，f_x(x,y),f_y(x,y)为由E_ext在X_i处对x,y的偏导数组成的矢量。

步骤4.4，通过迭代求解血管的中外膜轮廓。

将轮廓C视作s和时间t的函数C(s,t)，C_t(s,t)＝αC''(s,t)-βC''''(s,t)-▽E_ext，初值为：C(s,0)＝C₀(s)，迭代求解结果为：

$\left\{\begin{array}{c}{x}_i={(A+\mathrm{\γI})}^{-1}({\mathrm{\γx}}_{i-1}-f_x(x_{i-1},y_{i-1}))\\ y_i={(A+\mathrm{\γI})}^{-1}({\mathrm{\γy}}_{i-1}-f_y(x_{i-1},y_{i-1}))\end{array}\right.$

选取初始的Snake曲线作为方程的初解，当解收敛时，即C_i(s,t)≈0时得到最优曲线，该曲线即为血管的中外膜轮廓。

步骤5，利用全局最小化主动轮廓模型算法对血管斑块部分进行分割。

将步骤4得到的超声图像中外膜边缘作为超声图像ROI的临界线，选取步骤4中所分割的中外膜轮廓线以内的区域作为斑块分割的ROI，利用全局最小化主动轮廓模型算法对血管斑块部分进行分割。

步骤5.1，构造需要最小化的泛函。

基于全局最小化活动轮廓模型思想，为了得到血管内超声图像ROI的斑块轮廓，构造需要最小化的泛函：

$\begin{array}{c}{E}_1(u={\mathrm{\ψ}}_{\mathrm{\Ωc}},c_1,c_2,\mathrm{\λ})={\∫}_c\mathrm{gds}+\mathrm{\λ}{\∫}_{\mathrm{\Ω}}{((c_1-f\left(x\right))}^2-{(c_2-f\left(x\right))}^2{\mathrm{\ψ}}_{\mathrm{\Ωc}}\mathrm{dx}\\ ={\mathrm{TV}}_g\left({\mathrm{\ψ}}_{\mathrm{\Ωc}}\right)+\mathrm{\λ}{\∫}_{\mathrm{\Ω}}{r}_1(x,c_1,c_2){\mathrm{\ψ}}_{\mathrm{\Ωc}}\mathrm{dx}\end{array}$

其中，TV_g(u)是带有权重函数g(x)的函数u的总变分函数，g(x)是边缘检测函数，函数u是特征函数ψ_Ωc。

步骤5.2，构造区域信息项，表达式为：

r₁(x,c₁,c₂)＝(f_ROI(x)-c₁)²-(f_ROI(x)-c₂)²

其中，c₁为图像I(x)中u≥0.5区域的灰度均值，c₂为图像I(x)中u＜0.5区域的灰度均值，f_ROI(x)为待分割的血管内超声ROI图像。

步骤5.3，正则化泛函E₁(u＝ψ_Ωc,c₁,c₂,λ)得：

$E_2^r(u,c_1,c_2,\mathrm{\λ},a,\mathrm{\θ})={\mathrm{TV}}_g\left(u\right)+\frac{1}{2\mathrm{\θ}}\left|\right|u-v|{|}_{L^2}^2+{\∫}_{\mathrm{\Ω}}{\mathrm{\λr}}_1(x,c_1,c_2)v+\mathrm{aw}\left(v\right)\mathrm{dx}$

步骤5.4，固定v，搜索u作为

的解：

u＝v-θdivp

其中，p＝(p¹,p²)，由不动点方法解方程g(x)▽(θdivp-v)-(θdivp-v)p＝0得到：

p0＝0，且 $p^{n+1}=\frac{p^n+{\mathrm{\δ}}_t\▿({\mathrm{div}p}^n-v/\mathrm{\θ})}{1+\frac{{\mathrm{\δ}}_t}{g\left(x\right)}|\▿({\mathrm{div}p}^n-v/\mathrm{\θ})|}$

其中， $g=\frac{1}{1+\mathrm{\β}|\▿{(I*G_{\mathrm{\σ}})}^2|}.$

步骤5.5，固定u，搜索v作为 $\min \left\{\frac{1}{2\mathrm{\θ}}\right||u-v{\left|\right|}_{L^2}^2+{\∫}_{\mathrm{\Ω}}{\mathrm{\λr}}_1(x,c_1,c_2)v+\mathrm{aw}\left(v\right)\mathrm{dx}\}$ 的解：

v＝min{max{u(x)-θλr₁(x,c₁,c₂),0},1}

步骤5.6，通过迭代完成斑块分割。

初始化u₀＝0,_p0＝0，设定最大迭代次数i_max＝N，根据步骤5.2、5.4和5.5分别计算p_i，u_i，

v_i，当迭代次数达到设定的次数N时停止迭代，保存当前变量v_N即斑块分割结果。

与现有技术相比，本发明具有以下有益效果：

本发明提出了一种基于血管内超声图像的ROI分割方法，实现了血管ROI的管腔膜中外膜以及斑块的轮廓信息可视化，与基于统计学的IVUS图像分割方法相比，本发明摒弃了其复杂的统计建模过程且分割结果不受IVUS图像伪影和斑块特征的影响；与传统的基于活动轮廓模型的IVUS图像分割方法相比，本发明只需通过设置简单的矩形初始轮廓模板即可依次得到管腔区域轮廓从而定位Snake初始轮廓以进行中外膜的分割，省掉了对IVUS图像中外膜边缘进行初始轮廓的预分割这一步骤，提高了分割效率。本发明引入了血管斑块的分割，这是目前IVUS图像分割方法未涉及到的。

附图说明

图1是IVUS图像血管ROI分割组图：（a）为待分割的IVUS图像，（b）为血管管腔膜分割结果，（c）为二值化的管腔轮廓区域，（d）为血管中外膜分割结果，（e）为血管ROI图像，（f）为血管ROI的斑块轮廓显示；

图2是连续9帧血管内超声图像管腔膜分割效果图；

图3是连续9帧血管内超声图像中外膜分割效果图；

图4是连续9帧血管内超声图像ROI斑块轮廓图；

图5是本发明方法的流程图。

具体实施方式

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

一种基于血管内超声图像的血管ROI分割方法。首先，对待分割图像结合改进的水平集模型算法和窄带法实现对血管内超声图像管腔膜的分割，然后，定位管腔区域的中心并生成Snake模型的初始轮廓曲线进行迭代收敛从而得到血管中外膜的边缘，最后，选取中外膜以内区域为斑块的ROI，结合全局最小化主动轮廓模型的算法实现对血管ROI的斑块轮廓的分割。

上述基于血管内超声图像的血管ROI分割方法，包括下述步骤：

步骤1，利用血管内超声仪，以0.5mm/s的速度匀速回拉导管，获得人体冠状动脉的血管内超声视频影像；

步骤2，将步骤(1)得到的血管内超声视频影像导入计算机，从视频中截取连续的血管内超声图像序列作为实验图像，图像分辨率为384*384，以下简称为超声图像；

步骤3，对超声图像进行血管管腔膜的分割；

设待分割图像为I(x,y)，创建大小同样为384*384的空矩阵M_init，确立中心坐标为O(192,192)，令M_init(192-a:192+b,192-c:192+d)＝1，即创建了长为c+d宽为a+b的矩形初始化区域模板。

设M_init为初始轮廓模板矩阵，计算M_init中当前像素点与最近的非零像素点的欧氏距离（单位像素距离为1），设符号距离函数映射矩阵：

φ＝bwdist(M_init)-bwdist(1-M_init)-0.5（1）

其中，bwdist为距离变换函数，默认计算二值图中当前像素点与最近的非0像素点的距离，并返回与原二值图相同大小的结果矩阵。

设置窄带的宽度为2k，设窄带W_narrowband为满足-k≤φ≤k条件的像素点组成的集合。设点集P_inner为满足条件(φ≤0)的像素点的集合，点集P_outer为满足条件(φ＞0)的像素点的集合，并计算待分割图像I(x,y)在P_inner以及P_outer区域的灰度均值。

构造窄带区域内的图像数据力F_image：

F_image＝[I_Wnarrowband(x,y)-c_a]²-[I_Wnarrowband(x,y)-c_b]²（2）

对于窄带区域的点φ(x,y)(-k≤φ≤k)，采用八邻域的方法分别构造φ_x,φ_y,φ_xx,φ_yy,φ_xy

φ_x＝φ(x+1,y)-φ(x-1,y)（3）

φ_y＝φ(x,y+1)-φ(x,y-1)（4）

φ_xx＝φ(x-1,y)+φ(x+1,y)-2φ(x,y)（5）

φ_yy＝φ(x,y-1)+φ(x,y+1)-2φ(x,y)（6）

φx_y＝-0.25*φ(x-1,y-1)-0.25*φ(x+1,y+1)+0.25*φ(x+1,y-1)+0.25*φ(x-1,y+1)

（7）

通过下式计算得到曲率速度项：

$V_{\mathrm{curve}}=\frac{{\mathrm{\φ}}_{\mathrm{xx}}{\mathrm{\φ}}_y^2-2{\mathrm{\φ}}_y{\mathrm{\φ}}_x{\mathrm{\φ}}_{\mathrm{xy}}+{\mathrm{\φ}}_{\mathrm{yy}}{\mathrm{\φ}}_x^2}{{({\mathrm{\φ}}_x^2+{\mathrm{\φ}}_y^2)}^{\frac{3}{2}}}---\left(8\right)$

计算改进的水平集速度函数F：

$F=\frac{F_{\mathrm{image}}}{\max \left(\right|F_{\mathrm{image}}\left|\right)}+a\·\frac{{\mathrm{\φ}}_{\mathrm{xx}}{\mathrm{\φ}}_y^2-{2\mathrm{\φ}}_y{\mathrm{\φ}}_x{\mathrm{\φ}}_{\mathrm{xy}}+{\mathrm{\φ}}_{\mathrm{yy}}{\mathrm{\φ}}_x^2}{{({\mathrm{\φ}}_x^2+{\mathrm{\φ}}_y^2)}^{\frac{3}{2}}}---\left(9\right)$

其中光顺权重系数α为常数，通过窄带W_narrowband内的符号距离函数迭代公式

进行符号距离函数的迭代，其中N为迭代次数，时间步长设置为Δt₁＝0.45/max(F)。经过N次迭代，最终迭代结果φ即为血管管腔区域的二值矩阵，提取管腔膜边缘，记该二值矩阵为f_init，如图1（c）所示。

步骤4，实现超声图像中外膜的分割；

对于步骤3中所获得的血管管腔区域轮廓图像f_init，假设：

$m_{\mathrm{pq}}=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{i}^p{j}^p{f}_{\mathrm{init}}(i,j)---\left(10\right)$

m_pq为图像f_init的p+q阶矩，其中p和q都是非负整数，血管管腔的中心C_centre(x_c,y_c)的坐标按下式计算：

$\left\{\begin{array}{c}{x}_c=\frac{m_{10}}{m_{00}}\\ y_c=\frac{m_{01}}{m_{00}}\end{array}\right.---\left(11\right)$

通过下式得到Snake初始轮廓曲线C₀(s)上的点的坐标(x_snake,_ysnake)：

$\left\{\begin{array}{c}{x}_{\mathrm{snake}}=x_c+r\cos \mathrm{\θ}\\ y_{\mathrm{snake}}=y_c+r\sin \mathrm{\θ}\end{array}\right.---\left(12\right)$

其中，θ∈[0,2π]。

令其中s∈[0,1]。定义曲线C(s)＝(x(s),y(s))的能量为：

$E_{\mathrm{snake}}=\underset{0}{\overset{1}{\∫}}\left\{\frac{1}{2}\right[\mathrm{\α}\left|C^{\′}\left(s\right)\right|{]}^2+\mathrm{\β}\left|C^{\′\′}\left(s\right){|}^2\right]+E_{\mathrm{ext}}\left(C\left(s\right)\right)\}\mathrm{ds}---\left(13\right)$

图像能量E_ext为：

E_ext＝▽[G_σ*I(x,y)]²（14）

其中I(x,y)为待分割图像。对上式求导得到：

αC''(s)-βC''''(s)-▽E_ext＝0（15）

设空间步长为h，轮廓线上点的坐标为X_i＝(x_i,y_i)，将导数用差分近似得：

$X_i^{\′\′}=\frac{X_{i-1}-{2X}_i+X_{i+1}}{h^2}---\left(16\right)$

$X_i^{\′\′\′\′}=\frac{[X_i-{2X}_{i-1}+X_i]-2[X_{i-1}-{2X}_i+X_{i+1}]+[X_i-{2X}_{i+1}+X_{i+2}]}{h^4}---\left(17\right)$

上式可以化为如下X_i的矢量形式的线性方程组：

$\left\{\begin{array}{c}\mathrm{Ax}+f_x(x,y)=0\\ \mathrm{Ay}+f_y(x,y)=0\end{array}\right.---\left(18\right)$

其中A为由α,β,h构成的对称五对角循环矩阵，(x,y)为轮廓线上点的横坐标和纵坐标组成的矢量，f_x(x,y),f_y(x,y)为由E_ext在X_i处对x,y的偏导数组成的矢量。

将轮廓C视作s和时间t的函数C(s,t)，则：

C_t(s,t)＝αC''(s,t)-βC''''(s,t)-▽E_ext（19）

迭代求解结果为：

$\left\{\begin{array}{c}{x}_i={(A+\mathrm{\γI})}^{-1}({\mathrm{\γx}}_{i-1}-f_x(x_{i-1},y_{i-1}))\\ y_i={(A+\mathrm{\γI})}^{-1}({\mathrm{\γy}}_{i-1}-f_y(x_{i-1},y_{i-1}))\end{array}\right.---\left(20\right)$

取初值C(s,0)＝C₀(s)作为方程的初解，经过N次迭代后当解收敛时即C_i(s,t)≈0时得到最优曲线，该曲线即为血管的中外膜轮廓，轮廓线上点的坐标为(x_fin,y_fin)。

步骤5，载入超声图像的ROI图像，如图1（e）所示利用全局最小化主动轮廓模型算法，分割血管ROI斑块轮廓。

由步骤4所获得的中外膜轮廓曲线，选取中外膜轮廓曲线内部区域为血管ROI，得到超声图像的血管ROI图像f_ROI。为得到血管内超声图像ROI的斑块轮廓，基于全局最小化活动轮廓模型思想构造需要最小化的泛函

$\begin{array}{c}{E}_1(u={\mathrm{\ψ}}_{\mathrm{\Ωc}},c_1,c_2,\mathrm{\λ})={\∫}_c\mathrm{gds}+\mathrm{\λ}{\∫}_{\mathrm{\Ω}}{((c_1-f\left(x\right))}^2-{(c_2-f\left(x\right))}^2{\mathrm{\ψ}}_{\mathrm{\Ωc}}\mathrm{dx}\\ ={\mathrm{TV}}_g\left({\mathrm{\ψ}}_{\mathrm{\Ωc}}\right)+\mathrm{\λ}{\∫}_{\mathrm{\Ω}}{r}_1(x,c_1,c_2){\mathrm{\ψ}}_{\mathrm{\Ωc}}\mathrm{dx}\end{array}---\left(21\right)$

其中TV_g(u)是带有权重函数g(x)的函数u的总变分函数g(x)是边缘检测函数，函数u是特征函数ψ_Ωc。

构造区域信息项：

r₁(x,c₁,c₂)＝(f_ROI(x)-c₁)²-(f_ROI(x)-c₂)²（22）

其中c₁为图像I(x)中u≥0.5区域灰度均值，c₂为图像I(x)中u＜0.5区域灰度均值。f_ROI(x)为待分割的血管ROI图像。

正则化E₁(u＝ψ_Ωc,c₁,c₂,λ)得下式：

$E_2^r(u,c_1,c_2,\mathrm{\λ},a,\mathrm{\θ})={\mathrm{TV}}_g\left(u\right)+\frac{1}{2\mathrm{\θ}}\left|\right|u-v|{|}_{L^2}^2+{\∫}_{\mathrm{\Ω}}{\mathrm{\λr}}_1(x,c_1,c_2)v+\mathrm{aw}\left(v\right)\mathrm{dx}---\left(23\right)$

通过固定v搜索u作为

的解，解得：

u＝v-θdivp（24）

其中p＝(p¹,p²)由下式给出，

g(x)▽(θdivp-v)-(θdivp-v)p＝0（25）

上式可由不动点方法解出：

$p^{n+1}=\frac{p^n+{\mathrm{\δ}}_t\▿({\mathrm{div}p}^n-v/\mathrm{\θ})}{1+\frac{{\mathrm{\δ}}_t}{g\left(x\right)}|\▿({\mathrm{div}p}^n-v/\mathrm{\θ})|}---\left(26\right)$

其中p⁰＝0， $g=\frac{1}{1+\mathrm{\β}\left|{\▿(I*G_{\mathrm{\σ}})}^2\right|}.$

通过固定u搜索v作为 $\min \left\{\frac{1}{2\mathrm{\θ}}\right||u-v{\left|\right|}_{L^2}^2+{\∫}_{\mathrm{\Ω}}{\mathrm{\λr}}_1(x,c_1,c_2)v+\mathrm{aw}\left(v\right)\mathrm{dx}\}$ 的解，解得：

v＝min{max{u(x)-θλr₁(x,c₁,c₂),0},1}（27）

在斑块分割的迭代过程中，初始化u₀＝0,p₀＝0，设置迭代次数i，根据式（26）（24）和（22）（27）分别计算p_i，u_i，r₁ ⁱ(x,c₁,c₂)，v_i，当迭代次数达到设定的次数N时，停止迭代，保存当前变量v_N即斑块分割结果，如图1（f）所示。

Claims (2)

1.一种基于血管内超声图像的血管ROI分割方法，其特征在于，包括以下步骤：

步骤1，利用血管内超声仪，匀速回拉导管，获得人体冠状动脉的血管内超声视频影像；

步骤2，将步骤1得到的血管内超声视频影像导入计算机，从视频中截取连续的血管内超声图像序列作为实验图像，图像分辨率为2x₀*2y₀，以下简称为超声图像；

步骤3，对上述超声图像进行超声图像管腔膜的分割；

步骤3.1，创建与待分割血管内超声图像大小相同的初始管腔轮廓模板矩阵；

设待分割图像为I(x,y)，创建大小同样为2x₀*2y₀的空矩阵M_init，确立中心坐标为O(x₀,y₀)，令M_init(x₀-a:x₀+b,y₀-c:y₀+d)＝1，即创建了长为c+d宽为a+b的矩形初始化区域模板；

步骤3.2，初始化符号距离映射矩阵；

设M_init为步骤3.1中设置的初始轮廓模板矩阵，计算M_init中当前像素点与最近的非零像素点的欧氏距离，单位像素距离为1，设符号距离函数映射矩阵为：

φ＝bwdist(M_init)-bwdist(1-M_init)-0.5

其中，bwdist为距离变换函数，默认计算二值图中当前像素点与最近的非0像素点的距离，并返回与原二值图相同大小的结果矩阵；

步骤3.3，设置窄带的宽度；

得到符号距离函数映射矩阵φ之后，设置窄带W_narrowband的宽度为2k，窄带W_narrowband为满足-k≤φ≤k条件的像素点组成的集合；

步骤3.4，计算窄带区域内的图像数据力F_image，公式如下：

F_image＝[I_Wnarrowband(x,y)-c_a]²-[I_Wnarrowband(x,y)-c_b]²

其中c_a、c_b分别是待分割图像I(x,y)内外两个区域上的灰度均值；

步骤3.5，计算曲率速度项V_curve；

对于窄带区域的点φ(x,y)，-k≤φ≤k，通过八邻域的方法分别构造φ_x、φ_y、φ_xx、φ_yy和φ_xy：

φ_x＝φ(x+1,y)-φ(x-1,y)

φ_y＝φ(x,y+1)-φ(x,y-1)

φ_xx＝φ(x-1,y)+φ(x+1,y)-2φ(x,y)

φ_yy＝φ(x,y-1)+φ(x,y+1)-2φ(x,y)

φx_y＝-0.25*φ(x-1,y-1)-0.25*φ(x+1,y+1)

+0.25*φ(x+1,y-1)+0.25*φ(x-1,y+1)

曲率速度项V_curve为：

$V_{\mathrm{curve}}=\frac{{\mathrm{\φ}}_{\mathrm{xx}}{\mathrm{\φ}}_y^2-2{\mathrm{\φ}}_y{\mathrm{\φ}}_x{\mathrm{\φ}}_{\mathrm{xy}}+{\mathrm{\φ}}_{\mathrm{yy}}{\mathrm{\φ}}_x^2}{{({\mathrm{\φ}}_x^2+{\mathrm{\φ}}_y^2)}^{\frac{3}{2}}}$

步骤3.6，计算速度函数，公式如下：

$F=\frac{F_{\mathrm{image}}}{\max \left(\right|F_{\mathrm{image}}\left|\right)}+a\·\frac{{\mathrm{\φ}}_{\mathrm{xx}}{\mathrm{\φ}}_y^2-{2\mathrm{\φ}}_y{\mathrm{\φ}}_x{\mathrm{\φ}}_{\mathrm{xy}}+{\mathrm{\φ}}_{\mathrm{yy}}{\mathrm{\φ}}_x^2}{{({\mathrm{\φ}}_x^2+{\mathrm{\φ}}_y^2)}^{\frac{3}{2}}}$

其中，α为光顺权重系数；

步骤3.7，对窄带W_narrowband内符号距离函数φ进行迭代，迭代公式为：

${\mathrm{\φ}}_W^{n+1}={\mathrm{\φ}}_W^n+{\mathrm{\Δt}}_1\·F$

式中，N为迭代次数，Δt₁＝0.45/max(F)为时间步长；

步骤3.8，经过N次迭代，最终迭代结果φ即为血管管腔区域的二值矩阵，记做管腔区域轮廓图像为f_init，提取管腔膜边缘即得到管腔膜分割结果；

步骤4，实现超声图像的中外膜分割；

步骤4.1，确定管腔轮廓区域图像的中心；

计算由步骤3中得到的血管管腔区域轮廓图像f_init的p+q阶矩m_pq，公式如下：

$m_{\mathrm{pq}}=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{i}^p{j}^p{f}_{\mathrm{init}}(i,j)$

其中，p和q都是非负整数；

管腔轮廓区域图像的中心C_centre(x_c,y_c)的坐标为：

$\left\{\begin{array}{c}{x}_c=\frac{m_{10}}{m_{00}}\\ y_c=\frac{m_{01}}{m_{00}}\end{array}\right.$

其中，m₀₀表示图像f_init的0阶矩，m₁₀和m₀₁表示图像f_init的一阶矩；

步骤4.2，初始化Snake初始轮廓；

设Snake初始轮廓曲线为C0(s)，其坐标点为(x(s),y(s))，令 $\left\{\begin{array}{c}x\left(s\right)=x_c+r\cos \mathrm{\θ}\\ y\left(s\right)=y_c+r\sin \mathrm{\θ}\end{array}\right.,s=\frac{\mathrm{\θ}}{2\mathrm{\π}},$ s∈[0,1]，θ∈[0,2π]；

步骤4.3，求Snake轮廓曲线满足的方程；

定义C(s)＝(x(s),y(s))的能量：

$E_{\mathrm{snake}}=\underset{0}{\overset{1}{\∫}}\left\{\frac{1}{2}\right[\mathrm{\α}\left|C^{\′}\left(s\right)\right|{]}^2+\mathrm{\β}\left|C^{\′\′}\left(s\right){|}^2\right]+E_{\mathrm{ext}}\left(C\left(s\right)\right)\}\mathrm{ds}$

其中，E_ext＝▽[G_σ*I(x,y)]²为图像能量，I(x,y)为待分割图像；

对s求导得：

αC''(s)-βC''''(s)-▽E_ext＝0

设空间步长为h，轮廓线上点的坐标为X_i＝(x_i,y_i)，将导数用差分近似，

$X_i^{\′\′}=\frac{X_{i-1}-{2X}_i+X_{i+1}}{h^2},X_i^{\′\′\′\′}=\frac{[X_i-{2X}_{i-1}+X_i]-2[X_{i-1}-{2X}_i+X_{i+1}]+[X_i-{2X}_{i+1}+X_{i+2}]}{h^4}$

将上式化为X_i的矢量形式的线性方程组：

$\begin{array}{c}\mathrm{Ax}+f_x(x,y)=0\\ \mathrm{Ay}+f_y(x,y)=0\end{array}$

其中，A为由α,β,h构成的对称五对角循环矩阵，(x,y)为轮廓线上点的横坐标和纵坐标组成的矢量，f_x(x,y),f_y(x,y)为由E_ext在X_i处对x,y的偏导数组成的矢量；

步骤4.4，通过迭代求解血管的中外膜轮廓；

将轮廓C视作s和时间t的函数C(s,t)，C_t(s,t)＝αC''(s,t)-βC''''(s,t)-▽E_ext，初值为：

C(s,0)＝C₀(s)，迭代求解结果为：

$\left\{\begin{array}{c}{x}_i={(A+\mathrm{\γI})}^{-1}({\mathrm{\γx}}_{i-1}-f_x(x_{i-1},y_{i-1}))\\ y_i={(A+\mathrm{\γI})}^{-1}({\mathrm{\γy}}_{i-1}-f_y(x_{i-1},y_{i-1}))\end{array}\right.$

选取初始的Snake曲线作为方程的初解，当解收敛时，即C_i(s,t)≈0时得到最优曲线，该曲线即为血管的中外膜轮廓；

步骤5，利用全局最小化主动轮廓模型算法对血管斑块部分进行分割；

将步骤4得到的超声图像中外膜边缘作为超声图像ROI的临界线，选取步骤4中所分割的中外膜轮廓线以内的区域作为斑块分割的ROI，利用全局最小化主动轮廓模型算法对血管斑块部分进行分割。

2.根据权利要求1所述的一种基于血管内超声图像的血管ROI分割方法，其特征在于，所述步骤5利用全局最小化主动轮廓模型算法对血管斑块部分进行分割的方法包括以下步骤：

步骤5.1，构造需要最小化的泛函；

基于全局最小化活动轮廓模型思想，为了得到血管内超声图像ROI的斑块轮廓，构造需要最小化的泛函：

$\begin{array}{c}{E}_1(u={\mathrm{\ψ}}_{\mathrm{\Ωc}},c_1,c_2,\mathrm{\λ})={\∫}_c\mathrm{gds}+\mathrm{\λ}{\∫}_{\mathrm{\Ω}}{((c_1-f\left(x\right))}^2-{(c_2-f\left(x\right))}^2{\mathrm{\ψ}}_{\mathrm{\Ωc}}\mathrm{dx}\\ ={\mathrm{TV}}_g\left({\mathrm{\ψ}}_{\mathrm{\Ωc}}\right)+\mathrm{\λ}{\∫}_{\mathrm{\Ω}}{r}_1(x,c_1,c_2){\mathrm{\ψ}}_{\mathrm{\Ωc}}\mathrm{dx}\end{array}$

其中，TV_g(u)是带有权重函数g(x)的函数u的总变分函数，g(x)是边缘检测函数，函数u是特征函数ψ_Ωc；

步骤5.2，构造区域信息项，表达式为：

r₁(x,c₁,c₂)＝(f_ROI(x)-c₁)²-(f_ROI(x)-c₂)²

其中，c₁为图像I(x)中u≥0.5区域的灰度均值，c₂为图像I(x)中u＜0.5区域的灰度均值，f_ROI(x)为待分割的血管内超声ROI图像；

步骤5.3，正则化泛函E₁(u＝ψ_Ωc,c₁,c₂,λ)得：

$E_2^r(u,c_1,c_2,\mathrm{\λ},a,\mathrm{\θ})={\mathrm{TV}}_g\left(u\right)+\frac{1}{2\mathrm{\θ}}\left|\right|u-v|{|}_{L^2}^2+{\∫}_{\mathrm{\Ω}}{\mathrm{\λr}}_1(x,c_1,c_2)v+\mathrm{aw}\left(v\right)\mathrm{dx}$

步骤5.4，固定v，搜索u作为

的解：

u＝v-θdivp

其中，p＝(p¹,p²)，由不动点方法解方程g(x)▽(θdivp-v)-(θdivp-v)p＝0得到：

p⁰＝0，且 $p^{n+1}=\frac{p^n+{\mathrm{\δ}}_t\▿({\mathrm{div}p}^n-v/\mathrm{\θ})}{1+\frac{{\mathrm{\δ}}_t}{g\left(x\right)}|\▿({\mathrm{div}p}^n-v/\mathrm{\θ})|}$

其中， $g=\frac{1}{1+\mathrm{\β}\left|{\▿(I*G_{\mathrm{\σ}})}^2\right|};$

步骤5.5，固定u，搜索v作为 $\min \left\{\frac{1}{2\mathrm{\θ}}\right||u-v{\left|\right|}_{L^2}^2+{\∫}_{\mathrm{\Ω}}{\mathrm{\λr}}_1(x,c_1,c_2)v+\mathrm{aw}\left(v\right)\mathrm{dx}\}$ 的解：

v＝min{max{u(x)-θλr₁(x,c₁,c₂),0},1}

步骤5.6，通过迭代完成斑块分割；

初始化u₀＝0,p₀＝0，设定最大迭代次数i_max＝N，根据步骤5.2、5.4和5.5分别计算p_i，u_i，v_i，当迭代次数达到设定的次数N时停止迭代，保存当前变量v_N即斑块分割结果。

Images