基于HESSIAN增强和形态学尺度空间的视网膜血管分割

于　挥　王小鹏

(兰州交通大学电子与信息工程学院　甘肃 兰州 730070)



基于HESSIAN增强和形态学尺度空间的视网膜血管分割

于挥王小鹏

(兰州交通大学电子与信息工程学院甘肃 兰州 730070)

眼底视网膜血管的走向、弯曲度、分叉度等性状分析已成为医学上诊断全身血管性疾病的重要手段。采集到的眼底图像常存在光照不均匀等现象，利用传统的血管分割方法难以对微小血管进行检测。为此提出一种基于改进Hessian矩阵增强和形态学尺度空间的分割方法。首先利用高斯函数构建多尺度Hessian增强滤波器，采用新型的血管相似性函数对血管网络进行对比度增强，同时平滑图像以减轻噪声；然后利用改进的Top-hat变换尺度空间从背景中提取血管，并引入形态学重建方法进一步突出血管像素，消除伪边缘及孤立点噪声；最后使用二次阈值化方法实现血管的最终分割。仿真结果表明，改进的分割方法在保证大血管脉络准确分割的同时，能够较好地实现微小血管分割。

视网膜血管Hessian增强尺度空间形态学分割

0　引　言

近年来，由于临床诊断全身血管类疾病的需要，国内外学者针对视网膜血管的增强和分割提取进行了大量研究。Soares等[1]采用基于Gabor小波变换方法，利用大量预分割标准图像的特征信息进行分割；姚畅等[2]运用分布式遗传算法与Otsu相结合提取对比度较低的边界，通过训练算法实现分割；Oliveira等[3]提出了Hessian矩阵理论，对视网膜血管图像提取具有良好的效果；Frangi等[4]利用Hessian矩阵构造二维空间和三维空间的血管相似性函数以增强视网膜血管；Zhou等[5]在Frangi的基础上提出了一种基于Hessian矩阵的滤波方法，但由于其血管相似性函数设定的局限性，容易合并相邻的平行血管；Bauer等[6]针对这种情况引入了梯度矢量流进行分割，该方法较好地区分了相邻血管，并对微小血管有一定的检测能力。

针对视网膜血管呈现管状和线性的结构，Hessian算法对微小血管和低对比度血管有良好的检测和增强能力。经过传统的Hessian增强后的图像虽增大了血管图像的对比度，但仍存在灰度不均的现象。形态学是一种常用的图像分割方法，使用基于多尺度结构元素的Top-hat变换分割视网膜血管，在小血管等细节处理上有较好的效果[7,8]。本文提出一种基于新型血管相似性函数的Hessian矩阵增强，并结合改进的形态学尺度空间视网膜血管分割方法。首先对视网膜图像绿色分量进行改进的Hessian增强，然后利用改进的形态学尺度空间方法从增强后的图像中分割出视网膜血管。方法改进了传统的形态学Top-hat变换以消除灰度不均匀的背景，首次运用形态学重建去除分割结果中的噪声和伪边缘。最后进行阈值化得到最终的血管结构图像。

1　方法流程

在传统的分割基础上，本文将改进的Hessian算法与形态学尺度空间相结合，原始图像经过预处理后，首先构造改进的多尺度Hessian增强滤波器对视网膜血管进行对比度增强，然后利用改进的形态学尺度空间对增强后的血管图像实现分割。改进的Hessian增强的关键在于提出了一种新型的尺度因子集合及血管相似性函数；形态学尺度空间分割的关键在于提出了一种改进的Top-hat变换尺度空间方法，并加入形态学开重建，以达到眼底背景及小细节噪声的去除的目的。方法具体流程如图1所示。

图1　方法流程图

2　改进的Hessian矩阵血管增强

2.1Hessian矩阵原理

视网膜血管图像整体呈现出类似树杈的形状，其拓扑结构主要由线段构成，且具有很强的连续性。Hessian矩阵是一个由多元函数的二阶偏导数构成的方阵，其特征值和特征向量能很好地描述这种线性和管状结构。图像某一点A的局部特性由其泰勒展开式表示：

I(A+ΔA)≈I(A)+ΔAT▽I(A)+ΔATH(A)ΔA

(1)

在2D图像中，▽I(A)为图像在点A的梯度，H(X)为点A的Hessian矩阵：

(2)

其中fxx(A)、fyy(A)、fxy(A)、fyx(A)分别代表图像在点A处的二阶微分。由于图像元素是离散的，离散Hessian矩阵微分表现为差分运算。X方向上的二阶偏微分可表示为：

fxx=f(x-1,y)+f(x+1,y)-2f(x,y)

(3)

Y方向上的二阶偏微分为：

fyy=f(x,y-1)+f(x,y+1)-2f(x,y)

(4)

X，Y方向上的混合偏微分为：

fxy=f(x+1,y+1)+f(x,y)-f(x+1,y)-f(x,y+1)

(5)

对于连续函数f(x,y)，二阶偏导数的求导顺序没有区别，即fxy=fyx，因而H为实对称矩阵，具有两个特征值λi，i=1,2。其中较大特征值对应的特征向量与血管垂直，曲率最大；而较小特征值对应的特征向量与血管平行，是血管的真正走向，曲率也最小。Hessian矩阵的两个特征值λ1、λ2可由式(6)计算：

(6)

2.2多尺度Hessian算法

单一尺度的Hessian增强算法对视网膜血管直径变化较大的图像的增强效果较差。因此使用多尺度Hessian算法来实现增强，将高斯核函数引入其差分运算，改变高斯核函数的标准差δ以获得多尺度下的增强结果。将图像I与高斯核函数的二阶偏导数Gxx卷积可将fxx改写为空间导数Ixx：

(7)

其中δ为空间尺度因子，二维高斯函数表达式为：

(8)

同理Iyy、Ixy也由此方法求出。卷积后的结果不仅构造了基于δ的多尺度滤波器，由于高斯函数的性质，滤波器同时平滑了图像，消除了噪声的影响。空间导数Ixx、Iyy、Ixy的值与空间尺度因子δ的平方成反比。将二阶偏导数Gxx乘以δ2后再进行卷积，可实现比较多尺度下滤波器的输出值，则式(7)变化为：

(9)

(10)

2.3新型血管相似性函数

在二维图像中，Frangi等人[4]利用两个因子RB、S来构造视网膜血管的相似性函数：

(11)

其中参数β一般设为0.5，c一般为矩阵最大范数的一半，而：

(12)

式中λ1和λ2分别为Hessian矩阵的两个特征值，‖H‖F则表示矩阵的范数，D为图像的维数，此处D=2。

理想情况下，血管像素的λ1远大于λ2，输出响应可取得较大值；背景像素的λ1和λ2均很小，输出响应很弱。而Hessian矩阵特征值的输出响应对局部特性极其敏感，在实际临床中为保护患者，已逐渐使用眼底拍照彩色图像代替眼底造影图像，采集到的图像多有灰度不均及背景噪声环境复杂等情况。一般来说，噪声像素的灰度值起伏很大，RB因子在增强血管像素的同时，也对复杂背景有部分增强。对此提出一种新型的血管相似性函数对血管目标进行检测，表达式为：

(13)

该函数修改了RB因子，将特征值模的大小与矩阵范数相关联，利用范数调节两个特征值模之间的比值，并取消了参数β。这样，相似性函数较前具有更强的自适应性，且输出响应可减弱灰度不均匀及噪声的影响，增强效果较为理想。

2.4软件实现方法

改进算法可通过Matlab仿真工具的图像处理工具箱IPT(ImageProcessTool)实现，构建Hessian增强函数。先求解二维高斯核函数的二阶偏导，乘以尺度因子的平方，构建Hessian矩阵求出特征值，代入改进的相似性函数得到输出响应，核心代码如下：

Function Ih=HE(I,s)

Dxx(p,q)=1/(2*pi*sigma^4)*(x(p,q)^2/sigma^2-1)*exp(-(x(p,q)^2+y(p,q)^2)/(2*sigma^2));

%高斯核函数二阶偏导

Dsxx=Dxx*sigma^2;

%与尺度因子平方相乘

Ixx=imfilter(I,Dsxx,′replicate′);

%卷积求空间导数

H{i,j}=[Ixx(i,j) Ixy(i,j);Ixy(i,j) Iyy(i,j)];

%构建Hessian矩阵

E{i,j}=eig(H{i,j});

%计算H{m,n}矩阵的特征值

tzz=E{i,j};

%提取特征值

r1(i,j)=tzz(1,1);

%提取特征值r1

r2(i,j)=tzz(2,1);

%提取特征值r2

Ih(ic,jc)=exp((-1)*(abs(r1(ic,jc))-S(ic,jc))/(abs(r2(ic,jc))-S(ic,jc)))*(1-exp((-1)*S(ic,jc)^2/(2*c^2)));

%由相似性函数求输出响应

end

该HE函数的输出变量为Hessian增强输出响应，输入变量I为待增强图像，s为尺度因子，使用for循环对所有像素点进行逐个处理。多尺度增强通过多次调用该函数实现，代码如下：

Ih(:,:,i)=HE(Igc,se(i));

%尺度为se(i)的Hessian增强

Ihm=max(Ih,[],3);

%多尺度增强的最大值，各增强图像叠加

3　形态学尺度空间分割

形态学尺度空间是非线性的形态学运算，相对于单尺度的线性形态学运算，尺度空间运算能更有效地保留感兴趣的边缘像素，消除多余的微小细节，同时可以避免图像平滑过程中造成的边缘模糊和轮廓移位。现已有膨胀腐蚀尺度空间、开闭尺度空间、Top-hat变换空间和重建尺度空间等多种形态学尺度空间[9-11]，每种空间的使用都存在一定的优势和局限性。

因此，本文提出了一种结合两种尺度空间的分割方法。首先将改进的Top-hat变换与形态学尺度空间相结合以提取增强图像的血管目标，满足灰度不均、直径大小不一的血管目标提取要求；然后利用开重建尺度空间对提取的血管边缘图像进行重建，经过重建后的视网膜血管图像保留了微小血管轮廓，消除了孤立点噪声，同时保证了血管边缘轮廓的位置不发生偏移。

3.1改进的Top-hat变换尺度空间

基于形态学开运算的传统Top-hat变换[9]是一种高效的目标提取滤波器，定义为：

g(x,y)=f(x,y)-(f(x,y)∘b)

(14)

其中f(x,y)为输入图像，∘表示形态学开运算，b为结构元素。结构元素是Top-hat变换处理的关键，采用合理的结构元素可以很好地消除灰度不均匀的背景，提取目标图像。

传统Top-hat变换是基于形态学开运算的，经过开运算后图像的灰度值比原图灰度值低或不变，原图减去开运算图像后将检测出原图所有的灰度值变化点。然而血管图像的背景较为复杂，这些变化点不全是血管像素点，还有可能存在噪声点，视网膜图像的不均匀光照使此方法误检出噪声点的概率大大增加。为此，使用一种改进的 Top-hat变换，采用基于开运算和闭运算相结合的Top-hat变换，其表达式如下：

E(f)=f-min((f·Sc)∘S0;f)

(15)

其中Sc为闭运算的结构元素，So为开运算的结构元素，符号·和∘分别代表形态学闭运算和开运算。首先对图像进行闭运算，然后对闭运算后的图像做开运算，所得结果与原图进行对比，求出最小值。得到的结果图像只有在边缘像素点处与原图像不同，从原图中减去该图像即可得到图像的边缘目标。

同时引入尺度空间，克服了单尺度Top-hat变换对复杂视网膜血管特性检测能力较弱的缺点，其定义为：

Se(i)=strel(′disk′,i)i∈[imin,imax]

(16)

其中i表示结构元素的大小，imin、imax分别为尺度空间范围的最小值和最大值，空间步长为N。变换使用圆盘型结构元素，其具有的各向相同性、旋转不变性可以防止形态学滤波过程中图像的特征发生畸变。各尺度结构元素对图像处理后，利用最值法合并结果图像以作为尺度空间的输出图像。改进后的Top-hat变换尺度空间在检测图像边缘的同时，较好地消除了复杂的背景信息，对视网膜图像的分割与识别有很好的效果。

3.2开重建尺度空间

传统的血管分割并未涉及到开重建尺度空间，而经过改进Top-hat变换尺度空间后的图像分割效果虽比传统的分割效果明显，但仍含有部分孤立点噪声。血管的灰度值是连续变化的且跨度较大，大血管的灰度值高，小血管的灰度值低，而孤立点噪声的灰度值一般包含在血管的灰度级范围内，随后的阈值化很难在保留微小血管的前提下将其去除。形态学开重建[10]相对于开运算具有更准确的还原性，因此在阈值化之前，使用开重建尺度空间对Top-hat变换尺度空间的结果图像进行重建，达到保留微小血管并消除非血管的孤立点噪声的目的。

形态学开重建算法建立在测地膨胀理论之上，对于灰度图像I(x,y)和参考图像r(r取I-1)，形态学测地膨胀定义为：

(17)

(18)

3.3阈值化

本文使用基于Otsu的特殊阈值法进行二值化，通过两次阈值化和形态学重建进行实现。首先利用Otsu法求出最小类间方差作为level值，对血管图像进行阈值化。然后手动设定level值进行阈值化。最后对比两种阈值化的level值，采用小level值结果图像作为掩膜Mask，对大level值结果图像进行膨胀重建，表达式为：

H(I)=RT[l,255](I)[T[L,255](I)]

(19)

其中L和l分别代表较大和较小的level值，T[l,255]表示阈值化中保留大于l小于255的值，RT(I)表示用T(I)做掩膜进行重建操作。

3.4软件实现方法

形态学尺度空间分割通过Top-hat变换分割、形态学重建去噪、阈值化三个步骤实现：

(1) 构建改进的Top-hat函数，并进行多尺度分割，核心代码为：

function Itp=TopHat(I,se)

Ic=imclose(I,se);

%先闭运算

Ico=imopen(Ic,se);

%再开运算

Imin(i,j)=min(I(i,j),Ico(i,j));

%与原图比出较小值

Itp=imsubtract(I,Imin);

%图像减法运算

end

Itp(:,:,i)=TopHat(V,se(i));

%Top-hat分割

Itpm=max(Itp,[],3);

%多尺度取最大值

函数中se=strel(′disk′,i)为圆盘型结构元素。

(2) 分割完成后，使用OBR函数进行形态学开重建：

function Fobr=OBR(F,se)

Fe=imerode(F,se);

%结构元素为se的形态学腐蚀

Fobr=imreconstruct(Fe,F);

%重建

end

Iobr(:,:,i)=OBR(Itpmb,se(i));

Iobrm=max(Iobr,[],3);

%多尺度重建最大值

(3) 最后对重建后的血管图像进行二次阈值化操作，图像分别使用level1和level2为掩膜进行阈值化，再将两者结果图像进行膨胀重建得到最终分割结果图：

level=graythresh(Iobrm);

%使用Otsu法获取level值

BW1=im2bw(Iobrm,level);

%一次阈值化

level2=DL(Iobrm);

%使用DL函数手动设定level值

BW2=im2bw(Iobrm,level2);

%二次阈值化

if level>level2

BW=imreconstruct(BW1,BW2);

%膨胀重建

else

BW=imreconstruct(BW2,BW1);

end

4　仿真结果及分析

为验证方法有效性和分割精确度，选取Drive公共眼底图像库的两幅标准眼底测试图像进行仿真测试。仿真使用Intel Core i3处理器、2 GB内存的PC机，在Matlab 2014a平台上进行分析。源图像为Cannon CR5照相机拍摄的眼底彩色图像16_test和19_test，分辨率为584×565，包含一定程度的背景噪声及光照不均匀现象，视网膜血管由粗到细呈树状结构分布。

图2(a)为原始眼底图像，由于其绿色分量含噪声较少且对光照亮度相对不敏感，一般提取其绿色分量作为分割源图像。对绿色分量图像(如图2(b))取反并作双精度化预处理，以满足Hessian滤波增强的需要。血管直径从几像素到十几像素不等，为使大血管脉络和小血管分支均得到增强，选取尺度因子δi的最小值δmin为0.2，最大值δmin为2，量化灰度级N为10可达到全面增强的效果。血管的相似性函数中RB因子由特征值λ1和λ2计算获得，参数c的取值取决于图像的灰度范围，这里设为Hessian矩阵最大范数的一半，同样由特征值计算得出。图2(c)给出了尺度因子δ=1的单尺度增强效果图，图2(d)为本文方法的改进多尺度Hessian增强滤波图。可以看出，单尺度增强效果不佳，肉眼难以观察，多尺度因子对粗细不同的全局血管均有较好的增强效果。

图2　单尺度增强和多尺度增强对比图

图3(a)为Frangi的多尺度相似性函数增强结果，β取0.5，c取15；图3(b)为本文提出的新型的相似性函数增强结果，c取10。对比图像可以发现，本文所构建的相似性函数对细小血管网络有更好的增强性。

图3　相似性函数对增强效果的影响

图4(a)为原始图像经过提取绿色分量、双精度化预处理后的图像，进行多尺度Hessian增强(图4(b))后，利用改进的Top-hat变换尺度空间进行分割。设定圆形结构元素的大小i从1到20递增，步长为1，分割图像如图4(c)所示。图4(d)为直接对图4(c)进行阈值化，获得的图像含有大量的孤立点噪声，丢失了部分细小血管。图4(e)为对图4(c)进行本文提出的尺度空间开重建后，进一步增强血管网络的同时去除了伪血管的眼底边缘轮廓。图4(f)为最终分割图像，使用Otsu法计算得出的掩膜level值为0.45，根据图像的灰度范围设定level2值为0.2。与专家手动分割的金标准图4(g)相比，本文在精确检测出大血管轮廓的基础上，较好地检测出了微小血管。

图4　不同算法的分割结果图

使用图像配准算法将最终结果图像与原始图像进行配准，仿真结果如图5所示。(a)分别为两幅Drive测试图像的最终分割图像，(b)为原始图像的绿色分量，(c)为分割图与绿色分量图经配准后的结果图像。可以看出，良好的方法保持了血管原有位置，且未导致血管变粗致使直径发生变化。

在眼底图像中视网膜血管由线段组成，呈连续的树杈状由粗到细分布，这种图像特征是本文方法实现有效分割的关键。仿真中参数的设置也很重要，空间因子和结构元素尺度范围应与血管直径的大小相匹配，以达到全局血管的增强和有效分割。为定量分析方法的准确性，使用基于像素的统计测量评价方法，将专家手工勾画的血管网络图像作为标准参考图像。计算本文方法分割结果图像的真阳性TP(True Positive)、假阳性FP(False Positive)、真阴性TN(True Negative)、假阴性FN(False Negative)。TP表示血管点被正确识别为血管点的像素数；FN表示血管点被错误检测为非血管点的像素数；TN表示非血管点被正确识别为非血管点的像素数目；FP表示非血管点被错误检测为血管点的像素数。将它们组合为三个测度来评价分割的有效性，分别为精确度Acr(accuracy)、灵敏度Sns(sensitivity)和特异性Spc(specificity)[12]：

(20)

(21)

(22)

理想情况下，精确度最大为1代表所有像素被正确分类，灵敏度最大为1代表着所有血管像素被正确标记，特异性最大为1代表所有背景像素被正确标记。表1给出了两幅图像分别使用不同方法的精确度、灵敏度和特异性，本文方法与未增强的尺度空间分割和Hessian增强的Top-hat分割进行相比，精确度、灵敏度和特异性分别为0.913、0.887、0.906。对比数据可以看出，其结果图像与专家手工绘制的金标准图像最为接近。

表1　使用不同方法的图像统计测量值

5　结　语

视网膜血管分割是医学影像处理的重要分支，在计算机辅助下改善眼底图像质量、突出血管脉络，为临床医生进行疾病诊断提供了帮助。本文在传统的图像分割基础上提出改进的Hessian增强和形态学尺度空间相结合的方法，以达到准确分割视网膜血管图像的目的。利用Hessian矩阵对线状物体的敏感性，提出一种新型的相似性函数以实现血管像素增强；形态学尺度空间的引用，有效地从灰度不均图像背景中提取出视网膜血管结构，消除了噪声及伪边缘，使分割结果更加精确。仿真结果表明，该方法在准确分割大血管的同时，对微小血管有较好的分割检测能力，更大程度上接近眼科专家手工勾画的血管图像，提高了医生的疾病诊断效率，减少诊断错误的发生，并提供了重要的参考价值。

[1] Soares J B V,Leandro J J G,Cesar R M,et al.Retinal vessel segmentation using the 2-D gabor wavelet and supervised classification[J].IEEE,Transactions on Medical Imaging,2006,25(9):1214-1222.

[2] 姚畅,陈后金.一种新的视网膜血管网络自动分割方法[J].光电子·激光,2009,20(2):274-278.

[3] Oliveira W S,Tsang I R,Cavalcanti G D C.Retinal vessel segmentation using Average of Synthetic Exact Filters and Hessian matrix[C]//IEEE International Conference on:Image Processing (ICIP),2012(19):2017-2020.

[4] Alejandro F Frangi,Wiro J Niessen,Koen L Vincken,et al.Multiscale vessel enhancement filtering[J].Medical Image Computing and Computer-Assisted Intervention-Miccai98,1998,1496(1):130-137.

[5] Zhou Jinghao,Chang Sukmoon,Meraxas D,et al.Vascular structure segmentation and bifurcation detection[C]//IEEE International Symposium on Biomedical Imaging:From Nano to Macro,2007:872-875.

[6] Bauter C,Bischof H.A novel approach for detection of tubular objects and its application to medical image analysis[J].Pattern Computer Science,2008,5096(1):163-172.

[7] Bai Xiangzhi.Image enhancement using multi scale image features extracted by top-hat transform[J].Optics & Laser Technology,2012,44(2):328-226.

[8] Jinsung Oh,Heesoo Hwang.Feature enhancement of medical images using morphology-based homomorphic filter and different evolution algorithm[J].International Journal of Control,Automation and Systems,2010,8(4):857-861.

[9] 廖苗.一种新的视网膜血管图像增强方法[J].光电子·激光,2012,23(11):2237-2242.

[10] 刘剑秋,阮秋琦.形态学重建滤波器的研究与应用[J].通信学报,2002,23(1):116-121.

[11] 王爱华.基于Curvelet变换和形态学的视网膜血管分割[D].武汉:华中科技大学,2013.

[12] 张石,董建威,佘黎煌.医学图像分割算法的评价方法[J].中国图象图形学报,2009,14(9):1872-1880.

RETINAL VESSELS SEGMENTATION BASED ON HESSIAN ENHANCEMENT AND MORPHOLOGICAL SCALE SPACE

Yu HuiWang Xiaopeng

(SchoolofElectronicsandInformationEngineering,LanzhouJiaotongUniversity,Lanzhou730070,Gansu,China)

Characters analysis in regard to the trend, curvature and bifurcation of retinal vessels in fundus has become the important means of systemic vascular diseases diagnosis in medicine science. Because of most collected fundus images has the phenomenon of light unevenness, it is difficult to use traditional vessel segmentation methods to detect the micro vessels. Therefore we proposed a segmentation algorithm, it is based on the improved Hessian matrix enhancement and morphological scale space. First, by using Gauss function the algorithm constructs multi-scale Hessian enhanced filter, and uses a novel vascular similarity function to carry out the contrast enhancement on vascular network, while smoothes the image to weaken noise as well; then it extracts the vessels from background using an improved Top-hat transformation scale space, and introduces morphological reconstruction method to further highlight the vascular pixels and to eliminate the pseudo-edges and the noise of outliers; finally the algorithm uses secondary thresholding approach to realise final vessel segmentation. Simulation experimental results showed that while ensuring the accurate segmentation of great vessels and choroid, the improved segmentation method can better realise the segmentation of micro vessels.

Retinal vesselHessian enhancementScale spaceMorphological segmentation

2015-01-20。国家自然科学基金项目(61261029)；兰州市科技计划项目(2013-4-63)；金川公司预研基金项目(JCYY20130 09)。于挥，硕士生，主研领域：图像处理。王小鹏，教授。

TP391.4

A

10.3969/j.issn.1000-386x.2016.08.045