基于Hessian算子的多尺度视网膜血管增强滤波方法

基于Hessian算子的多尺度视网膜血管增强滤波方法
丘!立;蒋先刚;熊娟
【摘要】提出一种基于Hessian矩阵的多尺度视网膜血管集成增强滤波算法。

将Hessian矩阵的特征值和特征向量应用于血管特征的响应函数及形态学和非线性各向异性扩散处理的各个环节，选择合理的尺度空间范围及尺度空间增量和调节因子进而达到平滑非线状区域和锐化增强血管区域的效果。

与其他血管增强方法相比，该方法在同等准确率下具有较高的稳定鲁棒性。

%We propose a Hessian matrix-based multi-scale integrated enhancement filtering method for retinal vessels.It applies theeigenvalue and eigenvector of Hessian matrix to response function of vessel feature and morphology as well as all to aspects of nonlinearanisotropic diffusion processing,selects reasonable scale space range and scale space increment and adjustment factor so as to achieve theeffects of smoothing the non-linear region and sharpening and enhancing the vessel pared with other vessel enhancement methods,this one has higher stable robustness with same accuracy.
【期刊名称】《计算机应用与软件》
【年(卷),期】2014(000)009
【总页数】5页(P201-205)
【关键词】视网膜血管;海森矩阵;多尺度滤波;非线性扩散;数学形态学
【作者】丘!立;蒋先刚;熊娟
【作者单位】华东交通大学基础科学学院江西南昌330013;华东交通大学基础科学学院江西南昌330013;华东交通大学基础科学学院江西南昌330013
【正文语种】中文
【中图分类】TP391.41
0 引言
人体血管图像可以提供人体相关组织大量的信息，有助于医生更好地对疾病实施更可靠的临床诊断。

视网膜血管是人体唯一可非创伤性直接观察的微血管，其形态、分布和病变是就医者的高血压、糖尿病、动脉硬化、心血管疾病等的内在表征，视网膜血管的直径、分叉角度以及血管弯曲程度及其变化都是诊断相关疾病的重要指标。

视网膜血管图像对比度较低，而且通常伴有大量的随机噪声，这就需要探索对细节结构损失少的血管图像增强技术。

视网膜血管增强的目的是强调眼底图像中的血管结构形态而同时抑制非重要的图像背景和噪声。

Catte F等用非线性扩散的图像滤波方法在边缘和平坦区域采用不同的标量系数保证边缘不被模糊的前提下使图像得到平滑，Orkisz提出沿血管方向进行中值滤波等形态学的方法而使血管区域得到一定程度的增强。

但这些方法没有考虑到血管图像的局部对比非常弱且非常模糊，基于单一像素邻域的梯度计算的不准确性和形态学操作中的结构元素的有向控制问题，即这些方法只进行了单一尺度的操作，而没考虑多尺度情况下的算法的适应性[1，2]，Frangi等人提出了用基于Hessian矩阵的多尺度相似性测度的方法分析血管得到比较好的加强效果[3]。

在研究眼底图像特征的基础上，本文提出了一种基于Hessian矩阵的多尺度增强滤波的视网膜血管集成增强方法，将Hessian矩阵的特征值和特征向量应用于血
管特征的响应函数和形态学操作和非线性扩散，实验表明这种集成增强滤波算法可达到相当高的准确率并具备较高的鲁棒性，满足血管医学图像分析的临床需要。

1 视网膜血管增强滤波的集成方法
1.1 基于Hessian矩阵的多尺度线状增强滤波器的构造
Hessian方法是一种用高阶微分提取图像特征方向的方法，Hessian方法认为具有最大模的二阶方向导数的方向是垂直于图像特征的方向，与它垂直的方向被认为是平行于图像特征的方向。

对于具备高斯函数构造的线性模型，可以用与直线正交的绝对值较大的二阶导数、沿线方向的绝对值很小的二阶导数来表示，这恰好是二维Hessian矩阵的两个特征值所代表的几何意义。

下面介绍它的基本原理以及在二维图像中Hessian矩阵的应用原理。

视网膜血管都是以线状的形态出现，因此本文主要研究线状增强滤波器，将Hessian矩阵的差分运算与高斯函数结合，通过改变高斯函数的标准偏移量来获得不同尺度σ下的线形增强滤波[4]。

根据高斯函数的卷积性质，尺度空间导数Iab 由输入图像与高斯滤波器的二阶导数的卷积得到[5]：
Iab=I⊗
(1)
其中，高斯函数表达式为：
(2)
其中，σ是高斯滤波器的标准差，也称为空间尺度因子。

根据高斯函数构造的线形模型特点，对于线形结构元素，只有当尺度因子σ与血管的实际宽度最匹配时，此滤波器的输出最大。

一个理想的但考虑到医学血管图像的模糊性的血管结构，需要用比梯度计算更准确的血管边缘处的Hessian二阶算子来改进血管分析的精确性，我们用每一个像素的二维偏导数来构造每一个像素(x, y)的Hessian矩阵，它
的表达式为：
(3)
其中，fxx、fyy、fxy和fyx分别表示像素(x, y)的灰度的四个二阶偏导数。

fxy=fyx，即H是实对称矩阵。

本文的增强滤波器由Hessian矩阵的两个特征值
λ1和λ2来构造。

局部特性分析的窗口矩形的半宽一般取为3σ。

在当前尺度下，
血管的直径小于当前尺度相应的窗口矩形的宽和高，则管状血管的Hessian矩阵
的特征值满足：|λ1|≈0，|λ1|<<|λ2|。

用两个算子来定义这两个特征值的关系：
(4)
在当前尺度下，图像中的点p属于血管区域的响应度定义为：
(5)
其中，通过f(x,y)<ft直接将图中眼球外黑底背景去掉而省去后续的掩膜去背景干
扰运算，β用于调整线状和块状物的区别，c和γ是控制线状物整体平滑的参数，线状响应度函数的调节参数的选择范围是β∈[0.3,2]，c∈[10-5,10-6]，γ∈[3,6]，每一考查像素点p的矩形窗口的半宽一般取为3σ，当不同尺度对应的窗口与血管直径相匹配时，v(p)会得到最大的响应，在多尺度下的一点属于血管区域的响应函数为：
(6)
其中，v(p)的值将用于非线性扩散的扩散因子、形态学操作中的结构元素的选择。

形态学操作中的有方向的结构元素的选择还需考虑最大响应度下尺度σ和Hessian矩阵的特征向量表达的方向。

通过迭代尺度因子σ的选取而得到最大的响应度作为该点的适应度输出，迭代尺
度因子的变化范围为[σmin,σmax]，迭代步长增量为 step，对于定尺度的响应函数，尺度因子σ越小，对直径细小的血管的增强效果越好，而尺度因子σ越大，
对直径较大的血管的增强效果越好。

尺度因子的变化范围越大且迭代步长越小，就能越好地兼顾到各种直径的血管的加强效果。

为了测试尺度因子对增强效果的影响，实验中手动设置固定的尺度因子σ为 0.5，2，6，以及设置尺度因子的取值范围
是[0.5, 6]且迭代步长为0.4，并求4个相应的综合线状适应度。

如取各线状特征响应度的值与图像的灰度值成正比，对图像的各像素取反得到如图1所示的效果。

图1 单尺度和多尺度下Hessian线状特征响应效果对比
实验结果显示，特定尺度下只有相应对的度的部分血管的增强效果明显。

图1(a)
中σ=0.5使尺寸小的血管得到加强，图1(b)σ=2和图1(c)中σ=6使尺寸中等的
血管得到加强，而图1(d)则同时保证了尺寸较小以及中等的血管都得到了较好的
加强效果。

前3种尺度的线状加强都是针对特定尺寸的血管，而只有多尺度和小
步长的综合选取适应度峰值的方法才能使线状区域内的对象都得到适度甄别。

1.2 自适应灰度形态学处理方法
数学形态学处理技术是基于集合论的非线性图像处理方法。

它的最基本的运算算子包括膨胀、腐蚀、开运算和闭运算。

在灰度等级上运用这些算子和它们的组合可对灰度图像的形态和结构进行滤波、分割、填充、提脊等图像加强和分析处理[6，7]。

Tophat运算是形态学处理中一个重要的操作算子，它是将灰色原始图像与经形态学开运算处理后的图像进行差运算，即：
Tophat(f)=f-(f∘b)
(7)
其中f是原始图像，b为结构元素，(f∘b)为灰度开运算，这表示先对图像进行灰度腐蚀再进行灰度膨胀。

对已提取亮区目标的取反图像再进行Tophat运算可以在保
留更多的阴影信息的前提下提升亮区目标的灰度峰值，加强了亮区目标与阴影背景的对比度和增强亮区目标的平滑度。

结构元素的大小和类型与像素点区域的线状适应度取最大响应值时的Hessian矩阵的特征值、特征向量方向和尺度大小有关，
所以需要自适应选择形态学结构而进行处理。

结构元素的大小包括3×3，5×5和
7×7等，而其形状则包括四个主方向加强、四个对角加强和各分角加强等形状。

图2是对视网膜图管状响应图的各种灰度形态学处理结果和比较结果。

图2 视网膜图管状响应图的各种灰度形态学处理结果对比
1.3 非线性扩散处理方法
图像扩散增强的方法来源于热传递扩散的物理概念，扩散过程用Fick定律描述为：▽u|2)·▽u]
(8)
其中，u是图像的强度值，c是扩散张量，▽u表示强度值梯度。

各向异性扩散方
法既考虑扩散的大小又考虑扩散的方向。

扩散张量c由平行梯度和垂直梯度的两
个正交的特征向量e1和e2组成，这个张量的特征值λ1和λ2控制这两个方向的扩散量，当λ2较大时沿血管方向进行平滑和去噪，当λ1是较小正数时，以较小
的扩散程度来维护和保留血管的边界，当λ1是负数时，对血管起锐化作用。

尺度的扩散张量c的计算比梯度的计算更加适应模糊血管图像区域的分析。

从管状物的响应程度V(σ,p)及其相应的图像可以看出，非血管区域的V(σ,p)的值
相对较小，因此对V(σ,p)小于阈值vt的非血管区域用各向同性(λ1=λ2)的高斯滤
波进行平滑。

而在血管内部区域的V(σ,p)较大，像素点的各向异性扩散的扩散量
总是与V(σ,p)有关，而扩散方向由梯度方向的二阶方向导数确定。

对于V(σ,p)值
大于阈值Vt的血管内部的像素的平滑只在梯度正交方向(血管流动方向，即
λ1≈0,λ2≈1)进行，基于管状响应函数的特征向量的选取和血管图像平滑和锐化增
强的非线性各向异性扩散函数按下列公式计算：
λ1=1+(ω-1)v1/sλ2=1+(ε-1)v1/s
(9)
Ut =▽·(D▽U)
=▽▽U)
=▽
(10)
其中，λ1和λ2是对应于最大线状响应的特征值，e1和e2是最大响应函数尺度下的Hessian矩阵的正交的特征矢量，c是扩散张量，取ω∈[20,28]和
ε∈[0.0005,0.003]用于调整响应函数在各特征方向上的扩散程度作用，取s∈[4,7]为调整响应在各特征方向上的特征值的敏感程度。

上述的扩散函数表明通过各点的张量和梯度的相互作用逐步变化而实现有向的平滑扩散。

通过张量c在方向上的控制而实现非线性各向同性或异性扩散，在血管流向区域进行较大前向扩散，在血管边界方向进行后向锐化处理。

考虑到扩散主要在e1和e2两个方向进行，图像扩散方程的离散表达为：
Ut+1[I,J]= Ut[I,J]+Δt(Ut[I-1,J](ct[I,J]+
ct[I-1,J])/2+Ut[I+1,J](ct[I,J]+
ct[I+1,J])/2+Ut[I,J-1](ct[I,J]+
ct[I,J-1])/2+Ut[I,J+1](ct[I,J]+
ct[I,J+1])/2-Ut[I,J](4ct[I,J]+
ct[I-1,J]+ct[I+1,J]+ct[I,J-1]+
ct[I,J+1])/2)
(11)
其中，Ut[I,J]表示t次迭代时(I, J)处的灰度值，Δt为迭代时间增量，一般取为0.003～0.008，它控制每次迭代的扩散程度，上式表明通过主方向上的张量和梯
度的中心差分运算而使相邻的灰度值得到有向的平滑。

图3的非线性扩散处理的
结果比较，图3(b)是以梯度为扩散系数且迭代次数为20次的非线性各向同性扩散结果，在消除背景噪声的同时也使血管细部区域有所减弱，图3(c)、(d)是迭代次
数为10和20的非线性各向异性扩散增强图，它们在消除背景噪声的同时也使血
管区域特别是血管边界区域有所增强和平滑，同时不会随着迭代扩散次数的增加而破坏原有线型结构。

图3 眼底管状响应值经非线性扩散处理结果比较
2 线状增强滤波器算法的实现
将上述方法进行综合应用可实现基于Hessian矩阵的多尺度视网膜血管集成增强
方法。

对于眼球图像而言，尺度因子的范围与血管宽度相适应，步长设置得越小，血管细节越平滑，但耗时越多。

通过反复的实验，我们得出对于Drive的视网膜
血管图像设定尺度因子为[0.5, 6]，迭代步长为0.4效果较好。

考虑到血管原图的
绿色分量都比较低，通过这个特性可对血管图像进行快速的血管绿色对比加强预处理。

基于Hessian矩阵的多尺度视网膜血管滤波增强集成方法实现的步骤如下：
(1) 输入原始彩色图像I，根据图像的中绿色分量增强其血管部分；
(2) 设定尺度因子σ范围[σmin,σmax]和迭代步长Step，对于I的每一个像素Iij，重复(3)-(10)；
(3) 初始化空间尺度因子σ和管状响应V(σ,p)max；
(4) 若尺度σ满足极值则停止循环，跳转(9)；
(5) 计算元素Iij与高斯函数二阶微分的卷积；
(6) 生成Hessian矩阵H，并计算它的两个特征值λ1和λ2及特征向量e1和e2；
(7) 计算增强滤波的管状响应输出值V(σ,p)；
(8) 迭代尺度因子σ=σ+step，跳转(4)；
(9) 输出该像素的最大增强滤波输出值V(σ,p)max及相应的特征值和特征向量；
(10) 以V(σ,p)max及其特征向量等对血管图像进行非线性各向异性扩散；
(11) 以V(σ,p)max、对应的σ、特征值和特征向量选取不同方向和大小的结构元素进行自适应灰度形态学运算。

其中，非线性各向异性扩散通过源程序的调整使得求解过程得到优化的主要代码如下：
class procedure AnisotropicDiffusion4D(out ADstImg: TRealSpace2D; const ASrcImg: TGraySpace2D; const ACount: Integer; const AK, AT: Real);
var
// 定义局部变量…
begin
// 初始化局部变量…
for Index :=0 to ACount - 1 do
//ACount为扩散次数
begin
for I :=1 to W do
//计算X方向的扩散系数
for J :=1 to H do
//计算Y方向的扩散系数
Coefs[I,J]:=1/(1+AK*Sqr(RsltImg[RsltIndex][I,J]-RsltImg[RsltIndex][I-
1,J])+Sqr(RsltImg[RsltIndex][I,J]-RsltImg[RsltIndex][I,J-1]));
for I :=1 to W do
//计算X方向的扩散结果
for J :=1 to H do
//计算Y方向的扩散结果
RsltImg[RsltIndex xor
1][I,J]:=RsltImg[RsltIndex][I,J]+AT*((Coefs[I,J]+Coefs[I-
1,J])*RsltImg[RsltIndex][I-
1,J]/2+(Coefs[I,J]+Coefs[I+1,J])*RsltImg[RsltIndex][I+1,J]/2+(Coefs[I,J]+Coef s[I,J-1])*RsltImg[RsltIndex][I,J-
1]/2+(Coefs[I,J]+Coefs[I,J+1])*RsltImg[RsltIndex][I,J+1]/2-
(4*Coefs[I,J]+Coefs[I-1,J]+Coefs[I+1,J]+Coefs[I,J-
1]+Coefs[I,J+1])*RsltImg[RsltIndex][I,J]/2);
RsltIndex :=RsltIndex xor 1;
//结果下标赋值
end;
for I :=0 to W - 1 do
//将结果复制到输出图像
Move(RsltImg[RsltIndex][I+1,1],ADstImg[I,0],H*SizeOf(ADstImg[0,0]));
// 释放局部变量…
end;
而基于非线性各向异性扩散和形态学的双重调节的响应度也可以用源程序计算出最佳输出，具体实现如下：
procedure VesselFilter2D(out OIm, OScale, ODir: TRealSpace2D; const AIm: TGraySpace2D; const AOp: TVesselOptions);
var
// 定义局部变量…
begin
// 初始化局部变量…
for Index :=0 to LSigma - 1 do
begin
Sigma :=AOp.MinSigma+AOp.SigmaStep*Index;
if Sigma > AOp.MaxSigma then
//尺度限制
Sigma :=AOp.MaxSigma;
GSigma :=Sigma;
//设置高斯函数的参数
TGauss.GetHessian(Hxx, Hxy, Hyy);
Hxx.Filter(Dxx,AIm,mmConv);
// 计算XX卷积
Hxy.Filter(Dxy,AIm,mmConv);
// 计算XY卷积
Hyy.Filter(Dyy,AIm,mmConv);
// 计算YY卷积
for I :=0 to W - 1 do
for J :=0 to H - 1 do
begin
//卷积单量的调整
Dxx[I, J] :=Sqr(Sigma) * Dxx[I, J];
Dxy[I, J] :=Sqr(Sigma) * Dxy[I, J];
Dyy[I, J] :=Sqr(Sigma) * Dyy[I, J];
// 计算特征值和特征向量
GetEigen2D(EV1, EV2, EV1x, EV1y, EV2x, EV2y, Dxx[I, J], Dxy[I, J], Dyy[I, J]);
Angle :=ArcTan2(EV1x, EV1y);
//计算夹角
if EV1=0 then EV1 :=0.000001;
Rb :=Sqr(EV2 / EV1);
//计算特征值之比
S2 :=Sqr(EV1) + Sqr(EV2);
// 计算最大的响应度
if (AOp.BlackObject and (EV1 < 0)) or (not AOp.BlackObject and (EV1 > 0)) then Rslt :=0
else Rslt :=Exp(-Rb / B) * (1 - Exp(-S2 / C));
if OIm[I, J] < Rslt then
begin
OIm[I, J] :=Rslt;
//相应的响应值
OScale[I, J] :=Sigma;
//相应的Sigma
ODir[I, J] :=Angle * W * H;
//相应的方向
end;
end;
// 释放局部变量…
end;
end;
3 实验结果分析
本文中的实验在主机采用CPU Duo 8700，内存为2GB，32位Windows XP操
作系统上进行，开发环境为可视化编程平台Delphi 7。

为了验证本文提出的基于Hessian矩阵的视网膜图像增强方法的有效性，本文在Staal等人提供的Drive公共数据库里的眼底图像上进行实验，这些图像均为565×584的8位灰度图像。

本文实验使用20个视网膜图像样本作为测试用例，每个样本重复5次取平均值。

对于每个样本，该公共数据库有对应的一幅由专家手动分割的图像作为判断有效性的金标准的参考图像，图4是库中40号视网膜图像的处理过程和对比结果。

图4 视网膜血管提取过程及结果比较
尺度因子范围的选择应与视网膜图像的血管宽度相适应，在本文实验中的σ的范
围设置为[0.5, 6]，以包含所有可能的血管宽度，同时设置迭代步长为0.4，对应的窗口矩形的半宽设定为3σ。

取线状响应度函数的各调节参数β=1，γ=5，c=10-6，最小尺度参数σmin将确定可被加强的最细血管，如需消除处理后大直径血管中的断续细小空洞，可通过增加最大尺度参数σmax和进行后续的形态学等处理。

图
4(a)是视网膜血管原图，通过选择合理参数进行绿色对比加强可使血管部分的对比更加明显，图4(b)与图4(c)分别是对应的依绿色分量的增强图像和Hessian响应
度函数增强图，图4(d)是该眼球血管的金标准图像，图4(e)为本文集成方法提取
的血管融合图像图，对比血管原图可以看出，利用本文构造的滤波器可以使得非线形对象得到抑制，同时线形区域得到增强，并客观地反映了血管管径的逐级匀称变化。

本文对眼底图像加强后提取的血管形态区域比金标准图像更加细致地反映了血管管径变化和分布。

对照实验的素材和标准图都来自DRIVE数据库。

对比了多种血管增强方法和本文
采用的集成方法并计算了这些方法的准确性、敏感性及特异性的评价数据如表1
所示，对照方法的评价数据可参考文献[7，8]。

平均准确率表示按多个样本实验的总体平均来看，增强图像中正确识别的血管像素数与正确识别的背景像素数之和被
图像中所有像素相除，ROC 曲线下的面积Az是最常用的评价曲线特征的参数，
Az表示系统同时在敏感性与特异性上的性能表现，Az的值越大表示其结果越可靠。

表1是对DRIVE数据库上的图像进行血管图像增强效果的几种对比结果。

表1 几种血管图像增强效果的对比结果增强方法平均准确率可靠性
Azstaal0.94420.9520zana0.93770.8984非线性扩散0.90570.8767形态学方法0.87230.8223本文方法0.94450.9530
从上面结果分析可知，单独采用形态学、非线性扩散和Hessian特征加强的线状
目标提取都有一定的局限性，而将这些方法进行有效地集成和复合就具备相当的准确率和可靠性，即表明本文方法在敏感性和特异性上都具备较优性能。

4 结语
本文提出了一种基于Hessian矩阵的多尺度视网膜血管的集成增强滤波方法，对
血管区域绿色分量加强的预处理使需提取区域的特征更加明显，继而构造了基于Hessian矩阵的多尺度线形滤波器对视网膜图像的血管进行增强，将表达线状结构的Hessian矩阵特征矢量应用于线状结构的响应函数及自适应的灰度形态学操作
和非线性扩散。

通过中心差分方法而使图像各像素得到有向的平滑及增强处理比单独依靠梯度标量进行的各向同性扩散更适应血管图像的增强处理。

对这些方法的合理结合而获得了较高正确率的眼部血管形态区域的提取和加强，并在同等准确率下有较高的鲁棒性。

在和某医院的合作项目使用中，相比采用Matlab或Photoshop的相关处理功能，本文算法对视网膜的处理使得对某些眼睛疾病的判
断更为准确和快捷，并且已扩展到三维血管的提取和加强中。

进一步的研究工作将放在用模板表达的快速Hessian特征向量提取上。

参考文献
[1] 周琳,沈建新,廖文和,等.基于Hessian矩阵的视网膜血管中心线提取[J].中国制
造业信息化,2011,40(1):46-49.
[2] 李颖超,刘越,王涌天.基于多尺度Hessian矩阵和Gabor滤波的造影图像冠脉中心线提取[J].中国医学影像技术,2007,23(1):133-136.
[3] 刘晏丽,赵卫东,陈宇飞,等.基于Hessian矩阵和区域生长的肝血管树的分割算法研究[J].计算机与现代化,2011(1):113-116.
[4] 李光明,田捷,赵明昌,等.基于Hessian矩阵的中心路径提取算法[J].软件学
报,2003,14(12):2074-2081.
[5] Du Y P,Parker D L.Vessel enhancement filtering in three dimensional MR angiography[J].Journal of Magnetic Resonance Imaging,1995,5(3):353-359.
[6] Orkisz M M.Improved vessel visualization in MR angiography by nonlinear anisotropic filtering [J].Magnetic Resonance in
Medicine,1997,7(6):914-918.
[7] Frangi A F,Niessen W J,Vincken K L.Multiscale vessel enhancement filtering [C]//Medical Image Computing and Computer-Assisted Intervention,LNCS 1496.Berlin:Springer,1998:130-137.
[8] Lorenz C,Carlen I,Buzug T M.Multi-scale line segmentation with automatic estimation of width,contrast and tangential direction in 2D and 3D medical images[C]//CVRMed-MRCAS,LNCS
1205.Berlin:Springer,1997:233-242.
[9] Zana F,Klein J C.Segmentation of vessel-like patterns using mathematical morphology and curvature evaluation [J].IEEE Transactions on Image Processing,2001,10(7):1010-1019.
[10] Staal J J,Abramoff M D,Niemeijer M A,et al.Ridge based vessel segmentation in color images of the retina [J].IEEE Transactions on Medical
Imaging,2004,23(4):501-509.。