CN109978781B - 一种基于极值区域检测的血管内超声图像分割方法 - Google Patents

Abstract

一种基于极值区域检测的血管内超声图像分割方法属于基于区域检测的医学图像分割领域。该方法针对传统IVUS图像分割方法的建模复杂、运算量大、检测速度慢的缺点，由于提取的极值区域包含很多子区域，其中包含状态不稳定的、不能代表管腔和介质的面积较大或者较小的噪声区域。因此，获取极值区域之后，本发明继续设计了基于区域稳定性的极值区域筛选算法，最终提取出精确的管腔和介质区域，完成分割。本发明分割的内外膜与标准内外膜极为接近，内膜平均HD距离为0.28mm，外膜平均HD距离为0.41mm，内外膜性能比传统方法提取的内外膜的平均性能分别提高31.4％和48.2％，且算法运行时间为线性，不仅提取精度较高、速度较快，而且鲁棒性也较好。

Description

一种基于极值区域检测的血管内超声图像分割方法

技术领域

本发明属于基于区域检测的医学图像分割领域。

背景技术

血管内超声IVUS技术是目前最先进的心脑血管疾病诊断技术，该技术结合无创性的超声技术和有创性的导管技术，在冠心病诊断、治疗中发挥着重要作用。IVUS图像动脉壁边界分割不仅对血管壁和斑块特征的定量分析至关重要，而且在血管弹性定性分析和动脉三维重建模型等医学应用中也是必需的。当前研究背景下，IVUS图像分割方法主要分为两种：半自动方法和全自动方法。半自动方法包括两类：(1)先利用计算机自动获取近似边缘轮廓，再根据轮廓信息人工进行交互修正。(2)手动粗略勾画初始轮廓，再操作计算机自动修正。上述半自动方法对IVUS图像中的噪声比较敏感，通常不能为IVUS图像的边界检测提供理想的准确度，并且在进行手动修正时，操作者要有一定的IVUS图像处理经验。这不但要利用IVUS图像本身的灰度或者纹理信息，而且还要结合患者的生理特征以及操作者自身积累的病理知识。因此，当面对大量的IVUS图像时，半自动方法很耗时，同样带有操作者的主观性。相比之下，全自动方法在整个过程中不需要人工干预，只要提供IVUS图像输入，计算机就会自动提取内外膜边缘轮廓。

目前，IVUS图像的具体分割算法主要有四种：第一种为基于统计学的图像分割。该算法通过对IVUS图像的灰度分布进行统计学建模，把求取能量极小化的问题转化为求最大后验概率的问题，实现IVUS图像的分割。它的优点是允许更加灵活的把关于轮廓及噪声特性的模型应用于算法中，但是IVUS图像中的伪影以及斑块等较为复杂的图像特征的存在会大大降低统计建模的准确性，且该算法有运算量大、收敛速度慢等缺点。第二种是基于主动轮廓Snake模型及其相关的改进算法，如：(1)引入梯度矢量流场GVF为新的外力来定义GVF-Snake能量模型。该模型很好地解决了边界凹陷处的曲线收敛问题，但它对参数的设置非常的敏感，且梯度矢量流场的计算量大，耗费时间长。(2)B样条活动轮廓B-Snake。B-Snake是从内力角度对Snake模型进行改进。在Snake模型中，曲线的弹性和刚性用曲线的一阶导和二阶导来表示，B样条在每个节点左右两侧的一阶导和二阶导连续，且可以用较少的控制点来表示较多的数据点，因此可以用B样条取代内力，只用外力来控制整条曲线的进化。B样条虽然省去了选择弹性、刚性参数的过程，但控制点个数仍需要设定。(3)拓扑自适应活动轮廓T-Snake。其思想是把图像空间划分成网格状，在模型随内外作用力的演化过程中，将重新计算模型与网格的交点，并将其作为模型新的节点，从而找到拓扑结构更加复杂的图像边界。T-Snake模型是一个拓扑可变模型，受初始轮廓的影响较小，会随着图像的边界信息进行拓扑结构的改变，具有相对较强的自适应能力。第三种是水平集模型。该模型将二维的闭合曲线的演化问题转化为三维空间中水平集函数曲面演化的隐含方式求解问题，其思想是将曲面或曲线看成一个函数的零水平集，通过函数的水平集演化确定曲线和曲面的演化。其计算精度高、算法稳定。但是，由于水平集方法不仅要对零水平集的函数求解，还要对图像的所有水平集求解，大大影响了检测速度。第四种为人工神经网络ANN方法作为特征学习算法对IVUS中外膜和管腔内膜进行检测。该方法模型过于复杂，需要大量人工标注数据进行网络训练学习，在实际运用中往往会受到较多的限制。最近，最大稳定极值区域算法MSER作为一种新的轮廓检测方法，得到很多研究者的青睐。极值区域ER是指图片中一片相似区域的集合，此集合内所有的像素值都与外界的像素值存在着较大的差距。

由上述可知，目前的方法大都采用了一种能量最小化方法或者需要注释数据来训练分类器或ANN，而且对于IVUS中外膜和管腔内膜需要分别设计分割算法，是一种串行提取方法。然而，本发明提出了一种基于极值区域的方法，不仅避开繁琐的数据标注和网络训练，而且不使用任何变分方法或可变形模型，实现并行提取内外膜，在保证了内外膜轮廓精度的情况下，大大提高了运算效率。

发明内容

本发明的目的是要提供一种基于极值区域检测的IVUS图像分割方法。该方法针对传统IVUS图像分割方法的建模复杂、运算量大、检测速度慢的缺点，提出基于极值区域的IVUS图像分割方法。由于提取的极值区域包含很多子区域，其中包含状态不稳定的、不能代表管腔和介质的面积较大或者较小的噪声区域。因此，获取极值区域之后，本发明继续设计了基于区域稳定性的极值区域筛选算法，最终提取出精确的管腔和介质区域，完成分割。结果表明，本发明方法分割的内外膜与标准内外膜极为接近，内膜平均HD距离为0.28mm，外膜平均HD距离为0.41mm，内外膜性能比传统方法提取的内外膜的平均性能分别提高31.4％和48.2％，且算法运行时间为线性，不仅提取精度较高、速度较快，而且鲁棒性也较好。

具体的实现步骤如下：

Step1：对IVUS图像进行预处理：

Step1.1：对图像I(x,y)进行中值滤波，得到图像I_m(x,y)；

Step1.2：对中值滤波后的图像I_m(x,y)，进行ROF降噪处理，得到图像I_rof(x,y)。

Step2：对图像I_rof(x,y)提取极值区域ER：

Step2.1：利用Scharr滤波器，求出图像I_rof(x,y)的梯度

并在此基础上计算边缘点，得到二值图像M(x,y)。本发明中某一边缘点p(x,y)的评价标准为点p(x,y)的梯度幅值大于其邻域P内的所有点的梯度幅值，即：

其中，||·||为点的梯度幅值函数，P是点p(x,y)的半径r＝1的邻域集合。

而且，点p(x,y)的邻域P所有点的梯度幅值的均值不小于梯度

的幅值的α倍，即：

其中，E[·]为均值函数，α为常量阈值系数，用于抑制边缘点的数目，α越大，边缘点越多，α越小，边缘点越少，α的取值范围为[0,1]，本文取α＝0.5，这样边缘点的数目既不会太少影响结果检测的精确度，也不会太多而造成冗余的计算；

Step2.2：利用图像M(x,y)将所有边缘点分为内边缘点和外边缘点，即边缘点在预处理后的图像I_rof(x,y)中的像素值大于σ时，得到内边缘点图像M⁺(x,y)，边缘点在预处理后的图像I_rof(x,y)中的像素值大于Γ-σ时，得到外边缘点图像M^-(x,y)，计算公式为：

其中，阈值σ是所有边缘点像素值的中值，由于医用IVUS图像均为8位图像，故Γ＝256；

Step2.3：计算阈值图像T_θ。当图像I_rof(x,y)的灰度值小于阈值θ时，T_θ取1，当图像I_rof(x,y)的灰度值大于阈值θ时，T_θ取0，其计算公式为：

其中，θ在[0,255]中按顺序逐一取值；

Step2.4：阈值取θ时，计算阈值图像T_θ与二值图像M^-的交点和矢量N_θ，计算公式为：

N_θ＝Θ(T_θ(x,y)·M^-(x,y)) (6)

其中，Θ(·)是矩阵元素求和运算；

Step2.5：将N_θ与阈值图像T_θ中白色像素总数目变化的比率记为V_θ：

其中，Θ(·)是矩阵元素求和运算；

Step2.6：使用中心差异，计算V_θ的加权二阶导数即获得选取极值水平EL的全局标准Ψ_θ，计算公式为：

Step2.7：获取极值水平EL。所有θ中属于极值水平EL的集合为β邻域内Ψ_θ的局部最大值，β表示局部最大值选择过程中涉及的相邻灰度相差值。β越小，Ψ_θ选取局部最大值时对比的相邻灰度值的范围越窄，为了获得足够多数量的局部最大值以便后续筛选，本文将当前Ψ_θ与其前后一个相邻的灰度值对应的全局标准作比较，即取β＝1；

Step2.8：用EL的每一位元素对图像I_rof(x,y)进行阈值化，并使用连通域算法获得极值区域ER。

Step3：筛选获得代表管腔和介质的两个极值区域ER：

Step3.1：根据极值区域ER的面积大小，进行初步过滤，过滤掉面积大于A_max和面积小于A_min的极值区域，A_max和A_min取值如下：

A_max＝(R×C)/3 (9)

A_min＝(R×C)/50 (10)

其中，R、C分别为IVUS图像的行数和列数；

Step3.2：计算每个区域中每个像素的局部二值模式特征LBP，取平均值作为每个区域的LBP特征，记为B＝{B₁,B₂,…,B_n}，其中n为ER的数量；

Step3.3：将ER区域边界到区域质心的最长距离为半径的圆的区域减去ER区域作为该ER区域的邻域，记为ER*，利用JSD距离，计算每个ER区域与其邻域ER*的灰度差异CD，记为C＝{C₁,C₂,…,C_n}，其中其中n为ER的数量，CD计算公式为：

其中，i代表所有获得的区域顺序，从第1个遍历到第n个，C_i(·)为区域的直方图分布，j从直方图分布的第1块遍历到第b块，b为直方图总块数，DS_JSD(·||·)为JSD距离函数；

Step3.4：计算每个ER区域的边界轮廓长度，记为L＝{L₁,L₂,…,L_n}，其中n为ER的数量；

Step3.5：创建矢量K＝[B₁C₁L₁,B₂C₂L₂,…,B_nC_nL_n]，其中n为ER的数量；

Step3.6：计算筛选标准矢量

其中，i代表所有获得的区域顺序，为避开边界的计算，当i＝1或i＝n时，Ω_i＝0，n为ER的数量；

Step3.7：在Ω_i的取值中，选取第二个峰值对应的ER为管腔，最后一个峰值对应的ER为介质。

Step4：分别对筛选得到的代表管腔和介质的ER进行椭圆拟合化，描绘轮廓，完成分割。

本发明突破IVUS图像传统分割方法中先设计算法提取内膜，再提取外膜或先设计算法提取外膜，再设计算法提取内膜且将前一帧图像提取的最终轮廓作为后一帧的初始轮廓的串行提取方式，采用基于极值区域的并行提取方式，同时提取内外膜，而且时间复杂度也比较低。本方法只涉及当前帧，以便其它帧提取的结果不会影响当前帧的准确性，并添加了边缘点作为图像分割的先验信息，使得提取结果更为精确。本方法不但不需要训练，而且不使用任何变分方法或可变形模型，是一种基于极值区域检测的分割方法。

附图说明

图1是本发明方法的流程图；

图2是一帧IVUS原图像；

图3是本发明筛选标准取值示意图；

图4是本发明方法的结果图。(1)内膜(2)外膜。

具体实施方式

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

一种基于极值区域的IVUS图像分割方法。先对IVUS图像进行预处理：使用中值滤波后的图像作为ROF算法的输入，得到预处理之后的图像；再对预处理之后的图像提取ER区域；然后对ER进行筛选，筛选出最能代表真实管腔和介质的区域；最后描绘区域轮廓，并对轮廓进行拟合化。

上述基于极值区域的IVUS图像分割方法，包括下述步骤：

Step1：对IVUS图像进行预处理：

Step1.1：对图像I(x,y)进行中值滤波，得到图像I_m(x,y)；

Step1.2：对中值滤波后的图像I_m(x,y)，进行ROF降噪处理，得到图像I_rof(x,y)。

Step2：对图像I_rof(x,y)提取极值区域ER：

Step2.1：利用Scharr滤波器，求出图像I_rof(x,y)的梯度

并在此基础上计算边缘点，得到二值图像M(x,y)。本发明中某一边缘点p(x,y)的评价标准为点p(x,y)的梯度幅值大于其邻域P内的所有点的梯度幅值，即：

其中，||·||为点的梯度幅值函数，P是点p(x,y)的半径r＝1的邻域集合。

而且，点p(x,y)的邻域P所有点的梯度幅值的均值不小于梯度

的幅值的α倍，即：

其中，E[·]为均值函数，α为常量阈值系数，用于抑制边缘点的数目，α越大，边缘点越多，α越小，边缘点越少，α的取值范围为[0,1]，本文取α＝0.5，这样边缘点的数目既不会太少影响结果检测的精确度，也不会太多而造成冗余的计算；

Step2.2：利用图像M(x,y)将所有边缘点分为内边缘点和外边缘点，即边缘点在预处理后的图像I_rof(x,y)中的像素值大于σ时，得到内边缘点图像M⁺(x,y)，边缘点在预处理后的图像I_rof(x,y)中的像素值大于Γ-σ时，得到外边缘点图像M^-(x,y)，计算公式为：

其中，阈值σ是所有边缘点像素值的中值，由于医用IVUS图像均为8位图像，故Γ＝256；

Step2.3：计算阈值图像T_θ。当图像I_rof(x,y)的灰度值小于阈值θ时，T_θ取1，当图像I_rof(x,y)的灰度值大于阈值θ时，T_θ取0，其计算公式为：

其中，θ在[0,255]中按顺序逐一取值；

Step2.4：阈值取θ时，计算阈值图像T_θ与二值图像M^-的交点和矢量N_θ，计算公式为：

N_θ＝Θ(T_θ(x,y)·M^-(x,y)) (6)

其中，Θ(·)是矩阵元素求和运算；

Step2.5：将N_θ与阈值图像T_θ中白色像素总数目变化的比率记为V_θ：

其中，Θ(·)是矩阵元素求和运算；

Step2.6：使用中心差异，计算V_θ的加权二阶导数即获得选取极值水平EL的全局标准Ψ_θ，计算公式为：

Step2.7：获取极值水平EL。所有θ中属于极值水平EL的集合为β邻域内Ψ_θ的局部最大值，β表示局部最大值选择过程中涉及的相邻灰度相差值。β越小，ψ_θ选取局部最大值时对比的相邻灰度值的范围越窄，为了获得足够多数量的局部最大值以便后续筛选，本文将当前ψ_θ与其前后一个相邻的灰度值对应的全局标准作比较，即取β＝1；

Step2.8：用EL的每一位元素对图像I_rof(x,y)进行阈值化，并使用连通域算法获得极值区域ER。

Step3：筛选获得代表管腔和介质的两个极值区域ER：

Step3.1：根据极值区域ER的面积大小，进行初步过滤，过滤掉面积大于A_max和面积小于A_min的极值区域，A_max和A_min取值如下：

A_max＝(R×C)/3 (9)

A_min＝(R×C)/50 (10)

其中，R、C分别为IVUS图像的行数和列数；

Step3.2：计算每个区域中每个像素的局部二值模式特征LBP，取平均值作为每个区域的LBP特征，记为B＝{B₁,B₂,…,B_n}，其中n为ER的数量；

Step3.3：将ER区域边界到区域质心的最长距离为半径的圆的区域减去ER区域作为该ER区域的邻域，记为ER*，利用JSD距离，计算每个ER区域与其邻域ER*的灰度差异CD，记为C＝{C₁,C₂,…,C_n}，其中其中n为ER的数量，CD计算公式为：

其中，i代表所有获得的区域顺序，从第1个遍历到第n个，C_i(·)为区域的直方图分布，j从直方图分布的第1块遍历到第b块，b为直方图总块数，DS_JSD(·||·)为JSD距离函数；

Step3.4：计算每个ER区域的边界轮廓长度，记为L＝{L₁,L₂,…,L_n}，其中n为ER的数量；

Step3.5：创建矢量K＝[B₁C₁L₁,B₂C₂L₂,…,B_nC_nL_n]，其中n为ER的数量；

Step3.6：计算筛选标准矢量

其中，i代表所有获得的区域顺序，为避开边界的计算，当i＝1或i＝n时，Ω_i＝0，n为ER的数量；

Step3.7：在Ω_i的取值中，选取第二个峰值对应的ER为管腔，最后一个峰值对应的ER为介质。

Step4：分别对筛选得到的代表管腔和介质的ER进行椭圆拟合化，描绘轮廓，完成分割。

Claims (1)

1.一种基于极值区域的IVUS图像分割方法，其特征在于，包括以下步骤：

Step1：对IVUS图像进行预处理：

Step1.1：对图像I(x，y)进行中值滤波，得到图像I_m(x，y)；

Step1.2：对中值滤波后的图像I_m(x，y)，进行ROF降噪处理，得到图像I_rof(x，y)；

Step2：对图像I_rof(x，y)提取极值区域ER：

Step2.1：利用Scharr滤波器，求出图像I_rof(x，y)的梯度

并在此基础上计算边缘点，得到二值图像M(x，y)；某一边缘点p(x，y)的评价标准为点p(x，y)的梯度幅值大于其邻域P内的所有点的梯度幅值，即：

其中，||·||为点的梯度幅值函数，P是点p(x，y)的半径r＝1的邻域集合；

而且，点p(x，y)的邻域P所有点的梯度幅值的均值不小于梯度

的幅值的α倍，即：

其中，E[·]为均值函数，α为常量阈值系数，用于抑制边缘点的数目，取α＝0.5；

Step2.2：利用图像M(x，y)将所有边缘点分为内边缘点和外边缘点，即边缘点在预处理后的图像I_rof(x，y)中的像素值大于σ时，得到内边缘点图像M⁺(x，y)，边缘点在预处理后的图像I_rof(x，y)中的像素值大于Γ-σ时，得到外边缘点图像M^-(x，y)，计算公式为：

其中，阈值σ是所有边缘点像素值的中值，由于医用IVUS图像均为8位图像，故Γ＝256；

Step2.3：计算阈值图像T_θ；当图像I_rof(x，y)的灰度值小于阈值θ时，T_θ取1，当图像I_rof(x，y)的灰度值大于阈值θ时，T_θ取0，其计算公式为：

其中，θ在[0，255]中按顺序逐一取值；

Step2.4：阈值取θ时，计算阈值图像T_θ与二值图像M-的交点和矢量N_θ，计算公式为：

N_θ＝Θ(T_θ(x，y)·M^-(x，y)) (6)

其中，Θ(·)是矩阵元素求和运算；

Step2.5：将N_θ与阈值图像T_θ中白色像素总数目变化的比率记为V_θ：

其中，Θ(·)是矩阵元素求和运算；

Step2.6：使用中心差异，计算V_θ的加权二阶导数即获得选取极值水平EL的全局标准Ψ_θ，计算公式为：

Step2.7：获取极值水平EL；所有θ中属于极值水平EL的集合为β邻域内Ψ_θ的局部最大值，β＝1，

Step2.8：用EL的每一位元素对图像I_rof(x，y)进行阈值化，并使用连通域算法获得极值区域ER；

Step3：筛选获得代表管腔和介质的两个极值区域ER：

Step3.1：根据极值区域ER的面积大小，进行初步过滤，过滤掉面积大于A_max和面积小于A_min的极值区域，A_max和A_min取值如下：

A_max＝(R×C)/3 (9)

A_min＝(R×C)/50 (10)

其中，R、C分别为IVUS图像的行数和列数；

Step3.2：计算每个区域中每个像素的局部二值模式特征LBP，取平均值作为每个区域的LBP特征，记为B＝{B₁，B₂，...，B_n}，其中n为ER的数量；

Step3.3：将ER区域边界到区域质心的最长距离为半径的圆的区域减去ER区域作为该ER区域的邻域，记为ER*，利用JSD距离，计算每个ER区域与其邻域ER*的灰度差异CD，记为C＝{C₁，C₂，...，C_n}，其中其中n为ER的数量，CD计算公式为：

其中，i代表所有获得的区域顺序，从第1个遍历到第n个，C_i(·)为区域的直方图分布，j从直方图分布的第1块遍历到第b块，b为直方图总块数，DS_JSD(·||·)为JSD距离函数；

Step3.4：计算每个ER区域的边界轮廓长度，记为L＝{L₁，L₂，...，L_n}，其中n为ER的数量；

Step3.5：创建矢量K＝[B₁C₁L₁，B₂C₂L₂，...，B_nC_nL_n]，其中n为ER的数量；

Step3.6：计算筛选标准矢量

其中，i代表所有获得的区域顺序，为避开边界的计算，当i＝1或i＝n时，Ω_i＝0，n为ER的数量；

Step3.7：在Ω_i的取值中，选取第二个峰值对应的ER为管腔，最后一个峰值对应的ER为介质；

Step4：分别对筛选得到的代表管腔和介质的ER进行椭圆拟合化，描绘轮廓，完成分割。

Images