CN110717888B - 一种血管内光学相干层析成像血管壁内轮廓自动识别方法 - Google Patents

Abstract

本发明涉及一种血管内光学相干层析成像血管壁内轮廓自动识别方法，包括下列步骤：(1)对原始图像Ori进行预处理；(2)血管壁内轮廓识别，方法如下：第一步：腐蚀运算；第二步：面积排序；第三步：形态学方向与面积检测；第四步：导丝存在与否条件判断；第五步：导丝识别纠正；第六步：重建残缺轮廓。

Description

一种血管内光学相干层析成像血管壁内轮廓自动识别方法

技术领域

本发明涉及内窥光学相干层析成像(Intravascular Optical CoherenceTomography,IV-OCT)覆盖心血管、脑血管等内窥OCT成像领域，尤其涉及一种基于形态学方向特征提取的IV-OCT血管图像的管腔轮廓的全自动识别方法。

背景技术

血管内窥光学相干层析成像技术(IV-OCT)是一种基于导管的检查方法，它利用近黑外光，可获取分辨率为十倍于血管内超声(IVUS)的体内血管壁微结构的高分辨率成像。血管管腔轮廓信息对评价血管的狭窄程度，获取支架植入的最佳位置，评价支架的贴壁状况都有很大的应用价值。因此，血管管腔轮廓的准确识别及识别在冠状动脉疾病的治疗中显得尤为重要。

目前，血管管腔轮廓的识别及识别方法主要可以分为半自动识别和全自动识别。而全自动识别相对于半自动识别来说更为快速和便捷。近年来，不断有学者提出全自动识别血管管腔轮廓的方法，其中：Bourezak等人通过利用A-scan的强度图找到组织与背景有差异的地方，再利用组织的厚度信息确定出内边界的轮廓；Macedo等人通过Otsu二值化和数学形态学腐蚀膨胀运算相结合的方法划定感兴趣区域获取血管管腔边界；Zhao Wang等人通过动态规划算法实现了血管管腔及导丝的自动识别；Junedh等人通过利用K均值聚类的方法实现了管腔的轮廓检测及支架的识别；Kenji Sihan等人利用固定阈值的Canny算子对管腔进行成功识别；另外还有利用水平集等主动轮廓以及深度学习的方法识别管腔轮廓的相关报道。

上述方法对于较为规则的血管轮廓的识别效果较为显著，但对于血管轮廓不规则，导管与导丝距离血管管腔距离过近的IV-OCT血管管腔轮廓的识别和提取，出现受环境噪声干扰或者迭代次数过多，过程繁琐的问题。Huishuo Zhao等人提出了基于形态学特征的血管管腔轮廓的识别和提取方法，但该方法中，由于采用形态学闭运算造成容易造成导管、导丝与血管壁相连，在相连时无法对血管管腔轮廓识别；无法处理存在残余血液的血管管腔轮廓识别；难以对不规则管腔轮廓进行识别。基于上述问题，本专利提出了一种血管内光学相干层析成像血管壁内轮廓自动识别方法，相比于其他血管轮廓的识别方法，本方法的优点主要在于能够正确提取出低信噪比、形状不规则、受残余血液、导丝和导管影响较为严重而难以识别和提取血管轮廓。

发明内容

本发明的目的是提供一种IV-OCT血管图像的管腔轮廓全自动识别方法，以克服现有的IV-OCT血管图像的管腔轮廓识别和分割的过程中，因导丝和导管距离血管管腔过近、残余血液影响或者管腔轮廓过于复杂而提取的管腔轮廓出现识别困难的问题。本发明在整个识别过程使用了图像对比度增强、二值化、形态学开运算、连通域方向与面积特征筛选、动态切割方法等步骤能够正确提取出低信噪比、形状不规则、受残余血液、导丝和导管影响较为严重的血管轮廓。技术方案如下：

一种血管内光学相干层析成像血管壁内轮廓自动识别方法，所处理的原始图像为IVOCT极坐标下图像，记为Ori，图像左下角为极坐标零点，其中m为图像行数，n为图像列数，包括下列步骤：

(1)对原始图像Ori进行处理，方法如下：

第一步：图像对比度增强；

第二步：二值化处理，采用最佳阈值法确定二值化阈值。

将图像目标和背景点分布看成正态分布，选取阈值T的初始估计值，用T分割图像，得到两组像素，分别计算这两组像素的平均灰度值u1和u2，取新的T值为u1与u2的平均值，再重复以上步骤，迭代足够多的次数，直到逐次所得的T值之差小于10^-4，以此T值作为最佳阈值进行二值化处理，得到二值化图像BW；

第三步：进行形态学开运算；

选用半径为1的圆盘状结构元素se,对BW进行形态学开运算，去除分散的噪声，得到预处理后的结果图用PRE表示，选用大小为[5*5]的矩形结构元素SE,对PRE进行形态学腐蚀运算，得到腐蚀后的结果图用ER表示；

(2)血管壁内轮廓识别，方法如下：

第一步：面积排序

将ER中的所有连通区域按照面积从大到小的顺序进行排序，取前5个，得到面积排序后的图像S；

第二步：形态学方向与面积检测

对S中的所有连通区域进行标记，得到所有连通区域的标记矩阵L，从中选取连通区域的方向和面积属性，设定方向阈值和面积阈值，将方向属性大于方向阈值并且面积小于阈值的连通区域的像素值置为0，将经过方向与面积检测后的图像记为I；基于IVOCT图像噪声水平，选取方向阈值为89.5°，面积阈值0.05mm²。

第三步：动态切割去噪

对图像I中经过的残余血液噪声进行切割去除，方法如下：设n为图像列数，切割阈值th＝n/10，对该阈值周围的像素分布情况对自身进行动态调整，将矩阵列数小于该切割阈值的所有像素置为0，将处理后图像记为Inone。

第四步：导丝识别纠正

将图像Inone在0～n/5的列坐标范围内按照行数进行遍历，若同一行中的连续为1的像素数目小于20，则将这些像素置为0，从而去除导丝，同样将纠正后的图像记为IC；

第五步：第二次形态学开运算

再次对图像IC进行形态开操作，去除所有噪声，将修正后的图像记为bw。

第六步：重建残缺轮廓；

对于经过导丝识别纠正的图像bw，取每一行中第一个为1的点作为轮廓的边缘点，找到被狭长阴影分割的上下几部分的轮廓；使用双立方插值，重构出狭长阴影部分的轮廓，将提取出的完整轮廓记为Img；

第七步：误去除信号的寻回：

对于有可能受到开运算时误去除影响的轮廓Img，取图像I中在此轮廓附近的每一个连通域进行填充，得到的图像记为I2。

第八步：通过几何中心去除寻回噪声

通过噪声与信号部分不同的几何中心特征，消除掉其中的噪声，得到的图像记为I3，随后使用图像I3重复第六步，实现更完善的轮廓识别，该轮廓记为Img2；

(3)坐标变换

利用提取出的轮廓Img2对原始图像Ori进行分割，对分割后的图像进行坐标变换，得到笛卡尔坐标系下的输出图像Re。

与现有技术相比，本发明的有益效果是：

本发明克服了由于血管轮廓不规则，导管与导丝距离或血管管腔过近导致IVOCT血管图像中血管管腔轮廓难以识别和提取的问题。在不设定感兴趣区域，不对A-scan进行强度对比的情况下，使用了迭代次数较少的动态切割去噪与基于面积排序和形态学方向与面积特征提取的轮廓识别方法，通过二次识别更完善地去除噪声，修正结果，从而能够正确提取出低信噪比、形状不规则、受残余血液、导丝和导管影响较为严重而的血管轮廓。

附图说明

图1是一种IVOCT血管图像的管腔轮廓全自动识别方法的流程图；

图2经过图像对比度增强处理后IVOCT极坐标下图像

图3二值化后IVOCT极坐标下图像

图4(a)和(b)简单切割处理后可能会产生的IVOCT极坐标下图像

图5动态切割处理后IVOCT极坐标下图像

图6第二次开运算处理后IVOCT极坐标下图像

图7第一次轮廓重建后IVOCT极坐标图像(黑线表示血管轮廓)

图8为寻回误去除部分的IVOCT极坐标图像(黑点标识几何中心)

图9为利用几何中心去除其他噪声后IVOCT极坐标图像(黑点标识几何中心)

图10进行第二次轮廓重建结果(黑线表示血管轮廓)

图11坐标变换后笛卡尔坐标下血管轮廓识别结果(黑线表示血管轮廓)

具体实施方式

下面结合附图和具体实施例对本发明技术方案作进一步详细描述，所描述的具体实施例仅对本发明进行解释说明，并不用以限制本发明。

如图1所示，本发明提出的一种IVOCT血管图像的管腔轮廓全自动识别方法包括以下步骤：

本发明提供了一种基于形态学方向与面积特征提取的IV-OCT血管图像的管腔轮廓全自动识别方法，本发明克服了现有的IV-OCT血管图像的管腔轮廓识别和分割的过程中，因导丝和导管距离血管管腔过近，或者管腔轮廓过于复杂而提取的管腔轮廓出现识别困难的问题，本方法首先对原始极坐标系下的IV-OCT图像进行预处理；其次，对预处理后的图像进行形态学腐蚀运算，使用matlab中的工具箱函数regionprops，基于形态学的方向和面积特征对原始图像中的各个对象进行标记得到各个连通区域的方向及面积属性，去除导丝以及导管；然后通过双立方差值的方法重构了由于导丝的遮挡而缺失的血管管腔轮廓；最后通过坐标变换得到了笛卡尔坐标系下的输出图像。

为了解决上述技术问题，本发明提出了一种基于形态学方向特征提取的IV-OCT血管图像的管腔轮廓自动识别方法，其特征在于整个识别过程使用了形态学方向及面积特征提取运算，方法简单；并能够正确提取出形状不规则、受导丝和导管影响较为严重而难以识别和提取的血管轮廓，设所处理图像为IVOCT极坐标下图像，图像左下角为极坐标零点，其中m为图像行数，n为图像列数，所述方法包括以下步骤：

(1)对原始图像Ori进行对比度增强处理，方法如下：

步骤一：为了尽可能保留图像Ori的原始信息，将极坐标系下的16位的原始图像Ori进行对比度增强，将16位的图像转为无符号整型8位的图像，以使图像进行清晰的显示，方便后期的管腔轮廓的识别和提取；

步骤二：将中值滤波后的图像使用大津法进行二值化处理，具体步骤如下：

第一步，将图像目标和背景点分布看成正态分布，求图像的最大灰度A与最小灰度B，令阈值T的初始估计值为二者平均；

第二步，用T分割图像，得到两组像素：G1为所有灰度值>T的像素组成，G2为所有灰度值<T的像素组成；

第三步，对区域G1和G2中的所有像素计算平均灰度值u1和u2。

第四步，取新的T值为u1与u2的平均值。

第五步，重复第二步到第四步，直到逐次迭代所得的T值之差小于10^-4，以此T值作为最佳阈值进行二值化处理，得到二值化图像BW，便于后期导丝、导管的去除以及轮廓的寻找和提取。

步骤三：进行形态学开运算；

选用半径为1的圆盘状结构元素se,对BW进行形态学开运算，去除分散的噪声。得到预处理后的结果图用PRE表示；

(2)血管壁内轮廓识别，方法如下：

步骤一：提出了面积排序的方法，其主要思想为：对腐蚀运算后的二值化图像ER进行形态学标记，得到不同的连通区域，由于血管部分占绝大部分，因此被标记的连通区域中面积大的区域属于血管轮廓；而连通区域面积较小的部分则属于导管及导丝。通过统计分析，我们发现，这些图片的面积较大的血管的连通区域的个数不超过5个，因此我们将这些血管连通区域部分按照面积从大到小的顺序进行排序，取前5个，得到面积排序后的图像S。

步骤二：通过统计分析可得，导管及导丝主要位于极坐标图像的左半部分，形状比较规则(呈现直线型分布)，连通区域比较小；而血管组织信息主要位于图像的中部及右半部分，形状和轮廓呈现出不规则性，这种情况下导管和导丝不能通过直接选取有效区域的方法进行去除。因此我们提出了使用形态学方向与面积特征提取的方法，对导管以及导丝进行有效的去除。

通过使用matlab中的工具箱函数bwlabel标记面积排序后的图像S中的各个连通区域，得到包括各个目标的方向和面积结构属性的标记矩阵L。识别出的各个目标的方向具有不同性。其中，形状为直线型的导管的方向主要集中在90°左右，而血管管腔的方向比较分散，且主要集中在0-80°范围内。因此我们通过设定合理的阈值即可将导管和血管进行合理的区分，本专利选用的方向阈值为89.5°，即连通区域的方向若大于设定的方向阈值，便认定为导管，反之，则将此连通区域确定为血管。经过统计分析，发现有3个导管也满足大于方向阈值的条件，因此这三个血管的连通区域会被误判为导管。

为进一步降低错误率，对经过方向检测后的连通区域继续进行面积检测。可以看到小于面积阈值(0.05mm²)的都是导管，并且在面积检测中的3个被误判的血管的面积大于面积阈值。因此，根据方向检测的错误率进一步被降低。将方向属性大于方向阈值并且面积阈值小于0.05mm²的连通区域的像素值置为0，将经过方向与面积检测后的图像记为I。

步骤三：经过图像的预处理和形态学方向和面积检测之后，图像剩余噪声基本由于残余血液引起，且成为大面积的连通域形式。动态切割去噪是根据噪声的位置特点(靠近导管)进行去噪的方法。简单的切割去噪是在图像的列数中选定一个阈值(设图像总列数为n，约在n/10～n/5之间)，将该列之前的所有像素置为0。这种方法实现起来极其简单，但很多噪声并非竖直分布，有一定起伏，几乎不可能去除干净，如果阈值设置较大(约n/5)，还很有可能误去除信号部分。

动态切割方法是通过设置阈值为n/10，再以这个阈值为参考，对该行像素p的8邻域依次进行检索。当检索到的这些像素点存在1或已经通过检索的像素p存在，继续检索；当检索到的像素点全部为0，停止检索并将该行中小于检索终点列坐标的所有像素点置0。这种方法可以去除起伏不定的噪声连通域，并且在到达噪声与信号的分隔时自动停止。且原理简单，易于实现。将处理后的图像记为Inone。下图2为对比度增强后IVOCT极坐标下图像，在图2中可以明显观察到残留血液产生的大片连通域。图3为二值化后方法处理前IVOCT极坐标下图像，图4为简单切割处理后可能会产生的IVOCT极坐标下图像，图5为动态切割处理后IVOCT极坐标下图像。从图3中看出有大量残留血液和导管，通过切割方法可以去掉残留血液和导管，但简单切割方法如图4所示，找不到合适的阈值完全去除噪声并保留信号。采用动态切割的方法可以保留血管组织如图5所示，仍有部分小斑块噪声未被去除，并且由于信噪比不高，部分血管组织被去掉，如图2中线框所示部分。这些问题将在之后的步骤中被修正。

步骤四：导丝识别纠正

由于内窥OCT在进行冠状动脉的疾病检查时需要有导丝的指引，然而导丝的反射太强，在强度图中表现为同一个深度方向上导丝部位的强度值高，随着深度的增加，强度值突变近乎为0的值，呈现出梯度突变的情况。经过之前的处理，依然可能有不连通的导丝反光残留，导丝反光主要集中在图像列坐标的0～n/5附近，具有细长的特征，将图像Inone在0～n/5的列坐标范围内按照行数进行遍历，若同一行中的连续为1的像素数目小于20，则将这些像素置为0，从而去除导丝，同样将纠正后的图像记为IC；

步骤五：第二次形态学开运算

经过之前的处理，仍残留的噪声基本呈现小斑块形态，在绝大多数情况下，面积要小于信号区域，故设置开运算阈值为800pixel，再次进行形态开操作，去除所有噪声。将修正后的图像记为bw。至此去噪部分完成，如果出现误去除和去除失败等情况，在随后的步骤中依然可以修正；图6为第二次开运算结果。

步骤六：重构残缺轮廓，由于被导丝遮挡部分或信噪比低造成二值化有效信号被滤除的血管组织信息不能观察到，在OCT图像中会出现一条由于导丝的遮挡引起的一条长长的暗线。对后期支架的置入位置的确定及后期FFR的计算造成很大的影响。因此，填充出残缺的血管轮廓显得至关重要。

为了找到残缺的血管轮廓，通过分析预处理及初步形态学填充孔洞操作的图像可知，血管的有效信息大致可以分为三个部分，被导丝遮挡的部分，以及位于被导丝遮挡的部分的上下两个方向的部分。对于经过导丝识别纠正的图像bw，取每一行中第一个为1的点作为轮廓的边缘点，找到被狭长阴影分割的上下几部分的轮廓；使用双立方插值，重构出狭长阴影部分的轮廓，将提取出的完整轮廓记为Img。

步骤七：由于图像信噪比不高，信号强度不够，二值化时不能兼顾较高的信号强度与较低的噪声强度，故较弱的信号会发生较小程度的丢失，使血管壁区域出现间断。这些间断会分割本是一体的血管壁，分割出的较小区域会被处理过程中的图像开运算误去除。并且，使用上述动态切割去噪后，对于不连通的较大块噪声，依然会发生无法去除的现象。图7为采用之前算法的一次识别轮廓：

为了更好地解决以上问题，本算法使用了第一次轮廓重建后不完善的轮廓识别结果，之后根据识别到的轮廓与图像中连通域的相互关系补足缺漏。

得到不完善的轮廓后，根据识别出的轮廓，保留第二次开运算以前图像中轮廓附近的全部连通域，将不在轮廓附近的连通域全部置为0，得到的图像记为I2。这一步的主要目的是将之前步骤的误去除信号寻回；图8为寻回误去除部分的二值图像。

步骤八：在步骤七的处理中，由于不完善的轮廓可能经过噪声区域，也会将一部分噪声寻回，本算法利用噪声与信号的不同位置特征将寻回噪声去除，因噪声一般集中在导丝附近，即分布在图像列坐标的0～n/5，而有效信号大多分布在图像列坐标的n/5～n。首先，计算图像I2保留的所有k个连通域的几何中心：

其中，

p,q＝0,1,2,···

x_c和y_c分别是图像几何中心坐标，

M₀₀(p＝q＝0)为连通域的零阶矩；

M₁₀(p＝1,q＝0),M₀₁(p＝0,q＝1)为连通域的一阶矩；

将这些几何中心按照列坐标进行从大到小的排序，计算其中前k/2(向上取整)个连通域的几何中心列坐标平均值，并将每个连通域的几何中心与此值相比较，将所得差值大于120pixel的连通域置为0，得到的图像记为I3，这一步的主要目的是去除寻回误去除部分后产生的其他噪声。图9为利用几何中心去除其他噪声后的二值图像，几何中心已在图8与图9中标出。

处理之后，使用I2重复步骤六，再次识别出较完善的轮廓，记为Img2。效果如图10所示：

(3)坐标变换

将提取出的完整轮廓记为Img2对原始图像Ori进行分割，对分割后的图像进行坐标变换，得到笛卡尔坐标系下的输出图像Re。双线性插值的基本原理为：

对于N*M的极坐标图像，在与之对应于直角坐标下的目标显示图像中(大小为2N*2N)，每个像素点(x,y),对于图像扫描中心(N，N)也就是坐标零点需要转换为距离零点的距离和与初始角度的夹角，其对应关系为：

双线性内插法是通过极坐标图像中邻近的2*2像素点区域来计算x*y的值，其计算方法为：

f(x,y)＝(1-Δρ)(1-Δθ)f([ρ],[θ])+Δρ(1-Δθ)f([ρ]+1,[θ])

+Δθ(1-Δρ)f([ρ],[θ]+1)+ΔρΔθf([ρ]+1,[θ]+1)

(6)

经过本发明处理后的输出图像与未经过处理的原始图像相比，图像中血管的管腔轮廓得到正确的识别和提取。

IV-OCT血管图像的管腔轮廓的正确识别及对于评价心血管支架手术后血管的再狭窄，研究动脉粥样硬化特别有用，血管图像的轮廓识别是治疗冠状动脉疾病非常重要的基础性工作。现有动态规划，主动轮廓模型，深度学习等方法可以实现IV-OCT血管轮廓的识别和提取。但在识别和提取轮廓的过程中，可能会受到迭代次数过多，参数设置不合理等问题的限制。本专利提出整个识别过程使用了图像对比度增强、最佳阈值二值化、形态学开运算、连通域方向与面积特征筛选、动态切割方法等步骤能够正确提取出低信噪比、形状不规则、受残余血液、导丝和导管影响较为严重的血管轮廓。尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。

Claims (2)

1.一种血管内光学相干层析成像血管壁内轮廓自动识别方法，需要进行血管壁内轮廓识别的原始图像为IVOCT极坐标下图像，记为Ori，图像左下角为极坐标零点，其中m为图像行数，n为图像列数，包括下列步骤：

（1）对原始图像Ori进行处理，方法如下：

第一步：图像对比度增强；

第二步：二值化处理，采用最佳阈值法确定二值化阈值；

将图像目标和背景点分布看成正态分布，选取阈值T的初始估计值，用T分割图像，得到两组像素，分别计算这两组像素的平均灰度值u1和u2，取新的T值为u1与u2的平均值，再重复以上步骤，迭代足够多的次数，直到逐次所得的T值之差小于10^-4，以此T值作为最佳阈值进行二值化处理，得到二值化图像BW；

第三步：进行形态学开运算；

选用半径为1的圆盘状结构元素se, 对BW进行形态学开运算，去除分散的噪声，得到预处理后的结果图用PRE表示，选用大小为[5*5]的矩形结构元素SE, 对PRE进行形态学腐蚀运算，得到腐蚀后的结果图用ER表示；

（2）血管壁内轮廓识别，方法如下：

第一步：面积排序

将ER中的所有连通区域按照面积从大到小的顺序进行排序，取前5个，得到面积排序后的图像S；

第二步：形态学方向与面积检测

使用matlab中的函数bwlabel标记面积排序后的图像S中的各个连通区域，得到标记矩阵L，使用regionprops函数选取连通区域的方向和面积属性,设定方向阈值和面积阈值，将方向属性大于方向阈值并且面积小于阈值的连通区域的像素值置为0，将经过方向与面积检测后的图像记为I；基于IVOCT图像噪声水平；

第三步：动态切割去噪

动态确定每行的切割阈值，以该切割阈值作为切割划分，对图像I中在阈值列左侧的残余血液噪声进行切割去除，方法如下：设n为图像列数，切割阈值列数th=n/10，以每一行阈值列数th_i为起点，对该列数右侧的像素的分布情况进行检索，以检索结果对每行阈值自身进行动态调整，将每行切割阈值列数左侧的所有像素置为0，将处理后图像记为Inone；

第四步：导丝识别纠正

将图像Inone在0~n/5的列坐标范围内按照行数进行遍历，若同一行中的连续为1的像素数目小于20，则将这些像素置为0，从而去除导丝，同样将纠正后的图像记为IC；

第五步：第二次形态学开运算

再次对图像IC进行形态开操作，去除所有噪声，将修正后的图像记为bw；

第六步：重建残缺轮廓；

对于经过导丝识别纠正的图像bw，取每一行中第一个为1的点作为轮廓的边缘点，找到被狭长阴影分割的上下几部分的轮廓；使用双立方插值，重构出狭长阴影部分的轮廓，将提取出的完整轮廓记为Img；

第七步：误去除信号的寻回：

对于有可能受到开运算时误去除影响的轮廓Img，取图像I中在此轮廓附近的每一个连通域进行填充，得到的图像记为I2；

第八步：通过几何中心去除寻回噪声

通过噪声与信号部分不同的几何中心特征，消除掉其中的噪声，得到的图像记为I3，随后使用图像I3重复第六步，实现更完善的轮廓识别，该轮廓记为Img2；

（3）坐标变换

利用提取出的轮廓Img2对原始图像Ori进行分割，对分割后的图像进行坐标变换，得到笛卡尔坐标系下的输出图像Re。

2.根据权利要求1所述的方法，其特征在于，选取方向阈值为89.5°，面积阈值0.05mm²。

Images