CN116138760A - 一种自适应增强的激光散斑衬比血流成像方法 - Google Patents

一种自适应增强的激光散斑衬比血流成像方法 Download PDF

Info

Publication number
CN116138760A
CN116138760A CN202310190366.9A CN202310190366A CN116138760A CN 116138760 A CN116138760 A CN 116138760A CN 202310190366 A CN202310190366 A CN 202310190366A CN 116138760 A CN116138760 A CN 116138760A
Authority
CN
China
Prior art keywords
image
value
speckle
calculating
original
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.)
Pending
Application number
CN202310190366.9A
Other languages
English (en)
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.)
Tianjin Polytechnic University
Original Assignee
Tianjin Polytechnic University
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 Tianjin Polytechnic University filed Critical Tianjin Polytechnic University
Priority to CN202310190366.9A priority Critical patent/CN116138760A/zh
Publication of CN116138760A publication Critical patent/CN116138760A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
    • A61B5/026—Measuring blood flow
    • A61B5/0261—Measuring blood flow using optical means, e.g. infrared light
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0033—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • 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

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Biophysics (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Pathology (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Physics & Mathematics (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Physiology (AREA)
  • Quality & Reliability (AREA)
  • Hematology (AREA)
  • Cardiology (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种自适应增强的激光散斑衬比血流成像方法；利用激光散斑衬比成像系统探测待测组织表面，使用CMOS相机采集原始散斑图像序列；根据所述原始散斑图像序列，计算得到初始散斑衬比图像；对所述初始散斑衬比图像进行图像分割，得到分割图像；根据所述原始散斑图像序列和分割图像，计算得到自适应增强的散斑衬比图像；本发明提供的激光散斑衬比血流成像方法能够提高血管边缘和微小血管的成像效果，使血管的衬比图像更为清晰，反映更多细节；在较深层血管成像中，能够提高血管的可视化效果，具有更高的噪声衰减和更好的成像质量。

Description

一种自适应增强的激光散斑衬比血流成像方法
技术领域
本发明涉及医学图像处理技术领域，特别涉及一种自适应增强的激光散斑衬比血流成像方法。
背景技术
激光散斑衬比血流成像技术是一种实时的、非扫描式的全场血流动力学成像技术，具有非接触、无创伤、无需造影剂、高时空分辨率、仪器结构简单等优势；广泛应用于眼科、脑科学、微循环及术中监测等基础研究和临床应用中；使用激光散斑衬比血流成像技术能够反映出不同血管在同一时间的形态和流速，也能反映出同一血管在不同时刻的血流流速变化，具有长期、动态监测血流动力学的能力；心脑血管疾病、高血压等都与微循环状态密切相关，血流是衡量微循环状态的关键参数之一，实现血流的实时在体测量对疾病预防具有重要意义。
传统重建散斑衬比图像的方法是对空间窗、时间窗或时空窗内的所有像素计算光强的标准差和均值，来计算窗口中心点的散斑衬比值，以此来重建二维散斑衬比图像；在较深层血管成像中，使用传统的散斑衬比成像方法重建的图像含有较高的噪声，动态血管区域与静态组织区域对比度噪声比低；同时对血管边缘以及微小血管的成像质量很差，甚至无法成像出来，这不利于较深层微循环中对血流的监测。
发明内容
本发明旨在解决上述现有技术中的不足，提供一种自适应增强的激光散斑衬比血流成像方法。
本发明所采用的技术方案是：一种自适应增强的激光散斑衬比血流成像方法，所述方法包括如下步骤：
S100，利用激光散斑衬比成像系统测定待测生物组织，获得原始散斑图像序列；
S200，根据所述原始散斑图像序列，计算得到初始散斑衬比图像；
S300，根据所述初始散斑衬比图像，计算得到分割图像；
S400，根据所述原始散斑图像序列和分割图像，计算得到自适应增强的散斑衬比图像。
进一步的，在S100中，利用激光散斑衬比成像系统测定待测生物组织，获得原始散斑图像序列的方法为：
S101，设定所述激光散斑衬比成像系统的CMOS相机的曝光时间为Tms，即在每秒内采集
Figure BSA0000296668000000011
帧；
S102，设定在采集时间t秒内，设置所述激光散斑衬比成像系统的激光光源打开输出，经过扩束镜、反射镜后，照射在待测生物组织表面进行测定；
S103，所述激光散斑衬比成像系统的CMOS相机采集所述待测生物组织表面漫反射以及组织内散射的光束，形成原始散斑图像；
S104，设定在采集时间t秒内，所述CMOS相机总记录
Figure BSA0000296668000000012
帧所述原始散斑图像，由计算机进行灰度预处理后得到原始散斑图像序列。
进一步的，在S200中，根据所述原始散斑图像序列，计算得到初始散斑衬比图像的方法为：
S201，设定所述原始散斑图像序列为I(x，y，z)，其中I(x，y，z)表示在第z帧所述原始散斑图像在像素点坐标(x，y)处的灰度值，其中设定图像分辨率大小为[M，N]，x的取值范围为[1，M]，y的取值范围为[1，N]，z的取值范围为
Figure BSA0000296668000000021
S202，选取滑动空间窗口大小为V×V，滑动空间窗口遍历像素点坐标(x，y)的取值范围：x∈[1，M]，y∈[1，N]；根据原始散斑图像序列I(x，y，z)，选取滑动空间窗口下对应的原始散斑图像子帧数据IV(x，y，z)，其中x的取值范围为[1，V]，y的取值范围为[1，V]，z的取值范围为
Figure BSA0000296668000000022
计算空间散斑衬比图像Kz(x，y，z)：
Figure BSA0000296668000000023
其中
Figure BSA0000296668000000024
S203，根据空间散斑衬比图像Kz(x，y，z)，计算第z帧像素点(x，y)处四个方向上的梯度KN(x，y)z，KS(x，y)z，KW(x，y)z，KE(x，y)z，根据四个方向上的梯度计算四个方向的扩散系数dN(x，y)z，ds(x，y)z，dW(x，y)z，dE(x，y)z，计算得到初始散斑衬比图像K(x，y，z)(n+1)：
K(x，y，z)(n+1)＝K(x，y，z)n+λ×(N+S+WE+E)
上式为得到所述初始散斑衬比图像K(x，y，z)(n+1)的迭代式，n为迭代次数，λ为扩散系数；
其中：
N＝dN(x，y)z×KN(x，y)z，
Figure BSA0000296668000000025
KN(x，y)z＝Kz(x，y-1，z)-Kz(x，y，z)；
S＝dS(x，y)z×KS(x，y)z，
Figure BSA0000296668000000028
KS(x，y)z＝Kz(x，y+1，z)-Kz(x，y，z)；
WE＝dW(x，y)z×KW(x，y)z，
Figure BSA0000296668000000029
KW(x，y)z＝Kz(x-1，y，z)-Kz(x，y，z)；
E＝dE(x，y)z×KE(x，y)z)，
Figure BSA0000296668000000034
KE(x，y)z＝Kz(x+1，y，z)-K(x，y，z)；
其中Kthr为边缘参数。
进一步的，在S300中，根据所述初始散斑衬比图像，计算得到分割图像的方法为：
S301，根据所述初始散斑衬比图像K(x，y，z)(n+1)，使用分割算法分割动态区域、静态区域，得到分割图像Se(x，y)：
Figure BSA0000296668000000032
Se(x，y)＝Se(x1，y1)+Se(x2，y2)+Se(x3，y3)
其中Se(x1，y1)表示动态区域，Se(x2，y2)表示静态区域，Se(x3，y3)表示过渡区域。
进一步的，在S400中，根据所述原始散斑图像序列和分割图像，计算得到自适应增强的散斑衬比图像的方法为：
S401，设定选择的滑动空间窗口大小为W×W，根据所述原始散斑图像序列I(x，y，z)和所述分割图像Se(x，y)，滑动空间窗口遍历像素点坐标(x，y)的取值范围：x∈[1，M]，y∈[1，N]；
S402，选取滑动窗口下对应的原始散斑图像子帧数据IW(x，y，z)和分割图像子帧数据SeW(x，y)，其中x的取值范围为[1，W]，y的取值范围为[1，W]，z的取值范围为
Figure BSA0000296668000000033
S403，计算SeW(x，y)中像素点(x，y)为0的个数C0、像素点为1的个数C1、像素点为0.5的个数C0.5；若C0≠W2或C1≠W2或C0.5≠W2，则说明所述滑动窗口下对应的原始散斑图像子帧数据IW(x，y，z)处于边界区域，使用自适应窗口大小计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka(x，y，z)，即若SeW(x，y)中心点的值为1，则使用SeW(x，y)中等于1所对应的像素点IW1(x，y，z)的值计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka1(x，y，z)；同理，若SeW(x，y)中心点的值为0，则使用SeW(x，y)中等于0所对应的像素点IW0(x，y，z)的值计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka0(x，y，z)；若SeW(x，y)中心点的值为0.5，则使用SeW(x，y)中等于0.5所对应的像素点IW0.5(x，y，z)的值计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka0.5(x，y，z)；
若
Figure BSA0000296668000000041
则Ka(x，y，z)＝Ka1(x，y，z)
Figure BSA0000296668000000042
Figure BSA0000296668000000043
(1)(2)式中，i，j需满足SeW(i，j)＝1；
若
Figure BSA0000296668000000044
则Ka(x，y，z)＝Ka0(x，y，z)
Figure BSA0000296668000000045
Figure BSA0000296668000000046
(3)(4)式中，i，j需满足SeW(i，j)＝0；
若
Figure BSA0000296668000000047
则Ka(x，y，z)＝Ka0.5(x，y，z)
Figure BSA0000296668000000048
Figure BSA0000296668000000049
(5)(6)式中，i，j需满足SeW(i，j)＝0.5；
其中
Figure BSA00002966680000000410
表示对
Figure BSA00002966680000000411
四舍五入取整，z代表第z帧；
S404，若C0＝W2或C1＝W2或C0.5＝W2，则表示所述滑动窗口下对应的原始散斑图像子帧数据IW(x，y，z)不处于边界区域，采用方差判据计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka(x，y，z)；
其中，在S404中，采用方差判据计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka(x，y，z)的方法为：
S4041，计算IW(x，y，z)的窗口方差dall和4个方向(0°、45°、90°、135°)上的方差d0°，d45°，d90°，d135°，构成dn(z)：
dn(z)＝{d0°，d45°，d90°，d135°，dall}
Figure BSA0000296668000000051
Figure BSA0000296668000000052
Ij(z)表示第z帧对应方向上的像素光强值，
Figure BSA0000296668000000053
表示第z帧对应方向上像素光强的均值，
Figure BSA0000296668000000054
表示整个窗口的像素光强均值；
S4042，若C1＝W2说明该区域属于动态区域，采用方差最小方向上的像素点的值计算Ka(x，y，z)：
ddynamic(z)＝argdmin[var(dn(z))]
Figure BSA0000296668000000055
ddynamic(z)表示动态区域方向，也即方差最小方向，σ(ddynamic(z))表示方差最小方向上的像素点的值的标准差，
Figure BSA0000296668000000056
表示方差最小方向上的像素点的值的均值；
S4043，若C0＝W2说明该区域属于静态区域，采用方差最大方向上的像素点的值计算ka(x，y，z)：
dstatic(z)＝argdmax[var(dn(z))]
Figure BSA0000296668000000057
dstatic(z)表示静态区域方向，也即方差最大方向，σ(dstatic(z))表示方差最大方向上的像素点的值的标准差，
Figure BSA0000296668000000058
表示方差最大方向上的像素点的值的均值；
S4044，若C0.5＝W2说明该区域属于过渡区域，采用方差中值方向上的像素点的值计算Ka(x，y，z)：
dtransition(z)＝argdmedian[var(dn(z))]
Figure BSA0000296668000000059
dtransition(z)表示过渡区域方向，也即方差中值方向，σ(dtransition(z))表示方差中值方向上的像素点的值的标准差，
Figure BSA0000296668000000061
表示方差中值方向上的像素点的值的均值；
S405，滑动空间窗口遍历像素点坐标(x，y)的取值范围：x∈[1，M]，y∈[1，N]，即循环S402、S403、S404、S4041、S4042、S4043、S4044可得到新的散斑衬比图像Kz(x，y，z)。
如上所述，本发明所述的一种自适应增强的激光散斑衬比血流成像方法，具有以下有益效果：
(1)基于原始散斑图像序列和分割图像，通过对动静态边界区域自适应窗口大小处理，能够提高血管边缘和微小血管的成像效果，使血管的衬比图像更为清晰，反映更多细节；
(2)基于原始散斑图像序列和分割图像，分割出动静态区域，通过对动态、静态区域使用不同的方差判据进行成像，能够提高动静态区域的对比度；在较深层血管成像中，能够提高血管的可视化效果，具有更高的噪声衰减和更好的成像质量。
附图说明
通过对结合附图所示出的实施方案进行详细说明，本公开的上述及其他特征将更加明显，本公开附图中相同的参考标号表示相同或相似的元素，显而易见的，下面描述中的附图仅仅是本公开的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，在附图中：
图1是本发明的一种改进的激光散斑衬比血流成像方法于一实施例中的流程图；
图2是本发明的一种改进的激光散斑衬比血流成像方法于一实施例中的激光散斑衬比成像系统的硬件结构框图组成示意图；
图3是本发明的一种改进的激光散斑衬比血流成像方法于一实施例中基于公开的原始散斑图像采用空间散斑衬比成像方法与采用本发明成像方法的一个对比图；
图4是本发明的一种改进的激光散斑衬比血流成像方法于一实施例中基于仿体实验采集到得原始散斑图像采用空间散斑衬比成像方法与采用本发明成像方法的一个对比图。
具体实施方式
以下将结合实施例和附图对本公开的构思、具体结构及产生的技术效果进行清楚、完整的描述，以充分地理解本公开的目的、方案和效果；需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互结合；以下实施例中所提供的图示仅以示意方式说明本发明的基本构思，旨在用于解释本申请，而不能理解为对本申请的限制。
本发明所述的一种自适应增强的激光散斑衬比血流成像方法，通过对动静态边界区域自适应窗口大小处理，能够提高血管边缘和微小血管的成像效果；通过对动态、静态区域使用不同的方差判据进行成像，能够提高动静态区域的对比度；在较深层血管成像中，相比于传统的激光散斑衬比成像方法具有更高的噪声衰减和更好的成像质量。
如图1所示为根据本发明一种自适应增强的激光散斑衬比血流成像方法的流程图，下面结合图1来阐述根据本发明的实施方式的一种自适应增强的激光散斑衬比血流成像方法。
本公开提出一种自适应增强的激光散斑衬比血流成像方法，所述方法具体包括以下步骤：
S100，利用激光散斑衬比成像系统测定待测生物组织，获得原始散斑图像序列；
S200，根据所述原始散斑图像序列，计算得到初始散斑衬比图像；
S300，根据所述初始散斑衬比图像，计算得到分割图像；
S400，根据所述原始散斑图像序列和分割图像，计算得到自适应增强的散斑衬比图像。
进一步的，在S100中，利用激光散斑衬比成像系统测定待测生物组织，获得原始散斑图像序列的方法为：
S101，设定所述激光散斑衬比成像系统的CMOS相机的曝光时间为Tms，即在每秒内采集
Figure BSA0000296668000000071
帧；
优选地，T值取值范围为5ms～10ms；
S102，设定在采集时间t秒内，设置所述激光散斑衬比成像系统的激光光源打开输出，经过扩束镜、反射镜后，照射在待测生物组织表面进行测定；
S103，所述激光散斑衬比成像系统的CMOS相机采集所述待测生物组织表面漫反射以及组织内散射的光束，形成原始散斑图像；
S104，设定在采集时间t秒内，所述CMOS相机总记录
Figure BSA0000296668000000072
帧所述原始散斑图像，由计算机进行灰度预处理后得到原始散斑图像序列。
可选地，本具体实施例中的激光散斑衬比成像系统的硬件结构构成如图2所示，可见成像系统的硬件组成包括光源模块、透镜光路模块以及成像采集模块；
光源模块包括激光驱动和激光器，优选地激光器的中心波长为785nm；
透镜光路模块包括扩束镜、反射镜；激光通过扩束镜扩束后，再经过反射镜照射到待测生物组织；
成像采集模块包括相机和处理器，相机将采集到的信号传输给处理器，记录一定采集时间内得到的原始散斑图像序列。
进一步的，在S200中，根据所述原始散斑图像序列，计算得到初始散斑衬比图像的方法为：
S201，设定所述原始散斑图像序列为I(x，y，z)，其中I(x，y，z)表示在第z帧所述原始散斑图像在像素点坐标(x，y)处的灰度值，其中设定图像分辨率大小为[M，N]，x的取值范围为[1，M]，y的取值范围为[1，N]，z的取值范围为
Figure BSA0000296668000000073
S202，选取滑动空间窗口大小为V×V，滑动空间窗口遍历像素点坐标(x，y)的取值范围：x∈[1，M]，y∈[1，N]；根据原始散斑图像序列I(x，y，z)，选取滑动空间窗口下对应的原始散斑图像子帧数据IV(x，y，z)，其中x的取值范围为[1，V]，y的取值范围为[1，V]，z的取值范围为
Figure BSA0000296668000000074
计算空间散斑衬比图像Kz(x，y，z)：
Figure BSA0000296668000000075
其中
Figure BSA0000296668000000076
可选地，本实施例中选择V＝5；
S203，根据空间散斑衬比图像Kz(x，y，z)，计算第z帧像素点(x，y)处四个方向上的梯度KN(x，y)z，KS(x，y)z，KW(x，y)z，KE(x，y)z，根据四个方向上的梯度计算四个方向的扩散系数dN(x，y)z，dS(x，y)z，dW(x，y)z，dE(x，y)z，计算得到初始散斑衬比图像K(x，y，z)(n+1)：
K(x，y，z)(n+1)＝K(x，y，z)n+λ×(N+S+WE+E)
上式为得到所述初始散斑衬比图像K(x，y，z)(n+1)的迭代式，n为迭代次数，λ为扩散系数；
可选地，本实施例中选择n＝50，λ＝0.25；
其中：
N＝dN(x，y)z×KN(x，y)z，
Figure BSA0000296668000000081
KN(x，y)z＝Kz(x，y-1，z)-Kz(x，y，z)；
S＝dS(x，y)z×KS(x，y)z，
Figure BSA0000296668000000082
KS(x，y)z＝Kz(x，y+1，z)-Kz(x，y，z)；
WE＝dW(x，y)z×KW(x，y)z，
Figure BSA0000296668000000083
KW(x，y)z＝Kz(x-1，y，z)-Kz(x，y，z)；
E＝dE(x，y)z×KE(x，y)z)，
Figure BSA0000296668000000084
KE(x，y)z＝Kz(x+1，y，z)-Kz(x，y，z)；
其中Kthr为边缘参数。
可选地，本实施例中选择
Figure BSA0000296668000000085
Figure BSA0000296668000000086
为空间散斑衬比图像的均值。
进一步的，在S300中，根据所述初始散斑衬比图像，计算得到分割图像的方法为：
S301，根据所述初始散斑衬比图像K(x，y，z)(n+1)，使用分割算法分割动态区域、静态区域，得到分割图像Se(x，y)：
Figure BSA0000296668000000087
Se(x，y)＝Se(x1，y1)+Se(x2，y2)+Se(x3，y3)
其中Se(x1，y1)表示动态区域，Se(x2，y2)表示静态区域，Se(x3，y3)表示过渡区域。
可选地，本实施例中选择k-means聚类分割。
进一步的，在S400中，根据所述原始散斑图像序列和分割图像，计算得到自适应增强的散斑衬比图像的方法为：
S401，设定选择的滑动空间窗口大小为W×W，根据所述原始散斑图像序列I(x，y，z)和所述分割图像Se(x，y)，滑动空间窗口遍历像素点坐标(x，y)的取值范围：x∈[1，M]，y∈[1，N]；
可选地，本实施例中选择W＝5；
S402，选取滑动窗口下对应的原始散斑图像子帧数据IW(x，y，z)和分割图像子帧数据SeW(x，y)，其中x的取值范围为[1，W]，y的取值范围为[1，W]，z的取值范围为
Figure BSA0000296668000000091
S403，计算SeW(x，y)中像素点(x，y)为0的个数C0、像素点为1的个数C1、像素点为0.5的个数C0.5；若C0≠W2或C1≠W2或C0.5≠W2，则说明所述滑动窗口下对应的原始散斑图像子帧数据IW(x，y，z)处于边界区域，使用自适应窗口大小计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka(x，y，z)，即若SeW(x，y)中心点的值为1，则使用SeW(x，y)中等于1所对应的像素点IW1(x，y，z)的值计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka1(x，y，z)；同理，若SeW(x，y)中心点的值为0，则使用SeW(x，y)中等于0所对应的像素点IW0(x，y，z)的值计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka0(x，y，z)；若SeW(x，y)中心点的值为0.5，则使用SeW(x，y)中等于0.5所对应的像素点IW0.5(x，y，z)的值计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka0.5(x，y，z)；
若
Figure BSA0000296668000000092
则Ka(x，y，z)＝Ka1(x，y，z)
Figure BSA0000296668000000093
Figure BSA0000296668000000094
(1)(2)式中，i，j需满足SeW(i，j)＝1；
若
Figure BSA0000296668000000101
则Ka(x，y，z)＝Ka0(x，y，z)
Figure BSA0000296668000000102
Figure BSA0000296668000000103
(3)(4)式中，i，j需满足SeW(i，j)＝0；
若
Figure BSA0000296668000000104
则Ka(x，y，z)＝Ka0.5(x，y，z)
Figure BSA0000296668000000105
Figure BSA0000296668000000106
(5)(6)式中，i，j需满足SeW(i，j)＝0.5；
其中
Figure BSA0000296668000000107
表示对
Figure BSA0000296668000000108
四舍五入取整，z代表第z帧；
S404，若C0＝W2或C1＝W2或C0.5＝W2，则表示所述滑动窗口下对应的原始散斑图像子帧数据IW(x，y，z)不处于边界区域，采用方差判据计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka(x，y，z)；
其中，在S404中，采用方差判据计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka(x，y，z)的方法为：
S4041，计算IW(x，y，z)的窗口方差dall和4个方向(0°、45°、90°、135°)上的方差d0°，d45°，d90°，d135°，构成dn(z)：
dn(z)＝({0°，d45°，d90°，d135°，dall}
Figure BSA0000296668000000109
Figure BSA00002966680000001010
Ij(z)表示第z帧对应方向上的像素光强值，
Figure BSA0000296668000000111
表示第z帧对应方向上像素光强的均值，
Figure BSA0000296668000000112
表示整个窗口的像素光强均值；
S4042，若C1＝W2说明该区域属于动态区域，采用方差最小方向上的像素点的值计算Ka(x，y，z)
ddynamic(z)＝argdmin[var(dn(z))]
Figure BSA0000296668000000113
ddynamic(z)表示动态区域方向，也即方差最小方向，σ(ddynamic(z))表示方差最小方向上的像素点的值的标准差，
Figure BSA0000296668000000114
表示方差最小方向上的像素点的值的均值；
S4043，若C0＝W2说明该区域属于静态区域，采用方差最大方向上的像素点的值计算Ka(x，y，z)：
dstatic(z)＝argdmax[var(dn(z))]
Figure BSA0000296668000000115
dstatic(z)表示静态区域方向，也即方差最大方向，σ(dstatic(z))表示方差最大方向上的像素点的值的标准差，
Figure BSA0000296668000000116
表示方差最大方向上的像素点的值的均值；
S4044，若C0.5＝W2说明该区域属于过渡区域，采用方差中值方向上的像素点的值计算Ka(x，y，z)：
dtransition(z)＝argdmedian[var(dn(z))]
Figure BSA0000296668000000117
dtransition(z)表示过渡区域方向，也即方差中值方向，σ(dtransition(z))表示方差中值方向上的像素点的值的标准差，
Figure BSA0000296668000000118
表示方差中值方向上的像素点的值的均值；
S405，滑动空间窗口遍历像素点坐标(x，y)的取值范围：x∈[1，M]，y∈[1，N]，即循环S402、S403、S404、S4041、S4042、S4043、S4044可得到新的散斑衬比图像Ka(x，y，z)。
参照图3为一实施例中基于公开的原始散斑图像，左图为使用空间散斑衬比成像方法得到的伪彩色图像示意图，右图为使用本实例方法得到的自适应增强的散斑衬比血流图像伪彩色示意图；右图较左图显示了更多微血管细节，血管的清晰度得到有效提升；证明了本发明方法能够使血管的衬比图像更为清晰，能够提高血管边缘和微血管的成像效果，反映更多的细节。
参照图4为开展仿体实验，使用环氧树脂制作仿体模型，在环氧树脂中加入不同含量的二氧化钛和印度墨水模拟真皮层和表皮层，使用内径1mm、外径2mm的毛细玻璃管模拟血管，其中毛细玻璃管距离表面深度为400μm，通过注射泵以10mm/s的速度往毛细玻璃管中注射浓度为3％的脂肪乳溶液模拟血液，通过相机采集得到原始散斑图像序列；图4左图为使用空间散斑衬比成像方法得到的伪彩色图像示意图，右图为使用本实例方法得到的自适应增强的散斑衬比血流图像伪彩色示意图；右图相比于左图血管区域与背景区域对比度更高，血管区域更明显；证明了本发明方法在较深层血管成像中，能够提高血管的可视化效果，具有更高的噪声衰减和更好的成像质量。
尽管本发明的实施方案已公开如上，但本发明并不局限于上述实施方式，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节。

Claims (5)

1.一种自适应增强的激光散斑衬比血流成像方法，其特征在于，所述方法包括以下步骤：
S100，利用激光散斑衬比成像系统测定待测生物组织，获得原始散斑图像序列；
S200，根据所述原始散斑图像序列，计算得到初始散斑衬比图像；
S300，根据所述初始散斑衬比图像，计算得到分割图像；
S400，根据所述原始散斑图像序列和分割图像，计算得到自适应增强的散斑衬比图像。
2.根据权利要求1所述的一种自适应增强的激光散斑衬比血流成像方法，其特征在于，在S100中，利用激光散斑衬比成像系统测定待测生物组织，获得原始散斑图像序列的方法为：
S101，设定所述激光散斑衬比成像系统的CMOS相机的曝光时间为Tms，即在每秒内采集
Figure FSA0000296667990000011
帧；
S102，设定在采集时间t秒内，设置所述激光散斑衬比成像系统的激光光源打开输出，经过扩束镜、反射镜后，照射在待测生物组织表面进行测定；
S103，所述激光散斑衬比成像系统的CMOS相机采集所述待测生物组织表面漫反射以及组织内散射的光束，形成原始散斑图像；
S104，设定在采集时间t秒内，所述CMOS相机总记录
Figure FSA0000296667990000012
帧所述原始散斑图像，由计算机进行灰度预处理后得到原始散斑图像序列。
3.根据权利要求1所述的一种自适应增强的激光散斑衬比血流成像方法，其特征在于，在S200中，根据所述原始散斑图像序列，计算得到初始散斑衬比图像的方法为：
S201，设定所述原始散斑图像序列为I(x，y，z)，其中I(x，y，z)表示在第z帧所述原始散斑图像在像素点坐标(x，y)处的灰度值，其中设定图像分辨率大小为[M，N]，x的取值范围为[1，M]，y的取值范围为[1，N]，z的取值范围为
Figure FSA0000296667990000013
S202，选取滑动空间窗口大小为V×V，滑动空间窗口遍历像素点坐标(x，y)的取值范围：x∈[1，M]，y∈[1，N]；根据原始散斑图像序列I(x，y，z)，选取滑动空间窗口下对应的原始散斑图像子帧数据IV(x，y，z)，其中x的取值范围为[1，V]，y的取值范围为[1，V]，z的取值范围为
Figure FSA0000296667990000014
计算空间散斑衬比图像Ks(x，y，z)：
Figure FSA0000296667990000015
其中
Figure FSA0000296667990000016
S203，根据空间散斑衬比图像Ks(x，y，z)，计算第z帧像素点(x，y)处四个方向上的梯度KN(x，y)z，KS(x，y)z，KW(x，y)z，KE(x，y)z，根据四个方向上的梯度计算四个方向的扩散系数dN(x，y)z，dS(x，y)z，dW(x，y)z，dE(x，y)z，计算得到初始散斑衬比图像K(x，y，z)(n+1)：
K(x，y，z)(n+1)＝K(x，y，z)n+λ×(N+S+WE+E)
上式为得到所述初始散斑衬比图像K(x，y，z)(n+1)的迭代式，n为迭代次数，λ为扩散系数；
其中：
N＝dN(x，y)z×KN(x，y)z，
Figure FSA0000296667990000021
S＝dS(x，y)z×KS(x，y)z，
Figure FSA0000296667990000022
WE＝dW(x，y)z×KW(x，y)z，
Figure FSA0000296667990000023
E＝dE(x，y)z×KE(x，y)z)，
Figure FSA0000296667990000024
其中Kthr为边缘参数。
4.根据权利要求1所述的一种自适应增强的激光散斑衬比血流成像方法，其特征在于，在S300中，根据所述初始散斑衬比图像，计算得到分割图像的方法为：
S301，根据所述初始散斑衬比图像K(x，y，z)(n+1)，使用分割算法分割动态区域、静态区域，得到分割图像Se(x，y)：
Figure FSA0000296667990000025
Se(x，y)＝Se(x1，y1)+Se(x2，y2)+Se(x3，y3)
其中Se(x1，y1)表示动态区域，Se(x2，y2)表示静态区域，Se(x3，y3)表示过渡区域。
5.根据权利要求1所述的一种自适应增强的激光散斑衬比血流成像方法，其特征在于，在S400中，根据所述原始散斑图像序列和分割图像，计算得到自适应增强的散斑衬比图像的方法为：
S401，设定选择的滑动空间窗口大小为W×W，根据所述原始散斑图像序列I(x，y，z)和所述分割图像Se(x，y)，滑动空间窗口遍历像素点坐标(x，y)的取值范围：x∈[1，M]，y∈[1，N]；
S402，选取滑动窗口下对应的原始散斑图像子帧数据IW(x，y，z)和分割图像子帧数据SeW(x，y)，其中x的取值范围为[1，W]，y的取值范围为[1，W]，z的取值范围为
Figure FSA0000296667990000031
S403，计算SeW(x，y)中像素点(x，y)为0的个数C0、像素点为1的个数C1、像素点为0.5的个数C0.5；若C0≠W2或C1≠W2或C0.5≠W2，则说明所述滑动窗口下对应的原始散斑图像子帧数据IW(x，y，z)处于边界区域，使用自适应窗口大小计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka(x，y，z)，即若SeW(x，y)中心点的值为1，则使用SeW(x，y)中等于1所对应的像素点IW1(x，y，z)的值计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka1(x，y，z)；同理，若SeW(x，y)中心点的值为0，则使用SeW(x，y)中等于0所对应的像素点IW0(x，y，z)的值计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka0(x，y，z)；若SeW(x，y)中心点的值为0.5，则使用SeW(x，y)中等于0.5所对应的像素点IW0.5(x，y，z)的值计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka0.5(x，y，z)；
若
Figure FSA0000296667990000032
则Ka(x，y，z)＝Ka1(x，y，z)
Figure FSA0000296667990000033
Figure FSA0000296667990000034
(1)(2)式中，i，j需满足SeW(i，j)＝1；
Figure FSA0000296667990000035
则Ka(x，y，z)＝Ka0(x，y，z)
Figure FSA0000296667990000041
Figure FSA0000296667990000042
(3)(4)式中，i，j需满足SeW(i，j)＝0；
若
Figure FSA0000296667990000043
则Ka(x，y，z)＝Ka0.5(x，y，z)
Figure FSA0000296667990000044
Figure FSA0000296667990000045
(5)(6)式中，i，j需满足SeW(i，j)＝0.5；
其中
Figure FSA0000296667990000046
表示对
Figure FSA0000296667990000047
四舍五入取整，z代表第z帧；
S404，若C0＝W2或C1＝W2或C0.5＝W2，则表示所述滑动窗口下对应的原始散斑图像子帧数据IW(x，y，z)不处于边界区域，采用方差判据计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka(x，y，z)；
其中，在S404中，采用方差判据计算该窗口中心点对应所述原始散斑图像(x，y，z)处的衬比值Ka(x，y，z)的方法为：
S4041，计算IW(x，y，z)的窗口方差dall和4个方向(0°、45°、90°、135°)上的方差d0°，d45°，d90°，d135°，构成dn(z)：
dn(z)＝{d0°，d45°，d90°，d135°，dall}
Figure FSA0000296667990000048
Figure FSA0000296667990000051
Ij(z)表示第z帧对应方向上的像素光强值，
Figure FSA0000296667990000052
表示第z帧对应方向上像素光强的均值，
Figure FSA0000296667990000053
表示整个窗口的像素光强均值；
S4042，若C1＝W2说明该区域属于动态区域，采用方差最小方向上的像素点的值计算Ka(x，y，z)：
ddynamic(z)＝argdmin[var(dn(z))]
Figure FSA0000296667990000054
ddynamic(z)表示动态区域方向，也即方差最小方向，σ(ddynamic(z))表示方差最小方向上的像素点的值的标准差，
Figure FSA0000296667990000055
表示方差最小方向上的像素点的值的均值；
S4043，若C0＝W2说明该区域属于静态区域，采用方差最大方向上的像素点的值计算Ka(x，y，z)：
dstatic(z)＝argdmax[var(dn(z))]
Figure FSA0000296667990000056
dstatic(z)表示静态区域方向，也即方差最大方向，σ(dstatic(z))表示方差最大方向上的像素点的值的标准差，
Figure FSA0000296667990000057
表示方差最大方向上的像素点的值的均值；
S4044，若C0.5＝W2说明该区域属于过渡区域，采用方差中值方向上的像素点的值计算Ka(x，y，z)：
dtransition(z)＝argdmediam[var(dn(z))]
Figure FSA0000296667990000058
dtransition(z)表示过渡区域方向，也即方差中值方向，σ(dtransition(z))表示方差中值方向上的像素点的值的标准差，
Figure FSA0000296667990000059
表示方差中值方向上的像素点的值的均值；
S405，滑动空间窗口遍历像素点坐标(x，y)的取值范围：x∈[1，M]，y∈[1，N]，即循环S402、S403、S404、S4041、S4042、S4043、S4044可得到新的散斑衬比图像Ka(x，y，z)。
CN202310190366.9A 2023-03-02 2023-03-02 一种自适应增强的激光散斑衬比血流成像方法 Pending CN116138760A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310190366.9A CN116138760A (zh) 2023-03-02 2023-03-02 一种自适应增强的激光散斑衬比血流成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310190366.9A CN116138760A (zh) 2023-03-02 2023-03-02 一种自适应增强的激光散斑衬比血流成像方法

Publications (1)

Publication Number Publication Date
CN116138760A true CN116138760A (zh) 2023-05-23

Family

ID=86356220

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310190366.9A Pending CN116138760A (zh) 2023-03-02 2023-03-02 一种自适应增强的激光散斑衬比血流成像方法

Country Status (1)

Country Link
CN (1) CN116138760A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117814731A (zh) * 2024-01-11 2024-04-05 华中科技大学同济医学院附属同济医院 一种血流成像方法、系统、设备和存储介质

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117814731A (zh) * 2024-01-11 2024-04-05 华中科技大学同济医学院附属同济医院 一种血流成像方法、系统、设备和存储介质

Similar Documents

Publication Publication Date Title
CN107361791B (zh) 一种快速超分辨血流成像方法
JP5528083B2 (ja) 画像生成装置、画像生成方法、及び、プログラム
US9330461B2 (en) Image-based method for measuring elasticity of biological tissues and system thereof
Schwarz et al. Motion correction in optoacoustic mesoscopy
EP2237725B1 (en) Therapy assessment with ultrasonic contrast agents
Rajendran et al. Photoacoustic imaging aided with deep learning: a review
Cheng et al. High-resolution photoacoustic microscopy with deep penetration through learning
WO2007084771A2 (en) Quantitative optoacoustic tomography with enhanced contrast
CN107862724B (zh) 一种改进的微血管血流成像方法
CN114209278B (zh) 一种基于光学相干层析成像的深度学习皮肤病诊断系统
Huang et al. Improved ultrafast power Doppler imaging by using spatiotemporal non-local means filtering
Zhao et al. Multiscale vascular enhancement filter applied to in vivo morphologic and functional photoacoustic imaging of rat ocular vasculature
JP2022524163A (ja) 光音響データを解析するためのデバイスおよび方法、光音響システムならびにコンピュータプログラム
JP2022549669A (ja) 時空間データに基づき医用画像を解析するシステム及び方法
CN111436909B (zh) 一种活体组织的光学相干层析成像系统及方法
CN116138760A (zh) 一种自适应增强的激光散斑衬比血流成像方法
JP6946307B2 (ja) 情報取得装置および信号処理方法
Yang et al. Recent advances in deep-learning-enhanced photoacoustic imaging
WO2023039353A2 (en) Real-time super-resolution ultrasound microvessel imaging and velocimetry
US10548477B2 (en) Photoacoustic apparatus, information processing method, and storage medium
US10578588B2 (en) Photoacoustic apparatus, information processing method, and storage medium
JP6486056B2 (ja) 光音響装置および光音響装置の処理方法
CN110522438B (zh) 计算血流速度的方法、装置、介质及血流成像方法和系统
Han et al. A study on an improved laser speckle contrast blood flow imaging methodology
Chizari et al. Mitigation of Motion Artifacts in Handheld Laser Speckle Contrast Imaging Illustrated on Psoriasis Lesions

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