CN109584210B - 多光谱三维静脉显像系统 - Google Patents

Abstract

本发明公开了一种多光谱三维静脉显像系统，显像系统包括：两路CCD相机，用于在两个不同的波段上通过软件触发同时采集静脉图像；图像增强和分割模块，用于将两幅图像进行图像增强和图像分割处理之后，得到分割后图像；静脉图像立体匹配模块，用于将分割后图像进行图像立体匹配，得到静脉的着色深度图；多光谱静脉影像融合模块，基于局部特征残差与统计匹配，对采集到的双波段静脉图像经小波变换分解得到的低频系数和高频系数运用不同准则进行处理来确定高频融合系数，再通过小波逆变换得到融合的静脉影像。本发明可以快速、准确地找到静脉，方便医用定位静脉，进行有效地扎针。

Description

多光谱三维静脉显像系统

技术领域

本发明属于生物医学成像与辅助诊疗领域，具体涉及一种多光谱三维静脉显像系统。

背景技术

在生物医学成像与辅助诊疗领域，通过对多光谱医学图像信息的获取与分析对于可疑病灶部位的定位与定性具有重要意义，将有助于实现相关疾病的早期诊疗，降低误诊率。静脉显像技术可以快速、准确地找到静脉，再将静脉图像投影回手背的对应的静脉处，可方便医用定位静脉，进行有效地扎针。

传统的多光谱显像系统存在着集成度不高的缺点，少有小型化和微型化的成功范例；静脉显像系统也存在近红外图像中静脉分布不明显问题；传统的系统通过单相机单波段成像直接获取静脉图像，存在静脉图像轮廓不明显，信息不丰富，无法从中获取深度信息的问题。

发明内容

本发明的目的在于提供一种多光谱三维静脉显像系统。

实现本发明目的的技术解决方案为：一种多光谱三维静脉显像系统，包括：

两路CCD相机，用于在两个不同的波段上通过软件触发同时采集静脉图像；

图像增强和分割模块，用于对两幅静脉图像进行图像增强和图像分割处理，得到分割后图像；

静脉图像立体匹配模块，用于将分割后图像进行图像立体匹配，得到静脉的着色深度图；

多光谱静脉影像融合模块，对采集到的双波段静脉图像经小波变换分解得到低频系数和高频系数，基于局部特征残差与统计匹配确定高频融合系数，再通过小波逆变换得到融合的静脉影像。

一种多光谱三维静脉显像方法，包括以下步骤：

在两个不同的波段上通过软件触发同时采集静脉图像；

对两幅静脉图像进行图像增强和图像分割处理，得到分割后图像；

采用基于低纹理区域的静脉图像立体匹配算法，将分割后图像进行图像立体匹配，得到静脉的着色深度图；

对采集到的双波段静脉图像经小波变换分解得到低频系数和高频系数，基于局部特征残差与统计匹配确定高频融合系数，再通过小波逆变换得到融合的静脉影像。

与现有技术相比，本发明的显著优点为：(1)本发明采用基于低纹理区域的静脉图像立体匹配算法，得到静脉的着色深度图，可以快速、准确地找到静脉，方便医用定位静脉，进行有效地扎针；采用基于局部特征残差与统计匹配的多光谱静脉影像融合算法，得到融合的静脉影像，静脉图像轮廓明显，信息更加丰富；(2)本发明搭建了多光谱三维静脉显像系统，将传统的单相机显像替换成了双相机显像，并对获取的多波段的高质量的手背静脉图像进行后端计算；与传统的单相机显像相比，双相机显像得到着色深度图像拥有单相机显像没有的多光谱信息和三维信息，轮廓明显，信息丰富，并实现了成像视场和投影视场的重叠；在此基础上，针对不同成像区域，本发明中的快速切换滤光片装置可以更换不同波段的滤光片，选取最优波段的成像效果图作为后续嵌入式计算模块的输入。

附图说明

图1为多光谱三维静脉显像系统匹配算法流程图。

图2为双波段静脉融合结果图。

图3为投影的深度图。

图4为立体匹配算法流程图。

图5为基于像素位置的代价聚合窗口图。

图6为代价聚合示意图。

图7为融合流程图。

具体实施方式

本发明提出一种多光谱三维静脉显像系统，包括：

两路CCD相机，用于在两个不同的波段上通过软件触发同时采集静脉图像；

图像增强和分割模块，用于对两幅静脉图像进行图像增强和图像分割处理，得到分割后图像；

静脉图像立体匹配模块，用于将分割后图像进行图像立体匹配，得到静脉的着色深度图；

多光谱静脉影像融合模块，对采集到的双波段静脉图像经小波变换分解得到低频系数和高频系数，基于局部特征残差与统计匹配确定高频融合系数，再通过小波逆变换得到融合的静脉影像。

进一步的，静脉成像过程中，被照物反射LED光源的两束光，分别通过808nm窄带滤光片和940nm窄带滤光片的滤光进入两路CCD相机，在第一帧图像的获取时刻采用软件触发，同时控制两个相机进行图像采集，获得两个波段的静脉图像；808nm窄带滤光片中心波长为808nm，半带宽为30nm，940nm窄带滤光片中心波长为940nm，半带宽为30nm；

进一步的，两个相机通过连接杆固定在光学平板上，两个CCD相机等高，极线距离为20～40cm。

进一步的，图像增强采用OTSU算法设定阈值后，扩大静脉图像中静脉和皮肤的差别，抑制不感兴趣的特征；图像分割将静脉图像转换为灰度图，其中静脉为白色，其他为黑色。

进一步的，静脉图像立体匹配模块采用基于低纹理区域的静脉图像立体匹配算法，具体为：

通过对图像进行归一化和梯度计算，检测图像中的低纹理区域，所述低纹理区域的分辨率为256*256；利用本发明提出的基于像素颜色空间的匹配代价函数，比较像素间的差异；根据像素在聚合窗口中的位置来分配不同的权值并对匹配代价值进行聚合。

所述基于低纹理区域的静脉图像立体匹配算法，具体为：

1)基于图像梯度的低纹理区域检测

采用gamma校正法对采集的静脉图像进行颜色空间的归一化，归一化公式为：

I(x,y)＝I(x,y)^gamma (1)

式中I(x,y)为输入的图像，gamma为校正参数；

对归一化后的图像进行梯度计算，分别计算出水平方向的梯度I_g,x和垂直方向的梯度I_g,y；然后建立一个尺寸为5×5的窗口，若窗口内的所有像素都满足以下条件：

{I_g,x(p)≤0.1,I_g,y(p)≤0.1|p∈W} (2)

则该窗口内的像素都标记为低纹理区域，式中W为像素p所在的窗口；

2)基于像素颜色空间的匹配代价计算

匹配代价函数为：

式中p为参考图像中的像素，q为目标图像中的待匹配像素，f(p,q)表示像素p和像素q的匹配代价值，g_p＝[g_p ^r,g_p ^g,g_p ^b]为像素p的RGB三个通道的值，g_q＝[g_q ^r,g_q ^g,g_q ^b]为像素q的RGB三个通道的值；

对颜色内相关进行相应的变化:t¹为R分量到G分量的传递函数，t²为G分量到B分量的传递函数，t³为B分量到R分量的传递函数，具体定义如下：

将像素的颜色内相关用一个三维向量t＝[t¹,t²,t³]来表示，因为颜色内相关可以进一步反应像素之间的相似性，将匹配代价计算函数修正为：

f(p,q)＝1-cos(g_p,g_q)+α₀(1-cos(t_p,t_q)) (5)

式中α₀＝0.1，使得1-cos(t_p,t_q)在匹配代价计算函数f(p,q)中所占的权重很小，削弱了1-cos(t_p,t_q)对1-cos(g_p,g_q)的影响。当两个像素点利用颜色空间RGB向量计算的结果1-cos(g_p,g_q)相差较大时，利用颜色内相关向量1-cos(t_p,t_q)计算出来的微小差异不足以影响最终的视差选择；若两个像素点利用颜色空间RGB向量计算的结果1-cos(g_p,g_q)相差很小或者相等时，利用颜色内相关向量可以进一步将两个像素区别开来；即使利用颜色空间RGB向量计算的匹配代价值出现误差，利用颜色内相关向量计算的结果也可以进行纠正，从而选择最佳的视差值。两个像素颜色越接近，cos(g_p,g_q)和cos(t_p,t_q)越大，1-cos(g_p,g_q)和1-cos(t_p,t_q)就越小，从而匹配代价值f(p,q)也越小；其中t_p和t_q分别为参考图像、目标图像的颜色内相关。

3)基于像素窗口位置的代价聚合

利用聚合窗口中像素位置的不同来进行代价聚合，窗口尺寸为7×7，窗口中间三列像素组成了中间部分，其余的区域为边缘部分，在聚合过程中对这两个部分分配不同的权值；每行的代价聚合由三部分组成：

(1)中间部分：位置与中心像素点连贯且值一样的像素点，对这部分进行取平均处理；

(2)边缘部分：位于窗口边缘，与上述中间部分不连贯，但与中心像素点值相同，对这部分进行折中处理，取平均后乘以加权系数α(0＜α＜1)；

(3)边缘部分中，与中心像素点值不同的像素点，对这一部分进行弱化处理，取平均后乘以弱化系数β(0＜β＜α＜1)。

所以该行的代价聚合值为：

式中，val_middle为上述分析中中间部分各像素的代价值之和，val_edge为边缘部分各像素的代价值之和，val_else为边缘部分中与中心像素点值不同的像素点的代价值之和，三者之和构成了窗口内总的代价值val_window，如公式(7)。N_middle为中间部分像素的个数，N_edge为边缘部分像素的个数，N_else为其他部分像素的个数，三者之和构成了窗口内像素总数N_window，如式(8)；

val_window＝val_middle+val_edge+val_else (7)

N_window＝N_middle+N_edge+N_else (8)

进一步的，多光谱静脉影像融合模块进行图像融合的方法为：

首先，对采集到的双波段静脉图像经小波变换分解为低频系数和高频系数；

然后，对高低频系数采取不同的融合准则：对于低频系数采取基于局部特征残差准则(LFRS)，通过局部特征残差求取源图像的结构信息。对于高频系数采取基于局部统计匹配特征准则(LSMS)，由局部自匹配和互匹配来确定高频融合系数。

最后，通过小波逆变换得到融合的静脉影像。

基于局部特征残差与统计匹配的多光谱静脉影像融合算法为：

1)基于局部特征残差准则的低频图像融合准则

手部静脉图像经过小波变换后得到低频信息，低频信息丰富，能量集中，反映了原图像近似信息，研究表明，人眼对图像中具有显著特征变化的区域感兴趣；局部特征可以描述图像局部信息，比如均值、熵、平均梯度、标准差等特征。局部特征值越大，表示图像局部信息越丰富。因此，局部特征可以很好地描述图像局部信息，可以通过局部特征差异来选取融合系数。

当图像局部灰度变化显著时，中心与周边区域的特征值存在差异。所以，利用中心-周边区域的对比差异进行低频融合，这样可以有效地减少亮度差异带来的影响。本节中的对比差异就是特征残差。

低频系数的对比差异的计算步骤如下：

(a)计算源图像的局部特征向量，其包含了均值l、标准差d、熵s、平均梯度g。局部特征向量表达式为：

F_in(x,y)＝[l,d,g,s] (9)

F_out(x,y)＝[l,d,g,s] (10)

其中，Fin表示3×3大小的中心区域；Fout表示7×7大小的周边区域；x与y表示空间位置。

(b)计算中心-周边区域特征残差和：

R(x,y)＝∑|F_in(x,y)-F_out(x,y)| (11)

通过上述步骤计算出图像局部特征残差和，它描述了图像局部纹理特征。残差和越大，说明中心像素处于纹理区域，残差和越小说明中心像素处于平滑区域。

利用局部残差法可以有效地提取出源图纹理结构，手部静脉皮肤纹理。

低频系数包含大量的纹理信息。为了获取纹理清晰的融合图像，将图像局部纹理特征考虑到低频融合中。因此，将上述获取的局部纹理特征用于图像融合。

低频融合准则表达式如下：

其中，

为静脉影像A的低频系数；

为静脉影像B的低频系数；w为加权系数；j为分解层数。采用一种简化的系数函数来确定加权系数，从而实现自适应融合。简化的系数函数定义如下：

式(13)中，t＝RA/RB。RA为静脉影像A的特征残差和，RB为静脉影像B的特征残差和；当RA值小于RB时，则t值小于1，w小于0.5，低频融合系数主要来自静脉影像A的低频系数；当RA值与RB值近似时，则t值趋于1，w趋于0.5，低频融合系数为静脉影像A的低频系数与静脉影像B的系数加权平均；当RA值大于RB值时，则t值大于1，w大于0.5，低频融合系数主要来自静脉影像B的低频系数；当RA值远大于RB值时，则t值较大，w趋于1，低频融合系数来自静脉影像B的系数。

2)基于局部统计匹配特征准则(LSMS)的高频图像融合准则

手部静脉图像经过小波变换后得到高频系数，高频系数反映了手部及血管边缘、轮廓等细节信息。高频系数的大小表示灰度变化程度，高频系数的符号表示灰度变化的方向。

考虑源图像的高频系数包含共有的边缘信息和独有的边缘信息，提出了一种基于局部统计匹配的高频融合准则，通过互匹配来区分不同波段手部静脉共有的边缘信息与独有的边缘信息。双波段共有的边缘信息具有高匹配度，而某一波段独有的边缘信息匹配度较低。对于共同的边缘信息，采取加权平均准则；对于独有的边缘信息，利用它的物理意义尽可能地保留它。

2.1基于改进马氏距离的互匹配

利用改进的马氏距离来计算局部区域间的空间特性差异。马氏距离定义如下：

其中，Xi和Xj为行向量；S为两行向量的协方差矩阵；S^-1为矩阵S的逆矩阵。由于矩阵S可能是奇异矩阵，不能直接计算矩阵S的逆矩阵。根据矩阵原理，可以利用奇异值分解(SVD)来计算协方差矩阵的广义逆矩阵。

协方差矩阵的奇异值分解定义为：

其中，Σ＝diag(λ₁,λ₂...λ_r)，λ_i是矩阵S的奇异值，它反映了数据能量；r是矩阵S的秩；H为共轭转置符号，U为均值矩阵，V为方差矩阵。

在实验中可发现，最大奇异值是远远大于其它的奇异值(λ_max＞＞λ_i,max≠i)，因此只保留最大奇异值，忽略其它奇异值，并将它们当作零。

因此，矩阵的广义逆矩阵定义如下：

根据以上分析，改进的马氏距离见式(18)；如果马氏距离值越小，意味着空间特征差异较小，局部区域为共有信息；

2.2高频融合

高频融合步骤分为以下：

(a)选取3×3区域，将图像局部系数按列排成1×N²矩阵。

G_j,i＝(G_j,i(1),G_j,i(2),…,G_j,i(N²)) (19)

其中，j表示分解层数；i表示方向。

(b)计算矩阵

和矩阵

的互匹配度。如果匹配度D小于阈值T，T＝0.02，则意味着源图像的局部区域是共有的边缘信息；否则，是独有的边缘信息。

(c)根据高频融合准则进行融合。

其中，H表示局部向量投影值。

本发明还提供一种基于上述多光谱三维静脉显像系统的显像方法，包括以下步骤：

在两个不同的波段上通过软件触发同时采集静脉图像；

对两幅静脉图像进行图像增强和图像分割处理，得到分割后图像；

采用基于低纹理区域的静脉图像立体匹配算法，将分割后图像进行图像立体匹配，得到静脉的着色深度图；

对采集到的双波段静脉图像经小波变换分解得到低频系数和高频系数，基于局部特征残差与统计匹配确定高频融合系数，再通过小波逆变换得到融合的静脉影像。

上述方法的四个过程与上述显像系统中相机、图像增强和分割模块、静脉图像立体匹配模块、多光谱静脉影像融合模块的具体方法相同。

本发明通过对生物医学成像研究，建立了一种多光谱三维静脉显像系统，采用两路CCD相机分别在808nm窄带和940nm窄带两个波段上通过软件触发同时采集静脉图像，获取了高质量的静脉图像；将两个波段的图像进行匹配之后，提出了一种基于低纹理区域的静脉图像立体匹配算法和一种基于局部特征残差与统计匹配的多光谱静脉影像融合算法，通过对采集到的双波段静脉图像经小波变换分解得到的低频系数和高频系数运用不同准则进行处理来确定高频融合系数，再通过小波逆变换得到融合的静脉影像。

下面结合实施例对本发明进行详细说明。

实施例

本发明提出一种多光谱三维静脉显像系统。基于低纹理区域的静脉图像立体匹配算法，获得静脉图像的三维信息；基于局部特征残差与统计匹配的多光谱静脉影像融合算法，获得静脉图像的多光谱信息。

一、多光谱三维静脉显像系统

本发明搭建的多光谱三维静脉显像系统，可以获取高质量的手背静脉深度图像，实现成像视场和投影视场的重叠；多光谱三维静脉显像系统设计流程如图1所示。

1)波段优化

本发明中的静脉图像采集装置为两台CCD相机，为了提高系统的稳定性，减少外界对相机和镜头的影响，同时简化处理过程，使用连接杆将两台相机固定在光学平板上。为了避免可见光等其他光源对成像效果的影响，本发明创新性地在相机的镜头前放置窄带滤光片，使得只能通过近红外光，经过对不同波段的成像效果分析可知940、980nm波段对皮下静脉的成像效果较好。

2)多光谱融合

针对不同成像区域，本发明中的快速切换滤光片装置可以更换不同波段的滤光片，为静脉提供了丰富的光谱维信息，选取最优波段的成像效果图作为后续融合计算模块的输入。多光谱融合得到的轮廓清晰，图像信息丰富的静脉图像如图2所示。

3)高精度静脉深度投影方法

本发明将嵌入式计算模块输出的着色深度图进行投影，可以获得明显的静脉信息，且信息量更大更加稳定。通过改变硬件布局，采用手动选取特征点的方法对成像原图和着色深度图进行了图像配准，再加上自动梯形矫正，实现了成像视场和投影视场的重叠。投影的深度图如图3所示。

二、基于低纹理区域的静脉图像立体匹配算法

针对手部静脉图像中的低纹理区域，本发明提出了一种基于像素颜色空间和窗口位置的适用于低纹理区域的立体匹配算法，算法包括四个过程：图像预处理、匹配代价计算、代价聚合、视差优化，如图4所示。

1)基于图像梯度的低纹理区域检测

通常情况下静脉的低纹理区域的像素特征不够明显，所以匹配代价计算无法准确地区分出像素之间的差异。为了改善匹配结果，检测出图像中的低纹理区域并充分利用这部分区域的纹理是很有必要的。低纹理区域检测的算法包括图像归一化、梯度计算和区域标记三个过程。

本发明采用gamma校正法对图像进行颜色空间的归一化，通过对gamma参数的调整，可以改变图像的对比度同时减少光照因素的影响，提高检测的精度。归一化公式为：

I(x,y)＝I(x,y)^gamma (1)

式中I(x,y)为输入的图像，gamma为校正参数。

对归一化后的图像进行梯度计算，分别计算出水平方向的梯度I_g,x和垂直方向的梯度I_g,y。然后建立一个尺寸为5×5的窗口，若窗口内的所有像素都满足以下条件：

{I_g,x(p)≤0.1,I_g,y(p)≤0.1|p∈W} (2)

则该窗口内的像素都标记为低纹理区域，式中W为像素p所在的窗口。

2)基于像素颜色空间的匹配代价计算

传统的基于窗口的匹配代价函数如SAD、SSD、ZSAD等算法简单高效，对噪声较为敏感；NCC中包含开平方运算，计算复杂度高。这些算法都建立在灰度的基础上，可以通过RGB三个通道扩展到彩色空间，但是会增加计算量。Le Thanh SACH中提出了一种新的匹配代价函数,采用了改进的余弦距离来衡量两个像素之间的差异性。

式中p为参考图像中的像素，q为目标图像中的待匹配像素，f(p,q)表示像素p和像素q的匹配代价值，g_p＝[g_p ^r,g_p ^g,g_p ^b]为像素p的RGB三个通道的值，g_q＝[g_q ^r,g_q ^g,g_q ^b]为像素q的RGB三个通道的值。

本发明对颜色内相关进行了相应的变化:t¹为R分量到G分量的传递函数，t²为G分量到B分量的传递函数，t³为B分量到R分量的传递函数，具体的定义如下：

将像素的颜色内相关用一个三维向量t＝[t¹,t²,t³]来表示，因为颜色内相关可以进一步反应像素之间的相似性，因此本发明将匹配代价计算函数修正为：

f(p,q)＝1-cos(g_p,g_q)+α₀(1-cos(t_p,t_q)) (5)

式中α₀＝0.1，使得1-cos(t_p,t_q)在匹配代价计算函数f(p,q)中所占的权重很小，削弱了1-cos(t_p,t_q)对1-cos(g_p,g_q)的影响。当两个像素点利用颜色空间RGB向量计算的结果1-cos(g_p,g_q)相差较大时，利用颜色内相关向量1-cos(t_p,t_q)计算出来的微小差异不足以影响最终的视差选择；若两个像素点利用颜色空间RGB向量计算的结果1-cos(g_p,g_q)相差很小或者相等时，利用颜色内相关向量可以进一步将两个像素区别开来；即使利用颜色空间RGB向量计算的匹配代价值出现误差，利用颜色内相关向量计算的结果也可以进行纠正，从而选择最佳的视差值。两个像素颜色越接近，cos(g_p,g_q)和cos(t_p,t_q)越大，1-cos(g_p,g_q)和1-cos(t_p,t_q)就越小，从而匹配代价值f(p,q)也越小。

3)基于像素窗口位置的代价聚合

代价聚合的一个关键是聚合窗口的选择，目前已经提出了许多针对图像纹理特征进行自适应窗口选择的算法，保证不丢失图像细节同时可以改善匹配精度。Jean LucLotti提出了一种利用图像的边缘作为约束条件，以窗口中心像素的四个斜角作为窗口生长的方向，图像提取的边缘为结束条件，这四个方向所包围的区域就是中心像素的自适应窗口。如果对每个像素都进行这样的操作，那么会大大增加算法的计算量，且误匹配点主要分布在左右图像中不对称的区域，所以计算出来的自适应窗口误差较大。

本发明创新性地提出利用聚合窗口中像素位置的不同来进行代价聚合，窗口尺寸为7×7，窗口内的像素分配如图5所示。

图中(x,y)为窗口中心像素，窗口中间三列像素组成了中间部分，如图中红色区域，其余的蓝色区域为边缘部分，在聚合过程中对这两个部分分配不同的权值。取窗口中任意一行的7个像素为例，对本发明的代价聚合过程进行解释，如图6所示。

图6中(a)行表示经过低纹理区域检测的图像，以此凸显边界区域的纹理特性。(b)行表示当前像素的代价值。该行的代价聚合由三部分组成：

(1)中间部分：位置与中心像素点连贯且值一样的像素点，对这部分进行取平均处理；

(2)边缘部分：位于窗口边缘，与上述中间部分不连贯，但与中心像素点值相同，对这部分进行折中处理，取平均后乘以加权系数α(0＜α＜1)；

(3)其他部分：与中心像素点值不同的像素点，对这一部分进行弱化处理，取平均后乘以弱化系数β(0＜β＜α＜1)。

所以该行的代价聚合值为：

对于整个窗口下的一般情况，可以用如下公式来表示

式中，val_middle为上述分析中中间部分各像素的代价值之和，val_edge为边缘部分各像素的代价值之和，val_else为其他部分各像素的代价值之和，三者之和构成了窗口内总的代价值val_window，如公式(8)。N_middle为中间部分像素的个数，N_edge为边缘部分像素的个数，N_else为其他部分像素的个数，三者之和构成了窗口内像素总数N_window，如式(9)。需要注意的是，边缘部分和其他部分的像素个数可能为0，因此需要提前进行判断。

val_window＝val_middle+val_edge+val_else (8)

N_window＝N_middle+N_edge+N_else (9)

三、基于局部特征残差与统计匹配的多光谱静脉影像融合算法

传统静脉成像以单波段、多视角为主，故而缺少多波段、综合化的静脉影像融合提升算法。对获取到的多波段静脉影像经小波变换后得到低频信息和高频信息，传统的低频融合准则绝大部分采取简单的线性组合，未考虑纹理特征，使得融合后的静脉手部图像细节不丰富。而常用的系数绝对值选大融合准则，虽然该准则能够将边缘信息提取出来，但同时也会将干扰信息作为高频融合系数，为静脉融合结果引入了新的噪声。针对以上问题，提出一种基于局部特征残差与统计匹配的图像融合算法以获得轮廓清晰，图像信息丰富的手部静脉融合影像。融合流程如图7所示。

1)基于局部特征残差准则(LFRS)的低频图像融合准则

手部静脉图像经过小波变换后得到低频信息，低频信息丰富，能量集中，反映了原图像近似信息人，研究表明，人眼对图像中具有显著特征变化的区域感兴趣。局部特征可以描述图像局部信息，比如均值、熵、平均梯度、标准差等特征。局部特征值越大，表示图像局部信息越丰富。因此，局部特征可以很好地描述图像局部信息，可以通过局部特征差异来选取融合系数。

当图像局部灰度变化显著时，中心与周边区域的特征值存在差异。所以，利用中心-周边区域的对比差异进行低频融合，这样可以有效地减少亮度差异带来的影响。本节中的对比差异就是特征残差。

低频系数的对比差异的计算步骤如下：

(a)计算源图像的局部特征向量，其包含了均值(l)、标准差(d)、熵(s)、平均梯度(g)。局部特征向量表达式为：

F_in(x,y)＝[l,d,g,s] (10)

F_out(x,y)＝[l,d,g,s] (11)

其中，Fin表示3×3大小的中心区域；Fout表示7×7大小的周边区域；x与y表示空间位置。

(b)计算中心-周边区域特征残差和。

R(x,y)＝∑|F_in(x,y)-F_out(x,y)| (12)

通过上述步骤计算出图像局部特征残差和，它描述了图像局部纹理特征。残差和越大，说明中心像素处于纹理区域，残差和越小说明中心像素处于平滑区域。

利用局部残差法可以有效地提取出源图纹理结构，手部静脉皮肤纹理。

低频系数包含大量的纹理信息。为了获取纹理清晰的融合图像，将图像局部纹理特征考虑到低频融合中。因此，将上述获取的局部纹理特征用于图像融合与。

低频融合准则表达式如下：

其中，L_A为静脉影像A的低频系数；L_B为静脉影像B的低频系数；w为加权系数；j为分解层数。采用一种简化的系数函数来确定加权系数，从而实现自适应融合。简化的系数函数定义如下：

式(14)中，t＝RA/RB。当RA值小于RB时，则t值小于1，w小于0.5，低频融合系数主要来自静脉影像A的低频系数；当RA值与RB值近似时，则t值趋于1，w趋于0.5，低频融合系数为静脉影像A的低频系数与静脉影像B的系数加权平均；当RA值大于RB值时，则t值大于1，w大于0.5，低频融合系数主要来自静脉影像B的低频系数；当RA值远大于RB值时，则t值较大，w趋于1，低频融合系数来自静脉影像B的系数。

2)基于局部统计匹配特征准则(LSMS)的高频图像融合准则

手部静脉图像经过小波变换后得到高频系数，高频系数反映了手部及血管边缘、轮廓等细节信息。高频系数的大小表示灰度变化程度，高频系数的符号表示灰度变化的方向。

考虑源图像的高频系数包含共有的边缘信息和独有的边缘信息，提出了一种基于局部统计匹配的高频融合准则，通过互匹配来区分不同波段手部静脉共有的边缘信息与独有的边缘信息。双波段共有的边缘信息具有高匹配度，而某一波段独有的边缘信息匹配度较低。对于共同的边缘信息，采取加权平均准则；对于独有的边缘信息，利用它的物理意义尽可能地保留它。

2.1基于改进马氏距离的互匹配

利用改进的马氏距离来计算局部区域间的空间特性差异。马氏距离定义如下：

其中，Xi和Xj为行向量；S为两行向量的协方差矩阵；S-1为矩阵S的逆矩阵。由于矩阵S可能是奇异矩阵，不能直接计算矩阵S的逆矩阵。根据矩阵原理，可以利用奇异值分解(SVD)来计算协方差矩阵的广义逆矩阵。

协方差矩阵的奇异值分解定义为：

其中，Σ＝diag(λ₁,λ₂...λ_r)，λ_i是矩阵S的奇异值，它反映了数据能量；r是矩阵S的秩；H为共轭转置符号。在实验中可发现，最大奇异值是远远大于其它的奇异值(λ_max＞＞λ_i,max≠i)，因此只保留最大奇异值，忽略其它奇异值，并将它们当作零。

因此，矩阵的广义逆矩阵定义如下：

根据以上分析，改进的马氏距离见式(19)。如果马氏距离值越小，意味着空间特征差异较小，局部区域为共有信息。

2.2基于非负矩阵分解的自匹配

将具有低匹配度的手部静脉边缘信息作为独特的边缘信息。对于同一区域的独有的边缘信息，应选择显著的边缘信息作为高频融合系数。通过理论和实验分析发现，非负矩阵分解(NMF)的物理意义可以很好地描述图像边缘信息的显著性。非负矩阵分解的详细介绍见小节。从非负矩阵分解的物理意义角度出发，矩阵分解中的加权系数就是指原矩阵中列向量在基向量上的投影。投影值越大，意味着对应列向量包含更多信息。因此，采用基于非负矩阵分解的自匹配来判断边缘信息的显著性。

2.3高频融合

高频融合步骤分为以下：

(a)选取3×3区域，将图像局部系数按列排成1×N2矩阵。

G_j,i＝(G_j,i(1),G_j,i(2),…,G_j,i(N²)) (20)

其中，j表示分解层数；i表示方向。

(b)计算矩阵

和矩阵

的互匹配度。如果匹配度(D)小于阈值(T＝0.02)，则意味着源图像的局部区域是共有的边缘信息。否则，是独有的边缘信息。

(c)根据高频融合准则进行融合。

其中，H表示局部向量投影值。

本发明的具体实现步骤如下：

步骤1,建立多光谱三维静脉显像系统，静脉成像过程中，被照物反射led光源的两束光，分别通过808nm窄带滤光片和940nm窄带滤光片的滤光进入相机1和相机2中，同时在第一帧图像的获取时刻采用软件触发，同时控制两个相机进行图像采集，获得了两个波段的高质量的清晰静脉图像。在满足左相机与右相机等高的条件下，将左相机与右相机的极线距离拉大至30cm，使显像分辨精度得到了显著提升。利用上述多光谱三维静脉显像系统采集到的静脉图像，经过旋转平移的极线校正操作之后，作为立体匹配的输入。把808nm窄带滤光后的图像作为基准图像对另外一路图像进行极线校正。

步骤2，将经过808nm窄带滤光后的图像进行图像增强和图像分割处理之后，得到808nm分割后图像。

步骤3，将经过940nm窄带滤光后的图像进行图像增强和图像分割处理之后，得到940nm分割后图像。

步骤4，将步骤2和步骤3分别得到的808nm分割后图像与940nm分割后图像进行图像立体匹配，采用基于像素颜色空间和窗口位置的适用于低纹理区域的立体匹配算法，首先通过对图像进行归一化和梯度计算，检测图像中的低纹理区域；其次利用本发明提出的基于像素颜色空间的匹配代价函数，可以比较出像素间的微小差异；最后根据像素在聚合窗口中的位置来分配不同的权值并对匹配代价值进行聚合，最终得到静脉的着色深度图。

步骤5，通过基于局部特征残差与统计匹配的多光谱静脉影像融合算法进行图像融合：首先对采集得到的808nm和990nm波段的静脉图像经小波变换分解为低频系数和高频系数；然后，对于低频系数采取基于局部特征残差准则(LFRS)，通过局部特征残差求取源图像的结构信息；对于高频系数采取基于局部统计匹配特征准则(LSMS)，由局部自匹配和互匹配来确定高频融合系数；最后，通过小波逆变换得到轮廓清晰，图像信息丰富的静脉融合图像。

Claims (8)

1.一种多光谱三维静脉显像系统，其特征在于，包括：

两路CCD相机，用于在两个不同的波段上通过软件触发同时采集静脉图像；

图像增强和分割模块，用于对两幅静脉图像进行图像增强和图像分割处理，得到分割后图像；

静脉图像立体匹配模块，用于将分割后图像进行图像立体匹配，得到静脉的着色深度图；采用基于低纹理区域的静脉图像立体匹配算法，通过对图像进行归一化和梯度计算，检测图像中的低纹理区域，利用基于像素颜色空间的匹配代价函数，比较像素间的差异；根据像素在聚合窗口中的位置来分配不同的权值并对匹配代价值进行聚合；具体为：

1)基于图像梯度的低纹理区域检测

采用gamma校正法对采集的静脉图像进行颜色空间的归一化，归一化公式为：

I(x,y)＝I(x,y)^gamma (1)

式中I(x,y)为输入的图像，gamma为校正参数；

对归一化后的图像进行梯度计算，分别计算出水平方向的梯度I_g,x和垂直方向的梯度I_g,y；然后建立一个尺寸为5×5的窗口，若窗口内的所有像素都满足以下条件：

{I_g,x(p)≤0.1,I_g,y(p)≤0.1|p∈W} (2)

则该窗口内的像素都标记为低纹理区域，式中W为像素p所在的窗口；

2)基于像素颜色空间的匹配代价计算

匹配代价函数为：

式中p为参考图像中的像素，q为目标图像中的待匹配像素，f(p,q)表示像素p和像素q的匹配代价值，g_p＝[g_p ^r,g_p ^g,g_p ^b]为像素p的RGB三个通道的值，g_q＝[g_q ^r,g_q ^g,g_q ^b]为像素q的RGB三个通道的值；

对颜色内相关进行相应的变化:t¹为R分量到G分量的传递函数，t²为G分量到B分量的传递函数，t³为B分量到R分量的传递函数，具体定义如下：

将像素的颜色内相关用一个三维向量t＝[t¹,t²,t³]来表示，将匹配代价计算函数修正为：

f(p,q)＝1-cos(g_p,g_q)+α₀(1-cos(t_p,t_q)) (5)

3)基于像素窗口位置的代价聚合

利用聚合窗口中像素位置的不同来进行代价聚合，窗口尺寸为7×7，窗口中间三列像素组成了中间部分，其余的区域为边缘部分，在聚合过程中对这两个部分分配不同的权值；每行的代价聚合由三部分组成：

(1)中间部分：位置与中心像素点连贯且值一样的像素点，对这部分进行取平均处理；

(2)边缘部分中，位于窗口边缘，与上述中间部分不连贯，但与中心像素点值相同，对这部分进行折中处理，取平均后乘以加权系数α，0＜α＜1；

(3)边缘部分中，与中心像素点值不同的像素点，对这一部分进行弱化处理，取平均后乘以弱化系数β，0＜β＜α＜1；

该行的代价聚合值为：

式中，val_middle为上述分析中中间部分各像素的代价值之和，val_edge为边缘部分各像素的代价值之和，val_else为边缘部分中与中心像素点值不同的像素点的代价值之和，三者之和构成了窗口内总的代价值val_window，如公式(7)；N_middle为中间部分像素的个数，N_edge为边缘部分像素的个数，N_else为其他部分像素的个数，三者之和构成了窗口内像素总数N_window，如式(8)；

val_window＝val_middle+val_edge+val_else (7)

N_window＝N_middle+N_edge+N_else (8)

多光谱静脉影像融合模块，对采集到的双波段静脉图像经小波变换分解得到低频系数和高频系数，基于局部特征残差与统计匹配确定高频融合系数，再通过小波逆变换得到融合的静脉影像。

2.根据权利要求1所述的多光谱三维静脉显像系统，其特征在于，静脉成像过程中，被照物反射LED光源的两束光，分别通过808nm窄带滤光片和940nm窄带滤光片的滤光进入两路CCD相机，在第一帧图像的获取时刻采用软件触发，同时控制两个相机进行图像采集，获得两个波段的静脉图像；808nm窄带滤光片中心波长为808nm，半带宽为30nm，940nm窄带滤光片中心波长为940nm，半带宽为30nm。

3.根据权利要求2所述的多光谱三维静脉显像系统，其特征在于，两个相机通过连接杆固定在光学平板上，两个CCD相机等高，极线距离为20～40cm。

4.根据权利要求1所述的多光谱三维静脉显像系统，其特征在于，多光谱静脉影像融合模块采用基于局部特征残差与统计匹配的多光谱静脉影像融合算法进行图像融合，首先，对采集到的双波段静脉图像经小波变换分解为低频系数和高频系数；然后，对高低频系数采取不同的融合准则：对于低频系数采取基于局部特征残差准则，通过局部特征残差求取源图像的结构信息；对于高频系数采取基于局部统计匹配特征准则，由局部自匹配和互匹配来确定高频融合系数；最后，通过小波逆变换得到融合的静脉影像；

所述基于局部特征残差与统计匹配的多光谱静脉影像融合算法具体为：

1)基于局部特征残差准则的低频图像融合准则

(a)计算源图像的局部特征向量，包含均值l、标准差d、熵s、平均梯度g，局部特征向量表达式为：

F_in(x,y)＝[l,d,g,s] (9)

F_out(x,y)＝[l,d,g,s] (10)

其中，Fin表示3×3大小的中心区域，Fout表示7×7大小的周边区域，x与y表示空间位置；

(b)计算中心-周边区域特征残差和：

R(x,y)＝∑|F_in(x,y)-F_out(x,y)| (11)

低频融合准则表达式如下：

其中，

为静脉影像A的低频系数；

为静脉影像B的低频系数；w为加权系数，j为分解层数；采用简化的系数函数来确定加权系数，简化的系数函数定义如下：

式(13)中，t＝RA/RB，RA为静脉影像A的特征残差和，RB为静脉影像B的特征残差和；

2)基于局部统计匹配特征准则的高频图像融合准则

2.1基于改进马氏距离的互匹配

利用改进的马氏距离来计算局部区域间的空间特性差异，马氏距离定义如下：

其中，X_i和X_j为行向量，S为两行向量的协方差矩阵，S^-1为矩阵S的逆矩阵；利用奇异值分解计算协方差矩阵的广义逆矩阵；

协方差矩阵的奇异值分解定义为：

其中，Σ＝diag(λ₁,λ₂...λ_r)，λ_i是矩阵S的奇异值，r是矩阵S的秩，H为共轭转置符号，U为均值矩阵，V为方差矩阵；

最大奇异值远远大于其它的奇异值，λ_max＞＞λ_i,max≠i，因此只保留最大奇异值，忽略其它奇异值，并将它们当作零，有：

因此，矩阵的广义逆矩阵定义如下：

根据以上分析，改进的马氏距离为：

2.2高频融合

(a)选取3×3区域，将图像局部系数按列排成1×N²矩阵；

G_j,i＝(G_j,i(1),G_j,i(2),…,G_j,i(N²)) (19)

其中，j表示分解层数，i表示方向；

(b)计算矩阵

和矩阵

的互匹配度；如果匹配度D小于阈值，则源图像的局部区域是共有的边缘信息，否则，是独有的边缘信息；

(c)根据高频融合准则进行融合：

其中，K表示局部向量投影值。

5.一种基于权利要求1所述多光谱三维静脉显像系统的显像方法，其特征在于，包括以下步骤：

在两个不同的波段上通过软件触发同时采集静脉图像；

对两幅静脉图像进行图像增强和图像分割处理，得到分割后图像；

采用基于低纹理区域的静脉图像立体匹配算法，将分割后图像进行图像立体匹配，得到静脉的着色深度图，具体方法为：

1)基于图像梯度的低纹理区域检测

采用gamma校正法对采集的静脉图像进行颜色空间的归一化，归一化公式为：

I(x,y)＝I(x,y)^gamma (1)

式中I(x,y)为输入的图像，gamma为校正参数；

对归一化后的图像进行梯度计算，分别计算出水平方向的梯度I_g,x和垂直方向的梯度I_g,y；然后建立一个尺寸为5×5的窗口，若窗口内的所有像素都满足以下条件：

{I_g,x(p)≤0.1,I_g,y(p)≤0.1|p∈W} (2)

则该窗口内的像素都标记为低纹理区域，式中W为像素p所在的窗口；

2)基于像素颜色空间的匹配代价计算

匹配代价函数为：

式中p为参考图像中的像素，q为目标图像中的待匹配像素，f(p,q)表示像素p和像素q的匹配代价值，g_p＝[g_p ^r,g_p ^g,g_p ^b]为像素p的RGB三个通道的值，g_q＝[g_q ^r,g_q ^g,g_q ^b]为像素q的RGB三个通道的值；

对颜色内相关进行相应的变化:t¹为R分量到G分量的传递函数，t²为G分量到B分量的传递函数，t³为B分量到R分量的传递函数，具体定义如下：

将像素的颜色内相关用一个三维向量t＝[t¹,t²,t³]来表示，将匹配代价计算函数修正为：

f(p,q)＝1-cos(g_p,g_q)+α₀(1-cos(t_p,t_q)) (5)

3)基于像素窗口位置的代价聚合

利用聚合窗口中像素位置的不同来进行代价聚合，窗口尺寸为7×7，窗口中间三列像素组成了中间部分，其余的区域为边缘部分，在聚合过程中对这两个部分分配不同的权值；每行的代价聚合由三部分组成：

(1)中间部分：位置与中心像素点连贯且值一样的像素点，对这部分进行取平均处理；

(2)边缘部分：位于窗口边缘，与上述中间部分不连贯，但与中心像素点值相同，对这部分进行折中处理，取平均后乘以加权系数α，0＜α＜1；

(3)边缘部分中，与中心像素点值不同的像素点，对这一部分进行弱化处理，取平均后乘以弱化系数β，0＜β＜α＜1；

该行的代价聚合值为：

式中，val_middle为上述分析中中间部分各像素的代价值之和，val_edge为边缘部分各像素的代价值之和，val_else为边缘部分中与中心像素点值不同的像素点的代价值之和，三者之和构成了窗口内总的代价值val_window，如公式(7)；N_middle为中间部分像素的个数，N_edge为边缘部分像素的个数，N_else为其他部分像素的个数，三者之和构成了窗口内像素总数N_window，如式(8)；

val_window＝val_middle+val_edge+val_else (7)

N_window＝N_middle+N_edge+N_else (8)

对采集到的双波段静脉图像经小波变换分解得到低频系数和高频系数，基于局部特征残差与统计匹配确定高频融合系数，再通过小波逆变换得到融合的静脉影像。

6.根据权利要求5所述的方法，其特征在于，被照物反射LED光源的两束光，分别通过808nm窄带滤光片和940nm窄带滤光片的滤光进入两路CCD相机，在第一帧图像的获取时刻采用软件触发，同时控制两个相机进行图像采集，获得两个波段的静脉图像；808nm窄带滤光片中心波长为808nm，半带宽为30nm，940nm窄带滤光片中心波长为940nm，半带宽为30nm。

7.根据权利要求6所述的方法，其特征在于，两个相机通过连接杆固定在光学平板上，两个CCD相机等高，极线距离为20～40cm。

8.根据权利要求5所述的方法，其特征在于，采用基于局部特征残差与统计匹配的多光谱静脉影像融合算法进行图像融合，首先，对采集到的双波段静脉图像经小波变换分解为低频系数和高频系数；然后，对高低频系数采取不同的融合准则：对于低频系数采取基于局部特征残差准则，通过局部特征残差求取源图像的结构信息；对于高频系数采取基于局部统计匹配特征准则，由局部自匹配和互匹配来确定高频融合系数；最后，通过小波逆变换得到融合的静脉影像；

基于局部特征残差与统计匹配的多光谱静脉影像融合算法具体为：

1)基于局部特征残差准则的低频图像融合准则

(a)计算源图像的局部特征向量，包含均值l、标准差d、熵s、平均梯度g，局部特征向量表达式为：

F_in(x,y)＝[l,d,g,s] (9)

F_out(x,y)＝[l,d,g,s] (10)

其中，Fin表示3×3大小的中心区域，Fout表示7×7大小的周边区域，x与y表示空间位置；

(b)计算中心-周边区域特征残差和：

R(x,y)＝∑|F_in(x,y)-F_out(x,y)| (11)

低频融合准则表达式如下：

其中，

为静脉影像A的低频系数；

为静脉影像B的低频系数；w为加权系数，j为分解层数；采用简化的系数函数来确定加权系数，简化的系数函数定义如下：

式(13)中，t＝RA/RB，RA为静脉影像A的特征残差和，RB为静脉影像B的特征残差和；

2)基于局部统计匹配特征准则的高频图像融合准则

2.1基于改进马氏距离的互匹配

利用改进的马氏距离来计算局部区域间的空间特性差异，马氏距离定义如下：

其中，X_i和X_j为行向量，S为两行向量的协方差矩阵，S^-1为矩阵S的逆矩阵；利用奇异值分解计算协方差矩阵的广义逆矩阵；

协方差矩阵的奇异值分解定义为：

其中，Σ＝diag(λ₁,λ₂...λ_r)，λ_i是矩阵S的奇异值，r是矩阵S的秩，H为共轭转置符号，U为均值矩阵，V为方差矩阵；

最大奇异值远远大于其它的奇异值，λ_max ＞＞λ_i,max≠i，因此只保留最大奇异值，忽略其它奇异值，并将它们当作零，有：

因此，矩阵的广义逆矩阵定义如下：

根据以上分析，改进的马氏距离为：

2.2高频融合

(a)选取3×3区域，将图像局部系数按列排成1×N²矩阵；

G_j,i＝(G_j,i(1),G_j,i(2),…,G_j,i(N²)) (19)

其中，j表示分解层数，i表示方向；

(b)计算矩阵

和矩阵

的互匹配度；如果匹配度D小于阈值，则源图像的局部区域是共有的边缘信息，否则，是独有的边缘信息；

(c)根据高频融合准则进行融合：

其中，K表示局部向量投影值。

Images