CN104504711A - 一种基于圆形轮廓图极性的血管图像处理方法 - Google Patents

Abstract

本发明公开一种基于圆形轮廓图极性的血管图像处理方法，属于图像处理技术领域，提出了从几何结构分布来判断图像区域特点的两个圆形轮廓极性的指标，分别为圆形轮廓图极性指标一和指标二，能够鉴别血管图像中的噪声污染区域与血管区域，从而分别对噪声污染区域和血管区域采取不同的阈值化方法进行处理，对污染区域采取亮度值阈值化与结构阈值化两种处理，而对血管区域则仅采用结构阈值化处理，可以有效得去除血管图像中的噪声，得到高质量的血管提取结果。

Description

一种基于圆形轮廓图极性的血管图像处理方法

技术领域

本发明属于图像处理技术领域，具体涉及一种基于圆形轮廓图极性的血管图像处理方法。

背景技术

血管图像在医学研究和临床诊断上有着广泛应用，比如肿瘤血管生成的定量分析，糖尿病视网膜病变的诊断等。各种技术应用于生成血管图像，许多成像方法也在发展中。血管成像技术包括计算机断层扫描、磁共振血管造影技术、整装电镜和双光子荧光血管造影等，新型成像技术能够获取微米级的分辨率，能够更好地获取小血管(毛细血管)图像。这些优秀技术的发展使得我们能够一窥血管结构。由于血管在拓扑学和几何学上的复杂度，对血管进行形态学分析仍然具有很大挑战。这些分析的可靠性严重依赖于诸如去噪、阈值化等的图像预处理方法的结果。因此，如何提出和建立更高精度和效率的血管图像预处理算法具有十分重要的意义。

图像预处理的一个主要目标是去除噪声的同时更多地保留真实信号。对于血管图像，两种类型的噪声最为常见。在统计血管表征之前，需要特别处理这两种噪声。一种噪声是由于成像背景污染、血管漏洞或者荧光染料通过血管壁渗透扩散至周围组织的云状污染。另一种是染色剂或镜面反射造成的结构噪声，这种噪声比血管有更高的信号强度，即亮度值更大。当荧光染料渗透扩散时，污染和血管有着相近的信号强度，这使得一些小动脉、小静脉和毛细血管很容易被污染覆盖，令我们很难从污染区域中区分出血管。因此，在研究更有效的算法时，通常需要考虑那些更具血管鉴别力的特征(例如血管形状)。此外，在一些阈值化算法中，尤其是血管图像的阈值化算法，那些高强度的结构噪声很容易被误认为是真实信号，从而误识别为血管。因此，当前迫切需要那些能够恰当的处理图像污染和结构噪声的血管提取方法。

目前，血管图像的阈值化方法、分割和提取方法已经有大量的研究。根据所采用的图像特征，血管提取方法大致可分为三类：基于像素强度值的算法，基于形状的算法，考虑几何特征和像素强度值的综合性算法。其中，基于像素强度值的算法主要利用聚类、直方图、熵、像素强度值的局部统计和图像特征(例如边缘)来获取阈值。对于严重污染的图像，这些算法的结果很难让人满意。而且截至目前，没有一个算法能够去除结构噪声。基于形状的算法通常假定血管的形状与圆柱体、管状体、椭圆体相似，从而区分血管和背景、噪声。但是这些算法的参数化描述通常十分复杂(例如超椭圆体，三维B样条曲面，三维圆柱体模型)，而且复杂的血管分支模式，以及形状的近似误差，也对这些算法的效果有着较大的影响。

最后，尽管能够使用一些复杂的图像处理程序包，诸如Amira，VidaSuite，VesSeg，来代替上述算法，但是必须通过人工调整和修改大量的图像处理步骤(例如中值滤波、高斯滤波、图像增强、腐蚀、膨胀和间隙填充)。因为需要人工调整大量的参数，因此导致非常大的工作量。

发明内容

为了克服上述现有技术存在的缺陷，本发明的目的在于提供一种基于圆形轮廓图极性的血管图像处理方法，该方法同时考虑血管几何特征和像素强度，能够有效去除血管图像中的噪声，从而得到高质量的血管提取结果。

本发明是通过以下技术方案来实现：

一种基于圆形轮廓图极性的血管图像处理方法，包括：

首先，输入采集的血管图像，进行背景清理；

其次，对经过背景清理的血管图像通过窗口自适应处理，确定血管区域和噪声污染区域；

判断经窗口自适应处理后的血管图像的圆形轮廓图极性指标一是否小于设定的阈值t_h：若是，则判断为噪声污染区域，若否，则判断为血管区域；

最后，分别对血管区域和污染区域进行阈值化处理：

对噪声污染区域同时采取亮度值阈值化与结构阈值化处理；对血管区域仅采用结构阈值化处理；将噪声污染区域和血管区域的处理结果汇总，得到高质量的血管提取结果；

所述结构阈值化处理基于对血管图像圆形轮廓图极性指标二的处理；

所述的圆形轮廓图极性指标一为PCP₁：

其中，σ²表示该窗口内图像的方差，μ表示该窗口内图象的均值；

所述的圆形轮廓图极性指标二为PCP₂：

${\mathrm{PCP}}_2=E[\underset{m=\mathrm{1,2},...,M}{\max }f\left({\mathrm{\θ}}_m\right)-\underset{l=\mathrm{1,2},...,M}{\min }f\left({\mathrm{\θ}}_l\right)],$ E(·)表示数学期望函数，f(θ_m)和f(θ_l)分别表示圆周中θ_m和θ_l两个方向上的圆形轮廓函数的值。

具体包括以下步骤：

步骤一：图像背景处理

输入血管图像，清理血管图像的背景，统计血管图像中最亮的像素值，取该像素值的5％作为背景清理阈值，将图像中的像素值低于该背景清理阈值的像素置为零；

步骤二：窗口自适应

将整幅血管图像分割成多个子窗口，设置初始窗口尺寸值为w_ini，并定义阈值t_h，计算w_ini的圆形轮廓图极性指标一PCP₁；

如果计算得到的PCP₁的值小于阈值t_h，则判断该区域属于血管被噪声覆盖的噪声污染区域；

如果计算得到的PCP₁的值大于阈值t_h，则扩大窗口尺寸，再次计算PCP₁的值，并与上次计算所得的PCP₁的值进行比较，重复计算PCP₁的值与比较操作，直到PCP₁的值减小的幅度大于1为止，确定此时的窗口尺寸，这些区域属于血管区域；

步骤三：基于圆形轮廓图的血管阈值化处理

对判断属于噪声污染区域的图像区域，先进行亮度阈值化处理，再进行基于圆形轮廓图极性指标二PCP₂的几何结构阈值化处理，从而完成血管提取；

对判断属于血管区域的图像区域，进行基于PCP₂的几何结构阈值化处理；

将对所有子窗口进行的阈值化处理结果汇总整合，得到最终的血管提取结果。

步骤三所述的亮度阈值化处理采用三级或四级Otsu阈值化方法，并采用该阈值化方法中最高亮度作为亮度阈值，将大于等于该亮度阈值的像素保存，小于该亮度阈值的像素则置为0。

步骤三所述的基于圆形轮廓图极性指标二PCP₂的几何结构阈值化处理为：

当该图像区域中一个像素处所获得的PCP₂的值大于特定阈值t_p＝15时，判断该像素为血管像素并将其置1，反之则为噪声像素并将其置为0。

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

本发明公开的基于圆形轮廓图极性的血管图像处理方法，提出了从几何结构分布来判断图像区域特点的两个圆形轮廓极性的指标，分别为圆形轮廓图极性指标一和指标二，通过圆形轮廓图极性指标一能够区分血管图像中的污染区域和血管区域，而圆形轮廓图极性指标二则能够有效得区分血管像素，并去除结构噪声。本发明方法能够鉴别血管图像中的噪声污染区域与血管区域，从而分别对噪声污染区域和血管区域采取不同的阈值化方法进行处理，对污染区域采取亮度值阈值化与结构阈值化(基于圆形轮廓图极性指标二)两种处理，而对血管区域则仅采用结构阈值化处理，可以有效得去除血管图像中的噪声，得到高质量的血管提取结果。该方法无需复杂的数学模型，计算速度快，抗噪声能力强，血管图像提取精度高。

附图说明

图1a为包含云状污染的血管图像；

图1b为包含结构噪声的血管图像；

图2a为血管图像中的各种情形；

图2b-1为背景干净的血管像素的场景示意图；

图2b-2为背景干净的血管像素处的圆形轮廓图；

图2c-1为污染区域内像素处的场景示意图；

图2c-2为污染区域内像素处的圆形轮廓图；

图2d-1为污染区域边界像素处的场景示意图；

图2d-2为污染区域边界像素处的圆形轮廓图；

图2e-1为干净背景中结构噪声处像素的场景示意图；

图2e-2为干净背景中结构噪声处像素的圆形轮廓图；

图2f-1为污染区域中靠近血管像素的场景示意图；

图2f-2为污染区域中靠近血管像素的圆形轮廓图；

图2g-1为靠近血管的结构噪声处像素的场景示意图；

图2g-2为靠近血管的结构噪声处像素的圆形轮廓图；

图3为五种情景中的PCP₂范围图；

图4为发明的血管图像处理方法流程图；

图5a血管原图；

图5b本发明方法的血管提取结果图；

图5c Jiang算法的血管提取结果；

图5d Qian算法的血管提取结果；

图5e Otsu算法的血管提取结果；

图5f Huang算法的血管提取结果；

图5g Kittler算法的血管提取结果；

图6a本发明用于三维血管提取的原图；

图6b本发明用于三维血管提取的结果图。

具体实施方式

下面结合具体的实施例对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。

本发明分析了血管图像中的两类噪声的特点，为了在准确提取血管的同时去除两类噪声，需要设计一种新型的血管提取方法。该方法需要满足以下条件：1)能够同时处理云状污染和结构噪声；2)无需引入复杂和特殊形状约束的血管几何特征；3)具有高精度和高效率。

基于上述需求，本发明有如下两个考虑。首先，血管是管状物而污染和结构噪声不是；其次，“血管形状”能够使用数学或统计方法量化，而不是利用特殊的形状约束。一些先前的工作对“血管形状”已经使用相似的概念做了尝试。例如，海森矩阵方法考虑了局部区域的特征向量；不幸的是这种方法对噪声十分敏感。受圆形(循环)统计学的启发，本发明发现血管像素周围像素的空间分布和噪声像素周围像素的空间分布不同。基于此发现，本发明首先定义了一个名为“圆形轮廓图”的统计采样结构，然后通过探查像素圆形轮廓图(Polarityof Circular Profile，简写为PCP)的极性来判断该像素是真实血管信号还是噪声，本发明主要检测了每个像素周围像素的空间分布形式。通过大量的仿真实验和真实血管图像实验证明，对于血管提取，本发明中的方法在去除污染和结构噪声方面远优于其它现有方法。

首先，详细介绍圆形轮廓图的概念。血管图像中最为常见的两种噪声包括云状污染噪声(如图1a所示)和结构噪声(如图1b所示)，其中结构噪声指的是其亮度比血管像素还要亮的一些亮点，无法采用一般的阈值化处理方法去除。一个合理的血管检测算法应该能够在各种情况下检测出血管形状，无论血管拓扑结构(例如层级结构、分支形式)，几何结构(例如血管密度、形状、直径大小等)，噪声级别和分布有着多大的差异。为了实现这个功能，一个方案是在不同的方向取样多个点，然后总结这些采样点的特性。更具体地说，对于任意一个抽样点(称为探针)，将对所有位于以该抽样点为圆心的一个圆内的非背景像素的空间分布形式进行检测，检测结果称为圆形轮廓图。

令i(i＝1,2,...,N)表示第i个探针，θ(0≤θ＜2π)表示抽样方向和极轴之间的极角，j_k(k＝1,2,...,K)表示θ方向上的抽样点，K表示抽样点的数量，T_K(θ)是θ方向上抽样点集合{j₁,j₂,...,j_K}的统计量。考虑到可实现性，对每一个探针仅有有限的不同方向进行抽样(例如间隔特定角度的多个方向)，分别用θ_m(m＝1,2,...,M)表示。因此，探针i的圆形轮廓图可以用一个集合{T_K(θ_m):m＝1,2,...,M}_i定义。

为了更好地理解这个定义，本发明列出了一个有污染和结构噪声的图像(图2a)，同时标出了6个具有典型意义的血管图像中的场景(图2a中的点b、c、d、e、f和g)。此外，在所检测的像素上放置探针，在每个探针上选取16个方向(M＝16)，并在每个方向上取30个连续的像素点(K＝30)。简单地表示亮度非零像素的统计值，其中I_k表示像素j_k的强度，1(·)是一个指示函数。图2b到图2g分别给出了对应6个典型场景的抽样所得的圆形轮廓图，其中图2b-1到图2g-1是场景及抽样方法示意图，图2b-2到图2g-2是所得的圆形轮廓图。需要说明的是，为了更好的可视化，图2b到图2g已经被旋转，因此它们可能与图2a中的方向不同。

对于在干净背景中的血管上的探针(图2b)，T_K(θ)在血管轴向方向上达到最大值，而在垂直于血管轴向方向上取得最小值。因此，圆形轮廓图形状与一个圆有较大的偏差。对于在污染区域的探针，如果探针大部分在污染区域(图2c)或靠近血管(图2f)，T_K(θ)在各个方向上几乎有相同的值，这使得圆形轮廓图形状接近一个圆。当探针处于污染区域的边界时(图2d)，圆形轮廓图向污染区域隆起，不同于其他情形(图2b、图2c、图2f)的圆形轮廓图。如果探针位于结构噪声处(图2e)，或在结构噪声处且离血管较近(图2g)，圆形轮廓图的结构随机不规则，但是T_K(θ)的最大值和最小值的绝对差是一个很小的值。基于这些分析结果可见，圆形轮廓图的差别能够区分血管像素和噪声。

其次，详细介绍圆形轮廓图函数及其量化。如上所述，圆形轮廓极性指标(PCP)的定义应该能够反映一个圆形轮廓图是否向某个方向隆起。也就是说，一个大的PCP值表示T_K(θ)在不同θs方向上的差异较大，而一个小的PCP值代表T_K(θ)在不同方向上差异很小。因此，T_K的值能够量化PCP，可以通过考虑图2b到图2g的每个场景来处理这个问题。

令K表示某一方向上采样像素的数目，可以认为它表示了一个探针的半径，将从探针中心到半径尾部的线段命名为探针臂。假设探针位于像素i上，令I_i表示像素的强度，方向θ_m上的统计量T_K(θ)定义如下：

$T_K\left({\mathrm{\θ}}_m\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}1(I_k>0),m=\mathrm{1,2},...,M,$   公式(一)

其中，1(·)为指示函数，I_k表示像素j_k的强度，K表示一个方向上的采样点的数量，M表示采样的方向的数量，同时，令p_v和p_n分别表示血管和噪声的非背景像素的概率分布密度。圆形轮廓函数定义如下：

f(θ)＝E(T_K(θ))，

其中E(·)表示数学期望函数。对于图2b中的场景，图像中只有血管，假设探针位于血管，令d_v表示极坐标和上层血管壁的距离，那么圆形轮廓图函数可以定义为：

$f\left(\mathrm{\&theta;}\right)=\left\{\begin{array}{cc}{\mathrm{kp}}_v,& 0\&le;\mathrm{\&theta;}\&le;{\mathrm{\&theta;}}_t,\mathrm{\&pi;}-{\mathrm{\&theta;}}_t\&le;\mathrm{\&theta;}\&le;\mathrm{\&pi;}+{\mathrm{\&theta;}}_t,2\mathrm{\&pi;}-{\mathrm{\&theta;}}_t\&le;\mathrm{\&theta;}<2\mathrm{\&pi;\&theta;}\\ \frac{d_v{p}_v}{\sin \mathrm{\&theta;}},& {\mathrm{\&theta;}}_t\&le;\mathrm{\&theta;}\&le;\mathrm{\&pi;}-{\mathrm{\&theta;}}_t\\ \frac{(D-d_v)p_v}{\left|\sin \mathrm{\&theta;}\right|},& \mathrm{\&pi;}+{\mathrm{\&theta;}}_t<\mathrm{\&theta;}<2\mathrm{\&pi;}-{\mathrm{\&theta;}}_t\end{array}\right.,$   公式(二)

其中，θ_t表示探针半径末端与血管壁相交的夹角，D表示血管的直径，d_v表示极坐标和上层血管壁的距离。

对于图2c和图2d中的场景2，探针位于污染区域，并且远离血管。图2c和图2d的圆形轮廓图函数分别定义为：

f(θ)＝Kp_n，  公式(三)

$f\left(\mathrm{\&theta;}\right)=\left\{\begin{array}{cc}{\mathrm{Kp}}_n,& \mathrm{\&theta;}\&Element;{\mathrm{\&Theta;}}_{\mathrm{in}}\\ \frac{d_c{p}_n}{\left|\sin \mathrm{\&theta;}\right|},& (\mathrm{\&theta;}\&Element;{\mathrm{\&Theta;}}_{\mathrm{out}})\&cap;(\mathrm{\&theta;}\&NotEqual;0)\\ d_0,& (\mathrm{\&theta;}\&Element;{\mathrm{\&Theta;}}_{\mathrm{out}})\&cap;(\mathrm{\&theta;}=0)\end{array}\right..$   公式(四)

其中，Θ_in表示那些使探针半径完全在污染区域里的θs集合，Θ_out是Θ_in的补集，d_c表示从探针臂和污染边界的交点到极轴的距离，d₀表示从极点到(θ∈Θ_out)∩(θ＝0)的探针臂和污染边界交点的距离。

对于图2e的场景3，探针位于孤立的结构噪声处，我们定义与d_c相似的变量d_s，圆形轮廓图函数为：

$f\left(\mathrm{\&theta;}\right)=\left\{\begin{array}{cc}\frac{d_s{p}_n}{\left|\sin \mathrm{\&theta;}\right|},& \mathrm{\&theta;}\&NotEqual;0\\ d_0,& \mathrm{\&theta;}=0\end{array}\right.,$   公式(五)

其中d₀的定义同公式(四)中的定义。

对于图2f中的场景4，污染区域靠近血管，靠近血管的探针的圆形轮廓图函数为：

$f\left(\mathrm{\&theta;}\right)=\left\{\begin{array}{cc}{\mathrm{Kp}}_n,& \mathrm{\&theta;}\&Element;{\mathrm{\&Theta;}}_{\mathrm{in}}\\ {\mathrm{Kp}}_v-\frac{d_{\mathrm{cv}}(p_v-p_n)}{\left|\sin \mathrm{\&theta;}\right|},& (\mathrm{\&theta;}\&Element;{\mathrm{\&Theta;}}_{\mathrm{out}})\&cap;(\mathrm{\&theta;}\&NotEqual;0)\\ d_{0c}{p}_n+d_{0v}{p}_v,& (\mathrm{\&theta;}\&Element;{\mathrm{\&Theta;}}_{\mathrm{out}})\&cap;(\mathrm{\&theta;}=0)\end{array}\right.,$   公式(六)

其中，d_cv是从探针臂和血管壁交点到极轴的距离，d_0c表示当θ＝0时从极点到交点的距离，d_0v表示当θ＝0时从交点到探针臂末端的距离。

最后对于图2f中的场景5，结构噪声靠近血管，圆形轮廓图函数为：

$f\left(\mathrm{\&theta;}\right)=\left\{\begin{array}{cc}\frac{d_s{p}_n}{\left|\sin \mathrm{\&theta;}\right|},& (\mathrm{\&theta;}\&Element;{\mathrm{\&Theta;}}_{\mathrm{in}})\&cap;(\mathrm{\&theta;}\&NotEqual;0)\\ \frac{d_s{p}_n}{\left|\sin \mathrm{\&theta;}\right|}+(K-\frac{d_{\mathrm{sv}}}{\left|\sin \mathrm{\&theta;}\right|})p_v,& (\mathrm{\&theta;}\&Element;{\mathrm{\&Theta;}}_{\mathrm{out}})\&cap;(\mathrm{\&theta;}\&NotEqual;0)\\ d_{0s}{p}_n,& (\mathrm{\&theta;}\&Element;{\mathrm{\&Theta;}}_{\mathrm{in}})\&cap;(\mathrm{\&theta;}=0)\\ d_{0s}{p}_n+d_{0v}{p}_v,& (\mathrm{\&theta;}\&Element;{\mathrm{\&Theta;}}_{\mathrm{out}})\&cap;(\mathrm{\&theta;}=0)\end{array}\right.,$   公式(七)

其中d_sv表示交点到极轴的距离，d_0s和d_0v的定义类似于d_0c和d_0v。对于场景4和5，我们没有考虑探针臂穿透血管的情况，因为这些在外面的像素归于背景，它们对f(θ)没有影响。

再次，介绍基于圆形轮廓函数的极性指标的建立。圆形轮廓图极性指标二PCP₂定义为：

${\mathrm{PCP}}_2=E[\underset{m=\mathrm{1,2},...,M}{\max }f\left({\mathrm{\&theta;}}_m\right)-\underset{l=\mathrm{1,2},...,M}{\min }f\left({\mathrm{\&theta;}}_l\right)],$   公式(八)

其中，E(·)表示数学期望函数，该定义强调了圆形轮廓图中各个方向上的差异。方便起见，我们用S1表示图2b中的情景1，S2表示图2c和2d中的场景2，以此类推；对于不同场景的PCP₂值的变动范围如图3所示。

在上述推导中，假定p_v大于p_n，因为物理意义上，真实信号即血管像素的概率分布密度大于噪声像素的概率分布密度。即使p_n＞p_v，公式(二)到(七)仍然相同，以下推导结果也依旧成立。首先，场景1中的PCP₂范围和场景3、4和5中的PCP₂范围存在间隙，正如图中条纹区域所示范围。既然场景1中的探针像素在血管中，场景3到5的探针像素在污染或结构噪声上，那么间隙中的任意值都可以作为区分血管和噪声的阈值(当然，间隙的中间值是最优选择)。这个发现证明了使用PCP₂对血管提取的有效性。同时注意到当探针噪声像素离血管壁很近的时候，间隙会消失。对于与血管像素有相近亮度的像素，所有算法都很难处理。其次，场景2中的PCP₂范围可能与场景1中的PCP₁范围重叠，这使得公式(八)中PCP₂的定义对于场景2失效。

所以，本发明发现场景2中的探针像素都在污染区域，考虑到污染区域(云状)和血管(管状)在空间分布形式上的不同，定义了圆形轮廓极性指标一PCP₁，该定义强调了各向同向性，其具体形式如下：

${\mathrm{PCP}}_1=\frac{{\mathrm{\&sigma;}}^2}{\mathrm{\&mu;}},$   公式(九)

其中，σ²表示方差，μ表示均值。相似于PCP₂，一个大PCP₁值表示血管，一个小PCP₁值表示噪声。因此，可以通过PCP₁的值来判断一个图像区域是云状噪声区域，还是血管区域。对于血管区域，可以直接进行基于PCP₂的几何结构阈值化处理，而对于污染区域，则采用亮度阈值化处理和几何结构阈值化处理两种方式。

最后，介绍本发明的具体实施示例，其具体步骤如图4所示。

步骤一，输入血管图像，对其进行预处理清理背景，统计血管图像中最亮的像素值，取该值的5％作为背景清理阈值，将图像中的像素值低于该阈值的像素置为零；

步骤二，窗口自适应，确定血管区域和污染区域：将整幅血管图像分割成多个子窗口，设置初始窗口尺寸为某个较小的值，如20×20，计算其PCP₂(圆形轮廓图极性指标一，)的值，如果PCP₁的值大于t_h，则扩大窗口尺寸，再次计算PCP₁的值，如此反复，直到PCP₁的值开始减小，从而确定最后的窗口尺寸，这些区域属于血管区域，而PCP₁的值小于t_h的区域则属于噪声污染区域；

步骤三，分别对血管区域和污染区域进行阈值化处理：PCP₁的值小于t_h的图像区域，其中既存在血管，又存在大量噪声，先进行基于多级Otsu算法的亮度阈值化处理，再进行基于PCP₁(圆形轮廓图极性指标一)的几何结构阈值化处理，从而完成血管提取，在采用多级Otsu阈值化方法时，采用三级或四级Otsu阈值化方法，并采用其中最高的亮度值作为阈值，将大于或等于该亮度阈值的像素予以保留，而小于该阈值的像素置为0，从而保证能够尽量去除图像中的污染噪声，几何结构阈值化处理中，当一个像素处所获得的PCP₂的值大于特定阈值t_p＝15时，判断该像素为血管像素并将其置1，反之则为噪声像素并将其置为0；对PCP₁的值大于t_h的图像区域，其中可能只存在血管物体，因此只需进行基于PCP₂的几何结构阈值化处理，其处理方法同上，最后，将所有窗口的阈值化处理结果整合，形成最终的血管提取结果。

图5b和图6b分别给出了两组基于本发明的血管提取结果，图5a和图6a是血管原图，分别为二维血管图像和三维血管图像，图5c到图5g给出了一些其他方法的结果图作为比较。明显可见，本发明既可以有效地去除噪声，又可以有效地完成血管提取。

Claims (4)

1.一种基于圆形轮廓图极性的血管图像处理方法，其特征在于，包括：

首先，输入采集的血管图像，进行背景清理；

其次，对经过背景清理的血管图像通过窗口自适应处理，确定血管区域和噪声污染区域；

判断经窗口自适应处理后的血管图像的圆形轮廓图极性指标一是否小于设定的阈值t_h：若是，则判断为噪声污染区域，若否，则判断为血管区域；

最后，分别对血管区域和污染区域进行阈值化处理：

对噪声污染区域同时采取亮度值阈值化与结构阈值化处理；对血管区域仅采用结构阈值化处理；将噪声污染区域和血管区域的处理结果汇总，得到高质量的血管提取结果；

所述结构阈值化处理基于对血管图像圆形轮廓图极性指标二的处理；

所述的圆形轮廓图极性指标一为PCP₁：

其中，σ²表示该窗口内图像的方差，μ表示该窗口内图象的均值；

所述的圆形轮廓图极性指标二为PCP₂：

${\mathrm{PCP}}_2=E[\underset{m=\mathrm{1,2},...,M}{\max }f\left({\mathrm{\&theta;}}_m\right)-\underset{l=\mathrm{1,2},...,M}{\min }f\left({\mathrm{\&theta;}}_l\right)],$ E(·)表示数学期望函数，f(θ_m)和f(θ_l)分别表示圆周中θ_m和θ_l两个方向上的圆形轮廓函数的值。

2.根据权利要求1所述的一种基于圆形轮廓图极性的血管图像处理方法，其特征在于，包括以下步骤：

步骤一：图像背景处理

输入血管图像，清理血管图像的背景，统计血管图像中最亮的像素值，取该像素值的5％作为背景清理阈值，将图像中的像素值低于该背景清理阈值的像素置为零；

步骤二：窗口自适应

将整幅血管图像分割成多个子窗口，设置初始窗口尺寸值为w_ini，并定义阈值t_h，计算w_ini的圆形轮廓图极性指标一PCP₁；

如果计算得到的PCP₁的值小于阈值t_h，则判断该区域属于血管被噪声覆盖的噪声污染区域；

如果计算得到的PCP₁的值大于阈值t_h，则扩大窗口尺寸，再次计算PCP₁的值，并与上次计算所得的PCP₁的值进行比较，重复计算PCP₁的值与比较操作，直到PCP₁的值减小的幅度大于1为止，确定此时的窗口尺寸，这些区域属于血管区域；

步骤三：基于圆形轮廓图的血管阈值化处理

对判断属于噪声污染区域的图像区域，先进行亮度阈值化处理，再进行基于圆形轮廓图极性指标二PCP₂的几何结构阈值化处理，从而完成血管提取；

对判断属于血管区域的图像区域，进行基于PCP₂的几何结构阈值化处理；

将对所有子窗口进行的阈值化处理结果汇总整合，得到最终的血管提取结果。

3.根据权利要求2所述的一种基于圆形轮廓图极性的血管图像处理方法，其特征在于，步骤三所述的亮度阈值化处理采用三级或四级Otsu阈值化方法，并采用该阈值化方法中最高亮度作为亮度阈值，将大于等于该亮度阈值的像素保存，小于该亮度阈值的像素则置为0。

4.根据权利要求2所述的一种基于圆形轮廓图极性的血管图像处理方法，其特征在于，步骤三所述的基于圆形轮廓图极性指标二PCP₂的几何结构阈值化处理为：

当该图像区域中一个像素处所获得的PCP₂的值大于特定阈值t_p＝15时，判断该像素为血管像素并将其置1，反之则为噪声像素并将其置为0。