CN1897034A - 一种在数字减影血管造影图像中分割血管数据的方法 - Google Patents

Abstract

本发明公开了一种在数字减影血管造影图像中分割血管数据的方法，其步骤为：①将DSA图像划分成N个互不重叠的a1×a2大小的子块E₁、E₂、……、E_N，其中，d≤a1≤2d，d≤a2≤2d，d为血管直径的最大值，令l＝1；②判断各区域b_1l、b_2l、b_3l和b_4l中是否包含血管，b_1l、b_2l、b_3l、b_4l为大小为b1×b2且包含子块El在内的四个区域，其中，max(a1，a2)≤b1≤3max(a1，a2)，max(a1，a2)≤b2≤3max(a1，a2)；③对各区域的重叠子块进行分割；④令1＝1+1，对子块E_l重复步骤②和③，完成整幅DSA图像的分割。本发明方法能对DSA图像有效分割，并且操作简便，工作效率高。

Description

一种在数字减影血管造影图像中分割血管数据的方法

技术领域

本发明属于医学成像技术领域，具体涉及一种在数字减影血管造影图像中分割血管数据的方法。

背景技术

数字减影血管造影(DSA：Digital Subtraction Angiography)技术在临床已应用20多年，是心脑血管疾病无创诊断与介入治疗手术导航的重要依据。DSA图像处理中的一个关键任务就是进行图像分割，以使血管的生理特征能够更清楚的显示出来，结构分析、运动分析、三维可视化等后续操作，以及图像引导手术、肿瘤放射治疗、治疗评估等应用研究都是以图像分割为基础的。由于X射线经过的组织厚度以及血液中的造影剂浓度不均匀，经过剪影后的图像不能完全消除人体组织引起的噪声信号，背景和目标混杂在一起。另外，人体的组织结构和形状很复杂，而且人与人之间有相当大的差异。因此血管的分割是一项困难的任务。

C.Kirbas和F.Quek提出，现有的的血管分割技术大致可以分为6大类：模式识别方法、基于模型的方法、基于跟踪的方法、人工智能方法、神经网络方法、管状物体检测法，其中一些类型还可以细分(见“血管提取技术和算法的综述”，ACM Computing Surveys，36(2)：81-121，2004)。而阈值分割是最常见的直接检测区域的分割方法，其优点是简单，同时对于不同类的物体灰度值或其他特征值相差很大时，它能很有效的对图像进行分割，但它对图像中不存在明显灰度差异或灰度值范围有较大重叠的图像分割问题难以得到准确的结果。另外，由于它仅仅考虑图像的灰度信息而不考虑图像的空间信息，所以对噪声和灰度不均匀很敏感。

阈值分割方法分为全局阈值分割和局部阈值分割技术。早期用得较普遍的是全值阈值分割方法，但其并不适用于某些情况。以亮目标和暗背景的图像为例，目标的某些区域可能比背景的部分区域要暗，则任意一个确定的门限无法完全将目标和背景分离开来。也就是说，现有的全局分割技术很难在血管和背景灰度范围接近的局部地方取得好的分割效果，不适用于灰度直方图呈细窄的单峰分布的DSA图像，很难分得连续的血管段。

而另一方面，现有的局部阈值分割技术通常对图像进行分块后，对每个子块分别计算分割门限。由于未判断子块中是否存在血管，因此对不含血管的背景区域也将得到一个分割结果，这显然是不合理的。虽然直方图双峰法计算分割门限时，可利用直方图的双峰性在一定程度上对子块中是否存在血管作出判决，但其要求子块内所含血管点数和背景点数相当，这是十分苛刻很难满足的，且当子块较小，包含的象素数较少时，直方图双峰性的判断本身也存在困难。

发明内容

本发明的目的在于提供一种在数字减影血管造影图像中分割血管数据的方法，该方法能对DSA图像有效分割，并且操作简便，工作效率高。

本发明提供的一种在数字减影血管造影图像中分割血管数据的方法，其步骤为：

(1)将DSA图像划分成N个互不重叠的a1×a2大小的子块E₁、E₂、……、E_N，其中，d≤a1≤2d，d≤a2≤2d，d为血管直径的最大值，令l＝1；

(2)判断各区域b_1l、b_2l、b_3l、b_4l中是否包含血管，b_1l、b_2l、b_3l、b_4l为大小为b1×b2且包含子块E_l在内的四个区域，其中，max(a1，a2)≤b1≤3max(a1，a2)，max(a1，a2)≤b2≤3max(a1，a2)；判断过程为：

(2.1)计算各区域b_kl的方差值σ_kl，k＝1、2、3和4；

(2.2)将方差值σ_kl与血管存在性标准的门限σ_Tl进行比较，如果区域b_kl的方差σ_kl小于σ_Tl，则记该区域为背景；否则，计算该区域的分割门限，记为T_kl；

(3)对各区域的重叠子块进行分割：

如果四个区域b_kl均被判定为背景，则将这四个区域的重叠子块E_l判为背景；如果四个区域中只有一个区域b_λl被判定包含血管，1≤λ≤4，则将该区域的分割门限T_λl作为子块E_l的分割门限T_kl；如果四个区域中不只一个区域被判定包含血管，则对于针对不同区域b_kl已求出的分割门限T_kl，选取其中一个最优的门限作为E_l的分割门限T_l；

(4)令l＝l+1，对子块E_l重复步骤(2)和(3)，完成整幅DSA图像的分割。

上述步骤(2.2)中门限σ_Tl采用以下方法确定：

预先训练DSA图像中大量的b1×b2大小的区域，得到对应于方差的不同取值范围的概率曲线p(x)，其中x轴表示b1×b2大小的区域的方差的不同取值范围而p(x)表示方差在不同取值范围的概率，假设p₁(x)对应于仅含背景的b1×b2大小的区域而p₂(x)对应于同时包含血管和背景的b1×b2大小的区域，选择p₁(x)和p₂(x)曲线的交点处的x_T值作为血管存在性门限的范围，范围x_T的中值σ_Tl被用来作为区域E_l的血管存在性标准的门限。

上述步骤(3)采用以下方法确定子块E_l的分割门限T_l：

设共有m个区域被判定包含血管，2≤m≤4，对子块E_l，设分别对其取门限T_kl后所得的二值子块为E_kl，假设二值子块E_kl的象素值为E_kl(i，j)，计算E_kl的方差作为选取分割门限T_l的测度：

$V_{\mathrm{kl}}=\mathrm{\Σ}{[E_{\mathrm{kl}}(i,j)-\underset{i=1,j=1}{\overset{i=a1,j=a2}{\mathrm{\Σ}}}{E}_{\mathrm{kl}}(i,j)/(a1\×a2)]}^2/(a1\×a2)---\left(I\right)$

2≤k≤n；1≤i≤a1，1≤j≤a2；

选取V_kl中的最大值所对应的门限T_kl作为子块E_l的分割门限T_l。

本发明针对现有全局阈值分割技术以及部分局部阈值分割技术对DSA图像处理的缺陷，根据DSA图像的特点，提出一种基于血管直径先验知识以及重叠分块技术的局部阈值分割方法。本发明根据血管直径的先验知识将图像分成若干区域，所取区域大小保证在区域内绝对含有背景，因此区域只存在两种可能性：包含血管和只含背景。本发明采取一些判断准则预先判断区域内是否包含血管，然后对包含血管的区域采用一种阈值分割方法求取其分割门限，对只包含背景的区域直接将其标记为背景。同时本发明还采用重叠分块技术，综合考虑重叠区域象素信息来为区域的重叠子块选取最优门限，这样可以避免在血管所占比例较少的子块内局部阈值分割技术形同全局阈值分割技术，改善局部阈值分割的效果。

附图说明

图1示出了本方法的流程图；

图2示出了子块E_l与子块b_kl(k＝1，2，3，4)的位置示意图；

图3示出了对应于b1×b2大小的区域的方差的不同取值范围的概率曲线p(x)；

图4所示的反映了在本发明的典型实施例中使用上述方法与其他全局和局部阈值分割方法的一个比较结果，其中，图4a表示原始肝脏DSA图像；图4b表示本方法的分割结果；图4c表示基于非重叠分块的局部阈值分割结果；图4d表示用Bernsen的方法的分割结果；图4e用Niblack的方法的分割结果；图4f表示用全局阈值分割方法的结果。

具体实施方式

将参考附图来说明本发明的典型实施例。

下面的说明将假定图像含有暗目标和亮背景，以门限T进行分割可以将灰度值小于或等于T值的象素以及大于T值的象素分别划分为目标和背景，得到一个二值图像。相应的，在二值图像中，象素值为0对应目标，象素值为1对应背景。这里的1与0都是逻辑值。当然，本发明同样适用于亮目标和暗背景的情况。

下面结合图1详细说明本发明的步骤：

(1)将DSA图像划分成N个互不重叠的a1×a2大小的子块E_l，l＝1，2…N，其中，d≤a1≤2d，d≤a2≤2d，在此d被假设为血管直径的最大值，大小为b1×b2且包含子块E_l在内的区域为b_1l、b_2l、b_3l、b_4l，其中，max(a1，a2)≤b1≤3max(a1，a2)，max(a1，a2)≤b2≤3max(a1，a2)；

(2)判断各区域b_1l、b_2l、b_3l、b_4l中是否包含血管，其过程为：

(2.1)计算各区域b_1l、b_2l、b_3l、b_4l的方差值σ_kl；

(2.2)将方差值σ_kl与血管存在性标准的门限σ_Tl进行比较，如果区域b_kl的方差σ_kl小于σ_Tl，说明该区域内很均匀，则记该区域为背景；否则，采用一种分割方法(本发明采用的是Otsu法，即最大类间方差法，当然，也可采用其他阈值分割方法如矩不变法、直方图凹度分析法、最大熵法等)计算该区域的分割门限，记为T_kl。其中，门限σ_Tl通过DSA图像的先验知识得到，具体方法如下：

预先训练DSA图像中大量的b1×b2大小的区域，得到如图3所示的对应于方差的不同取值范围的概率曲线p(x)，其中x轴表示b1×b2大小的区域的方差的不同取值范围而p(x)表示方差在不同取值范围的概率。假设p₁(x)对应于仅含背景的b1×b2大小的区域而p₂(x)对应于同时包含血管和背景的b1×b2大小的区域。选择p₁(x)和p₂(x)曲线的交点处的x_T值作为血管存在性门限的范围，以使判定一个给定区域是仅包含背景还是同时包含血管和背景时的平均出错率降至最小。范围x_T的中值σ_Tl被用来作为区域E_l的血管存在性标准的门限。

(3)综合各重叠区域的分割情况对各区域的重叠子块进行分割，具体方法如下：

如果四个区域b_kl(k＝1，2，3，4)均被判定为背景，则将这四个区域的重叠子块E_l也判为背景；如果四个区域中只有一个区域b_λ1(1≤λ≤4)被判定包含血管，则将该区域的分割门限T_λ1作为子块E_l的分割门限；如果四个区域中不只一个区域被判定包含血管，则对于针对不同区域b_kl已求出的门限T_kl，选取其中最优的一个门限作为E_l的分割门限T_l。

下面介绍对于多个备选门限选取分割门限T_l的方法：

设共有m个区域被判定包含血管(2≤m≤4)，采用可以表示区域内部均匀性的方差作为选取分割门限T_l的测度：对子块E_l，设分别对其取门限T_kl后所得的二值子块为E_kl，假设二值子块E_kl的象素值为E_kl(i，j)，计算E_kl的方差作为选取分割门限T_l的测度：

$V_{\mathrm{kl}}=\mathrm{\Σ}{[E_{\mathrm{kl}}(i,j)-\underset{i=1,j=1}{\overset{i=a1,j=a2}{\mathrm{\Σ}}}{E}_{\mathrm{kl}}(i,j)/(a1\×a2)]}^2/(a1\×a2)---\left(1\right)$

在此，2≤k≤n；1≤i≤a1，1≤j≤a2。

V_kl越大表示区域内灰度越不均匀，即包含的目标血管越多，背景越少。本发明选取V_kl中的最大值所对应的门限T_kl作为子块E_l的分割门限。

对所有子块E_l重复(2)(3)过程，即完成了对整幅DSA图像的分割。

对于亮目标和暗背景的情况所采取的步骤同上，只不过对子块求得分割门限T时，将灰度值小于或等于T值的象素以及大于T值的象素分别划分为背景和目标，得到的二值图像中，象素值为1对应目标，象素值为0对应背景，这里的1与0都是逻辑值。

图4所示的反映了在本发明的典型实施例中使用上述方法以及与其他全局和局部阈值分割方法的一个比较结果。从中可以看出，对整幅图像选一个全局门限很难在血管和背景的灰度范围彼此接近的局部区域将血管段从背景中分割出来，部分微细血管难以被有效保留，甚至出现不连续的血管段。而本发明因为采取了一些判断准则来去除背景点以及基于血管直径的先验知识的重叠分块技术，在有效去除背景点的同时尽可能保留了有用信息，和其他局部阈值分割技术(比如Bernsen的方法、Niblack的方法、基于非重叠分块的局部阈值分割方法等)相比，可以达到更好的二值分割结果。同时本发明具有执行效率高的特点。

根据本发明的典型实施例，用于实现本发明的计算机系统可以包括，特别是，中央处理器(CPU)、存储器和输入/输出(I/O)接口。计算机系统通常通过I/O接口与显示器和诸如鼠标和键盘此类的各种输入设备相连，配套电路可以包括像高速缓存、电源、时钟电路和通信总线这样的电路。存储器可以包括随机存储器(RAM)、只读存储器(ROM)、磁盘驱动器、磁带机等，或它们的组合。计算机平台还包括操作系统和微指令代码。此处所述各种过程和功能可以是通过操作系统执行的微指令代码或应用程序(或它们的组合)的一部分。此外，各种其他外围设备可以连接到该计算机平台，如附加数据存储设备和打印设备。

还应理解，因为附图中所述的某些构成系统的组件和方法步骤可以软件形式来实现，所以系统组件(或过程步骤)之间的实际连接可能有所不同，具体视本发明的编程方式而定。基于此处提出的本发明原理，相关领域的普通专业人员可以设想本发明的这些以及类似实施方案或配置。

Claims (3)

1、一种在数字减影血管造影图像中分割血管数据的方法，其步骤为：

(1)将DSA图像划分成N个互不重叠的a1×a2大小的子块E₁、E₂、......、E_N，其中，d≤a1≤2d，d≤a2≤2d，d为血管直径的最大值，令l＝1；

(2)判断各区域b_1l、b_2l、b_3l、b_4l中是否包含血管，b_1l、b_2l、b_3l、b_4l为大小为b1×b2且包含子块E_l在内的四个区域，其中，max(a1，a2)≤b1≤3max(a1，a2)，max(a1，a2)≤b2≤3max(a1，a2)；判断过程为：

(2.1)计算各区域b_kl的方差值σ_kl，k＝1、2、3和4；

(2.2)将方差值σ_kl与血管存在性标准的门限σ_Tl进行比较，如果区域b_kl的方差σ_kl小于σ_Tl，则记该区域为背景；否则，计算该区域的分割门限，记为T_kl；

(3)对各区域的重叠子块进行分割：

如果四个区域b_kl均被判定为背景，则将这四个区域的重叠子块E_l判为背景；如果四个区域中只有一个区域b_λl被判定包含血管，1≤λ≤4，则将该区域的分割门限T_λl作为子块E_l的分割门限T_kl；如果四个区域中不只一个区域被判定包含血管，则对于针对不同区域b_kl已求出的分割门限T_kl，选取其中一个最优的门限作为E_l的分割门限T_l；

(4)令l＝l+1，对子块E_l重复步骤(2)和(3)，完成整幅DSA图像的分割。

2、根据权利要求1所述的方法，其特征在于：门限σ_Tl采用以下方法确定：

预先训练DSA图像中大量的b1×b2大小的区域，得到对应于方差的不同取值范围的概率曲线p(x)，其中x轴表示b1×b2大小的区域的方差的不同取值范围而p(x)表示方差在不同取值范围的概率，假设p₁(x)对应于仅含背景的b1×b2大小的区域而p₂(x)对应于同时包含血管和背景的b1×b2大小的区域，选择p₁(x)和p₂(x)曲线的交点处的x_T值作为血管存在性门限的范围，范围x_T的中值σ_Tl被用来作为区域E_l的血管存在性标准的门限。

3、根据权利要求1或2所述的方法，其特征在于：步骤(3)采用以下方法确定子块E_l的分割门限T_l：

设共有m个区域被判定包含血管，2≤m≤4，对子块E_l，设分别对其取门限T_kl后所得的二值子块为E_kl，假设二值子块E_kl的象素值为E_kl(i，j)，计算E_kl的方差作为选取子块E_l的分割门限T_l的测度：

$V_{\mathrm{kl}}=\mathrm{\Σ}{[E_{\mathrm{kl}}(i,j)-\underset{i=1,j=1}{\overset{i=a1,j=a2}{\mathrm{\Σ}}}{E}_{\mathrm{kl}}(i,j)/(a1\×a2)]}^2/(a1\×a2)$

2≤k≤n；1≤i≤a1，1≤j≤a2；

选取V_kl中的最大值所对应的门限T_kl作为子块E_l的分割门限T_l。

Images