CN115067911A - 一种基于gpu实时处理的octa图像优化方法与装置 - Google Patents

一种基于gpu实时处理的octa图像优化方法与装置 Download PDF

Info

Publication number
CN115067911A
CN115067911A CN202210600548.4A CN202210600548A CN115067911A CN 115067911 A CN115067911 A CN 115067911A CN 202210600548 A CN202210600548 A CN 202210600548A CN 115067911 A CN115067911 A CN 115067911A
Authority
CN
China
Prior art keywords
signal
blood flow
octa
classifier
signals
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
CN202210600548.4A
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.)
Zhejiang University ZJU
Jiaxing Research Institute of Zhejiang University
Original Assignee
Zhejiang University ZJU
Jiaxing Research Institute of Zhejiang 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 Zhejiang University ZJU, Jiaxing Research Institute of Zhejiang University filed Critical Zhejiang University ZJU
Priority to CN202210600548.4A priority Critical patent/CN115067911A/zh
Publication of CN115067911A publication Critical patent/CN115067911A/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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/404—Angiography
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/412—Dynamic
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/424—Iterative
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/428—Real-time

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Biophysics (AREA)
  • General Health & Medical Sciences (AREA)
  • Physiology (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Pathology (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Cardiology (AREA)
  • Hematology (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

本发明公开了一种基于GPU实时处理的OCTA图像优化方法与装置。光学相干探测方法，用于收集三维空间内散射信号样本的OCT散射信号；实时血流优化成像方法，基于GPU构建OCT散射信号的信噪比倒数与去相关系数之间的二维特征空间的优化分类器，利用优化分类器对动态血流信号与静态组织信号进行实时优化分类，进而优化获得OCTA图像；通过构建基于图形处理单元(GPU)加速计算的数据处理框架，实现OCTA数据高速处理与实时成像；并实时计算OCTA图像质量；通过当前图像质量反馈调整数据处理与成像参数，实时获得质量最优的OCTA图像。本发明通过GPU加速计算，实现OCTA数据的实时处理和成像，以及OCTA血流投影图像质量的实时自优化。

Description

一种基于GPU实时处理的OCTA图像优化方法与装置
技术领域
本发明大体涉及生物医学成像领域，且更具体地涉及与光学相干层析成像技术(Optical Coherence Tomography，OCT)和血流成像(OCT Angiography，OCTA)相关联的实时图像处理与图像质量优化方法。
背景技术
血液动力学响应是衡量生理功能和病理状态的重要参数，而光学相干层析血流成像(OCTA)技术可以内源性的血流运动代替传统的外源荧光标记物，对生物组织内的微血管网络进行非侵入性、高分辨率的无标记三维成像，可作为观察血液动力学响应的有效手段。在临床实验中，实时的OCTA成像与高质量的血流造影图是分析血液动力学响应的可靠前提。
为了获取OCTA血流图像，通常需要在生物组织的每个空间位置、以一定的时间间隔进行重复采样(重复的A线扫描或B帧扫描)，每个信号处的运动强度通过分析OCT散射信号的时间动态来进行量化，根据量化得到的运动强度来对血流信号和静态组织信号进行分类。目前已报道的OCTA的血流分类，主要是基于相邻的A线扫描间(或相邻B扫描帧间)的差分、方差或去相关计算。其中基于去相关计算的OCTA血流分类，由于其对于窗口内的多个信号的统计特性的充分利用，因此理论上分类结果的可靠性更高。同时，由于去相关衡量的是相邻B扫描帧间的相似度，因此受整体光源强度变化的影响小。
但是，去相关对于运动对比度的量化效果，对原始的OCT散射信号的噪声水平具有显著的依赖性。随着信号强度的衰减(如在深层组织区域)，随机性噪声将逐渐占据主要成分，也将产生较大的去相关值，带来去相关伪影。基于去相关运算生成的运动对比度无法区分噪声的随机性和血红细胞的运动导致的去相关，因此信噪比较弱的区域容易被误判为血流信号区域，严重影响血流图像的对比度。已有的基于信噪比倒数-去相关系数(inverseSNR-decorrelation，ID)特征空间的方法以静态信号在信噪比倒数-去相关系数特征空间分布的3σ边界作为固定的分类边界。ID-OCTA算法虽然能去除大部分静态区域，但是将在信噪比倒数-去相关系数特征空间和静态信号重叠的动态信号一同去除，不能达到动态信号的准确分类，生成的OCTA图像质量存在进一步提升空间。同时，在CPU框架下的ID-OCTA方法因其无法对所有OCT信号进行并行处理，处理速度较慢。
发明内容
本发明针对现有技术的不足，提出了一种基于GPU实时处理的OCTA图像优化方法与装置，方法中利用GPU构建处理框架，通过构建图像质量评价模块实时计算生成血流投影图像的质量参数，根据血流投影图质量反馈调节信噪比倒数-去相关系数空间中线性分类器阈值，从而实现OCTA信号实时处理与图像质量优化。本发明能显著加速OCTA信号处理速度，优化血流投影图像质量，提高血流和背景对比度，改善血管连续性。
本发明的目的是通过如下技术方案实现的：
一、一种基于GPU实时处理的OCTA图像优化方法：
一种光学相干探测方法，用于收集三维空间内散射信号样本的OCT散射信号；
一种实时血流优化成像方法，基于GPU构建OCT散射信号的信噪比倒数与去相关系数之间的二维特征空间的优化分类器，利用优化分类器对动态血流信号与静态组织信号进行实时优化分类，进而优化获得OCTA图像。
所述的一种光学相干探测方法是采用以下的一种：
通过扫描改变参考臂光程的时间域OCT成像方法；
或者利用光谱仪记录光谱干涉信号的光谱域OCT成像方法；
利用探测器记录光谱干涉信号的扫频OCT成像方法。
所述的一种实时血流优化成像方法中，使用GPU实时优化步骤具体包括：
S1、采用一阶自协方差和零阶自协方差分析对OCT散射信号计算分析得到各OCT散射信号的信噪比倒数和去相关系数的两个特征(21)；
S2、构建OCT散射信号的信噪比倒数-去相关系数特征空间(22)，基于多元时间序列模型和莱伊达准则在信噪比倒数-去相关系数特征空间中确定静态组织信号分布的上界和动态血流信号分布的下界，根据上界和下界将信噪比倒数-去相关系数特征空间划分成三个区域，分别为动态信号为主的动态区域、静态信号为主的静态区域和动静态信号混合的过渡区域(23)；
S3、在过渡区域中寻找优化分类线，使得该分类线提取的动态信号生成的血流造影图质量最高。根据优化分类线再从过渡区域进一步分离出信号，再与动态区域的信号结合得到动态OCTA信号，利用动态OCTA信号生成优化的血流造影图(24)。
所述步骤S1具体为：计算各OCT散射信号的信噪比倒数，对同一空间位置T个不同时间点扫描得到的OCT散射信号的幅度部分或对包含幅度和相位的OCT散射信号计算去相关系数，得到各OCT散射信号的信噪比倒数和去相关系数两个特征。
所述步骤S2具体为：根据多元时间序列模型获得OCT散射信号在信噪比倒数-去相关系数特征空间的渐近分布，借助仿真分析得到渐近分布的方差，确定静态信号分布的上界及动态信号分布的下界。
所述多元时间序列模型的具体定义为：
采集到的OCT散射信号X(s，t)由生物组织样本的后向散射光信号A(s，t)和复数高斯白噪声n(s，t)叠加而成，其中后向散射光信号为A(s，t)＝αAd(s，t)+(1-α)As(s，t)，Ad(s，t)和As(s，t)分别为后向散射光信号中的动态和静态成分，α调整两者的比例，s，t分别表示体素的三维空间和时间。一般地，后向散射光信号A和噪声n满足以下条件：1、n(s，t)是时间和空间维度相互独立的白噪声，并且和后向散射光A(s，t)无关；2、在时间维度，后向散射光的动态部分Ad(s，t)遵循独立同分布，静态部分As(s，t)为一个常矢量。
通过多元时间序列模型分析，对于运动引起动态程度为α的OCT散射信号，其在信噪比倒数-去相关系数特征空间的渐近分布满足：
Figure BDA0003669090620000031
其中，→表示收敛，iSNR表示信噪比倒数，D表示去相关系数，α表示OCT散射信号的动态程度。
所述借助仿真分析得到渐近分布的方差，具体为：当计算去相关系数的平均核大小有限时，进一步通过仿真模拟得到对于当α＝0下的静态信号在信噪比倒数-去相关系数特征空间的分布方差σs，满足下式：
Figure BDA0003669090620000032
其中，G为仿真分析得到的常数，N为计算去相关系数的平均核大小，N＝S×T，S、T分别表示体素的三维空间尺寸和时间尺寸。
所述确定静态信号分布的上界及动态信号分布的下界，具体为：按照以下公式获得静态信号分布的上界Ds1和动态信号分布的下界Dd2：
Figure BDA0003669090620000033
Figure BDA0003669090620000041
式中，α0为显著动态的下限。
所述S2中，在信噪比倒数-去相关系数特征空间中，将去相关值大于静态信号分布的上界Ds1的区域作为以动态信号为主的动态区域Rd，将介于动态信号分布的下界Dd2和静态信号分布的上界Ds1之间的区域作为动静态信号混合的过渡区域Ri，将去相关值小于动态信号分布的下界Dd2的区域作为以静态信号为主的静态区域Rs。
所述步骤S3具体为：
S31、根据两个分类器阈值分别确定两条分界线，对信噪比倒数分别小于两个分类器阈值对应的两条分界线的信号进行提取获得两种动态OCTA信号，利用两种动态OCTA信号各自生成一幅血流造影图，获得两幅血流造影图(24)；
S32、计算两幅血流造影图(24)的质量评价参数AQI；
S33、比较两幅血流造影图(24)的质量评价参数AQI大小迭代更新分类器阈值，利用更新分类器阈值回到步骤S31中重新获得两幅血流造影图，再进行步骤S32中重新获得两幅血流造影图(24)的质量评价参数AQI，由此不断重复步骤迭代直至质量评价参数AQI最大，则取质量评价参数AQI最大时所对应的两个分类器阈值；
S34、利用质量评价参数AQI最大时所对应的两个分类器阈值求平均值，以平均值确定的分界线对过渡区域进行分割，在过渡区域中将去相关值大于分界线的区域作为分离出的信号，将分离出的信号和动态区域的信号叠加得到动态OCTA信号。
上述过程中采用GPU加速计算，从而实现血流造影图质量实时优化。
所述步骤S31中，以步骤S2获得的静态组织信号分布的上界和动态血流信号分布的下界作为初始的两条分界线，在信噪比倒数-去相关系数特征空间中以两条初始分界线和去相关系数坐标轴之间的夹角角度为初始的两个分类器阈值。
对信噪比倒数分别小于初始的两个分类器阈值对应的初始的两条分界线的信号进行提取，分别获得两种动态OCTA信号，利用两种动态OCTA信号各自生成一幅初始血流造影图，获得两幅初始血流造影图。
所述步骤S32中，使用GPU实时计算两幅血流造影图(24)的血管密度和图像对比度的参数，将两者的乘积作为血流投影图的质量评价参数AQI。
所述计算血流造影图的血管密度及图像对比度参数，具体为：
对二维血流投影图A进行分析，记血流投影图的边长为L，在血流投影图行列
Figure BDA0003669090620000051
和
Figure BDA0003669090620000052
处分别选取两两相互垂直的4条直线，计算每条直线上各相邻像素的灰度强度差值，组成像素强度梯度曲线，记4条直线中第i条直线的强度梯度曲线与横坐标轴(像素位置)的交点个数为Vi，得到血流投影图的血管密度VD：
Figure BDA0003669090620000053
二维血流投影图像在血管位置处呈现锐利的灰度强度变化，这样变化的频率被所提出的血管密度VD捕捉。
血流投影图像的图像对比度用均方根对比度CRMS来表示：
Figure BDA0003669090620000054
其中，NA代表血流投影图A中总的像素数目，I(x，y)代表血流投影图A中坐标在(x，y)处的像素灰度强度，
Figure BDA0003669090620000055
代表血流投影图A的平均灰度强度。
最后由血管密度VD和均方根对比度CRMS的乘积得到血流投影图A的质量评价参数AQI：
AQI＝VD×CRMS。
所述步骤S33中，迭代更新分类器阈值，具体为：
两个分类器阈值分别为第一分类器阈值γ1与第一分类器阈值γ2，第一分类器阈值γ1与第一分类器阈值γ2对应获得的血流投影图的质量评价参数AQI分别为第一质量评价参数AQI1与第二质量评价参数AQI2，并进行以下判断：
若第一质量评价参数AQI1<第二质量评价参数AQI2，则下次迭代的第一分类器阈值γ1与第一分类器阈值γ2分别设为
Figure BDA0003669090620000056
与AQI2；否则下次迭代的第一分类器阈值γ1与第一分类器阈值γ2分别设为AQI1与
Figure BDA0003669090620000057
所述步骤S33中，某次迭代结束后第一分类器阈值γ1与第一分类器阈值γ2之差小于1，认为此时质量评价参数AQI最大。
二、基于GPU实时处理的光学相干血流造影装置：
一个OCT光学相干层析探测装置，用于对三维空间内的散射信号样本的OCT散射信号进行采集；
一个或多个基于GPU加速的信号处理器，用于获取并分析OCT散射信号的信噪比倒数和去相关系数，从而构建信噪比倒数与去相关系数之间的二维特征空间的优化分类器，利用优化分类器对动态血流信号与静态组织信号进行实时优化分类，进而优化获得OCTA图像。
具体实施中内部还构建流投影图质量评价模块，用于信噪比倒数-去相关系数特征空间中动态信号线性分类器阈值的反馈调节和更新。
所述的一OCT光学相干层析装置是采用以下的一种：
包括低相干光源、干涉仪和探测器；
或者包括低相干光源、干涉仪和光谱仪；
或者包括扫频宽光谱光源、干涉仪和探测器。
所述的一个OCT光学相干层析探测装置中选择地配置一个可见光指示装置，用于指示OCT探测光束的位置，指导探测目标的放置位置。
所述可见光指示装置主要由可见光指示光源和准直透镜组成。
本发明通过构建基于图形处理单元(GPU)加速计算的数据处理框架，实现OCTA数据高速处理与实时成像；并实时计算OCTA图像质量；通过当前图像质量反馈调整数据处理与成像参数，实时获得质量最优的OCTA图像。
本发明通过GPU加速计算，实现OCTA数据的实时处理和成像，以及OCTA血流投影图像质量的实时自优化。
本发明的有益效果和创新点如下：
本发明通过计算OCT散射信号的信噪比倒数以及去相关系数，构成二维信噪比倒数-去相关系数特征空间，结合时间序列模型，将信噪比倒数-去相关系数特征空间分为动态区域，过渡区域及静态区域。设计血流投影图像质量评价模块，基于血管密度及图像对比度反馈优化调整信噪比倒数-去相关系数特征空间中的线性分类器阈值参数，从而提取动态血流信号。构建GPU计算框架，用于加速过程中所有信号处理流程，从而实现OCTA信号的实时处理与血流投影图像质量优化。相比于现有方法，本发明能显著加速OCTA信号处理速度，提升动态信号分类的准确率，优化OCTA图像质量。
本发明的图像质量优化方法为数据驱动式的，从数据本身的血管密度和对比度出发，对于每组数据均能自适应地优化至其最佳的图像质量。
本发明对比已有技术具有以下显著优点：
1.基于去相关计算的OCTA，由于OCT散射信号的去相关系数和信号强度间的依赖关系，低信噪比区域随机噪声引入的去相关伪影无法和血流运动引入的去相关区分。常用解决办法是设置经验性分类器阈值进行强度掩膜，相当于在强度-去相关系数特征空间中对信号用去相关分类器阈值区分后，再用一个信噪比倒数分类器阈值进行区分，而信号的去相关系数和信噪比间更复杂的依赖关系使得实际血流信号和其它信号的分界与强度分类器阈值直线有很大差异，造成误分类。但本发明提出的分类器基于对信噪比倒数-去相关系数特征空间中静态信号的时间序列模型分析，具有信噪比自适应的优势，此外通过血流投影图像质量评价模块，反馈调整线性分类器阈值，从而自优化分类效果。
2.本发明提出的基于血流投影图像质量评价的图像质量优化方法，能自适应地调整信噪比倒数-去相关系数特征空间中线性分类器的分类器阈值参数，从而有效地对过渡区域中的动静态信号进行细分类，在抑制静态信号的同时，有效地提取出了更充分的动态信号。
3.对比现有的方法，本发明构建的基于GPU加速的计算框架，可以实现实时的OCTA信号处理与图像质量优化，其处理速度表现显著优于其他方法；对信噪比倒数-去相关系数特征空间线性分类器进行基于图像质量的自适应优化方法，经过大量样本验证，在图像信噪比、对比噪声比、血管连续性参数的表现上显著优于传统方法。
附图说明
图1为本发明方法的示意图；
图2为本发明装置的示意图；
图3为本发明实施例的装置的示意图；
图4为本发明实施例的基于GPU加速的信号处理框架流程图；
图5为本发明示例性实施例的图像质量评价参数反馈调整分类分类器阈值结果图；
图6为本发明示例性实施例的活体小鼠视网膜血流成像实验结果图；
图7为本发明对小鼠视网膜造影结果与传统方法结果的分析对比图。
图中：1-一种光学相干探测方法；2-一种实时血流优化成像方法；21-计算OCT信号的信噪比倒数和去相关系数两个特征；22-构建OCT散射信号的信噪比倒数-去相关系数信噪比倒数-去相关系数特征空间；23-计算静态信号分布的上下界，将信噪比倒数-去相关系数特征空间分为动态区域，过渡区域，静态区域三部分；24-在过渡区域寻找最优分类线提取动态信号，生成血流造影图。
具体实施方式
下面将结合附图对本发明的具体实施方式做详细说明，附图形成本发明的一部分。需要注意的是，这些说明及示例仅仅为示例性，不能被理解为限制了本发明的范围，本发明的保护范围由随附的权利要求书限定，任何在本发明权利要求基础上的改动都是本发明的保护范围。
本发明的实施例如下：
为了便于理解本发明的实施例，将各操作描述成多个离散的操作，但是描述的顺序不代表实施操作的顺序。
本描述中针对样品测量空间采用基于空间方向的x-y-z三维坐标表示。这种描述仅仅用于促进讨论，而不意欲限制本发明的实施例的应用。其中：深度z方向为沿入射光轴的方向；x-y平面为垂直于光轴的平面，其中x与y正交，且x表示OCT横向快扫描方向，y表示慢扫描方向。
上述α，A，N等表示变量，仅仅用于促进讨论，而不意欲限制本发明的实施例的应用，可以是1,2,3等任一数值。
本发明方法如图1所示，信号采集部分，对组织样本进行OCT三维扫描成像，相同或相邻空间位置在T个不同时间点重复采样1。
信号分类部分，结合散射信号强度与时间去相关性构建二维特征空间，实现动态血流信号与静态组织背景的分类2。
其具体步骤是：首先通过采集器对小鼠视网膜样本进行时空采样，对采集到的原始光谱信号进行OCT重构，具体包括插值，色散校正，快速傅里叶变换(FFT)，帧间相关配准，如在图4的流程图中OCT部分所展示的。之后对重建得到的OCT复数信号矩阵进行OCTA处理，具体包括分层，构建信噪比倒数-去相关系数特征空间分类器提取血流信号，如在图4的流程图中OCTA部分所展示的。再设计图像质量评价模块对上一部中生成的OCTA投影图像进行质量评价，根据评价参数反馈调整OCTA部分中信噪比倒数-去相关系数特征空间线性分类器的分类器阈值，如在图4的流程图中血流投影图像质量评价与线性ID分类器阈值优化部分所展示的。最后保存三维OCT/OCTA处理结果，生成自优化的OCTA血流投影图。以上步骤中的信号处理流程均在所搭建的基于GPU加速计算的数据处理框架中进行，以实现实时的OCTA信号处理与血流投影图像质量优化。
构建信噪比倒数-去相关系数特征空间所用的去相关系数计算基于分析血流信号和周围组织的相对运动，得到各OCT散射信号的相关程度信息。去相关运算包括对T个不同时间点扫描得到的OCT散射信号的幅度部分、或对包含幅度和相位的复数OCT散射信号计算去相关系数，其中对复数信号的相关计算理论上可以获得更高的运动对比度，相邻B扫描帧(x-z平面)间的相关系数可以由下式计算：
Figure BDA0003669090620000081
其中，CorrMapi(z，x)表示在(z，x)处计算得到体素的相关系数，i表示第i帧B扫描帧，z表示深度方向，x表示线扫描方向，Ai(，)表示复数信号，*表示复数的共轭，
Figure BDA0003669090620000091
和
Figure BDA0003669090620000092
是相邻两B扫描帧的复数强度。P和Q分别表示进行相关系数计算时的所取窗口的z方向和x方向的大小，这里窗口大小选定为5×1(z×x)像素，p和q分别表示相应的窗口内的像素位置。计算得到的相关系数的范围是0到1，用1减相关系数可将其转化为去相关系数。计算过程中，每个窗口对应一个去相关系数值，移动窗口使其遍历整个z-x平面，可以得到整个平面内所有位置的去相关值。同一扫描位置相邻各帧的去相关值取平均作为结果。
在OCT系统中，噪声来源主要为散粒噪声，视为在整个扫描体积内近似不变，可以通过扫描空气的数据或事先对系统进行标定而确定。
然后采用以下公式计算每个体素的信噪比倒数iSNR，定义如下：
Figure BDA0003669090620000093
其中，s2是OCT系统的噪声水平，I表示每个体素的OCT信号强度值。
通过计算OCT散射信号中数据空白区域的强度或者针对OCT系统事先标定得到噪声水平。至此已完成OCT信号去相关系数与信噪比倒数的计算21。
结合上述信号处理得到的OCT信号信噪比倒数和基于去相关计算得到的OCTA血流信息建立ID二维特征空间，并将OCT散射信号投影在特征空间中22。根据通过多元时间序列模型分析的结果，对于运动引起的OCT散射信号，其在信噪比倒数-去相关系数特征空间的渐近分布满足：
Figure BDA0003669090620000094
其中，→表示收敛，iSNR表示信噪比倒数，α表示OCT散射信号的动态程度，D表示去相关系数。
当平均核大小有限时，进一步通过仿真模拟得到对于静态信号(α＝0)在信噪比倒数-去相关系数特征空间的分布方差σs，满足下式：
Figure BDA0003669090620000095
其中，G≈1.5为仿真分析得到的常数，N＝S×T为高斯平均核大小，S、T分别表示体素的三维空间尺寸和时间尺寸。
之后根据多元时间序列模型及仿真分析，结合经验值在信噪比倒数-去相关系数特征空间定义静态信号分布的上界Ds1和动态信号分布的下界Dd2，定义如下：
Figure BDA0003669090620000101
Figure BDA0003669090620000102
式中，α0为显著动态的下限，实例中其值为0.4。
根据定义的(5)式和(6)式可以将信噪比倒数-去相关系数特征空间分成三个部分23：以静态信号为主的静态区域Rs，以动态信号为主的动态区域Rd，动静态信号混合的过渡区域Rt，如图5(a)中所示。
对信噪比倒数-去相关系数特征空间中的动态信号进行初始分类时，所采用的分类器阈值即为Ds1与Dd2，以保证在之后的分类器阈值迭代优化过程中，线性分类器一直在过渡区域，从而在过渡区域中提取适当的动态信号。以Ds1与Dd2为分类器对过渡区域的混合信号进行分类，提取出其中的动态信号，生成对应的血流造影图像24。
基于血流造影图像的血管密度与图像对比度构建OCTA图像质量评价参数，如图5(b)中各二维血流图所示，对二维血流投影图A进行分析，记血流投影图的边长为L，在血流投影图行列
Figure BDA0003669090620000103
和
Figure BDA0003669090620000104
处分别选取两两相互垂直的4条直线，计算每条直线上各相邻像素的灰度强度差值，组成像素强度梯度曲线，如图5(c)所示，记4条直线中第i条直线的强度梯度曲线与横坐标轴(像素位置)的交点个数为Vi，可得到血流投影图的血管密度VD：
Figure BDA0003669090620000105
二维血流投影图像在血管位置处呈现锐利的灰度强度变化，这样变化的频率被所给出的血管密度VD捕捉，可以反映血流投影图像的血管数量。
血流投影图像的图像对比度用均方根对比度CRMS来表示：
Figure BDA0003669090620000106
其中，NA代表二维血流投影图A中总的像素数目，I(x，y)代表坐标在(x，y)处的像素灰度强度，
Figure BDA0003669090620000107
代表血流投影图的平均灰度强度。
则可由血管密度VD和均方根对比度CRMS的乘积得到血流投影图的质量评价参数AQI：
AQI＝VD×CRMS (9)
图5(b)展示了同一OCT数据经过不同线性ID分类器γ分类生成的血流图像，其图像质量AQI值。在信噪比倒数-去相关系数特征空间的过渡区域中，有某一分类器阈值对应的线性ID分类器可以较完整地提取动态信号，生成图像质量最优的OCTA投影图。同一数据由不同ID分类器分类生成的血流投影图，其血管密度VD，均方根对比度CRMS，血流投影图质量评价参数AQI与其对应分类器的分类器阈值关系如图5(d)所示，血管密度VD随分类器阈值γ的上升，先由于更多的动态信号被提取而上升，随着动态信号被提取完毕，由于静态组织和噪声的影响在其最大值附近波动；而均方根对比度CRMS则随分类器阈值γ的上升呈先上升后下降的趋势。由两者相乘得到的血流投影图质量评价参数AQI也因此随分类器阈值γ的上升呈先上升后下降的趋势。
同一数据条件下，其血流投影图质量评价参数AQI最大时对应的ID分类器阈值可以通过基于的迭代搜索法确定。具体为：记某次迭代中两个起始分类器阈值为γ1与γ2，对应血流投影图的质量评价参数AQI分别为AQI1与AQI2，若AQI1＜AQI2，则下次迭代的起始分类器阈值为
Figure BDA0003669090620000111
与AQI2，否则为AQI1与
Figure BDA0003669090620000112
直至某次迭代结束后的起始分类器阈值之差小于1，此时即找到血流投影图质量评价参数AQI最大时对应的ID分类器阈值，从而实现OCTA图像质量的自优化。
对信噪比倒数-去相关系数特征空间按照基于血流投影图质量评价反馈优化的分类器阈值进行分类及动态信号提取，可再从图形角度进行高斯滤波或Hessian滤波等处理，从而进一步得到血流造影结果。
上述信号处理过程均在如图4所示的基于GPU加速计算的处理流程框架中进行，实例中所用GPU型号为RTX 2080Ti(32GB显存)，达到133kHz的线处理速度，对于线采集速度为120kHz的采集器(1)可以实现实时的OCTA信号处理与图像质量优化。
图2示出的是本发明中基于GPU实时处理的OCT血流成像技术的采集装置结构示意图。该装置的低相干干涉测量部分的主体结构为一干涉仪，由11～23构成，其中光源11发出的光被分束器12分成两部分光束：其中的一束光经一偏振控制器13进入到干涉仪的参考臂，通过参考臂准直镜14照射于平面反射镜15上；另一束光经另一偏振控制器13进入到样品臂，经过准直透镜16和扫描装置光路聚焦到待测样品21上。其中扫描装置光路中，光束经过二维扫描振镜组17、18、“4f”透镜组54、55和二向色镜19的反射后，经过聚焦物镜20聚焦在待测样品21上，透镜组54、55的设计是为了确保扫描时二维扫描振镜镜面的光束中心和二相色镜反射面的光束中心固定不变化，使得OCT样品臂中的光束在扫描时不影响物镜的成像性质。而后参考臂和样品臂各自反射回的光发生干涉后由干涉信号探测装置22接收，干涉信号探测装置22再连接到信号处理器模块与计算单元23。对于光纤型光路，采用偏振控制器13调整光束的偏振态，最大化信号干涉效果。
具体实施还设置有可见光指示装置，可见光指示装置包括低功率可见光光源25，准直透镜24和滤光片52，用于指示的可见光经过准直透镜24、二向色镜19和聚焦物镜20后到待测样品21。
依据低相干干涉探测信号的不同方式，图2所示的一种基于GPU实时处理的OCT血流成像技术的采集装置具体包括：
1)时间域测量装置。光源11采用宽带低相干光，平面反射镜15可沿光轴方向移动，干涉信号探测装置22为一点探测器。通过移动平面反射镜15改变参考臂光程，两臂的干涉信号由点探测器22探测到，对某一空间深度的z方向的散射信号的低相干干涉探测，从而得到深度空间维度的采样体。
2)光谱域测量装置。光源11采用宽带低相干光，平面反射镜15固定不动，干涉信号探测装置22采用光谱仪。干涉信号经过光谱仪中的线阵相机同时记录干涉光谱。采用傅里叶分析方法分析干涉光谱信号，并行获取深度z方向的散射信息，从而得到深度空间维度的采样体。
3)扫频测量装置。光源11采用扫频光源，平面反射镜15固定不动，干涉信号探测装置22采用点探测器。点探测器分时记录扫频光源的低相干干涉光谱。采样傅里叶分析干涉光谱信号，并行获取深度z方向的散射信息，从而得到深度空间维度的采样体。
图3示出的是利用本发明的一个示例性实施例。基于特征空间的三维血流成像系统，包括宽带低相干光源26、光环形器27、分光比为50:50的光纤耦合器28、第一偏振控制器29、第一光纤准直器件30、聚焦透镜36、平面反射镜37、第二偏振控制器38、第二光纤准直器件39、二维扫描振镜组合40和41、二向色镜42、聚焦物镜43、第三光纤准直器件45、光栅46、聚焦透镜47、高速线阵相机48、信号处理器模块与计算单元49、可见光指示光源50、准直透镜51、“4f”透镜组56和57，其中宽带低相干光源26采用中心波长为1325nm、带宽为100nm的超发光二极管光源，聚焦物镜43采用焦距为30mm的消色差双胶合透镜，高速线阵相机48采用由2048像素单元组成的线阵扫描相机；其中由本发明装置所使用的低相干宽带光源26发出的光，经过光环行器27后进入到分光比为50:50的光纤耦合器28，从光纤耦合器28出射的光被分成两部分子光束：其中一束光通过光纤经过第一偏振控制器29连接至参考臂中的第一光纤准直器件30，经过准直和聚焦透镜36后照射到平面反射镜37；另一束光通过光纤经过第二偏振控制器38连接至样品臂部分的第二光纤准直器件39，准直后经过两个扫描振镜40、41、“4f”透镜组56、57和二向色镜42反射后，由聚焦物镜43聚焦到被测样品44上，其中透镜组56、57的设计是为了确保扫描时二维扫描振镜镜面的光束中心和二相色镜反射面的光束中心固定不变化。由参考臂中平面反射镜37反射的光与样品臂中被测样品背向散射的光在光纤耦合器28处干涉，干涉光经过光谱仪45～48探测并被记录，而后由信号处理器模块与计算单元49采集并作信号分析处理。光谱仪包括依次连接的器件45～48，其中器件45为光纤耦合器，46为光栅，47为汇聚透镜，将光栅色散分出的光聚焦到48所示的线阵探测器上。
具体实施还设置有可见光指示装置，可见光指示装置包括可见光指示光源50、准直透镜51，可见光指示光源50发出用于指示的可见光经过准直透镜51、二向色镜42和聚焦物镜43后到待测样品44。
图6为活体小鼠视网膜血流成像实验结果图。本发明提出的一种基于GPU实时处理的OCTA图像优化方法与装置可以准确对视网膜各结构层分层，生成高质量的视网膜各层血流投影图像，包括视网膜全层(RE)，视网膜浅层(SVP)，视网膜中层(ICP)，视网膜深层(DCP)与脉络膜层(CH)。
图7为本发明对小鼠视网膜造影结果与传统方法结果的分析对比图。传统方法选择的是采用静态信号在信噪比倒数-去相关系数特征空间分布的上界作为分类分类器阈值的ID-OCTA方法。本发明提出的基于图像质量评价的分类器阈值自适应优化方法(图7(a))和ID-OCTA(图7(b))都能得到较低的背景噪声，而本发明提取的血流信号更为丰富，相对背景有更高的对比度，并且血管连续性更好，说明本发明能够在抑制过渡区域噪声信号的同时有效地提取出动态信号。而在ID-OCTA中，由于采用静态信号在信噪比倒数-去相关系数特征空间分布的上界作为分类分类器阈值，在保留完全动态信号的同时也将过渡区域中可能的动态信号完全剔除，导致部分血流信号被除去，影响了血管的连续性。图7(c)反映了本发明在信噪比倒数-去相关系数特征空间中比ID-OCTA方法提取了更多的OCTA信号。
进一步地，为了量化基于图像质量评价的分类器阈值自适应优化方法相对ID-OCTA方法带来的图像质量提升，本发明基于图像信噪比SNR，对比噪声比CNR和血管连续性CON来表征两种方法的性能，噪声比CNR由下式计算得到：
Figure BDA0003669090620000141
其中，
Figure BDA0003669090620000142
分别是血管区域和背景噪声区域的去相关系数的均值，σn是背景噪声区域的标准差。
血管连续性CON由下式计算得到：
Figure BDA0003669090620000143
其中，Nc为投影图骨架化后连接点的数目，Ns为投影图骨架化后所有像素点的数目。
表1为对10组数据的小鼠视网膜全层血流投影图的定量评估结果。对SNR,CNR,CON各参数采用平均值±标准差的方式描述，Δ表示参数数值的提升。
表1本发明与ID-OCTA血流分类方法对于血流成像性能的定量比较
Figure BDA0003669090620000144
表1的定量分析结果表明，本发明所提出的基于图像质量评价的分类器阈值自适应优化方法得到的血流造影结果中，生成血流投影图的信噪比，对比噪声比，血管连续性均优于ID-OCTA方法。该对比结果充分说明：利用本发明所涉及的基于GPU实时处理的OCTA图像优化方法，能够提升血流信号分类的准确率，实现血流对比度的有效增强与信号处理速度的提升，具有其突出显著的技术效果。

Claims (10)

1.一种基于GPU实时处理的OCTA图像优化方法，其特征在于包括：
一种光学相干探测方法(1)，用于收集三维空间内散射信号样本的OCT散射信号；
一种实时血流优化成像方法(2)，基于GPU构建OCT散射信号的信噪比倒数与去相关系数之间的二维特征空间的优化分类器，利用优化分类器对动态血流信号与静态组织信号进行实时优化分类，进而优化获得OCTA图像。
2.根据权利要求1所述的一种基于GPU实时处理的OCTA图像优化方法，其特征在于：所述的一种光学相干探测方法(1)是采用以下的一种：
通过扫描改变参考臂光程的时间域OCT成像方法；
或者利用光谱仪记录光谱干涉信号的光谱域OCT成像方法；
利用探测器记录光谱干涉信号的扫频OCT成像方法。
3.根据权利要求1所述的一种基于GPU实时处理的OCTA图像优化方法，其特征在于：所述的一种实时血流优化成像方法(2)中，具体包括：
S1、采用一阶自协方差和零阶自协方差分析对OCT散射信号计算分析得到各OCT散射信号的信噪比倒数和去相关系数的两个特征(21)；
S2、构建OCT散射信号的信噪比倒数-去相关系数特征空间(22)，在信噪比倒数-去相关系数特征空间中确定静态组织信号分布的上界和动态血流信号分布的下界，根据上界和下界将信噪比倒数-去相关系数特征空间划分成三个区域，分别为动态信号为主的动态区域、静态信号为主的静态区域和动静态信号混合的过渡区域(23)；
S3、在过渡区域中寻找优化分类线，根据优化分类线再从过渡区域进一步分离出信号，再与动态区域的信号结合得到动态OCTA信号，利用动态OCTA信号生成优化的血流造影图(24)。
4.根据权利要求3所述的一种基于GPU实时处理的OCTA图像优化方法，其特征在于：所述S2中，在信噪比倒数-去相关系数特征空间中，将去相关值大于静态信号分布的上界Ds1的区域作为以动态信号为主的动态区域Rd，将介于动态信号分布的下界Dd2和静态信号分布的上界Ds1之间的区域作为动静态信号混合的过渡区域Ri，将去相关值小于动态信号分布的下界Dd2的区域作为以静态信号为主的静态区域Rs。
5.根据权利要求3所述的一种基于GPU实时处理的OCTA图像优化方法，其特征在于：所述步骤S3具体为：
S31、根据两个分类器阈值分别确定两条分界线，对信噪比倒数分别小于两条分界线的信号进行提取获得两种动态OCTA信号，利用两种动态OCTA信号各自生成一幅血流造影图，获得两幅血流造影图(24)；
S32、计算两幅血流造影图(24)的质量评价参数AQI；
S33、比较两幅血流造影图(24)的质量评价参数AQI大小迭代更新分类器阈值，利用更新分类器阈值回到步骤S31中重新获得两幅血流造影图，再进行步骤S32中重新获得两幅血流造影图(24)的质量评价参数AQI，由此不断重复步骤迭代直至质量评价参数AQI最大，则取质量评价参数AQI最大时所对应的两个分类器阈值；
S34、利用质量评价参数AQI最大时所对应的两个分类器阈值求平均值，以平均值确定的分界线对过渡区域进行分割，在过渡区域中将去相关值大于分界线的区域作为分离出的信号，将分离出的信号和动态区域的信号叠加得到动态OCTA信号。
6.根据权利要求5所述的一种基于GPU实时处理的OCTA图像优化方法，其特征在于：所述步骤S31中，以步骤S2获得的静态组织信号分布的上界和动态血流信号分布的下界作为初始的两条分界线，在信噪比倒数-去相关系数特征空间中以两条初始分界线和去相关系数坐标轴之间的夹角角度为初始的两个分类器阈值。
7.根据权利要求5所述的一种基于GPU实时处理的OCTA图像优化方法，其特征在于：所述步骤S32中，使用GPU实时计算两幅血流造影图(24)的血管密度和图像对比度的参数，将两者的乘积作为血流投影图的质量评价参数AQI。
8.根据权利要求5所述的一种基于GPU实时处理的OCTA图像优化方法，其特征在于：所述步骤S33中，迭代更新分类器阈值，具体为：
两个分类器阈值分别为第一分类器阈值γ1与第一分类器阈值γ2，第一分类器阈值γ1与第一分类器阈值γ2对应获得的血流投影图的质量评价参数AQI分别为第一质量评价参数AQI1与第二质量评价参数AQI2，并进行以下判断：
若第一质量评价参数AQI1＜第二质量评价参数AQI2，则下次迭代的第一分类器阈值γ1与第一分类器阈值γ2分别设为
Figure FDA0003669090610000021
与AQI2；否则下次迭代的第一分类器阈值γ1与第一分类器阈值γ2分别设为AQI1与
Figure FDA0003669090610000022
9.用于实施权利要求1～8任一所述方法的基于GPU实时处理的光学相干血流造影装置，包括：
一个OCT光学相干层析探测装置，用于对三维空间内的散射信号样本的OCT散射信号进行采集；
一个或多个基于GPU加速的信号处理器，用于获取并分析OCT散射信号的信噪比倒数和去相关系数，从而构建信噪比倒数与去相关系数之间的二维特征空间的优化分类器，利用优化分类器对动态血流信号与静态组织信号进行实时优化分类，进而优化获得OCTA图像。
10.根据权利要求9所述的基于GPU实时处理的光学相干血流造影装置，其特征在于：
所述的一OCT光学相干层析装置是采用以下的一种：
包括低相干光源、干涉仪和探测器；
或者包括低相干光源、干涉仪和光谱仪；
或者包括扫频宽光谱光源、干涉仪和探测器。
所述的一个OCT光学相干层析探测装置中选择地配置一个可见光指示装置，用于指示OCT探测光束的位置，指导探测目标的放置位置。
所述可见光指示装置主要由可见光指示光源和准直透镜组成。
CN202210600548.4A 2022-05-30 2022-05-30 一种基于gpu实时处理的octa图像优化方法与装置 Pending CN115067911A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210600548.4A CN115067911A (zh) 2022-05-30 2022-05-30 一种基于gpu实时处理的octa图像优化方法与装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210600548.4A CN115067911A (zh) 2022-05-30 2022-05-30 一种基于gpu实时处理的octa图像优化方法与装置

Publications (1)

Publication Number Publication Date
CN115067911A true CN115067911A (zh) 2022-09-20

Family

ID=83249811

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210600548.4A Pending CN115067911A (zh) 2022-05-30 2022-05-30 一种基于gpu实时处理的octa图像优化方法与装置

Country Status (1)

Country Link
CN (1) CN115067911A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116701920A (zh) * 2023-08-08 2023-09-05 北京理工大学 一种提取oct色散失配系数的方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116701920A (zh) * 2023-08-08 2023-09-05 北京理工大学 一种提取oct色散失配系数的方法
CN116701920B (zh) * 2023-08-08 2023-10-31 北京理工大学 一种提取oct色散失配系数的方法

Similar Documents

Publication Publication Date Title
CN109907731B (zh) 基于特征空间的光学相干层析的三维血流造影方法
CN108670239B (zh) 一种基于特征空间的三维血流成像方法与系统
US7995814B2 (en) Dynamic motion contrast and transverse flow estimation using optical coherence tomography
CN106073700B (zh) 图像生成方法和图像生成装置
CN107595250B (zh) 基于运动与图形混合对比度的血流成像方法与系统
CN107788950B (zh) 基于自适应阈值分割的血流成像方法与系统
US10588572B2 (en) Bulk motion subtraction in optical coherence tomography angiography
CN112057049B (zh) 一种基于多维度特征空间的光学相干血流造影方法与系统
CN107862724B (zh) 一种改进的微血管血流成像方法
CN106137134B (zh) 多角度复合的血流成像方法及系统
CN113331809B (zh) 基于mems微型振镜的腔道内三维血流成像方法和装置
JP2014188373A (ja) ３次元光コヒーレンス断層撮影インターフェログラムデータから２次元画像を生成する装置および方法
CN112136182A (zh) 基于Gabor光学相干层析术血流成像的系统和方法
CN115067911A (zh) 一种基于gpu实时处理的octa图像优化方法与装置
CN112396622B (zh) 基于多维特征空间的微血流图像分割量化方法和系统
CN111568373A (zh) 一种重复扫描的octa毛细血管网成像方法
CN113017593B (zh) 血流信号强度分层滤波的血管尾部伪影去除方法与系统
CN111543971B (zh) 时空自适应样本系综去相关运算的血流量化方法与系统
CN114894793A (zh) 基于消除伪影的成像方法、成像系统及服务器
CN111543937A (zh) 一种快速成像相干光层析扫描检眼镜装置
CN113706567A (zh) 一种结合血管形态特征的血流成像量化处理方法与装置
CN113712527A (zh) 一种基于幅度去相关的三维血流成像方法与系统
Chen et al. B-scan-sectioned dynamic micro-optical coherence tomography for bulk-motion suppression
Wijanto et al. Research on Dispersion Compensation of FD-OCT System via Pix2Pix GAN Technique
JP6992031B2 (ja) 画像生成装置、画像生成方法及びプログラム

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