CN103606152A - 基于sift特征点聚类及布尔差运算的dsa血管图像分割方法 - Google Patents

Abstract

本发明属于医学图像处理领域，特别涉及一种基于SIFT特征点聚类及布尔差运算的DSA血管图像分割方法。该方法包括下述步骤：输入相对应的DSA蒙片图像和活片图像；在相同阈值条件下，利用SIFT算法对蒙片图像和活片图像进行几何特征点的提取；采用基于灰度梯度值的方法进行几何特征点局部位置的调整；进行基于欧氏距离的几何特征点聚类及布尔差运算；采用基于Otsu阈值的方法进行血管图像分割。利用该方法得到患者DSA减影中血管图像的过程，具有鲁棒性较强、计算过程方便等优点，能够为基于DSA介入治疗的临床手术提供准确的血管图像数据。

Description

基于SIFT特征点聚类及布尔差运算的DSA血管图像分割方法

技术领域

本发明属于医学图像处理领域，特别涉及基于SIFT特征点聚类及布尔差运算的DSA血管图像分割方法。

背景技术

数字减影（Digital Subtraction Angiography，DSA）技术是一种使X光射线序列图片中的血管可视化的强大技术，是常规血管造影术和电子计算机图像处理技术相结合的产物，已在临床应用20多年，是血管疾病无创诊断与介入治疗手术导航的重要依据，广泛应用于X光射线序列成像中血管的可视化系统中。其基本原理简言之为：X光射线穿过人体各解剖结构形成荧光影像，经过影像增强器增强后为电视摄像管采集而形成视频影像；再经对数增幅和模/数转换形成数字影像；这些数字信息输入计算机处理后，再经蒙片图像（打入造影剂前的图像）和活片图像（打入造影剂后的图像）的减影、对比度增强和数/模转换，产生数字减影图像。DSA技术主要用于血管疾病的诊断、术中观察，如：颈段动脉狭窄或闭塞、动静脉畸形、颅内动脉瘤、动脉闭塞和血管发育异常等，以及颅内肿瘤供血动脉的观察及定性定位诊断等。

但是，由于患者头部微小转动及吞咽等动作的存在，致使蒙片图像和活片图像存在空间上的微小错位，导致直接减影后的图像中存在大量的椒盐噪声及运动伪影，严重影响了血管病变诊断与治疗的准确性。因此，根据蒙片图像和活片图像的数据信息，尽可能地提取血管组织的准确图像已成为当前基于DSA技术的血管诊疗中亟待解决的问题。

针对基于DSA技术的血管病变诊疗中血管图像的提取问题，国内学者提出了一些解决方法，发表的文献主要包括：《深圳大学学报》的“基于拉伸空间变换算法在DSA配准中的应用”；《中国医疗器械杂志》的“基于仿射变换模型的外周DSA的自动配准”。此类技术主要利用自动栅格线或者角点检测的方法提取DSA图像中的标记点，然后利用线型/非线性的方法进行空间变换，进而将蒙片图像和活片图像进行配准，最后提取出血管图像。此类方法的血管图像精度相比传统的减影方法有所提高，但其计算过程较为繁琐，速度较慢。另外，当蒙片图像和活片图像之间的位移变化较大时，此类方法的效果很差。

国外学者针对DSA图像中血管提取的问题所提出的方法包括：在《ComputerVision and Image Understanding》的“A 3-D space–time motion detection for aninvariant image registration approach in digital subtraction angiography”和《PatternRecognition》的“An invariant approach for image registration in digital subtractionangiography”。此类技术的基本原理同样是从图像的空间变换出发，以蒙片图像和活片图像之间的逆向运动估计作为计算目标，从而实现蒙片图像和活片图像的位置配准，以达到提取血管图像的目的。利用此类方法所得到的减影图像中运动伪影可进一步减少，但其在运动估计过程中的计算稳定性较低，运算过程复杂。同时，对DSA图像原始数据的测试环境也提出了苛刻的要求。

发明内容

本发明的目的正是为了克服现有技术的不足，解决如何准确利用患者蒙片图像和活片图像，从分割角度进行患者血管图像的提取问题，而提出一种基于SIFT特征点聚类及布尔差运算的DSA血管图像分割方法。利用该方法得到患者DSA减影中血管图像的过程，具有鲁棒性较强、计算过程方便等优点，能够为基于DSA介入治疗的临床手术提供准确的血管图像数据。

一种基于SIFT特征点聚类及布尔差运算的DSA血管图像分割方法，所述方法包括以下步骤：

1）将对应的蒙片图像和活片图像进行直接减影，获得一张直接减影图像。在相同阈值条件下，利用SIFT算法分别在相对应的蒙片图像与活片图像中确定几何特征点，同时构造两个几何特征点集，分别为“蒙片图像几何特征点集”和“活片图像几何特征点集”，并将“活片图像几何特征点集”复制映射到直接减影图像上，生成“直接减影图像几何特征点集”；

2）在直接减影图像上根据几何特征点附近的灰度梯度进行几何特征点位置的局部调整，以便将几何特征点移动至直接减影图像中血管的内部，具体步骤如下：

①判断“直接减影图像几何特征点集”是否为空。若非空，则从该特征点集中取出一个几何特征点，转向②步。否则几何特征点位置调整的计算过程结束；

②比较该几何特征点的灰度值与其8邻域内各像素的灰度值的大小。若该几何特征点的灰度值最小，则将该几何特征点从特征点集中移除，并转向①步。否则转向③步；

③计算该几何特征点与Kirsch算子8个模板卷积得到的8个方向上的梯度值，记为

（i为Kirsch算子8个模板中的第i个模板，1≤i≤8），使该几何特征点向梯度值最大的方向

移动一个像素单位，将移动后的几何特征点重新加入特征点集中，转向①步。

3）将“蒙片图像几何特征点集”和“直接减影图像几何特征点集”中的几何特征点根据欧氏距离进行聚类，然后进行几何特征点的类间布尔差运算，以防止处在毛细血管等精细结构上的几何特征点被误减，具体步骤如下：

①针对蒙片图像和直接减影图像分别构造两个几何特征点聚类集合，即“蒙片图像几何特征点聚类集合”和“直接减影图像几何特征点聚类集合”，并初始化这两个聚类集合为空；

②针对“蒙片图像几何特征点集”和“直接减影图像几何特征点集”，分别判断几何特征点集中是否有未被聚类的特征点。若有，则任选一个几何特征点，转向③步。若所有几何特征点都已被聚类处理，则转向④步；

③在蒙片图像/直接减影图像中，依次遍历“蒙片图像几何特征点聚类集合”/“直接减影图像几何特征点聚类集合”中的所有聚类，计算该几何特征点与某个聚类中的某个几何特征点之间的欧氏空间距离dis(p₀,p)＝sqrt((x_p0-x_p)²+(y_p0-y_p)²)，其中p₀,p分别为待聚类几何特征点和当前聚类中的某个几何特征点，x,y分别表示几何特征点的横纵坐标。当这两点间距离小于阈值dis_threshold时，将该几何特征点p₀归入点p所在聚类。若几何特征点p₀与所有聚类中的所有几何特征点都不满足欧氏空间距离小于阈值dis_threshold，则将该几何特征点作为新的聚类加入到“蒙片图像几何特征点聚类集合”/“直接减影图像几何特征点聚类集合”中，转向②步；

④针对“直接减影图像几何特征点聚类集合”中的每一个聚类

计算“蒙片图像几何特征点聚类集合”中所有聚类

与该聚类

重叠的面积，然后再将所有重叠部分面积进行并集运算

（假设直接减影图像中共有W个聚类，蒙片图像中共有T个聚类），计算得到的总面积

与该聚类

的面积

的比值，若该比值大于某一阈值

则删除该聚类，即将

内所有几何特征点从直接减影图像上删去。

4）以类间布尔差运算后的几何特征点作为种子点集合，通过Otsu算法确定分割阈值，以邻域区域生长的方式进行血管图像的分割，具体步骤如下：

①初始化种子点集，将所有种子点标记为未生长状态，并计算当前直接减影图像的Otsu分割阈值；

②判断种子点集中的种子点是否全部为已生长状态。若不是，任取一个种子点，标记为已生长状态，转向③步。否则算法结束；

③考察该种子点8邻域内各像素（忽略种子点集中的点）的灰度值，取灰度值最小的前三个像素，判断它们与Otsu阈值的关系。若小于Otsu阈值，则将这些像素加入到种子点集中，并标记为未生长状态。若大于Otsu阈值，则舍去。转向②步。

附图说明

图1为本发明DSA血管图像分割方法的流程图。

图2为本发明实施例中所涉及的DSA蒙片图像。

图3为本发明实施例中所涉及的DSA活片图像。

图4为本发明实施例中蒙片图像几何特征点集（小圆圈）的示意图。

图5为本发明实施例中直接减影图像几何特征点集（小圆圈）的示意图。

图6为本发明实施例中根据灰度梯度进行特征点局部位置调整的示意图。

图7为本发明实施例中蒙片图像几何特征点聚类集合（横、竖线段范围内）的示意图。

图8为本发明实施例中直接减影图像几何特征点聚类集合（横、竖线段范围内）的示意图。

图9为本发明实施例中直接减影图像经类间布尔差运算后几何特征点集（小圆圈）的示意图。

图10为本发明实施例中利用Otsu阈值进行血管分割后的示意图。

具体实施方式

下面结合附图对本发明的较佳实施例进行详细阐述，以使本发明的优点和特征能够更易于被本领域技术人员理解，从而对本发明的保护范围做出更为清楚明确的界定。

请参考图1，本发明提供一种基于SIFT特征点聚类及布尔差运算的DSA血管图像分割方法，该方法包括下述步骤：输入相对应的DSA蒙片图像和活片图像；在相同阈值条件下，利用SIFT算法对蒙片图像和活片图像进行几何特征点的提取；采用基于灰度梯度值的方法进行几何特征点局部位置的调整；进行基于欧氏距离的几何特征点聚类及布尔差运算；采用基于Otsu阈值的方法进行血管图像分割。

请参考图2～10，下面将具体介绍该基于SIFT特征点聚类及布尔差运算的DSA血管图像分割方法的实施过程。

输入由常规医用DSA机所采集到的病灶区域DSA图像（需保证蒙片图像与活片图像相对应）。

首先，将对应的蒙片图像和活片图像进行直接减影，获得一张直接减影图像，具体步骤为：

因为直接减影图像的文件存储结构和所占内存大小与活片图像相同，所以可以开辟一块同样大小的内存，并拷贝活片图像的数据到该内存块中，以此存储直接减影图像的数据；

将蒙片图像和活片图像的对应像素直接做差，并将每一个像素对的差值结果赋值给已开辟的直接减影图像文件的对应像素，待全部的对应像素做差完毕，即可得到一张直接减影初始化图像。该图像的背景为黑色，而血管区域为灰白色，不符合临床DSA减影图像的使用习惯，因此将该图像每一个像素的灰度值用“255-灰度值”进行重新赋值，这样即可得到一张背景为白色，血管区域为灰黑色的直接减影图像。

其次，在相同阈值条件下，利用SIFT算法分别在相对应的蒙片图像和活片图像上计算并提取若干个几何特征点，同时构造两个几何特征点集，分别为“蒙片图像几何特征点集”和“活片图像几何特征点集”，用于存储已提取到的几何特征点。并将“活片图像几何特征点集”的特征点像素位置坐标复制映射到直接减影图像上，以此生成“直接减影图像几何特征点集”。

然后，在直接减影图像上根据几何特征点附近像素的灰度梯度值进行几何特征点位置的局部调整，以便将“直接减影图像几何特征点集”中的几何特征点移动到直接减影图像中血管的内部区域（因为利用SIFT算法所提取到的几何特征点一般处于角点、边缘等位置，这样不利于后续基于种子点生长的血管区域填充分割方法的实施，容易导致错误分割），具体步骤为：

判断“直接减影图像几何特征点集”是否为空。若为空，说明所有的几何特征点均已被局部调整完毕，则该计算过程结束；若不为空，则从该特征点集中任取出一个几何特征点，比较该几何特征点的灰度值与其周围8个像素的灰度值的大小，若该几何特征点的灰度值最小，则将该几何特征点从特征点集中移除，并重新在特征点集中选取新的几何特征点进行邻域灰度值比较；若该几何特征点的灰度值并非最小，则计算该几何特征点与Kirsch算子8个模板卷积得到的8个方向上的梯度值，记为

（i为Kirsch算子8个模板中的第i个模板，1≤i≤8），令该几何特征点向梯度值最大的方向移动一个像素单位，将移动后的几何特征点重新加入到特征点集中。并继续在特征点集中选取新的几何特征点进行邻域灰度值比较。

而后，为了尽可能萃取“直接减影图像几何特征点集”中的有效几何特征点（血管区域内的几何特征点），剔除无效几何特征点（非血管区域内的几何特征点），需将“蒙片图像几何特征点集”和“直接减影图像几何特征点集”进行做差处理。将当前“蒙片图像几何特征点集”和“直接减影图像几何特征点集”中的几何特征点根据欧氏距离进行聚类，这样便于进行几何特征点的类间布尔差运算，从而防止处在毛细血管等精细结构上的几何特征点被误减，具体步骤为：

针对蒙片图像和直接减影图像分别构造两个几何特征点聚类集合，即“蒙片图像几何特征点聚类集合”和“直接减影图像几何特征点聚类集合”，同时将这两个聚类集合初始化为空；针对当前“蒙片图像几何特征点集”和“直接减影图像几何特征点集”，分别判断其中是否有未被聚类的几何特征点；

若存在未被聚类的几何特征点，则从当前“蒙片图像几何特征点集”/“直接减影图像几何特征点集”中选取一个几何特征点，并依次遍历“蒙片图像几何特征点聚类集合”/“直接减影图像几何特征点聚类集合”中的所有聚类，计算该几何特征点与某个聚类中的某个几何特征点之间的欧氏空间距离dis(p₀,p)＝sqrt((x_p0-x_p)²+(y_p0-y_p)²)，其中p₀,p分别为待聚类几何特征点和当前聚类中的某个几何特征点，x,y分别表示几何特征点的横纵坐标。当这两点间的欧氏空间距离小于阈值dis_threshold（本例中取dis_threshold=20）时，将该几何特征点p₀归入点p所在聚类。如果几何特征点p₀与所有聚类中的所有几何特征点都不满足欧氏空间距离小于阈值dis_threshold（本例中取dis_threshold=20），则将该几何特征点作为新的聚类加入到“蒙片图像几何特征点聚类集合”/“直接减影图像几何特征点聚类集合”中。并从当前“蒙片图像几何特征点集”/“直接减影图像几何特征点集”中重新开始选取新的几何特征点，进行聚类处理；

若所有几何特征点都已被聚类处理，则针对“直接减影图像几何特征点聚类集合”中的每一个聚类

计算“蒙片图像几何特征点聚类集合”中所有聚类

与该聚类

重叠的面积，然后再将所有重叠部分的面积进行并集求和运算

（假设直接减影图像中共有W个聚类，蒙片图像中共有T个聚类），计算得到的总面积与该聚类的面积的比值，如果该比值大于阈值

（本例中取α_threshold=0.3），则删除该聚类，也就是将

内所有几何特征点从直接减影图像上删去。

最后，以类间布尔差运算后的几何特征点作为种子点集合，通过Otsu算法确定分割阈值，以邻域区域生长（亦称作种子点填充）的方式进行血管图像的分割，具体步骤为：

初始化种子点集，将所有种子点标记为未生长状态。假设待分割的血管区域图像为前景图像，其他区域为背景图像，通过穷举搜索灰度值t（取值范围为0～255）而找到前景图像和背景图像之间的方差g(t)的最大值，g(t)＝w₀(t)*(u₀(t)-u)²+w₁(t)*(u₁(t)-u)²（式中各参数含义为：当前景图像和背景图像的分割阈值为灰度值t时，w₁(t)表示前景图像所占整幅图像的像素数量比例，w₀(t)表示背景图像所占整幅图像的像素数量比例，u₁(t)表示前景图像所有像素的灰度值平均值，u₀(t)表示背景图像所有像素的灰度值平均值，u表示整幅图像中所有像素的灰度平均值），当g(t)最大时所对应的灰度值t即为当前直接减影图像中待分割血管区域的Otsu分割阈值；

判断种子点集中的种子点是否全部为已生长状态。如果是，则结束整个血管图像分割过程；如果不是，则从种子点集中任取一个种子点，并标记为已生长状态，考察种子点8邻域内各像素（忽略种子点集中的点）的灰度值，取灰度值最小的前三个像素，判断它们与Otsu阈值t的关系；如果小于Otsu阈值t，则将这些像素加入到种子点集中，并标记为未生长状态；如果大于Otsu阈值t，则停止生长，并舍去；重新从种子点集中选取新的种子点继续进行分割操作。

据此，整个计算过程之后便得到了直接减影图像中的目标血管区域图像。

本发明所述的基于SIFT特征点聚类及布尔差运算的DSA血管图像分割方法，在普通个人电脑上即可实施，无需特殊设备的支持，这使得本发明很容易推广。以上所述实施例仅表达了本发明的实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (1)

1.基于SIFT特征点聚类及布尔差运算的DSA血管图像分割方法，其特征包括以下步骤：

1）将对应的蒙片图像和活片图像进行直接减影，获得一张直接减影图像；在相同阈值条件下，利用SIFT算法分别在相对应的蒙片图像与活片图像中确定几何特征点，同时构造两个几何特征点集，分别为“蒙片图像几何特征点集”和“活片图像几何特征点集”，并将“活片图像几何特征点集”复制映射到直接减影图像上，生成“直接减影图像几何特征点集”；

2）在直接减影图像上根据几何特征点附近的灰度梯度进行几何特征点位置的局部调整，将几何特征点移动至直接减影图像中血管的内部，步骤如下：

①判断“直接减影图像几何特征点集”是否为空；若非空，则从该特征点集中取出一个几何特征点，转向②，否则几何特征点位置调整的计算过程结束；

②比较该几何特征点的灰度值与其8邻域内各像素的灰度值的大小；若该几何特征点的灰度值最小，则将该几何特征点从特征点集中移除，并转向①步；否则转向③步；

③计算该几何特征点与Kirsch算子8个模板卷积得到的8个方向上的梯度值，记为

使该几何特征点向梯度值最大的方向

移动一个像素单位，将移动后的几何特征点重新加入特征点集中，转向①步；其中，i为Kirsch算子8个模板中的第i个模板，1≤i≤8；

3）将“蒙片图像几何特征点集”和“直接减影图像几何特征点集”中的几何特征点根据欧氏距离进行聚类，然后进行几何特征点的类间布尔差运算，以防止处在毛细血管等精细结构上的几何特征点被误减，步骤如下：

①针对蒙片图像和直接减影图像分别构造两个几何特征点聚类集合，即“蒙片图像几何特征点聚类集合”和“直接减影图像几何特征点聚类集合”，并初始化这两个聚类集合为空；

②针对“蒙片图像几何特征点集”和“直接减影图像几何特征点集”，分别判断几何特征点集中是否有未被聚类的特征点；若有，则任选一个几何特征点，转向③步；若所有几何特征点都已被聚类处理，则转向④步；

③在蒙片图像/直接减影图像中，依次遍历“蒙片图像几何特征点聚类集合”/“直接减影图像几何特征点聚类集合”中的所有聚类，计算该几何特征点与某个聚类中的某个几何特征点之间的欧氏空间距离dis(p₀,p)＝sqrt((x_p0-x_p)²+(y_p0-y_p)²)，其中p₀,p分别为待聚类几何特征点和当前聚类中的某个几何特征点，x,y分别表示几何特征点的横纵坐标；当这两点间距离小于阈值dis_threshold时，将该几何特征点p₀归入点p所在聚类；若几何特征点p₀与所有聚类中的所有几何特征点都不满足欧氏空间距离小于阈值dis_threshold，则将该几何特征点作为新的聚类加入到“蒙片图像几何特征点聚类集合”/“直接减影图像几何特征点聚类集合”中，转向②步；

④针对“直接减影图像几何特征点聚类集合”中的每一个聚类

计算“蒙片图像几何特征点聚类集合”中所有聚类与该聚类

重叠的面积，然后再将所有重叠部分面积进行并集运算

（假设直接减影图像中共有W个聚类，蒙片图像中共有T个聚类），计算得到的总面积

与该聚类

的面积的比值，若该比值大于某一阈值

则l删d除该聚类，即将内所有几何特征点从直接减影图像上删去；

4）以类间布尔差运算后的几何特征点作为种子点集合，通过Otsu算法确定分割阈值，以邻域区域生长的方式进行血管图像的分割，步骤如下：

①初始化种子点集，将所有种子点标记为未生长状态，并计算当前直接减影图像的Otsu分割阈值；

②判断种子点集中的种子点是否全部为已生长状态；若不是，任取一个种子点，标记为已生长状态，转向③步，否则算法结束；

③考察该种子点8邻域内各像素的灰度值，取灰度值最小的前三个像素，判断它们与Otsu阈值的关系；若小于Otsu阈值，则将这些像素加入到种子点集中，并标记为未生长状态；若大于Otsu阈值，则舍去；转向②步。

Images