一种近红外静脉特征提取方法研究

杨 刚 徐学军

(长沙理工大学物理与电子科学学院 湖南 长沙 410000)

0 引 言

无论是静脉穿刺,还是中医针灸,临床上,医护人员主要通过目测或经验定位血管,但对部分血管过细、脂肪较厚或血管塌陷的特殊患者(例如儿童、体胖者、老人等),血管难以辨别,严重影响静脉穿刺成功率[1]。而在中医针灸中,医生需要将针灸针一根一根施于患者的相应穴位中,在施针的时候,医生用双眼很难精确地去确定每个穴位的具体位置,同时还得避免扎错地方,若是扎到了血管或其他要害部位,还会危害患者的健康[2-3]。为了使静脉能够被准确识别,人们通过光照显影、红外显像、静脉投影等辅助方式,提高静脉精准定位。而近年来根据相关研究,主要利用近红外成像技术,来区分静脉和其他组织,有较好的实施效果和适用性[1],其静脉图像提取系统如图1所示。

图1 近红外静脉图像提取系统

然而,基于近红外成像技术采集的静脉图像,由于受到光照的不均匀性和个体差异性的影响,导致近红外图像出现对比度低、边缘模糊和细节信息丢失等缺陷,造成了图像分割后出现很多虚假的静脉信息。因此,如何增强静脉图像的对比度以及将静脉纹路与皮肤背景分割开来成为了研究的热点。鲁周迅等[4]提出了基于子直方图均衡化算法,结合NiBlack分割方法对图像进行静脉提取。林剑等[5]提出了基于Hessian矩阵的静脉图像分割方法,通过Hessian矩阵的迹初次滤除了非静脉区的像素点。王定汉等[6]提出基于静脉灰度值特征的图像分割与对比度增强算法。Malutan等[7]提出了将LLBP描述子与形态学算子相结合用于提取静脉特征,有效地提取出了静脉纹理特征。以上算法都取得了一定的效果,但是也存在一些不足。因此,探索新的静脉图像的对比度增强和分割方法是非常必要的。

由于人体静脉特征网络具有复杂性和方向高度随机性,本文算法在基于方向滤波器的基础上,首先提出多方向Gabor滤波算法进行静脉图像增强,然后将每个方向所得的静脉图像进行重组得到增强后的图像。最后利用马尔可夫随机场分割方法对静脉网络进行分割,该静脉特征提取流程如图2所示。实验结果表明,该方法能够得到可靠有效的静脉网络特征。

1 静脉图像增强算法

静脉没有虹膜那么丰富的纹理,也不像指纹具有较多的细节点特征,其主要体现在具有复杂变化关系的血管网络结构中[8]。人体静脉特征通常要通过静脉图像区域分割来提取,而静脉图像分割的好坏会直接影响后续的静脉特征提取质量。因此,本文提出一种基于多方向的Gabor滤波器的静脉图像增强算法。

Gabor滤波是一种众所周知的特征提取方法,已在机器视觉领域得到广泛研究和应用[9]。Gabor 滤波器对于方向和角度的表达类似于人类视觉系统对于方向和角度的表达[1],而在本文中引入的多方向Gabor滤波器目的是能够检测不同方向的静脉血管。在空间域中,一般二维(2D)Gabor滤波器被定义为包含复正弦信号的高斯函数[10]。其表达式如下:

(1)

(2)

式中:G(x,y)是像素的坐标;σx、σy分别为沿着x和y方向的高斯包络的标准偏差,xθ、yθ表示二维Gabor滤波器旋转θ角后更新后的x、y坐标;*表示卷积运算符;f0代表二维Gabor滤波器的空间中心频率。

Gabor滤波函数又分为实部和虚部,经研究发现,Gabor滤波器的实部也叫做偶对称 Gabor滤波器,特别适合图像中脊线的检测,且增强了Gabor滤波器处理时间的有效性[1]。偶对称Gabor滤波器如式(3)所示。

(3)

基于前人的经验可知, 当滤波器的方向角θ个数超过8个时, 对静脉的增强效果很差,故在选择方向个数的时候, 通常会选择4个或6个方向,一般来说, 利用Gabor滤波器的多个方向可以挖掘出更多的静脉信息。然而, 选择方向过多时会产生更多的冗余信息,这会大大降低Gabor 滤波器的性能,选择方向较少时又不能挖掘出更多的静脉信息[8]。所以本文选择6个方向,其中i(i=1,2,3,4,5,6)表示6个方向的索引值,对应的6个方向θi为π/6、2π/6、3π/6、4π/6、5π/6。fi表示根据第i个方向的偶数对称Gabor滤波器的空间中心频率,其值的选择决定了输出滤波图像的质量。图3显示了6个方向的偶对称Gabor滤波器的空间响应。

图3 6个方向的偶对称Gabor滤波器的空间响应

设I(x,y)为一幅近红外成像设备获得的静脉图像,将每一方向对应的偶对称Gabor滤波函数与原图像进行卷积运算,获得滤波后的静脉图像Oi(x,y)如下:

(4)

由于近红外静脉图像中静脉线比皮肤区域暗,因此,选择6幅图像中每个位置像素点灰度值最低的点构成滤波后的静脉图像,如图4所示。

图4 滤波后的静脉图像

2 MRF图像分割

众所周知,静脉图像分割提供了具有分离的静脉区域的图案分割图像,这对于特征提取很重要[11]。对于静脉图像而言,分割使得静脉纹路和皮肤背景区分开来,以便于特征提取等后续处理。传统的图像分割方法主要分为基于像素的方法、基于区域的方法、基于边界的方法、基于模型的方法和多种组合算法[12]。而本文在使用Gabor滤波器静脉图像增强后,采用MRF随机场用于静脉图像分割。实验表明,该算法在静脉图像分割中取得很好的分割效果。

基于马尔可夫随机场理论是指图像中的像素点与其相邻像素相关,即相邻像素的强度值变化不大。任何数字图像都是由离散的像素集合构成的,这些像素可以被建模为一个随机场。图像中的每个像素都是一个站点,每个站点都有一个标签X,它通常表示一个像素的强度值[13]。

由上文可知静脉纹路暗,皮肤背景亮,可分别记为标签“0”和“1”,这样将大小为m×n的整幅图像分为了两类,表示为:

W={w1,w2,…,wm×n}

式中:W为所得点集;w(w=0,1)表示图像各点像素值为s所属类。

为了使得分割图像更精准,我们需要将像素值更好地分为所属类。P(w|s)表示在已知像素信息s的条件下属于一个标签的概率,也叫后验概率。因此,该过程就是通过最大化后验概率来寻找像素点s所属最优标签w。由于后验分布是一个条件分布,我们可以用Bayes规则来估计它,根据这个规则可得后验概率表达式[14]:

(5)

式中:P(w)是先验概率,在此认为其符合MRF模型;P(s|w)称为似然函数;P(s)是一个常量。

根据Hammersley-Clifford[15]定理,用 Gibbs 随机场的概率密度函数描述马尔可夫随机场,则先验概率P(w)表示为:

(6)

(7)

本文采用双点势团的 8 邻域系统,像素点t为像素点邻域s内的点。其中β为耦合系数,控制着势团内两点相互作用的强弱,一般取值在[0.5,1.0]区间之间,耦合系数β越小,对图像的分割越细腻,随着β的增大,对细碎纹理的检测能力逐渐降低[16]。

而P(s|w)为似然函数,表示在已知分类标签w情况下,点像素灰度值为s的概率。因此,我们需要将图像随机初始化分类为w,并将属于w类的点找出来。一般认为每一类里面的所有点分布是具有不同均值和标准差的高斯函数[16],对每一类像素点建立高斯密度函数。对于静脉图像而言,我们得到两类高斯密度函数。

那么已知像素点灰度值则能计算出该像素点分别属于“0”“1”类的概率值P(s|w0)、P(s|w1)。由式(5)可知,当P(s|w)×P(w)达到最大值时,图像分割效果最理想。具体分割算法描述如下:

(1) 获取Gabor滤波增强后的静脉图像。

(2) 使用k-means聚类算法初始化图像标签“0”和“1”。

该矿体为浸染(星点)状辉钼矿，矿体呈灰白色—浅肉红色—褐红色，地表出露长度为142 m，厚度1.72～2.86 m，产状310°∠69°。矿体规模不大，沿走向露头不好。矿石金属矿物主要有黄铁矿、褐铁矿、辉钼矿、钼华，蚀变主要为硅化。脉石矿物主要为二长花岗岩及碎裂石英岩。围岩为二长花岗岩，矿体与围岩界线不明显，沿矿体两侧围岩产生钾化、硅化蚀变。矿石钼品位0.038%～0.137%。

(3) 统计每个标号场8领域的相同标签个数,获得图像标号场的先验概率。

(4) 计算似然函数P(s|w)。

(5) 求取最大化后验概率,并更新标号场。

(6) 判断是否达到最大迭代次数,若没有,重复步骤(2)-步骤(4),否则,算法结束。

本文在 Windows 10、64位操作系统下,通过 MATLAB 2014a实现上述算法。相关函数来自MATLAB自带函数库,部分主要代码如下:

Class_num=2;

%设置分类数

Max_iter=20;

%最大迭代次数

Label=kmeans(img(:),Class_num);

for i=1:Class_num

Label_i=i*ones(size(label));

Temp=～(label_i-label_u)+～(label_i-

Label_d)……;

p_c=(i,:)=temp(:)/8;

%计算最大后验概率

end

for j=i:Class_num

MU=repmat(mu(j),size(img,1)*

size(img,2),1);

p_sc(j,:)=1/sqrt(2*pi*sigma(j))*

exp(-(img(:)-MU).^2/2/sigma(j));

end

%计算似然概率P(s|w)

[M,label]=max(log(p_c)+log(p_sc));

%寻找最大后验概率

分割后的静脉图像如图5所示。

图5 分割后的静脉图像

3 实验结果分析

为了验证本文方法的有效性,实验采用Intel(R)Core(TM)i5-5200U CPU @2.20 GHz 2.20 GHz,内存为8 GB的PC机,在MATLAB 2014 a环境下完成的。

本实验选择了3幅实验室自行设计的近红外静脉成像设备采集到的噪声未经处理的静脉纹路高度随机性的分辨率为640×480像素的静脉图像进行测验。在进行Gabor滤波时,为避免图像欠增强和产生更多的冗余信息,因此将这3幅低复杂度背景的静脉图像选用六方向分别为30°、60°、90°、120°、150°、180°,x、y方向的高斯包络偏差σx、σy均为2的Gabor滤波器进行静脉图像增强,如图6所示,经实验所得静脉图像得到了显著增强,克服了GHE等图像增强算法出现的局部过增强和欠增强等问题,同时也很好地抑制了处理过程中产生的噪声、边缘伪影等信息,为后续进行静脉网络特征提取提供了好的静脉环境。

图6 实验所得静脉图像

其次,得到增强后的静脉图像用马尔可夫随机场进行分割,耦合系数β越小,对静脉图像的纹理与背景分离越明显,随着β的增大,许多静脉细节信息丢失。在进行中值滤波去毛刺处理时,经多次实验对比,选择窗口为5×5的滤波模板时去毛刺效果最好。将本文方法与其他三种常用静脉分割方法进行比较,结果如图7所示。Otsu分割算法对静脉图像做分割时,分割后的图像仍然保留了过多的噪声。使用Niblack分割算法得到的二值图像中对于皮肤背景信息较复杂的静脉图像而言,容易引进较多的伪静脉信息。K-means聚类算法对于较高质量的静脉图像能够提取到很好的静脉信息,且残留的噪声较少,对于复杂背景的静脉图像分割依旧有噪声存在。

图7 各种算法对比

本文提出的静脉图像分割方法对低背景复杂度的静脉图像做分割时,能够提取到更好地静脉信息,没有噪声、伪静脉信息的干扰。而对于背景信息相对复杂的静脉图像而言,能够更好地抑制噪声信息的产生。采用主观评价法通过人眼来比较不同的分割结果图像,其稳定性较差,难以定量描述,在应用中受到很大限制。为了更好地评价本文提出的静脉特征提取模型性能,因此选择无监督评价指标中的区域内一致性指标对图像分割性能进行客观评价。而区域一致性指标又分为最大对比度评价和方差评价一致性,本文使用最大对比度评价(Zeb)[18],结果如表1所示。可以看出使用本文算法所得Zeb值相比几种传统分割算法而言,其值更小,所得静脉分割区域内一致性更好,所得分割图像效果更好。因此本文方法提取出来的静脉较为完整,提取的静脉结构非常适合原始图像的血管拓扑结构。

表1 不同分割算法的最大对比度评价指标

4 结 语

本文为了解决利用近红外成像设备采集到的静脉所存在的问题,提出一种基于多方向Gabor滤波器的静脉图像增强方法,系统有效地解决了静脉信息的噪声大、对比度低、图像清晰度不够的实际困难,主要采取了两种技术手段。(1) 根据静脉的纹理特征,采用了六维度Gabor滤波器进行滤波,获得了满意的增强静脉图像;(2) 为了实现静脉特征提取和针灸穴位精准定位,用马尔可夫随机场对增强后的静脉进行二值分割,并采用中滤波取毛刺处理。通过范例实验结果表明,本文算法能够有效精准地提取皮下近红外静脉特征和针灸穴位,可靠地增强静脉图案质量,为中医针灸、西医穿刺、生物特征识别等领域提供了新思路。