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

Abstract

本发明公开了一种在数字减影血管造影图像中分割血管数据的方法，其步骤为：①将DSA图像划分成N个互不重叠的a1×a2大小的子块；d≤a1、a2≤2d，b_1l、b_2l、b_3l、b_4l为大小为b1×b2且包含子块E_l在内的四个区域，其中，max(a1，a2)≤b1、b2≤3max(a1，a2)；②对各区域进行分割，然后利用繁忙度作为测度来判断对区域分割是否合理；③如果四个区域均被判定为背景，则将这四个区域的重叠子块E_l判为背景；如果四个区域中只有一个区域b_λl被判定包含血管，其中1≤λ≤4，则将该区域的分割门限T_λl作为子块E_l的分割门限；如果四个区域中不只一个区域被判定包含血管，则对于针对不同区域b_kl已求出的门限T_kl，选取其中最优的一个门限作为E_l的分割门限T_l。本发明方法能对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为血管直径的最大值；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)，令l＝1；

(2)对区域b_1l、b_2l、b_3l、b_4l进行分割，然后利用繁忙度作为测度来判断对区域分割是否合理，其过程为：

(2.1)对上述各区域进行分割，得到各区域的分割门限T_kl；

假设用b_kl(i，j)表示区域b_kl中象素的灰度值，在对区域b_kl分割后所得的二值图对应的象素强度由f_kl(i，j)来表示：

$f_{\mathrm{kl}}(i,j)=\left\{\begin{array}{ccc}0& \mathrm{if}& b_{\mathrm{kl}}(i,j)\&GreaterEqual;T_{\mathrm{kl}}\\ 1& \mathrm{if}& b_{\mathrm{kl}}(i,j)<T_{\mathrm{kl}}\end{array}\right.---\left(I\right)$

其中，1≤i≤b1，1≤j≤b2

(2.2)采用繁忙度作为判断对区域分割是否合理的测度，来接受或者拒绝对各区域的假设门限，其过程为：

(a)利用公式(II)计算上述各区域的繁忙度busy(b_kl)；

设区域b_kl有r行c列，r＝b1，c＝b2，繁忙度计算公式如下：

$\mathrm{busy}\left(b_{\mathrm{kl}}\right)=\underset{i=1,j=1}{\overset{i=r,j=c}{\mathrm{\&Sigma;}}}{t}_{\mathrm{kl}}(i,j)---\left(\mathrm{II}\right)$

其中，1≤i≤r，1≤j≤c，

(b)利用公式(III)计算作为血管存在性标准的繁忙度门限H_σ；

H_σ＝w·busy_max+(1-w)·busy_min    (III)

其中，busy_max、busy_min表示理想情况下b1×b2大小的区域繁忙度的最大或最小值，w是一个取值范围为[0，1]的系数；

(c)如果busy(b_kl)＞H_σ，则判定该区域b_kl只包含背景，将二值区域各象素值f_kl(i，j)置为背景所对应的逻辑值；否则，即判定该区域包含血管，保留二值区域的结果，即f_kl(i，j)值保持不变；

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

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

上述步骤(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{\&Sigma;}{[E_{\mathrm{kl}}(i,j)-\underset{i=1,j=1}{\overset{i=a1,j=a2}{\mathrm{\&Sigma;}}}{E}_{\mathrm{kl}}(i,j)/(a1\&times;a2)]}^2/(a1\&times;a2)---\left(\mathrm{IV}\right)$

2≤k≤m；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示出了只含背景(血管)的子块以及对其分别取门限后所得的二值子块，其中，图3a和图3b分别代表只含背景的子块和包含血管的子块，图3a′和图3b′分别表示对图3a和图3b取门限后所得的二值子块；

图4示出了对应于繁忙度最大(最小)值的二值图像，其中，图4a表示对应于繁忙度最大值的二值图像，图4b表示对应于繁忙度最小值的二值图像；

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

具体实施方式

下面结合附图来说明本发明的典型实施例。

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

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

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

(2)对区域b_1l、b_2l、b_3l、b_4l进行分割，然后利用繁忙度作为测度来判断对该区域分割是否合理，其过程为：

(2.1)对上述各区域进行分割，得到各区域的分割门限T_kl；本发明采用的是Otsu法，即最大类间方差法，当然，也可采用其他阈值分割方法如矩不变法、直方图凹度分析法、最大熵法等。

假设用b_kl(i，j)表示区域b_kl中象素的灰度值，在对区域b_kl分割后所得的二值图对应的象素强度由f_kl(i，j)来表示：

$f_{\mathrm{kl}}(i,j)=\left\{\begin{array}{ccc}0& \mathrm{if}& b_{\mathrm{kl}}(i,j)\&le;T_{\mathrm{kl}}\\ 1& \mathrm{if}& b_{\mathrm{kl}}(i,j)>T_{\mathrm{kl}}\end{array}\right.---\left(1\right)$

其中，1≤i≤b1，1≤j≤b2。

(2.2)采用繁忙度作为判断对区域分割是否合理的测度，来接受或者拒绝对各区域的假设门限。

(2.2.1)繁忙度作为测度的数学原理；

图3a和图3b分别代表只含背景的区域和包含血管的区域，图3a′和图3b′分别表示对图3a和图3b取门限后所得的二值区域。对于一个区域而言，在对其进行分割以后，所得的二值区域中，如果其象素值变迁频繁(0→1，1→0的转换频繁)，即该区域背景和目标分布较均匀，则该区域本来只包含背景；反之，因为血管本身的空间连续性，在对包含血管的区域进行二值分割后，所得区域的二值分布不会很均匀，会出现局部全0，局部全1的情况，即区域象素值变迁并不频繁。显然，图3a′中区域的象素值变迁很频繁，而3b′中区域的象素值变迁的频率就小得多。数学原理说明如下：

假设原始的二值图像区域为I_ori，I_x表示象素x的强度值，即I_x为0或1(逻辑值)，象素y为象素x的邻域象素并且I_y表示象素y的强度值。P_t(I_x＝α，I_y＝β)表示象素x的强度值为I_x同时y的强度值为I_y的概率。而一个同时包含血管和背景的区域以及一个仅含背景的区域的繁忙度分别被表示为busy′(b_kl)和busy″(b_kl)。在仅含背景的区域中象素x和y的强度值彼此独立，而在同时包含血管和背景的区域中象素y的强度值相关于象素x的值。这样，得到如下等式：

p_t(I_x＝α，I_y＝β)＝p_t(I_x＝α)·P_t(I_y＝β)  (x∈I_ori，y∈N(x)，α，β∈{0，1})(2)

${\mathrm{busy}}^{\&prime;}\left(b_{\mathrm{kl}}\right)=\underset{x\&Element;I_{\mathrm{ori}},\mathrm{\&alpha;}\&NotEqual;\mathrm{\&beta;}}{\mathrm{\&Sigma;}}{p}_t(I_x=\mathrm{\&alpha;})\&CenterDot;p_t(I_y=\mathrm{\&beta;})---\left(3\right)$

p_t(I_x＝α，I_y＝β)＝p_t(I_x＝α)·p_t(I_y＝β|I_x＝α)(x∈I_orl，y∈N(x)，α，β∈{0，1})(4)

${\mathrm{busy}}^{\&prime;\&prime;}\left(b_{\mathrm{kl}}\right)=\underset{x\&Element;I_{\mathrm{ori}},\mathrm{\&alpha;}\&NotEqual;\mathrm{\&beta;}}{\mathrm{\&Sigma;}}{p}_t(I_x=\mathrm{\&alpha;})\&CenterDot;p_t(I_y=\mathrm{\&beta;}|I_x=\mathrm{\&alpha;})---\left(5\right)$

很明显p_t(I_y＝β)＞p_t(I_y＝β|I_x＝α)。即使仅含背景的区域中象素x和y的强度值并不完全独立，同时包含血管和背景的区域中象素x和y的强度值的相关性也高于仅含背景的区域中，从而得到busy′(b_kl)＞busy″(b_kl)。

(2.2.2)利用公式(6)计算上述各区域的繁忙度busy(b_kl)；

设区域b_kl有r行c列(在此，r＝b1，c＝b2)，定义繁忙度计算公式如下：

$\mathrm{busy}\left(b_{\mathrm{kl}}\right)=\underset{i=1,j=1}{\overset{i=r,j=c}{\mathrm{\&Sigma;}}}{t}_{\mathrm{kl}}(i,j)---\left(6\right)$

其中，1≤i≤r，1≤j≤c。

(2.2.3)计算作为血管存在性标准的繁忙度门限；

取一个合适的门限H_σ如下：

H_σ＝w·busy_max+(1-w)·busy_min    (7)

在此，busy_max(busy_min)表示理想情况下b1×b2大小的区域繁忙度的最大(最小)值，w是一个取值范围为[0，1]的系数。

图4示出了理想情况下0对应于繁忙度最大(最小)值的二值图像。图4a表示对应于繁忙度最大值的二值图像；图4b表示对应于繁忙度最小值的二值图像。这里，采用3×3的区域作为样图，其中黑色和白色的单元格分别表示一个血管象素和一个背景象素，而这一状况也同样适用于b1×b2的区域。假设busy_max(busy_min)表示b1×b2大小的区域繁忙度的最大(最小)值。理想情况下有busy_max＝b1×b2且busy_min＝2min(b1，b2)，本发明中k＝0.2是一个令人满意的选择。

(2.2.4)判断各区域是否本来包含血管，作相应处理。

如果busy(b_kl)＞H_σ，则判定该区域b_kl只包含背景，将二值区域各象素值f_kl(i，j)置为背景所对应的逻辑值；否则，即判定该区域包含血管，保留二值区域的结果，即f_kl(i，j)值保持不变。

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

如果四个区域b_kl(k＝1，2，3，4)均被判定为背景，则将这四个区域的重叠子块E_l也判为背景；如果四个区域中只有一个区域b_λl(1≤λ≤4)被判定包含血管，则将该区域的分割门限T_λl作为子块E_l的分割门限；如果四个区域中不只一个区域被判定包含血管，则对于针对不同区域b_kl已求出的门限T_kl，选取其中最优的一个门限作为E_l的分割门限T_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{\&Sigma;}{[E_{\mathrm{kl}}(i,j)-\underset{i=1,j=1}{\overset{i=a1,j=a2}{\mathrm{\&Sigma;}}}{E}_{\mathrm{kl}}(i,j)/(a1\&times;a2)]}^2/(a1\&times;a2)---\left(8\right)$

2≤k≤m；1≤i≤a1，1≤j≤a2。

V_kl越大表示区域内灰度越不均匀，即包含的目标血管越多，背景越少。本发明选取V_kl中的最大值所对应的门限T_kl作为子块E_l的分割门限。对所有子块E_l重复此过程，即完成了对整幅DSA图像的分割。

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

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

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

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

Claims (2)

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

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

(2)对区域b_1l、b_2l、b_3l、b_4l进行分割，然后利用繁忙度作为测度来判断对区域分割是否合理，其过程为：

(2.1)对上述各区域进行分割，得到各区域的分割门限T_kl；

假设用b_kl(i，j)表示区域b_kl中象素的灰度值，在对区域b_kl分割后所得的二值图对应的象素强度由f_kl(i，j)来表示：

$f_{\mathrm{kl}}=\left\{\begin{array}{ccc}0& \mathrm{if}& b_{\mathrm{kl}}(i,j)\&GreaterEqual;T_{\mathrm{kl}}\\ 1& \mathrm{if}& b_{\mathrm{kl}}(i,j)<T_{\mathrm{kl}}\end{array}\right.---\left(I\right)$

其叶中，1≤i≤b1，1≤j≤b2

(2.2)采用繁忙度作为判断对区域分割是否合理的测度，来接受或者拒绝对各区域的假设门限，其过程为：

(a)利用公式(II)计算上述各区域的繁忙度busy(b_kl)；

设区域b_kl有r行c列，r＝b1，c＝b2，繁忙度计算公式如下：

$\mathrm{busy}\left(b_{\mathrm{kl}}\right)=\underset{i=1,j=1}{\overset{i=r,j=c}{\mathrm{\&Sigma;}}}{t}_{\mathrm{kl}}(i,j)---\left(\mathrm{II}\right)$

其中，1≤i≤r，1≤j≤c，

(b)利用公式(III)计算作为血管存在性标准的繁忙度门限H_σ；

H_σ＝w·busy_max+(1-w)·busy_min    (III)

其中，busy_max、busy_min表示理想情况下b1×b2大小的区域繁忙度的最大或最小值，w是一个取值范围为[0，1]的系数；

(c)如果busy(b_kl)＞H_σ，则判定该区域b_kl只包含背景，将二值区域各象素值f_kl(i，j)置为背景所对应的逻辑值；否则，即判定该区域包含血管，保留二值区域的结果，即f_kl(i，j)值保持不变；

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

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

2.根据权利要求1所述的方法，其特征在于：步骤(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{\&Sigma;}{[E_{\mathrm{kl}}(i,j)-\underset{i=1,j=1}{\overset{i=a1,j=a2}{\mathrm{\&Sigma;}}}{E}_{\mathrm{kl}}(i,j)/(a1\&times;a2)]}^2/(a1\&times;a2)$

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

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

Images