[go: up one dir, main page]

CN114359357B - 非刚性图像配准方法及其系统和用途 - Google Patents

非刚性图像配准方法及其系统和用途

Info

Publication number
CN114359357B
CN114359357B CN202111642914.4A CN202111642914A CN114359357B CN 114359357 B CN114359357 B CN 114359357B CN 202111642914 A CN202111642914 A CN 202111642914A CN 114359357 B CN114359357 B CN 114359357B
Authority
CN
China
Prior art keywords
image
sub
feature points
registration method
feature
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
CN202111642914.4A
Other languages
English (en)
Other versions
CN114359357A (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.)
Peking University
Original Assignee
Peking 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 Peking University filed Critical Peking University
Priority to CN202111642914.4A priority Critical patent/CN114359357B/zh
Publication of CN114359357A publication Critical patent/CN114359357A/zh
Application granted granted Critical
Publication of CN114359357B publication Critical patent/CN114359357B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; 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 OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/06—Physical realisation, i.e. hardware implementation of neural networks, neurons or parts of neurons
    • G06N3/067—Physical realisation, i.e. hardware implementation of neural networks, neurons or parts of neurons using optical means
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/02—Affine transformations
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/10—Image enhancement or restoration using non-spatial domain filtering
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G06T2207/20056—Discrete and fast Fourier transform, [DFT, FFT]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Artificial Intelligence (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Neurology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Analysis (AREA)

Abstract

本发明涉及非刚性图像配准方法及系统，特别是基于三角剖分和分段仿射变换的非刚性图像配准方法及系统，属于图像处理和分析领域。所述非刚性图像配准方法包括：将特征图像进行剖分，得到三角形图像块，对三角形图像块进行分段仿射变换，以及合并所有分段仿射变换后的三角形图像块。本发明将分段配准的三角形图像块组合起来，高效和高准确度地实现了跨多个时程的神经环路成像配准。

Description

非刚性图像配准方法及其系统和用途
技术领域
本发明涉及非刚性图像配准方法及其系统和用途，特别是基于三角剖分和分段仿射变换的非刚性图像配准方法及其系统和用途，属于图像处理和分析领域。
背景技术
探索神经环路的特性是了解大脑认知能力的关键。神经环路成像技术可以在多日内以单神经元分辨率监测数百至数千个神经元的活动。长期记录的多时程神经成像是研究神经环路功能机制的重要技术手段，准确有效的非刚性图像配准对于长期记录的多时程神经环路成像至关重要。但是，在不同时程的成像过程中，由于实验操作过程中显微镜重定位的位置误差，以及大脑的非刚性变形，不同时程中的视野通常无法精确匹配。
在分析多时程神经环路成像数据之前，需要将不同时程中的相同特征进行位置配准。图像配准广泛用于许多图像分析任务，例如计算机视觉和医学图像。典型的图像配准过程包括四个关键步骤：特征检测、特征匹配、映射函数设计和图像变换。在神经环路成像配准中，发挥重要作用的关键图像特征包括神经元、轴突、树突和血管。大多数已经建立的神经图像配准工具都是为了解决同一时程成像结果中小变形的配准问题而开发的，跨多个时程的神经环路成像配准仍然是一个挑战。传统的图像刚性仿射变换具有完善的理论基础，但不适合直接解决大脑的非刚性变形问题。而且，多天多时程神经环路成像具有较大的复杂变形，传统的刚性变形配准方法和适用于小范围非刚性变形配准的方法，不适用于多天内多时程的神经环路成像配准，不适合直接解决大脑的非刚性变形问题。
发明内容
为改善上述技术问题，本发明提供一种非刚性图像配准方法及其系统和用途，特别是基于三角剖分和分段仿射变换的非刚性图像配准方法及其系统和用途。
根据本发明的实施方案，所述非刚性图像配准方法包括：将特征图像进行剖分，得到三角形图像块，对三角形图像块进行分段仿射变换，以及合并所有分段仿射变换后的三角形图像块。
根据本发明的实施方案，所述非刚性图像配准方法包括以下步骤：
S1：从多时程神经环路成像结果中提取特征图像；
S2：从特征图像中识别局部特征点；
S3：将局部特征点进行配对；
S6：向外扩展图像，增加扩展图像边界特征点；
S7：根据配对的特征点，将特征图像进行三角剖分，得到三角形图像块；
S8：对每一个三角形图像块进行分段仿射变换；
S9：合并所有分段仿射变换后的三角形图像块，得到最终的非刚性图像配对结果。
根据本发明的实施方案，所述配准方法的步骤S2和/或S3可以自动进行或非自动进行。
根据本发明的实施方案，当所述配准方法的步骤S2自动进行时，所述配准方法还可以包括如下步骤S4：
S4：手动调整局部特征点，包括删除、增加局部特征点。
根据本发明的实施方案，当所述配准方法的步骤S3自动进行时，所述配准方法还可以包括如下步骤S5：
S5：手动调整局部特征点对，包括删除、增加局部特征点对。
根据本发明的实施方案，步骤S1可以适用于多种不同的神经环路成像方式。例如，神经环路的成像方式可以选自台式双光子成像、小型头戴式双光子成像、小型头戴式单光子成像、三光子成像等。
优选地，步骤S1还包括将神经元活动转化为光信号。其中，将神经元活动转化为光信号的方式可以选自OGB、Ga520、病毒、转基因等。
优选地，步骤S1针对的对象可以为动物体，例如选自小鼠、大鼠、猴、树鼩等。作为实例，可以采用多时程的方式对同一动物体的特定神经环路进行多次观察，包括一天内不同时间段的多次观察，和/或多天的观察。其中，可以将每一个时程的神经环路成像结果得到一个图像序列，对每一个时程的成像图像序列提取一个特征图像。
根据本发明的实施方案，从图像序列中提取特征图像的方式可以选自下列中的一种或多种：Z轴平均，Z轴最大投射，Z轴方差等。
根据本发明的实施方案，对两个时程的特征图像，可以将其中一张特征图像作为静态图像IM，另外一张图像作为动态图像im。在图像配准过程中，静态图像保持不变，动态图像向静态图像做非刚性变换，实现图像配准。
根据本发明的实施方案，在步骤S2从特征图像中识别局部特征点的过程中，关键的局部图像特征包括但不限于选自神经元、轴突、树突和血管等的局部图像特征。
例如，可以采用多种方式从特征图像中识别局部特征点，包括选自下列中的一种或多种：BRISK特征、MSER特征、ORB特征、SURF特征、KAZE特征、基于FAST算法的cornerPoints特征、基于Harris–Stephens算法的cornerPoints特征、基于minimumeigenvalue算法的cornerPoints特征等。
根据本发明的实施方案，较好的局部图像特征应具有旋转不变性、尺度不变性和对噪声的鲁棒性等。
根据本发明的实施方案，在步骤S3的局部特征点配对过程中，尤其是自动配对的过程中，可以采用多尺度的二维傅里叶变换(2D FFT)的方式实现特征点配对。其中，静态图像中的特征点表示为Pk＝(Xk,Yk)，动态图像中潜在的与之对应的特征点表示为pk。优选地，Pk与pk处于不同图像中，但是表示同一个对象，原因在于图像变形可能导致二者的空间位置存在差异。
根据本发明的实施方案，步骤S3可包括步骤S3-1、S3-2、S3-3和S3-4：
S3-1：设定一级二维傅里叶变换特征点配对参数；
S3-2：提取图像块对，计算互相关矩阵，筛选潜在的特征点对；
S3-3：设定二级二维傅里叶变换特征点配对参数；
S3-4：提取图像块对，计算互相关矩阵，确定最终的特征点对。
根据本发明的实施方案，在步骤S3-1设定一级二维傅里叶变换特征点配对参数的过程中，可根据需要设置图像块大小的取值S、以及潜在的配对特征点互相关阈值CorrThre1。
根据本发明的实施方案，在步骤S3-2提取图像块对，计算互相关矩阵，筛选潜在的特征点对的过程中，分别以Pk与pk为中心点，从静态和动态图像提取图像块IMk1和imk1，图像块的大小为S×S：
IMk1(x,y)＝IM(Xk+x-S/2,Yk+y-S/2)
imk1(x,y)＝im(Xk+x-S/2,Yk+y-S/2)
其中x,y＝1,2,…,S；
计算图像块IMk1和imk1之间的互相关矩阵R1，计算互相关矩阵R1最大峰值与次大峰值之间的差值DiffTMax1；
比较差值DiffTMax1与配对特征点互相关阈值CorrThre1之间的关系：
如果Flag1＝0，特征点Pk与pk被判断为不配对的特征点；
如果Flag1＝1，特征点Pk与pk被判断为潜在的配对特征点，然后计算互相关矩阵R1最大峰值距离0点的偏移(Δx1,Δy1)，根据偏移量(Δx1,Δy1)调整动态图像中特征点pk的位置；
pk＝(Xk-Δx1,Yk-Δy1)。
根据本发明的实施方案，在步骤S3-3设定二级二维傅里叶变换特征点配对参数的过程中，需要设置图像块大小的取值T、以及潜在的配对特征点互相关阈值CorrThre2；
作为实例，步骤S3-1的图像块尺寸S大于步骤S3-3的图像块尺寸T，作为优选，T/S的取值范围为[0.4 0.8]。
根据本发明的实施方案，在步骤S3-4提取图像块对，计算互相关矩阵，确定最终的特征点对的过程中，开展精细尺度的2D FFT；
分别以Pk与更新后的pk为中心点，从静态和动态图像提取图像块IMk2和imk2，图像块的大小为T×T：
IMk2(x,y)＝IM(Xk+x-T/2,Yk+y-T/2)
imk2(x,y)＝im(Xk+x-Δx1-T/2,Yk+y-Δy1-T/2)
其中x,y＝1,2,…,T；
计算图像块IMk2和imk2之间的互相关矩阵R2，计算互相关矩阵R2最大峰值与次大峰值之间的差值DiffTMax2；
比较差值DiffTMax2与配对特征点互相关阈值CorrThre2之间的关系：
如果Flag2＝0，特征点Pk与pk被判断为不配对的特征点；
如果Flag2＝1，特征点Pk与pk被确定为配对特征点，然后计算互相关矩阵R2最大峰值距离0点的偏移(Δx2,Δy2)，根据偏移量(Δx1,Δy1)与(Δx2,Δy2)确定动态图像中特征点pk的最终位置，记为
xk＝Xk-Δx1-Δx2
yk＝Yk-Δy1-Δy2
分别位于静态图像和动态图像中的特征点Pk与为最终的配对特征点对。
根据本发明的实施方案，在步骤S4手动调整局部特征点，包括删除、增加局部特征点的过程中，可以通过人机交互的方式，手动添加和删除局部特征点；
根据本发明的实施方案，在步骤S5手动调整局部特征点对，包括删除、增加局部特征点对的过程中，可以通过人机交互的方式，手动为增加的局部特征点进行配对，以及删除已经存在的局部特征点对；
根据本发明的实施方案，在步骤S6向外扩展图像，增加扩展图像边界特征点的过程中，以原始的静态图像和动态图像为中心，将图像向外扩展，扩展的参数为z，原始图像的尺寸为M×N，扩展后的图像尺寸为(M*z)×(N*z)；作为优选，扩展参数z的取值范围为[1.11.5]。扩展得到的图像部分可以用多种方式进行填充，如全0填充、全1填充、任意给定值填充等。
优选地，在扩展后的图像边界上，按照等距离的方式设定特征点；
优选地，扩展后的静态图像和动态图像的边界上的特征点，具有相同的空间位置。
根据本发明的实施方案，在步骤S7根据配对的特征点，将特征图像进行三角剖分的过程中，采用基于Delaunay三角剖分方法，将静态图像和动态图像进行三角剖分。
优选地，在对静态图像和动态图像进行三角剖分的过程中，二者具有相同的三角形网格对应顺序；即对于静态图像中三角剖分之后的任意一个三角形图像块，在动态图像三角剖分的结果中存在另一个三角形图像块，这两个三角形图像块对应三组配对的特征点对；静态图像IM与动态图像im被分解为N个三角形图像块：
其中imPatchn和imPatchn分别表示静态和动态图像中的第n个三角形图像块；
根据本发明的实施方案，在步骤S8对每一个三角形图像块进行分段仿射变换的过程中，对所有的IMPatchn和imPatchn三角形图像块对进行配准。具体操作为imPatchn向IMPatchn配准，配准后的imPatchn记为imPatchRegn。
imPatchn向IMPatchn配准的变换矩阵记为TF：
对每个配对的特征点对，Pk＝(Xk,Yk)，
[Xk Yk 1]＝[xk yk 1]·TF
优选地，每一组三角形图像块对imPatchn向IMPatchn具有3组配对的特征点，其中imPatchn的特征点记为(x1n,y1n),(x2n,y2n),(x3n,y3n)，IMPatchn的特征点记为(X1n,Y1n),(X2n,Y2n),(X3n,Y3n)：
优选地，求解上述方程组得到变换矩阵记TF，并根据仿射变换得到配准后的图像块imPatchRegn。
imPatchRegn＝warp(imPatchn,TF)
根据本发明的实施方案，在步骤S9合并所有分段仿射变换后的三角形图像块，完成非刚性图像配配准的过程中，将所有配准后的三角形图像块imPatchRegn组合起来，得到配准的动态图像imReg：
根据本发明的实施方案，还提供一种用于实现所述配准方法的非刚性图像配准系统，包括适用于所述配准方法的图像输入装置、图像提取和特征点识别装置、图像剖分装置、剖分图像处理装置和结果输出装置。
本发明还提供所述配准方法或所述系统的用途，其用于跨多个时程的神经环路成像配准。
有益效果
本发明的技术方案是一种基于三角剖分和分段仿射变换的非刚性图像配准方法及系统，采用多尺度的二维傅里叶变换(2D FFT)的方式实现特征点配对，根据配对的特征点，将静态图像与动态图像进行三角剖分，以及对三角剖分后的三角形图像块，开展分段仿射变换，实现图像配准。并且，本发明将分段配准的三角形图像块组合起来，高效和高准确度地实现了跨多个时程的神经环路成像配准。
附图说明
图1为本发明基于三角剖分和分段仿射变换的非刚性图像配准方法的步骤示意图；
图2为本发明非刚性图像配准方法步骤3-1至3-4的示意图；
图3为实施例1中对存在空间位置变形的多时程神经环路成像结果进行图像配准的结果；
图4为实施例2中从静态图像和动态图像中识别局部特征点的结果；
图5为实施例3中基于二维傅里叶变换的特征点配对识别效果；
图6为实施例4中将三角剖分后的三角形图像块进行分段仿射变换，并合并所有分段仿射变换后的三角形图像块，完成非刚性图像配配准的结果；
图7为实施例5中将三角剖分后的三角形图像块进行分段仿射变换，并合并所有分段仿射变换后的三角形图像块，完成非刚性图像配配准的结果；
图8为实施例6中图像配准前后的效果。
具体实施方式
下文将结合具体实施例对本发明的技术方案做更进一步的详细说明。应当理解，下列实施例仅为示例性地说明和解释本发明，而不应被解释为对本发明保护范围的限制。凡基于本发明上述内容所实现的技术均涵盖在本发明旨在保护的范围内。
除非另有说明，以下实施例中使用的仪器均为市售商品。
实施例1
本实施例提供利用本发明所述的一种基于三角剖分和分段仿射变换的非刚性图像配准方法及系统，对存在空间位置变形的多时程神经环路成像结果进行图像配准。
本实施例中，使用双光子(two-photon,2p)神经环路成像技术，记录小鼠PFC皮层的神经环路活动。所有动物协议均经北京大学机构动物护理和使用委员会(IACUC)批准。实验程序是根据实验动物护理和使用指南进行的。所有小鼠均饲养在北京大学实验动物中心的动物设施中，25±2℃，12小时光照/黑暗循环，直至用于2p成像实验。使用双光子显微镜(B-Scope，Thor Labs)和飞秒激光(MaiTai BB DS-OL，Spectra-Physics)监测小鼠中前额叶皮层(mPFC)中皮层神经元的荧光瞬变。激光波长设置为920nm。使用16倍水浸物镜(0.8NA，尼康)，图像视场为517.77×517.77μm，像素为768×768。扫描频率为20.4Hz。AAV9.syn.GCaMP6f.WPRE.SV40(Penn Vector Core)用于标记皮层神经元的钙信号。在头部固定手术后7天开始2p钙成像实验，这段时间记为第0天；2p成像记录的时程1-5的实验在第0、3、6、9和12天实施。每个时程记录1300帧图像序列。
如图3所示，其中静态图像(Session 1)是第一个时程神经环路图像序列的特征图像，动态图像(Session N)是第N个时程神经环路图像序列的特征图像，融合图像是将静态图像和动态图像融合在一起的结果，从融合图像可以看出来，参考图像和目标图像之间存在较大的错位。
实施例2
本实施例提供利用本发明所述的一种基于三角剖分和分段仿射变换的非刚性图像配准方法及系统，从静态图像和动态图像中识别局部特征点的结果。
如图4所示，该实施例中，特征点识别方式采用最大稳定极值区域(MSER-Maximally Stable Extremal Regions，MSER)特征识别。MSER的基本原理是对一幅灰度图像(灰度值为0～255)取阈值进行二值化处理，阈值从0到255依次递增。阈值的递增类似于分水岭算法中的水面的上升，随着水面的上升，有一些较矮的丘陵会被淹没，如果从天空往下看，则大地分为陆地和水域两个部分，这类似于二值图像。在得到的所有二值图像中，图像中的某些连通区域变化很小，甚至没有变化，则该区域就被称为最大稳定极值区域。这类似于当水面持续上升的时候，有些被水淹没的地方的面积没有变化。MSER具有以下特点：对图像灰度具有仿射变换的不变性；稳定性，具有相同阈值范围内所支持的区域才会被选择；无需任何平滑处理就可以实现多尺度检测，即小的和大的结构都可以被检测到。
该实施例中，基于MATLAB R2021a软件平台实现MSER特征识别，采用detectMSERFeatures函数实现特征识别，参数设置为RegionAreaRange：[40 100]；MaxAreaVariation：0.1；ThresholdDelta：1。对检测到的原始MSER特征做进一步筛选，得到最终的MSER特征；筛选方式为：计算所有检测到的MSER特征两两之间的交并比(Intersection over Union，IoU)，对于交并比大于阈值(0.2)的两个MSER特征，保留面积较大的MSER特征，删除面积较小的MSER特征。
实施例3
本实施例提供利用本发明所述的一种基于三角剖分和分段仿射变换的非刚性图像配准方法及系统，实现基于二维傅里叶变换的特征点配对识别效果。
如图5所示，在实施例中，一级二维傅里叶变换特征点配对参数为：图像块大小的取值S＝36；潜在的配对特征点互相关阈值CorrThre1＝0.5。对于图示的3个MSER特征点，在静态图像中提取图像块P1,P2,P3,以及在动态图像中对应的图像块p1,p2,p3。基于MATLABR2021a软件平台的fft2函数实现二维快速傅里叶变换，计算图像块对之间的互相关矩阵，并计算互相关矩阵中最大峰值与次大峰值之间的差值DiffTMax。对于图像I1与图像I2，计算其互相关矩阵的方式为：
FI1＝fft2(I1)；
FI2＝fft2(I2)；
FR＝FI1.*conj(FI2)；％calculating correlation
R＝ifft2(FR)；
R＝fftshift(R)；
图像块对P1与p1，P2与p2，P3与p3之间的DiffTMax值分别为3.45，1.30，0.39。根据设定的互相关阈值CorrThre1(0.5)，因此P1与p1，P2与p2被判定为潜在的特征点对，P3与p3被判定为非配对的特征点。
后续进一步设定二级二维傅里叶变换特征点配对参数为：图像块大小的取值T＝24；潜在的配对特征点互相关阈值CorrThre2＝0.5，对于潜在的特征点对P1与p1，P2与p2，进行精细的特征点配对筛查。
实施例4
本实施例提供利用本发明所述的一种基于三角剖分和分段仿射变换的非刚性图像配准方法及系统，根据配对的特征点对将图像进行三角剖分的结果。
如图6所示，在该实施例中，其中参考图像和目标图像中的原始的配对特征点分别用圆形和加号标记。原始的特征点一共有8对；将原始的参考图像和目标图像进行扩展，扩展参数为z＝1.2，并对扩展区域填充0；在扩展的图像边界上建立8对特征点，包括4对角点和4对图像边界中点；使用基于MATLAB R2021a软件平台实现Delaunay三角剖分，具体实现函数为delaunay；一共得到22个三角形图像块。
实施例5
本实施例提供利用本发明所述的一种基于三角剖分和分段仿射变换的非刚性图像配准方法及系统，将三角剖分后的三角形图像块进行分段仿射变换，并合并所有分段仿射变换后的三角形图像块，完成非刚性图像配配准的结果，该结果如图7所示。
在该实施例中，针对图6所示的所有三角剖分后得到的三角形图像块对，对每一对三角形图像块对进行仿射变换。仿射变换的结果是：将动态图像中的每一个三角形图像块，向对应的静态图像三角形图像块的位置进行变换；即将动态图像的内容，映射到静态图像对应的位置，从而实现每一个三角形图像块的配准。使用基于MATLAB R2021a软件平台实现每对三角形图像块的仿射变换，具体的实现函数为estimateGeometricTransform，输入为三角形图像块的3对位置坐标，方法为affine；计算得到变换矩阵TF.
实施例6
本实施例提供利用本发明所述的一种基于三角剖分和分段仿射变换的非刚性图像配准方法及系统，图像配准前后的效果。
如图8所示，配准前，静态图像与动态图像的融合图像中，图像内容之间存在较大的空间位置错位。配准后，静态图像与配准后的融合图像中，图像内容之间的不匹配得到消除。
以上对本发明技术方案的实施方式进行了示例性的说明。应当理解，本发明的保护范围不拘囿于上述实施方式。凡在本发明的精神和原则之内，本领域技术人员所做的任何修改、等同替换、改进等，均应包含在本申请权利要求书的保护范围之内。

Claims (24)

1.一种非刚性图像配准方法，其中所述方法包括：将特征图像进行剖分，得到三角形图像块，对三角形图像块进行分段仿射变换，以及合并所有分段仿射变换后的三角形图像块；其中所述非刚性图像配准方法包括以下步骤：
S1：从多时程神经环路成像结果中提取特征图像；
S2：从特征图像中识别局部特征点；
S3：将局部特征点进行配对；
S6：向外扩展图像，增加扩展图像边界特征点；
S7：根据配对的特征点，将特征图像进行三角剖分，得到三角形图像块；
S8：对每一个三角形图像块进行分段仿射变换；
S9：合并所有分段仿射变换后的三角形图像块，得到最终的非刚性图像配对结果；步骤S1针对的对象为动物体；采用多时程的方式对同一动物体的特定神经环路进行多次观察，包括一天内不同时间段的多次观察，和/或多天的观察；将每一个时程的神经环路成像结果得到一个图像序列，对每一个时程的成像图像序列提取一个特征图像；
对两个时程的特征图像，将其中一张特征图像作为静态图像IM，另外一张图像作为动态图像im；在图像配准过程中，静态图像保持不变，动态图像向静态图像做非刚性变换，实现图像配准。
2.如权利要求1所述的配准方法，其中所述配准方法的步骤S2和/或S3自动进行或非自动进行；
当所述配准方法的步骤S2自动进行时，所述配准方法还包括如下步骤S4：
S4：手动调整局部特征点，包括删除、增加局部特征点；
和/或
当所述配准方法的步骤S3自动进行时，所述配准方法还包括如下步骤S5：
S5：手动调整局部特征点对，包括删除、增加局部特征点对。
3.如权利要求1所述的配准方法，其中步骤S1适用于选自下列的神经环路成像方式：台式双光子成像、小型头戴式双光子成像、小型头戴式单光子成像、三光子成像。
4.如权利要求1所述的配准方法，步骤S1还包括将神经元活动转化为光信号，将神经元活动转化为光信号的方式选自OGB、病毒、转基因。
5.如权利要求1所述的配准方法，从图像序列中提取特征图像的方式选自下列中的一种或多种：Z轴平均，Z轴最大投射，Z轴方差。
6.如权利要求1-5任一项所述的配准方法，其中在步骤S2从特征图像中识别局部特征点的过程中，局部图像特征包括但不限于选自神经元、轴突、树突和血管的局部图像特征；
从特征图像中识别局部特征点选自下列中的一种或多种：BRISK特征、MSER特征、ORB特征、SURF特征、KAZE特征、基于FAST算法的cornerPoints特征、基于Harris–Stephens算法的cornerPoints特征、基于minimum eigenvalue算法的cornerPoints特征。
7.如权利要求1-5任一项所述的配准方法，其中在步骤S3的局部特征点配对过程中，采用多尺度的二维傅里叶变换(2D FFT)的方式实现特征点配对；其中，静态图像中的特征点表示为Pk＝(Xk，Yk)，动态图像中潜在的与之对应的特征点表示为pk。
8.如权利要求7所述的配准方法，Pk与pk处于不同图像中，但是表示同一个对象。
9.如权利要求7所述的配准方法，其中，步骤S3包括如下步骤S3-1、S3-2、S3-3和S3-4：
S3-1：设定一级二维傅里叶变换特征点配对参数；
S3-2：提取图像块对，计算互相关矩阵，筛选潜在的特征点对；
S3-3：设定二级二维傅里叶变换特征点配对参数；
S3-4：提取图像块对，计算互相关矩阵，确定最终的特征点对。
10.如权利要求9所述的配准方法，在步骤S3-1设定一级二维傅里叶变换特征点配对参数的过程中，可根据需要设置图像块大小的取值S、以及潜在的配对特征点互相关阈值CorrThre1。
11.如权利要求9所述的配准方法，在步骤S3-2提取图像块对，计算互相关矩阵，筛选潜在的特征点对的过程中，分别以Pk与pk为中心点，从静态和动态图像提取图像块IMk1和imk1，图像块的大小为S×S：
IMk1(x,y)＝IM(Xk+x-S/2,Yk+y-S/2)
imk1(x,y)＝im(Xk+x-S/2,Yk+y-S/2)
其中x,y＝1,2,…,S；Xk、Yk分别为中心点PK的坐标点(Xk、Yk)的横坐标点值和纵坐标点值；
计算图像块IMk1和imk1之间的互相关矩阵R1，计算互相关矩阵R1最大峰值与次大峰值之间的差值DiffTMax1；
比较差值DiffTMax1与配对特征点互相关阈值CorrThre1之间的关系：
如果Flag1＝0，特征点Pk与pk被判断为不配对的特征点；
如果Flag1＝1，特征点Pk与pk被判断为潜在的配对特征点，然后计算互相关矩阵R1最大峰值距离0点的偏移(Δx1，Δy1)，根据偏移量(Δx1，Δy1)调整动态图像中特征点pk的位置；
pk＝(Xk-Δx1，Yk-Δy1)。
12.如权利要求9所述的配准方法，其中，在步骤S3-3设定二级二维傅里叶变换特征点配对参数的过程中，需要设置图像块大小的取值T、以及潜在的配对特征点互相关阈值CorrThre2。
13.如权利要求12所述的配准方法，步骤S3-1的图像块尺寸S大于步骤S3-3的图像块尺寸T，T/S的取值范围为[0.40.8]。
14.如权利要求11所述的配准方法，在步骤S3-4提取图像块对，计算互相关矩阵，确定最终的特征点对的过程中，开展精细尺度的2DFFT；
分别以Pk与更新后的pk为中心点，从静态和动态图像提取图像块IMk2和imk2，图像块的大小为T×T：
IMk2(x，y)＝IM(Xk+x-T/2，Yk+y-T/2)
imk2(x，y)＝im(Xk+x-Δx1-T/2，Yk+y-Δy1-T/2)
其中x，y＝1，2，...，T；
计算图像块IMk2和imk2之间的互相关矩阵R2，计算互相关矩阵R2最大峰值与次大峰值之间的差值DiffTMax2；
比较差值DiffTMax2与配对特征点互相关阈值CorrThre2之间的关系：
如果Flag2＝0，特征点Pk与pk被判断为不配对的特征点；
如果Flag2＝1，特征点Pk与pk被确定为配对特征点，然后计算互相关矩阵R2最大峰值距离0点的偏移(Δx2，Δy2)，根据偏移量(Δx1，Δy1)与(Δx2，Δy2)确定动态图像中特征点pk的最终位置，记为
xk＝Xk-Δx1-Δx2
yk＝Yk-Δy1-Δy2
分别位于静态图像和动态图像中的特征点Pk与为最终的配对特征点对。
15.如权利要求1-5任一项所述的配准方法，在步骤S4手动调整局部特征点，包括删除、增加局部特征点的过程中，通过人机交互的方式，手动添加和删除局部特征点。
16.如权利要求1-5任一项所述的配准方法，在步骤S5手动调整局部特征点对，包括删除、增加局部特征点对的过程中，通过人机交互的方式，手动为增加的局部特征点进行配对，以及删除已经存在的局部特征点对。
17.如权利要求1-5任一项所述的配准方法，在步骤S6向外扩展图像，增加扩展图像边界特征点的过程中，以原始的静态图像和动态图像为中心，将图像向外扩展，扩展的参数为z，原始图像的尺寸为M×N，扩展后的图像尺寸为(M*z)×(N*z)。
18.如权利要求17所述的配准方法，扩展参数z的取值范围为[1.11.5]。
19.如权利要求17所述的配准方法，在扩展后的图像边界上，按照等距离的方式设定特征点；扩展后的静态图像和动态图像的边界上的特征点，具有相同的空间位置。
20.如权利要求1-5任一项所述的配准方法，在步骤S7根据配对的特征点，将特征图像进行三角剖分的过程中，采用基于Delaunay三角剖分方法，将静态图像和动态图像进行三角剖分。
21.如权利要求14所述的配准方法，在对静态图像和动态图像进行三角剖分的过程中，二者具有相同的三角形网格对应顺序；即对于静态图像中三角剖分之后的任意一个三角形图像块，在动态图像三角剖分的结果中存在另一个三角形图像块，这两个三角形图像块对应三组配对的特征点对；静态图像IM与动态图像im被分解为N个三角形图像块：
其中IMPatchn和imPatchn分别表示静态和动态图像中的第n个三角形图像块。
22.如权利要求21所述的配准方法，在步骤S8对每一个三角形图像块进行分段仿射变换的过程中，对所有的IMPatchn和imPatchn三角形图像块对进行配准：imPatchn向IMPatchn配准，配准后的imPatchn记为imPatchRegn；
imPatchn向IMPatchn配准的变换矩阵记为TF：
对每个配对的特征点对，Pk＝(Xk,Yk)，
[Xk Yk 1]＝[xk yk 1]·TF
23.如权利要求22所述的配准方法，每一组三角形图像块对imPatchn向IMPatchn具有3组配对的特征点，其中imPatchn的特征点记为(x1n,y1n),(x2n,y2n),(x3n,y3n)，IMPatchn的特征点记为(X1n,Y1n),(X2n,Y2n),(X3n,Y3n)：
求解方程组得到变换矩阵记TF，并根据仿射变换得到配准后的图像块imPatchRegn；
imPatchRegn＝warp(imPatchn,TF)。
24.如权利要求22所述的配准方法，在步骤S9合并所有分段仿射变换后的三角形图像块，完成非刚性图像配配准的过程中，将所有配准后的三角形图像块imPatchRegn组合起来，得到配准的动态图像imReg：
CN202111642914.4A 2021-12-29 2021-12-29 非刚性图像配准方法及其系统和用途 Active CN114359357B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111642914.4A CN114359357B (zh) 2021-12-29 2021-12-29 非刚性图像配准方法及其系统和用途

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111642914.4A CN114359357B (zh) 2021-12-29 2021-12-29 非刚性图像配准方法及其系统和用途

Publications (2)

Publication Number Publication Date
CN114359357A CN114359357A (zh) 2022-04-15
CN114359357B true CN114359357B (zh) 2025-11-18

Family

ID=81103743

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111642914.4A Active CN114359357B (zh) 2021-12-29 2021-12-29 非刚性图像配准方法及其系统和用途

Country Status (1)

Country Link
CN (1) CN114359357B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105869153A (zh) * 2016-03-24 2016-08-17 西安交通大学 融合相关块信息的非刚性人脸图像配准方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2284795A2 (de) * 2001-03-08 2011-02-16 Universite Joseph Fourier Quantitative Analyse, Visualisierung und Bewegungskorrektur in dynamischen Prozessen
CN106934761A (zh) * 2017-02-15 2017-07-07 苏州大学 一种三维非刚性光学相干断层扫描图像的配准方法
CN110599478B (zh) * 2019-09-16 2023-02-03 中山大学 一种图像区域复制粘贴篡改检测方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105869153A (zh) * 2016-03-24 2016-08-17 西安交通大学 融合相关块信息的非刚性人脸图像配准方法

Also Published As

Publication number Publication date
CN114359357A (zh) 2022-04-15

Similar Documents

Publication Publication Date Title
Li et al. Online tracking of migrating and proliferating cells imaged with phase-contrast microscopy
Serradell et al. Non-rigid graph registration using active testing search
CN109767454B (zh) 基于时-空-频显著性的无人机航拍视频运动目标检测方法
Alba et al. Automatic cardiac LV segmentation in MRI using modified graph cuts with smoothness and interslice constraints
CN113608663A (zh) 一种基于深度学习和k-曲率法的指尖跟踪方法
CN115423851B (zh) 一种基于os-sift的可见光-sar图像配准算法
CN109300113A (zh) 一种基于改进凸包方法的肺结节辅助检测系统及方法
CN105574527A (zh) 一种基于局部特征学习的快速物体检测方法
Jaiswal et al. Tracking virus particles in fluorescence microscopy images using multi-scale detection and multi-frame association
CN110782428A (zh) 一种用于构建临床脑部ct图像roi模板的方法及系统
Siddiqui et al. Clustering techniques for image segmentation
Hazarika et al. A comparative study on different skull stripping techniques from brain magnetic resonance imaging
Mahapatra et al. MRF-based intensity invariant elastic registration of cardiac perfusion images using saliency information
CN119445005A (zh) 一种基于视觉的点云图像融合方法
Feng et al. Retinal mosaicking with vascular bifurcations detected on vessel mask by a convolutional network
CN109766850B (zh) 基于特征融合的指纹图像匹配方法
CN113780421B (zh) 基于人工智能的脑部pet影像识别方法
Tang et al. Retinal image registration based on robust non-rigid point matching method
Fang et al. Images crack detection technology based on improved K-means algorithm
Chen et al. Image segmentation based on mathematical morphological operator
CN114359357B (zh) 非刚性图像配准方法及其系统和用途
Rohr et al. Tracking and quantitative analysis of dynamic movements of cells and particles
Tanács et al. Estimation of linear deformations of 2D and 3D fuzzy objects
Palraj et al. Retinal fundus image registration via blood vessel extraction using binary particle swarm optimization
Mohammed et al. A Review of Image Segmentation Strategies from Classical Methods to Deep Learning

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