CN113011333B - 基于近红外图像获取最佳静脉穿刺点和方向的系统及方法 - Google Patents

基于近红外图像获取最佳静脉穿刺点和方向的系统及方法 Download PDF

Info

Publication number
CN113011333B
CN113011333B CN202110295569.5A CN202110295569A CN113011333B CN 113011333 B CN113011333 B CN 113011333B CN 202110295569 A CN202110295569 A CN 202110295569A CN 113011333 B CN113011333 B CN 113011333B
Authority
CN
China
Prior art keywords
vein
image
point
value
pixel
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110295569.5A
Other languages
English (en)
Other versions
CN113011333A (zh
Inventor
赵悦
宋健
吴昊
刘健行
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN202110295569.5A priority Critical patent/CN113011333B/zh
Publication of CN113011333A publication Critical patent/CN113011333A/zh
Application granted granted Critical
Publication of CN113011333B publication Critical patent/CN113011333B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V40/00—Recognition of biometric, human-related or animal-related patterns in image or video data
    • G06V40/10—Human or animal bodies, e.g. vehicle occupants or pedestrians; Body parts, e.g. hands
    • G06V40/12—Fingerprints or palmprints
    • G06V40/1347—Preprocessing; Feature extraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/20—Image enhancement or restoration by the use of local operators
    • G06T5/70—
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/13—Edge detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/20—Image preprocessing
    • G06V10/25—Determination of region of interest [ROI] or a volume of interest [VOI]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/44—Local feature extraction by analysis of parts of the pattern, e.g. by detecting edges, contours, loops, corners, strokes or intersections; Connectivity analysis, e.g. of connected components
    • G06V10/443—Local feature extraction by analysis of parts of the pattern, e.g. by detecting edges, contours, loops, corners, strokes or intersections; Connectivity analysis, e.g. of connected components by matching or filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V40/00—Recognition of biometric, human-related or animal-related patterns in image or video data
    • G06V40/10—Human or animal bodies, e.g. vehicle occupants or pedestrians; Body parts, e.g. hands
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10048—Infrared image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20024—Filtering details
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20024—Filtering details
    • G06T2207/20028—Bilateral filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20092—Interactive image processing based on input by user
    • G06T2207/20104—Interactive definition of region of interest [ROI]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30101—Blood vessel; Artery; Vein; Vascular
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V40/00—Recognition of biometric, human-related or animal-related patterns in image or video data
    • G06V40/10—Human or animal bodies, e.g. vehicle occupants or pedestrians; Body parts, e.g. hands
    • G06V40/14—Vascular patterns

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Multimedia (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Human Computer Interaction (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

基于近红外图像获取最佳静脉穿刺点和方向的系统及方法，它属于人体手背静脉血管识别技术领域。本发明解决了采用现有方法难以找出最佳静脉穿刺点及穿刺方向的问题。本发明采用双边滤波方法对采集的手背近红外图像进行降噪，以尽可能保护静脉边缘信息的同时减少噪声，并依据手背轮廓特征动态选取ROI感兴趣区域，减少计算量以提高运算效率；设计基于Hessian矩阵的静脉特征增强算法，利用二阶微分算子与ROI图像进行微分以求取Hessian矩阵，借助其特征值和特征向量构造静脉增强滤波函数，获取静脉局部特征，提取出完整静脉网络结构，提高静脉位置识别的精度；开发静脉宽度特征可视化算法，筛选最佳静脉穿刺点以及穿刺方位。本发明可以应用于手背静脉血管识别。

Description

基于近红外图像获取最佳静脉穿刺点和方向的系统及方法
技术领域
本发明属于人体手背静脉血管识别技术领域，具体涉及一种基于近红外图像获取最佳静脉穿刺点和方向的系统及方法。
背景技术
皮下静脉穿刺已经成为目前临床治疗中最为普遍的医疗手段，包括静脉输液、采血以及输血等。静脉穿刺的成功与否直接影响到患者的治疗、抢救是否及时、有效。因此，对于医护人员来说，能够准确且熟练地为患者进行手背静脉穿刺是非常重要的。然而现实生活中，由于患者个体皮肤颜色、静脉深浅以及脂肪厚度不同等因素的影响，经常会造成静脉穿刺失误。统计数据表明：普通成人的静脉穿刺首次失败率为28％；儿童的静脉穿刺首次失败率高达44％，其中需尝试穿刺3次以上才能成功进行注射的比例为43％。另一方面，从医护人员的角度来说，静脉访问成功率是评价其职业技术水平的标准之一。在临床操作中，如果未能及时、准确地访问患者静脉，很容易会使医护人员产生沮丧、急躁、缺乏自信等不利于身心健康以及工作效率的负面情绪。这些心理压力会影响医护人员的进一步发挥，从而影响到医院整体的医疗水平。
随着人工智能技术的发展，智能采血系统得到了越来越广的关注。智能采血系统主要依赖图像引导采血机器人实现。图像引导静脉穿刺的首要条件是能够采集到较为清晰的静脉血管图像并识别出静脉血管的位置。现有的静脉访问的辅助设备通常以近红外成像技术为基础，采集静脉图像，而后通过图像特征增强来处理改善图像中静脉血管模糊、背景不均匀的特性。目前，市场上在售的静脉显示仪主要有两种：一种是利用人体肌肉组织和静脉血管对可见红光的透过吸收率不同来提高护士目视识别静脉走向，分辨能力较强，操作简单。另一种是近红外静脉显像仪：其将近红外光投影至待检测部位，用红外摄像模块来接收而后成像在显示器上，医护人员通过观察显示器来获取静脉信息。但是由于实际采集到的近红外图像中噪声干扰较多，对比度偏低，图像清晰度不够，而这些因素又都会对采集图像中的静脉血管位置识别精度产生影响，因此，采用现有方法对静脉血管位置识别的精度仍然较低，进而难以找出最佳的静脉穿刺点及穿刺方向，通过配合图像增强算法改进图像清晰度并去除图像中噪声，便于医护人员进行静脉访问，对于寻找最佳静脉穿刺点及穿刺方向是十分必要的。
发明内容
本发明的目的是为解决由于采用现有方法对静脉血管位置识别的精度低，导致难以找出最佳的静脉穿刺点及穿刺方向的问题，而提出了一种基于近红外图像获取最佳静脉穿刺点和方向的系统及方法。
本发明为解决上述技术问题所采取的技术方案是：
基于本发明的一个方面，一种基于近红外图像获取最佳静脉穿刺点和方向的系统，所述系统包括近红外图像采集模块、图像预处理模块、ROI区域提取模块、特征增强模块、穿刺点及穿刺方向定位模块、图像配准模块；其中：
所述近红外图像采集模块用于采集人体手背的近红外图像，并将采集的图像发送给图像预处理模块；
所述图像预处理模块用于对接收的图像进行灰度化和滤波降噪处理，获得预处理后的图像，并将预处理后的图像发送给ROI区域提取模块；
所述ROI区域提取模块用于提取和拟合手背图像轮廓，根据提取和拟合的手背图像轮廓定位ROI区域位置；并将ROI区域的图像发送给特征增强模块；
所述特征增强模块用于对ROI区域图像中的静脉进行特征增强，获得特征增强后的图像，并对特征增强后的图像进行灰度反转，获得灰度反转后的图像；
所述穿刺点及穿刺方向定位模块用于提取灰度反转后图像的静脉边缘，根据提取的静脉边缘检测静脉横向宽度和静脉骨架位置，并定位穿刺点及穿刺方向；
所述图像配准模块用于将近红外图像中ROI区域的静脉位置信息配准到自然光学图像中，得到在自然光学图像中的最佳静脉穿刺点以及穿刺方向。
基于本发明的另一个方面，一种基于近红外图像获取最佳静脉穿刺点和方向的方法，所述方法具体包括以下步骤：
步骤一、采用近红外静脉采集器采集人体手背的近红外图像，并对采集的图像进行图像灰度化处理，获得灰度化处理后的图像；再使用双边滤波器对灰度化处理后的图像进行滤波降噪，获得滤波降噪后的图像；
从滤波降噪后的图像中提取出手背图像轮廓，并计算出手背图像轮廓的矩特征，根据矩特征计算出手背图像轮廓的中心；
对手背图像轮廓进行拟合，计算出手背图像轮廓的宽度width和下边界点纵坐标ydown；并根据手背图像轮廓的中心的横坐标、宽度width以及下边界点纵坐标ydown选取ROI区域；
步骤二、提取ROI区域内的静脉结构，获取静脉特征增强图像；并对获得的静脉特征增强图像进行灰度反转，得到灰度反转后的图像；
步骤三、利用Canny边缘检测方法提取出灰度反转后图像中的静脉边缘，再根据提取出的静脉边缘进行静脉横向宽度和静脉骨架位置检测；
根据检测出的静脉骨架位置计算方向修正角度，利用方向修正角度对静脉横向宽度进行修正，获得修正后的静脉横向宽度；
根据修正后的静脉横向宽度进行最佳静脉穿刺点以及穿刺方向定位；
步骤四、利用步骤二的方法对人体手背的自然光学图像进行ROI区域的提取，再将近红外图像ROI区域中的静脉位置信息配准到自然光学图像中，得到在自然光学图像中的最佳静脉穿刺点以及穿刺方向。
本发明的有益效果是：本发明提出了一种基于近红外图像获取最佳静脉穿刺点和方向的系统及方法，本发明采集手背近红外图像，采用双边滤波的方法进行降噪，在尽可能保护静脉边缘信息的前提下减少噪声，而后依据手背轮廓特征，动态选取ROI感兴趣区域，减少计算量，以提高运算效率；设计基于Hessian矩阵的静脉特征增强算法，利用二阶微分算子与ROI图像进行微分以求取Hessian矩阵，借助其特征值和特征向量构造静脉增强滤波函数，获取静脉局部特征，提取出完整的静脉网络结构，提高了对静脉位置识别的精度；开发静脉宽度特征可视化算法，由此筛选出最佳静脉穿刺点以及穿刺方位，并进行穿刺安全性检验；将上述获得的近红外图像中的静脉信息配准到自然光学图片中，增强可观性，为医护人员提供可视化的信息辅助指导。
附图说明
图1为本发明方法的流程图；
图2为二维图像线状结构的Hessian矩阵特征值分布图；
图3a)为血管模型图；
图3b)为血管横切面灰度值分布图；
图4为静脉方向修正示意图；
图5为手背近红外静脉图像；
图6为手背静脉轮廓图；
图7为手背ROI区域图像；
图8为切割ROI图像；
图9为反转后多尺度静脉增强图像；
图10为静脉边缘提取效果图；
图11a)为静脉横向宽度示意图图1；
图11b)为对图11a)进行宽度修正后的实际静脉宽度；
图11c)为静脉横向宽度示意图图2；
图11d)为对图11c)进行宽度修正后的实际静脉宽度；
图中标注：白点为静脉骨架位置，数据依次为静脉宽度、静脉方向(弧度制)，静脉方向取向下为初始方向，顺时针为正方向；
图12为穿刺方案示意图；
图中标注：白线为筛选出的适于穿刺的静脉片段，白点为推荐静脉访问位置，数据依次为静脉访问推荐指数排名、推荐静脉穿刺方向(弧度制)，静脉方向取向上为初始方向，逆时针为正方向；
图13a)为穿刺方案在近红外图像下的示意图；
图13b)为穿刺方案在自然光学图像下的示意图；
图14a)为GUI进程——ROI区域的示意图；
图14b)为GUI进程——宽度可视化的示意图；
图14c)为GUI进程——穿刺点及穿刺方向定位算法的示意图；
图14d)为GUI进程——近红外图像与自然光学图像配准的示意图；
图15为静脉宽度特征可视化的流程图；
图16为静脉分割的流程图。
具体实施方式
具体实施方式一、本实施方式所述的基于近红外图像获取最佳静脉穿刺点和方向的系统，所述系统包括近红外图像采集模块、图像预处理模块、ROI区域提取模块、特征增强模块、穿刺点及穿刺方向定位模块、图像配准模块；其中：
所述近红外图像采集模块用于采集人体手背的近红外图像，并将采集的图像发送给图像预处理模块；
所述图像预处理模块用于对接收的图像进行灰度化和滤波降噪处理，获得预处理后的图像，并将预处理后的图像发送给ROI区域提取模块；
所述ROI区域提取模块用于提取和拟合手背图像轮廓，根据提取和拟合的手背图像轮廓定位ROI区域位置；并将ROI区域的图像发送给特征增强模块；
所述特征增强模块用于对ROI区域图像中的静脉进行特征增强，获得特征增强后的图像，并对特征增强后的图像进行灰度反转，获得灰度反转后的图像；
所述穿刺点及穿刺方向定位模块用于提取灰度反转后图像的静脉边缘，根据提取的静脉边缘检测静脉横向宽度和静脉骨架位置，并定位穿刺点及穿刺方向；
所述图像配准模块用于将近红外图像中ROI区域的静脉位置信息配准到自然光学图像中，得到在自然光学图像中的最佳静脉穿刺点以及穿刺方向。
本发明已经初步实现了基于近红外图像的静脉血管识别的功能。由此开发的静脉宽度特征可视化算法和最佳静脉穿刺点及穿刺方向自动定位算法为医护工作人员提供了更加便捷、可视化的信息辅助手段，有助于提高静脉访问成功率，同时也为智能采血系统的开发打下了一定的基础。
具体实施方式二：本实施方式与具体实施方式一不同的是，所述近红外图像采集模块采集图像时，光源频段为850nm的近红外波段。
具体实施方式三：本实施方式与具体实施方式二不同的是，所述近红外图像采集模块采集图像时采用的设备为近红外静脉采集器，所述近红外静脉采集器由中央的图像采集孔以及在图像采集孔的周围呈环形均匀分布的12个LED组成，所述12个LED包括11个红外LED(肉眼不可见)和一个指示用的绿光LED。
具体实施方式四、本实施方式所述的基于近红外图像获取最佳静脉穿刺点和方向的方法，所述方法具体包括以下步骤：
步骤一、采用近红外静脉采集器采集人体手背的近红外图像，采集到的图像如图5所示，并对采集的图像进行图像灰度化处理，获得灰度化处理后的图像；再使用双边滤波器对灰度化处理后的图像进行滤波降噪，获得滤波降噪后的图像；
从滤波降噪后的图像中提取出手背图像轮廓，手背静脉轮廓如图6所示，并计算出手背图像轮廓的矩特征，根据矩特征计算出手背图像轮廓的中心；
对手背图像轮廓进行拟合，计算出手背图像轮廓的宽度width和下边界点纵坐标ydown；并根据手背图像轮廓的中心的横坐标、宽度width以及下边界点纵坐标ydown选取ROI区域；选取的ROI区域如图7所示；
步骤二、提取ROI区域内的静脉结构，获取静脉特征增强图像；并对获得的静脉特征增强图像进行灰度反转，得到灰度反转后的图像；反转后的静脉特征增强图像如图9所示；
步骤三、利用Canny边缘检测方法提取出灰度反转后图像中的静脉边缘，结果如图10所示，再根据提取出的静脉边缘进行静脉横向宽度和静脉骨架位置检测；
根据检测出的静脉骨架位置计算方向修正角度，利用方向修正角度对静脉横向宽度进行修正，获得修正后的静脉横向宽度；效果如图11a)至图11d)所示；
根据修正后的静脉横向宽度进行最佳静脉穿刺点以及穿刺方向定位；效果如图12所示；
步骤四、利用步骤二的方法对人体手背的自然光学图像进行ROI区域的提取，再将近红外图像ROI区域中的静脉位置信息配准到自然光学图像中，得到在自然光学图像中的最佳静脉穿刺点以及穿刺方向。配准效果如图13a)和13b)所示。
具体实施方式五：本实施方式与具体实施方式四不同的是，所述步骤一中，对采集的图像进行图像灰度化处理，获得灰度化处理后的图像，其具体过程为：
对采集的图像的R、G、B三个分量进行加权处理：
Gray(x,y)＝WR*R(x,y)+WG*G(x,y)+WB*B(x,y)
其中，R(x,y)为采集的图像在像素点(x,y)的R分量，G(x,y)为采集的图像在像素点(x,y)的G分量，B(x,y)为采集的图像在像素点(x,y)的B分量，WR为R(x,y)的加权值，WG为G(x,y)的加权值，WB为B(x,y)的加权值，WB为B(x,y)的加权值，WR的取值为0.299，WG的取值为0.587，WB的取值为0.114。
具体实施方式六：本实施方式与具体实施方式五不同的是，所述双边滤波器表示为：
Figure BDA0002984216610000061
其中，g(i,j)为双边滤波器，ω(i,j,k,l)为双边滤波器的模板权值，f(k,l)为中心点位置灰度值；
所述双边滤波器的模板权值ω(i,j,k,l)为：
ω(i,j,k,l)＝ωd(i,j,k,l)*ωr(i,j,k,l)
其中，ωd(i,j,k,l)为空间域核模板权值，ωr(i,j,k,l)为值域核模板权值；
所述空间域核模板权值ωd(i,j,k,l)为：
Figure BDA0002984216610000062
其中，(i,j)为当前点位置，(k,l)为中心点位置，δd为空间域标准差；
所述值域核模板权值ωr(i,j,k,l)为：
Figure BDA0002984216610000063
其中，f(i,j)为当前点位置灰度值，δr为值域标准差，||·||代表范数。
具体实施方式七：本实施方式与具体实施方式四不同的是，所述步骤一中，从滤波降噪后的图像中提取出手背图像轮廓，并计算出手背图像轮廓的矩特征，根据矩特征计算出手背图像轮廓的中心；其具体过程为：
轮廓的零阶矩m00为：
Figure BDA0002984216610000071
其中，f(i′,j′)为手背图像轮廓内的像素点(i′,j′)的灰度值，M为手背图像轮廓的宽度，N为手背图像轮廓的高度，i′＝1,2,…,M，j′＝1,2,…,N；
轮廓的一阶矩m10和m01为：
Figure BDA0002984216610000072
Figure BDA0002984216610000073
则手背图像轮廓的中心
Figure BDA0002984216610000074
为：
Figure BDA0002984216610000075
具体实施方式八：本实施方式与具体实施方式七不同的是，所述ROI区域的横坐标范围为
Figure BDA0002984216610000076
纵坐标范围为
Figure BDA0002984216610000077
Figure BDA0002984216610000078
其中，xcenter为手背图像轮廓的中心的横坐标。
具体实施方式九：本实施方式与具体实施方式八不同的是，所述步骤二中，提取ROI区域内的静脉结构，获取静脉特征增强图像；其具体过程为：
步骤二一、对于ROI区域内的任意一个像素(x,y)，在像素(x,y)处的Hessian(黑塞)矩阵H为：
Figure BDA0002984216610000079
其中，fxx(x,y)为ROI区域内的图像f在x方向上的二阶偏导数，fxy(x,y)为ROI区域内的图像f在两个垂直方向的混合偏导数，fyy(x,y)为ROI区域内的图像f在y方向上的二阶偏导数；
Figure BDA0002984216610000081
Figure BDA0002984216610000082
Figure BDA0002984216610000083
其中，f(x,y)是像素(x,y)的灰度值，G(x,y；σ)为高斯卷积核；
若Hessian矩阵多项式中一次项的系数小于等于t，则像素(x,y)不属于静脉像素，若一次项的系数大于t，则计算出Hessian矩阵的特征值λ1和λ2后，其中λ2＞λ1，判断是否满足λ2/λ1＞k，若不满足，则像素(x,y)不属于静脉像素，若满足，则执行步骤二二；
步骤二二、定义Hessian滤波器如下：
Figure BDA0002984216610000084
其中，U(λ)为像素(x,y)在尺度因子σ下的特征增强值，参数β和c为调节系数，RB和S为中间变量；
Figure BDA0002984216610000085
Figure BDA0002984216610000086
若在ROI区域内的图像中，静脉宽度范围是[r1,r2]，则尺度因子σ的取值范围为[σmin,σmax]，其中，σmin＝r1/4，σmax＝r2/4，将区间[σmin,σmax]平均分割为N个量，再分别计算出像素)x,y)在每个量下的特征增强值，并从中选取出最大的特征增强值：
Figure BDA0002984216610000087
将选取出的最大U值作为像素(x,y)的特征增强值；
步骤二三、同理，对ROI区域内的每一个像素均进行步骤二一至步骤二二的处理后，获得静脉特征增强图像。
参数t和参数k的取值可根据实际情况调整，在本发明中t的取值为0.01，k的取值为2。
具体实施方式十：结合图15说明本实施方式。本实施方式与具体实施方式九不同的是，所述步骤三中，根据提取出的静脉边缘进行静脉横向宽度和静脉骨架位置检测，其具体过程为：
步骤三一、对于灰度反转后图像中的第L行像素，判断处于静脉边缘上的像素点是否为静脉边界点；
在第L行像素中，从左至右对各个处于静脉边缘上的像素点进行遍历，当首次遇到某个像素点的左右两侧像素灰度呈下降趋势时，即位于该像素点左侧的像素点的灰度值比位于该像素点右侧的像素点的灰度值高，则将该像素点记为左边界点x1；
步骤三二、继续寻找x1后的第一个边界点，若寻找出的边界点的左侧像素点的灰度值比右侧像素点的灰度值低，则将寻找出的边界点记为右边界点x2；执行步骤三三；
否则，将寻找出的边界点记为左边界点x1；并重复步骤三二；
步骤三三、若|x2-x1|＞r2或|x2-x1|＜r1，则舍弃边界点x1和x2；否则r1≤|x2-x1|≤r2，则|x2-x1|为静脉横向宽度，
Figure BDA0002984216610000091
为静脉骨架位置；
直至第L行的全部边界点遍历完成，得到第L行包含的全部静脉的横向宽度以及静脉骨架位置；将第i个静脉的横向宽度记为Wi，骨架位置记为ci；
步骤三四、将灰度反转后图像中的第L-d行像素以及第L+d行像素作为辅助行，分别对第L-d行像素以及第L+d行像素进行步骤三一至步骤三三的处理，寻找出与ci对应的骨架位置ci-d和ci+d，ci-d为ci在第L-d行对应的骨架位置，ci+d为ci在第L+d行对应的骨架位置；
将ci-d和ci+d的连线作为第i个静脉在第L行上的静脉方向；
方向修正角度通过以下公式计算：
Figure BDA0002984216610000092
Figure BDA0002984216610000093
其中，
Figure BDA0002984216610000094
和
Figure BDA0002984216610000095
分别为骨架位置ci-d和ci+d的横坐标值，θi为方向修正角度；
则第i个静脉的修正后横向宽度W′i为：
W′i＝Wi*cosθi。
具体实施方式十一：结合图16说明本实施方式。本实施方式与具体实施方式十不同的是，所述步骤三中，根据修正后的静脉横向宽度进行最佳静脉穿刺点以及穿刺方向定位；其具体过程为：
步骤1、静脉分割
由上至下遍历静脉骨架图像，对于静脉骨架图像的首行，分别创建每个静脉骨架点从属的静脉片段集合，对于静脉骨架图像的第L行，判断第L行上的各个静脉骨架点是否从属于已有的静脉片段集合，若是，则将静脉骨架点加入已有的静脉片段集合，否则，将创建新的静脉片段集合；
所述判断第L行上的各个静脉骨架点是否从属于已有的静脉片段集合，其具体过程为：
对于第L行上的任意一个静脉骨架点，计算该静脉骨架点与第L-1行中静脉骨架点的横向最小距离，若计算出的横向最小距离小于等于设定的阈值a，则该静脉骨架点从属于已有的静脉片段集合，否则，该静脉骨架点不从属于已有的静脉片段集合；
阈值a可根据实际使用中，红外摄像模组规格以及成像情况自主调节，本次实验中a的取值为10个像素，a的取值根据静脉的最大宽度r2来设定；切割后图像如图8所示。
步骤2、删除微小静脉
即删除静脉长度小于阈值b的静脉；
本发明中阈值b的取值设置为50个像素；
步骤3、筛选静脉片段
对于每条静脉，遍历其上的每一个静脉片段，静脉片段的长度大于阈值c，分别计算出每个静脉片段上的宽度均值、宽度方差以及方向修正角度均值、方向修正角度方差；
根据宽度方差和方向修正角度方差1：4的权重对每条静脉上的静脉片段的访问价值进行排序，排序越靠前的静脉片段，越适于进行静脉穿刺，用于提供静脉穿刺方案；
对于一个静脉片段，在该静脉片段上的穿刺方案为：穿刺点为该静脉片段的靠近手指侧的端点，穿刺方向为该静脉片段的方向修正角度均值；
本发明中阈值c的取值设置为50个像素；
步骤4、静脉穿刺安全性检验
分别对每个静脉穿刺方案进行验证，只有进针区域全部位于目标静脉片段范围内，且针距离静脉片段边界的距离大于阈值d时，对应的穿刺方案才被保留；否则，舍弃对应的穿刺方案。
本发明中阈值d的取值设置为3个像素；
具体实施方式十二：本实施方式与具体实施方式十一不同的是，所述步骤四中，将近红外图像ROI区域中的静脉位置信息配准到自然光学图像中所采用的方法是透视变换方法。
实施例
下面结合附图进一步说明本发明方法的整个流程。如图1所示为本发明方法的流程图，具体步骤如下：
步骤一：感兴趣区域(ROI)的动态确定
1)图像灰度化
采集图像时光源频段为850nm的近红外波段不可见光。为了避免可见光的干扰，图像采集环境需是封闭的、屏蔽可见光的。但由于我们的图像采集设备很难做到完全屏蔽可见光，因此采集的图像中依然包含着一定量的颜色信息。为了减少后续步骤的计算量，需先将图像灰度化，以消除彩色(R、G、B)信息。
图像灰度化原理：依据人肉眼对R、G、B三个光源颜色分量的敏感程度不同(由高到低依次为绿、红、蓝)，因此将R、G、B三个分量进行加权处理，加权值分别为WR＝0.299,WG＝0.587,WB＝0.114可得到最佳的灰度变换图像：
Gray(x,y)＝0.299*R(x,y)+0.587*G(x,y)+0.114*B(x,y)
2)滤波降噪
这里采用了双边滤波，在降噪的同时，保护并且增强静脉纹路特征。
双边滤波是一种非线性滤波方法，是结合了图像的空间邻近度和像素值相似度的一种折中处理，同时考虑空域信息和灰度相似性，在滤除噪声的同时可以保留原图的边缘信息。具有简单、非迭代、局部的特点。双边滤波器之所以能够做到在平滑去噪的同时还能够很好的保存边缘，是由于其滤波器的核由两个函数生成：空间域核函数(又称定义域核，空间系数或空间域)和值域核函数(又称像素范围域)。
空间域核是由像素位置欧式距离决定的模板权值ωd，公式如下：
Figure BDA0002984216610000111
其中，(i,j)为当前点位置，(k,l)为中心点位置，δd为空间域标准差。
值域核是由像素值的差值决定的模板权值ωr,公式如下：
Figure BDA0002984216610000112
其中，f(i,j)为当前点位置灰度值，f(k,l)为中心点位置灰度值，δr为值域标准差。
将空间域核和值域核两个模板相乘，即可得双边滤波器的模板权值，公式如下：
Figure BDA0002984216610000113
化简之后的双边滤波器可以表示为：
Figure BDA0002984216610000121
对于高斯滤波，仅用空间距离的权值系数核与图像卷积后，确定中心点的灰度值。即认为离中心点越近的点，其权重系数越大。双边滤波中加入了对灰度信息的权重，即在邻域内，灰度值越接近中心点灰度值的点的权重更大，灰度值相差大的点权重越小。此权重大小，则由值域高斯函数确定。借助opencv中提供的bilateralFilter()函数我们可以非常简便地实现双边滤波操作。
4)动态确定感兴趣区域(ROI)：
采集到的图像中包含我们不需要的手背边缘以及背景信息，如果保留的话会极大地增加我们的工作量，使算法的运算速度大幅度降低，降低算法效率。因此，我们有必要动态确定一个只包含手背皮肤和静脉区域的感兴趣区域(ROI)，以便于进行后续处理。
本发明提出了一种动态确定ROI区域的方法，该方法能有效地避免手背平移对静脉识别地影响，具有很好的鲁棒性。在本发明的实验中，采集到的图片中大部分区域为手背区域，所含背景区域较小。因此，图像中最大的轮廓即为手背轮廓。借助提取到的图像的轮廓，计算图像轮廓的矩特征，由此计算图像的中心点。
轮廓零阶矩:
Figure BDA0002984216610000122
当图像为二值图时，m00就是这个图像上白色区域的总和，因此，m00可以用来求二值图像(轮廓，连通域)的面积。
轮廓一阶矩：
Figure BDA0002984216610000123
Figure BDA0002984216610000124
当图像为二值图时,m10就是白色像素关于x坐标的累加和，而m01则是y坐标的累加和。由此，可获得图像的中心：
Figure BDA0002984216610000125
对轮廓进行拟合，计算轮廓的宽度width和下边界点纵坐标ydown。由轮廓的中心横坐标xcenter和轮廓宽度进行ROI区域的选取(因为每次采集图像时手腕进入图像的长度以及手掌位置具有不确定性，因此这里舍弃中心纵坐标xcenter和轮廓高度，进行后续计算)。
经反复验证，当ROI区域设定为：
a)横坐标范围为
Figure BDA0002984216610000131
b)纵坐标范围为
Figure BDA0002984216610000132
所确定的ROI区域效果最好。我们整个ROI区域的动态确定的算法是基于有限的手背近红外图像实现的，因此，我们有必要进行算法的鲁棒性实验，以保证算法在应用于实际场景时的有效性。我们发现，在保证待检测手背以与视野上边界成近似90度的角度进入视野，且进入视野的小臂与手掌尽量保持同一直线，算法的鲁棒性较好，可以满足预期要求。
步骤二：基于Hessian矩阵的多尺度静脉特征增强
1)Hessian矩阵
Hessian矩阵是一个由多元函数的二阶偏导数所构成的方阵，它描述了多元函数的局部曲率。对于多元函数f(x0,x1,…,xn)，若其所有二阶偏导数均存在，则其Hessian矩阵为：
H(f(x))ij＝DiDjf(x)
其中，x为向量(x0,x1,…,xn)，即：
Figure BDA0002984216610000133
在图像边缘检测领域中，要判断某一个像素点是否具有边缘特征，需要考虑该点在输入图像中的局部特征(包括局部极值和极值方向)。因此，可利用Hessian矩阵的特征值以及特征向量获得其所包含的图像边缘信息。
根据线性尺度空间理论(LOG)，对一个函数求导，等于求其与高斯函数导数的卷积。当我们构建一幅图像的尺度空间时，遵循的基本思想是：将原始图像的特征信息，通过卷积的运算来加载到一个受尺度参数σ控制的渐进的偏导信号(高斯核)中。对于二维信号f(x)，它的尺度空间L(x,σ)表达式为：
L(x，y，；σ)＝σγf(x，y)*G(x，y；σ)
其中，σ为空间尺度因子，G(x,y；σ)为高斯卷积核，定义为：
Figure BDA0002984216610000141
因此，图像I在x方向上的二阶方向偏导数为：
Figure BDA0002984216610000142
在y方向上的二阶方向偏导数为：
Figure BDA0002984216610000143
在两个垂直方向的混合偏导数为：
Figure BDA0002984216610000144
综上可得图像I在(x,y)处的Hessian矩阵为：
Figure BDA0002984216610000145
假设Hessian矩阵的第k个特征值为λk，νk为对应的特征向量，由矩阵的特征值和特征向量的定义可知：
Hvk＝Akvk
vk THvk＝λk
Hessian矩阵对图像中不同的结构类型响应不同(主要表现为特征值不同)。表1为二维和三维图像中几种典型的几何结构及其与Hessian矩阵特征值的关系(假设三维图像中特征值的关系为：|λ1|＜|λ2|＜|λ3|；二维图像中特征值的关系为：|λ1|＜|λ2|)。其中，特征值λ1近似为0，其特征向量v1变化程度最小，代表静脉的方向；λ1、λ3的绝对值很大且符号相同，他们的特征向量v2、v3变化程度较大，v2、v3构成的平面垂直于静脉方向。
表1二维和三维图像中几种典型的几何结构及其与Hessian矩阵特征值的关系(H为high，L为low，N为噪声，+/-为特征值符号)
Figure BDA0002984216610000146
Figure BDA0002984216610000151
而本发明中考察的是二维静脉图像，其几何结构包括：圆形结构、线状结构以及噪声。实际上，除了节点处的形状不容易分辨外，静脉图像中其余的几何结构都能大致看作为线性结构。在二维图像中，像素点的Hessian矩阵特征值与其指示的线状结构的关系如图2所示。
2)静脉增强滤波器
由于Hessian矩阵的不同特征值对应不同的结构特性。根据这一原理构造出了一种静脉相似性函数(也称为Hessian滤波器)，用于判断像素点是否具有静脉特征。Frangi定义的二维血管Hessian滤波器如下：
Figure BDA0002984216610000152
其中，
Figure BDA0002984216610000153
Figure BDA0002984216610000154
其中，参数β、c为调节系数，调整RB、S的大小。RB用来分辨局部结构为线状结构还是球状结构，在|λ1|趋近于0时，RB趋近于0，该区域认定为线状结构，即我们所默认的静脉血管区域。实际中，由于背景及噪声产生的随机波动，很可能会产生难以预估的滤波响应。所以，为了减少随机噪声对滤波结果造成的影响，我们引进参数S来剔除背景噪声。由表1可知，当局部区域像素特征为噪声信号时，其所对应的特征值都很小。因此，测度S可以有效地减少滤波结果中的噪声影响。一般来说，我们设置β的取值为0.5，c的取值为相应的Hessian矩阵最大范数的二分之一。结合以上各参数进行分析，可知：当U(λ)的取值为最大值时，我们所考察的局部区域为线状结构，相比之下，其他结构的响应U(λ)值均趋近于0。这样很容易就可以达到抑制背景噪声，增强静脉血管区域的目标。
3)多尺度融合技术
实际应用时，由于静脉血管尺度并不统一，有粗有细，有大有小，如果直接应用上述所提出的滤波器对静脉特征进行增强的话，滤波效果并不理想。这是因为上述所提出的滤波器仅仅只是单一尺度的，其无法对图像中包含的所有尺度的静脉血管都产生强烈的响应。因此，为了改善图像静脉特征增强效果，我们需要采用一种多尺度融合技术。经过反复验证，我们发现：仔细挑选一个合适的尺度范围对于采用多尺度相融合的静脉特征增强非常重要，一个不合适的尺度范围会很大程度地降低图像静脉特征的算法处理效果。对于像静脉血管这样的线状结构来说，如果二阶高斯导数的尺度因子σ与真实的静脉宽度σ0匹配度最高时，静脉增强滤波器的输出响应也将达到最大值。设待检测的静脉图像中，静脉宽度范围是[r1,r2]，那么，理论上来说，尺度因子σ的最小值σmin为r1/4，其最大值σmax取值为r2/4，取值范围为[σmin,σmax]。为了实现对ROI图像中静脉尺度范围内的所有静脉的增强，我们设置尺度因子的步长为Δσ，将区间[σmin,σmax]平均分割为N个量，分别计算N次，从而得到N个尺度上的增强图像，之后选取所有尺度中响应最大的点作为最终的输出：
Figure BDA0002984216610000161
4)快速Hessian算法
上面所述传统的Frangi管状结构滤波方法计算量大且耗时长，这里使用了一种基于快速Hessian矩阵的滤波算法进行静脉增强滤波。通过对Hessian矩阵多项式系数的条件分析，预识别出不满足条件的元素为非血管元素，可省略对该部分元素的Hessian矩阵特征值求解，从而降低计算量，减少滤波耗时。
对于二维图像中任意像素点而言，满足：
Figure BDA0002984216610000162
其中，λ是待求解的特征值，β1、β2是系数，满足：
β1＝λ1+λ2
β2＝λ1λ2
特征值的关系为：|λ1|＜|λ2|。对于静脉血管处像素点而言，其特征值满足：λ1≈0，λ2＞＞λ1。因此，可通过设定对Hessian矩阵多项式系数的条件分析，预识别出不满足条件的元素为非血管元素，省略对该部分元素的Hessian矩阵特征值求解，这里我们设定第一判定条件为β1＞0.01。在计算出剩余部分图像的Hessian矩阵特征值后，利用第二判定条件λ2/λ1＞2，进一步对非血管元素进行筛选排除。
实验结果表明，在保证有效评价图像中管状结构相似度的情况下，该方法可将滤波耗时平均降低至传统Frangi算法的40％以下。
步骤三：静脉宽度特征可视化
1)静脉切面模型
血管模型图如图3a)所示，图3b)为图3a)对应的静脉切面数学模型图。可见，静脉血管灰度分布特征可近似视为高斯分布，这一特征主要体现在：位于静脉血管边缘的像素点，其梯度具有明显的增大趋势，且对于静脉中心线上的像素，其灰度值为局部最小值。
2)静脉横向宽度检测
这里首先要确定的就是血管边缘和血管骨架的位置。鉴于静脉血管灰度分布特征为：位于静脉血管边缘的像素点，其梯度具有明显的增大趋势，且对于静脉中心线上的像素，其灰度值为局部最小值。因此，这里将通过Canny边缘检测的原理进行静脉边缘的提取。
Canny边缘检测可以通过opencv库Canny函数完成，其算法步骤为：
(1)使用高斯滤波器对初始图像进行滤波降噪；
(2)用一阶有限差分近似代替偏导数，计算图像梯度强度和方向；
(3)非极大值抑制(NMS)：寻找像素局部最大值点。沿着梯度方向，比较它前后的梯度值。在沿其边缘方向上的邻域梯度幅值最大，保留；否则，抑制。
(4)双阈值的选取、边缘连接：选取低阈值TH和高阈值TL，取出非极大值抑制后的图像中的最大梯度幅值，将小于TL的点抛弃，赋0；将大于TL的点立即标记(这些点就是边缘点)，赋1。将大于TL，小于TL的点使用8连通区域确定(即：只有与TL像素连接时才会被接受，成为边缘点，赋1)。
而后利用我们所获得的静脉边缘计算静脉宽度和静脉骨架位置，算法步骤如下(以待检测行L为例)：
(1)静脉左边界点检测：对于待检测行L，由左至右进行边缘判定。当遇到边界点时，基于对应点两侧像素点灰度梯度判断，若灰度呈下降趋势，则判断为可能的左边界点，记为x1。
(2)静脉边界预判：寻找像素点x1后的第一个边界点，判断其两侧像素点灰度是否呈上升趋势。若是，则判定为右边界点，记为x2。此两点判断为可能的静脉边界。否则，将此边界点记为新的x1，重新进行本步骤。
(3)静脉边界判别：判断血管左右边界的准确性，删除横向宽度过大或者过小的结果。宽度过大的将会被视作伪静脉考虑；横向宽度过小，是因为这是由于图像噪声造成的，另一方面来说，宽度过小的静脉也不适合进行静脉穿刺访问。
(4)静脉宽度和静脉骨架位置的计算：将|x1-x2|记为静脉横向宽度，将
Figure BDA0002984216610000171
记为静脉骨架位置。
(5)通过以上操作可检测到L行上所有存在的静脉宽度信息，包括静脉骨架位置ci和相应的宽度Wi，我们从左向右依次保存静脉数据。
3)方向修正的静脉宽度
实际中，手背静脉方向并非总是严格的由上至下的，大多数都存在一定的角度差异。因此，我们有必要对上述获得的静脉横向宽度进行方向修正，从而获得真实的静脉宽度特征。
这里我们采用与检测L行相同的宽度信息检测方法，对两辅助行(L-d行和L+d行)进行检测，寻找辅助行上与L静脉骨架点ci对应的骨架点ci-d和ci+d，借助这些静脉骨架点寻找静脉血管的方向信息。两辅助行中骨架点ci-d和ci+d的连接线，即可视为该静脉血管在检测行L上的静脉方向，如图4所示。
方向修正角度通过以下公式计算：
Figure BDA0002984216610000172
Figure BDA0002984216610000173
其中，
Figure BDA0002984216610000181
和
Figure BDA0002984216610000182
分别为辅助行静脉骨架点像素横坐标值。
因此，修正后静脉宽度W′i为：
W′i＝Wi*cosθi
4)最佳静脉穿刺点及穿刺方向定位。
在获取了手背静脉的全部的宽度信息和方向信息之后，我们有必要进一步地提高静脉特征的可视化程度。前文算法支持用户通过鼠标点击图像中某一横行以获取当前行的静脉宽度和方向信息。然而，在实际应用中，对于医护人员而言，这种方式始终不够直观、便捷。因此，我们有必要设计一种针对手背静脉的、可以筛选最佳穿刺点以及穿刺方位的自动定位算法。
我们根据静脉穿刺的特点可以分析得出，适合进行访问的静脉主要需要满足3个条件：①静脉具有一定的长度；②静脉弯曲度较小；③静脉具有一定的宽度。基于此3点的要求，我们可以对静脉依次进行分割、提取，而后借助上述算法处理所得到的静脉纹路的位置和方向信息，筛选出粗壮且位置较浅的静脉，设计定位算法，以实现对适合穿刺的静脉片段的自动定位。最后，在经过穿刺安全性检验之后，将筛选出的质量评价靠前的定位方案提供给医护人员进行选择。
在手背近红外静脉图像的静脉特征增强图像的基础上，我们有必要对手背静脉进行分割，同时依次去除静脉节点、删除微小静脉，然后进一步地选取最有利于穿刺的静脉片段，实现最佳静脉穿刺点及穿刺方向自动定位算法的开发。
算法步骤如下：
(1)静脉分割：上一章节我们获得了静脉骨架位置，由上至下遍历静脉骨架图像，对于首行，由相应的的静脉骨架点创建其从属的静脉集合。对于待检测行L，设定阈值，判断L行上的静脉骨架点是否从属于已有的静脉集合。若是，则加入；否则，将创建新的静脉片段集合。到此，我们完成了对静脉的分割、提取。
(2)删除微小静脉：设定阈值，删除长度过短的、不适合进行静脉访问的静脉。
(3)筛选每条静脉上最适合进行穿刺的片段：对于每条静脉，遍历其上每一个静脉片段(设定阈值，以保证片段长度足够进行静脉穿刺)，计算静脉片段的宽度均值、方差以及静脉方向的均值、方差。依据宽度方差和角度方差1：4的权重对每条静脉上的静脉片段的访问价值进行排序，筛选出每条静脉上最适合进行穿刺的静脉片段。
(4)依据筛选出的静脉片段对静脉的穿刺价值进行排序，将穿刺方案提供给医护人员。
(5)静脉穿刺安全性检验：为保证静脉访问方案的安全性，我们还需要进行穿刺安全性检验。对于算法提供的每一种穿刺方案，我们需要重新验证进针区域是否全部位于目标静脉范围内，且距离静脉边界尚有足够的安全距离，以防止医护人员操作不当时刺破血管，保证算法的安全性，提升患者的治疗体验。如果满足要求，则保留该穿刺方案；否则，舍弃。
步骤四：近红外图像与自然光学图像的配准。
通常来说，在近红外图像与自然光学图像中，同一物体的灰度分布情况通常存在较大的差异，因此，如果我们使用常规的直接基于图像灰度的提取角点的方式，会使得近红外图像中的特征点在自然光学图像中完全找不到对应，极端的情况下甚至大部分不存在同名点。这时候，采用点对点配准的方法正确率将会非常低。如果我们仔细观察对应的近红外图像和自然光学图像，就可以看出它们的边缘特征相似、灰度梯度分布及其相似。由此，我们可以依据手背的轮廓特征对近红外图像与自然光学图像进行配准。
由于近红外图像与自然光学图像拍摄位置、角度相同，因此在两幅图像中，手背的轮廓特征基本相同。在这种情况下，我们可以借助步骤二所述动态确定ROI区域的方法，分别对近红外图像与自然光学图像进行ROI区域提取，所获得的将会是同一个ROI区域。
在确定了近红外图像与自然光学图像的ROI区域后，需要将近红外图像ROI区域中的静脉位置信息配准到自然光学图像中。这里采用了透视变换的方法。
透视变换是一种二维坐标(x,y)到二维坐标(u,v)的线性变换，其数学表达式形式如下：
Figure BDA0002984216610000191
透视变换矩阵A为：
Figure BDA0002984216610000192
令a33＝1，则有：
Figure BDA0002984216610000193
展开：
Figure BDA0002984216610000194
分别代入近红外图像与自然光学图像ROI区域4对角点即可求得变换矩阵，而后依据变换矩阵将静脉信息配准到自然光学图像中。
步骤五：GUI界面设置
GUI界面分为图像区、交互区以及功能区三个部分。图像区用于显示图像，同时在特定功能中支持鼠标按键操作以显示更多的静脉信息；交互区用于显示当前进程，向用户发布功能相关信息；功能区则结合特定进程显示不同功能按键，支持用户按键操作确认参数、切换图像等。几个典型GUI界面如图14a)至图14d)所示。
本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (8)

1.基于近红外图像获取最佳静脉穿刺点和方向的方法，其特征在于，所述方法具体包括以下步骤：
步骤一、采用近红外静脉采集器采集人体手背的近红外图像，并对采集的图像进行图像灰度化处理，获得灰度化处理后的图像；再使用双边滤波器对灰度化处理后的图像进行滤波降噪，获得滤波降噪后的图像；
从滤波降噪后的图像中提取出手背图像轮廓，并计算出手背图像轮廓的矩特征，根据矩特征计算出手背图像轮廓的中心；
对手背图像轮廓进行拟合，计算出手背图像轮廓的宽度width和下边界点纵坐标ydown；并根据手背图像轮廓的中心的横坐标、宽度width以及下边界点纵坐标ydown选取ROI区域；
步骤二、提取ROI区域内的静脉结构，获取静脉特征增强图像；并对获得的静脉特征增强图像进行灰度反转，得到灰度反转后的图像；
步骤三、利用Canny边缘检测方法提取出灰度反转后图像中的静脉边缘，再根据提取出的静脉边缘进行静脉横向宽度和静脉骨架位置检测；
根据检测出的静脉骨架位置计算方向修正角度，利用方向修正角度对静脉横向宽度进行修正，获得修正后的静脉横向宽度；
根据修正后的静脉横向宽度进行最佳静脉穿刺点以及穿刺方向定位；
所述根据修正后的静脉横向宽度进行最佳静脉穿刺点以及穿刺方向定位；其具体过程为：
步骤1、静脉分割
由上至下遍历静脉骨架图像，对于静脉骨架图像的首行，分别创建每个静脉骨架点从属的静脉片段集合，对于静脉骨架图像的第L行，判断第L行上的各个静脉骨架点是否从属于已有的静脉片段集合，若是，则将静脉骨架点加入已有的静脉片段集合，否则，将创建新的静脉片段集合；
所述判断第L行上的各个静脉骨架点是否从属于已有的静脉片段集合，其具体过程为：
对于第L行上的任意一个静脉骨架点，计算该静脉骨架点与第L-1行中静脉骨架点的横向最小距离，若计算出的横向最小距离小于等于设定的阈值a，则该静脉骨架点从属于已有的静脉片段集合，否则，该静脉骨架点不从属于已有的静脉片段集合；
步骤2、删除微小静脉
即删除静脉长度小于阈值b的静脉；
步骤3、筛选静脉片段
对于每条静脉，遍历其上的每一个静脉片段，静脉片段的长度大于阈值c，分别计算出每个静脉片段上的宽度均值、宽度方差以及方向修正角度均值、方向修正角度方差；
根据宽度方差和方向修正角度方差1：4的权重对每条静脉上的静脉片段的访问价值进行排序，排序越靠前的静脉片段，越适于进行静脉穿刺，用于提供静脉穿刺方案；
对于一个静脉片段，在该静脉片段上的穿刺方案为：穿刺点为该静脉片段的靠近手指侧的端点，穿刺方向为该静脉片段的方向修正角度均值；
步骤4、静脉穿刺安全性检验
分别对每个静脉穿刺方案进行验证，只有进针区域全部位于目标静脉片段范围内，且针距离静脉片段边界的距离大于阈值d时，对应的穿刺方案才被保留；否则，舍弃对应的穿刺方案；
步骤四、利用步骤二的方法对人体手背的自然光学图像进行ROI区域的提取，再将近红外图像ROI区域中的静脉位置信息配准到自然光学图像中，得到在自然光学图像中的最佳静脉穿刺点以及穿刺方向。
2.根据权利要求1所述的基于近红外图像获取最佳静脉穿刺点和方向的方法，其特征在于，所述步骤一中，对采集的图像进行图像灰度化处理，获得灰度化处理后的图像，其具体过程为：
对采集的图像的R、G、B三个分量进行加权处理：
Gray(x,y)＝WR*R(x,y)+WG*G(x,y)+WB*B(x,y)
其中，R(x,y)为采集的图像在像素点(x,y)的R分量，G(x,y)为采集的图像在像素点(x,y)的G分量，B(x,y)为采集的图像在像素点(x,y)的B分量，WR为R(x,y)的加权值，WG为G(x,y)的加权值，WB为B(x,y)的加权值，WR的取值为0.299，WG的取值为0.587，WB的取值为0.114。
3.根据权利要求2所述的基于近红外图像获取最佳静脉穿刺点和方向的方法，其特征在于，所述双边滤波器表示为：
Figure FDA0003763536180000021
其中，g(i,j)为双边滤波器，ω(i,j,k,l)为双边滤波器的模板权值，f(k,l)为中心点位置灰度值；
所述双边滤波器的模板权值ω(i,j,k,l)为：
ω(i,j,k,l)＝ωd(i,j,k,l)*ωr(i,j,k,l)
其中，ωd(i,j,k,l)为空间域核模板权值，ωr(i,j,k,l)为值域核模板权值；
所述空间域核模板权值ωd(i,j,k,l)为：
Figure FDA0003763536180000031
其中，(i,j)为当前点位置，(k,l)为中心点位置，δd为空间域标准差；
所述值域核模板权值ωr(i,j,k,l)为：
Figure FDA0003763536180000032
其中，f(i,j)为当前点位置灰度值，δr为值域标准差，||·||代表范数。
4.根据权利要求1所述的基于近红外图像获取最佳静脉穿刺点和方向的方法，其特征在于，所述步骤一中，从滤波降噪后的图像中提取出手背图像轮廓，并计算出手背图像轮廓的矩特征，根据矩特征计算出手背图像轮廓的中心；其具体过程为：
轮廓的零阶矩m00为：
Figure FDA0003763536180000033
其中，f(i′,j′)为手背图像轮廓内的像素点(i′,j′)的灰度值，M为手背图像轮廓的宽度，N为手背图像轮廓的高度，i′＝1,2,…,M，j′＝1,2,…,N；
轮廓的一阶矩m10和m01为：
Figure FDA0003763536180000034
Figure FDA0003763536180000035
则手背图像轮廓的中心
Figure FDA0003763536180000036
为：
Figure FDA0003763536180000037
5.根据权利要求4所述的基于近红外图像获取最佳静脉穿刺点和方向的方法，其特征在于，所述ROI区域的横坐标范围为
Figure FDA0003763536180000038
纵坐标范围为
Figure FDA0003763536180000039
其中，xcenter为手背图像轮廓的中心的横坐标。
6.根据权利要求5所述的基于近红外图像获取最佳静脉穿刺点和方向的方法，其特征在于，所述步骤二中，提取ROI区域内的静脉结构，获取静脉特征增强图像；其具体过程为：
步骤二一、对于ROI区域内的任意一个像素(x,y)，在像素(x,y)处的Hessian矩阵H为：
Figure FDA0003763536180000041
其中，fxx(x,y)为ROI区域内的图像f在x方向上的二阶偏导数，fxy(x,y)为ROI区域内的图像f在两个垂直方向的混合偏导数，fyy(x,y)为ROI区域内的图像f在y方向上的二阶偏导数；
Figure FDA0003763536180000042
Figure FDA0003763536180000043
Figure FDA0003763536180000044
其中，f(x,y)是像素(x,y)的灰度值，G(x,y；σ)为高斯卷积核；
若Hessian矩阵多项式中一次项的系数小于等于t，则像素(x,y)不属于静脉像素，若一次项的系数大于t，则计算出Hessian矩阵的特征值λ1和λ2后，其中λ2＞λ1，判断是否满足λ2/λ1＞k，若不满足，则像素(x,y)不属于静脉像素，若满足，则执行步骤二二；
步骤二二、定义Hessian滤波器如下：
Figure FDA0003763536180000045
其中，U(λ)为像素(x,y)在尺度因子σ下的特征增强值，参数β和c为调节系数，RB和S为中间变量；
Figure FDA0003763536180000046
Figure FDA0003763536180000047
若在ROI区域内的图像中，静脉宽度范围是[r1,r2]，则尺度因子σ的取值范围为[σmin,σmax]，其中，σmin＝r1/4，σmax＝r2/4，将区间[σmin,σmax]平均分割为N个量，再分别计算出像素(x,y)在每个量下的特征增强值，并从中选取出最大的特征增强值：
Figure FDA0003763536180000048
将选取出的最大U值作为像素(x,y)的特征增强值；
步骤二三、同理，对ROI区域内的每一个像素均进行步骤二一至步骤二二的处理后，获得静脉特征增强图像。
7.根据权利要求6所述的基于近红外图像获取最佳静脉穿刺点和方向的方法，其特征在于，所述步骤三中，根据提取出的静脉边缘进行静脉横向宽度和静脉骨架位置检测，其具体过程为：
步骤三一、对于灰度反转后图像中的第L行像素，判断处于静脉边缘上的像素点是否为静脉边界点；
在第L行像素中，从左至右对各个处于静脉边缘上的像素点进行遍历，当遇到某个像素点的左右两侧像素灰度呈下降趋势时，即位于该像素点左侧的像素点的灰度值比位于该像素点右侧的像素点的灰度值高，则将该像素点记为左边界点x1；
步骤三二、继续寻找x1后的第一个边界点，若寻找出的边界点的左侧像素点的灰度值比右侧像素点的灰度值低，则将寻找出的边界点记为右边界点x2；执行步骤三三；
否则，将寻找出的边界点记为左边界点x1；并重复步骤三二；
步骤三三、若|x2-x1|＞r2或|x2-x1|＜r1，则舍弃边界点x1和x2；否则r1≤|x2-x1|≤r2，则|x2-x1|为静脉横向宽度，
Figure FDA0003763536180000051
为静脉骨架位置；
直至第L行的全部边界点遍历完成，得到第L行包含的全部静脉的横向宽度以及静脉骨架位置；将第i个静脉的横向宽度记为Wi，骨架位置记为ci；
步骤三四、将灰度反转后图像中的第L-d行像素以及第L+d行像素作为辅助行，分别对第L-d行像素以及第L+d行像素进行步骤三一至步骤三三的处理，寻找出与ci对应的骨架位置ci-d和ci+d，ci-d为ci在第L-d行对应的骨架位置，ci+d为ci在第L+d行对应的骨架位置；
将ci-d和ci+d的连线作为第i个静脉在第L行上的静脉方向；
方向修正角度通过以下公式计算：
Figure FDA0003763536180000052
Figure FDA0003763536180000053
其中，
Figure FDA0003763536180000054
和
Figure FDA0003763536180000055
分别为骨架位置ci-d和ci+d的横坐标值，θi为方向修正角度；
则第i个静脉的修正后横向宽度Wi′为：
Wi′＝Wi*cosθi。
8.根据权利要求7所述的基于近红外图像获取最佳静脉穿刺点和方向的方法，其特征在于，所述步骤四中，将近红外图像ROI区域中的静脉位置信息配准到自然光学图像中所采用的方法是透视变换方法。
CN202110295569.5A 2021-03-19 2021-03-19 基于近红外图像获取最佳静脉穿刺点和方向的系统及方法 Active CN113011333B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110295569.5A CN113011333B (zh) 2021-03-19 2021-03-19 基于近红外图像获取最佳静脉穿刺点和方向的系统及方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110295569.5A CN113011333B (zh) 2021-03-19 2021-03-19 基于近红外图像获取最佳静脉穿刺点和方向的系统及方法

Publications (2)

Publication Number Publication Date
CN113011333A CN113011333A (zh) 2021-06-22
CN113011333B true CN113011333B (zh) 2022-11-18

Family

ID=76403126

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110295569.5A Active CN113011333B (zh) 2021-03-19 2021-03-19 基于近红外图像获取最佳静脉穿刺点和方向的系统及方法

Country Status (1)

Country Link
CN (1) CN113011333B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113487570B (zh) * 2021-07-06 2024-01-30 东北大学 基于改进的yolov5x网络模型的高温连铸坯表面缺陷检测方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104123703A (zh) * 2014-07-09 2014-10-29 广州中国科学院先进技术研究所 一种保持表皮原色的静脉显像方法
CN107194928A (zh) * 2017-06-15 2017-09-22 华中科技大学同济医学院附属协和医院 一种基于视觉的静脉采血扎针点自动提取方法
CN107749049A (zh) * 2017-09-07 2018-03-02 广州中国科学院先进技术研究所 一种静脉分布显示方法及装置
CN110147769A (zh) * 2019-05-22 2019-08-20 成都艾希维智能科技有限公司 一种手指静脉图像匹配方法
CN112085802A (zh) * 2020-07-24 2020-12-15 浙江工业大学 一种基于双目摄像头获取三维指静脉图像的方法
CN112102227A (zh) * 2020-06-05 2020-12-18 哈尔滨工业大学 一种基于血管尺寸特征的血管最优穿刺点选取方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080298642A1 (en) * 2006-11-03 2008-12-04 Snowflake Technologies Corporation Method and apparatus for extraction and matching of biometric detail
CN104091145B (zh) * 2013-06-02 2018-05-15 广东智冠实业发展有限公司 人体掌脉特征图像采集方法
CN111820919A (zh) * 2020-06-05 2020-10-27 哈工大机器人(中山)无人装备与人工智能研究院 一种采血穿刺控制方法、装置及存储介质

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104123703A (zh) * 2014-07-09 2014-10-29 广州中国科学院先进技术研究所 一种保持表皮原色的静脉显像方法
CN107194928A (zh) * 2017-06-15 2017-09-22 华中科技大学同济医学院附属协和医院 一种基于视觉的静脉采血扎针点自动提取方法
CN107749049A (zh) * 2017-09-07 2018-03-02 广州中国科学院先进技术研究所 一种静脉分布显示方法及装置
CN110147769A (zh) * 2019-05-22 2019-08-20 成都艾希维智能科技有限公司 一种手指静脉图像匹配方法
CN112102227A (zh) * 2020-06-05 2020-12-18 哈尔滨工业大学 一种基于血管尺寸特征的血管最优穿刺点选取方法
CN112085802A (zh) * 2020-07-24 2020-12-15 浙江工业大学 一种基于双目摄像头获取三维指静脉图像的方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
基于MSCT的冠脉血管快速增强和斑块检测方法研究;王丽瑶;《中国优秀硕士学位论文全文数据库（电子期刊） 医药卫生科技辑》;20200615;21-26,33-37 *
基于基准点和 NMI 的手背静脉识别算法研究;屈冰广;《中国优秀硕士学位论文全文数据库（电子期刊） 信息科技辑》;20170415;7-9,12-13 *
基于手背静脉的生物特征识别关键算法研究;耿宏雨;《中国优秀硕士学位论文全文数据库（电子期刊） 信息科技辑》;20160815;15-16 *
近红外静脉图像特征增强及可视化算法;任泓宇;《中国优秀硕士学位论文全文数据库（电子期刊） 信息科技辑》;20150215;10-14,28-47 *

Also Published As

Publication number Publication date
CN113011333A (zh) 2021-06-22

Similar Documents

Publication Publication Date Title
CN109859203B (zh) 基于深度学习的缺陷牙齿图像识别方法
Grau et al. Automatic localization of cephalometric landmarks
Mamonov et al. Automated polyp detection in colon capsule endoscopy
CN105760841B (zh) 一种身份识别方法及系统
KR102206621B1 (ko) 딥러닝 알고리즘을 이용한 근감소증 분석 프로그램 및 애플리케이션
CN111933275A (zh) 一种基于眼动与面部表情的抑郁评估系统
CN114170201B (zh) 基于边缘光流信息的非接触式呼吸率检测方法及系统
CN110298273A (zh) 一种基于多光谱图像的3d指静脉提取方法及系统
CN116630762B (zh) 一种基于深度学习的多模态医学图像融合方法
CN114022554A (zh) 一种基于yolo的按摩机器人穴位检测与定位方法
CN111339828B (zh) 基于红外影像和超声多普勒结合的静脉显影识别方法
CN113011333B (zh) 基于近红外图像获取最佳静脉穿刺点和方向的系统及方法
CN104715459B (zh) 血管图像增强方法
Mohammed et al. Digital medical image segmentation using fuzzy C-means clustering
Tavakoli et al. Automated optic nerve head detection in fluorescein angiography fundus images
Wei et al. Automatic recognition of major fissures in human lungs
JP5740403B2 (ja) 網膜の異常を検出するためのシステムおよび方法
CN108416792B (zh) 基于活动轮廓模型的医学计算机断层扫描图像分割方法
Li et al. A novel method for low-contrast and high-noise vessel segmentation and location in venipuncture
Tan et al. A real-time image analysis system for computer-assisted diagnosis of neurological disorders
Waghulde et al. Detection of skin cancer lesions from digital images with image processing techniques
Kang et al. Segmentation of coronary arteries based on transition region extraction
Yan et al. Segmentation of pulmonary parenchyma from pulmonary CT based on ResU-Net++ model
Luo et al. An Optic Disc Segmentation Method Based on Active Contour Tracking.
Zamani et al. A Fast and Reliable Three-Dimensional Centerline Tracing: Application to Virtual Cochlear Implant Surgery

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant