CN101236647A - 一种融合上下文信息的数字血管造影图像增强方法 - Google Patents

Abstract

本发明公开了一种融合上下文信息的数字血管造影图像增强方法，本发明首先通过Gabor滤波器将图像分解到不同的频率和方位，再利用血管局部空间的结构一致性，基于一种共圆规则融合血管上下文的信息，并通过一个迭代过程使得一致性的空间成分得到增强，最后通过各向异性扩散平滑去噪，得到较为满意的视觉效果。心血管造影图像通常充满噪声，传统增强技术难以将血管，尤其是细小的血管，从对比度低且背景复杂的造影图像中增强出来。本发明根据目标的全局特征来增强我们感兴趣的成分，将血管这一空间结构性较好的成分从造影图像中增强出来，具有良好的抗局部干扰性。尤其对级数较高、对比度较低的血管末梢表现出好的增强性能。

Description

一种融合上下文信息的数字血管造影图像增强方法

技术领域

本发明属于医学成像技术领域，具体涉及一种融合上下文信息的数字血管造影图像增强方法。

背景技术

图像增强作为一种基本的图像处理技术，其目的是对图像进行加工，以得到对具体应用来说视觉效果更好，更有用的图像。由于具体应用的目的和要求不同，所需要的具体增强技术也可以大不相同。从根本上说，并没有图像增强的通用标准，观察者是某种增强技术优劣的最终判断者。

血管造影技术是一种历史悠久的血管可视化技术。数字减影血管造影(Digital Subtraction Angiography，简称DSA)是在此基础上发展的一种广泛使用的血管可视化技术，在临床已应用20多年，是血管疾病无创诊断与介入治疗手术导航的重要依据，广泛应用于X射线序列成像中的血管的可视化***中。但是由于种种原因，某些部位的数字血管造影图像不适合做剪影处理，因此，血管造影图像在医学图像的诊断分析中扮演着重要的角色。但由于实际造影图像中存在的造影剂的分布不均、衰减，X射线的曝光不均等复杂情况，以及成像环境中各种组织以及相关噪声的影响，其图像质量一般都比较差。因此，血管造影图像处理中的一个关键任务就是进行图像增强，以使血管的生理特征能够更清楚的显示出来，便于医生的分析与诊断。

常用的增强技术根据其处理所进行的空间不同，可分为空间域方法和频域方法。与后者相比，前者是直接在图像所在空间进行的。进一步而言，前者根据对图像的每次处理是对单个像素进行的还是对小的子图像进行的，又可分为基于像素点的和基于模板的方法。基于模板的方法也常称为空间滤波。

目前针对数字血管造影图像的增强技术大多在图像空间直接进行。相对于频域方法，空域方法能更直观地针对图像特点制定改善视觉效果的增强策略。常用的方法有各向异性扩散，哈勃(Gabor)滤波，匹配滤波，形态学方法等等。但是这些方法多是从局部区域的灰度信息出发，而未考虑血管的空间结构以及上下文之间的相互影响，因而增强效果有限。尤其是对于级数较高，对比度极低的血管末梢，很难达到理想的视觉效果，给医生的诊断造成困难。

发明内容

本发明的目的在于提供一种融合上下文信息的数字血管造影图像增强方法，该方法具有良好的抗局部干扰性，能对血管造影图像有效增强，获得良好的视觉效果。

本发明提供的融合上下文信息的数字血管造影图像增强方法，其步骤为：

(1)用哈勃滤波器将灰度反转的数字血管造影图像分解至K个方向L个尺度，K的取值范围为4～16，L根据最大血管宽度W_max与最小血管宽度W_min确定；

(2)将哈勃滤波器的输出作为增强模型的初始输入，并对每个方向和尺度的哈勃滤波响应实施基于共圆约束的增强；其过程为：

(2.1)将中心位于(x，y)的像素在方向α、中心频率ω分量上受到来自环境的增强作用F(x，y，α，ω)表示为：

$F(x,y,\mathrm{\α},\mathrm{\ω})=\underset{(x^{\′},y^{\′})\∈A_e}{\mathrm{\Σ}}\underset{\mathrm{\β}}{\mathrm{\Σ}}{W}_c(x^{\′},y^{\′},\mathrm{\β};x,y,\mathrm{\α})W_d(x^{\′},y^{\′};x,y)R(x^{\′},y^{\′},\mathrm{\β},\mathrm{\ω})$

式中A_e表示满足增强特性的区域，R(x′，y′，β，ω)表示位于(x′，y′)的像素在中心频率为ω、方向为β上的响应，式中的 $W_c(x^{\′},y^{\′},\mathrm{\β};x,y,\mathrm{\α})=\exp (-\frac{Q^2}{{\mathrm{\σ}}_c^2})$ 为曲率加权函数，Q为曲率，且

$Q=\frac{1}{r}=\frac{2}{d}\sin \left(\mathrm{\&phi;}\right)=\left\{\begin{array}{cc}\frac{2}{d}\sin \left|\frac{\mathrm{\&beta;}-\mathrm{\&alpha;}}{2}\right|& 0\&le;2\mathrm{\&gamma;}-\mathrm{\&alpha;}<\mathrm{\&pi;}\\ \frac{2}{d}\cos \left|\frac{\mathrm{\&beta;}-\mathrm{\&alpha;}}{2}\right|& 2\mathrm{\&gamma;}-\mathrm{\&alpha;}<0\mathrm{or}|2\mathrm{\&gamma;}-\mathrm{\&alpha;}\&GreaterEqual;\mathrm{\&pi;}\end{array}\right.$

式中 $d=\sqrt{{(x^{\&prime;}-x)}^2+{(y^{\&prime;}-y)}^2}$

$W_d(x^{\&prime;},y^{\&prime;};x,y)=\exp [-\frac{{(x^{\&prime;}-x)}^2+{(y^{\&prime;}-y)}^2}{{\mathrm{\&sigma;}}_d^2}]$ 为距离加权函数；

(2.2)通过S形变换进行迭代处理：

式中s表示增强作用的上限饱和系数，取值范围为20～80，

的初始值为图像经过哈勃滤波器的输出值；

(2.3)选每个像素所有方向和尺度上的响应值之中最大的一个作为模型的输出，即：

$R(x,y)=\underset{\mathrm{\&omega;},\mathrm{\&alpha;}}{\max }\left[R(x,y,\mathrm{\&alpha;},\mathrm{\&omega;})\right]$

(3)运用各向异性扩散去噪方法对步骤(2)的结果进行平滑去噪处理，得到增强后的图像。

血管作为一种特定的空间结构，在图像中具有一定的空间连续性致性。利用这种空间全局性特征，通过上下文的相互影响来增强血管显然是一种合理的方法。本发明针对现有数字血管造影图像增强技术的缺陷，根据数字血管造影图像的特点，提出一种融合上下文信息的数字血管造影图像增强方法。本发明根据血管的空间结构特点，利用上下文信息的相互影响增强图像。我们根据目标的全局特征来增强感兴趣的成分，将血管这一空间结构性较好的成分从造影图像中增强出来，具有良好的抗局部干扰性。尤其对级数较高，对比度较低的血管末梢具有好的增强效果。

附图说明

图为本发明方法的流程图；

图2为采用哈勃(Gabor)滤波器将一幅血管造影图像分解到4个方向4个尺度上的响应图；

图3(a)为轴特性与模特性；图3(b)为共圆约束增强的空间几何关系；

图4(a)为用本发明的方法增强后，综合各方向和尺度响应的重建图像；图4(b)示出了各向异性平滑去噪的结果；

图5为本发明的实施例中使用上述方法与其他血管造影图像增强方法的一个比较结果，其中，图5(a)表示原始冠脉造影图像；图5(b)传统为各向异性扩散的结果；图5(c)表示基于对比度的各向异性扩散结果；图5(d)表示用模糊形态学方法处理的结果；图5(e)表示用本发明的方法增强的结果。

具体实施方式

下面结合附图和实例对本发明作进一步详细的说明。

实例：

如图1所示，本实例的过程为：

(1)用哈勃(Gabor)滤波器将灰度反转的数字血管造影图像分解至K个方向L个尺度，K的取值范围通常为4～16，尺度的个数L可以根据最大血管宽度W_max与最小血管宽度W_min来确定。

(1.1)本实例中K取值为12，L取值为4，对应的角度

$\mathrm{\&theta;}=\frac{\mathrm{\&pi;}\&times;(k-1)}{K},k=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,K$

每一个角度都代表有可能的增强方向。

(1.2)不同尺度可以通过调整Gabor滤波器的中心频率来获得。在中心频率为ω的尺度上宽度

的血管具有最佳响应。所选择的尺度范围应能使图像中所有的血管都能获得最佳响应。

(1.3)Gabor滤波器的数学表达式为h_e(x，y，ω，θ)＝g(x′，y′)cos(2πωx′)式中：(x′，y′)＝(xcosθ+ysinθ，-xsinθ+ycosθ)，θ表示滤波器的方向，通过在x-y平面旋转可使滤波器达到任何所期望的方位。ω表示中心频率，它决定了滤波器带通区域中心在频域上的位置，通过选取不同的ω可以得到不同的尺度。g(x，y)为如下高斯函数：

$g(x,y)=\frac{1}{2\mathrm{\&pi;}{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y}\exp \{-\frac{1}{2}[{\left(\frac{x}{{\mathrm{\&sigma;}}_x}\right)}^2+{\left(\frac{y}{{\mathrm{\&sigma;}}_y}\right)}^2\left]\right\}$

式中σ_x和σ_y分别是沿着x和y轴的高斯包络的空间常量。σ_x、σ_y分别与频率带宽B_f和方位带宽B_θ有关，关系如下：

${\mathrm{\&sigma;}}_x=\sqrt{\frac{\ln 2}{2}}\frac{1}{\mathrm{\&pi;\&omega;}}\frac{2^{B_f}+1}{2^{B_f}-1}$

${\mathrm{\&sigma;}}_y=\sqrt{\frac{\ln 2}{2}}\frac{1}{\mathrm{\&pi;\&omega;}}\frac{1}{\tan (B_{\mathrm{\&theta;}}/2)}$

在本发明中B_f＝1.25，B_θ＝π/6。

(2)将Gabor滤波器的输出作为增强模型的输入对每个方向和尺度的Gabor滤波响应图实施基于共圆约束的增强。

共圆约束体现了上下文之间相互影响的两种特性：轴特性(axialspecificity)与模特性(modular specificity)。所谓轴特性就是对中心像素具有增强作用的像素应该位于中心像素上下文方向的轴方向；模特性是指这些具有增强作用的像素与中心像素应呈现相似的上下文方向。图3(a)显示了这两种特性。

图3(b)给出了基于共圆约束的曲率计算，A、B分别表示中心像素和环境中一点，其空间坐标分别为(x，y)，(x′，y′)。γ为A、B连线的方位且0≤γ＜π。当B的上下文方位为β时，如果它和A满足共圆约束能形成一个较一致的空间结构，此方位对A的增强作用贡献最大。如果B的主方位与β越接近，那么β方位的响应幅值越大(相对其它方位)，反之，越小。因此我们可以只考虑B在β方位的响应，其幅值的大小反应了它满足共圆约束的程度。

曲率是决定自然一致性结构可探测性的一个重要因素，低曲率结构的显著性通常更高，因此产生的共圆轮廓的曲率越小，其增强作用表现越强烈。

(2.1)中心位于(x，y)的像素在方向α、中心频率ω分量上受到来自环境的增强作用F(x，y，α，ω)为：

$F(x,y,\mathrm{\&alpha;},\mathrm{\&omega;})=\underset{(x^{\&prime;},y^{\&prime;})\&Element;A_e}{\mathrm{\&Sigma;}}\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}{W}_c(x^{\&prime;},y^{\&prime;},\mathrm{\&beta;};x,y,\mathrm{\&alpha;})W_d(x^{\&prime;},y^{\&prime;};x,y)R(x^{\&prime;},y^{\&prime;},\mathrm{\&beta;},\mathrm{\&omega;})$

式中A_e表示满足增强特性的区域(本发明中为与中心的偏离角不超过π/4的区域)，R(x′，y′，β，ω)表示位于(x′，y′)的像素在中心频率为ω，方向β上产生的响应。

(2.1.1)曲率的计算方法如下：A、B分别表示中心像素和环境中一点，其空间坐标分别为(x，y)，(x′，y′)。A的上下文方位为α(0≤α＜π)，B的上下文方位为β时曲率Q表示如下：

$Q=\frac{1}{r}=\frac{2}{d}\sin \left(\mathrm{\&phi;}\right)=\left\{\begin{array}{cc}\frac{2}{d}\sin \left|\frac{\mathrm{\&beta;}-\mathrm{\&alpha;}}{2}\right|& 0\&le;2\mathrm{\&gamma;}-\mathrm{\&alpha;}<\mathrm{\&pi;}\\ \frac{2}{d}\cos \left|\frac{\mathrm{\&beta;}-\mathrm{\&alpha;}}{2}\right|& 2\mathrm{\&gamma;}-\mathrm{\&alpha;}<0\mathrm{or}2\mathrm{\&gamma;}-\mathrm{\&alpha;}\&GreaterEqual;\mathrm{\&pi;}\end{array}\right.$

式中 $d=\sqrt{{(x^{\&prime;}-x)}^2+{(y^{\&prime;}-y)}^2}.$

(2.1.2) $W_c(x^{\&prime;},y^{\&prime;},\mathrm{\&beta;};x,y,\mathrm{\&alpha;})=\exp (-\frac{Q^2}{{\mathrm{\&sigma;}}_c^2})$ 为曲率加权函数，Q为曲率。

$W_d(x^{\&prime;},y^{\&prime;};x,y)=\exp [-\frac{{(x^{\&prime;}-x)}^2+{(y^{\&prime;}-y)}^2}{{\mathrm{\&sigma;}}_d^2}]$ 为距离加权函数。

本发明中曲率和距离加权函数中的高斯标准差分别为σ_c＝0.25，σ_d＝12。

(2.2)通过一个迭代方程实现模型的动态描述，将增强作用向周围扩散，从而实现造影血管图像上血管结构的增强。通过如下S形(Sigmoid)变换实现：

式中s代表神经元的上限饱和系数，本实例中s＝50。图像哈勃(Gabor)滤波器的输出值作

的初始值。

增强模型的迭代次数根据实际数字血管造影图像的质量而定，图像对比度越低，噪声越多，需要迭代的次数越多。

(2.3)选每个像素所有方向和尺度上的响应值之中最大的一个作为模型的输出，重建增强后的图像，也就是：

$R(x,y)=\underset{\mathrm{\&omega;},\mathrm{\&alpha;}}{\max }\left[R(x,y,\mathrm{\&alpha;},\mathrm{\&omega;})\right]$

(3)运用各向异性扩散去噪方法对步骤(2)的结果进行平滑去噪处理，得到增强后的图像。Perona和Malik提出的各向异性的扩散方程：

$\frac{\&PartialD;I(x,t)}{\&PartialD;t}=\mathrm{div}[c(x,y.t)\&dtri;I],I(x,0)=I_0$

式中▽为梯度算子。传导函数c(x，y，t)在扩散的过程中能够随着图像中灰度的变化而变化，自适应地控制扩散的速度。本发明采用以下形式的传导函数：

$c(x,y,t)=\frac{1}{1+{\left[\right||\&dtri;I(x,y,t)||/T]}^2}$

其中T称为传导参数，‖▽I(x，y，t)‖为梯度的幅值。

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

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

图5所示的反映了在本发明的典型实施例中使用上述方法与其他血管造影图像增强方法的一个比较结果，其中，图5(a)表示原始冠脉造影图像；图5(b)传统为各向异性扩散的结果(P Perona，J Malik.Scale-space and edgedetection using anisotropic diffusion[J].IEEE Tansactions on PAMI，1990，12(7)：629-639)；图5(c)表示基于对比度的各向异性扩散结果(李抱朴，桑农等.一种新的血管造影图像增强方法.电子学报，2006，34(4)：695-697)；图5(d)表示用模糊形态学方法处理的结果(K Sun，N Sang，et al.Extraction ofVascular Tree on Angiogram with Fuzzy Morphological Method.InternationalJournal of Information Technology，2005，11(9)：119-127)；图5(e)表示用本发明的方法增强的结果。

Claims (1)

1、一种融合上下文信息的数字血管造影图像增强方法，其步骤为：

(1)用哈勃滤波器将灰度反转的数字血管造影图像分解至K个方向L个尺度，K的取值范围为4～16，L根据最大血管宽度W_max与最小血管宽度W_min确定；

(2)将哈勃滤波器的输出作为增强模型的初始输入，并对每个方向和尺度的哈勃滤波响应实施基于共圆约束的增强；其过程为：

(2.1)将中心位于(x，y)的像素在方向α、中心频率ω分量上受到来自环境的增强作用F(x，y，α，ω)表示为：

$F(x,y,\mathrm{\&alpha;},\mathrm{\&omega;})=\underset{(x^{\&prime;},y^{\&prime;})\&Element;A_e}{\mathrm{\&Sigma;}}\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}{W}_c(x^{\&prime;},y^{\&prime;},\mathrm{\&beta;};x,y,\mathrm{\&alpha;})W_d(x^{\&prime;},y^{\&prime;};x,y)R(x^{\&prime;},y^{\&prime;},\mathrm{\&beta;},\mathrm{\&omega;})$

式中A_e表示满足增强特性的区域，R(x′，y′，β，ω)表示位于(x′，y′)的像素在中心频率为ω、方向为β上的响应，式中的 $W_c(x^{\&prime;},y^{\&prime;},\mathrm{\&beta;};x,y,\mathrm{\&alpha;})=\exp (-\frac{Q^2}{{\mathrm{\&sigma;}}_c^2})$ 为曲率加权函数，Q为曲率，且

$Q=\frac{1}{r}=\frac{2}{d}\sin \left(\mathrm{\&phi;}\right)=\left\{\begin{array}{cc}\frac{2}{d}\sin \left|\frac{\mathrm{\&beta;}-\mathrm{\&alpha;}}{2}\right|& 0\&le;2\mathrm{\&gamma;}-\mathrm{\&alpha;}<\mathrm{\&pi;}\\ \frac{2}{d}\cos \left|\frac{\mathrm{\&beta;}-\mathrm{\&alpha;}}{2}\right|& 2\mathrm{\&gamma;}-\mathrm{\&alpha;}<0\mathrm{or}2\mathrm{\&gamma;}-\mathrm{\&alpha;}\&GreaterEqual;\mathrm{\&pi;}\end{array}\right.$

式中 $d=\sqrt{{(x^{\&prime;}-x)}^2+{(y^{\&prime;}-y)}^2}$

$W_d(x^{\&prime;},y^{\&prime;};x,y)=\exp [-\frac{{(x^{\&prime;}-x)}^2+{(y^{\&prime;}-y)}^2}{{\mathrm{\&sigma;}}_d^2}]$ 为距离加权函数；

(2.2)通过S形变换进行迭代处理：

式中s表示增强作用的上限饱和系数，取值范围为20～80，

的初始值为图像经过哈勃滤波器的输出值；

(2.3)选每个像素所有方向和尺度上的响应值之中最大的一个作为模型的输出，即：

$R(x,y)=\underset{\mathrm{\&omega;},\mathrm{\&alpha;}}{\max }\left[R(x,y,\mathrm{\&alpha;},\mathrm{\&omega;})\right]$

(3)运用各向异性扩散去噪方法对步骤(2)的结果进行平滑去噪处理，得到增强后的图像。

Images