CN111710012B - 一种基于两维复合配准的octa成像方法与装置 - Google Patents
一种基于两维复合配准的octa成像方法与装置 Download PDFInfo
- Publication number
- CN111710012B CN111710012B CN202010534378.5A CN202010534378A CN111710012B CN 111710012 B CN111710012 B CN 111710012B CN 202010534378 A CN202010534378 A CN 202010534378A CN 111710012 B CN111710012 B CN 111710012B
- Authority
- CN
- China
- Prior art keywords
- registration
- octa
- strip
- stripe
- registered
- 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
- 239000002131 composite material Substances 0.000 title claims abstract description 44
- 238000003384 imaging method Methods 0.000 title claims abstract description 44
- FCKYPQBAHLOOJQ-UHFFFAOYSA-N Cyclohexane-1,2-diaminetetraacetic acid Chemical compound OC(=O)CN(CC(O)=O)C1CCCCC1N(CC(O)=O)CC(O)=O FCKYPQBAHLOOJQ-UHFFFAOYSA-N 0.000 title claims abstract 38
- 238000000034 method Methods 0.000 claims abstract description 87
- 230000017531 blood circulation Effects 0.000 claims abstract description 65
- 230000033001 locomotion Effects 0.000 claims abstract description 60
- 210000004204 blood vessel Anatomy 0.000 claims abstract description 55
- 238000001514 detection method Methods 0.000 claims abstract description 21
- 238000002601 radiography Methods 0.000 claims abstract description 21
- 238000000605 extraction Methods 0.000 claims abstract description 13
- 238000012014 optical coherence tomography Methods 0.000 claims description 162
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 claims description 36
- 230000009466 transformation Effects 0.000 claims description 36
- 239000011159 matrix material Substances 0.000 claims description 32
- 238000013519 translation Methods 0.000 claims description 22
- 238000012545 processing Methods 0.000 claims description 18
- 238000006073 displacement reaction Methods 0.000 claims description 16
- 238000013507 mapping Methods 0.000 claims description 15
- 238000012216 screening Methods 0.000 claims description 13
- 230000003068 static effect Effects 0.000 claims description 12
- 230000002457 bidirectional effect Effects 0.000 claims description 10
- 238000011478 gradient descent method Methods 0.000 claims description 9
- 239000000126 substance Substances 0.000 claims description 8
- 238000005516 engineering process Methods 0.000 claims description 7
- 230000011218 segmentation Effects 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 5
- 238000010276 construction Methods 0.000 claims description 4
- 230000001419 dependent effect Effects 0.000 claims description 4
- 230000000717 retained effect Effects 0.000 claims description 3
- 230000002829 reductive effect Effects 0.000 abstract description 10
- 230000004927 fusion Effects 0.000 abstract description 5
- 230000000694 effects Effects 0.000 description 23
- 230000006870 function Effects 0.000 description 23
- 238000010586 diagram Methods 0.000 description 20
- 238000011156 evaluation Methods 0.000 description 9
- 230000003287 optical effect Effects 0.000 description 8
- 230000008569 process Effects 0.000 description 8
- 230000002207 retinal effect Effects 0.000 description 7
- 230000002792 vascular Effects 0.000 description 6
- 238000002583 angiography Methods 0.000 description 4
- 210000001525 retina Anatomy 0.000 description 4
- 230000003595 spectral effect Effects 0.000 description 4
- 230000000052 comparative effect Effects 0.000 description 3
- 238000012937 correction Methods 0.000 description 3
- 239000000284 extract Substances 0.000 description 3
- 230000004424 eye movement Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000000007 visual effect Effects 0.000 description 3
- 230000000670 limiting effect Effects 0.000 description 2
- 238000002310 reflectometry Methods 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 208000009795 Microphthalmos Diseases 0.000 description 1
- 206010044565 Tremor Diseases 0.000 description 1
- 238000001574 biopsy Methods 0.000 description 1
- 239000008280 blood Substances 0.000 description 1
- 210000004369 blood Anatomy 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000002591 computed tomography Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 206010020718 hyperplasia Diseases 0.000 description 1
- 230000003902 lesion Effects 0.000 description 1
- 239000012528 membrane Substances 0.000 description 1
- 201000010478 microphthalmia Diseases 0.000 description 1
- 210000004088 microvessel Anatomy 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000003607 modifier Substances 0.000 description 1
- 206010029864 nystagmus Diseases 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 238000011158 quantitative evaluation Methods 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000008439 repair process Effects 0.000 description 1
- 230000004599 slow eye movement Effects 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
- 238000000844 transformation Methods 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Images
Classifications
-
- 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/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/40—Image enhancement or restoration by the use of histogram techniques
-
- G06T5/80—
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
- G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR 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/10072—Tomographic images
- G06T2207/10101—Optical tomography; Optical coherence tomography [OCT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30101—Blood vessel; Artery; Vein; Vascular
- G06T2207/30104—Vascular flow; Blood flow; Perfusion
Abstract
本发明公开了一种基于两维复合配准的OCTA成像方法与装置。本发明通过一种两维复合配准的方法，融合多张初步血流造影图，在提取出OCTA有效条带后，复合使用基于大血管骨架配准度和SIFT特征点提取的全局配准方法以及基于SIFT流目标函数检测的局部配准方法来对OCTA有效条带对进行两维层面上的配准，进而消除运动伪影。本发明得到的融合OCTA结果图像的像素平均误差小、质量好、与原OCTA血流造影图的结构相似度高，且包含的运动伪影大幅减少甚至消失不见，实现高质量的OCTA血流造影成像。
Description
技术领域
本发明大体涉及生物医学成像领域，且更具体地涉及与光学相干层析成像(Optical Coherence Tomography，OCT)和光学相干血流成像(Optical CoherenceTomography Angiography，OCTA)相关联的运动伪影纠正与两维配准方法。
背景技术
自1991年时域Optical Coherence Tomography(OCT)被首次提出后，OCT技术作为一种非侵入性诊断工具逐渐被广泛应用在生物医学成像研究领域。不久之后，研究人员敏锐地感知到了OCT技术在视网膜检测领域的应用前景，为了得到高分辨率视网膜血流体积数据进而对视网膜血管网成像，Optical Coherence Tomography Angiography(OCTA)技术应运而生。作为一种新兴的活体组织成像技术，OCTA通过检测在视网膜相同轴向位置重复扫描得到的OCT信号反射率，并进行逐像素计算去相关信号来可视化视网膜血管网。多方研究结果显示，细节多、可视视场大的视网膜血管网成像非常有利于在临床医学等领域来诊断眼科疾病，辅助治疗。
但细节多、可视视场大的OCTA图像并不能简单地通过单次或多次扫描得到。在实际应用中，OCTA常常使用横向扫描(B-scan扫描)来提取OCT信号反射率的波动，进而得到一组可比较的双向预测帧(B-frame)。正是因为横向扫描速度的限制，在活体检测中不可避免的眼球运动会在一定程度上影响组织微血管的可视化和定量。这种非自愿性的眼球轴向运动和横向运动可以大致分为以下三类：相对缓慢的震颤、漂移和大范围快速移动的微眼跳。这三类眼球运动都会导致静态组织中去相关信号的增加，降低相邻组织之间的内在联系。其中，眼球震颤、漂移等较温和缓慢的眼球运动可以通过体部运动信号减运算和横向扫描的预配准来进行高效的抵消。但是，还没有成熟的方法来抵消微眼跳这类快速剧烈的运动，它们往往会给OCTA图像带来肉眼可见的运动伪影，严重影响视网膜血管网的构建。
为了消除处理掉这些肉眼可见的运动伪影，目前主要有基于硬件的和基于软件的两大类解决方案。基于硬件的方案往往会给OCT成像系统带来复杂性增加和探测效率大大降低的问题。所以，基于软件的解决方案是目前研究的主流方向，这一类方案大多通过两维或三维上的OCTA图像配准来人工构建拟合出新的OCTA数据或者用变形后的同一区域的无运动伪影的有效条带来填补运动伪影，形成细节更丰富，视场范围更大的OCTA图像。
因为两维配准方法要比三维配准方法更为简便、效率更高，所以我们主要针对使用横向扫描方式得到的初步血流造影图进行两维配准层面上的研究。近年来，已经有不少两维配准方法被提出，但都存在着局限性。在2017年Morgan Heisler等人提出了带状配准方法来配准OCTA图像。为了保证他们方法的稳定性，需要提前获取一张完整的无运动伪影的完美的初步血流造影图来作为模板贯穿整个配准过程。因为这种完美的初步血流造影图在实际应用中无法得到，所以他们配准方法受模板的质量限制，对模板的要求比较高。在2018年Ang Li等人提出一种基于微血管的检测和破口填塞的单一OCTA图像修补方法，该方法避免了复杂的图像采集过程，但此方法不能修补原图像数据中没有出现过的细节数据，局部细节很有可能存在缺失。
发明内容
本发明针对现有技术的不足，提出一种两维复合配准方法和装置来对使用横向扫描方式得到的包含运动伪影的初步血流造影图进行配准。
本发明的技术方案是：
一、一种基于两维复合配准的OCTA成像方法：
一种信号采集方法，利用OCT技术收集三维空间内散射信号样本；
一种运动信号分类方法，用于实现动态血流信号与静态组织的分类，生成初步血流造影图；
一种复合配准方法，用于融合多张初步血流造影图，消除运动伪影，获得高质量的OCTA血流造影成像。
所述一种信号采集方法，利用OCT技术收集三维空间内散射信号样本，包括：
对散射信号样本进行三维空间的OCT扫描成像，对于相同/相近空间区域及其附近位置，在快轴扫描方向进行n次重复扫描，在慢轴扫描方向进行m次重复扫描采样共计完成n+m次扫描。
所述对散射信号样本进行三维空间的OCT扫描成像，采用以下方式之一：通过扫描改变OCT扫描成像的参考臂光程的时间域OCT成像方法；利用光谱仪记录光谱干涉信号的光谱域OCT成像方法；利用扫频光源记录光谱干涉信号的扫频OCT成像方法。
所述一种运动信号分类方法，用于实现动态血流信号与静态组织的分类，生成初步血流造影图，包括：
结合OCT散射信号的信噪比与去相关系数构建二维信噪比倒数-去相关性系数特征空间，进行动态血流信号与静态组织的分类，获得初步血流造影图。
所述一种复合配准方法，用于融合多张初步血流造影图，消除运动伪影，具体包括：
S1：沿快轴扫描方向，从初步血流造影图中提取无运动伪影的OCTA有效条带；
S2：采用基于大血管骨架配准度和SIFT特征点提取的全局配准方法对待配准条带处理获得全局配准条带；
S21、从剩余其他未选择过的保留的OCTA有效条带中选取一个与基准条带有固定重合面积的条带作为待配准条带，以待配准条带和基准条带共同构成一对OCTA有效条带对，初始情况先随机选一个OCTA有效条带作为基准条带，也和任一选取的一个与基准条带有固定重合面积的OCTA有效条带共同构成初始OCTA有效条带对；
S22、针对OCTA有效条带对提取大血管骨架；
S23、提取OCTA有效条带对中成功匹配的SIFT特征点对；
S231、提取OCTA有效条带对中的SIFT特征点；
S232、使用最近/次近距离算法来匹配SIFT特征点对；
S233、利用尺度和方向信息筛选出成功匹配的SIFT特征点对；
S24、通过成功匹配的SIFT特征点对和大血管骨架配准度来进行全局配准；
S241、利用RANSAC随机数算法构建出与成功匹配的SIFT特征点对相对应的仿射变换矩阵，利用仿射变换矩阵对待配准条带n2进行初步刚性变形；
S242、处理获得大血管骨架配准度，并利用大血管骨架配准度来修正仿射变换矩阵的平移参数；
S243、利用修正后的仿射变换矩阵对待配准条带进行再次刚性变形来完成全局配准，得到全局配准条带；
S3：采用基于SIFT流目标函数检测的局部配准方法对全局配准条带处理获得局部配准条带；
S31、建立OCTA有效条带对中的SIFT流目标函数；(此处的OCTA有效条带对指基准条带和全局配准条带，但在S22和S23中指基准条带和待配准条带)
S32、通过SIFT流目标函数来改变OCTA图像中局部像素位置对全局配准条带进行局部配准得到局部配准条带；
S4：融合基准条带和局部配准条带来作为新的基准条带，回到步骤S2重复迭代进行S2至S3的步骤，直至所有OCTA有效条带均被选择处理过，以最后获得的基准条带作为高质量的OCTA血流造影成像。
所述S1具体为：
S11、计算单个利用横向扫描B-scan扫描得到双向预测帧B-frame的灰度；
S12、根据双向预测帧B-frame的灰度使用灰度峰值检测算法剔除OCTA无效条带，保留获得OCTA有效条带。
所述S22中，针对OCTA有效条带对提取大血管骨架，具体包括：
S221、对初步血流造影图先进行直方图均衡化和阈值分割预处理，得到预处理后OCTA图像；
S222、使用连通区域像素数筛选算法针对预处理后OCTA图像中OCTA有效条带对进行处理获得基准条带和待配准条带各自的大血管骨架
和
大血管骨架是指为一个图像区域，其中
和
分别表示基准条带和待配准条带的大血管骨架中的像素点(i,j)的像素值，连通区域像素数筛选算法具体如下：先计算预处理后OCTA图像的中OCTA有效条带连通区域的数目，如果数目达到预设数量阈值，则进一步计算并剔除像素数目少于k的连通区域，保留剩下连通区域；否则返回S221修改阈值分割的阈值来使连通区域的数目达到预设数量阈值。
所述S23，具体为：
S231：在OCTA有效条带对的基准条带或者待配准条带的图像区域范围内提取关键点描述子，然后把以关键点描述子为中心的连通邻近区域均匀划分为4×4的16块小区域，然后每个小区域生成8个方向的梯度直方图，最终所有梯度直方图组成具有16×8＝128维信息的SIFT特征点；SIFT特征点包含有尺度信息Scale(i,j)n、方向信息Orientation(i,j)n和二维位置信息(i,j)n；
S232：从基准条带中取一个SIFT特征点，然后再找到其与待配准条带中欧式距离最近的前两个SIFT特征点，按照以下方式判断：在这前两个SIFT特征点中，只有当欧式距离最近的SIFT特征点对除以欧式距离次近的SIFT特征点对少于预设距离倍数x，则接受欧式距离最近的SIFT特征点对作为初步匹配的SIFT特征点对；否则不作为初步匹配的SIFT特征点对；
S233：
利用SIFT特征点的尺度信息与方向信息对初步匹配的特征点对进行筛选，采用以下公式计算尺度因子和方向因子：
其中，
为基准条带n1和待配准条带n2之间的尺度因子，
为基准条带n1和待配准条带n2之间的方向因子；
然后将只有在同时符合尺度因子在[L1,L2]范围内以及方向因子在[θ1,θ2]范围内的特征点对，保留下来作为成功匹配的SIFT特征点对，L1,L2分别表示预期尺度因子的上限值和下限值，θ1,θ2分别表示预期方向因子的上限值和下限值。
所述S24，具体为：
S241：
利用RANSAC随机数算法对成功匹配的SIFT特征点对进行拟合运算得到仿射变换矩阵M，仿射变换矩阵M表示如下：
其中，α1、α2分别表示第一、第二刚性变形参数，β1、β2分别表示第三、第四刚性变形参数，γ1、γ2分别表示第一、第二平移参数；
然后将仿射变换矩阵M与待配准条带n2进行叉乘得到刚性变形后的待配准条带n2′；
S242：
构建大血管骨架配准度Cvessel如下所示：
Cvessel＝0～1
其中，
表示逻辑与运算函数，当基准条带n1和待配准条带n2对应位置像素点(i,j)处的大血管骨架重合相同时取为1，其余情况均取0；
表示刚性变形后的待配准条带n2′中的大血管骨架，
表示基准条带n1中的大血管骨架；Ck(i,j)表示像素点(i,j)的大血管骨架权重，k表示像素点(i,j)的序号，n表示像素点总数目；
然后使用大血管骨架配准度Cvessel来对第一、第二平移参数γ1、γ2进行修正，进而对仿射变换矩阵M的修正：把大血管骨架配准度Cvessel作为代价函数，第一、第二平移参数γ1和γ2是自变量，使用梯度下降法找到代价函数的全局最优解，刚性变形参数用α1、α2、β1、β2保持不变，并利用实施变化后的包含第一、第二平移参数γ1和γ2的仿射变换矩阵M来更新计算刚性变形后的待配准条带n2′中的大血管骨架，进而更新计算大血管骨架配准度Cvessel；通过梯度下降法中取得大血管骨架配准度Cvessel的最大值，将大血管骨架配准度Cvessel最大值对应的第一、第二平移参数γ1和γ2作为最优第一、第二平移参数γ1′和γ2′，最后把原来仿射变换矩阵M中的第一、第二平移参数γ1和γ2替换为最优第一、第二平移参数γ1′和γ2得到修改后的仿射变换矩阵M′。
S243：
把修正后的仿射变换矩阵M′与待配准条带n2再次进行叉乘获得全局配准条带n2″，达到全局配准的目的。
所述S3，具体为：
S31：
建立以下SIFT流目标函数E(ω)，具体表示如下：
ω(p)＝(u(p)+v(p))
其中，n1(p)、n2″(p)、n2″(p+ω(p))分别代表基准条带、全局配准条带和待求的局部配准条带，p是对应的OCTA有效条带中的像素点，p＝(i，j)，ω(p)代表在像素点p处的流向量，也就是像素p处的位移；u(p)和v(p)分别代表了二维空间中的像素点p处的水平流向量和竖直流向量，这里u(p)和v(p)只能取整数，因为是进行的是像素大小的位移；ε代表像素点p的邻域，是二维空间区域，q是邻域ε中的像素点；min(||n1(p)-n2″(p+ω(p))||，t)为数据项，t表示数据项的阈值；η(|u(p)+v(p)|)为微位移项，η表示微位移程度；min(α|u(p)-u(q)|，d)+min(α|v(p)-v(q)|，d)为平滑度项，d表示平滑度项的阈值，α表示平滑程度；min()表示求取最小值函数，|| ||、| |均表示取绝对值；
在数据项和平滑度项中分别有阈值t和d来控制配准失败点对和不连续度，微位移项中的η控制了像素p的微位移程度，平滑度项中的α控制了像素p相对于相邻像素q的位移的平滑程度。
本发明的局部配准借助SIFT流目标函数来专门修正全局配准条带中局部区域内的非刚性变形。
S32：
利用SIFT流目标函数E(ω)对全局配准条带进行处理，获得局部配准条带，具体为：使用梯度下降法计算SIFT流目标函数E(ω)对应应变量的流向量ω(p)的全局最优解，此时SIFT流目标函数E(ω)取得全局极小值，对应获得的n2″(p+ω(p))即为局部配准条带。
二、一种基于两维复合配准的OCTA成像装置，其特征在于：
一OCT光学相干层析探测装置，用于对三维空间内的散射信号样本的OCT散射信号进行采集；
一个或多个信号处理器，耦合到所述OCT光学相干层析扫描装置并适用于使所述OCT光学相干层析扫描装置：
实现动态血流信号与静态组织的分类；
生成初步血流造影图；
提取无运动伪影的OCTA有效条带；
基于大血管骨架配准度和SIFT特征点提取实现全局配准；
基于SIFT流目标函数检测实现局部配准；
融合多张初步血流造影图，消除运动伪影，获得高质量的OCTA血流造影成像。
本发明是通过一种两维复合配准的方法，融合多张初步血流造影图，消除运动伪影，获得高质量的OCTA血流造影成像。在从初步血流造影图中提取出OCTA有效条带后，使用两维复合配准方法对其进行配准。
该两维复合配准方法包括两部分：一部分是基于Scale Invariant FeatureTransform(SIFT)特征点提取和大血管骨架配准度的全局配准方法；另一部分是基于SIFT流目标函数检测的局部配准方法。配准时采用一系列限制条件对配准因子进行修正，并把融合后的新的基准条带作为下一次的配准模板，融合所有OCTA有效条带得到最终的无运动伪影的高质量的OCTA血流造影成像。对配准结果进行人眼观察和客观评价体系双重评估，结果表明：使用本发明得到的OCTA血流造影成像平均误差小、质量好、与原OCTA血流造影图的结构相似度高，运动伪影大幅减少甚至消失不见。
本发明对比已有技术具有以下显著优点：
(1)无需完美的无运动伪影的初步血流造影图模板。本论文从一组有运动伪影的初步血流造影图中直接提取出OCTA有效条带然后逐次进行配准融合，而不是采用常见的方法即：把条带与初步血流造影图模板进行配准。因为使用二维复合配准方法单次配准出的基准条带运动伪影少，不会影响后续配准过程的进行，所以突破了以往必须提前获取OCTA图像模板的两维配准方法的限制。在实际应用中考虑到效率的原因，医疗研究方只会对活体(如人眼)提取2到3组初步血流造影图，为了得到完美无运动伪影的初步血流造影图模板而去重复采样是不太现实的。所以本发明在一定程度上具有普遍性与实用性。
(2)在配准过程中融入了大血管配准度这一评价体系。对应于我们在观察OCTA图像时主要被视网膜大血管网络所吸引的视觉习惯，本发明构建了大血管骨架配准度的评价参数并把它融入到全局配准方法中用来修正仿射变换矩阵的平移参数，使得全局配准出的全局配准条带更加适合人眼观察。
(3)在面对大区域重合的图像的配准时，仿射变换以其简单和快速的计算效率闻名，但在处理有着复杂非刚性形变的图像的配准时得到的效果很差。对应的本发明的处理对象——初步血流造影图不仅包含刚性形变，还包含强烈的非刚性形变。本发明采用局部配准方法处理了非刚性形变，所用的时间是只采用仿射变换配准的两倍左右，同时像素水平上的误差下降了一半，减少了更多可见的运动伪影。综合来说，本发明具有优秀的配准质量和配准效率。
本发明得到的OCTA血流造影成像的像素平均误差小、质量好、与原OCTA血流造影图的结构相似度高，且包含的运动伪影大幅减少甚至消失不见。
附图说明
通过下面的详细描述连同附图，将更容易理解本发明的实施例程。在所有附图中，通过示例的方式而非限制的方式示出本发明的实施例。
图1为本发明方法的示意图；
图2中示出的是本发明提出的两维复合配准方法的示例性实施例的流程框图；
图3示意性地示出了横向扫描方式的说明图；
图4中给出的是双向预测帧是否位于运动伪影处的比较说明图；
图5中给出的是灰度峰值检测算法的处理效果对比图；
图6中示意性地示出了提取OCTA大血管骨架的效果对比图；
图7中示意性地示出了在OCTA中提取SIFT特征点的说明图；
图8中示出的是根据成功匹配的SIFT特征点对来对待配准条带进行仿射变换的得到的全局配准条带的效果图；
图9中示意性地示出了局部变形极其严重的待配准条带经历全局配准后得到的全局配准条带的效果对比图；
图10中示出的是对全局配准条带进行局部配准后的得到的局部配准条带的效果图；
图11为使用本发明提出的二维复合配准方法得到的实验效果图。
图中：1—为一种信号采集方法，利用OCT技术收集三维空间内散射信号样本；2—为一种运动信号分类方法，用于实现动态血流信号与静态组织的分类，生成初步血流造影图；3—为一种复合配准方法，用于融合多次OCTA成像结果，消除运动伪影。
更详细地，3包含：S1—沿快轴扫描方向，从初步血流造影图中提取无运动伪影的OCTA有效条带；S2—采用基于大血管骨架配准度和SIFT特征点提取的全局配准方法对待配准条带处理获得全局配准条带；S3—采用基于SIFT流目标函数检测的局部配准方法对全局配准条带处理获得局部配准条带；S4—融合基准条带和局部配准条带来作为新的基准条带，回到步骤S2重复迭代进行S2至S3的步骤，直至所有OCTA有效条带均被选择处理过，以最后获得的基准条带作为融合OCTA结果图像。
具体实施方式
下面将结合附图对本发明的具体实施例做详细说明，附图形成本发明的一部分。需要注意的是，这些说明及实施例仅仅为示例性，不能被理解为限制了本发明的范围，本发明的实施例的保护范围由附属的权利要求及其等同物限定，任何在本发明权利要求基础上的改动都是本发明的保护范围。
为了促进对本公开的实施例的审查，提供下面的特定术语的解释：
基准条带：在每次配准周期中不直接处理而作为基准的OCTA条带。在一周期的配准后通过融合基准条带与局部配准条带来更新，并在最后作为完美的无运动伪影的高质量的OCTA血流造影成像输出。
待配准条带：在每次配准周期中通过本发明提出的二维复合配准方法以基准条带作为基准来进行变形配准的OCTA条带。
全局/局部配准条带：在每次配准周期得到的经过全局/局部配准方法处理过的待配准/全局配准条带。
OCTA无效条带：长边沿快轴扫描方向的含大量运动伪影的OCTA条带。
OCTA有效条带：长边沿快轴扫描方向的不含运动伪影的OCTA条带。
运动伪影条带：在图像中因为目标物与相机之间的相对运动而导致的重叠阴影区域。因为横向扫描中沿快轴的扫描次数远多于慢轴的扫描次数，所以绝大多数运动伪影呈灰色条带状，且条带的长边沿快轴扫描方向。
全局配准：相对来说大范围的整体图像配准。在本发明中特指借助仿射变换来专门修正待配准条带中大范围区域内的刚性变形。
局部配准：相对来说小范围的局部图像配准。在本发明中特指借助SIFT流目标函数来专门修正全局配准条带中局部区域内的非刚性变形。
横向扫描B-scan：也即为截面断层扫描，通过单方向的轴向扫描获得。
双向预测帧B-frame：在横向扫描过程中获得的与光轴垂直的x-y截面。
大血管骨架：本发明是从图像学上来定义的，即血流造影图中出现的枝干粗、平均灰度高、最容易被观察到的视网膜血管骨架图像区域。
本发明的实施例如下：
将各操作描述成多个离散的操作有助于理解本发明的实施例；然而，描述的顺序不应被解释成隐含这些操作是依赖于顺序的。
各操作步骤中对应于计算因子、阈值因子等的数字比值不应被认定为本发明的实施例的唯一有效数值，而是作为实施例的操作示范数值。具体的数字比值应根据具体实施例的情况而进行规定才能发挥最好的效果。
实施例中可使用诸如“不含”、“有少量”等描述性的修饰词来表示对象间的比例关系。这种描述仅仅用于帮助理解本发明的实施例的操作对象的属性，不应被认定为表示绝对含义的数字量化比例。此外，如关于本发明的实施例所使用的术语如“包含”、“由……组成”等等均为同义的。
本描述中针对横向扫描样品空间采用基于空间方向的x-y-z三维坐标表示。这种描述仅仅用于促进讨论，而不意欲限制本发明的实施例的应用。其中：深度z方向为入射光轴的方向；x-y平面为垂直于光轴的平面。更详细地，x与y正交，且OCT快轴扫描方向沿x轴，慢轴扫描方向沿y轴。
本发明方法如图1所示，信号采集部分1表示对散射信号样本进行三维空间的OCT扫描成像，即对于相同空间区域及其附近位置，在快轴扫描方向x方向进行n次重复扫描，在慢轴扫描方向y方向进行m次重复扫描采样，共计完成n+m次扫描。采用的方式可以是以下任意一种：通过扫描改变OCT扫描成像的参考臂光程的时间域OCT成像方法；利用光谱仪记录光谱干涉信号的光谱域OCT成像方法；利用扫频光源记录光谱干涉信号的扫频OCT成像方法。
运动信号分类部分2表示用于实现动态血流信号与静态组织的分类，用于生成初步血流造影图。更详细地，第2部分具体包括：结合OCT散射信号的信噪比与去相关系数构建二维信噪比倒数-去相关性系数特征空间，进行动态血流信号与静态组织的分类，获得初步血流造影图。
复合配准方法3的目标是融合多张初步血流造影图，消除运动伪影。其包含三大操作步骤：1、沿快轴扫描方向，提取无运动伪影的OCTA有效条带；2、基于大血管骨架配准度和SIFT特征点提取的全局配准；3、基于SIFT流目标函数检测的局部配准。在上述操作步骤依次完成后，融合所有配准后的OCTA条带得到一个完整无运动伪影的OCTA图像。
图2中示出的是本发明提出的两维复合配准方法的示例性实施例的流程框图，实施例包含的所有操作步骤将在下面给出细致的说明。
S1：
图3中示意性地示出了横向扫描方式的说明图。在301中沿着光轴方向z前进，通过扫描相同位置处的x-y平面即得到一系列可比较的双向预测帧。而如302所示，对于同一处的双向预测帧有快轴扫描方向x和慢轴扫描方向y两种正交的扫描采样过程。但是由于沿快轴扫描方向x的扫描次数远多于慢轴扫描方向y的扫描次数，所以出现的运动伪影条带的长边基本是沿着快轴扫描方向的。
S11：
图4中给出的是双向预测帧是否位于运动伪影处的比较说明图，包含位于运动伪影处的双向预测帧图像401和不位于运动伪影处的双向预测帧图像402。由于双向预测帧是通过单方向行进的横向扫描方式得到的逐帧扫描两维截面图，如果在前进扫描时某处出现强烈的相对运动，相应位置处获得的双向预测帧中像素的平均灰度会有明显上升。通过401与402的比较，可以发现401的平均灰度远大于402，401中出现大片高灰度连通区域，符合上述说明。
S12：
进一步地，可以通过双向预测帧的平均灰度的变化幅度大小来判断OCTA无效条带位置并予以剔除。
本发明采用灰度峰值检测算法而不是传统的阈值检测方法，从而避免出现图像局部区域灰度突然上升而导致传统的阈值检测方法失效的情况。
图5中给出的是灰度峰值检测算法的处理效果对比图，包含初步血流造影图301以及处理效果图302。在301中的所有OCTA无效条带都被检测出来并予以剔除。二维x-y坐标图303展示了灰度峰值检测算法的核心思想，它以双向预测帧的平均灰度作为y轴纵坐标，横向扫描方向的前进位移量作为x轴横坐标。提取二维x-y坐标图的峰值点及其邻近帧作为OCTA无效条带并予以剔除，剩下的双向预测帧位置作为OCTA有效条带予以保留。采用峰值与距离峰值最近的极小值的比值来作为筛选峰值点的依据，针对实施例选用的初步血流造影图，本发明使用的峰值与距离峰值最近的极小值的比值为1.5。
另外，为了防止OCTA有效条带中出现高比例的运动伪影并且方便后续SIFT特征点的提取，如304方框中所示的宽度小于10个像素点的OCTA有效条带应予以剔除。
以上属于本发明提出的二维复合配准方法的操作步骤，目的是对OCTA初步成像结果进行预处理，提取OCTA有效条带、剔除OCTA无效条带。后面本发明将对二维复合配准方法的操作步骤，即全局配准方法进行说明。
S2：
S21、从剩余其他未选择过的保留的OCTA有效条带中选取一个与基准条带有一定的重合面积的条带作为待配准条带，以待配准条带和基准条带共同构成一对OCTA有效条带对，初始情况先随机选一个OCTA有效条带作为基准条带，也和任一选取的一个与基准条带有一定的重合面积的OCTA有效条带共同构成初始OCTA有效条带对；
S22：
提取OCTA有效条带对的大血管骨架的具体方法包含两大步骤：
S221、对初步血流造影图先进行直方图均衡化和阈值分割预处理，得到预处理后OCTA图像；
S222、使用连通区域像素数筛选算法针对预处理后OCTA图像中OCTA有效条带对进行处理获得基准条带和待配准条带各自的大血管骨架
和
大血管骨架是指为一个图像区域，其中
和
分别表示基准条带和待配准条带的大血管骨架中的像素点(i,j)的像素值，连通区域像素数筛选算法具体如下：先计算预处理后OCTA图像的中OCTA有效条带连通区域的数目，如果数目达到预设数量阈值，则进一步计算并剔除像素数目少于k的连通区域，保留剩下连通区域；否则返回S221修改阈值分割的阈值来使连通区域的数目达到预设数量阈值。
针对实施例选用的初步血流造影图，连通区域像素数筛选算法中使用的连通区域像素数的阈值为100。
图6中示意性地示出了提取OCTA大血管骨架的效果对比图，包含了初步血流造影图601、采用传统图像处理方法即反复迭代图像开闭运算得到的大血管骨架效果图602、以及采用本发明提出的方法得到的大血管骨架效果图603。可以清晰地看到，603的大血管骨架的细节更多、边缘轮廓更清晰。
基于上述实施例方法，本发明在此处提取大血管骨架是后续全局配准方法的基础，目的是减少肉眼可见的运动伪影，使得全局配准出的全局配准条带更加适合人眼观察。
S23：
S231：
在OCTA有效条带对的图像区域范围内进行提取SIFT特征点，图7中示意性地示出了在OCTA中提取SIFT特征点的说明图。
SIFT特征点的提取方法为：在OCTA有效条带对的基准条带或者待配准条带的图像区域范围内提取关键点描述子，然后把以关键点描述子为中心的连通邻近区域均匀划分为4×4的16块小区域，然后每个小区域生成8个方向的梯度直方图，最终所有梯度直方图组成具有16×8＝128维信息的SIFT特征点；SIFT特征点包含有尺度信息Scale(i,j)n、方向信息Orientation(i,j)n和二维位置信息(i,j)n。
具体实施的8个方向分别为上下左右的四个方向以及左上、左下、右上、右下的四个方向。
把提取出的SIFT特征点在原OCTA图像中标示出如效果图701所示，再从中随机筛选出SIFT特征点进行示例性呈现，得到效果图702后再对局部方形区域进行放大得到示意图703。正如703中标示的，SIFT特征点含有三类特征信息，以基准条带n1为例，由下列符号进行一一表示：
(1)
其中，
表示SIFT特征点在基准条带n1中的尺度信息710，其中i,j分别代表SIFT特征点在图像上的横、纵坐标。
(2)
其中，
表示SIFT特征点在基准条带n1中的方向信息720。
(3)
其中，
表示SIFT特征点在基准条带n1中的二维位置信息730。
S232：
在提取出OCTA有效条带对各自的SIFT特征点后，进一步地，需要对基准条带与待配准条带进行特征点匹配，匹配成功的SIFT特征点对才可作为后续全局配准的依据。首先使用最近/次近距离来对两条带进行SIFT特征点匹配，从基准条带中取一个SIFT特征点，然后再找到其与待配准条带中欧式距离最近的前两个SIFT特征点，按照以下方式判断：在这两个SIFT特征点中，只有当欧式距离最近的SIFT特征点对除以欧式距离次近的SIFT特征点对少于预设距离倍数x，则接受欧式距离最近的SIFT特征点对作为初步匹配的SIFT特征点对；否则不作为初步匹配的SIFT特征点对；
针对实施例选用的初步血流造影图，本发明采用x为0.8时得到的效果较好，每次配准能提取出大约70对特征点对。
S233：
在上述步骤获得的初步匹配的SIFT特征点对基础上，进一步剔除错误匹配的SIFT特征点对，本发明利用SIFT特征点的尺度信息710与方向信息720对初步匹配的特征点对进行筛选，计算尺度因子和方向因子：
(4)
其中，
为基准条带n1和待配准条带n2之间的尺度因子，
为基准条带n1和待配准条带n2之间的方向因子；
然后将只有在同时符合尺度因子在[L1,L2]范围内以及方向因子在[θ1,θ2]范围内的特征点对，保留下来作为成功匹配的SIFT特征点对，L1,L2分别表示符合预期的尺度因子的上限值和下限值，θ1,θ2分别表示符合预期的方向因子的上限值和下限值。
针对实施例选用的初步血流造影图，本发明采用[L1,L2]＝[0.5,2]，[θ1,θ2]＝[-π/16,π/16]得到的效果较好，每次配准过程能筛选出大约30对成功匹配的SIFT特征点对。
S24：
图8中示出的是根据成功匹配的SIFT特征点对来对待配准条带进行仿射变换的得到的全局配准条带的效果图，并针对有无使用大血管骨架配准度对仿射变换矩阵进行修正作了一个比较说明。其中，811与812分别为未使用大血管骨架配准度修正得到的融合伪彩色图和灰度图；821与822分别为使用大血管骨架配准度修正后得到的融合伪彩色图和灰度图。在使用大血管骨架配准度修正后，虚线方框区域内不仅重合叠影大大减少，而且视网膜血管网络的细节更加丰富，大幅提升了全局配准效果。
S241：
在获得基准条带与待配准条带之间的成功匹配的SIFT特征点对801后，利用RANSAC随机数算法对成功匹配的SIFT特征点对进行拟合运算得到仿射变换矩阵M，然后利用仿射变换矩阵M对待配准条带进行配准，仿射变换矩阵M表示如下：
(5)
其中，α1、α2分别表示第一、第二刚性变形参数，β1、β2分别表示第三、第四刚性变形参数，γ1、γ2分别表示第一、第二平移参数；使用α1、α2、β1、β2来对待配准条带进行刚性变形(包括旋转、剪切形变和尺度缩放等)，使用γ1、γ2来对待配准条带进行平移。
然后将仿射变换矩阵M与待配准条带n2进行叉乘以进行配准得到刚性变形后的待配准条带n2′；
S242：
构建大血管骨架配准度Cvessel如下所示：
(7)
Cvessel＝0～1
其中，
表示逻辑与运算函数，当基准条带n1和待配准条带n2应像素点(i，j)的大血管骨架重合相同时取为1，其余情况均取0；
表示刚性变形后的待配准条带n2′中的大血管骨架，
表示待配准条带n2中的大血管骨架；Ck(i，j)表示像素点(i，j)的大血管骨架权重，k表示像素点(i，j)的序号，n表示像素点总数目，数值不同表示在像素点(i，j)位置处已累积起来的大血管骨架数目不同；如若想得到更好的视觉效果，在像素点(i，j)位置处重合相同的大血管骨架数目越多的情况下，赋给Ck(i，j)的权值越高，Ck(i，j)的权值设置视实施例的具体情况而定。
然后使用大血管骨架配准度Cvessel来对第一、第二平移参数γ1、γ2进行修正，进而对仿射变换矩阵M的修正：把大血管骨架配准度Cvessel作为代价函数，第一、第二平移参数γ1和γ2是自变量，使用梯度下降法找到代价函数的全局最优解，刚性变形参数用α1、α2、β1、β2保持不变，并利用实施变化后的包含第一、第二平移参数γ1和γ2的仿射变换矩阵M来更新计算刚性变形后的待配准条带n2′中的大血管骨架，进而更新计算大血管骨架配准度Cvessel；通过梯度下降法中取得大血管骨架配准度Cvessel的最大值，将大血管骨架配准度Cvessel最大值对应的第一、第二平移参数γ1和γ2作为最优第一、第二平移参数γ1′和γ2′，最后把原来仿射变换矩阵M中的第一、第二平移参数γ1和γ2替换为最优第一、第二平移参数γ1′和γ2′得到修改后的仿射变换矩阵M′。
每配准完一对OCTA条带即要计算一次Cvessel，相应地，Ck矩阵也要更新一次。
针对实施例选用的初步血流造影图，具体实施在当Cvessel的值均小于0.75时，证明待配准条带的局部变形极其严重，进而需要对这个条带数据进行切分或者直接舍弃此条带的配准。
图9中示意性地示出了局部变形极其严重的待配准条带经历全局配准后得到的全局配准条带的效果对比图，包含大血管骨架融合伪彩色图901与融合伪彩色图902。虚线方框内严重的局部形变导致血管“增生”，但其他区域的全局配准效果基本未受影响。这种严重的局部形变极少出现，所以直接舍弃此待配准条带是一种高效率的选择。
S243：
把修正后的仿射变换矩阵M′与待配准条带n2再次进行叉乘获得全局配准条带n2″，达到全局配准的目的。
S3：
在得到的全局配准条带后，全局配准完毕。后面本发明将对二维复合配准方法的步骤3，即局部配准方法进行说明。
S31：
建立OCTA有效条带对中的SIFT流目标函数E(ω)，包含三部分：数据项、微位移项和平滑度项，具体形式表示如下：
(8)
ω(p)＝(u(p)+v(p))
其中，n1(p)、n2″(p)、n2″(p+ω(p))分别代表基准条带、全局配准条带和待求的局部配准条带，p是对应的OCTA有效条带中的像素点，p＝(i，j)，ω(p)代表在像素点p处的流向量，也就是像素p处的位移；u(p)和v(p)分别代表了二维空间中的像素点p处的水平流向量和竖直流向量，这里u(p)和v(p)只能取整数，因为是进行的是像素大小的位移；ε代表像素点p的邻域，是二维空间区域，q是邻域ε中的像素点；
min(||n1(p)-n2″(p+ω(p))||，t)为数据项，t表示数据项的阈值；
η(|u(p)+v(p)|)为微位移项，η表示微位移程度；
min(α|u(p)-u(q)|，d)+min(α|v(p)-v(q)|，d)为平滑度项，d表示平滑度项的阈值，α表示平滑程度；
在数据项和平滑度项中分别有阈值t和d来控制配准失败点对和不连续度，微位移项中的η控制了像素p的位移程度，平滑度项中的α控制了像素p相对于相邻像素q的位移的平滑程度。
本发明的局部配准借助SIFT流目标函数来专门修正全局配准条带中局部区域内的非刚性变形。
S32：
利用SIFT流目标函数E(ω)对全局配准条带进行处理，获得局部配准条带，具体为：使用梯度下降法计算E(ω)对应的应变量ω(p)的全局最优解，此时E(ω)取得全局极小值，对应的n2″(p+ω(p))即为局部配准条带。
图10中示出的是对全局配准条带进行局部配准后的得到的局部配准条带的效果图,1011、1012、1013分别为局部配准条带与基准条带的融合伪彩色图、大血管骨架伪彩色图和灰度图，1021、1022、1023分别为局部配准条带与基准条带的融合伪彩色图、大血管骨架伪彩色图和灰度图。通过观察效果图可得，经过局部配准后虚线方框内的重叠伪影大大减少，血管网络错位被修正，对比度上升。
至此，本发明提出的二维复合配准方法的三个主要操作步骤介绍完毕，并完整地经历了一个配准周期。
S4：融合基准条带和局部配准条带来作为新的基准条带，回到步骤S2重复迭代进行S2至S3的步骤，直至所有OCTA有效条带均被选择处理过，直至所有OCTA有效条带都遍历，以最后获得的基准条带作为成像结果，得到最后的目标——无运动伪影的融合OCTA结果图像。
本发明实施例的示例性实施情况：
使用本发明所描述的二维复合配准方法对一组人眼视网膜进行配准，目的是消除运动伪影，获得高质量的OCTA血流造影成像。实验在搭载Intel(R)Core(TM)i7-7700HQCPU@2.80GHz 2.80GHz处理器和8GB内存的电脑上运行，程序版本为MATLAB R2016a，单次全局配准的平均时间为3.157s，单次局部配准的平均时间为6.554s。
图11为使用本发明提出的二维复合配准方法得到的实验效果图。采用图2所示的流程，不断重复进行上述说明的配准周期，对所有OCTA有效条带依次进行配准可以得到效果图。1101、1102、1103作为一组初步血流造影图，选出了7个OCTA有效条带共经历了6个配准周期，经由二维复合配准方法处理后得到融合OCTA血流造影成像1110。可以清楚地观察到，相比于初步血流造影图，融合OCTA血流造影成像中所有的运动伪影已基本消失。
进一步地，为了定量描述二维复合配准系统的配准匹配度，并证明其强健性，本发明使用三个客观评价参数来对三个例示性实施例进行评估。其中，三个例示性实施例分别为：1、仅使用全局配准方法；2、仅使用局部配准方法；3、使用二维复合配准方法；每个实施例都采用相同的初步血流造影图，条带的配准顺序也完全一致。更详细地，客观评价参数分别为：Mean Squared Error(MSE)、Peak Signal to Noise Ratio(PSNR)和StructuralSimilarity Index for measuring image quality(SSIM)。
使用上述三个客观评价参数来对三个例示性实施例进行量化评估，得到的结果如表1所示：(其中所有客观评价参数均为平均值)
表1：例示性实施例的量化评估结果
例示性实施例3的MSE相比于其他例示性实施例显著减少，说明使用二维复合配准方法配准后，一对OCTA有效条带在配准后像素水平上的误差有较为明显的减少。更详细地，例示性实施例3的PSNR较高，表明图像的峰值信噪比有所提升，从侧面印证了噪声的减少，OCTA血流造影成像的对比度提升。另外，例示性实施例3的SSIM较高表示基准条带与局部配准条带的结构相似度高，适合人眼观察。
虽然各项结果表明二维复合配准方法是一种可以有效消除初步血流造影图中大量运动伪影的有效方案，但我们也承认这个方案在配准后得到的融合OCTA血流造影成像的质量上存在局限性。造成这种局限性有两方面原因，一方面是采集数据的限制，如果所有采集数据都在某一区域的质量偏低、细节偏少(如局部区域的病变)，那么融合图像的质量也很难有所提高，只有通过再采集几组数据来增强局部区域的对比度。另一方面是图像变形的参数估计的限制，如果采集到的初步血流造影图在局部区域的变形过于强烈，那么在全局配准时部分区域就会出现严重的运动伪影进而影响后面整体的配准步骤。虽然这两方面限制只在极少部分情况下出现，但考虑到二维复合配准方法的实际应用，我们还是不能忽视此方法的局限性。
上面所阐释的本公开包含多个实施例。虽然这些实施例中的每一个已经以它优选的形式公开，但是如本发明所公开和示出的特定实施例不以限定的方式看待，因为大量变化是可能的。但总的来说，本发明描述的基于两维复合配准的OCTA成像方法与装置是一种可以有效消除运动伪影并得到高质量的OCTA血流造影成像的发明。
Claims (8)
1.一种基于两维复合配准的OCTA成像方法，包括：
一种信号采集方法(1)，利用OCT技术收集三维空间内散射信号样本；
一种运动信号分类方法(2)，用于实现动态血流信号与静态组织的分类，生成初步血流造影图；
一种复合配准方法(3)，用于融合多张初步血流造影图，消除运动伪影，获得高质量的OCTA血流造影成像；
所述一种复合配准方法(3)，用于融合多张初步血流造影图，消除运动伪影，具体包括：
S1：沿快轴扫描方向，从初步血流造影图中提取无运动伪影的OCTA有效条带；
S2：采用基于大血管骨架配准度和SIFT特征点提取的全局配准方法对待配准条带处理获得全局配准条带；
S21、从剩余其他未选择过的OCTA有效条带中选取一个与基准条带有重合面积的条带作为待配准条带，以待配准条带和基准条带共同构成一对OCTA有效条带对，初始情况先随机选一个OCTA有效条带作为基准条带；
S22、针对OCTA有效条带对提取大血管骨架；
S23、提取OCTA有效条带对中成功匹配的SIFT特征点对；
S231、提取OCTA有效条带对中的SIFT特征点；
S232、使用最近/次近距离算法来匹配SIFT特征点对；
S233、利用尺度和方向信息筛选出成功匹配的SIFT特征点对；
S24、通过成功匹配的SIFT特征点对和大血管骨架配准度来进行全局配准；
所述S24，具体为：
S241：
利用RANSAC随机数算法对成功匹配的SIFT特征点对进行拟合运算得到仿射变换矩阵M，仿射变换矩阵M表示如下：
(5)
其中，α1、α2分别表示第一、第二刚性变形参数，β1、β2分别表示第三、第四刚性变形参数，γ1、γ2分别表示第一、第二平移参数；
然后将仿射变换矩阵M与待配准条带n2进行叉乘得到刚性变形后的待配准条带n2′；
S242：
构建大血管骨架配准度Cvessel如下所示：
(7)
Cvessel＝0～1
其中，表示逻辑与运算函数，当基准条带n1和待配准条带n2对应位置像素点(i，j)处的大血管骨架重合相同时取为1，其余情况均取0；表示刚性变形后的待配准条带n2′中的大血管骨架，表示基准条带n1中的大血管骨架；Ck(i，j)表示像素点(i，j)的大血管骨架权重，k表示像素点(i，j)的序号，n表示像素点总数目；
然后使用大血管骨架配准度Cvessel来对第一、第二平移参数γ1、γ2进行修正，进而对仿射变换矩阵M的修正：把大血管骨架配准度Cvessel作为代价函数，第一、第二平移参数γ1和γ2是自变量，使用梯度下降法找到代价函数的全局最优解，并利用实施变化后的包含第一、第二平移参数γ1和γ2的仿射变换矩阵M来更新计算刚性变形后的待配准条带n2′中的大血管骨架，进而更新计算大血管骨架配准度Cvessel；通过梯度下降法中取得大血管骨架配准度Cvessel的最大值，将大血管骨架配准度Cvessel最大值对应的第一、第二平移参数γ1和γ2作为最优第一、第二平移参数γ1′和γ2′，最后把原来仿射变换矩阵M中的第一、第二平移参数γ1和γ2替换为最优第一、第二平移参数γ1′和γ2′得到修改后的仿射变换矩阵M′；
S243：把修正后的仿射变换矩阵M′与待配准条带n2再次进行叉乘获得全局配准条带n2″，达到全局配准的目的；
S3：采用基于SIFT流目标函数检测的局部配准方法对全局配准条带处理获得局部配准条带；
S31、建立OCTA有效条带对中的SIFT流目标函数；
S32、通过SIFT流目标函数来改变OCTA图像中局部像素位置对全局配准条带进行局部配准得到局部配准条带；
S4：融合基准条带和局部配准条带来作为新的基准条带，回到步骤S2重复迭代进行S2至S3的步骤，直至所有OCTA有效条带均被选择处理过，以最后获得的基准条带作为高质量的OCTA血流造影成像。
2.根据权利要求1所述的一种基于两维复合配准的OCTA成像方法，其特征在于：所述一种信号采集方法(1)，利用OCT技术收集三维空间内散射信号样本，包括：对散射信号样本进行三维空间的OCT扫描成像，对于相同/相近空间区域，在快轴扫描方向进行n次重复扫描，在慢轴扫描方向进行m次重复扫描采样。
3.根据权利要求1所述的一种基于两维复合配准的OCTA成像方法，其特征在于：所述一种运动信号分类方法(2)，用于实现动态血流信号与静态组织的分类，生成初步血流造影图，包括：
结合OCT散射信号的信噪比与去相关系数构建二维信噪比倒数和去相关性系数特征空间，进行动态血流信号与静态组织的分类，获得初步血流造影图。
4.根据权利要求1所述的一种基于两维复合配准的OCTA成像方法，其特征在于：所述S1具体为：
S11、计算单个利用横向扫描B-scan扫描得到双向预测帧B-frame的灰度；
S12、根据双向预测帧B-flame的灰度使用灰度峰值检测算法剔除OCTA无效条带，保留获得OCTA有效条带。
5.根据权利要求1所述的一种基于两维复合配准的OCTA成像方法，其特征在于：所述S22中，针对aCTA有效条带对提取大血管骨架，具体包括：
S221、对初步血流造影图先进行直方图均衡化和阈值分割预处理，得到预处理后OCTA图像；
S222、使用连通区域像素数筛选算法针对预处理后OCTA图像中OCTA有效条带对进行处理获得基准条带和待配准条带各自的大血管骨架和其中和分别表示基准条带和待配准条带的大血管骨架中的像素点(i，j)的像素值，连通区域像素数筛选算法具体如下：先计算预处理后OCTA图像的中OCTA有效条带连通区域的数目，如果数目达到预设数量阈值，则进一步计算并剔除像素数目少于k的连通区域，k表示预设的阈值，保留剩下连通区域；否则返回S221修改阈值分割的阈值来使连通区域的数目达到预设数量阈值。
6.根据权利要求1所述的一种基于两维复合配准的OCTA成像方法，其特征在于：所述S23，具体为：
S231：在OCTA有效条带对的基准条带或者待配准条带的图像区域范围内提取关键点描述子，然后把以关键点描述子为中心的连通邻近区域均匀划分为4×4的16块小区域，然后每个小区域生成8个方向的梯度直方图，最终所有梯度直方图组成具有16×8＝128维信息的SIFT特征点；SIFT特征点包含有尺度信息Scale(i，j)n、方向信息Orientation(i，j)n和二维位置信息(i，j)n；
S232：从基准条带中取一个SIFT特征点，然后再找到其与待配准条带中欧式距离最近的前两个SIFT特征点，按照以下方式判断：在这前两个SIFT特征点中，只有当欧式距离最近的SIFT特征点对除以欧式距离次近的SIFT特征点对少于预设距离倍数x，则接受欧式距离最近的SIFT特征点对作为初步匹配的SIFT特征点对；否则不作为初步匹配的SIFT特征点对；
S233：
利用SIFT特征点的尺度信息与方向信息对初步匹配的特征点对进行筛选，采用以下公式计算尺度因子和方向因子：
其中，为基准条带n1和待配准条带n2之间的尺度因子，为基准条带n1和待配准条带n2之间的方向因子；
然后将只有在同时符合尺度因子在[L1，L2]范围内以及方向因子在[θ1，θ2]范围内的特征点对，保留下来作为成功匹配的SIFT特征点对，L1，L2分别表示预期尺度因子的上限值和下限值，θ1，θ2分别表示预期方向因子的上限值和下限值。
7.根据权利要求1所述的一种基于两维复合配准的OCTA成像方法，其特征在于：所述S3，具体为：
S31：
建立以下SIFT流目标函数E(ω)，具体表示如下：
(8)
ω(p)＝(u(p)+v(p))
其中，n1(p)、n2″(p)、n2″(p+ω(p))分别代表基准条带、全局配准条带和待求的局部配准条带，p是对应的OCTA有效条带中的像素点，p＝(i，j)，ω(p)代表在像素点p处的流向量，u(p)和v(p)分别代表了二维空间中的像素点p处的水平流向量和竖直流向量，ε代表像素点p的邻域，q是邻域ε中的像素点；min(||n1(p)-n2″(p+ω(p))||，t)为数据项，t表示数据项的阈值；η(|u(p)+v(p)|)为微位移项，η表示微位移程度；min(α|u(p)-u(q)|，d)+min(α|v(p)-v(q)|，d)为平滑度项，d表示平滑度项的阈值，α表示平滑程度；min()表示求取最小值函数，|| ||、| |均表示取绝对值；
S32：
利用SIFT流目标函数E(ω)对全局配准条带进行处理，获得局部配准条带，具体为：使用梯度下降法计算SIFT流目标函数E(ω)对应应变量的流向量ω(p)的全局最优解，此时SIFT流目标函数E(ω)取得全局极小值，对应获得的n2″(p+ω(p))即为局部配准条带。
8.用于实施权利要求1～7任一所述方法的一种基于两维复合配准的OCTA成像装置，包括：
一OCT光学相干层析探测装置，用于对三维空间内的散射信号样本的OCT散射信号进行采集；
一个或多个信号处理器，耦合到所述OCT光学相干层析扫描装置并适用于使所述OCT光学相干层析扫描装置：
实现动态血流信号与静态组织的分类；
生成初步血流造影图；
提取无运动伪影的OCTA有效条带；
基于大血管骨架配准度和SIFT特征点提取实现全局配准；
基于SIFT流目标函数检测实现局部配准；
融合多张初步血流造影图，消除运动伪影，获得高质量的OCTA血流造影成像。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010534378.5A CN111710012B (zh) | 2020-06-12 | 2020-06-12 | 一种基于两维复合配准的octa成像方法与装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010534378.5A CN111710012B (zh) | 2020-06-12 | 2020-06-12 | 一种基于两维复合配准的octa成像方法与装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111710012A CN111710012A (zh) | 2020-09-25 |
| CN111710012B true CN111710012B (zh) | 2023-04-14 |
Family
ID=72540242
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010534378.5A Active CN111710012B (zh) | 2020-06-12 | 2020-06-12 | 一种基于两维复合配准的octa成像方法与装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111710012B (zh) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112396622B (zh) * | 2020-11-24 | 2023-10-31 | 浙江大学 | 基于多维特征空间的微血流图像分割量化方法和系统 |
| CN112819867A (zh) * | 2021-02-05 | 2021-05-18 | 苏州大学 | 基于关键点匹配网络的眼底图像配准方法 |
| CN113469986A (zh) * | 2021-07-13 | 2021-10-01 | 深圳市中科微光医疗器械技术有限公司 | 图像处理的方法、装置、电子设备及存储介质 |
| CN114066889B (zh) * | 2022-01-12 | 2022-04-29 | 广州永士达医疗科技有限责任公司 | 一种oct主机的成像质量检测方法及装置 |
| CN115496917B (zh) * | 2022-11-01 | 2023-09-26 | 中南大学 | 一种GPR B-Scan图像中的多目标检测方法及装置 |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101604444A (zh) * | 2009-07-09 | 2009-12-16 | 复旦大学 | 用图像质量评估作为相似性测度的血管减影图像配准方法 |
| CN101822545A (zh) * | 2010-05-11 | 2010-09-08 | 河南大学 | 一种数字减影血管造影运动伪影消除方法及其系统 |
| CN102722890A (zh) * | 2012-06-07 | 2012-10-10 | 内蒙古科技大学 | 基于光流场模型的非刚性心脏图像分级配准方法 |
| US8401276B1 (en) * | 2008-05-20 | 2013-03-19 | University Of Southern California | 3-D reconstruction and registration |
| CN103606152A (zh) * | 2013-11-15 | 2014-02-26 | 大连理工大学 | 基于sift特征点聚类及布尔差运算的dsa血管图像分割方法 |
| CN106934761A (zh) * | 2017-02-15 | 2017-07-07 | 苏州大学 | 一种三维非刚性光学相干断层扫描图像的配准方法 |
| CN107886508A (zh) * | 2017-11-23 | 2018-04-06 | 上海联影医疗科技有限公司 | 差分减影方法和医学图像处理方法及系统 |
| CN111260543A (zh) * | 2020-01-19 | 2020-06-09 | 浙江大学 | 一种基于多尺度图像融合和sift特征的水下图像拼接方法 |
-
2020
- 2020-06-12 CN CN202010534378.5A patent/CN111710012B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8401276B1 (en) * | 2008-05-20 | 2013-03-19 | University Of Southern California | 3-D reconstruction and registration |
| CN101604444A (zh) * | 2009-07-09 | 2009-12-16 | 复旦大学 | 用图像质量评估作为相似性测度的血管减影图像配准方法 |
| CN101822545A (zh) * | 2010-05-11 | 2010-09-08 | 河南大学 | 一种数字减影血管造影运动伪影消除方法及其系统 |
| CN102722890A (zh) * | 2012-06-07 | 2012-10-10 | 内蒙古科技大学 | 基于光流场模型的非刚性心脏图像分级配准方法 |
| CN103606152A (zh) * | 2013-11-15 | 2014-02-26 | 大连理工大学 | 基于sift特征点聚类及布尔差运算的dsa血管图像分割方法 |
| CN106934761A (zh) * | 2017-02-15 | 2017-07-07 | 苏州大学 | 一种三维非刚性光学相干断层扫描图像的配准方法 |
| CN107886508A (zh) * | 2017-11-23 | 2018-04-06 | 上海联影医疗科技有限公司 | 差分减影方法和医学图像处理方法及系统 |
| CN111260543A (zh) * | 2020-01-19 | 2020-06-09 | 浙江大学 | 一种基于多尺度图像融合和sift特征的水下图像拼接方法 |
Non-Patent Citations (2)
| Title |
|---|
| 王安娜 等.基于SIFT 特征提取的非刚性医学图像配准算法研究.生物医学工程学杂志.2010,第27卷(第4期),P763-P784. * |
| 赵明.基于改进SIFT 特征的红外与可见光图像配准方法.光电工程.2011,第38卷(第9期),P130-P136. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111710012A (zh) | 2020-09-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111710012B (zh) | 一种基于两维复合配准的octa成像方法与装置 | |
| Gao et al. | A deep learning based approach to classification of CT brain images | |
| CN110163809A (zh) | 基于U-net生成对抗网络DSA成像方法及装置 | |
| Ricco et al. | Correcting motion artifacts in retinal spectral domain optical coherence tomography via image registration | |
| Mahapatra et al. | Retinal vasculature segmentation using local saliency maps and generative adversarial networks for image super resolution | |
| US9089280B2 (en) | Image processing apparatus, image processing method, and program storage medium | |
| CN110717956A (zh) | 一种有限角投影超像素引导的l0范数最优化重建方法 | |
| CN112164043A (zh) | 一种用于多眼底图像拼接的方法及系统 | |
| CN112562058B (zh) | 一种基于迁移学习的颅内血管模拟三维模型快速建立方法 | |
| CN114092405A (zh) | 一种针对黄斑水肿oct图像的视网膜层自动分割方法 | |
| CN108665474B (zh) | 一种基于b-cosfire的眼底图像视网膜血管分割方法 | |
| Xie et al. | Speckle denoising of optical coherence tomography image using residual encoder–decoder CycleGAN | |
| Niemeijer et al. | Automated localization of the optic disc and the fovea | |
| CN110033496B (zh) | 时间序列三维视网膜sd-oct图像的运动伪差校正方法 | |
| Lin et al. | Res-UNet based optic disk segmentation in retinal image | |
| CN113160050B (zh) | 基于时空神经网络的小目标识别方法及系统 | |
| CN115393239A (zh) | 一种多模态眼底图像配准和融合方法及系统 | |
| CN113012184B (zh) | 基于Radon变换与多类型图像联合分析的微血管瘤检测方法 | |
| CN115205241A (zh) | 一种用于视细胞密度的计量方法及系统 | |
| Pohle-Fröhlich et al. | Estimation of muscle fascicle orientation in ultrasonic images | |
| Kumar et al. | Semiautomatic method for segmenting pedicles in vertebral radiographs | |
| CN114240893A (zh) | 一种外像图像中脊柱Cobb角的测算方法 | |
| Ahdi et al. | A hybrid method for 3D mosaicing of OCT images of macula and optic nerve head | |
| CN113112473A (zh) | 一种人体扩张型心肌病的自动诊断系统 | |
| CN112102327B (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 |