1. 引言
皮下静脉穿刺是常见的医疗手段之一,在静脉采血、输血、注射药物之前都需要对患者进行皮下静脉的穿刺。但是,许多因素会影响到皮下静脉穿刺的成功率,如儿童静脉纤细、深色皮肤等。手背静脉显像技术可以帮助医护人员找到血管,提高皮下静脉穿刺的成功率。主流的手背静脉成像技术基于近红外成像技术,利用了静脉血管中的血红蛋白对于700 nm~900 nm的近红外光的吸收率高于水、脂肪等周围的皮肤组织的特性 [1] 。由于光强、光照角度、手背皮肤散射、个体手背厚薄等因素,静脉血管和周围的皮肤组织之间的分界不够明显,图像的对比度较低且噪声较多。噪声会影响后续的血管分割提取,导致静脉血管细节的模糊或缺失,因此需要在保证图像中静脉血管的细节特征不受影响下对图像中的噪声进行滤除。现如今,双边滤波由于其能够为图像去噪的同时保留图像的边缘信息而被越来越多地使用,该方法非常适合于手背静脉图像的去噪处理。考虑到手背静脉显像的实时性要求较高,双边滤波用软件实现比较慢,因此采用硬件实现能达到更快的处理速度和显示效果。随着近年来芯片工艺的发展,FPGA的性能不断提高,使用成本大幅降低。用户可以根据需求使用软件设计出硬件逻辑,并利用FPGA中集成的资源帮助图像采集及运算,降低了开发难度。利用FPGA的并行以及流水线处理的特点,可以有效地提高双边滤波算法的计算速度 [2] 。本文利用双边滤波的优点,在FPGA上实现了双边滤波的算法,能够有效地对手背静脉图像进行去噪,提高手背静脉图像的清晰度。
2. 双边滤波算法
双边滤波是一种结合了空间邻近关系和亮度相似度的滤波器。空间邻近度滤波器通过平均滤波器窗口中心点邻域的像素值,起到低通滤波器的作用,对图像进行平滑。而亮度相似度滤波器只允许平均相似的像素值,假如滤波器窗口内的像素与中心像素的像素值相差较大,该像素则会被忽略,起到保持边缘的作用。因此,双边滤波器既平滑了图像,又保持了图像的边缘。
双边滤波由下式定义:
(1)
定义m = (x, y)代表像素的坐标,m0 = (x0, y0)和
分别代表滤波前后中心像素的坐标。因此,式中
为滤波后中心像素的像素值,
和
分别为滤波窗口F中邻域像素和中心像素的像素值。
(2)
(3)
式(2)和(3)分别称之为亮度相似度算子和空域邻近度算子 [3] 。式中,σph和σc定义了高斯曲线的宽度,决定了双边滤波器的性能 [4] 。亮度相似度算子比较每个邻域像素与中心像素的灰度值,基于σph因子计算出相应权重系数。当像素值差值的绝对值越大,相应的权重系数越小,反之亦然。空域邻近度算子是一个低通滤波器,其权重与中心像素与邻域像素的空间距离成反比。两种权重相结合,每个像素由其邻近点的加权平均代替。
(4)
式(4)对加权后的像素值进行归一化,保证滤波后的图像灰度级保持在原来的范围。
3. 双边滤波算法的FPGA实现
3.1. 双边滤波处理步骤
双边滤波算法实现的总体框图如图1所示,图像逐个像素从SDRAM中读入。为了形成N × N大小的窗口,输入的像素需要先经过N − 1个图像宽度的行缓存,再逐列输出到N × N的窗口寄存器中。窗口寄存器中的每个像素点分别输入到各自的亮度相似度算子,并行计算出亮度相似度权重系数,并与已知的空域邻近度权重系数相乘。计算出权重系数后,每个像素与相应权重系数相乘,然后分别对加权像素值与权重系数求和,二者相除得到双边滤波后的图像。
Figure 1. Flow chart of algorithm process
图1. 算法处理步骤流程图
3.2. 亮度相似度权重系数计算
根据式(2),计算亮度相似度权重系数需要复杂的指数运算。在FPGA中实现指数运算的常用函数有CORDIC算法、分段线性逼近法等,然而其计算范围小、精度差,而且需要耗费大量硬件资源和运算周期 [5] 。而查表法是根据精度要求计算出结果,并将结果预先写入查找表(LUT)中,然后在计算时通过查表在1个时钟周期可以直接得到结果。因此使用查表法查找权重系数预先的计算结果代替权重系数的直接计算,可以有效地减少了资源的消耗和计算的时间。
根据式(2),两个像素值的差值可以作为查找表中的相应权重系数的查找地址,且权重系数的数量是有限的,并受到下列三个参数的影响:
1) 像素值的位宽N;
2) 参数σph;
3) 权重系数的位宽W;
提升图像的像素值位宽,会导致像素值差值的数量增多,预先存入查找表中的权重系数的数量会增加。参数σph会决定高斯曲线的陡峭程度,与权重系数的位宽共同影响计算精度。根据指数的性质,亮度相似度权重与空间邻近度权重的取值范围为
与
,因此使用定点数代替浮点数来表示权重系数。当像素值差值越大,负指数向零收敛,此时只有有限数量的权重系数大于零。在权重系数的位宽W下,对权重系数采取四舍五入,并将最后一个大于0的权重系数所对应的像素值差值定为极限值。假设像素值位宽N = 8,像素值差值的可取0~255,权重系数总共有256个值,查找表深度为256,权重系数位宽W = 8,消耗256 × 8 bit的存储器资源。如图2所示,此时像素值差值90为极限值,只有90个权重系数会被存储到查找表,消耗的存储器资源随之减少。
Figure 2. Limitation of the number of weight coefficient
图2. 数量有限的权重系数
亮度相似度权重系数计算步骤如图3所示,窗口寄存器中的邻域像素
中心像素
相减,得到像素值差值来查找基于特定σph计算出的亮度相似度权重系数表,得到像素m对应的亮度相似度权重系数p(m)。若像素值差值超过极限值,将跳过查表,输出结果0。
Figure 3. Brightness similarity operato
图3. 亮度相似度算子
3.3. 空间邻近度权重系数计算
根据式(3),计算空间邻近度权重系数与亮度相似度权重系数相似,需要进行复杂的指数运算。空间邻近度权重系数对于特定大小的窗口是恒定的,而且在水平和垂直方向上对称,可预先计算好结果。如图4所示,对于3 × 3大小的窗口,只需计算空间距离为0、1、2的权重系数,根据像素的空间位置直接查表可得相应的空间邻近度权重系数g(m),并与上一步的亮度相似度权重系数相乘得到每个像素的综合权重系数w(m)。
3.4. 归一化
根据式(1)和式(4),加权后的像素值和权重系数分别需要进行累加,二者相除进行归一化。如图5与图6所示,采用累加树的形式对加权像素值和权重系数分别累加,对于3 × 3大小的窗口,需要消耗4个时钟周期得到累加结果,最后二者相除得到双边滤波后的像素值。
Figure 6. Normalization of weight pixel value
图6. 加权像素值的归一化
4. 实验结果
采用Intel Cyclone IV E系列的EP4CE22F17C8芯片实现算法,在752 × 480的图像分辨率,及5 × 5窗口大小的参数下,双边滤波模块总共消耗了3038个逻辑寄存器,79440位RAM存储资源,及72个乘法器资源。使用Matlab对8位的752 × 480分辨率大小的手背静脉图像添加σnoise = 0.5的高斯噪声,用ModelSim对图像进行双边滤波仿真,双边滤波器中参数分别为:σph = 25,σc = 25,W = 8,另使用相同参数的高斯滤波对噪声图像进行处理作为比较。实验结果如图7所示,经过双边滤波后图像得到平滑,双边滤波的去噪效果良好,且相较于高斯滤波这一常见的图像去噪方法,保持了手背静脉血管的结构和边缘。由于流水线设计的特点,在25 MHz的时钟频率下,FPGA计算一帧图像需要耗费14 ms,每秒能计算超过60帧图像,保证了系统的实时性。
Figure 7. Images filtered by bilateral filter and gauss filter
图7. 双边滤波与高斯滤波处理后图像
5. 结论
本文采用双边滤波算法对手背静脉图像进行滤波,与传统的高斯滤波方法相比,能滤除图像中噪声的同时保持了血管的结构和边缘,为后续的静脉血管分割和提取奠定良好基础。在FPGA上实现的双边滤波算法,利用FPGA并行和流水线计算的优点,提升了算法运算速度,满足了手背静脉显像系统实时性的需求,并提高了图像的清晰度。该算法能广泛应用于图像处理领域,具有良好的应用前景。
不足之处在于,目前本文在FPGA上实现的双边滤波算法不能动态地调节参数,只能作算法的验证,未来需要加入窗口大小以及权重系数结果进行动态调整。
致谢
非常感谢2018国家级大学生创新创业计划资助项目(2018101450201)对本文的资助。
NOTES
*通讯作者。