本发明涉及一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，该方法包括如下步骤：S1、将光源接入3D内窥镜的n根导光束，并控制每个导光束依次通光形成一个照明周期，使用3D内窥镜的两个相机的近红外通道同步采集原始散斑图像，一个照明周期采集2n张原始散斑图像；S2、对于任意一个像素位置分别构建初始散射光强矩阵；S3、获取像素点的单次散射光强二阶中心矩、多次散射光强二阶中心矩以及多次散射光强一阶中心矩；S4、遍历所有像素点；S5、获取各像素点的单次散射衬比度值和多次散射衬比度值，按照像素点位置排列组成单次散射衬比度图像和多次散射衬比图像。与现有技术相比，本发明能够得到更为精准的表层和深部血流信息。
1.一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，其特征在于，该方法包括如下步骤：S1、将近红外相干激光光源接入3D内窥镜的n根导光束，并控制每个导光束依次通光形成一个照明周期，使用3D内窥镜的两个相机的近红外通道同步采集原始散斑图像，一个照明周期采集N＝2n张原始散斑图像；
S2、对于任意一个像素位置，采用T个照明周期采集的原始散斑图像的像素点灰度数据构建N×T大小的初始散射光强矩阵；
S3、基于初始散射光强矩阵获取对应像素点的单次散射光强二阶中心矩、多次散射光强二阶中心矩以及多次散射光强一阶中心矩；
S4、依次遍历原始散斑图像中的所有像素点，重复步骤S2～S3；
S5、基于单次散射光强二阶中心矩获取各个像素点的单次散射衬比度值，基于多次散射光强二阶中心矩以及多次散射光强一阶中心矩获取各个像素点的多次散射衬比度值，并按照像素点位置排列组成表征表层相对血流速度的单次散射衬比度图像以及表征深部相对血流速度的多次散射衬比图像。
2.根据权利要求1所述的一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，其特征在于，步骤S1具体为：将近红外相干激光光源接入3D内窥镜的n根导光束，3D内窥镜前端以固定位置分布n根导光束的出光口，不同导光束发光时与两个相机以及成像组织形成不同角度关系，在一个照明周期内导光束依次轮流发光，使用两个相机的近红外通道同步连续采集n帧原始散斑图像，单个照明周期记录到的总帧数为N＝2n。
3.根据权利要求1所述的一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，其特征在于，步骤S2中构建初始散射光强矩阵的具体方式包括：对于原始散斑图像的每个像素位置，将单个照明周期的N个成像角度的N张原始散斑图像的像素点灰度数据转化为一个列向量；
按时间顺序依次排列T个照明周期获得的列向量，构建出初始散射光强矩阵SH，SH的大小为N×T。
4.根据权利要求1所述的一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，其特征在于，步骤S3具体包括：S31、将初始散射光强矩阵转换为协方差矩阵；
S33、利用最小和次小的特征值以及 分布计算多次散射光强二阶中心矩；
S34、基于初始散射光强矩阵的二阶中心矩和多次散射光强二阶中心矩计算单次散射光强二阶中心矩；
S35、依据单次散射光强负指数分布特点求解单次散射光强一阶中心矩，进而求解多次散射光强一阶中心矩。
5.根据权利要求4所述的一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，其特征在于，步骤S31协方差矩阵通过下式获得：其中，MH为协方差矩阵，SH为初始散射光强矩阵， 矩阵与SH维度相同， 中的元素值为SH矩阵中对应位置元素所在行的均值， 为 的转置矩阵。
6.根据权利要求4所述的一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，其特征在于，步骤S32对协方差矩阵进行特征值分解表示为：T
其中，MH为协方差矩阵，Γ为对角阵，U和V是特征向量矩阵，Γ矩阵主对角线上的元素为协方差矩阵MH的特征值。
7.根据权利要求4所述的一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，其特征在于，步骤S33中多次散射光强二阶中心矩通过下式获得：其中， 为多次散射光强二阶中心矩，λN为协方差矩阵进行特征值分解得到的最小特征值，λN‑1为协方差矩阵进行特征值分解得到的次小特征值，α和β为正实数系数，α+β＝1，Q＝N/T。
8.根据权利要求4所述的一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，其特征在于，步骤S34中单次散射光强二阶中心矩通过下式获得：其中， 为单次散射光强二阶中心距， 为多次散射光强二阶中心矩， 为初始散射光强矩阵的二阶中心矩，γ为正实数系数，γ小于1。
9.根据权利要求4所述的一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，其特征在于，步骤S35单次散射光强一阶中心矩和多次散射光强一阶中心矩通过下式获得：其中， 为单次散射光强一阶中心矩， 为 的平方根， 为单次散射光强二阶中心距， 为多次散射光强一阶中心矩，μ为初始散射光强矩阵中像素点灰度数据的平均值。
10.根据权利要求1所述的一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，其特征在于，步骤S5中各个像素点的单次散射衬比度值和多次散射衬比度值通过下式获得：其中，KS为单次散射衬比度值， 为 的平方根， 为单次散射光强二阶中心距，KM为多次散射衬比度值， 为 的平方根， 为多次散射光强二阶中心矩， 为多次散射光强一阶中心矩。
基于多角度散射随机矩阵的3D内窥镜用血流成像方法
技术领域
[0001] 本发明涉及光学成像技术领域，尤其是涉及一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法。
背景技术
[0002] 3D内窥镜是外科手术中的一种新兴工具，可以提供实时深度感知。通过三维重建人体内部组织器官或病灶部位的轮廓数据，3D内窥镜有效提高了当前诊疗过程的准确性、可行性和科学性。目前，在腹腔、颅腔、消化内科等多种科室都已经引入使用了3D内窥镜。
[0003] 目前临床已经使用静脉注射吲哚菁绿染料并结合荧光成像获得肿瘤和脉管系统血流的可视化影像，但该技术受制于染料注射代谢周期影响，且易出现假阳性的结果。激光散斑衬比成像(laser speckle contrast imaging，LSCI)作为一种简单且低成本的方法，可以获取整个视野的2D血流灌注图，并实时提供血流变化的动态描述，因此可以用作术中血流监测工具。与荧光成像相比，LSCI无需造影剂，因此可以在手术期间的任何时候按需使用。反射式激光散斑系统记录到的信号同时包含单次散射和多次散射光子，分别对应术中对象表层和较深层血管的血流信息，目前传统的LSCI方法还无法实现单次和多次散射光强的统计分离。
[0004] 大部分医学临床诊断治疗不仅有对立体化可视的需求，也有对不同深度血管血流等生理信息显影的需求，如果能实现基于LSCI技术的单次和多次散射光强统计分离并与3D内窥镜有效结合，将显著提高3D内窥镜在精准手术、微创手术等应用中的作用。
发明内容
[0005] 本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法。
[0006] 本发明的目的可以通过以下技术方案来实现：
[0007] 一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，该方法包括如下步骤：
[0008] S1、将近红外相干激光光源接入3D内窥镜的n根导光束，并控制每个导光束依次通光形成一个照明周期，使用3D内窥镜的两个相机的近红外通道同步采集原始散斑图像，一个照明周期采集N＝2n张原始散斑图像；
[0009] S2、对于任意一个像素位置，采用T个照明周期采集的原始散斑图像的像素点灰度数据构建N×T大小的初始散射光强矩阵；
[0010] S3、基于初始散射光强矩阵获取对应像素点的单次散射光强二阶中心矩、多次散射光强二阶中心矩以及多次散射光强一阶中心矩；
[0011] S4、依次遍历原始散斑图像中的所有像素点，重复步骤S2～S3；
[0012] S5、基于单次散射光强二阶中心矩获取各个像素点的单次散射衬比度值，基于多次散射光强二阶中心矩以及多次散射光强一阶中心矩获取各个像素点的多次散射衬比度值，并按照像素点位置排列组成表征表层相对血流速度的单次散射衬比度图像以及表征深部相对血流速度的多次散射衬比图像。
[0013] 优选地，步骤S1具体为：
[0014] 将近红外相干激光光源接入3D内窥镜的n根导光束，3D内窥镜前端以固定位置分布n根导光束的出光口，不同导光束发光时与两个相机以及成像组织形成不同角度关系，在一个照明周期内导光束依次轮流发光，使用两个相机的近红外通道同步连续采集n帧原始散斑图像，单个照明周期记录到的总帧数为N＝2n。
[0015] 优选地，步骤S2中构建初始散射光强矩阵的具体方式包括：
[0016] 对于原始散斑图像的每个像素位置，将单个照明周期的N个成像角度的N张原始散斑图像的像素点灰度数据转化为一个列向量；
[0017] 按时间顺序依次排列T个照明周期获得的列向量，构建出初始散射光强矩阵SH，SH的大小为N×T。
[0018] 优选地，步骤S3具体包括：
[0019] S31、将初始散射光强矩阵转换为协方差矩阵；
[0020] S32、对协方差矩阵进行特征值分解；
[0021] S33、利用最小和次小的特征值以及 分布计算多次散射光强二阶中心矩；
[0022] S34、基于初始散射光强矩阵的二阶中心矩和多次散射光强二阶中心矩计算单次散射光强二阶中心矩；
[0023] S35、依据单次散射光强负指数分布特点求解单次散射光强一阶中心矩，进而求解多次散射光强一阶中心矩。
[0024] 优选地，步骤S31协方差矩阵通过下式获得：
[0025]
[0026]
[0027]
[0028] 其中，MH为协方差矩阵，SH为初始散射光强矩阵， 矩阵与SH维度相同， 中的元素值为SH矩阵中对应位置元素所在行的均值， 为 的转置矩阵。
[0029] 优选地，步骤S32对协方差矩阵进行特征值分解表示为：
[0030] MH＝UΓVT
[0031] 其中，MH为协方差矩阵，Γ为对角阵，U和V是特征向量矩阵，Γ矩阵主对角线上的元素为协方差矩阵MH的特征值。
[0032] 优选地，步骤S33中多次散射光强二阶中心矩通过下式获得：
[0033]
[0034] 其中， 为多次散射光强二阶中心矩，λN为协方差矩阵进行特征值分解得到的最小特征值，λN‑1为协方差矩阵进行特征值分解得到的次小特征值，α和β为正实数系数，α+β＝1，Q＝N/T。
[0035] 优选地，步骤S34中单次散射光强二阶中心矩通过下式获得：
[0036]
[0037] 其中， 为单次散射光强二阶中心距， 为多次散射光强二阶中心矩， 为初始散射光强矩阵的二阶中心矩，γ为正实数系数，γ小于1。
[0038] 优选地，步骤S35单次散射光强一阶中心矩和多次散射光强一阶中心矩通过下式获得：
[0039]
[0040]
[0041] 其中， 为单次散射光强一阶中心矩， 为 的平方根， 为单次散射光强二阶中心距， 为多次散射光强一阶中心矩，μ为初始散射光强矩阵中像素点灰度数据的平均值。
[0042] 优选地，步骤S5中各个像素点的单次散射衬比度值和多次散射衬比度值通过下式获得：
[0043]
[0044]
[0045] 其中，Ks为单次散射衬比度值， 为 的平方根， 为单次散射光强二阶中心距，KM为多次散射衬比度值， 为 的平方根， 为多次散射光强二阶中心矩， 为多次散射光强一阶中心矩。
[0046] 与现有技术相比，本发明具有如下优点：
[0047] (1)本发明对散射矩阵的特征描述更为准确，可以实现一阶中心矩和二阶中心矩的统计分离，从而实现单次散射和多次散射光强的统计分离，联合一阶矩和二阶矩的信息，可以更精确描述单次散射和多次散射的整体强度分布，并进一步获得更加准确的血流、组织灌注等功能信息。
[0048] (2)本发明可获取更为纯净的表层血管统计特性和血流信息，呈现更好的成像效果(单次散射部分)；可以呈现组织内较深部血管统计特性和血流信息，同时有效去除表面反光点(多次散射部分)；由于随机矩阵对刚性运动不敏感，因此本发明可以消除或降低部分生理性运动(呼吸、心跳等)产生的成像伪迹。
附图说明
[0049] 图1为本发明基于多角度散射随机矩阵的3D内窥镜用血流成像方法的流程框图；
[0050] 图2为本发明使用的3D内窥镜前端示意图。
具体实施方式
[0051] 下面结合附图和具体实施例对本发明进行详细说明。注意，以下的实施方式的说明只是实质上的例示，本发明并不意在对其适用物或其用途进行限定，且本发明并不限定于以下的实施方式。
[0052] 实施例
[0053] 如图1所示，本实施例提供一种基于多角度散射随机矩阵的3D内窥镜用血流成像方法，该方法包括如下步骤：
[0054] S1、将近红外相干激光光源接入3D内窥镜的n根导光束，并控制每个导光束依次通光形成一个照明周期，使用3D内窥镜的两个相机的近红外通道同步采集原始散斑图像，一个照明周期采集N＝2n张原始散斑图像。具体地：将近红外相干激光光源接入3D内窥镜的n根导光束，3D内窥镜前端以固定位置分布n根导光束的出光口，不同导光束发光时与两个相机以及成像组织形成不同角度关系，在一个照明周期内导光束依次轮流发光，使用两个相机的近红外通道同步连续采集n帧原始散斑图像，单个照明周期记录到的总帧数为N＝2n。本实施例中，近红外相干激光中心波长为830nm，功率为20mW，3D内窥镜设有4根导光束，并经过内窥镜的前端(如图2所示，呈对称分布)在主动控制下依次轮流将激光照射在被观察的术中对象上，来自于散射粒子(如红细胞)的背向散射光被内窥镜前端的两个(双目)相机的近红外通道同步连续采集(相机的采集帧速率和曝光时间与导光束依次照明时间同步)，单个照明周期(4根导光束依次轮流照明一次)，两个相机的近红外通道共采集N＝8帧原始散斑图像。
[0055] S2、对于任意一个像素位置，采用T个照明周期采集的原始散斑图像的像素点灰度数据构建N×T大小的初始散射光强矩阵，步骤S2中构建初始散射光强矩阵的具体方式包括：
[0056] 对于原始散斑图像的每个像素位置，将单个照明周期的N个成像角度的N张原始散斑图像的像素点灰度数据转化为一个列向量；
[0057] 按时间顺序依次排列T个照明周期获得的列向量，构建出初始散射光强矩阵SH，SH的大小为N×T；
[0058] 本实施例中，T＝24。
[0059] S3、基于初始散射光强矩阵获取对应像素点的单次散射光强二阶中心矩、多次散射光强二阶中心矩以及多次散射光强一阶中心矩；
[0060] 步骤S3具体包括：
[0061] S31、将初始散射光强矩阵转换为协方差矩阵，协方差矩阵通过下式获得：
[0062]
[0063]
[0064]
[0065] 其中，MH为协方差矩阵，SH为初始散射光强矩阵， 矩阵与SH维度相同， 中的元素值为SH矩阵中对应位置元素所在行的均值， 为 的转置矩阵。
[0066] S32、对协方差矩阵进行特征值分解，表示为：
[0067] MH＝UΓVT
[0068] 其中，MH为协方差矩阵，Γ为对角阵，U和V是特征向量矩阵，Γ矩阵其主对角线上的元素为协方差矩阵MH的特征值。
[0069] S33、利用最小和次小的特征值以及 分布计算多次散射光强二阶中心矩。由于3D内窥镜的两个(双目)相机采集到的背向散射光同时包含单次散射光和多次散射光，所以协方差矩阵MH可以表示为下式：
[0070] MH＝Ms+MM
[0071] 上式中，矩阵MS代表单次散射部分，矩阵MM代表多次散射部分。
[0072] 由于光在组织中多次散射的特性，矩阵MM为Wishart随机矩阵，符合分布，当T/N＝Q≥1时，其特征值的概率密度函数如下式：
[0073]
[0074] 上式中，λ±是矩阵WM特征值的上下边界，其与矩阵WM的方差 的关系如下式：
[0075]
[0076] 如果矩阵WM的四阶矩是有限的，当帧数T趋于无穷时，矩阵WM的最小特征值与下边界λ_近似相等；当不同成像角度的像素数N和帧数T均趋于无穷时，矩阵WH的最小特征值和矩阵WM的最小特征值近似相等。考虑到3D内窥镜的两个(双目)相机与导光束窗存在对称关系，需要利用最小和次小特征值和 分布计算多次散射光强二阶中心矩，进而多次散射光强二阶中心矩通过下式获得：
[0077]
[0078] 其中， 为多次散射光强二阶中心矩，λN为协方差矩阵进行特征值分解得到的最小特征值，λN‑1为协方差矩阵进行特征值分解得到的次小特征值，α和β为正实数系数，α+β＝1，Q＝N/T，本实施例中α＝0.5，β＝0.5，Q＝3。
[0079] S34、基于初始散射光强矩阵的二阶中心矩和多次散射光强二阶中心矩计算单次散射光强二阶中心矩。对矩阵MH的迹进行分解，如下式：
[0080]
[0081] 二阶中心矩和矩阵的迹关系如下式：
[0082]
[0083] 当不同成像角度的像素数N和帧数T均趋于无穷时，有如下关系式：
[0084]
[0085] 进而，可以单次散射光强二阶中心矩通过下式获得：
[0086]
[0087] 其中， 为单次散射光强二阶中心距， 为多次散射光强二阶中心矩， 为初始散射光强矩阵的二阶中心矩，γ为正实数系数，γ小于1。
[0088] S35、依据单次散射光强负指数分布特点求解单次散射光强一阶中心矩，进而求解多次散射光强一阶中心矩，具体通过下式获得：
[0089]
[0090]
[0091] 其中， 为单次散射光强一阶中心矩， 为 的平方根， 为单次散射光强二阶中心距， 为多次散射光强一阶中心矩，μ为初始散射光强矩阵中像素点灰度数据的平均值。
[0092] S4、依次遍历原始散斑图像中的所有像素点，重复步骤S2～S3；
[0093] S5、基于单次散射光强二阶中心矩获取各个像素点的单次散射衬比度值，基于多次散射光强二阶中心矩以及多次散射光强一阶中心矩获取各个像素点的多次散射衬比度值，并按照像素点位置排列组成表征表层相对血流速度的单次散射衬比度图像以及表征深部相对血流速度的多次散射衬比图像。其中，各个像素点的单次散射衬比度值和多次散射衬比度值通过下式获得：
[0094]
[0095]
[0096] 其中，KS为单次散射衬比度值， 为 的平方根， 为单次散射光强二阶中心距，KM为多次散射衬比度值， 为 的平方根， 为多次散射光强二阶中心矩， 为多次散射光强一阶中心矩。
[0097] 本发明可获取更为纯净的表层血管统计特性和血流信息，呈现更好的成像效果(单次散射部分)；可以呈现组织内较深部血管统计特性和血流信息，同时有效去除表面反光点(多次散射部分)；由于随机矩阵对刚性运动不敏感，因此本发明可以消除或降低部分生理性运动(呼吸、心跳等)产生的成像伪迹。相比传统基于单个相机的单次和多次散射光强统计分离方法，对散射矩阵的特征描述更为准确，可以实现一阶矩和二阶中心矩的统计分离，进而，联合一阶矩和二阶矩的信息，可以更精确描述单次散射和多次散射的整体强度分布，并进一步获得更加准确的血流、组织灌注等功能信息。
[0098] 上述实施方式仅为例举，不表示对本发明范围的限定。这些实施方式还能以其它各种方式来实施，且能在不脱离本发明技术思想的范围内作各种省略、置换、变更。
