基于Hessian特征的视网膜血管图像的增强滤波算法

蒋先刚，熊 娟，丘赟立，范德营

（华东交通大学基础科学学院，江西南昌 330013）

视网膜血管的直径、分叉角度以及血管弯曲程度及其变化都是诊断相关疾病的重要指标。视网膜血管图像对比度较低，同时伴有大量的随机噪声，这就需要探索对细节结构损失少的血管图像的增强技术。视网膜血管增强的目的是强调眼底图像中的视网膜血管结构形态并同时抑制非重要的图像背景和噪声，从而加强图像的判别和识别。Orkisz提出沿血管方向进行中值滤波等数学形态学的方法使血管区域得到一定程度的增强，这些方法只是基于单一像素邻域的梯度计算和形态学操作中的结构元素的有向控制，而没有考虑多尺度情况下算法的适应性［1-2］。Frangi等［3］人提出了基于Hessian矩阵的多尺度相似性测度的方法，采用这种方法分析血管能得到较好的加强效果。

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

1 视网膜血管增强滤波的集成方法

1.1 基于Hessian矩阵的多尺度线状增强滤波器的构造

Hessian是一种用高阶微分提取图像特征方向的方法。Hessian方法认为具有最大模的二阶方向导数的方向是垂直于图像特征的方向，与它垂直的方向被认为是平行于图像特征的方向。对于用高斯函数构造的线性模型，可以用与直线正交的绝对值较大的二阶导数和沿线方向的绝对值很小的二阶导数来表示。将Hessian矩阵的差分运算与高斯函数结合，通过改变高斯函数的标准偏移量就可获得不同尺度σ下的线形增强滤波［4］。由高斯函数的卷积性质可知，尺度空间导数Ixy是由输入图像与高斯滤波器的二阶导的卷积得到［5］

式中：σ是高斯滤波器的标准差，也就是空间尺度因子。只有当尺度因子σ与血管的实际宽度最匹配时，滤波器的输出才最大。对于一个理想的血管结构并且要兼顾到医学血管图像的模糊性，需要用比梯度方法更准确的处理血管边缘处的Hessian二阶算子来改进血管分析的精确性，我们用每一个像素的二阶偏导数来构造每一个像素(x，y)的Hessian矩阵，其表达式为

式中：H表示Hessian矩阵；fxx，fyy，fxy和fyx分别表示二维图像f( )x，y在指定区域内的4个二阶偏导数，其中fxy=fyx。由于Hessian是实对称矩阵，所以，可以用Hessian矩阵的两个特征值λ1和λ2来构造我们的增强滤波器。

Hessian矩阵的2个特征值λ1和λ2中，幅值最大的特征值对应的特征向量代表着该点曲率最大的方向，而幅值最小的特征值对应的特征向量代表该点曲率最小的方向，表1总结了二维数据中Hessian矩阵特征值和探测形状之间的对应关系，表1中的H，L分别表示特征值大与小，+/-表示特征值的正、负。

表1 二维情况下各种可能形状结构的Hessian矩阵特征值的关系表Tab.1 Hessianm atrix eigenvalues relationship of every possible shape under two-dimension conditions

对于一幅包含各种形状的灰色图像用Hessian矩阵的特征值而对不同的几何要素滤波效果如图1所示。

图1 按基本配比的λ1和λ2特征值的滤波效果Fig.1 The filtering effectsof eigenvalues λ1 and λ2 atbasic ratio

图1（a）是图像原图，图1（b），（c），（d）是取尺度σ=1.0 时的点状、线状和块状物的增强滤波处理效果。由于滤波过程中的尺度固定且参数选择只考虑了依据表1中的一般阈值限定条件而作，自然地对块状滤波过程中仍然保持了点、线状物体的一些特性，故对某一形状的滤波效果需要取决于尺度大小、λ1和λ2的合理权重及调和函数的选配。对于眼部血管的特殊情况，如取局部特性分析的窗口矩形的半宽一般取3σ，在当前尺度下，血管的直径小于当前尺度对应的窗口矩形的半宽，则管状血管的Hessian矩阵的特征值满足：|λ1|≈ 0，|λ1|≪|λ2|。下面我们用两个算子RB，S来定义这两个特征值的关系

式中：λ1和λ2表示Hessian矩阵的两个特征值；‖H‖F表示Hessian矩阵的F范数。在当前尺度下，图像中属于血管区域的点p的响应度可定义为

式中：V(σ，p)是点p的响应度；β是用于调整线状物和块状物区别的参数；c，γ是控制整体平滑的参数，通过f(x，y)＜ft直接将图中眼球外黑底背景去掉，这样可以省去后续的掩膜去除背景的干扰运算，线状响应度函数的调节参数β，γ，c的选择范围分别是0.2～3，2～6，10-5～10-6，每一考查像素点p的矩形窗口的半宽一般取3σ，当不同尺度对应的窗口与血管直径相匹配时，v(p)会得到最大的响应，在多尺度下的某一点属于血管区域的响应函数为

v(p)的值将用于非线性扩散的扩散因子和形态学操作中的结构元素的选择。形态学操作中的带有方向和大小的结构元素的选择还需考虑最大响应度下的尺度因子σ和Hessian矩阵的特征向量表达的方向。

图2 单尺度和多尺度下的Hessian矩阵线状特征响应效果对比Fig.2 Responseeffect contrastof Hessian linearmatrix under single-scaleandmulti-scale conditions

由实验结果可知，在相同尺度下只有相应宽度的部分血管的增强效果明显，图2（a）中，σ=1.0使直径小的血管得到加强，图2（b）中，σ=3和图2（c）中，σ=5使直径中等的血管得到加强，而前3种尺度的线状加强都是针对相应尺寸的血管的线状加强，而图2（d）的效果表明使用多尺度、步长小的综合方法来取的适应度峰值才能使线状区域内的对象都得到适度的增强。

1.2 非线性各向异性扩散处理方法

在非线性热扩散方程计算中，通过设计合适的扩散系数来控制扩散方程的扩散行为，可使得在平滑图像的同时能够保留甚至增强图像的特征信息。将Hessian矩阵的特征值和特征向量应用到非线性扩散中，能得到更好地增强效果。扩散过程用Fick定律描述为：

式中：div是散度算子；u是图像像素的灰度值；D是扩散张量；∇u表示强度值梯度；t表示物理中热扩散时的时间，程序设计中用迭代步长表示。各向异性扩散方法既考虑扩散的大小又考虑扩散的方向。扩散张量D由平行梯度和垂直梯度的两个正交的特征向量e1，e2组成，这个张量的特征值λ1，λ2控制这两个方向的扩散量，当λ2较大时，将沿血管方向进行平滑和去噪，当λ1是较小的正数时，以小的扩散量可维护和保留血管的边界，当λ1是负数时，对血管壁进行锐化处理。

从管状响应值V(σ，p)和管状响应图可以看出，在非血管区域的像素点的V(σ，p)值很小，对于V(σ，p)值小于阀值vt的非血管区域实施各向同性（λ1=λ2）的高斯滤波平滑。而在血管内部区域的像素点的V（σ，p）值较大，像素点的各向异性扩散的扩散量总与V(σ，p)有关，而扩散方向由梯度方向的二阶方向导数确定。对于V(σ，p)值大于阀值Vt的血管内部的点的平滑只在梯度正交（血管流动方向）方向（λ1≈0，λ2≈1）进行，基于管状响应函数的特征向量选取和血管图像平滑和锐化增强的非线性各向异性扩散函数按下列公式计算

式中:λ1，λ2是对应于最大线状响应的特征值，e1，e2是正交的两个特征向量，对应于最大响应函数尺度下的Hessian矩阵的特征向量，D是扩散张量，ω，ε用于调整响应函数在各特征方向上的扩散程度，ω一般取值为20～30，ε一般取值为0.000 8～0.00 2，s为调整响应在各特征方向上的特征值的敏感程度，它的取值范围为3～8。上式的扩散函数表明通过各点的张量和梯度的相互作用逐步变化而实现有向的平滑扩散。图3是非线性扩散处理的结果比较，图3（b）是以梯度为扩散系数且迭代次数为20的非线性各向同性扩散结果，在消除背景噪声的同时也使血管细部区域有所消弱，图3（c）,（d）是迭代次数分别为13和22的非线性各向异性扩散图，它们在消除背景噪声的同时也使血管区域特别是血管边界区域有所增强和平滑，同时不会随着迭代扩散次数的增加而使线型结构模糊。下图是对Drive图库中的2号图做的眼底管状响应值经非线性扩散处理结果对比图。

图3 眼底管状响应值经非线性扩散处理结果比较Fig.3 Thenonlinear processing comparison resultsof fundus tubular response value diffusion

2 线状增强滤波器的算法实现

将上述方法进行综合应用可实现基于Hessian特征的视网膜血管图像的增强滤波算法。对于眼底图像而言，尺度因子的范围与血管宽度相适应。通过反复的实验，我们得出对于Drive的视网膜血管图像设定尺度因子为［1.0，5］，迭代步长为0.5的效果最好。考虑到血管原图的绿色分量都比较低，由这个特性可先对血管图像进行快速的血管绿色对比加强预处理。基于Hessian矩阵的多尺度线状增强视网膜血管增强滤波集成方法实现的步骤如下：

Ⅰ输入原始彩色视网膜血管图像I，设定尺度因子σ范围[σmin，σmax] 和迭代步长Step，对于I的每一个像素Iij，重复Ⅱ～Ⅶ；

Ⅱ初始化空间尺度因子σ和管状响应V(σ，p)max；

Ⅲ若尺度σ满足停止条件，则跳转9）；

Ⅳ计算元素Iij与高斯函数二阶微分的卷积，生成Hessian矩阵H，并计算其两个特征值λ1，λ2及其特征向量e1，e2；

Ⅴ计算增强滤波的管状响应输出值V(σ，p)；

Ⅵ迭代尺度因子σ=σ+step，跳转Ⅳ；

Ⅶ尺度因子迭代结束，输出该像素的最大增强滤波输出值V(σ，p)max、特征值和特征向量；

Ⅷ图像像素遍历结束；

Ⅸ以V(σ，p)max特征向量等对血管图像进行非线性各向异性扩散，以V(σ，p)max对应的σ、特征值和特征向量选取不同方向和大小的结构元素进行自适应灰度形态学运算。

3 实验结果分析

图像滤波实验在Windows XP操作系统上进行，主机采用CPU Intel I3，内存为2GB，开发环境为可视化编程平台Delphi7。为了验证本文提出的基于Hessian矩阵的多尺度线状加强的视网膜图像增强方法的有效性，系统在由Staal等人提供的Drive公共数据库的眼底图像上进行实验。该公共数据库的每一眼底视网膜图都有对应的一幅由专家手动分割的图像作为判断有效性的金标准的参考图像，实验原始图像的尺寸为565×584，图4是库中2号和21号视网膜图像的处理过程和对比结果。

图4 视网膜血管提取过程及结果比较Fig.4 Theextraction processand the comparison resultsof retinalvascular

尺度因子的选择应与视网膜图像的血管宽度相适应，在本文实验中的σ的范围设置为［1.0，5］，包含了所有可能的血管宽度，同时设置迭代步长为0.5，对应的窗口矩形的半宽设定为3σ。取线状响应度函数的各调节参数β=1.5，γ=6，c=10-6，最小尺度参数σmin将确定可被加强的最细血管，如需消除处理后大直径血管中的断续细小空洞，可通过增加最大尺度参数σmax和进行后续的形态学和各向异性扩散等处理。形态学处理的结构元素的大小和类型与像素点区域的线状适应度最大响应值的Hessian矩阵的特征值、特征向量方向和尺度大小有关，需要自适应选择形态学结构进行处理。两组图中的图4（a）是视网膜血管原图，通过选择合理参数进行绿色对比加强可使血管部分的对比更加明显，图4（b）与图4（c）分别是对应的依绿色分量的增强图像和Hessian响应度函数增强图，图4（d）是该眼底血管的金标准参考图像，图4（e）为本文集成方法提取的血管融合图，与血管原图的对比可以看出，利用本文构造的滤波器可以使得非线形对象得到抑制的同时线形区域得到增强，并且客观地反映了血管管径的逐级匀称变化。

对照实验的素材和金标准都是采用DRIVE数据库的样本和参考数据，比较可仅在图像的脊线上进行。评价体系中的平均准确率表示按多个样本实验的总体平均来看，增强图像中正确识别的血管像素数与正确识别的背景像素数之和与被图像中所有像素数相除，ROC曲线下的面积Az是最常用的评价曲线特征的参数，Az表示系统同时在敏感性与特异性上的性能表现，Az的值越大表示其结果越可靠。表2是对DRIVE数据库上的图像进行血管图像增强效果的几种对比结果。

从上面结果分析可知，单独采用形态学、非线性扩散和Hessian特征加强的线状目标提取都有一定的局限性，效果也不是很好，而将这些方法进行有效地集成和复合就具备相当的准确率和可靠性。由图4中的原图可看出眼底血管原图图像成橙红色，它对我们需要增强的视网膜血管有很大的影响，我们采用将形态学、非线性扩散和Hessian特征进行有效的结合有效的克服了这方面的影响，并且大大提高了视网膜血管增强的准确率。由此可以看出通过这些方法合理有序的操作以及实施而获得了较高正确率的视网膜血管形态区域的加强和提取，并在同等准确率下能得到较强的鲁棒性。

表2 几种血管图像增强方法的效果对比Tab.2 The contrastof severalenhancementmethods for vascular images

4 结束语

提出了一种基于Hessian特征的视网膜血管图像的增强滤波算法，通过对Hessian矩阵特征值的管理利用而获得不同形状的增强方法，并构造了基于Hessian矩阵的多尺度线形滤波器对视网膜图像的血管进行集成增强，将表达线状结构的Hessian矩阵特征向量应用于线状结构的响应函数及自适应的灰度形态学操作和非线性各向异性扩散，通过中心差分法而使图像各点的张量和梯度得到有向的平滑增强。通过这些方法的合理有序的操作以及实施而获得了较高正确率的视网膜血管形态区域的加强和提取，并在同等准确率下能得到较强的鲁棒性。

［1］ 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.

［2］YOSHINOBU S，CARL F，ABHIR B，etal.Tissue classification based on 3d local intensity structure for volume rendering［J］.Visualization and ComputerGraphics，IEEETransactions，2000，6（2）：160-180.

［3］ FRANGIAF，NIESSENW J，VINCKEN K L．Multiscale vessel enhancement filtering［C］//Medical Image Computing and Computer-Assisted Intervention，LNCS 1496，Berlin：Springer，1998：130-137.

［4］游嘉，陈波，基于Hessian 矩阵的多尺度视网膜增强算法［J］.计算机应用，2011，31（6）：113-116.

［5］罗菁，林树忠，倪建云，等，基于Hessian矩阵的指纹细节点提取方法［J］.光电工程，2008，35（11）：134-138.

［6］李光明，田捷，赵明昌，等，基于Hessian矩阵的中心路径提取算法［J］.软件学报，2003，14（12）：2074-2081.

［7］许燕，胡广书，尚丽华，等，基于Hessian矩阵的冠状动脉中心线的跟踪算法［J］.清华大学学报，2007，47（6）：889-892.

［8］郑轶，蔡体健.稀疏表示的人脸识别及其优化算法［J］.华东交通大学学报，2012，29（1）：10-14.

［9］蒋先刚，许伦伦，赵莹.基于三维各向异性扩散的图像平滑及三维重构效果分析［J］.华东交通大学学报，2010，27（3）：78-82.