CN111260701B - 多模态视网膜眼底图像配准方法及装置 - Google Patents

多模态视网膜眼底图像配准方法及装置 Download PDF

Info

Publication number
CN111260701B
CN111260701B CN202010024197.8A CN202010024197A CN111260701B CN 111260701 B CN111260701 B CN 111260701B CN 202010024197 A CN202010024197 A CN 202010024197A CN 111260701 B CN111260701 B CN 111260701B
Authority
CN
China
Prior art keywords
feature
point set
characteristic
difference matrix
acquiring
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
CN202010024197.8A
Other languages
English (en)
Other versions
CN111260701A (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.)
South China University of Technology SCUT
Original Assignee
South China University of Technology SCUT
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 South China University of Technology SCUT filed Critical South China University of Technology SCUT
Priority to CN202010024197.8A priority Critical patent/CN111260701B/zh
Publication of CN111260701A publication Critical patent/CN111260701A/zh
Application granted granted Critical
Publication of CN111260701B publication Critical patent/CN111260701B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G06F18/253—Fusion techniques of extracted features
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
    • G06V10/462—Salient features, e.g. scale invariant feature transforms [SIFT]
    • G06V10/464—Salient features, e.g. scale invariant feature transforms [SIFT] using a plurality of salient features, e.g. bag-of-words [BoW] representations
    • 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/30041—Eye; Retina; Ophthalmic
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T10/00—Road transport of goods or passengers
    • Y02T10/10—Internal combustion engine [ICE] based vehicles
    • Y02T10/40—Engine management systems

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Software Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Medical Informatics (AREA)
  • Health & Medical Sciences (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • General Health & Medical Sciences (AREA)
  • Quality & Reliability (AREA)
  • Multimedia (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Computing Systems (AREA)
  • Image Analysis (AREA)
  • Eye Examination Apparatus (AREA)

Abstract

本申请公开了一种多模态视网膜眼底图像配准方法及装置，所述方法包括：根据尺度不变特征转换算法，提取浮动图像中的第一特征点集参考图像中的第二特征点集；根据第一特征点集以及第二特征点集的第二形状上下文特征，获取第一特征差异矩阵，并根据第一特征点集和第二特征点集的第二纹理特征，获取第二特征差异矩阵；根据期望最大化算法，对第一特征差异矩阵和第二特征差异矩阵通过高斯混合模型和贝叶斯定律求解进行计算，获取第一特征差异矩阵与第二特征差异矩阵基于混合高斯模型的贝叶斯规则的后验概率矩阵后，根据后验概率矩阵进行计算，获取点集坐标，直至期望最大化算法的计算结果收敛或达到预设的迭代次数；根据点集坐标，获取配准后的图像。

Description

多模态视网膜眼底图像配准方法及装置
技术领域
本申请涉及图像处理技术领域，尤其涉及一种多模态视网膜眼底图像配准方法及装置。
背景技术
眼底的视网膜血管是人体中可看见的血管，医生把它当作了解其它脏器血管情况的窗口，其变化在一定程度上反映了一些器官的改变程度。为能够更好的观测视网膜血管，通常采用图像配准的方法来对眼底彩色图像和眼底荧光造影图像所产生的图像进行处理，使两幅图像的信息互补，以得到更为全面的信息。
现有的图像配准方法，是通过对两幅图像进行Harris角点检测，得到包含分支点和交叉点的角点，判断角点的属性并记录眼底彩色图像及眼底荧光造影图像中的特征点并构成两个特征点集，根据两个特征点集进行计算得到进行图像配准的仿射变换参数后，以眼底荧光造影图像作为参考图像，以眼底彩色图像作为浮动图像，通过仿射变换参数将两者进行配准，得到配准后的图像，从而获取更为全面的信息。但在实现本申请的过程中，发现当两幅图像的局部结构相似时，仅使用单一的特征描述符去评估点集之间的对应关系，会导致大量的误判，导致最终得到的图像信息不准确。
发明内容
本申请实施例提供一种多模态视网膜眼底图像配准方法及装置，解决现有的图像配准方法进行配准时获取到的图像精度不足的问题。
为解决上述问题，本申请实施例提供一种多模态视网膜眼底图像配准方法，适于在计算设备中执行，至少包括如下步骤：
获取多模态的浮动图像和参考图像；
根据尺度不变特征转换算法，提取所述浮动图像中的多个第一特征点组成第一特征点集，以及提取所述参考图像中的多个第二特征点组成第二特征点集；
根据所述第一特征点集的第一形状上下文特征，以及所述第二特征点集的第二形状上下文特征，获取第一特征差异矩阵，并根据所述第一特征点集的第一纹理特征和所述第二特征点集的第二纹理特征，获取第二特征差异矩阵；
根据期望最大化算法，在E步处理中，对所述第一特征差异矩阵和所述第二特征差异矩阵通过高斯混合模型和贝叶斯定律求解进行计算，获取所述第一特征差异矩阵与所述第二特征差异矩阵基于混合高斯模型的贝叶斯规则的后验概率矩阵后，在M步处理中，根据所述后验概率矩阵进行计算，获取点集坐标，直至期望最大化算法的计算结果收敛或达到预设的迭代次数；
根据所述点集坐标，获取配准后的图像。
进一步的，所述根据尺度不变特征转换算法，提取所述浮动图像中的多个第一特征点组成第一特征点集，以及提取所述参考图像中的多个第二特征点组成第二特征点集，包括：
以预设的阈值，根据所述尺度不变特征转换算法对所述浮动图像和所述参考图像进行遍历，获取所述阈值下的各第一特征点和各第二特征点后，迭代地基于预设的初始值更新所述阈值，并根据每次迭代后的阈值，获取每次迭代后的阈值下的各所述第一特征点和各所述第二特征点，直至所述阈值的迭代次数达到预设次数后，生成包含所有所述第一特征点的第一特征点集，以及包含所有所述第二特征点的第二特征点集。
进一步的，所述E步处理包括：
在D维空间中，根据
Figure BDA0002360070790000021
获取第i个所述第二特征点由第j个所述第一特征点生成的生成概率；其中，yj为第j个所述第一特征点，xi为第i个第二特征点，各所述第一特征点的先验概率以及各向同性协方差σ2I相同；
为所述第一特征点和所述第一特征点添加预设权重，基于所述预设权重和所述生成概率，获取所述高斯混合模型的概率密度函数后，根据贝叶斯定理，获取第一特征点集与所述第二特征点集的欧式距离特征对应关系；
将所述第一特征差异矩阵与所述第二特征差异矩阵通过欧式距离特征对应关系进行计算，得到所述后验概率矩阵。
进一步的，所述根据所述点集坐标，获取配准后的图像，包括：
将所述点集坐标通过反推法进行图像配准，获取配准后的图像。
进一步的，还提供一种多模态视网膜眼底图像配准装置，包括：
图像获取模块，用于获取多模态的浮动图像和参考图像；
特征提取模块，用于根据尺度不变特征转换算法，提取所述浮动图像中的多个第一特征点组成第一特征点集，以及提取所述参考图像中的多个第二特征点组成第二特征点集；
差异矩阵获取模块，用于根据所述第一特征点集的第一形状上下文特征，以及所述第二特征点集的第二形状上下文特征，获取第一特征差异矩阵，并根据所述第一特征点集的第一纹理特征和所述第二特征点集的第二纹理特征，获取第二特征差异矩阵；
特征匹配模块，用于根据期望最大化算法，在E步处理中，对所述第一特征差异矩阵和所述第二特征差异矩阵通过高斯混合模型和贝叶斯定律求解进行计算，获取所述第一特征差异矩阵与所述第二特征差异矩阵基于混合高斯模型的贝叶斯规则的后验概率矩阵后，在M步处理中，根据所述后验概率矩阵进行计算，获取点集坐标，直至期望最大化算法的计算结果收敛或达到预设的迭代次数；
图像配准模块，用于根据所述点集坐标，获取配准后的图像。
进一步的，所述特征提取模块具体用于：
以预设的阈值，根据所述尺度不变特征转换算法对所述浮动图像和所述参考图像进行遍历，获取所述阈值下的各第一特征点和各第二特征点后，迭代地基于预设的初始值更新所述阈值，并根据每次迭代后的阈值，获取每次迭代后的阈值下的各所述第一特征点和各所述第二特征点，直至所述阈值的迭代次数达到预设次数后，生成包含所有所述第一特征点的第一特征点集，以及包含所有所述第二特征点的第二特征点集。
进一步的，所述E步处理包括：
在D维空间中，根据
Figure BDA0002360070790000041
获取第i个所述第二特征点由第j个所述第一特征点生成的生成概率；其中，yj为第j个所述第一特征点，xi为第i个第二特征点，各所述第一特征点的先验概率以及各向同性协方差σ2I相同；
为所述第一特征点和所述第一特征点添加预设权重，基于所述预设权重和所述生成概率，获取所述高斯混合模型的概率密度函数后，根据贝叶斯定理，获取第一特征点集与所述第二特征点集的欧式距离特征对应关系；
将所述第一特征差异矩阵与所述第二特征差异矩阵通过欧式距离特征对应关系进行计算，得到所述后验概率矩阵。
进一步的，所述图像配准模块具体用于：
将所述点集坐标通过反推法进行图像配准，获取配准后的图像。
实施本申请实施例，具有如下有益效果：
本实施例采用形状上下文特征和纹理特征进行科学计算的方式，形成全局和局部特征的混合特征框架，以使其在配准过程中形成互补关系，从而避免产生大量误判，从而提高图像配准的精准度。
附图说明
图1是多种病变类型在眼底图像上的具体表现的示例图；
图2本申请的一个实施例提供的多模态视网膜眼底图像配准方法的流程示意图；
图3本申请实施例的图像输出结果与现有技术的图像输出结果的对比结果图；
图4是本申请的一个实施例提供的多模态视网膜眼底图像配准装置的结构示意图。
具体实施方式
下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
眼底检查十分重要，许多疾病都可以从眼底上反映出来。眼底的视网膜血管是人体中可看见的血管，医生把它当作了解其它脏器血管情况的窗口。因此，它的变化在一定程度上反映了一些器官的改变程度。如高血压病人眼底可见到视网膜动脉硬化，糖尿病病人眼底可见毛细血管瘤、小的出血点和渗出物，具体如图1所示，其在一定程度上反映了全身的血管改变情况。
现有的图像配准方法，是通过对两幅图像进行Harris角点检测，得到包含分支点和交叉点的角点，判断角点的属性并记录眼底彩色图像及眼底荧光造影图像中的特征点并构成两个特征点集，根据两个特征点集进行计算得到进行图像配准的仿射变换参数后，以眼底荧光造影图像作为参考图像，以眼底彩色图像作为浮动图像，通过仿射变换参数将两者进行配准，得到配准后的图像，从而获取更为全面的信息。但当两幅图像存在较大的非刚性畸变，即浮动图像存在病变时，Harris固定阈值的特征点提取，会导致提取的内点数量不足，冗余点的存在会影响配准的精度，且当图像间的局部结构相似时，仅使用单一的特征描述符去评估点集之间的对应关系会导致大量的误判。
为解决上述问题，参见图2，是本申请的一个实施例提供的多模态视网膜眼底图像配准方法的流程示意图，各步骤具体如下：
步骤S11，获取多模态的浮动图像和参考图像。
在本实施例中，多模态的浮动图像和多模态的参考图像，包括多视角和多时间所捕获的浮动图像以及多视角和多时间所捕获的参考图像。
步骤S12，根据尺度不变特征转换算法，提取浮动图像中的多个第一特征点组成第一特征点集，以及提取参考图像中的多个第二特征点组成第二特征点集。
具体的，以预设的阈值，根据尺度不变特征转换算法对浮动图像和参考图像进行遍历，获取阈值下的各第一特征点和各第二特征点后，迭代地基于预设的初始值更新阈值，并根据每次迭代后的阈值，获取每次迭代后的阈值下的各第一特征点和各第二特征点，直至阈值的迭代次数达到预设次数后，生成包含所有第一特征点的第一特征点集，以及包含所有第二特征点的第二特征点集。
在本实施例中，首先采用一个较小的阈值分别对浮动图像和参考图像进行特征点提取，完成提取后，为阈值加上一个初始值，从而更新阈值，并基于更新后的阈值，重新通过尺度不变特征转换算法对浮动图像和参考图像进行特征点提取。其中，从浮动图像中提取到的特征点为第一特征点，即源点，从参考图像上提取到的特征点为第二特征点，即目标点。随着迭代的进行，阈值逐渐增加，提取的点的数量逐渐减少，此时内点的数量增加，冗余点逐渐减少。直至达到预先设置的最大迭代次数后，停止计算，并将所有提取到的第一特征点作为第一特征点集，即源点集，将所有提取到的第二特征点作为第二特征点集，即目标点集。
步骤S13，根据第一特征点集的第一形状上下文特征，以及第二特征点集的第二形状上下文特征，获取第一特征差异矩阵，并根据第一特征点集的第一纹理特征和第二特征点集的第二纹理特征，获取第二特征差异矩阵。
在本实施例中，形状上下文特征为图像的全局结构特征，纹理特征为图像的局部结构特征，通过第一形状上下文特征和第二形状上下文特征，计算获得两个点集之间的第一特征差异矩阵SC，通过第一纹理特征和第二纹理特征，计算获得两个点集之间的第二特征差异矩阵TF。
步骤S14，根据期望最大化算法，在E步处理中，对第一特征差异矩阵和第二特征差异矩阵通过高斯混合模型和贝叶斯定律求解进行计算，获取第一特征差异矩阵与第二特征差异矩阵基于混合高斯模型的贝叶斯规则的后验概率矩阵后，在M步处理中，根据后验概率矩阵进行计算，获取点集坐标，直至期望最大化算法的计算结果收敛或达到预设的迭代次数。
其中，E步处理包括：
在D维空间中，根据
Figure BDA0002360070790000071
获取第i个第二特征点由第j个第一特征点生成的生成概率；其中，yj为第j个第一特征点，xi为第i个第二特征点，各第一特征点的先验概率以及各向同性协方差σ2I相同；为第一特征点和第一特征点添加预设权重，基于预设权重和生成概率，获取高斯混合模型的概率密度函数后，根据贝叶斯定理，获取第一特征点集与第二特征点集的欧式距离特征对应关系；将第一特征差异矩阵与第二特征差异矩阵通过欧式距离特征对应关系进行计算，得到后验概率矩阵。
在本实施例中，通过期望最大化算法，实现特征点集匹配，其包括两个交替迭代的步骤：对应关系评估和空间变换更新。而在本实施例中，将两个点集的对应关系评估看作是概率密度估计问题，即通过高斯混合模型(GMM)和贝叶斯定律求解所得的后验概率来充当模糊对应。在D维空间中，令源点yj为第j个高斯成员的质心，目标点xi为第i个数据。则第i个数据由第j个高斯成员生成的概率为：
Figure BDA0002360070790000072
所有高斯成员均被分配以相等的各向同性协方差σ2I和先验概率。其中，I是单位矩阵。由于存在内点和冗余点，因此给它们添加权重
Figure BDA0002360070790000073
后，得到GMM的概率密度函数：
Figure BDA0002360070790000074
根据贝叶斯定理，得到后验概率：
Figure BDA0002360070790000075
这属于单一的欧式距离特征对应关系矩阵。将第一特征差异矩阵SC和第二特征差异矩阵TF代入上式，得到基于混和高斯模型的后验概率：
Figure BDA0002360070790000081
其中参数σ和β分别控制特征的强度，由M步通过最小化完整的负对数似然函数的期望找到的参数值进行迭代更新。
即E步(E-step)用于猜测计算基于混合模型的贝叶斯规则的后验概率分布的参数值，即求解后验概率矩阵CE，M-step通过最小化完整的负对数似然函数的期望来更新参数值。
为了更完美的完成M-step，本实施例在能量方程中加入全局-局部几何结构约束：
Figure BDA0002360070790000082
其中，第二项和第三项表示约束项,λ和η分别表示约束的强度控制。κ表示高斯核κ(yj,yi)＝exp(-||yi-yj||2/2ψ2)，常数ψ是高斯核的宽度。通过求导就可以得到括号里的参数，每一次更新的都带入E步中。就这样两个步骤不断迭代，直至能量方程收敛或达到设置的最大迭代次数，从而得到点集坐标T(Y,W)＝Y+κW。
步骤S15，根据点集坐标，获取配准后的图像。
具体的，将点集坐标通过反推法进行图像配准，获取配准后的图像。
在本实施例中，采用反推法(Backward approach)实现图像配准，其中映射函数采用薄板样条(Thin plate spline,TPS)，双三插值用于提高转换图像的平滑度。
获取配准后的图像后，使用交替显示的配准后的图像和参考图像棋盘格，能够很直观的看出配准结果，如图3所示，其中R表示参考图像，S表示浮动图像，ours表示本实施例的对比结果(棋盘格为参考图像与配准后的图像交替显示的结果)，PR-GLS表示现有技术的对比结果。
与现有技术相比，本实施例采用形状上下文特征和纹理特征进行科学计算的方式，形成全局和局部特征的混合特征框架，以使其在配准过程中形成互补关系，从而避免产生大量误判，从而提高图像配准的精准度。
除此之外，通过全局和局部的双约束去指导点集移动，以充分利用利用各自优势，使局部尺度下的任意点具有一定程度的运动自由，从而进一步提升算法在空间变换更新步骤的精度。
进一步的，参见图4，是本申请的一个实施例提供的多模态视网膜眼底图像配准装置的结构示意图。包括：
图像获取模块101，用于获取多模态的浮动图像和参考图像。
特征提取模块102，用于根据尺度不变特征转换算法，提取浮动图像中的多个第一特征点组成第一特征点集，以及提取参考图像中的多个第二特征点组成第二特征点集。
在本实施例中，特征提取模块102具体用于，以预设的阈值，根据尺度不变特征转换算法对浮动图像和参考图像进行遍历，获取阈值下的各第一特征点和各第二特征点后，迭代地基于预设的初始值更新阈值，并根据每次迭代后的阈值，获取每次迭代后的阈值下的各第一特征点和各第二特征点，直至阈值的迭代次数达到预设次数后，生成包含所有第一特征点的第一特征点集，以及包含所有第二特征点的第二特征点集。
差异矩阵获取模块103，用于根据第一特征点集的第一形状上下文特征，以及第二特征点集的第二形状上下文特征，获取第一特征差异矩阵，并根据第一特征点集的第一纹理特征和第二特征点集的第二纹理特征，获取第二特征差异矩阵。
特征匹配模块104，用于根据期望最大化算法，在E步处理中，对第一特征差异矩阵和第二特征差异矩阵通过高斯混合模型和贝叶斯定律求解进行计算，获取第一特征差异矩阵与第二特征差异矩阵基于混合高斯模型的贝叶斯规则的后验概率矩阵后，在M步处理中，根据后验概率矩阵进行计算，获取点集坐标，直至期望最大化算法的计算结果收敛或达到预设的迭代次数。
其中，E步处理包括：
在D维空间中，根据
Figure BDA0002360070790000101
获取第i个第二特征点由第j个第一特征点生成的生成概率；其中，yj为第j个第一特征点，xi为第i个第二特征点，各第一特征点的先验概率以及各向同性协方差σ2I相同；为第一特征点和第一特征点添加预设权重，基于预设权重和生成概率，获取高斯混合模型的概率密度函数后，根据贝叶斯定理，获取第一特征点集与第二特征点集的欧式距离特征对应关系；将第一特征差异矩阵与第二特征差异矩阵通过欧式距离特征对应关系进行计算，得到后验概率矩阵。
在本实施例中，通过期望最大化算法，实现特征点集匹配，其包括两个交替迭代的步骤：对应关系评估和空间变换更新。而在本实施例中，将两个点集的对应关系评估看作是概率密度估计问题，即通过高斯混合模型(GMM)和贝叶斯定律求解所得的后验概率来充当模糊对应。在D维空间中，令源点yj为第j个高斯成员的质心，目标点xi为第i个数据。则第i个数据由第j个高斯成员生成的概率为：
Figure BDA0002360070790000102
所有高斯成员均被分配以相等的各向同性协方差σ2I和先验概率。其中，I是单位矩阵。由于存在内点和冗余点，因此给它们添加权重
Figure BDA0002360070790000103
后，得到GMM的概率密度函数：
Figure BDA0002360070790000104
根据贝叶斯定理，得到后验概率：
Figure BDA0002360070790000105
这属于单一的欧式距离特征对应关系矩阵。将第一特征差异矩阵SC和第二特征差异矩阵TF代入上式，得到基于混和高斯模型的后验概率：
Figure BDA0002360070790000111
其中参数σ和β分别控制特征的强度，由M部通过最小化完整的负对数似然函数的期望找到的参数值进行迭代更新。
即E部(E-step)用于猜测计算基于混合模型的贝叶斯规则的后验概率分布的参数值，即求解后验概率矩阵CE，M-step通过最小化完整的负对数似然函数的期望来更新参数值。
为了更完美的完成M-step，本实施例在能量方程中加入全局-局部几何结构约束：
Figure BDA0002360070790000112
其中，第二项和第三项表示约束项,λ和η分别表示约束的强度控制。κ表示高斯核，κ(yj,yi)＝exp(-||yi-yj||2/2ψ2)，常数ψ是高斯核的宽度。通过求导就可以得到括号里的参数，每一次更新的都带入E步中。就这样两个步骤不断迭代，直至能量方程收敛或达到设置的最大迭代次数，从而得到点集坐标T(Y,W)＝Y+κW。
图像配准模块105，用于根据点集坐标，获取配准后的图像。
在本实施例中，图像配准模块105具体用于，将所述点集坐标通过反推法进行图像配准，获取配准后的图像。
与现有技术相比，本实施例采用形状上下文特征和纹理特征进行科学计算的方式，形成全局和局部特征的混合特征框架，以使其在配准过程中形成互补关系，从而避免产生大量误判，从而提高图像配准的精准度。
本申请的又一实施例还提供了一种多模态视网膜眼底图像配准终端设备，包括处理器、存储器以及存储在所述存储器中且被配置为由所述处理器执行的计算机程序，所述处理器执行所述计算机程序时实现如上述实施例所述的多模态视网膜眼底图像配准方法。
以上所述是本申请的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本申请原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本申请的保护范围。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random AccessMemory，RAM)等。

Claims (6)

1.一种多模态视网膜眼底图像配准方法，其特征在于，包括：
获取多模态的浮动图像和参考图像；
根据尺度不变特征转换算法，提取所述浮动图像中的多个第一特征点组成第一特征点集，以及提取所述参考图像中的多个第二特征点组成第二特征点集；
根据所述第一特征点集的第一形状上下文特征，以及所述第二特征点集的第二形状上下文特征，获取第一特征差异矩阵，并根据所述第一特征点集的第一纹理特征和所述第二特征点集的第二纹理特征，获取第二特征差异矩阵；
根据期望最大化算法，在E步处理中，对所述第一特征差异矩阵和所述第二特征差异矩阵通过高斯混合模型和贝叶斯定律求解进行计算，获取所述第一特征差异矩阵与所述第二特征差异矩阵基于混合高斯模型的贝叶斯规则的后验概率矩阵后，在M步处理中，根据所述后验概率矩阵进行计算，获取点集坐标，直至期望最大化算法的计算结果收敛或达到预设的迭代次数；
根据所述点集坐标，获取配准后的图像；
其中，所述E步处理包括：
在D维空间中，根据
Figure FDA0004246057640000011
获取第i个所述第二特征点由第j个所述第一特征点生成的生成概率；其中，yj为第j个所述第一特征点，xi为第i个第二特征点，各所述第一特征点的先验概率以及各向同性协方差σ2I相同；
为所述第一特征点和所述第一特征点添加预设权重，基于所述预设权重和所述生成概率，获取所述高斯混合模型的概率密度函数后，根据贝叶斯定理，获取第一特征点集与所述第二特征点集的欧式距离特征对应关系；
将所述第一特征差异矩阵与所述第二特征差异矩阵通过欧式距离特征对应关系进行计算，得到所述后验概率矩阵；
所述M步处理包括：
在能量方程中加入全局-局部几何结构约束：
Figure FDA0004246057640000021
其中，第二项和第三项表示约束项，λ和η分别表示约束的强度控制；κ表示高斯核κ(yj,yi)＝exp(-||yi-yj||2/2ψ2)，常数ψ是高斯核的宽度，W、σ和β分别需要优化求解的参数，CEji为后验概率矩阵，i,j示点集索引，X表示源点集，Y表示目标点集，L(T)为局部约束，G(T)为全局约束；
通过求导得到括号里的参数，每一次更新的都带入E步中；两个步骤不断迭代，直至能量方程收敛或达到设置的最大迭代次数，从而得到点集坐标T(Y,W)＝Y+κW。
2.根据权利要求1所述的多模态视网膜眼底图像配准方法，其特征在于，所述根据尺度不变特征转换算法，提取所述浮动图像中的多个第一特征点组成第一特征点集，以及提取所述参考图像中的多个第二特征点组成第二特征点集，包括：
以预设的阈值，根据所述尺度不变特征转换算法对所述浮动图像和所述参考图像进行遍历，获取所述阈值下的各第一特征点和各第二特征点后，迭代地基于预设的初始值更新所述阈值，并根据每次迭代后的阈值，获取每次迭代后的阈值下的各所述第一特征点和各所述第二特征点，直至所述阈值的迭代次数达到预设次数后，生成包含所有所述第一特征点的第一特征点集，以及包含所有所述第二特征点的第二特征点集。
3.根据权利要求1所述的多模态视网膜眼底图像配准方法，其特征在于，所述根据所述点集坐标，获取配准后的图像，包括：
将所述点集坐标通过反推法进行图像配准，获取配准后的图像。
4.一种多模态视网膜眼底图像配准装置，其特征在于，包括：
图像获取模块，用于获取多模态的浮动图像和参考图像；
特征提取模块，用于根据尺度不变特征转换算法，提取所述浮动图像中的多个第一特征点组成第一特征点集，以及提取所述参考图像中的多个第二特征点组成第二特征点集；
差异矩阵获取模块，用于根据所述第一特征点集的第一形状上下文特征，以及所述第二特征点集的第二形状上下文特征，获取第一特征差异矩阵，并根据所述第一特征点集的第一纹理特征和所述第二特征点集的第二纹理特征，获取第二特征差异矩阵；
特征匹配模块，用于根据期望最大化算法，在E步处理中，对所述第一特征差异矩阵和所述第二特征差异矩阵通过高斯混合模型和贝叶斯定律求解进行计算，获取所述第一特征差异矩阵与所述第二特征差异矩阵基于混合高斯模型的贝叶斯规则的后验概率矩阵后，在M步处理中，根据所述后验概率矩阵进行计算，获取点集坐标，直至期望最大化算法的计算结果收敛或达到预设的迭代次数；
图像配准模块，用于根据所述点集坐标，获取配准后的图像；
其中，所述E步处理包括：
在D维空间中，根据
Figure FDA0004246057640000031
获取第i个所述第二特征点由第j个所述第一特征点生成的生成概率；其中，yj为第j个所述第一特征点，xi为第i个第二特征点，各所述第一特征点的先验概率以及各向同性协方差σ2I相同；
为所述第一特征点和所述第一特征点添加预设权重，基于所述预设权重和所述生成概率，获取所述高斯混合模型的概率密度函数后，根据贝叶斯定理，获取第一特征点集与所述第二特征点集的欧式距离特征对应关系；
将所述第一特征差异矩阵与所述第二特征差异矩阵通过欧式距离特征对应关系进行计算，得到所述后验概率矩阵；
所述M步处理包括：
在能量方程中加入全局-局部几何结构约束：
Figure FDA0004246057640000041
其中，第二项和第三项表示约束项，λ和η分别表示约束的强度控制；κ表示高斯核κ(yj,yi)＝exp(-||yi-yj||2/2ψ2)，常数ψ是高斯核的宽度，W、σ和β分别需要优化求解的参数，CEji为后验概率矩阵，i,j示点集索引，X表示源点集，Y表示目标点集，L(T)为局部约束，G(T)为全局约束；
通过求导得到括号里的参数，每一次更新的都带入E步中；两个步骤不断迭代，直至能量方程收敛或达到设置的最大迭代次数，从而得到点集坐标T(Y,W)＝Y+κW。
5.根据权利要求4所述的多模态视网膜眼底图像配准装置，其特征在于，所述特征提取模块具体用于：
以预设的阈值，根据所述尺度不变特征转换算法对所述浮动图像和所述参考图像进行遍历，获取所述阈值下的各第一特征点和各第二特征点后，迭代地基于预设的初始值更新所述阈值，并根据每次迭代后的阈值，获取每次迭代后的阈值下的各所述第一特征点和各所述第二特征点，直至所述阈值的迭代次数达到预设次数后，生成包含所有所述第一特征点的第一特征点集，以及包含所有所述第二特征点的第二特征点集。
6.根据权利要求4所述的多模态视网膜眼底图像配准装置，其特征在于，所述图像配准模块具体用于：
将所述点集坐标通过反推法进行图像配准，获取配准后的图像。
CN202010024197.8A 2020-01-08 2020-01-08 多模态视网膜眼底图像配准方法及装置 Active CN111260701B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010024197.8A CN111260701B (zh) 2020-01-08 2020-01-08 多模态视网膜眼底图像配准方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010024197.8A CN111260701B (zh) 2020-01-08 2020-01-08 多模态视网膜眼底图像配准方法及装置

Publications (2)

Publication Number Publication Date
CN111260701A CN111260701A (zh) 2020-06-09
CN111260701B true CN111260701B (zh) 2023-06-30

Family

ID=70950322

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010024197.8A Active CN111260701B (zh) 2020-01-08 2020-01-08 多模态视网膜眼底图像配准方法及装置

Country Status (1)

Country Link
CN (1) CN111260701B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112308068B (zh) * 2020-10-28 2023-09-26 深圳大学 形状差异分析方法、装置、服务器及存储介质
CN112819867A (zh) * 2021-02-05 2021-05-18 苏州大学 基于关键点匹配网络的眼底图像配准方法
CN113095385B (zh) * 2021-03-31 2023-04-18 安徽工业大学 一种基于全局和局部特征描述的多模图像匹配方法
CN113298742A (zh) * 2021-05-20 2021-08-24 广东省人民医院 基于图像配准的多模态视网膜图像融合方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106485741A (zh) * 2016-10-19 2017-03-08 哈尔滨工业大学 一种保留局部结构的非刚点集配准的方法
CN109544610A (zh) * 2018-10-15 2019-03-29 天津大学 一种基于卷积神经网络的图像配准方法
CN109949349A (zh) * 2019-01-24 2019-06-28 北京大学第三医院（北京大学第三临床医学院） 一种多模态三维图像的配准及融合显示方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106485741A (zh) * 2016-10-19 2017-03-08 哈尔滨工业大学 一种保留局部结构的非刚点集配准的方法
CN109544610A (zh) * 2018-10-15 2019-03-29 天津大学 一种基于卷积神经网络的图像配准方法
CN109949349A (zh) * 2019-01-24 2019-06-28 北京大学第三医院（北京大学第三临床医学院） 一种多模态三维图像的配准及融合显示方法

Also Published As

Publication number Publication date
CN111260701A (zh) 2020-06-09

Similar Documents

Publication Publication Date Title
CN111260701B (zh) 多模态视网膜眼底图像配准方法及装置
CN111292314B (zh) 冠脉分段方法、装置、图像处理系统和存储介质
Annamalai et al. An early prediction of tumor in heart by cardiac masses classification in echocardiogram images using robust back propagation neural network classifier
Liu et al. Automatic whole heart segmentation using a two-stage u-net framework and an adaptive threshold window
JP7326499B2 (ja) 画像検出モデルのトレーニング方法及びその関連装置、機器、記憶媒体
CN111105424A (zh) 淋巴结自动勾画方法及装置
US7773789B2 (en) Probabilistic minimal path for automated esophagus segmentation
CN111383259B (zh) 图像分析方法、计算机设备和存储介质
WO2021136368A1 (zh) 钼靶图像中胸大肌区域自动检测方法及装置
US20230394670A1 (en) Anatomically-informed deep learning on contrast-enhanced cardiac mri for scar segmentation and clinical feature extraction
CN111932495B (zh) 一种医学图像检测方法、装置及存储介质
WO2020168647A1 (zh) 图像识别方法及相关设备
Tamoor et al. Automatic segmentation of medical images using a novel Harris Hawk optimization method and an active contour model
CN111127532B (zh) 基于深度学习特征光流的医学图像形变配准方法及系统
CN113327225A (zh) 用于提供气道信息的方法
Tang et al. Retinal image registration based on robust non-rigid point matching method
Furuzuki et al. A detection method for liver cancer region based on faster R-CNN
Lu et al. AutoMPR: Automatic detection of standard planes in 3D echocardiography
CN113962957A (zh) 医学图像处理方法、骨骼图像处理方法、装置、设备
CN115439423B (zh) 一种基于ct图像的识别方法、装置、设备及存储介质
CN115937163A (zh) 一种spect肺灌注成像的目标区域提取方法与系统
CN113177953B (zh) 肝脏区域分割方法、装置、电子设备及存储介质
CN112336365B (zh) 心肌血流量分布图像获取方法、系统、介质及电子设备
Paris et al. Shape constraint function for artery tracking in ultrasound images
CN111640126A (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