CN119579725B - 一种相干光宽场照明下单帧散斑图像重建方法 - Google Patents
一种相干光宽场照明下单帧散斑图像重建方法Info
- Publication number
- CN119579725B CN119579725B CN202411733896.4A CN202411733896A CN119579725B CN 119579725 B CN119579725 B CN 119579725B CN 202411733896 A CN202411733896 A CN 202411733896A CN 119579725 B CN119579725 B CN 119579725B
- Authority
- CN
- China
- Prior art keywords
- moment
- scattering component
- image
- calculating
- multiple scattering
- 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
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; 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 OR CALCULATING; 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/10068—Endoscopic image
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Medical Informatics (AREA)
- Radiology & Medical Imaging (AREA)
- Quality & Reliability (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
- Image Processing (AREA)
Abstract
本发明涉及一种相干光宽场照明下单帧散斑图像重建方法,方法包括以下步骤:对成像目标进行照明;收集成像目标背向散射光,得到散斑图像,并记录曝光时间;计算一阶矩、二阶矩、三阶矩和四阶矩;计算单次散射分量的一阶矩和二阶矩以及多次散射分量的一阶矩和二阶矩,得到对应的图像,若不要求实时成像,可以采用极大似然估计的方法得到对应的图像。与现有技术相比,本发明具有减少多帧重建造成的伪迹影响,提高重建精度等优点。
Description
技术领域
本发明涉及光学成像领域,尤其是涉及一种相干光宽场照明下单帧散斑图像重建方法。
背景技术
光学成像技术中,显微镜成像、内窥镜成像后都需要进行图像的重建,现有技术中,论文Miao,P.,Zhang,Y.,Wang,C.,&Tong,S.(2022).Random matrix description ofdynamically backscattered coherent waves propagating in a wide-field-illuminated random medium.Applied Physics Letters,120(4).建立了一种在相干光宽场照明下使用多张图像以生成一张单次散射一阶矩图像、一张单次散射二阶矩图像、一张多次散射一阶矩图像、一张多次散射二阶矩图像的方法,并应用于脑皮层成像(Zhang,Y.,Wang,C.,Tong,S.,&Miao,P.(2022).Separating single-and multiple-scatteringcomponents in laser speckle contrast imaging of tissue blood flow.BiomedicalOptics Express,13(5),2881-2895.)、内窥镜成像(Guo,Y.,Weng,Y.,Zhang,Y.,Tong,S.,Liu,Y.,Lu,Z.,&Miao,P.(2023).Random matrix-based laser speckle contrastimaging enables quasi-3D blood flow imaging in laparoscopicsurgery.Biomedical Optics Express,14(4),1480-1493.),同时拓展到3D内窥镜中(ZL202111142881.7)上述方法需要相机对同一位置连续拍照多张图像以生成一张单次散射一阶矩图像、一张单次散射二阶矩图像、一张多次散射一阶矩图像、一张多次散射二阶矩图像进行图像重建,但是这个过程费时费力;在临床应用中由于器官自主运动、患者呼吸心跳影响以及医生手持内窥镜的抖动效应,连续拍照多帧进行重建的方法往往存在很大的伪迹影响,无法准确反映组织浅层和深层的组织吸收和血流等信息。
发明内容
本发明的目的就是为了减少多帧重建造成的伪迹影响,提高重建精度而提供的一种相干光宽场照明下单帧散斑图像重建方法。
本发明的目的可以通过以下技术方案来实现:
一种相干光宽场照明下单帧散斑图像重建方法,方法包括以下步骤:
S1、对成像目标进行照明;
S2、收集成像目标背向散射光,得到散斑图像,并记录曝光时间τcG;
S3、对散斑图像中每一个区域所有灰度值数据计算一阶矩μ1、二阶矩μ2、三阶矩μ3和四阶矩μ4,若要求实时成像,则执行S4,若不要求实时成像,则执行S17;
S4、判断曝光时间τcG是否满足τcM≤τcG≤0.1τcS,其中,τcS和τcM分别为单次散射分量和多次散射分量散斑退相干时间,若是,则执行S5,反之执行S9;
S5、计算单次散射分量的一阶矩μS;
S6、计算单次散射分量的二阶矩
S7、计算多次散射分量的一阶矩μM;
S8、计算多次散射分量的二阶矩然后执行S16;
S9、计算协方差矩阵W;
S10、并对协方差矩阵W进行特征值分解,得到最小特征值为λmin;
S11、基于S3得到的二阶矩μ2、三阶矩μ3和四阶矩μ4求解无偏参数
S12、计算多次散射分量的二阶矩
S13、计算单次散射分量的二阶矩
S14、计算单次散射分量的一阶矩μS;
S15、计算多次散射分量的一阶矩μM,然后执行S16;
S16、在整个散斑图像上,使用的空间窗进行滑动,对将每个空间窗数据作为一个区域执行S3~S15,得到单次散射分量的一阶矩图像、单次散射分量的二阶矩图像、多次散射分量的一阶矩图像和多次散射分量的二阶矩图像;
S17、采用极大似然估计的方法,得到单次散射分量的一阶矩图像、单次散射分量的二阶矩图像、多次散射分量的一阶矩图像和多次散射分量的二阶矩图像。
进一步地,S5中单次散射分量的一阶矩μS的计算为:
其中,μ2表示二阶矩,μ3表示三阶矩。
进一步地,S6中单次散射分量的二阶矩的计算为:
其中,μ2表示二阶矩,μ3表示三阶矩。
进一步地,S7中多次散射分量的一阶矩μM的计算为:
其中,μ1表示一阶矩。
进一步地,S8中多次散射分量的二阶矩的计算为:
进一步地,无偏参数为:
其中,λmin为最小特征值,a1、a2、a3的取值为-1~1的实数,p1和p2的取值为0.5~3.5之间的实数,μ4表示四阶矩。
进一步地,S12中多次散射分量的二阶矩的计算为:
其中,c=M/N,<y>为-5~5之间的常数,M和N表示图像的区域的行数和列数。
进一步地,S17的具体步骤为:
构建关于单次散射分量的一阶矩、多次散射分量的一阶矩、多次散射分量的二阶矩的代价函数,最小化代价函数,设置初始值,然后得到代价函数最小时的单次散射分量的一阶矩、多次散射分量的一阶矩、多次散射分量的二阶矩和单次散射分量的一阶矩。
进一步地,所述代价函数为:
其中,l表示代价函数,erfc为高斯误差函数,zi表示图像区域内所有灰度值。
进一步地,所述照明采用和相机曝光时间上同步的脉冲式照明或采用连续波模式相干光的宽敞照明。
与现有技术相比,本发明具有以下有益效果:
本申请从单帧图像中单次散射和多次散射光强分量的统计特性入手,基于单帧图像光强的一阶矩、二阶矩、三阶矩、四阶矩等信息实现提取单次散射一阶矩图像、单次散射二阶矩图像、多次散射一阶矩图像、多次散射二阶矩图像。基于单帧图像提取的图像进行重建,减少多帧重建造成的伪迹影响,提高重建精度。
附图说明
图1为本发明的流程图。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施,给出了详细的实施方式和具体的操作过程,但本发明的保护范围不限于下述的实施例。
本发明提出一种相干光宽场照明下单帧散斑图像重建方法,方法包括以下步骤:
S1、对成像目标进行照明;
S2、收集成像目标背向散射光,得到散斑图像,并记录曝光时间τcG;
S3、对散斑图像中每一个区域所有灰度值数据计算一阶矩μ1、二阶矩μ2、三阶矩μ3和四阶矩μ4,若要求实时成像,则执行S4,若不要求实时成像,则执行S17;
S4、判断曝光时间τcG是否满足τcM≤τcG≤0.1τcS,其中,τcS和τcM分别为单次散射分量和多次散射分量散斑退相干时间,若是,则执行S5,反之执行S9;
S5、计算单次散射分量的一阶矩μS;
S6、计算单次散射分量的二阶矩
S7、计算多次散射分量的一阶矩μM;
S8、计算多次散射分量的二阶矩然后执行S16;
S9、计算协方差矩阵W;
S10、并对协方差矩阵W进行特征值分解,得到最小特征值为λmin;
S11、基于S3得到的二阶矩μ2、三阶矩μ3和四阶矩μ4求解无偏参数
S12、计算多次散射分量的二阶矩
S13、计算单次散射分量的二阶矩
S14、计算单次散射分量的一阶矩μS;
S15、计算多次散射分量的一阶矩μM,然后执行S16;
S16、在整个散斑图像上,使用的空间窗进行滑动,对将每个空间窗数据作为一个区域执行S3~S15,得到单次散射分量的一阶矩图像、单次散射分量的二阶矩图像、多次散射分量的一阶矩图像和多次散射分量的二阶矩图像;
S17、采用极大似然估计的方法,得到单次散射分量的一阶矩图像、单次散射分量的二阶矩图像、多次散射分量的一阶矩图像和多次散射分量的二阶矩图像。
其中S1~S16的流程图如图1所示。
本发明的具体步骤为:
1、使用连续波(CW)模式相干光(激光)对成像目标进行宽敞照明
2、使用镜头(Lens)或物景(Objectives)收集成像目标背向散射光。相机通过曝光时间τcG记录背向散射光形成的散斑图像I0(x,y)。设τcS和τcM分别为单次散射分量和多次散射分量散斑退相干时间。
3、对散斑图像中每一个M×N的区域(即M行N列,且M≤N),利用M×N区域的所有灰度值数据计算一阶矩(均值)μ1、二阶矩(方差)μ2、三阶矩μ3(偏度)和四阶矩μ4(峰度)。
4、当τcG满足τmM≤τcG≤0.1τcS执行步骤4~8的计算。当τcG不满足上述条件,执行步骤9~15的计算。
5、利用公式计算单次散射分量的一阶矩。
6、利用公式计算单次散射分量的二阶矩。
7、利用公式计算多次散射分量的一阶矩。
8、利用公式计算多次散射分量的二阶矩。
9、计算协方差矩阵W=(I-μJ)(I-μJ)T,其中I为图像I0中大小为M×N的光强矩阵,J为全1矩阵。
10、对W矩阵进行特征值分解,其最小特征值为λmin。
11、利用二阶矩(方差)μ2、三阶矩μ3(偏度)和四阶矩μ4(峰度)求解无偏参数 其中a1、a2、a3的取值为-1~1的实数。p1和p2的取值为0.5~3.5之间的实数。
12、利用Marchenko–Pastur分布和Tracy-Widom分布计算多次散射分量的二阶矩
其中:c=M/N,<y>为-5~5之间的常数。
13、利用公式计算单次散射分量的二阶矩。
14、利用公式μS=σS计算单次散射分量的一阶矩。
15、利用公式μM=μ1-μs计算多次散射分量的一阶矩。
16、在整个图像I0上,使用M×N的空间窗进行滑动,滑动间隔为1像素。每个M×N空间窗数据(设其中心区域像素位置为x0,y0),重复步骤3~15的计算,计算结果μS放入单次散射分量一阶矩图像(x0,y0)位置处;计算结果放入单次散射分量二阶矩图像(x0,y0)位置处;计算结果μM放入多次散射分量一阶矩图像(x0,y0)位置处;计算结果放入多次散射分量一阶矩图像(x0,y0)位置处。
步骤1的激光照明模式也可使用脉冲式照明,此时要求照明和相机曝光在时间上要同步。
对于M≥N的情况,可将图像矩阵I转置后按照同样的流程进行计算。
如果不要求实时成像,也可使用如下极大似然估计的方法来取代步骤4-15,重建单次和多次散射分量的一阶矩和二阶矩图像。M×N区域的所有灰度值{zi},其中i=1…M×N。
S17的具体步骤为:
(a)构建关于单次散射分量的一阶矩、多次散射分量的一阶矩、多次散射分量的二阶矩的代价函数l(μS,μM,σM)为:
其中erfc为高斯误差函数。
(b)优化方法:为最小化代价函数,可使用梯度下降算法或者遗传算法等方法。
(c)初始值:使用随机值。μM和的初始值也可使用{zi}数据的一阶矩和二阶矩。
(d)利用公式计算单次散射分量的一阶矩。
本发明中,单帧图像的计算结果往往信噪比较差,导致分辨率下降。为提升成像质量,可将M×N区域旋转不同的角度进行数据提取并计算不同角度的计算结果进行叠加平均。
本发明的应用场景包括:
(1)神经外科手术显微镜成像:
神经外科手术显微镜引入连续波模式785nm激光照明,使用单色相机记录785nm背向散射光图像,帧速率高于100帧每秒,曝光时间小于1毫秒。对于相机记录的每一帧,使用步骤3~8计算获得一张单次散射一阶矩图像、一张单次散射二阶矩图像、一张多次散射一阶矩图像、一张多次散射二阶矩图像。当需要实时成像,可使用GPU加速。
(2)内窥镜成像
医疗内窥镜引入脉冲模式842nm激光照明,使用RGBIR相机的近红外通道记录842nm背向散射光图像,帧速率为60帧每秒,曝光时间为16.67毫秒。对于相机记录的每一帧,使用步骤3和9~15计算获得一张单次散射一阶矩图像、一张单次散射二阶矩图像、一张多次散射一阶矩图像、一张多次散射二阶矩图像。当需要实时成像,可使用FPGA加速。
以上详细描述了本发明的较佳具体实施例。应当理解,本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此,凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案,皆应在由权利要求书所确定的保护范围内。
Claims (10)
1.一种相干光宽场照明下单帧散斑图像重建方法,其特征在于,方法包括以下步骤:
S1、对成像目标进行照明;
S2、收集成像目标背向散射光,得到散斑图像,并记录曝光时间τcG;
S3、对散斑图像中每一个区域所有灰度值数据计算一阶矩μ1、二阶矩μ2、三阶矩μ3和四阶矩μ4,若要求实时成像,则执行S4,若不要求实时成像,则执行S17;
S4、判断曝光时间τcG是否满足τcM≤τcG≤0.1τcS,其中,τcS和τcM分别为单次散射分量和多次散射分量散斑退相干时间,若是,则执行S5,反之执行S9;
S5、计算单次散射分量的一阶矩μS;
S6、计算单次散射分量的二阶矩
S7、计算多次散射分量的一阶矩μM;
S8、计算多次散射分量的二阶矩然后执行S16;
S9、计算协方差矩阵W;
S10、并对协方差矩阵W进行特征值分解,得到最小特征值为λmin;
S11、基于S3得到的二阶矩μ2、三阶矩μ3和四阶矩μ4求解无偏参数
S12、计算多次散射分量的二阶矩
S13、计算单次散射分量的二阶矩
S14、计算单次散射分量的一阶矩μS;
S15、计算多次散射分量的一阶矩μM,然后执行S16;
S16、在整个散斑图像上,使用的空间窗进行滑动,对将每个空间窗数据作为一个区域执行S3~S15,得到单次散射分量的一阶矩图像、单次散射分量的二阶矩图像、多次散射分量的一阶矩图像和多次散射分量的二阶矩图像;
S17、采用极大似然估计的方法,得到单次散射分量的一阶矩图像、单次散射分量的二阶矩图像、多次散射分量的一阶矩图像和多次散射分量的二阶矩图像。
2.根据权利要求1所述的一种相干光宽场照明下单帧散斑图像重建方法,其特征在于,S5中单次散射分量的一阶矩μS的计算为:
其中,μ2表示二阶矩,μ3表示三阶矩。
3.根据权利要求2所述的一种相干光宽场照明下单帧散斑图像重建方法,其特征在于,S6中单次散射分量的二阶矩的计算为:
其中,μ2表示二阶矩,μ3表示三阶矩。
4.根据权利要求3所述的一种相干光宽场照明下单帧散斑图像重建方法,其特征在于,S7中多次散射分量的一阶矩μM的计算为:
其中,μ1表示一阶矩。
5.根据权利要求4所述的一种相干光宽场照明下单帧散斑图像重建方法,其特征在于,S8中多次散射分量的二阶矩的计算为:
6.根据权利要求1所述的一种相干光宽场照明下单帧散斑图像重建方法,其特征在于,无偏参数为:
其中,λmin为最小特征值,a1、a2、a3的取值为-1~1的实数,p1和p2的取值为0.5~3.5之间的实数,μ4表示四阶矩。
7.根据权利要求6所述的一种相干光宽场照明下单帧散斑图像重建方法,其特征在于,S12中多次散射分量的二阶矩的计算为:
其中,c=M/N,<y>为-5~5之间的常数,M和N表示图像的区域的行数和列数。
8.根据权利要求1所述的一种相干光宽场照明下单帧散斑图像重建方法,其特征在于,S17的具体步骤为:
构建关于单次散射分量的一阶矩、多次散射分量的一阶矩、多次散射分量的二阶矩的代价函数,最小化代价函数,设置初始值,然后得到代价函数最小时的单次散射分量的一阶矩、多次散射分量的一阶矩、多次散射分量的二阶矩和单次散射分量的一阶矩。
9.根据权利要求8所述的一种相干光宽场照明下单帧散斑图像重建方法,其特征在于,所述代价函数为:
其中,l表示代价函数,erfc为高斯误差函数,zi表示图像区域内所有灰度值。
10.根据权利要求1所述的一种相干光宽场照明下单帧散斑图像重建方法,其特征在于,所述照明采用和相机曝光时间上同步的脉冲式照明或采用连续波模式相干光的宽敞照明。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202411733896.4A CN119579725B (zh) | 2024-11-29 | 2024-11-29 | 一种相干光宽场照明下单帧散斑图像重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202411733896.4A CN119579725B (zh) | 2024-11-29 | 2024-11-29 | 一种相干光宽场照明下单帧散斑图像重建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN119579725A CN119579725A (zh) | 2025-03-07 |
| CN119579725B true CN119579725B (zh) | 2025-10-17 |
Family
ID=94815165
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202411733896.4A Active CN119579725B (zh) | 2024-11-29 | 2024-11-29 | 一种相干光宽场照明下单帧散斑图像重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN119579725B (zh) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115705632A (zh) * | 2021-08-02 | 2023-02-17 | 深圳大学 | 白光照明的散射成像方法、装置、计算机设备和存储介质 |
| CN117974525A (zh) * | 2024-03-29 | 2024-05-03 | 华侨大学 | 基于二阶自相关函数计算的激光散斑衬比血流成像方法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20220039679A1 (en) * | 2018-11-27 | 2022-02-10 | ContinUse Biometrics Ltd. | System and method for remote monitoring of biomedical parameters |
| CN113729593B (zh) * | 2021-09-28 | 2022-11-01 | 上海交通大学 | 基于多角度散射随机矩阵的3d内窥镜用血流成像方法 |
| CN113729672A (zh) * | 2021-09-28 | 2021-12-03 | 上海交通大学 | 一种手术显微镜用术中血流成像方法、装置 |
| CN114587635B (zh) * | 2022-03-16 | 2025-11-25 | 苗鹏 | 一种基于实时血流成像的显微手术术中组织灌注自动重建方法 |
-
2024
- 2024-11-29 CN CN202411733896.4A patent/CN119579725B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115705632A (zh) * | 2021-08-02 | 2023-02-17 | 深圳大学 | 白光照明的散射成像方法、装置、计算机设备和存储介质 |
| CN117974525A (zh) * | 2024-03-29 | 2024-05-03 | 华侨大学 | 基于二阶自相关函数计算的激光散斑衬比血流成像方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN119579725A (zh) | 2025-03-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11170258B2 (en) | Reducing noise in an image | |
| US11141044B2 (en) | Method and apparatus for estimating the value of a physical parameter in a biological tissue | |
| JP6521845B2 (ja) | 心拍に連動する周期的変動の計測装置及び計測方法 | |
| JP2021532891A (ja) | マルチスペクトル情報を用いた観血的治療における拡張画像化のための方法およびシステム | |
| JP2021532881A (ja) | マルチスペクトル情報を用いた拡張画像化のための方法およびシステム | |
| US20150126875A1 (en) | Method and system for screening of atrial fibrillation | |
| JP2003529432A (ja) | イメージセンサおよびそれを使用した内視鏡 | |
| US11206991B2 (en) | Systems and methods for processing laser speckle signals | |
| Cheng et al. | Dilated residual learning with skip connections for real-time denoising of laser speckle imaging of blood flow in a log-transformed domain | |
| US20240346632A1 (en) | Medical imaging | |
| WO2021163603A1 (en) | Systems and methods for processing laser speckle signals | |
| CN114587635B (zh) | 一种基于实时血流成像的显微手术术中组织灌注自动重建方法 | |
| Duncan et al. | Spatio-temporal algorithms for processing laser speckle imaging data | |
| CN119579725B (zh) | 一种相干光宽场照明下单帧散斑图像重建方法 | |
| CN113367675A (zh) | 基于激光散斑成像的血流动态检测方法、系统和介质 | |
| CN113729593B (zh) | 基于多角度散射随机矩阵的3d内窥镜用血流成像方法 | |
| CN115063595B (zh) | 一种散斑噪声图像处理装置 | |
| Cheng et al. | Motion-robust anterior–posterior imaging ballistocardiography for non-contact heart rate measurements | |
| CN115908190A (zh) | 一种用于视频图像画质增强的方法及系统 | |
| CN113729672A (zh) | 一种手术显微镜用术中血流成像方法、装置 | |
| CN117315068A (zh) | 基于三维卷积神经网络的无散斑光学相干层析成像方法 | |
| CN118052757B (zh) | 用于观测心肌血流灌注规律的方法及装置 | |
| CN119136732A (zh) | 生物信息测量装置、生物信息测量方法以及生物信息测量程序 | |
| CN119405262B (zh) | 术中监测准周期运动心肌血流的方法及装置 | |
| Akaeva et al. | Methods of digital processing of medical optical images to improve image quality and disease diagnostics accuracy: review |
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 |