CN113706567A - 一种结合血管形态特征的血流成像量化处理方法与装置 - Google Patents

Abstract

本发明公开了一种结合血管形态特征的血流成像量化处理方法与装置。收集散射信号样品的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)特征空间的血管网络量化检测方法。本发明能显著抑制随机噪声的干扰，提高动静态信号的分类精度，改善二值化血管网络的连续性。

本发明的目的是通过如下技术方案实现的：

一、一种结合血管形态特征的血流成像量化处理方法：

一种散射信号采集方式，基于光学相干层析成像技术(OCT)采集三维空间内散射信号样本的OCT散射信号；

一种血流图像分割方法，结合形态特征、OCT散射信号的信噪比倒数、去相关系数构建三维特征空间，实现动态血流信号与静态组织信号的分类，得到二值化血管网络图像；

一种微血管形态量化处理方法，根据二值化血管网络图像进行血流骨架和轮廓提取获得血流骨架图，进而根据血流骨架图计算出反映血流形态的多种量化参数，多种量化参数包括血流平均管径、血流面积密度、血流单位面积长度和血流单位面积周长。

所述的散射信号样本为生物组织样本，生物组织样本例如可以为人和其他动物的皮肤、脑组织、眼睛。

所述一种散射信号采集方式，包括：对散射信号样本进行三维空间的OCT扫描成像，相同空间位置或其附近位置在T个不同时间点重复采样，采用以下方法之一：通过扫描改变参考臂光程的时间域OCT成像方法；利用光谱仪记录光谱干涉信号的光谱域OCT成像方法；利用扫频光源记录光谱干涉信号的扫频OCT成像方法。

所述的一种血流图像分割方法，具体包括：

S1、采用一阶和零阶自协方差对OCT散射信号计算分析，得到各OCT散射信号的信噪比倒数和去相关系数两个特征，所得到的信噪比倒数和去相关系数进一步在三维空间、时间、角度及偏振态等多个维度上做滑动平均或高斯平均，利用平均处理后的信噪比倒数和去相关系数的两个特征构建OCT散射信号的信噪比倒数-去相关系数(inverse SNR-decorrelation，ID)特征空间；

S2、基于形态特征、信噪比倒数、去相关系数多维度特征空间对信号进行分类，包括：遍历获得信噪比倒数-去相关系数特征空间中经过原点的两条线性分类边界，结合形态特征对三维空间的血管网络图像进行二值化处理获得二值化体数据，计算二值化体数据的结构相似度值，遍历所有线性分类边界的角度组合后，选取最小的结构相似度值对应的二值化结果作为最终的二值化血管网络。

所述S2具体为：

遍历信噪比倒数-去相关系数特征空间中经过原点的每两条分割阈值线，通过两条分割阈值线将信噪比倒数-去相关系数特征空间划分为动态区域、中间区域和静态区域；分割阈值线都从原点出发，通过两条分割阈值线将信噪比倒数-去相关系数特征空间分割为三部分区域，靠近去相关系数所在坐标轴的区域作为动态区域，靠近信噪比倒数所在坐标轴的区域作为静态区域，动态区域和静态区域之间的为中间区域；

由动态区域和中间区域共同作形态滤波，形态滤波的结果利用预设阈值二值化形成形态掩膜，借助形态特征构建形态掩膜提取出中间区域的动态信号；

根据对动静态信号的分类结果对三维空间的血管网络图像进行二值化处理获得二值化体数据，根据对动静态信号的分类结果，在三维空间中计算二值化体数据的结构相似度值BVSIM；

三维空间的血管网络图像通常是由OCT散射信号计算出的去相关系数构建图像获得。

遍历每两条分割阈值线形成各种两条分割阈值线的角度组合后，选取最小的结构相似度值BVSIM对应的两条分割阈值线作为两条线性分类边界，根据两条线性分类边界并结合形态掩膜区分中间区域的动静态信号，生成二值化血管网络。

所述的分割阈值线为在信噪比倒数-去相关系数二维特征空间中的一条过原点的直线，原点是指信噪比倒数-去相关系数特征空间中去相关系数和信噪比倒数均为零的位置，去相关系数和信噪比倒数都是非负数，分割阈值线和信噪比倒数-去相关系数二维特征空间的坐标系横坐标轴的夹角为分割阈值线角度。

所述S2具体为：

S21、随机在信噪比倒数-去相关系数特征空间中建立经过原点的每两条分割阈值线，结合形态掩膜对信号实现初步的分类，分为初步静态信号和初步动态信号；

S22、先在信噪比倒数-去相关系数特征空间中对初步动态信号生成一系列经过原点的多条分割线，一系列分割线和去相关系数所在坐标轴之间的夹角逐渐增大，每两条分割线间包含1/n的信噪比倒数-去相关系数特征空间中的总体素数目，利用分割线对动态区域进行二值化分割后获得一系列二值化体数据，将各个二值化体数据按照分割线的角度递增顺序组成一个序列作为初步动态信号的二值化体数据序列，计算初步动态区域内的体数据间的结构相似度，具体为：

首先，按照以下公式处理获得每个二值化体数据在位置(z，x，y)处的结构向量

其中，B(α，z+h，x+i，y+j)表示二值化体数据中坐标(z+h，x+i，y+j)处的值，α为二值化体数据对应的分割阈值线相对于去相关系数所在坐标轴的角度，k表示结构向量的窗口大小，h、i和j表示窗口内像素的三个坐标的索引，(h，i，j)表示一个三维矢量，三维矢量的大小和方向由h、i和j决定；

然后，按照以下公式计算各个二值化体数据的局部的结构差异值之和，作为整个区域的结构相似度值：

其中，m、l分别表示区域内二值化体数据序列中的二值化体数据的序号，V表示区域内的每两个二值化体数据之间的图像结构相似度总和，即区域的结构相似度值，Δv(m，l)表示第m个二值化体数据和第l个二值化体数据之间的结构差异度，|*|表示欧几里德距离，Z、X和Y分别是OCT深度方向、快扫描方向和慢扫描方向的总像素数；

S23、按照和S22同样处理方式计算初步静态区域内的体数据间的结构相似度；

S24、综合动静态区域内的体数据间的结构相似度，得到最终的二值化体数据的结构相似度值BVSIM，具体公式为：

其中，V_d代表动态区域的结构相似度，V_s代表静态区域的结构相似度，n_d和n_s分别代表动态区域、静态区域内的二值化体数据的数量，

表示从n_d个元素中选取2个元素的所有组合的个数，

表示从n_s个元素中选取2个元素的所有组合的个数。

当结构相似度值BVSIM最小时，以当前的线性分类边界的角度作为阈值，对应的动态信号(包括动态区域的信号以及采用形态掩膜在中间区域提取出的动态信号)作为血流信号，生成二值化微血管网络。

所述一种微血管形态量化处理方法中，根据二值化血管网络图像进行血流骨架和轮廓提取获得血流骨架图，具体为：在二值化血管网络图像中沿水平面建立横向和竖向两个方向，在横向和竖向两个方向上分别对每两个相邻像素进行差分运算，进而得到血流边缘图；再迭代删除二值化血管网络图像中的血流区域外部像素，直到得到宽度为单像素的三维血流骨架，得到血流骨架图。

二、基于多维特征空间的微血流图像分割量化系统：

一OCT光学相干层析探测装置，用于对三维空间内的散射信号样本的OCT散射信号进行采集；

一个图像处理器，用于用于获取并分析OCT散射信号的信噪比倒数和和去相关系数，并结合形态特征，进行动态血流信号与静态组织信号的分类，得到二值化血管网络图像；

一个数据处理器，根据二值化血管网络图像进行血流骨架和轮廓提取获得血流骨架图，进而根据血流骨架图计算出反应血流形态的多种量化参数，多种量化参数包括血流平均管径、血流面积密度、血流单位面积长度和血流单位面积周长。

所述的一OCT光学相干层析探测装置是采用以下的一种：

包括低相干光源、干涉仪和探测器；

或者包括低相干光源、干涉仪和光谱仪；

或者包括扫频宽光谱光源、干涉仪和探测器。

所述的一OCT光学相干层析探测装置中选择地配置一个可见光指示装置，用于指示OCT探测光束的位置，指导探测目标的放置位置。

方法中以形态、信噪比倒数和去相关系数三个特征描述OCT散射信号，构建了基于多维度特征空间的分类器。然后处理中，以二值化体数据的结构相似度值(BVSIM)衡量两个体数据的结构相似程度。根据采用不同分割边界线组合情况下的BVSIM值，二值化体数据的结构相似度值自动确定最优分割边界线组合阈值，进行血流图像分割。

本发明方法的二值化阈值是随信噪比自适应的变化的，所以可以有效保留低信噪比处血流信号，得到更好的血流图像。

本发明的有益效果和创新点如下：

对比已有技术，本发明利用OCT散射信号的三个特征(血管形态、信噪比倒数和去相关系数)，结合二值化微血管网络的结构相似性，建立信噪比自适应的分类器，可以更好地保留信噪比较低的血流信号，得到更准确的二值化微血管网络。同时提出了量化微血流形态特征的方法，可用于检测与血管形态变化关联的疾病。

本发明对比已有技术具有以下显著优势：

1.基于去相关计算的OCTA，由于OCT散射信号的去相关系数和信噪比间的依赖关系，低信噪比区域随机噪声引入的去相关伪影无法和血流运动引入的去相关区分。常用解决办法是设置经验性阈值进行强度掩膜，相当于在ID特征空间中用一个强度(信噪比)阈值来去除所有低信噪比信号，而信号的去相关系数和信噪比间更复杂的依赖关系使得实际血流信号和其它信号的边界与强度阈值直线有很大差异，造成较高的误分类率。但是本发明提出的分类器基于对ID空间的定量分析，具有信噪比自适应的优势，此外分类器进一步结合了血管形态特征，构建了多维度特征分类器。

2.本发明提出了具有自适应形态阈值的形态掩膜，利用形态特征对ID特征空间中重叠在中间区域的动静态信号进行分类，在有效抑制中间区域的静态信号的同时，提取出了中间区域的动态信号；

3.对比现有的方法，本发明建立的分类器更加可靠；同时由于去除了大部分的静态和噪声区域，能够提升血管造影图在所有信噪比下的可视度和整体对比度，经过大量样本验证表现显著优于传统方法。

4.血流分割阈值线仅由图像处理器自动搜索，无需其他复杂的对于系统其他参数的标定或对相关算法的复杂矫正；

5.由于微血管形态变化可反映多种疾病的发展，所以对微血管的形态特征进行量化分析，有助于提前发现疾病，辅助临床诊断。

附图说明

图1为本发明方法的示意图；

图2为本发明装置的示意图；

图3为本发明实施例的装置的示意图；

图4为本发明实施例中分类方法的原理示意图及流程图；

图4(a)为ID空间的划分结果，ID空间被随机选取的经过原点的两条分割阈值线划分成三个区域：以静态信号为主的静态区域R_s，以动态信号为主的动态区域R_d以及动静态信号混合的中间区域R_i；

图4(b)为λ₃＝0.2，0.5，1时血管形态评估函数输出值随|λ₁|、|λ₂|的变化，其中曲面上的等高线表示阈值敏感参数η＝0.4，0.6，0.8时的阈值。

图4(c)为方法的流程图，①：将OCT散射信号投影到ID空间后，将其划分成R_s、R_d和R_i三个区域，并分别取出R_d+R_i、R_i和R_d三个组分；②：将去除静态信号的三维去相关值(R_d+R_i对应的去相关值)作为血管形态评估函数的输入，采用自适应形态阈值，得到形态掩膜；③：用形态掩膜提取中间区域的动态信号；④：将动态区域的信号与中间区域中提取出的动态信号叠加，得到最终的分类结果。

具体实施方式

下面将结合附图对本发明的具体实施方式做详细说明，附图形成本发明的一部分。需要注意的是，这些说明及示例仅仅为示例性，不能被理解为限制了本发明的范围，本发明的保护范围由随附的权利要求书限定，任何在本发明权利要求基础上的改动都是本发明的保护范围。

本发明的实施例如下：

为了便于理解本发明的实施例，将各操作描述成多个离散的操作，但是描述的顺序不代表实施操作的顺序。

本描述中针对样品测量空间采用基于空间方向的x-y-z三维坐标表示。这种描述仅仅用于促进讨论，而不意欲限制本发明的实施例的应用。其中：深度z方向为沿入射光轴的方向；x-y平面为垂直于光轴的平面，其中x与y正交，且x表示OCT横向快扫描方向，y表示慢扫描方向。

上述

V，m，l等表示变量，仅仅用于促进讨论，而不意欲限制本发明的实施例的应用，可以是1，2，3等任一数值。为了表述简便，此处略去对于OCT系统波长、角度及偏振等维度进行平均的讨论，仅以时空维度为示例。其实际执行步骤与下述在时空维度上的操作相同。

本发明方法如图1所示，首先为信号采集部分，对组织样本进行OCT三维扫描成像，相同或相邻空间位置在T个不同时间点重复采样。其次为信号分类部分，结合形态特征、OCT散射信号的信噪比倒数以及去相关系数构建三维特征空间，生成二值化血管网络。最后根据二值化血管网络计算血管量化参数。

其具体步骤是：

1)分析血流和周围组织的相对运动，得到各OCT散射信号的信噪比倒数和去相关系数特征21。

去相关系数是对OCT散射信号通过去相关运算处理获得，去相关运算包括对T个不同时间点扫描得到的包含幅度和相位的复数OCT散射信号计算，计算获得去相关系数。其中对复数信号的相关计算理论上可以获得更高的运动对比度。

对于血流和周围组织中的某一局部区域，针对每个体素用其相邻的T个OCT扫描的B扫描帧(x-z平面)进行平均计算(即使用高维平均核进行卷积)获得每个体素的一阶与零阶自协方差和去相关系数：

其中，C表示一阶自协方差，I表示零阶自协方差，即强度，D表示去相关系数，作为OCTA血流信息；*表示复数的共轭，X(s，t)是某一个空间位置(z，x，y)在时刻t的复数信号；S表示进行去相关系数计算时所取高维平均核在x-y-z空间的总数；s表示进行去相关系数计算时所取高维平均核在x-y-z空间的序数；T表示进行去相关系数计算时所取高维平均核在时间维度t上的总数，即为OCT扫描中同一空间位置B扫描帧的帧数量；t表示进行去相关系数计算时所取高维平均核在时间维度上的序数。C表示一阶自协方差，I表示零阶自协方差。

计算过程中，先对于每个体素采用上述公式计算其一阶和零阶自协方差，在时间和空间等各个维度上进行平均从而得到整个散射信号样品的扫描体积内所有体素的去相关值，能够提升运算速度。

2)在OCT系统中，噪声来源主要为散粒噪声，视为在整个扫描体积内近似不变，可以通过计算断层图中顶部空气区域及底部噪声区域的OCT信号均值获得。

然后采用以下公式计算每个体素的信噪比倒数iSNR，定义如下：

其中，s²是OCT系统的噪声水平。I表示零阶自协方差。

3)结合上述OCT探测得到的信噪比倒数和基于去相关计算得到的OCTA血流信息建立ID二维特征空间，并将OCT散射信号投影在特征空间中22。之后计算任意分割阈值线角度组合下的BVSIM值，具体为：

任取角度组合α₁和α₂(0°＜α₁＜α₂＜90°，步长越短阈值计算精度越高，可根据需求调整，本发明中为方便描述设定步长为1°)，根据如下形式定义分割边界线：

D_1/2＝cot(α_1/2)·iSNR (5)

其中，D_1/2表示两条分割阈值线的去相关值D和信噪比倒数iSNR的关系，α_1/2表示两条分割阈值线的角度；

据此将ID特征空间划分成动态区域R_d、中间区域R_i和静态区域R_s，继而结合自适应管状掩膜技术对中间的信号进行区分，具体为：

将分布在静态区域R_s的静态信号去除后，采用血管形态评估函数对每个体素进行评估。函数的输入为将静态区域R_s的去相关值置零后的三维去相关系数值，输出为血管测度v(vesselness)，定义为：

式中，R_A、R_B为第一、第二几何比结构测度，R_C为区分背景像素的测度，λ₁、λ₂、λ₃分别为黑塞矩阵(Hessian Matrix)的第一、第二、第三特征根；a、b、c为第一、第二、第三函数灵敏度参数，e表示自然常数；

进一步提出基于形态的自适应形态阈值v_T，定义为：

其中，η为阈值整体水平参数，具体为介于0和1之间的常数，事先选定，用于控制阈值的整体水平，后续以0.6为例。

基于上述提出的具有自适应形态阈值的形态掩膜(将v≥v_T的体素识别为动态，其余为静态)，对中间区域的混合信号进行分类，提取出其中的动态信号232。

动静态信号的BVSIM按照以下获得：

首先按照以下公式处理获得二值化体数据在位置(z，x，y)处的结构向量

其中B(α，z+h，x+i，y+j)表示二值化体数据中坐标(z+h，x+i，y+j)处的值，α为二值化数据对应的分割阈值线的角度，k表示计算结构向量的窗口大小，h、i和j表示窗口内像素的索引，(h，i，j)表示一个三维矢量，矢量的大小和方向由h、i和j决定。

然后按照以下公式计算组内各个二值化体数据的结构差异值之和，作为每个组的结构相似度值：

V＝∑_m，lΔv(m，l) (9)

其中，m、l分别表示某一组内二值化体数据序列的序号，V表示组内的每两幅二值化图像之间的图像结构相似度总和，Δv(m，l)表示第m个二值化体数据和第l个二值化体数据之间的结构差异度，|*|表示欧几里德距离，Z、X和Y分别是OCT数据深度方向、快扫描和慢扫描方向的像素数；

进而综合动静态各自组内体数据相似度处理获得整体的结构相似程度：

其中，V_d代表动态信号的组内结构相似度，V_s代表静态信号的组内结构相似度，n_d和n_s分别代表动态、静态信号组内二值化体数据的数量，

表示组合数，即从n_d/s个元素中选取2个元素的所有组合的个数。

4)遍历所有角度后，选取BVSIM最小的组合，根据其分类结果生成最终的二值化血管网络。

5)利用二值化血管网络图像进行血流形态量化包括：

在上述方法得到的二值化血管网络图像中，迭代删除血流外围像素，得到血流区域为单像素宽度的骨架图。再进行正面投影得到血流骨架图。在三维血流图上进行骨架提取的优势在于：更容易将将深度方向上重合的血流区分开。在二值化血管网络图像中，在横向和竖向两个方向上进行相邻两个像素的差分运算，得到血流边缘图。

计算的血流形态量化参数如下：

其中，n表示二值化血管网络图像的宽度和高度，(x，y)表示图像中的索引，A表示二值化血管网络图像，S表示血流骨架图，P表示作为血流边缘图。VDI反映了图像中血流平均管径。VSD为血流骨架图中血流所占的长度与总面积的比值，反应血流单位面积长度。VAD计算为血流面积与图像总面积的比值，反应血流面积密度。VPI为血流周长与图像总面积的比值，反应血流单位面积周长。

图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示出的是利用本发明的一个示例性实施例。结合血管形态特征的血流成像量化处理装置，包括宽带低相干光源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。

上述实验对比结果充分说明：利用本发明所涉及的基于多维度特征空间的光学相干血流造影方法，能够提升血流信号分类的准确率，实现血流对比度的有效增强，和血流图像质量的提升，具有其突出显著的技术效果。

由此本发明能显著抑制随机噪声的干扰，提高动静态信号的分类精度，改善二值化血管网络的连续性。

Claims (9)

1.一种结合血管形态特征的血流成像量化处理方法，其特征在于包括：

一种散射信号采集方式(1)，基于光学相干层析成像技术(OCT)采集三维空间内散射信号样本的OCT散射信号；

一种血流图像分割方法(2)，结合形态特征、OCT散射信号的信噪比倒数、去相关系数构建三维特征空间，实现动态血流信号与静态组织信号的分类，得到二值化血管网络图像；

一种微血管形态量化处理方法(3)，根据二值化血管网络图像进行血流骨架和轮廓提取获得血流骨架图，进而根据血流骨架图计算出反映血流形态的多种量化参数。

2.根据权利要求1所述的一种结合血管形态特征的血流成像量化处理方法，其特征在于：所述一种散射信号采集方式(1)，包括：对散射信号样本进行三维空间的OCT扫描成像，相同空间位置或其附近位置在T个不同时间点重复采样，采用以下方法之一：通过扫描改变参考臂光程的时间域OCT成像方法；利用光谱仪记录光谱干涉信号的光谱域OCT成像方法；利用扫频光源记录光谱干涉信号的扫频OCT成像方法。

3.根据权利要求1所述的一种结合血管形态特征的血流成像量化处理方法，其特征在于：所述的一种血流图像分割方法(2)，具体包括：

S1、采用一阶和零阶自协方差对OCT散射信号计算分析，得到各OCT散射信号的信噪比倒数和去相关系数两个特征，所得到的信噪比倒数和去相关系数进一步在三维空间、时间、角度及偏振态等多个维度上做滑动平均或高斯平均(21)，利用平均处理后的信噪比倒数和去相关系数的两个特征构建OCT散射信号的信噪比倒数-去相关系数(inverse SNR-decorrelation，ID)特征空间(22)；

S2、基于形态特征、信噪比倒数、去相关系数多维度特征空间对信号进行分类(23)，包括：遍历获得信噪比倒数-去相关系数特征空间中经过原点的两条线性分类边界，结合形态特征对三维空间的血管网络图像进行二值化处理获得二值化体数据，计算二值化体数据的结构相似度值，遍历所有线性分类边界的角度组合后，选取最小的结构相似度值对应的二值化结果作为最终的二值化血管网络(25)。

4.根据权利要求3所述的一种结合血管形态特征的血流成像量化处理方法，其特征在于：所述S2具体为：

遍历信噪比倒数-去相关系数特征空间中经过原点的每两条分割阈值线，通过两条分割阈值线将信噪比倒数-去相关系数特征空间划分为动态区域、中间区域和静态区域(231)；

借助形态特征构建形态掩膜提取出中间区域的动态信号(232)；

根据对动静态信号的分类结果，在三维空间中计算二值化体数据的结构相似度值(233)；

遍历每两条分割阈值线后，选取最小的结构相似度值对应的两条分割阈值线作为两条线性分类边界(24)，根据两条线性分类边界并结合形态掩膜区分中间区域的动静态信号，生成二值化血管网络(25)。

5.根据权利要求3或4所述的一种结合血管形态特征的血流成像量化处理方法，其特征在于：所述S2具体为：

S21、随机在信噪比倒数-去相关系数特征空间中建立经过原点的每两条分割阈值线，结合形态掩膜对信号实现初步的分类，分为初步静态信号和初步动态信号；

S22、先在信噪比倒数-去相关系数特征空间中对初步动态信号生成一系列经过原点的分割线，一系列分割线和去相关系数所在坐标轴之间的夹角逐渐增大，每两条分割线间包含1/n的总体素数目，利用分割线对动态区域进行二值化分割后获得一系列二值化体数据，将各个二值化体数据按照分割线的角度递增顺序组成一个序列作为初步动态信号的二值化体数据序列，计算初步动态区域内的体数据间的结构相似度，具体为：

首先，按照以下公式处理获得每个二值化体数据在位置(z，x，y)处的结构向量

其中，B(α，z+h，x+i，y+j)表示二值化体数据中坐标(z+h，x+i，y+j)处的值，α为二值化体数据对应的分割阈值线相对于去相关系数所在坐标轴的角度，k表示结构向量的窗口大小，h、i和j表示窗口内像素的三个坐标的索引，(h，i，j)表示一个三维矢量，三维矢量的大小和方向由h、i和j决定；

然后，按照以下公式计算各个二值化体数据的结构差异值之和，作为整个区域的结构相似度值：

其中，m、l分别表示区域内二值化体数据序列中的二值化体数据的序号，V表示区域内的每两个二值化体数据之间的图像结构相似度总和，即区域的结构相似度值，Δv(m，l)表示第m个二值化体数据和第l个二值化体数据之间的结构差异度，|*|表示欧几里德距离，Z、X和Y分别是OCT深度方向、快扫描方向和慢扫描方向的总像素数；

S23、按照和S22同样处理方式计算初步静态区域内的体数据间的结构相似度；

S24、综合动静态区域内的体数据间的结构相似度，得到最终的二值化体数据的结构相似度值，具体公式为：

其中，V_d代表动态区域的结构相似度，V_s代表静态区域的结构相似度，n_d和n_s分别代表动态区域、静态区域内的二值化体数据的数量，

表示从n_d个元素中选取2个元素的所有组合的个数，

表示从n_s个元素中选取2个元素的所有组合的个数。

6.根据权利要求1所述的一种微血管形态量化处理方法，其特征在于：

所述一种微血管形态量化处理方法(3)中，根据二值化血管网络图像进行血流骨架和轮廓提取获得血流骨架图，具体为：在二值化血管网络图像中沿水平面建立横向和竖向两个方向，在横向和竖向两个方向上分别对每两个相邻像素进行差分运算，进而得到血流边缘图；再迭代删除二值化血管网络图像中的血流区域外部像素，直到得到宽度为单像素的三维血流骨架，得到血流骨架图。

7.用于实施权利要求1～6任一所述方法的基于多维特征空间的微血流图像分割量化系统，包括：

一OCT光学相干层析探测装置，用于对三维空间内的散射信号样本的OCT散射信号进行采集；

一个图像处理器，用于用于获取并分析OCT散射信号的信噪比倒数和和去相关系数，并结合形态特征，进行动态血流信号与静态组织信号的分类，得到二值化血管网络图像；

一个数据处理器，根据二值化血管网络图像进行血流骨架和轮廓提取获得血流骨架图，进而根据血流骨架图计算出反应血流形态的多种量化参数，多种量化参数包括血流平均管径、血流面积密度、血流单位面积长度和血流单位面积周长。

8.根据权利要求7所述的基于多维特征空间的微血流图像分割量化系统，其特征在于：所述的一OCT光学相干层析探测装置是采用以下的一种：

包括低相干光源、干涉仪和探测器；

或者包括低相干光源、干涉仪和光谱仪；

或者包括扫频宽光谱光源、干涉仪和探测器。

9.根据权利要求7或8所述的基于多维特征空间的微血流图像分割量化系统，其特征在于：

所述的一OCT光学相干层析探测装置中选择地配置一个可见光指示装置，用于指示OCT探测光束的位置，指导探测目标的放置位置。

Images