动脉粥样硬化是一种血管壁疾病，是心脏病或中风等心血管疾病的主要产生原因(Frostegård，2005)。血管内超声是一种基于导管的血管内成像技术，已经被广泛应用于评价冠状动脉粥样硬化特点及程度(Vykoukal等，2014)。它是将发射频率为20~40 MHz的微小超声源镶嵌在导管的顶端，放置于冠状动脉血管腔内(Mendizabal-Ruiz等，2013)。超声源发出的超声束通过不同层物质的反射和散射以及导管的匀速回拉而获得的血管壁及血管腔的序列切面图像(Nissen和Yock，2001)，其结构如图 1所示。IVUS图像包含3个区域：管腔区域、介质区域和血管外膜层区域。这3个区域被两个闭合轮廓分割：管腔内膜，对应于管腔壁边界，位于管腔和介质之间；外膜，位于介质和血管外膜层之间。IVUS图像分割的研究内容即如何准确地获取管腔内膜和外膜轮廓。

IVUS图像分割有半自动方法和全自动方法。半自动方法在提取内外膜时需要手工交互进行轮廓初始化。该方法对IVUS图像噪声比较敏感，提取精度不高。且手工参与的工作量较大，不同专家的手工检测结果存在主观差异。例如，Hossain等人(2013)和Sofian等人(2015)提出的算法在提取轮廓前，需在每一帧IVUS图像中手动创建一组初始轮廓点，严重影响了分割速度。

相比之下，全自动方法在整个过程中不需要人为干预，由计算机自动分析IVUS图像特征，可快速提取内外膜边缘(Mendizabal-Ruiz等，2012)。例如，Ciompi等人(2012)、Mendizabal-Ruiz等人(2008)、Jodas等人(2017)以及袁绍锋等人(2018a)的方法均为全自动方法，分割速度较快，但也各自存在不足之处：Ciompi等人(2012)将外膜视为整个血管的一部分，不仅考虑血管中的所有组织区域，而且建立外膜与血管其他组织区域之间的关系，从而实现稳健的分割。但该自动方法只能用于检测外膜。此外，Mendizabal-Ruiz等人(2008)提出的轮廓检测算法中，轮廓的平滑度受到IVUS图像中导管和导丝伪影的干扰，影响了分割精度，且该方法只能用于检测管腔内膜。Jodas等人(2017)将K均值算法和平均圆度应用于IVUS图像分割，并提出了一种识别和消除分岔伪影的算法，但也只能检测到管腔内膜。袁绍锋等人(2018a)构建了多分类随机决策森林模型实现边界检测，可以较好克服IVUS临床图像中声影和硬斑块干扰等问题，但也易受超声阴影和血管分支的干扰，这类影响甚至导致其检测失败。

根据轮廓分割原理，目前IVUS图像分割算法主要分为4种：基于统计学的分割、基于主动轮廓模型(Snake)的分割、基于水平集模型的分割以及基于机器学习的分割。但这4种方法模型的建立均过于复杂。

1) 基于统计学的分割，Mendizabal-Ruiz等人(2013)根据血液像素和非血液像素的概率定义了一个概率场，通过最小化概率成本函数使参数曲线变形。该种方法是对IVUS图像的灰度分布进行统计学建模完成分割，但伪影的存在大大降低了统计建模的准确性。

2) 基于Snake模型的分割，高重阳(2010)对Snake模型进行改进，将相似性测度信息应用到相邻图像的轮廓点上，用一致性测度构造Snake模型的图像力，让曲线演化到边缘。Zhu等人(2011)提出了一种基于梯度向量流模型(GVF-Snake)的方法，在GVF场上利用非线性滤波改变曲线的形态结构以完成分割。林慕丹等人(2016)结合轮廓的先验形状信息，剔除血管区域中不满足曲线曲率要求的部分从而提取关键点，最后结合梯度和相位信息，采用Snake模型进一步细化轮廓，提高了对巧化、纤维斑块及声影区域的识别能力，但该方法因有监督分类误判的血管区轮廓符合冠脉血管的凹凸性要求而导致算法无法校正误差，从而无法正确检测目标轮廓。基于Snake模型的方法虽然能够得到良好的分割结果，但依赖初始轮廓线的选取，而且将上一帧获取的最终轮廓作为当前帧的初始轮廓进行迭代，这样做的缺点是当某一帧图像受到复杂伪影的干扰而检测错误时，后续帧都将受到影响，误差较大。

3) 基于水平集的分割，Lazrag和Naceur(2012)提出基于水平集模型和Contourlet变换相结合的方法。该方法先利用Contourlet变换将原始图像分解为低通和带通定向带，然后利用水平集函数获得最终轮廓。方婷(2018)针对灰度不均匀图像分割效果不理想，以及对初始曲线敏感问题，设计了结合局部熵的混合水平集分割模型，进一步提高对局部灰度信息的利用率，达到分割的目的。但水平集方法要对所有水平集求解，严重降低了检测速度。

4) 基于机器学习的分割，Su等人(2017)用图像的空间和邻域特征作为数据输入，然后基于机器学习的方式，通过编码器和分类器区分不同的血管层。Yang等人(2018)提出IVUS-Net的全卷积神经网络，为管腔或介质预测二进制掩模，然后进行轮廓提取。袁绍锋等人(2018b)使用两个深度全卷积网络(DFCN)对图像和不同的标注信息进行感兴趣和非感兴趣区域分割训练，然后结合心血管先验形状信息对边界进行优化，提高了各类斑块和伪影的识别能力，能准确、可重复地检测出目标边界，但是该算法针对内外膜有不同的边界优化算法的组合，且不同的算法组合会对边界检测结果造成不同程度的影响。袁绍锋等人(2019)提出一种结合堆叠沙漏网络(SHGN)和有条件生成对抗网络(C-GAN)的IVUS内膜与中-外膜检测的改进方法，可以实现跨数据集的分割。但基于机器学习的算法，需要大量人工标注数据进行网络训练，在实际运用中往往会受到较多的限制。

最大稳定极值区域算法(MSER)作为一种新的轮廓检测方法，得到了很多研究者的青睐。极值区域(ER)(Faraji等，2015)是指图像中一片相似区域的集合，此集合内所有的像素值都与外界的像素值存在着较大的差异。例如，为了得到目标物轮廓，Tabassum和Dhondse(2015)和Yan等人(2018)利用最大稳定极值区域算法先提取候选区域，然后在候选区域中定位目标。MSER算法虽然定位准确度较高，但稳定性标准的选择对区域选择有较大影响，而且区域清理步骤复杂。

因此，IVUS图像分割仍面临诸多挑战：1)内膜和外膜的分割不能同时实现，实时性不高：传统方法中，由于IVUS图像内膜和外膜所处的“环境”不同，需要针对内膜和外膜的不同特点分别设计算法进行提取，该种串行提取方式不仅算法复杂，而且时间效率低下；2)IVUS图像中存在复杂的伪影，如分岔、侧血管、导丝伪影、阴影伪影等影响分割精度(Balocco等，2014)；3)建模复杂：需要根据不同的图像信息建立复杂的模型完成分割。针对IVUS图像分割面临的挑战，本文采用新的极值区域检测算法，结合了MSER算法的不足，不仅避开了繁琐的数据标注和网络训练，而且不使用任何复杂的可变模型，实现并行提取内外膜，在保证内外膜轮廓精度的情况下，也提高了算法鲁棒性。本文的主要贡献如下：

1) 将IVUS图像的边缘检测问题放宽到区域检测和区域选择问题，无需建立复杂的模型。

2) 本文设计的算法既适应于内膜提取，又适应于外膜提取，所以可以实现并行提取内外膜。

3) 利用局部二值模式(LBP)特征(Ojala等，1994)、灰度差异和边缘周长等信息设计区域筛选矢量，以便从候选区域中获取准确的管腔区域和介质区域，同时这些信息的使用提高了算法的鲁棒性。

算法流程如图 2所示，主要分为3大步骤：极值区域检测、极值区域筛选以及轮廓拟合。极值区域检测用于获得单帧IVUS图像的所有极值区域以供后续步骤筛选，主要包含边缘点检测、计算极值水平标准、获取极值水平和提取极值区域4部分；极值区域筛选的目的是获得与管腔和介质最接近的两个极值区域，包含面积筛选以及矢量筛选。最后对筛选得到的代表管腔和介质的两个极值区域进行轮廓拟合，完成分割。

1) 边缘点检测。某一边缘点$p(x, y)$评价标准为该点的梯度幅值大于其邻域$\mathit{\boldsymbol{P}}$内的所有像素点的梯度幅值(Faraji等，2015)，即

$ \left\| {\nabla \mathit{\boldsymbol{I}}\left({p\left({x, y} \right)} \right)} \right\| > \left\| {\nabla \mathit{\boldsymbol{I}}\left(\mathit{\boldsymbol{P}} \right)} \right\| $

(1)

式中，$\nabla \mathit{\boldsymbol{I}}$是图像$\mathit{\boldsymbol{I}}$的梯度，||·||为求点的梯度幅值函数，$\mathit{\boldsymbol{P}}$是半径$r$=1的邻域像素点集合。此外，点$p(x, y)$的邻域$\mathit{\boldsymbol{P}}$内所有点的梯度幅值的均值不小于梯度$\nabla \mathit{\boldsymbol{I}}$幅值的$\alpha $倍，即

$ E\left[ {\left\| {\nabla \mathit{\boldsymbol{I}}\left(\mathit{\boldsymbol{P}} \right)} \right\|} \right] \ge \alpha E\left[ {\left\| {\nabla \mathit{\boldsymbol{I}}} \right\|} \right] $

(2)

式中，$E\left[ \cdot \right]$为均值函数，$\alpha $为常量阈值系数，用于抑制边缘点的数目，$\alpha $越大，边缘点的数目越少，反之越多，其取值一般为[0, 1]，本文取$\alpha = 0.5$，适中的边缘点数量既不会因为太少影响结果检测的精确度，也不会因为太多而造成冗余的计算。

将所有边缘点划分为两部分，一部分为从暗区域到亮区域过渡的边缘点，称为内边缘点$\left({{M^ + }\left({x, y} \right)} \right)$，另一部分为从亮区域到暗区域过渡的边缘点，称为外边缘点$\left({{M^ - }\left({x, y} \right)} \right)$，分别定义为

$ {M^ + }\left({x, y} \right) = \left\{ {\begin{array}{*{20}{l}} 1&{M\left({x, y} \right) \times I\left({x, y} \right) \ge \sigma }\\ 0&{其他} \end{array}} \right. $

(3)

$ {M^ - }\left({x, y} \right) = \left\{ {\begin{array}{*{20}{l}} 1&{M\left({x, y} \right) \times I\left({x, y} \right) \ge \mathit{\Gamma } - \sigma }\\ 0&{其他} \end{array}} \right. $

(4)

式中，$M\left({x, y} \right)$为标记所有边缘点的二值图像中的像素，对于$m$位的图像，$\mathit{\Gamma } = {2^m}$，$\sigma $取所有边缘点像素的中值。假设${P_i}$为某一边缘点灰度值的概率，为保证边缘点划分为均匀，即$\sum\limits_{i = 0}^\sigma {{P_i} \ge 1/2} $与$\sum\limits_{i = 0}^{\sigma - 1} {{P_i} < 1/2} $相等，则$\sigma $为所有边缘点像素值的中位数。

2) 计算极值水平标准。对于$m$位图像，灰度取值范围为$\left[ {0, {2^m} - 1} \right]$，若直接用图像每一个可能的灰度值分别对图像进行阈值化，则工作量十分庞大，效率低下。因此，需要一个能够舍弃没用阈值的全局矢量。首先分别以图像$\mathit{\boldsymbol{I}}$的每一个灰度值作为阈值$\theta $，计算阈值图像，即

$ {T_\theta }\left({x, y} \right) = \left\{ {\begin{array}{*{20}{l}} 1&{I\left({x, y} \right) \le \theta }\\ 0&{其他} \end{array}} \right. $

(5)

式中，$\theta $在图像灰度范围$\left[ {0, {2^m} - 1} \right]$中逐一取值。

然后在灰度阈值$\theta $下，计算${\mathit{\boldsymbol{T}}_\theta }$与外边缘点图像${\mathit{\boldsymbol{M}}^ - }$的交点和矢量${\mathit{\boldsymbol{N}}_\theta }$，即二者相交的点的数量为

$ {\mathit{\boldsymbol{N}}_\theta } = \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}\left({{\mathit{\boldsymbol{T}}_\theta } \cdot {\mathit{\boldsymbol{M}}^ - }} \right) $

(6)

式中，$\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}$(·)是矩阵元素求和运算。

为获得稳定的阈值图像，除了考虑与阈值图像相交的边缘点的数量变化是否稳定外，还应考虑阈值图像中包含的所有像素总数量的变化是否稳定。因此，将矢量${\mathit{\boldsymbol{V}}_\theta}$定义为${\mathit{\boldsymbol{N}}_\theta}$与当前阈值图像中像素变化量的比率，即

$ {\mathit{\boldsymbol{V}}_\theta } = \frac{{{\mathit{\boldsymbol{N}}_\theta }}}{{1 + \mathit{\boldsymbol{ \boldsymbol{\varTheta} }}\left({{\mathit{\boldsymbol{T}}_\theta }} \right)}}, 0 \le \theta \le {2^m} - 1 $

(7)

所有灰度阈值$\theta $下计算得到的${\mathit{\boldsymbol{V}}_\theta}$组成矢量$\mathit{\boldsymbol{V}}$, 根据矢量$\mathit{\boldsymbol{V}}$，全局标准$\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}$的每一个元素定义为

$ {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_\theta } = \frac{{{\mathit{\boldsymbol{V}}_\theta }}}{{\left| {{\mathit{\boldsymbol{V}}_{\theta + 1}} - {\mathit{\boldsymbol{V}}_{\theta - 1}}} \right| + \varepsilon }} $

(8)

式中，为防止分母为0，$\varepsilon $是一个足够小的常量。由定义可以看出，${\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_\theta }$取值越大，${\mathit{\boldsymbol{V}}_\theta}$的变化越小，则阈值图像中总像素数量的变化越稳定，表明该阈值$\theta $下获得的阈值图像越稳定。

3) 获取极值水平。由以上分析可知，极值水平就是$\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}$取极大值时灰度阈值$\theta $的集合。为了计算矢量$\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}$的极大值，每个${\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_\theta }$需分别与其前后$\beta $个值作比较，若该${\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_\theta }$在该$\beta $邻域内为极大值，则该灰度阈值$\theta $就属于极值水平。另外邻域窗口半径越大，矢量$\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}$的极大值就越少，极值水平的数量也就越少，本文取$\beta $=1，以便获得最多数量的极值水平以供后续分析。

4) 提取极值区域。用极值水平中的每一个值分别对图像做二值化处理，之后对每一个二值图像做形态学闭运算，即先膨胀后腐蚀以填充区域内的细小空洞且不明显改变面积。然后，使用连通域标记算法(He等，2017)标记图像的连通域，此时的连通域即是极值区域。为便于后续处理，需要获取极值区域的相关参数，如区域的轮廓曲线坐标、区域内的像素值和像素坐标等。

为减少后续矢量筛选的工作量，上述获取的极值区域首先需做面积筛选，过滤掉面积过大(大于${A_{\max }}$)或者过小的(小于${A_{\min }}$)不能代表管腔和介质的极值区域以得到候选区域。根据图像尺寸以及内外膜大致面积范围，设置${A_{\max }} = \left({R{\rm{ \times }}C} \right)/3$，${A_{\min }} = \left({R{\rm{ \times }}C} \right)/50$，其中，$R$、$C$分别代表IVUS图像的高和宽。

然后，设计极值区域的筛选标准矢量，以便在候选区域中找到分别代表管腔和介质的极值区域。构建矢量$\mathit{\boldsymbol{K}}$(其中每一个元素都对应一个极值区域)，为了在矢量$\mathit{\boldsymbol{K}}$中找到对应的最稳定的管腔和介质区域，仍用求局部最大值的方法。为此，将区域的LBP特征值(Ojala等，1994)、灰度差异和极值区域轮廓周长作为矢量$\mathit{\boldsymbol{K}}$的组成部分。LBP特征是用来描述图像纹理特征的算子，真实的管腔和介质区域的平均LBP特征值应当足够稳定。灰度差异是指区域内的灰度值和其相邻区域的差异，两个区域的灰度差异用JSD距离(Majtey等，2005)度量。真实内外膜与其邻域有一定的灰度差异，且灰度差异越大，区域越稳定。极值区域轮廓周长是对区域边界的度量，目的是找到边界变化相对稳定的区域。

1) LBP(local binary pattern)特征值。原始的LBP算子(Ojala等，1994)定义在3×3的邻域窗口内，以窗口中心像素为阈值，并与相邻8个像素的灰度值进行比较，若周围像素值大于中心像素值，该对应像素点被标记为1，否则为0，则该8位二进制数转为十进制数后即为该窗口中心像素点的LBP值。

改进的LBP算子${L_{{\rm{BP}}}}\left({P, R} \right)$(Ahonen等，2004)可以计算不同邻域半径和不同像素点数的特征值，其中$P$表示周围像素点个数，$R$表示邻域半径。本文使用改进的LBP特征算法(Ahonen等，2004)计算每个区域中每个像素的LBP特征，取平均值作为每个极值区域的LBP特征，记为$\mathit{\boldsymbol{B = }}\left\{ {{B_1}, {B_2}, \cdots, {B_i}, \cdots, {B_n}} \right\}$，其中$i$表示第$i$个区域，$n$为候选极值区域数量。为保证亮度的同时，让图像纹理更精细，取$R$=1.5，$P$=12。

2) 灰度差异。本文将极值区域轮廓到区域质心的最长距离为半径的圆与极值区域的差作为该极值区域的邻域$\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{R}}^*}$，则利用JSD(Jensen-Shannon divergence)距离(Majtey等，2005)，计算每个极值区域与其邻域的灰度差异，记为$\mathit{\boldsymbol{C = }}\left\{ {{C_1}, {C_2}, \cdots, {C_i}, \cdots, {C_n}} \right\}$,

$ {C_i} = \sqrt {D{S_{JSD}}\left({{C_i}\left({\mathit{\boldsymbol{ER}}} \right)\left\| {{C_i}\left({\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{R}}^ * }} \right)} \right.} \right)} $

(9)

式中，${C_i}(\cdot)$为第$i$个区域的直方图分布，$D{S_{JSD}}(\cdot \left\| \cdot \right.)$为JSD距离函数(Majtey等，2005)，定义为

$ \begin{array}{*{20}{c}} {D{S_{JSD}}\left({{C_i}\left({\mathit{\boldsymbol{ER}}} \right)\left\| {{C_i}\left({\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{R}}^ * }} \right)} \right.} \right) = }\\ {1/2\left({D{S_{KID}}\left({{C_i}\left({\mathit{\boldsymbol{ER}}} \right)\left\| {{C_i}\left({\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{R}}_a}} \right)} \right.} \right) + } \right.}\\ {D{S_{KID}}\left({{C_i}\left({\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{R}}^ * }} \right)\left\| {{C_i}\left({\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{R}}_a}} \right)} \right.} \right)} \end{array} $

(10)

式中，${C_i}\left({\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{R}}_a}} \right) = 1/2\left({{C_i}(\mathit{\boldsymbol{ER}}) + {C_i}\left({\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{R}}^*}} \right)} \right)$，$D{S_{KLD}}(\cdot \left\| \cdot \right.)$为两个区域的KLD(Kullback-Leibler divegonce)距离(Klein和Frintrop，2011)，即

$ \begin{array}{*{20}{c}} {D{S_{KID}}\left({{C_i}\left({\mathit{\boldsymbol{ER}}} \right)\left\| {{C_i}\left({\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{R}}_a}} \right)} \right. = } \right.}\\ {{C_i}\left({\mathit{\boldsymbol{ER}}} \right)\log \frac{{{C_i}\left({\mathit{\boldsymbol{ER}}} \right)}}{{{C_i}\left({\mathit{\boldsymbol{E}}{\mathit{\boldsymbol{R}}_a}} \right)}}} \end{array} $

(11)

3) 极值区域轮廓周长。在提取极值区域过程中，由于使用了连通域算法，因此可以直接获取每个区域的边界长度、区域面积、质心坐标等参数。将每个候选区域轮廓周长表示为$\mathit{\boldsymbol{L = }}\left\{ {{L_1}, {L_2}, \cdots, {L_i}, \cdots, {L_n}} \right\}$。

4) 矢量筛选。针对每个候选区域，将以上3个特征的乘积作为矢量$\mathit{\boldsymbol{K}}$的元素，即

$ \mathit{\boldsymbol{K}} = \left[ {{B_1}\;{C_1}\;{L_2}, {B_2}\;{C_2}\;{L_2}, \cdots, {K_i}, \cdots, {B_n}\;{C_n}\;{L_n}} \right] $

(12)

从矢量$\mathit{\boldsymbol{K}}$的构成可以看出，$\mathit{\boldsymbol{K}}$的变化代表区域纹理信息、灰度差异，以及区域边界长度的变化。$\mathit{\boldsymbol{K}}$中某个元素变化越稳定，说明该元素对应的极值区域越稳定，成为管腔或介质区域的可能性也就越大。因此建立一个代表$\mathit{\boldsymbol{K}}$变化情况的矢量$\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$，表示管腔和介质极值区域的稳定程度，其定义为

$ {\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i} = \frac{{{\mathit{\boldsymbol{K}}_i}}}{{\varepsilon + \left({{\mathit{\boldsymbol{K}}_{i + 1}} - {\mathit{\boldsymbol{K}}_{i - 1}}} \right)}} $

(13)

式中，$\varepsilon $为一较小的常量。${\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}$取值越大，则${\mathit{\boldsymbol{K}}_i}$变化越少，即第$i$个极值区域就越稳定。所以从${\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}$的峰值中寻找代表管腔和介质的极值区域。图 3为选取的一帧IVUS图像的前3个，中间3个以及后3个峰值对应的极值区域轮廓。由于提取的极值区域范围由内膜附近一直扩展到外膜附近，且是相互嵌套的区域，所以矢量$\mathit{\boldsymbol{K}}$是递增的，也就是说管腔区域在前几个峰值对应的极值区域中，而介质区域在后几个峰值对应的极值区域中，中间的峰值代表的极值区域对应于内外膜之间的部分。当IVUS图像没有伪像时，${\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}$的峰值较多，而且较为明显，此时选择${\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}$的第2个峰值对应的极值区域作为管腔区域，倒数第2个峰值对应的极值区域作为介质区域。当${\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}_i}$的峰值较少时，直接将最后一个峰值对应的极值区域作为介质区域，如图 4所示为一帧IVUS图像极值区域的峰值示意图。

根据内外膜的固有形状特征，本文采用椭圆对轮廓进行拟合(Szpak等，2012)。将不规则的极值区域拟合为椭圆区域，需要求出椭圆的长轴$l$、短轴$\omega $，以及和水平方向的夹角$\theta $这3个参数。

根据待拟合区域$\xi $的重心位置$\left({{x_c}, {y_c}} \right)$, 计算该区域的中心二阶矩矩阵，$\mathit{\boldsymbol{U}} = \left[ {\begin{array}{*{20}{l}} {{\mu _{20}}}&{{\mu _{11}}}\\ {{\mu _{11}}}&{{\mu _{02}}} \end{array}} \right]$，即

$ {\mu _{20}} = \sum\limits_\xi {{{\left({x - {x_c}} \right)}^2}} I\left({x, y} \right) $

(14)

$ {\mu _{11}} = \sum\limits_\xi {\left({x - {x_c}} \right)} \left({y - {y_c}} \right)I\left({x, y} \right) $

(15)

$ {\mu _{02}} = \sum\limits_\xi {{{\left({y - {y_c}} \right)}^2}} I\left({x, y} \right) $

(16)

因此，椭圆的3个参数可以计算为

$ \theta = \frac{1}{2}\arctan \left({\frac{{2{\mu _{11}}}}{{{\mu _{20}} + {\mu _{02}}}}} \right) $

(17)

$ l = \sqrt {\frac{1}{{{\lambda _{\min }}}}} $

(18)

$ \omega = \sqrt {\frac{1}{{{\lambda _{\max }}}}} $

(19)

式中，${\lambda _{\min }}$和${\lambda _{\max }}$分别是二阶矩矩阵的最小、最大特征值。

鉴于现有方法多数采用标准公开数据集中20 MHz的图像展示分割效果，所以为了公平地与已有算法进行比较，实验使用的是包含326帧20 MHz的IVUS B模式图像的标准公开数据集，每帧图像的分辨率是384×384像素(Balocco等，2013)。数据集中还提供了两位专家手工标注的内膜和外膜的轮廓坐标。其中一位专家在标注完成一周之后进行了重复标注，以提高数据标注的准确性。另外，使用从北京安贞医院获取的临床数据集以验证本文算法在临床数据中的分割能力，该数据集包含432幅图像，图像分辨率仍为384×384像素。

实验首先选取标准数据集中不同类型的IVUS图像分别进行定性和定量分析，以验证算法可行性。实验把每帧图像提取的极值区域和拟合后的最终轮廓同时显示在原始图像上。图 5为选取9幅包含伪影和不包含伪影的IVUS图像提取的管腔极值区域轮廓和椭圆拟合化后的最终轮廓；图 6同样选取图 5中的9幅IVUS图像，描绘出介质极值区域轮廓和椭圆拟合化后的最终轮廓。其中蓝色曲线代表提取的极值区域的轮廓，白色曲线是椭圆拟合化之后的最终轮廓。另外，图 7是本文方法的实验结果和数据集提供的专家标注的轮廓之间的对比，其中白色曲线为标准的内、外膜轮廓，红色曲线为本文方法提取的内、外膜轮廓。图 5-图 7的(a)-(b)无伪影，(c)-(i)均有不同程度伪影。

可以看出，即使IVUS图像在外膜以外区域存在大量伪影，极值区域边缘与最终内外膜边缘仍然相当接近。这是由于虽然外膜以外的部分受到严重干扰，但极值区域仍然能够被检测到，这受益于极值区域选择过程中的面积筛选，能基本过滤掉不符合内外膜区域面积的极值区域，最大程度降低了外膜以外的区域质量所带来的影响，在一定程度上增加了算法的鲁棒性。另外，可以看出，外膜极值区域的边缘变化小于内膜，这是由于管腔区域内图像整体灰度值偏低，内膜边缘不清晰，极值区域边缘定位浮动大。但是经过椭圆拟合化之后，内膜边缘就能较为准确地被定位。因此，椭圆拟合化的主要作用有两点：1)边缘不清晰或部分边缘缺失的IVUS图像提取极值区域产生的误差在该步得到弥补；2)提取的极值区域具有凹凸不平的边界，这与实际管腔内膜和外膜不符，椭圆拟合可以平滑边缘。

本文使用DC(dice coefficient)、JI(jaccard index)、PAD(percentage of area difference)指标以及HD(hausdorff distance) 4个指标定量评估算法准确性。设${A_{{\rm{auto}}}}$是本文方法获得的轮廓包围的面积，${A_{{\rm{man}}}}$是数据集中人工标注轮廓包围的面积，则这4个指标如下所述。

1) DC系数。DC系数是两个区域相似性的度量，该指标的定义为

$ D = \frac{{2\left| {{A_{{\rm{auto}}}} \cap {A_{{\rm{man}}}}} \right|}}{{\left| {{A_{{\rm{auto}}}}} \right| + \left| {{A_{{\rm{man}}}}} \right|}} $

(20)

2) JI指标。JI指标也被称为交并比，用于比较本文自动分割的区域面积和人工标注的轮廓包围的面积的相似性，即

$ J = \frac{{\left| {{A_{{\rm{auto}}}} \cap {A_{\rm{man}}}} \right|}}{{\left| {{A_{{\rm{auto}}}} \cup {A_{\rm{man}}}} \right|}} $

(21)

3) PAD指标。PAD指标是区域差异性的度量，定义为

$ P = \frac{{\left| {{A_{{\rm{auto}}}} - {A_{\rm{man}}}} \right|}}{{{A_{\rm{man}}}}} $

(22)

4) HD距离。HD距离度量了两个点集间的最大不匹配程度：假设本文自动分割方法区域轮廓坐标集合${\mathit{\boldsymbol{C}}_{{\rm{auto }}}} = \left\{ {{a_1}, {a_2}, \cdots, {a_p}} \right\}$，数据集中专家标注的轮廓坐标集合${\mathit{\boldsymbol{C}}_{{\rm{man }}}} = \left\{ {{b_1}, {b_2}, \cdots, {b_p}} \right\}$，则这两个集合之间的HD距离定义为

$ H\left({{\mathit{\boldsymbol{C}}_{{\rm{auto}}}}, {\mathit{\boldsymbol{C}}_{\rm{man}}}} \right) = \max \left({h\left({{\mathit{\boldsymbol{C}}_{{\rm{auto}}}}, {\mathit{\boldsymbol{C}}_{\rm{man}}}} \right), h\left({{\mathit{\boldsymbol{C}}_{\rm{man}}}, {\mathit{\boldsymbol{C}}_{{\rm{auto}}}}} \right)} \right. $

(23)

$ h\left({{\mathit{\boldsymbol{C}}_{{\rm{auto}}}}, {\mathit{\boldsymbol{C}}_{\rm{man}}}} \right) = \mathop {\max }\limits_{a \in {\mathit{\boldsymbol{C}}_{{\rm{auto}}}}} \mathop {\min }\limits_{b \in {\mathit{\boldsymbol{C}}_{{\rm{man}}}}} \left\| {a - b} \right\| $

(24)

$ h\left({{\mathit{\boldsymbol{C}}_{{\rm{man}}}}, {\mathit{\boldsymbol{C}}_{{\rm{auto}}}}} \right) = \mathop {\max }\limits_{b \in {\mathit{\boldsymbol{C}}_{{\rm{man}}}}} \mathop {\min }\limits_{a \in {\mathit{\boldsymbol{C}}_{{\rm{auto}}}}} \left\| {b - a} \right\| $

(25)

式中，||·||为两点之间的欧氏距离。

首先，从标准数据集中不同专家描绘的可变性来看，将本文全自动方法提取的结果分别与两位专家手动描绘的结果比较，以及两位专家之间的比较，计算的DC系数、JI指标、PAD指标、HD距离如表 1所示。若两位专家手动描绘的轮廓完全相同，DC系数和JI指标应该为1，PAD指标和HD距离应该为0。表 1中专家1与专家2手动描绘的DC系数和JI指标不等于1，PAD指标和HD距离不等于0，所以二者手动描绘的结果有所差异，而且即使同一位专家前后两次描绘的结果也不完全相同。而本文方法与专家1以及本文方法与专家2测量的内膜的4个指标分别相差0.01、0.01、0.03、0.03 mm，外膜分别相差0.01、0.03、0.06、0.08 mm。说明本文方法与专家1和专家2手动描绘的结果具有一致性。

其次，根据测试图像是否存在伪影以及伪影的种类，表 2对标准数据集的鲁棒性和泛化测试结果进行了展示，并绘出了表 2的柱状图，如图 8所示。当图像中存在伪影时，各项指标的整体性能明显低于无伪像时的整体性能。由表 2和图 8可知，本文方法对分岔伪影最为敏感，此种情况下性能最低。在无伪影时，所提方法内膜的HD距离达到0.22 mm，外膜的HD距离达到0.23 mm。除此之外，表 2和图 8中数据显示，除分岔伪像外(内膜HD距离为0.37 mm，外膜HD距离为0.45 mm)，其他类型的伪像和无伪像时的数据相差不大，侧血管和阴影伪影的内外膜平均HD距离均不超过0.30 mm，因此本文算法鲁棒性较好。另外，本文在极值区域筛选算法中运用的LBP特征是纹理信息，能有效地消除亮度对图像的影响，所以一定程度上缓解了亮度变化的问题。灰度差异是利用图像灰度值这一先验信息计算区域和其邻域的灰度差。另外，区域周长是极值区域边界变化的度量，利用该信息能够得到边界变化较为稳定的区域。这3种图像信息的使用均一定程度上提高了算法的鲁棒性。

最后，表 3定量显示了本文方法和5个现有方法的比较结果，其柱状图如图 9所示。本文方法检测的内膜的DC系数为0.94±0.02，JI指标为0.90±0.04，PAD指标为0.05±0.05，HD距离为0.28±0.14 mm，外膜的DC系数为0.91±0.07，JI指标为0.87±0.11，PAD指标为0.11±0.11，HD距离为0.41±0.31 mm。虽然已有算法在不断地改进，效果也越来越好，但本文方法与这些方法相比仍然具有明显的优势，比如，与Yang等人(2018)的算法(内膜平HD距离为0.26±0.25 mm，外膜HD距离为0.48±0.44 mm)相比，本文内膜的平均HD距离(0.28 mm)虽然多0.02 mm，但是方差比较小，结果更加稳定；外膜的平均HD距离(0.41 mm)比Yang等人(2018)的算法少0.07 mm，且方差较小。综合来看，本文算法优于Yang等人(2018)的算法。

此外，图 10随机选取了9帧北京安贞医院跨数据集的不同类型的图像进行展示，其中红色曲线为本文方法提取的内外膜轮廓，白色曲线为临床专家手动勾画的内外膜轮廓。结果表明，本文方法提取的内外膜轮廓和临床专家的手动提取结果存在很大相似性。所以本文方法在临床数据中也能得到较好的分割结果，进一步验证了本文方法的泛化能力。

尽管本文方法总体性能较好，但仍有一些局限性：由于极值区域是相互嵌套的，因此极值区域筛选矢量是严格递增的，选取第2个峰值作为管腔区域，倒数第2个或最后1个峰值作为介质区域，但在这些峰值之外，可能还有更加接近真实内外膜的峰值，因此极值区域筛选标准矢量还有待改进。

在之后的工作中，找到一个统一的极值区域标准让结果更精确并使其对每一帧IVUS图像都适用仍是未来研究工作的重点。

IVUS图像内外膜的分割是快速识别和量化动脉粥样硬化的重要步骤，在传统方法中，IVUS图像分割只能分别设计算法串行提取内外膜，而本文算法摒弃复杂的数学建模和网络训练过程，将IVUS图像分割放宽到区域检测，并设计筛选标准矢量，将筛选标准矢量的第2个峰值对应的极值区域作为管腔，倒数第2个或最后1个峰值对应的极值区域作为介质，实现并行分割IVUS图像。本文算法的主要特色有：1)整个过程不需要任何用户交互，而且当前帧提取的结果不依赖于其他帧，只需对单帧图像操作，以保证其他帧提取的结果不会影响到当前帧；2)将图像边缘信息点作为先验信息，并利用区域的纹理信息、灰度差异度和边缘长度等信息设计筛选矢量，实现并行提取；3)相比其他多数方法，本文方法在保证精度的同时，鲁棒性较好。