CN112057049B - 一种基于多维度特征空间的光学相干血流造影方法与系统 - Google Patents

一种基于多维度特征空间的光学相干血流造影方法与系统 Download PDF

Info

Publication number
CN112057049B
CN112057049B CN202010962092.7A CN202010962092A CN112057049B CN 112057049 B CN112057049 B CN 112057049B CN 202010962092 A CN202010962092 A CN 202010962092A CN 112057049 B CN112057049 B CN 112057049B
Authority
CN
China
Prior art keywords
signal
signals
dynamic
oct
blood flow
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
CN202010962092.7A
Other languages
English (en)
Other versions
CN112057049A (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.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
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 filed Critical Zhejiang University ZJU
Priority to CN202010962092.7A priority Critical patent/CN112057049B/zh
Publication of CN112057049A publication Critical patent/CN112057049A/zh
Application granted granted Critical
Publication of CN112057049B publication Critical patent/CN112057049B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/7264—Classification of physiological signals or data, e.g. using neural networks, statistical classifiers, expert systems or fuzzy systems
    • 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
    • A61B5/0035—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 adapted for acquisition of images from more than one imaging mode, e.g. combining MRI and optical tomography
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0059—Measuring for diagnostic purposes; Identification of persons using light, e.g. diagnosis by transillumination, diascopy, fluorescence
    • A61B5/0062—Arrangements for scanning
    • A61B5/0066—Optical coherence imaging

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Surgery (AREA)
  • Public Health (AREA)
  • Pathology (AREA)
  • Veterinary Medicine (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Radiology & Medical Imaging (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Fuzzy Systems (AREA)
  • Mathematical Physics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physiology (AREA)
  • Psychiatry (AREA)
  • Signal Processing (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

本发明公开了一种基于多维度特征空间的光学相干血流造影方法与系统。收集三维空间内散射信号样本的OCT散射信号；通过分类器结合形态特征、OCT散射信号的信噪比倒数、去相关系数构建三维特征空间，实现动态血流信号与静态组织信号的分类；包括：采用一阶和零阶自协方差对OCT散射信号处理得到信噪比倒数和去相关系数；构建信噪比倒数‑去相关系数特征空间；确定静态信号分布的上界及动态信号分布的下界，据此将信噪比倒数‑去相关系数特征空间分成三个区域；采用形态掩膜提取过渡区域中的动态信号；结合动态区域的信号和过渡区域提取出的动态信号得到分类结果，生成血流造影结果。本发明能显著抑制噪声带来的伪影问题，提高血流和背景对比度，改善血管连续性。

Description

一种基于多维度特征空间的光学相干血流造影方法与系统
技术领域
本发明大体涉及生物医学成像领域，且更具体地涉及与光学相干层析成像技术(Optical Coherence Tomography，OCT)和血流成像(OCT Angiography，OCTA)相关联的血流造影和基于形态特征、OCT散射信号的信噪比倒数、去相关系数三维特征空间的分类器。
背景技术
血流灌注是衡量生理功能和病理状态的重要指标，目前在临床上常用的血管成像技术需要静脉注射外源性标记物，可能引发的副作用使其不适于对人体血流进行长期的、频繁的跟踪检测。近年来，以光学相干层析技术为基础发展起来的血管造影技术OCTA，以内源性的血流运动代替传统的外源荧光标记物，其非侵入性、无标记的特点，以及对生物组织内的微血管网络进行清晰、可靠的三维成像的能力，使得该技术在被发明以来得到了很快的发展，并在眼底成像和脑皮层血管成像的研究中得到了应用。
为了获取OCTA血流图像，通常需要在生物组织的每个空间位置、以一定的时间间隔进行重复采样(重复的A线扫描或B帧扫描)，每个信号处的运动强度通过分析OCT散射信号的时间动态来进行量化，根据量化得到的运动强度来对血流信号和静态组织信号进行分类。目前已报道的OCTA的血流分类，主要是基于相邻的A线扫描间(或相邻B扫描帧间)的差分、方差或去相关计算。其中基于去相关计算的OCTA血流分类，由于其对于窗口内的多个信号的统计特性的充分利用，因此理论上分类结果的可靠性更高。同时，由于去相关衡量的是相邻B扫描帧间的相似度，因此受整体光源强度变化的影响小。
但是，去相关对于运动对比度的量化效果，对原始的OCT散射信号的噪声水平具有显著的依赖性。随着信号强度的衰减(如在深层组织区域)，随机性噪声将逐渐占据主要成分，也将产生较大的去相关值，带来去相关伪影。基于去相关运算生成的运动对比度无法区分噪声的随机性和血红细胞的运动导致的去相关，因此信噪比较弱的区域容易被误判为血流信号区域，严重影响血流图像的对比度。常见的解决办法是设置一个经验性的强度阈值，生成强度掩膜来去掉所有低信噪比的信号。但是，由于去相关系数和信号强度间存在着复杂的依赖关系，简单的强度掩膜会导致高分类错误率和低运动对比度。
已有的基于信噪比倒数-去相关系数(inverse SNR-decorrelation，ID)特征空间的方法以静态信号在ID空间分布的3σ边界作为分类边界。ID-OCTA算法虽然能去除大部分静态区域，但是将在ID空间和静态信号重叠的动态信号一同去除，影响了血流相对背景噪声的对比度以及血管的连续性。
发明内容
为了解决背景技术中存在的问题，针对现有技术的不足，本发明提出了一种基于形态特征、OCT散射信号的信噪比倒数、去相关系数(shape-inverse SNR-decorrelation，SID)特征空间的光学相干血流造影方法(SID-OCTA)，能够精确区分OCTA图像中的血流区域和混有噪声的静态区域。本发明能显著抑制噪声带来的伪影问题，提高血流和背景对比度，改善血管连续性。
本发明的目的是通过如下技术方案实现的：
一、一种基于多维度特征空间的光学相干血流造影方法：
一种采集器，用于收集三维空间内散射信号样本的OCT散射信号；
一种分类器，结合形态特征、OCT散射信号的信噪比倒数、去相关系数构建三维特征空间，实现动态血流信号与静态组织信号的分类。
所述用于收集三维空间内散射信号样本的OCT散射信号，包括：对散射信号样本进行三维空间的OCT扫描成像，相同空间位置或其附近位置在T个不同时间点重复采样，或进一步结合其他维度，如：波长、角度及偏振等多维度进行并行采样。
所述对散射信号样本进行三维空间的OCT扫描成像，相同空间位置或其附近位置在T个不同时间点重复采样，采用以下方法之一：通过扫描改变参考臂光程的时间域OCT成像方法；利用光谱仪记录光谱干涉信号的光谱域OCT成像方法；利用扫频光源记录光谱干涉信号的扫频OCT成像方法。
所述结合形态特征、OCT散射信号的信噪比倒数、去相关系数构建三维特征空间，实现动态血流信号与静态组织信号的分类，具体包括：
S1、采用一阶和零阶自协方差对OCT散射信号计算分析得到各OCT散射信号的信噪比倒数和去相关系数的两个特征；
S2、构建OCT散射信号的信噪比倒数-去相关系数(inverse SNR-decorrelation，ID)特征空间，确定静态信号分布的上界及动态信号分布的下界并将特征空间划分成三个区域；
采用具有自适应形态阈值的形态掩膜对三个区域中过渡区域的动静态信号进行分类，将从过渡区域分离出的动态信号与动态区域的信号结合得到最终的分类结果并生成最终的血流造影图。
所述步骤S1具体为：采用一阶和零阶自协方差对OCT散射信号计算分析，在三维空间、时间和其他维度(包括波长、角度及偏振等)上做平滑或高斯平均，得到各OCT散射信号的信噪比倒数和去相关系数两个特征。
所述去相关系数的特征包括：对同一空间位置或其附近位置T个不同时间点扫描得到的OCT散射信号的幅度部分或对包含幅度和相位的OCT散射信号计算去相关系数，作为各OCT散射信号的OCTA血流对比度，即作为OCTA血流信息。
所述步骤S2具体为：
构建信噪比倒数-去相关系数(inverse SNR-decorrelation，ID)特征空间，将OCT散射信号的体素映射至信噪比倒数-去相关系数特征空间，特征空间实际为一张特征图；
结合多元时间序列模型及仿真分析，确定静态信号分布的上界及动态信号分布的下界，并根据静态信号分布的上界及动态信号分布的下界将ID特征空间划分成三个区域：动态信号为主的动态区域、静态信号为主的静态区域和动静态信号混合的过渡区域；
将静态区域去除后，利用血管形态评估函数对OCT散射信号的每个体素计算血管测度v(vesselness)，对于血管测度v大于自适应形态阈值vT的体素归为动态，进而生成形态掩膜；
利用形态掩膜对过渡区域的动静态混合信号进行分类；
将过渡区域提取出的动态信号和动态区域的信号结合，得到分类器最终的分类结果，并得到血流造影图。
所述确定静态信号分布的上界及动态信号分布的下界具体为：
根据多元时间序列模型获得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)为一个常矢量。
具体实施中，根据多元时间序列模型，计算得到去相关系数和信噪比倒数的渐近关系。在仿真分析过程中，根据通常认为的瑞利分布对A的强度和均匀分布对其相位进行模拟，并设定白噪声为高斯分布的复随机噪声，以此生成数张模拟图片，将得到的最终模拟信号投影在ID特征空间上。通过计算模拟，拟合静态信号在ID空间的分布方差和信噪比倒数、平均核大小的关系。
具体是通过多元时间序列模型分析，对于运动引起动态程度为α的OCT散射信号，其在ID空间的渐近分布满足：
Figure BDA0002680901580000041
其中，→表示收敛，iSNR表示信噪比倒数，D表示去相关系数，α为与流速/流量成正相关关系的系数；
当计算去相关系数的平均核大小有限时，进一步通过仿真模拟得到对于当α＝0下的静态信号在ID空间的分布方差σs，满足下式：
Figure BDA0002680901580000042
其中，G为仿真分析得到的常数，N为高斯平均核大小，N＝S×T，S表示进行去相关系数计算时所取高维平均核在x-y-z坐标空间的总数，其中：深度z方向为沿入射光轴的方向；x表示OCT横向快扫描方向，y表示慢扫描方向，T表示OCT进行去相关系数计算时所取高维平均核在时间维度t上的数目；
按照以下公式处理获得静态信号分布的上界Ds1和动态信号分布的下界Dd2：
Figure BDA0002680901580000043
Figure BDA0002680901580000044
式中，α0为显著动态的下限；
在信噪比倒数-去相关系数特征空间中，将去相关值大于静态信号分布的上界Ds1的区域作为以动态信号为主的动态区域Rd，将介于动态信号分布的下界Dd2和静态信号分布的上界Ds1之间的区域作为动静态信号混合的过渡区域Rt，将去相关值小于动态信号分布的下界Dd2的区域作为以静态信号为主的静态区域Rs。
所述的形态掩膜具体为：对去除静态区域的三维OCTA造影结果进行黑塞矩阵(Hessian Matrix)分析，得到Hessian矩阵的特征值并构造血管评估函数计算血管测度，将血管测度大于预先设定的自适应形态阈值的体素设置为1，否则设置为0，从而生成形态掩膜。
在信噪比倒数-去相关系数特征空间中，将分布在静态区域Rs的静态信号去除后，采用以下公式的血管形态评估函数对每个体素计算血管测度v(vesselness)：
Figure BDA0002680901580000051
RA＝|λ2|/|λ3|
Figure BDA0002680901580000052
Figure BDA0002680901580000053
|λ1|＜|λ2|＜|λ3|
式中，RA、RB为几何比结构测度，RC为区分背景像素的测度，λ1、λ2、λ3分别为Hessian矩阵的第一、第二、第三特征根，分别代表局部三维空间主曲率的大小(特征根数值大小)和方向(特征根正负)；a，b，c为第一、第二、第三函数灵敏度参数，e表示自然常数；
所述的自适应形态阈值vT具体采用以下公式获得：
Figure BDA0002680901580000054
其中，η为控制阈值参数。
所述基于分类器的分类结果生成血流造影图，包括：
将分类结果中血流信号的部分标记置为1，组织噪声信号的部分标记置为0作为第一种血流造影图，或将第一种血流造影图作为掩膜加在去相关血流图上作为优化后的第一种血流造影图；或将每一OCT散射信号被认为是血流信号的概率值赋为灰度值，灰度值为概率值乘以灰度值区间的计算结果，作为第二种血流造影图，或将第二种血流造影图作为掩膜加在去相关血流图上作为优化后的第二种血流造影图。
二、一种基于多维度特征空间的三维血流成像系统：
一OCT光学相干层析探测装置，用于对三维空间内的散射信号样本的OCT散射信号进行采集；
一个或多个信号处理器，用于获取并分析散射信号的信噪比倒数信息和OCTA血流信息，并综合形态特征，OCT散射信号的信噪比倒数信息和去相关系数特征对动态血流信号和静态组织信号进行分类。
所述的OCT光学相干层析探测装置是采用以下的一种：
包括低相干光源、干涉仪和探测器；
或者包括低相干光源、干涉仪和光谱仪；
或者包括扫频宽光谱光源、干涉仪和探测器。
所述的一OCT光学相干层析探测装置中选择地配置一个可见光指示装置，用于指示OCT探测光束的位置，指导探测目标的放置位置。
所述可见光指示装置主要由可见光指示光源和准直透镜组成。
本发明基于光学相干层析成像(OCT)的无标记、三维、血流运动造影技术。首先基于OCT散射信号的信噪比倒数和去相关系数计算建立ID二维特征空间，结合多元时间序列模型和仿真得到的动态、静态信号分布的边界线，并将ID空间划分成动态区域、静态区域以及动静态信号混合的过渡区域。采用具有自适应形态阈值的形态掩膜对过渡区域的信号进行分类，最终将过渡区域中提取出的动态信号和动态区域的信号结合得到最终的血流造影结果。相比于现有的方法，显著提高了血流和背景的对比度和血管的连续性。
本发明的有益效果和创新点如下：
对比已有技术，本发明结合时间序列模型得到OCT散射信号在ID空间的渐近分布，结合仿真分析定义静态信号分布的上界和动态信号分布的下界。在此基础上，结合具有自适应形态阈值的形态掩膜对过渡区域的信号进行分类，显著提高了血流和背景的对比度和血管的连续性。
本发明的分类器的决定参数仅由系统噪声水平决定，无需其他复杂的对于系统其他参数的标定。
本发明对比已有技术具有以下显著优势：
1.基于去相关计算的OCTA，由于OCT散射信号的去相关系数和信噪比间的依赖关系，低信噪比区域随机噪声引入的去相关伪影无法和血流运动引入的去相关区分。常用解决办法是设置经验性阈值进行强度掩膜，相当于在ID特征空间中用一个强度(信噪比)阈值来分类血流信号和其他信号，而信号的去相关系数和信噪比间更复杂的依赖关系使得实际血流信号和其它信号的分界与强度阈值直线有很大差异，造成较高的误分类率。但是本发明提出的分类器基于对ID空间的定量分析，具有信噪比自适应的优势，此外分类器进一步结合了形态信息，构建了多维度特征分类器。
2.本发明提出具有自适应形态阈值的形态掩膜，利用形态特征对在ID特征空间重叠的过渡区域的动静态信号进行分类，有效地在抑制静态信号的同时，提取出了过渡区域的动态信号；
3.对比现有的方法，本发明提出建立的分类器更加可靠；同时由于去除了大部分的静态和噪声区域，能够提升血管造影图在所有信噪比下的可视度和整体对比度，经过大量样本验证表现显著优于传统方法。
附图说明
图1为本发明方法的示意图；
图2为本发明装置的示意图；
图3为本发明实施例的装置的示意图；
图4为本发明实施例的方法的原理示意图及流程图；
图4(a)为ID空间的划分结果，ID空间被静态信号分布的上界Ds1和动态信号分布的下界Dd2划分成三个区域：以静态信号为主的静态区域Rs，以动态信号为主的动态区域Rd以及动静态信号混合的过渡区域Rt；
图4(b)为λ3＝0.2，0.5，1时血管形态评估函数输出值随|λ1|、|λ2|的变化，其中曲面上的等高线表示阈值敏感参数η＝0.4，0.6，0.8时的阈值。
图4(c)为方法的流程图，①：将OCT散射信号投影到ID空间后，将其划分成Rs、Rd和Rt三个区域，并分别取出Rd+Rt、Rt和Rd三个组分；②：将去除静态信号的三维去相关值(Rd+Rt对应的去相关值)作为血管形态评估函数的输入，采用自适应形态阈值，得到形态掩膜；③：用形态掩膜提取过渡区域的动态信号；④：将动态区域的信号与过渡区域中提取出的动态信号叠加，得到最终的分类结果；
图5为不同方法用于活体小鼠视网膜血流成像实验结果图；
图5(a)为提出的SID-OCTA方法得到的血流最大值投影图及断面图，其中第一道三行分别为浅层、中层和深层的的血流最大值投影图，第四行为断面图，右侧为局部放大图；
图5(b)为ID-OCTA方法得到的血流最大值投影图及断面图；
图5(c)为cmOCT-3σ(阈值为均值加3倍标准差)得到的血流最大值投影图及断面图；
图5(d)为cmOCT-6σ(阈值为均值加6倍标准差)得到的血流最大值投影图及断面图；
图5(e)表示所有动静态信号在ID空间的分布情况；
图5(e1)表示中层低信噪比区域(中层中深色的方框划定的区域)中动静态信号在ID空间的分布情况；
图5(e2)表示中层正常信噪比区域(中层中浅色的方框划定的区域)中动静态信号在ID空间的分布情况；
图6为本发明对小鼠视网膜造影结果及相比传统方法的量化分析对比图；
图6(a)为提出的SID-OCTA方法用于小鼠视网膜时得到的血流造影投影图；
图6(b)为人工标记的中层血管(去除表层大血管的尾影)，作为定量分析的标准；
图6(c)为采用SID-OCTA方法在小鼠视网膜中层的血流造影结果得到的血管与噪声区域的统计直方图；
图6(d)为采用ID-OCTA方法在小鼠视网膜中层的血流造影结果得到的血管与噪声区域的统计直方图；
图6(e)为采用cmOCT-3σ在小鼠视网膜中层的血流造影结果得到的血管与噪声区域的统计直方图；
图6(f)为采用cmOCT-6σ在小鼠视网膜中层的血流造影结果得到的血管与噪声区域的统计直方图；
图中：1-一种采集器；2-一种分类器；21-计算OCTA信号特征；22-将OCT散射信号投影在ID二维特征空间中；23-构建分类器，在形态特征，信噪比倒数，去相关系数三维特征空间对OCT散射信号进行分类；231-根据动静态信号分布的边界，将ID空间分成静态区域、动态区域及过渡区域；232-采用具有自适应形态阈值的形态掩膜提取过渡区域的动态信号；24-将动态区域信号及过渡区域提取出的动态信号叠加构成最终作为血流信号作为最终的分类结果；25-根据分类器的分类结果生成血流造影图。
具体实施方式
下面将结合附图对本发明的具体实施方式做详细说明，附图形成本发明的一部分。需要注意的是，这些说明及示例仅仅为示例性，不能被理解为限制了本发明的范围，本发明的保护范围由随附的权利要求书限定，任何在本发明权利要求基础上的改动都是本发明的保护范围。
本发明的实施例如下：
为了便于理解本发明的实施例，将各操作描述成多个离散的操作，但是描述的顺序不代表实施操作的顺序。
本描述中针对样品测量空间采用基于空间方向的x-y-z三维坐标表示。这种描述仅仅用于促进讨论，而不意欲限制本发明的实施例的应用。其中：深度z方向为沿入射光轴的方向；x-y平面为垂直于光轴的平面，其中x与y正交，且x表示OCT横向快扫描方向，y表示慢扫描方向。
上述s，t，A，n等表示变量，仅仅用于促进讨论，而不意欲限制本发明的实施例的应用，可以是1,2,3等任一数值。为了表述简便，此处略去对于OCT系统波长、角度及偏振等维度进行平均的讨论，仅以时空维度为示例。其实际执行步骤与下述在时空维度上的操作相同。
本发明方法如图1所示，首先为信号采集部分，对组织样本进行OCT三维扫描成像，相同或相邻空间位置在T个不同时间点重复采样1。其次为信号分类部分，结合形态特征、OCT散射信号的信噪比倒数以及去相关系数构建三维特征空间，实现动态血流信号与静态组织背景的分类2。
其具体步骤是：
1)分析血流和周围组织的相对运动，得到各OCT散射信号的信噪比倒数和去相关系数特征21。
去相关系数是对OCT散射信号通过去相关运算处理获得，去相关运算包括对T个不同时间点扫描得到的包含幅度和相位的复数OCT散射信号计算，计算获得去相关系数。其中对复数信号的相关计算理论上可以获得更高的运动对比度。
对于血流和周围组织中的某一局部区域，针对每个体素用其相邻的T个OCT扫描的B扫描帧(x-z平面)进行平均计算(即使用高维平均核进行卷积)获得每个体素的一阶与零阶自协方差和去相关系数：
Figure BDA0002680901580000091
Figure BDA0002680901580000092
Figure BDA0002680901580000093
其中，C表示一阶自协方差，I表示零阶自协方差，即强度，D表示去相关系数，作为OCTA血流信息；*表示复数的共轭，X(s，t)是某一个空间位置(z，x)在时刻t的复数信号；S表示进行去相关系数计算时所取高维平均核在x-y-z空间的总数；s表示进行去相关系数计算时所取高维平均核在x-y-z空间的序数；T表示进行去相关系数计算时所取高维平均核在时间维度t上的总数，即为OCT扫描的相邻B扫描帧的帧数量；t表示进行去相关系数计算时所取高维平均核在时间维度上的序数。
计算过程中，先对于每个体素采用上述公式计算其一阶和零阶自协方差，在时间和空间等各个维度上进行平均从而得到整个散射信号样品的扫描体积内所有体素的去相关值，能够提升运算速度。
2)在OCT系统中，噪声来源主要为散粒噪声，视为在整个扫描体积内近似不变，可以通过扫描空气的数据或事先对系统进行标定而确定。
然后采用以下公式计算每个体素的信噪比倒数iSNR，定义如下：
Figure BDA0002680901580000101
其中，s2是OCT系统的噪声水平。
通过计算OCT散射信号中数据空白区域的强度或者针对OCT系统事先标定得到噪声水平。
3)结合上述OCT探测得到的信噪比倒数和基于去相关计算得到的OCTA血流信息建立ID二维特征空间，并将OCT散射信号投影在特征空间中22。
根据通过多元时间序列模型分析的结果，对于运动引起动态程度为α的OCT散射信号，其在ID空间的渐近分布满足：
Figure BDA0002680901580000102
其中，→表示收敛，iSNR表示信噪比倒数，平均核大小充分大时能实现收敛。平均体素只能代表x-y-z空间的平均核数目，平均核包括体素数目和重复的B-scan数目。当平均核大小有限时，进一步通过仿真模拟得到对于静态信号(α＝0)在ID空间的分布方差σs，满足下式：
Figure BDA0002680901580000103
其中G≈1.5为仿真分析得到的常数，N＝S×T为高斯平均核大小。
4)之后根据多元时间序列模型及仿真分析，结合经验值在ID空间定义静态信号分布的上界Ds1和动态信号分布的下界Dd2，定义如下：
Figure BDA0002680901580000104
Figure BDA0002680901580000105
式中，α0为显著动态的下限，这里以0.4为例。
根据定义的(7)式和(8)式可以将ID空间分成三个部分231：以静态信号为主的静态区域Rs，以动态信号为主的动态区域Rd，动静态信号混合的过渡区域Rt。
5)将分布在静态区域Rs的静态信号去除后，采用血管形态评估函数对每个体素进行评估。函数的输入为将静态区域Rs的去相关值置零后的三维去相关系数值，输出为血管测度v(vesselness)，定义为：
Figure BDA0002680901580000111
式中，RA＝|λ2|/|λ3|、
Figure BDA0002680901580000112
为几何比结构测度，
Figure BDA0002680901580000113
Figure BDA0002680901580000114
为区分背景像素的测度，λ1、λ2、λ3分别为黑塞矩阵(Hessian Matrix)的第一、第二、第三特征根；a，b，c为第一、第二、第三函数灵敏度参数，e表示自然常数；
进一步提出基于形态的自适应形态阈值vT，定义为：
Figure BDA0002680901580000115
其中，η为介于0和1之间的常数，事先选定，用于控制阈值的整体水平，后续以0.6为例。
基于上述提出的具有自适应形态阈值的形态掩膜(将v≥vT的体素识别为动态，其余为静态)，对过渡区域的混合信号进行分类，提取出其中的动态信号232。
6)将动态区域Rd的信号以及采用形态掩膜识别出的过渡区域Rt中的动态信号叠加获得最终的分类结果24，并进一步得到血流造影结果25。
图2示出的是本发明中基于形态特征、OCT散射信号信噪比倒数以及去相关系数特征空间的OCT血流造影技术的采集装置结构示意图。该装置的低相干干涉测量部分的主体结构为一干涉仪，由11～23构成。光源11连接到分束器12一侧的输入端，光源11发出的光被分束器12分成两部分光束：其中的一束光经一偏振控制器13进入到干涉仪的参考臂，通过参考臂准直镜14照射于参考的平面反射镜15上；另一束光经另一偏振控制器13进入到样品臂，具体经过准直透镜16和扫描装置光路聚焦到待测样品21上。其中扫描装置光路中，光束经过二维扫描振镜组17、18，“4f”透镜组54、55和二向色镜19的反射后，经过物镜20聚焦在待测样品21上。其中透镜组54、55是由两个透镜54、55同光轴布置组成，透镜组54、55的设计是为了确保扫描时二维扫描振镜镜面的光束中心和二相色镜反射面的光束中心固定不变，使得OCT样品臂中的光束在扫描时不影响物镜的成像特性。
然后，参考臂和样品臂各自反射回的光经原路返回到分束器12输出，发生干涉后由干涉信号探测装置22接收，干涉信号探测装置22再连接到信号处理器模块与计算单元23。对于光纤型光路，采用两个偏振控制器13调整光束的偏振态，最大化信号干涉效果。
具体实施还设置有可见光指示装置，可见光指示装置包括低功率可见光光源25、准直透镜24和滤光片52，用于指示的可见光依次经过准直透镜24、滤光片52、二向色镜19和聚焦物镜20后到待测样品21。
依据低相干干涉探测信号的不同方式，图2所示的一种基于多维度特征空间的三维血流成像系统装置具体包括：
1)时间域测量装置。光源11采用宽带低相干光，平面反射镜15可沿光轴方向移动，干涉信号探测装置22为一点探测器。通过移动平面反射镜15改变参考臂光程，两臂的干涉信号由点探测器22探测到，对某一空间深度的z方向的散射信号的低相干干涉探测，从而得到深度空间维度的采样体。
2)光谱域测量装置。光源11采用宽带低相干光，平面反射镜15固定不动，干涉信号探测装置22采用光谱仪。干涉信号经过光谱仪22中的线阵相机同时记录干涉光谱。采用傅里叶分析方法分析干涉光谱信号，并行获取深度z方向的散射信息，从而得到深度空间维度的采样体。
3)扫频测量装置。光源11采用扫频光源，平面反射镜15固定不动，干涉信号探测装置22采用点探测器。点探测器22分时记录扫频光源的低相干干涉光谱。采样傅里叶分析干涉光谱信号，并行获取深度z方向的散射信息，从而得到深度空间维度的采样体。
图3示出的是利用本发明的一个示例性实施例。基于多维度特征空间的三维OCT血管造影系统，包括宽带低相干光源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。
图5为不同方法对小鼠视网膜进行血流造影的结果及信号在ID空间的分布情况。SID-OCTA(图5(a))和ID-OCTA(图5(b))都能得到较低的背景噪声，而SID-OCTA的血流信号相对背景有更高的对比度，并且血管连续性更好，说明SID-OCTA能够在抑制过渡区域噪声信号的同时有效地提取出动态信号。而在cmOCT中，采用一个强度阈值将低信噪比的信号去除，当采用低阈值时(图5(c))，背景噪声较大，降低了血流信号相对背景的对比度；采用高阈值时(图5(d))，很多血流信号被除去，影响了血管的连续性。
SID-OCTA对血管连续性的提升在低信噪比区域(图5中层视网膜中深色的方框区域)更加明显。在低信噪比情况下，大量的动态信号会偏移到过渡区域，被ID-OCTA及cmOCT-6σ划分为静态信号，故而严重影响了血管的连续性，而采用SID-OCTA时，过渡区域的动态信号可以被有效地提取出来，提升了血管的连续性。
进一步地，为了量化所提出的SID-OCTA方法带来的对比度的提升，本发明基于对直方图统计结果的分析，并采用血管和噪声的重叠程度(用归一化后的频率衡量)及对比度-噪声比(contrast-noise-ratio,CNR)来表征不同算法的性能，CNR由下式计算得到：
Figure BDA0002680901580000131
其中，
Figure BDA0002680901580000141
分别是血管区域和背景噪声区域的去相关系数的均值，σn是背景噪声区域的标准差。
图6为对小鼠视网膜中层的血流造影结果的直方图统计分析及定量评估。图6(a)为采用所提出的SID-OCTA对小鼠视网膜造影的投影图结果，图6(b)为人工划定的中层血管(浅层的大血管的尾影被去除)。对比图6(c)(SID-OCTA)和图6(d)(ID-OCTA)，两者噪声的直方图分布结果类似，而SID-OCTA得到的血管的直方图分布在高去相关值位置的频率更高，说明SID-OCTA能在抑制背景噪声的同时有效地提取出血流信号。
定量分析结果表明，所提出的SID-OCTA方法得到的血流造影结果中，动静态的重叠更少(SID-OCTA：22.3％，ID-OCTA：47.9％；cmOCT-3σ：49.3％，cmOCT-6σ：55.7％)。此外，SID-OCTA显著提高了CNR值(4.12)，和ID-OCTA(3.04)相比有36％的提升，和cmOCT-3σ(2.11)相比有95％的提升，和cmOCT-6σ(2.40)相比有71％的提升。
上述实验对比结果充分说明：利用本发明所涉及的基于多维度特征空间的光学相干血流造影方法，能够提升血流信号分类的准确率，实现血流对比度的有效增强，和血流图像质量的提升，具有其突出显著的技术效果。

Claims (9)

1.一种基于多维度特征空间的光学相干血流造影方法，其特征在于包括：
一种采集器(1)，用于收集三维空间内散射信号样本的OCT散射信号；
一种分类器(2)，结合形态特征、OCT散射信号的信噪比倒数、去相关系数构建三维特征空间，实现动态血流信号与静态组织信号的分类；
所述结合形态特征、OCT散射信号的信噪比倒数、去相关系数构建三维特征空间，实现动态血流信号与静态组织信号的分类，具体包括：
S1、采用一阶和零阶自协方差对OCT散射信号计算分析得到各OCT散射信号的信噪比倒数和去相关系数的两个特征(21)；
S2、构建OCT散射信号的信噪比倒数-去相关系数(inverse SNR-decorrelation，ID)特征空间(22)，确定静态信号分布的上界及动态信号分布的下界并将特征空间划分成三个区域(231)；
采用具有自适应形态阈值的形态掩膜对三个区域中过渡区域的动静态信号进行分类(232)，将从过渡区域分离出的动态信号与动态区域的信号结合得到最终的分类结果(24)并生成最终的血流造影图(25)。
2.根据权利要求1所述的一种基于多维度特征空间的光学相干血流造影方法，其特征在于：所述用于收集三维空间内散射信号样本的OCT散射信号，包括：对散射信号样本进行三维空间的OCT扫描成像，相同空间位置或其附近位置在T个不同时间点重复采样。
3.根据权利要求2所述的一种基于多维度特征空间的光学相干血流造影方法，其特征在于：所述对散射信号样本进行三维空间的OCT扫描成像，相同空间位置或其附近位置在T个不同时间点重复采样，采用以下方法之一：通过扫描改变参考臂光程的时间域OCT成像方法；利用光谱仪记录光谱干涉信号的光谱域OCT成像方法；利用扫频光源记录光谱干涉信号的扫频OCT成像方法。
4.根据权利要求1所述的一种基于多维度特征空间的光学相干血流造影方法，其特征在于：所述步骤S1具体为：采用一阶和零阶自协方差对OCT散射信号计算分析，在三维空间、时间、波长、角度及偏振上做平滑或高斯平均，得到各OCT散射信号的信噪比倒数和去相关系数两个特征(21)。
5.根据权利要求1所述的一种基于多维度特征空间的光学相干血流造影方法，其特征在于：所述去相关系数的特征包括：对同一空间位置或其附近位置T个不同时间点扫描得到的OCT散射信号的幅度部分或对包含幅度和相位的OCT散射信号计算去相关系数，作为各OCT散射信号的OCTA血流对比度，即作为OCTA血流信息。
6.根据权利要求1所述的一种基于多维度特征空间的光学相干血流造影方法，其特征在于：所述步骤S2具体为：
构建信噪比倒数-去相关系数特征空间(22)，将OCT散射信号的体素映射至信噪比倒数-去相关系数特征空间(22)；
结合多元时间序列模型及仿真分析，确定静态信号分布的上界及动态信号分布的下界，并根据静态信号分布的上界及动态信号分布的下界将ID特征空间划分成三个区域：动态信号为主的动态区域、静态信号为主的静态区域和动静态信号混合的过渡区域(231)；
将静态区域去除后，利用血管形态评估函数对OCT散射信号的每个体素计算血管测度v，将血管测度v大于自适应形态阈值vT的体素归为动态，进而生成形态掩膜；
利用形态掩膜对过渡区域的动静态混合信号进行分类(232)；
将过渡区域提取出的动态信号和动态区域的信号结合，得到分类器最终的分类结果(24)，并得到血流造影图(25)。
7.根据权利要求1所述的一种基于多维度特征空间的光学相干血流造影方法，其特征在于：所述确定静态信号分布的上界及动态信号分布的下界具体为：
根据多元时间序列模型获得OCT散射信号在信噪比倒数-去相关系数特征空间(22)的渐近分布，借助仿真分析得到渐近分布的方差，分别确定静态信号分布的上界及动态信号分布的下界。
8.根据权利要求1所述的一种基于多维度特征空间的光学相干血流造影方法，其特征在于：所述的形态掩膜具体为：对去除静态区域的三维OCTA造影结果进行黑塞矩阵(Hessian Matrix)分析，得到Hessian矩阵的特征值并构造血管评估函数计算血管测度，将血管测度大于预先设定的自适应形态阈值的体素设置为1，否则设置为0，从而生成形态掩膜。
9.一种实施权利要求1～8任一所述方法的基于多维度特征空间的光学相干血流造影系统，包括：
一OCT光学相干层析探测装置，用于对三维空间内的散射信号样本的OCT散射信号进行采集；
一个或多个信号处理器，用于获取并分析散射信号的信噪比倒数信息和OCTA血流信息，并综合形态特征、OCT散射信号的信噪比倒数信息和OCTA血流信息对动态血流信号和静态组织信号进行分类。
CN202010962092.7A 2020-09-14 2020-09-14 一种基于多维度特征空间的光学相干血流造影方法与系统 Active CN112057049B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010962092.7A CN112057049B (zh) 2020-09-14 2020-09-14 一种基于多维度特征空间的光学相干血流造影方法与系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010962092.7A CN112057049B (zh) 2020-09-14 2020-09-14 一种基于多维度特征空间的光学相干血流造影方法与系统

Publications (2)

Publication Number Publication Date
CN112057049A CN112057049A (zh) 2020-12-11
CN112057049B true CN112057049B (zh) 2021-08-10

Family

ID=73695600

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010962092.7A Active CN112057049B (zh) 2020-09-14 2020-09-14 一种基于多维度特征空间的光学相干血流造影方法与系统

Country Status (1)

Country Link
CN (1) CN112057049B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2022144476A (ja) * 2021-03-19 2022-10-03 株式会社トプコン グレード評価装置、眼科撮影装置、プログラム、記録媒体、およびグレード評価方法
CN113706567A (zh) * 2021-07-19 2021-11-26 浙江大学 一种结合血管形态特征的血流成像量化处理方法与装置
CN113712527A (zh) * 2021-07-19 2021-11-30 浙江大学 一种基于幅度去相关的三维血流成像方法与系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2010106913A1 (ja) * 2009-03-19 2010-09-23 富士フイルム株式会社 光立体構造計測装置及びその構造情報処理方法
CN104523233A (zh) * 2014-12-29 2015-04-22 浙江大学 基于复数互相关的微血管光学造影及抖动补偿方法与系统
CN105286779A (zh) * 2015-10-30 2016-02-03 温州医科大学 一种在体视网膜血流动力学的成像与绝对流速测量方法
CN109907731A (zh) * 2019-01-31 2019-06-21 浙江大学 基于特征空间的光学相干层析的三维血流造影方法及系统
CN110415216A (zh) * 2019-07-01 2019-11-05 南京理工大学 基于sd-oct和octa视网膜图像的cnv自动检测方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10327635B2 (en) * 2016-05-03 2019-06-25 Oregon Health & Science University Systems and methods to compensate for reflectance variation in OCT angiography

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2010106913A1 (ja) * 2009-03-19 2010-09-23 富士フイルム株式会社 光立体構造計測装置及びその構造情報処理方法
CN104523233A (zh) * 2014-12-29 2015-04-22 浙江大学 基于复数互相关的微血管光学造影及抖动补偿方法与系统
CN105286779A (zh) * 2015-10-30 2016-02-03 温州医科大学 一种在体视网膜血流动力学的成像与绝对流速测量方法
CN109907731A (zh) * 2019-01-31 2019-06-21 浙江大学 基于特征空间的光学相干层析的三维血流造影方法及系统
CN110415216A (zh) * 2019-07-01 2019-11-05 南京理工大学 基于sd-oct和octa视网膜图像的cnv自动检测方法

Also Published As

Publication number Publication date
CN112057049A (zh) 2020-12-11

Similar Documents

Publication Publication Date Title
CN109907731B (zh) 基于特征空间的光学相干层析的三维血流造影方法
CN107788950B (zh) 基于自适应阈值分割的血流成像方法与系统
CN112057049B (zh) 一种基于多维度特征空间的光学相干血流造影方法与系统
CN108670239B (zh) 一种基于特征空间的三维血流成像方法与系统
CN107595250B (zh) 基于运动与图形混合对比度的血流成像方法与系统
US7995814B2 (en) Dynamic motion contrast and transverse flow estimation using optical coherence tomography
WO2017133083A1 (zh) 基于全空间调制谱分割角度复合的微血管造影方法与系统
US8721077B2 (en) Systems, methods and computer-readable medium for determining depth-resolved physical and/or optical properties of scattering media by analyzing measured data over a range of depths
US20070025642A1 (en) Methods, systems and computer program products for analyzing three dimensional data sets obtained from a sample
CN105026879B (zh) 光学成像设备和用于对样本成像的方法
CN110693457B (zh) 一种基于光学相干技术的组织活性检测的方法与系统
CN106137134B (zh) 多角度复合的血流成像方法及系统
CN106073700A (zh) 图像生成方法和图像生成装置
CN107862724B (zh) 一种改进的微血管血流成像方法
CN109963494A (zh) 具有改进的图像质量的光相干断层成像系统
JP2014188373A (ja) ３次元光コヒーレンス断層撮影インターフェログラムデータから２次元画像を生成する装置および方法
US20160367146A1 (en) Phase Measurement, Analysis, and Correction Methods for Coherent Imaging Systems
CN112136182A (zh) 基于Gabor光学相干层析术血流成像的系统和方法
CN113017593B (zh) 血流信号强度分层滤波的血管尾部伪影去除方法与系统
CN111568373A (zh) 一种重复扫描的octa毛细血管网成像方法
CN111543971B (zh) 时空自适应样本系综去相关运算的血流量化方法与系统
CN115067911A (zh) 一种基于gpu实时处理的octa图像优化方法与装置
CN112396622B (zh) 基于多维特征空间的微血流图像分割量化方法和系统
CN113706567A (zh) 一种结合血管形态特征的血流成像量化处理方法与装置
CN113712527A (zh) 一种基于幅度去相关的三维血流成像方法与系统

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