CN110363738A - 一种具有仿射不变性的视网膜图像配准方法及其装置 - Google Patents

一种具有仿射不变性的视网膜图像配准方法及其装置 Download PDF

Info

Publication number
CN110363738A
CN110363738A CN201810305736.8A CN201810305736A CN110363738A CN 110363738 A CN110363738 A CN 110363738A CN 201810305736 A CN201810305736 A CN 201810305736A CN 110363738 A CN110363738 A CN 110363738A
Authority
CN
China
Prior art keywords
point
image
region
candidate feature
intersection
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.)
Granted
Application number
CN201810305736.8A
Other languages
English (en)
Other versions
CN110363738B (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.)
Central South University
Original Assignee
Central South 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 Central South University filed Critical Central South University
Priority to CN201810305736.8A priority Critical patent/CN110363738B/zh
Publication of CN110363738A publication Critical patent/CN110363738A/zh
Application granted granted Critical
Publication of CN110363738B publication Critical patent/CN110363738B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/136—Segmentation; Edge detection involving thresholding
    • 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/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30041—Eye; Retina; Ophthalmic

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • Quality & Reliability (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Image Analysis (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种具有仿射不变性的视网膜图像配准方法及其装置，涉及图像处理技术领域。与传统技术整体处理相比，本发明的一种具有仿射不变性的视网膜图像配准方法，在视网膜骨骼化网络图像的基础上，通过连通域内相互连通的两个交点来截取局部图像结构(即候选特征区域)，并将局部图像结构变换(旋转、高斯滤波、压缩)到圆形高斯滤波器适用的形式下，解决了视角和尺度变化问题，提高了多幅图像的配准精度，有利于构建完整的视网膜图像信息，减少了视网膜病变给诊断带来的干扰，提高了诊断的准确性。

Description

一种具有仿射不变性的视网膜图像配准方法及其装置
技术领域
本发明属于图像处理技术领域，尤其涉及一种具有仿射不变性的视网膜图像配准方法及其装置。
背景技术
20世纪以来医学影像技术发展迅速，各种影像设备层出不穷，先后出现了X射线影像、计算机断层摄影(CT)、数字减影血管造影(DR)、彩色眼底图(CFI)、核磁共振成像(MRI)、光学相干断层扫描(OCT)和正电子发射体层成像(PET)等，在临床诊断分析中越来越倚重于医学图像的分析。不同的影像设备因成像原理不同，对同一器官的成像提供的信息也各不相同，如何将这些图像进行适当的配准得到更多的信息，成为了临床医生的迫切需求。因此，对医学图像进行配准的研究受到了众多学者的青睐，逐渐成为了热门课题。
图像配准是图像处理中极其重要的一门技术，在计算机视觉、模式识别、医学影像分析、遥感图像处理以及导弹的地形与地图匹配等方面得到广泛的应用。医学图像配准是医学影像分析的基本课题，具有重要理论研究和临床应用价值。医学图像配准是指对一幅医学图像寻求一种(或一系列)空间变换，使它与另一幅医学图像上的对应点达到空间上的一致。
研究表明，青光眼已经成为了仅次于白内障的第二大常见致盲性疾病。治疗青光眼最有效的方法是早期的筛查和诊断。彩色眼底图和视网膜OCT图像(后面对两种图像简称为视网膜图像)是眼科医生诊断青光眼的重要依据，患有青光眼一般会导致视网膜在形态上发生变化。在临床诊疗的过程中，为了获得准确、有效的视网膜图像，眼科医生往往需要对同一患者在不同时期进行多次的拍摄。由于成像角度和成像环境的不同，多次拍摄往往会造成视网膜图像出现细微的差异和变化。而这种差异变化会对病变的检测带来难以排除的干扰，从而影响诊断的准确性。另外，由于眼底相机的视角问题，每幅拍摄的内容有限，为了更容易的观测眼部是否出现病变区域，就不得不将多幅图像进行融合。因此，对视网膜眼底图像进行配准，为后续的视网膜眼底图像融合以及眼底病变检测来做准备是极其重要的。总而言之，视网膜眼底图像配准对视网膜病变诊断的起到了重要作用，所以视网膜眼底图像的配准研究对眼底疾病的早期发现、辅助诊断有着重要的意义。
发明内容
针对现有技术的不足，本发明提供一种具有仿射不变性的视网膜图像配准方法及其装置，根据高斯滤波器的形状和大小要与图像结构相适应的原理，利用视网膜图像血管树的交点来截取局部图像结构，并将局部图像结构变换到圆形高斯滤波器适用的形式下，以解决视角和尺度变化问题；为了保证图像变换的正确性，采用旋转压缩的方式将各向异性的图像结构变换为各向同性的图像结构，使得在各向同性的图像结构上提取出的尺度不变特征变换(SIFT)特征点是完全仿射不变的，使得局部特征结构不仅仅满足尺度和旋转不变，也使得局部特征结构完全仿射不变，大大提高配准的精度，保证算法在有较大仿射变换的视网膜配准中有良好的效果。
本发明是通过如下的技术方案来解决上述技术问题的：一种具有仿射不变性的视网膜图像配准方法，包括如下步骤：
步骤1：获取血管主要区域图像；
输入两幅视网膜图像，一幅作为参考图像Isrc，另一幅作为待配准图像Idst，然后分别对两幅图像依次进行提取绿色通道、高斯滤波以及多次底帽变换处理，最后将所有底帽变换所得结果叠加在一起即得到血管主要区域图像IEng；
步骤2：对血管主要区域图像IEng进行骨骼化处理，得到骨骼化网络I(x,y)，再对骨骼化网络I(x,y)进行去噪处理；
步骤3：判断各个像素点为交点还是端点；
计算骨骼化网络I(x,y)中各个像素点八连通区域内的通路数CN；
根据CN的取值将各个像素点划分为普通点POrd或者交点PCross或者端点PEnd，并形成所有交点集合和所有端点集合；
计算各个端点到与之相连通的交点的距离，如果该距离小于设定阈值，则删除该端点到与之相连通的交点所形成的分支，该分支是端点到与之相连通的交点所成的连通域，并分别在端点集合和交点集合中去除该端点和该交点；
设定阈值优选为10个像素；
图像在处理过程中大小保持始终不变，在骨骼化网络中提取的交点，在血管主要区域图像中位置不变；
步骤4：候选特征区域的截取；
依次在血管主要区域图像IEng中找到相互连通的两个交点，并截取这两个交点对应的连通域作为候选特征区域，
候选特征区域大小为RS＝LD*2*D，其中，LD表示相互连通的两个交点间的距离，D为连通域内的点到相互连通的两个交点所成直线的最大距离；
步骤5：对候选特征区域依次进行旋转、高斯滤波以及压缩处理；
旋转方向为顺时针，旋转角度为θ，其中，θ为候选特征区域内对应的相互连通的两个交点所成直线的倾斜角，如果获得最大距离的点在直线的上方，那么θ为正，如果获得最大距离的点在直线的下方，那么θ为负；
高斯滤波沿着X轴方向进行；
压缩处理是将候选特征区域的矩形区域变成正方形区域，在具体操作时，候选特征区域矩形区域的短边保持不变，将长边压缩至与短边相等；
步骤6：提取特征点，形成特征矩阵；
对压缩后的候选特征区域提取特征点，将特征点坐标变回原始视网膜图像坐标；将各个候选特征区域提取的特征点归并为特征矩阵，分别形成参考图像Isrc特征矩阵Fsrc、待配准图像Idst特征矩阵Fdst；
步骤7：对Fsrc与Fdst进行特征点匹配，并滤除误匹配的点；
步骤8：显示配准结果；
将参考图像Isrc和待配准图像Idst按照特征点匹配后得到的变换矩阵重新采样到同一坐标系，其中参考图像Isrc保持不变，待配准图像Idst作半透明处理，并使用棋盘格的方式进行配准结果显示。
进一步地，所述步骤1中，多次底帽变换采用线性结构元素，底帽变换的过程由公式表示为：
Ibot＝IGau-IGau×S
其中，Ibot表示底帽变换所得结果，IGau表示经高斯滤波后的图像，S表示线性结构元素。
进一步地，所述线性结构元素为一组具有不同尺度l和不同方向θl的线性结构元素，其中，l为整数，且取值范围为：2≤l≤15，k＝0,1,…,l-1；
针对每一个尺度，首先把该尺度所有方向上的底帽变换所得结果叠加在一起，然后通过OSTU阈值分割方法得到该尺度上的血管分割结果，最后将所有尺度上的血管分割结果叠加在一起得到血管主要区域图像IEng。
进一步地，所述步骤2中，去噪处理的具体操作为：
在八连通区域模式下，对骨骼化网络I(x,y)的连通域进行标记，然后统计各个连通域内像素点的个数，如果连通域内像素点的个数小于设定值，则去除该连通域，达到去噪的目的；
设定值的取值为100。
进一步地，所述步骤3的具体操作过程为：
设P为待检测像素点，P在八连通区域内的点记为P1,P2,L,P8，P1,P2,L,P8的排列顺序是以P1为起始点，顺时针依次排列，且P1为P正上方的点，CN的计算公式表示为：
如果则P为端点；如果则P为普通点；如果则P为交点。
根据CN的取值判断待检测像素点P是端点、交点还是普通点的依据是骨骼化网络中线的宽度为一个像素点，而且像素的取值为0和1两种，一条通路在边界处会产生差值(就像一条马路边的水渠，路面是平的，而路与水渠连接处是凹陷的)，因此，如果一个点两边是有变化的，那么这个点就是一条通路。CN表示的就是某个像素点八连通区域内的通路数，该判断方法的优势在于只需要知道八连通区域内的通路数，而不用判断具体的通路方向，计算简单，易于实现。
进一步地，所述步骤4的具体操作如下：
首先，在血管主要区域图像IEng中对所有交点进行连通域检测，如果这两个交点是相互连通的，则归为一类，并记录连通域Con，直到所有交点分类完毕，得到K个连通域，记为：Conk为第k类相互连通的交点对应的连通域；
其次，找到连通域Conk内相互连通的两个交点Pc1(x1,y1),Pc2(x2,y2)，并计算所述两个交点所成直线
再次，找到连通域Conk内到直线lc1c2距离最大的点Pmax(x3,y3)，计算最大距离其中，A、B、C是直线lc1c2的直线方程对应的系数；
最后，截取两个交点Pc1(x1,y1),Pc2(x2,y2)对应的连通域Conk作为候选特征区域，候选特征区域大小为RS＝LD*2*D。
进一步地，所述步骤5中的旋转操作是以候选特征区域左下角的交点作为中心点，旋转角度θ的正负是由Pmax(x3,y3)的位置决定，如果点Pmax(x3,y3)在直线lc1c2的上方，那么θ为正，如果点Pmax(x3,y3)在直线lc1c2的下方，那么θ为负；
高斯滤波核函数的标准差是由候选特征区域对应矩形的高与宽的比值来决定，即标准差其中，H表示候选特征区域对应矩形的高度，W表示候选特征区域对应矩形的宽度。
进一步地，所述步骤6中，特征点的提取采用SIFT(Scale Invariant FeatureTransform)算法。
进一步地，所述步骤7中，特征点的匹配采用FLANN(Fast Library forApproximate Nearest Neighbors)算法，滤除误匹配采用串行RANSAC(RANdom SAmpleConsensus)算法。
进一步地，一种具有仿射不变性的视网膜图像配准装置，包括：
图像输入模块，用于输入视网膜图像，一幅作为参考图像，另一幅作为待配准图像；
区域图像处理生成模块，用于分别对两幅图像依次进行绿色通道提取、高斯滤波以及多次底帽变换处理，并将所有底帽变换所得结果叠加在一起生成血管主要区域图像；
网络图像处理生成模块，用于对血管主要区域图像进行骨骼化处理，并生成骨骼化网络；
交点判定模块，用于根据骨骼化网络中各个像素点八连通区域内的通路数判断各个像素点是交点还是端点，并形成所有交点集合和所有端点集合；还用于计算各个端点到与之相连通的交点的距离，如果该距离小于设定阈值，则删除该端点到与之相连通的交点所形成的分支，并分别在端点集合和交点集合中去除该端点和该交点；
设定阈值优选为10个像素；
区域截取模块，用于依次在血管主要区域图像中找到相互连通的两个交点，并截取这两个交点对应的连通域作为候选特征区域，候选特征区域大小为RS＝LD*2*D，
其中，LD表示相互连通的两个交点间的距离，D为连通域内的点到相互连通的两个交点所成直线的最大距离；
区域处理模块，用于对候选特征区域依次进行顺时针旋转θ度，沿X轴方向高斯滤波，以及压缩至候选特征区域的矩形区域变成正方形区域，其中，θ为候选特征区域内对应的相互连通的两个交点所成直线的倾斜角；
特征矩阵生成模块，用于对压缩后的候选特征区域提取特征点，并将各个候选特征区域提取的特征点归分并为特征矩阵，分别形成参考图像特征矩阵和待配准图像特征矩阵；
特征点匹配模块，用于对两个特征矩阵进行特征点匹配，并滤除误匹配点；
配准结果显示模块，用于将参考图像和待配准图像按照特征点匹配后得到的变换矩阵重新采样到同一坐标系，其中参考图像保持不变，待配准图像作半透明处理，并使用棋盘格的方式进行配准结果显示；
所述的图像输入模块、区域图像处理生成模块、网络图像处理生成模块、交点判定模块、区域截取模块、区域处理模块、特征矩阵生成模块、特征点匹配模块以及配准结果显示模块的具体处理过程均是采用所述的一种具有仿射不变性的视网膜图像配准方法来实现的。
进一步地，在所述区域图像处理生成模块中，底帽变换处理单元采用线性结构元素，底帽变换处理的过程由公式表示为：
Ibot＝IGau-IGau×S
其中，Ibot表示底帽变换所得结果，IGau表示经高斯滤波后的图像，S表示线性结构元素；S为一组具有不同尺度l和不同方向θl的线性结构元素，其中，l为整数，且取值范围为：2≤l≤15，k＝0,1,…,l-1；
针对每一个尺度，首先把该尺度所有方向上的底帽变换所得结果叠加在一起，然后通过OSTU阈值分割方法得到该尺度上的血管分割结果，最后将所有尺度上的血管分割结果叠加在一起得到血管主要区域图像IEng。
进一步地，所述网络图像处理生成模块，还用于对骨骼化网络进行去噪处理；去噪处理单元在八连通区域模式下，对骨骼化网络的连通域进行标记，然后统计各个连通域内像素点的个数，如果连通域内像素点的个数小于设定值，则去除该连通域；
设定值优选为100。
进一步地，在所述交点判定模块中，设P为待检测像素点，P在八连通区域内的点记为P1,P2,L,P8，P1,P2,L,P8的排列顺序是以P1为起始点，顺时针依次排列，且P1为P正上方的点，各个像素点八连通区域内的通路数CN的计算公式表示为：
如果则P为端点；如果则P为普通点；如果则P为交点。
进一步地，所述区域截取模块，包括：
交点聚类子模块，用于在血管主要区域图像IEng中对所有交点进行连通域检测，如果这两个交点是相互连通的，则归为一类，并记录连通域Con，直到所有交点分类完毕，得到K个连通域，记为：Conk为第k类相互连通的交点对应的连通域；
直线生成子模块，用于查找连通域Conk内相互连通的两个交点Pc1(x1,y1),Pc2(x2,y2)，并计算所述两个交点所成直线
最大距离计算子模块，用于查找连通域Conk内到直线lc1c2距离最大的点Pmax(x3,y3)，计算最大距离其中，A、B、C是直线lc1c2的直线方程对应的系数；
候选特征区域截取子模块，用于截取两个交点Pc1(x1,y1),Pc2(x2,y2)对应的连通域Conk作为候选特征区域，候选特征区域大小为RS＝LD*2*D。
进一步地，在所述区域处理模块中，顺时针旋转θ度操作单元是以候选特征区域左下角的交点为中心点进行旋转操作的；
高斯滤波单元核函数的标准差是由候选特征区域对应矩形的高与宽的比值来决定，即标准差其中，H表示候选特征区域对应矩形的高度，W表示候选特征区域对应矩形的宽度。
进一步地，在所述特征矩阵生成模块中，特征点提取单元采用SIFT算法进行特征点的提取。
进一步地，在所述特征点匹配模块中，特征点匹配单元采用FLANN算法进行特征点匹配，误匹配点滤除单元采用串行RANSAC算法进行误匹配点滤除。
有益效果：
1、本发明提供了一种具有仿射不变性的视网膜图像配准方法及其装置，与传统技术整体处理相比，在视网膜骨骼化网络图像的基础上，通过连通域内相互连通的两个交点来截取局部图像结构(即候选特征区域)，并将局部图像结构变换(旋转、高斯滤波、压缩)到圆形高斯滤波器适用的形式下，解决了视角和尺度变化问题，提高了多幅图像的配准精度，有利于构建完整的视网膜图像信息，减少了视网膜病变给诊断带来的干扰，提高了诊断的准确性。
2、本发明根据CN的取值对骨骼化网络中的像素点是交点还是端点进行判断，即通过某个像素点附近的通路数判断该像素点是交点还是端点，该判断方法仅仅需要知道八连通领域内的通路数即可，不用判断具体的通路方向，计算简单，易于实现。
3、本发明对局部图像结构采用了旋转、高斯滤波以及压缩的处理方式，使各向异性的图像结构变换为各向同性的图像结构，在各向同性的图像结构上提取出的尺度不变特征变换特征点是完全仿射不变的，使得局部特征结构不仅满足了尺度和旋转不变，而且保证局部特征结构完全仿射不变，大大提高了配准的精度，保证算法在有较大仿射变换的视网膜配准中有良好的效果。
附图说明
为了更清楚地说明本发明的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一个实施例，对于本领域普通技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明的流程图；
图2A是参考图像Isrc；
图2B是待配准图像Idst；
图3是根据参考图像Isrc处理后得到的血管主要区域图像IEng；
图4是根据血管主要区域图像IEng计算得到的骨骼化网络I(x,y)；
图5是根据CN的取值划分的交点集合；
图6是根据交点截取的候选特征区域；
图7是Fsrc与Fdst进行特征点匹配的结果；
图8是配准结果。
具体实施方式
下面结合本发明实施例中的附图，对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，本发明所提供的一种具有仿射不变性的视网膜图像配准方法，包括如下步骤：
步骤1：获取血管主要区域图像；
输入两幅视网膜图像(分别对两幅图像进行步骤1-6的操作)，一幅作为参考图像Isrc(如图2A所示)，另一幅作为待配准图像Idst(如图2B所示)，然后分别对两幅图像依次进行提取绿色通道(G通道)、高斯滤波以及多次底帽变换处理，最后将所有底帽变换所得结果合并在一起即得到血管主要区域图像IEng，如图3所示。
设参考图像Isrc(x,y)＝[r(x,y)；g(x,y)；b(x,y)]，提取绿色通道后变为Isrc(x,y)＝[g(x,y)]，采用线性结构元素对经高斯滤波后的图像进行底帽变换，底帽变换的过程由公式表示为：
Ibot＝IGau-IGau×S (1)
其中，Ibot表示底帽变换所得结果，IGau表示经高斯滤波后的图像，S表示线性结构元素；
S是一组具有不同尺度l和方向θl的结构元素，其中，l为整数，其取值范围为：2≤l≤15，k＝0,1,…,l-1；
针对每一个尺度，首先把该尺度所有方向上的底帽变换所得结果叠加在一起，即然后通过OSTU阈值分割方法得到该尺度上的血管分割结果，最后将所有尺度上的血管分割结果叠加在一起得到血管主要区域图像IEng，
步骤2：对血管主要区域图像IEng进行骨骼化处理，得到骨骼化网络I(x,y)(如图4所示)，再对骨骼化网络I(x,y)进行去噪处理。
本实施例采用Zhang并行快速细化算法对血管主要区域图像IEng进行骨骼化处理得到骨骼化网络I(x,y)，具体的细化删除条件为：
设P为待检测像素点，P在八连通区域内的点记为P1,P2,L,P8，P1,P2,L,P8的排列顺序是以P1为起始点，顺时针依次排列，且P1为P正上方的点；如果满足下列四个条件，则该待检测像素点P可以删除，并反复迭代，直至所有的像素点都被遍历为止：
a、2<＝N(P)<＝6，N(P)为P在八连通区域内为提取点的数目；
b、CN(P)＝1，
c、迭代次数为奇数P1*P3*P4＝0，迭代次数为偶数P1*P3*P7＝0、P2*P4*P6＝0或者P2*P4*P8＝0；
d、迭代次数为奇数P3*P5*P7＝0，迭代次数为偶数P1*P5*P7＝0、P4*P6*P8＝0或者P2*P6*P8＝0。
去噪处理的具体操作为：
在八连通区域模式下，对骨骼化网络的连通域进行标记，然后统计各个连通域内像素点的总个数，如果连通域内像素点的个数小于设定值(本实施例中，设定值为100)，则为较小的连通域，去除该连通域。
例如，对骨骼化网络的连通域进行标记，得到了1～50号连通域，第1号连通域内包含的像素点的总个数只有25个，第2号连通域内包含的像素点的总个数有200个……，第1号连通域包含像素点个数25<100，那第1号就是较小的连通域，删除第1号连通域，第2号连通域包含像素点个数200>100，不是较小的连通域，不删除第2号连通域。
步骤3：判断各个像素点为交点还是端点；
计算骨骼化网络I(x,y)中各个像素点八连通区域内的通路数CN；
根据CN的取值将各个像素点划分为普通点POrd或者交点PCross或者端点PEnd，并形成所有交点集合和所有端点集合；
计算各个端点到与之相连通的交点的距离，如果该距离小于设定阈值(本实施例中，设定阈值为10个像素)，则删除该端点到与之相连通的交点所形成的分支，并分别在端点集合和交点集合中去除该端点和该交点。
设P为待检测像素点，P在八连通区域内的点记为P1,P2,L,P8，P1,P2,L,P8的排列顺序是以P1为起始点，顺时针依次排列，且P1为P正上方的点，CN的计算公式表示为：
如果则P为端点；如果则P为普通点；如果则P为交点，划分的交点如图5所示。
步骤4：候选特征区域的截取；
依次在血管主要区域图像IEng中找到相互连通的两个交点，并截取这两个交点对应的连通域作为候选特征区域(如图6所示)，
候选特征区域大小为RS＝LD*2*D，其中，LD表示相互连通的两个交点间的距离，D为连通域内的点到相互连通的两个交点所成直线的最大距离，两个交点所成直线与截取区域对应的宽相互平行，与高相互垂直。
步骤4的具体操作过程如下：
首先，依次在血管主要区域图像IEng中对所有交点进行连通域检测，如果这两个交点是相互连通的，则归为一类，并记录连通域Con，直到所有交点分类完毕，得到K个连通域，记为：Conk为第k类相互连通的交点对应的连通域；
其次，找到连通域Conk内相互连通的两个交点Pc1(x1,y1),Pc2(x2,y2)，并计算所述两个交点所成直线
再次，找到连通域Conk内到直线lc1c2距离最大的点Pmax(x3,y3)，计算最大距离其中，A、B、C是直线lc1c2的直线方程对应的系数；
最后，截取两个交点Pc1(x1,y1),Pc2(x2,y2)对应的连通域Conk作为候选特征区域，候选特征区域大小为RS＝LD*2*D。
步骤5：对候选特征区域依次进行旋转、高斯滤波以及压缩处理；
旋转方向为顺时针，旋转角度为θ，其中，θ为候选特征区域内对应的相互连通的两个交点所成直线的倾斜角；如果点Pmax(x3,y3)在直线lc1c2的上方，那么θ为正，如果点Pmax(x3,y3)在直线lc1c2的下方，那么θ为负；
高斯滤波沿着X轴方向进行；
压缩处理是将候选特征区域的矩形区域变成正方形区域，在具体操作时，候选特征区域矩形区域的短边保持不变，将长边压缩至与短边相等。
旋转操作是以候选特征区域左下角的交点作为中心点，(x,y)为旋转前的坐标，(x',y')为旋转后的坐标，用矩阵表示为：
并进行高斯滤波，标准差为σ的二维高斯核函数定义为：
其中，其标准差σ是由候选特征区域对应矩形的高与宽的比值来决定，即高斯函数标准差其中，H表示候选特征区域对应矩形的高度，W表示候选特征区域对应矩形的宽度；然后进行压缩，压缩后的图像变为：
则，候选特征区域的点经过旋转、高斯滤波及压缩处理后变换为：
步骤6：提取特征点，形成特征矩阵；
对压缩后的候选特征区域采用SIFT算法提取特征点，将特征点坐标变回原始视网膜图像坐标；将各个候选特征区域提取的特征点归并为特征矩阵，分别形成参考图像Isrc特征矩阵Fsrc、待配准图像Idst特征矩阵Fdst。
步骤7：匹配特征点；对Fsrc与Fdst采用FLANN算法进行特征点匹配(如图7所示)，并采用串行RANSAC算法滤除误匹配的点。
步骤8：最终结果显示；
将参考图像Isrc和待配准图像Idst按照特征点匹配后得到的变换矩阵重新采样到同一坐标系，其中参考图像Isrc保持不变，待配准图像Idst作半透明处理，并使用棋盘格的方式进行配准结果显示，如图8所示。
为了处理具有视角和尺度变化的视网膜图像，本发明方法提出一种具有仿射不变的视网膜图像配准方法，对局部特征进行提取。根据高斯滤波器的形状和大小要与图像结构相适应的原理，本发明选取视网膜图像中最稳定的血管交点作为仿射不变量，通过寻找相连通两个交点所对应的连通域，并截取该连通域对应的最小包围矩形(候选特征区域)得到局部图像结构，并将局部图像结构变换到圆形高斯滤波器适用的形式下，以解决视角和尺度变化问题；本发明对局部图像结构进一步采用旋转、高斯滤波和压缩的方式将各向异性的图像结构变换为各向同性的图像结构，保证了图像变换的正确性；最后在各向同性的图像结构上提取尺度不变特征变换(SIFT)特征点，并将SIFT特征点的坐标变回原始视网膜图像坐标完成匹配；本发明方法提取的局部特征满足高斯尺度空间理论，其对视角变化、仿射变换、噪声保持一定程度的稳定性，独特性好，信息量丰富，适用于在有视角和尺度变化的视网膜图像中进行快速、准确的配准。
一种具有仿射不变性的视网膜图像配准装置，包括：
图像输入模块，用于输入视网膜图像，一幅作为参考图像，另一幅作为待配准图像。
区域图像处理生成模块，用于分别对两幅图像依次进行绿色通道提取、高斯滤波以及多次底帽变换处理，并将所有底帽变换所得结果叠加在一起生成血管主要区域图像；
在所述区域图像处理生成模块中，底帽变换处理单元采用线性结构元素，底帽变换处理的过程由公式表示为：
Ibot＝IGau-IGau×S
其中，Ibot表示底帽变换所得结果，IGau表示经高斯滤波后的图像，S表示线性结构元素；S为一组具有不同尺度l和不同方向θl的线性结构元素，其中，l为整数，且取值范围为：2≤l≤15，k＝0,1,…,l-1；
针对每一个尺度，首先把该尺度所有方向上的底帽变换所得结果叠加在一起，然后通过OSTU阈值分割方法得到该尺度上的血管分割结果，最后将所有尺度上的血管分割结果叠加在一起得到血管主要区域图像IEng。
网络图像处理生成模块，用于对血管主要区域图像进行骨骼化处理，并生成骨骼化网络；
所述网络图像处理生成模块，还用于对骨骼化网络进行去噪处理；去噪处理单元在八连通区域模式下，对骨骼化网络的连通域进行标记，然后统计各个连通域内像素点的个数，如果连通域内像素点的个数小于设定值(设定值优选为100)，则去除该连通域。
交点判定模块，用于根据骨骼化网络中各个像素点八连通区域内的通路数判断各个像素点是交点还是端点，并形成所有交点集合和所有端点集合；还用于计算各个端点到与之相连通的交点的距离，如果该距离小于设定阈值(设定阈值为10个像素)，则删除该端点到与之相连通的交点所形成的分支，并分别在端点集合和交点集合中去除该端点和该交点；
在所述交点判定模块中，设P为待检测像素点，P在八连通区域内的点记为P1,P2,L,P8，P1,P2,L,P8的排列顺序是以P1为起始点，顺时针依次排列，且P1为P正上方的点，各个像素点八连通区域内的通路数CN的计算公式表示为：
如果则P为端点；如果则P为普通点；如果则P为交点。
区域截取模块，用于依次在血管主要区域图像IEng中找到相互连通的两个交点，并截取这两个交点对应的连通域作为候选特征区域，候选特征区域大小为RS＝LD*2*D，
其中，LD表示相互连通的两个交点间的距离，D为连通域内的点到相互连通的两个交点所成直线的最大距离。
所述区域截取模块，包括：
交点聚类子模块，用于在血管主要区域图像IEng中对所有交点进行连通域检测，如果这两个交点是相互连通的，则归为一类，并记录连通域Con，直到所有交点分类完毕，得到K个连通域，记为：Conk为第k类相互连通的交点对应的连通域；
直线生成子模块，用于查找连通域Conk内相互连通的两个交点Pc1(x1,y1),Pc2(x2,y2)，并计算所述两个交点所成直线
最大距离计算子模块，用于查找连通域Conk内到直线lc1c2距离最大的点Pmax(x3,y3)，计算最大距离其中，A、B、C是直线lc1c2的直线方程对应的系数；
候选特征区域截取子模块，用于截取两个交点Pc1(x1,y1),Pc2(x2,y2)对应的连通域Conk作为候选特征区域，候选特征区域大小为Rs＝LD*2*D。
区域处理模块，用于对候选特征区域依次进行顺时针旋转θ度，沿X轴方向高斯滤波，以及压缩至候选特征区域的矩形区域变成正方形区域，其中，θ为候选特征区域内对应的相互连通的两个交点所成直线的倾斜角；
在所述区域处理模块中，顺时针旋转θ度操作单元是以候选特征区域左下角的交点为中心点进行旋转操作的；如果获得最大距离的点在直线的上方，那么θ为正，如果获得最大距离的点在直线的下方，那么θ为负；
高斯滤波单元核函数的标准差是由候选特征区域对应矩形的高与宽的比值来决定，即标准差其中，H表示候选特征区域对应矩形的高度，W表示候选特征区域对应矩形的宽度。
特征矩阵生成模块，用于对压缩后的候选特征区域提取特征点，并将各个候选特征区域提取的特征点归分并为特征矩阵，分别形成参考图像特征矩阵和待配准图像特征矩阵；特征点提取单元采用SIFT算法进行特征点的提取。
特征点匹配模块，用于对两个特征矩阵进行特征点匹配，并滤除误匹配点；特征点匹配单元采用FLANN算法进行特征点匹配，误匹配点滤除单元采用串行RANSAC算法进行误匹配点滤除。
配准结果显示模块，用于将参考图像和待配准图像按照特征点匹配后得到的变换矩阵重新采样到同一坐标系，其中参考图像保持不变，待配准图像作半透明处理，并使用棋盘格的方式进行配准结果显示。
图像输入模块、区域图像处理生成模块、网络图像处理生成模块、交点判定模块、区域截取模块、区域处理模块、特征矩阵生成模块、特征点匹配模块以及配准结果显示模块的具体处理过程均是采用一种具有仿射不变性的视网膜图像配准方法来实现的。
以上所述的具体实施方式，对本发明的技术领域、背景、目的、方案和有益效果做了进一步的详细说明，所应理解的是，本实施方式仅为本发明的优选方式而已，并不用于限制本发明，凡在本发明的精神和原则之内，做出的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种具有仿射不变性的视网膜图像配准方法，其特征在于，包括如下步骤：
步骤1：获取血管主要区域图像；
输入两幅视网膜图像，一幅作为参考图像Isrc，另一幅作为待配准图像Idst，然后分别对两幅图像依次进行提取绿色通道、高斯滤波以及多次底帽变换处理，最后将所有底帽变换所得结果叠加在一起即得到血管主要区域图像IEng；
步骤2：对血管主要区域图像IEng进行骨骼化处理，得到骨骼化网络I(x,y)，再对骨骼化网络I(x,y)进行去噪处理；
步骤3：判断各个像素点为交点还是端点；
计算骨骼化网络I(x,y)中各个像素点八连通区域内的通路数CN；
根据CN的取值将各个像素点划分为普通点POrd或者交点PCross或者端点PEnd，并形成所有交点集合和所有端点集合；
计算各个端点到与之相连通的交点的距离，如果该距离小于设定阈值，则删除该端点到与之相连通的交点所形成的分支，并分别在端点集合和交点集合中去除该端点和该交点；
步骤4：候选特征区域的截取；
依次在血管主要区域图像IEng中找到相互连通的两个交点，并截取这两个交点对应的连通域作为候选特征区域，
候选特征区域大小为RS＝LD*2*D，其中，LD表示相互连通的两个交点间的距离，D为连通域内的点到相互连通的两个交点所成直线的最大距离；
步骤5：对候选特征区域依次进行旋转、高斯滤波以及压缩处理；
旋转方向为顺时针，旋转角度为θ，其中，θ为候选特征区域内对应的相互连通的两个交点所成直线的倾斜角；
高斯滤波沿着X轴方向进行；
压缩处理是将候选特征区域的矩形区域变成正方形区域；
步骤6：提取特征点，形成特征矩阵；
对压缩后的候选特征区域提取特征点，将特征点坐标变回原始视网膜图像坐标；将各个候选特征区域提取的特征点归并为特征矩阵，分别形成参考图像Isrc特征矩阵Fsrc、待配准图像Idst特征矩阵Fdst；
步骤7：对Fsrc与Fdst进行特征点匹配，并滤除误匹配的点；
步骤8：显示配准结果；
将参考图像Isrc和待配准图像Idst按照特征点匹配后得到的变换矩阵重新采样到同一坐标系，其中参考图像Isrc保持不变，待配准图像Idst作半透明处理，并使用棋盘格的方式进行配准结果显示。
2.根据权利要求1所述的一种具有仿射不变性的视网膜图像配准方法，其特征在于，所述步骤1中，多次底帽变换采用线性结构元素，底帽变换的过程由公式表示为：
Ibot＝IGau-IGau×S
其中，Ibot表示底帽变换所得结果，IGau表示经高斯滤波后的图像，S表示线性结构元素。
3.根据权利要求2所述的一种具有仿射不变性的视网膜图像配准方法，其特征在于，所述线性结构元素为一组具有不同尺度l和不同方向θl的线性结构元素，其中，l为整数，且取值范围为：2≤l≤15，
针对每一个尺度，首先把该尺度所有方向上的底帽变换所得结果叠加在一起，然后通过OSTU阈值分割方法得到该尺度上的血管分割结果，最后将所有尺度上的血管分割结果叠加在一起得到血管主要区域图像IEng。
4.根据权利要求1所述的一种具有仿射不变性的视网膜图像配准方法，其特征在于，所述步骤2中，去噪处理的具体操作为：
在八连通区域模式下，对骨骼化网络I(x,y)的连通域进行标记，然后统计各个连通域内像素点的个数，如果连通域内像素点的个数小于设定值，则去除该连通；
设定值的取值为100。
5.根据权利要求1所述的一种具有仿射不变性的视网膜图像配准方法，其特征在于，所述步骤3的具体操作过程为：
设P为待检测像素点，P在八连通区域内的点记为P1,P2,L,P8，P1,P2,L,P8的排列顺序是以P1为起始点，顺时针依次排列，且P1为P正上方的点，CN的计算公式表示为：
如果则P为端点；如果则P为普通点；如果则P为交点。
6.根据权利要求1所述的一种具有仿射不变性的视网膜图像配准方法，其特征在于，所述步骤4的具体操作如下：
首先，在血管主要区域图像IEng中对所有交点进行连通域检测，如果这两个交点是相互连通的，则归为一类，并记录连通域Con，直到所有交点分类完毕，得到K个连通域，记为：Conk为第k类相互连通的交点对应的连通域；
其次，找到连通域Conk内相互连通的两个交点Pc1(x1,y1),Pc2(x2,y2)，并计算所述两个交点所成直线
再次，找到连通域Conk内到直线lc1c2距离最大的点Pmax(x3,y3)，计算最大距离其中，A、B、C是直线lc1c2的直线方程对应的系数；
最后，截取两个交点Pc1(x1,y1),Pc2(x2,y2)对应的连通域Conk作为候选特征区域，候选特征区域大小为RS＝LD*2*D。
7.根据权利要求1所述的一种具有仿射不变性的视网膜图像配准方法，其特征在于，所述步骤5中的旋转操作是以候选特征区域左下角的交点作为中心点；
高斯滤波核函数的标准差是由候选特征区域对应矩形的高与宽的比值来决定，即标准差其中，H表示候选特征区域对应矩形的高度，W表示候选特征区域对应矩形的宽度。
8.根据权利要求1所述的一种具有仿射不变性的视网膜图像配准方法，其特征在于，所述步骤6中，特征点的提取采用SIFT算法。
9.根据权利要求1所述的一种具有仿射不变性的视网膜图像配准方法，其特征在于，所述步骤7中，特征点的匹配采用FLANN算法，滤除误匹配采用串行RANSAC算法。
10.一种具有仿射不变性的视网膜图像配准装置，其特征在于，包括：
图像输入模块，用于输入视网膜图像，一幅作为参考图像，另一幅作为待配准图像；
区域图像处理生成模块，用于分别对两幅图像依次进行绿色通道提取、高斯滤波以及多次底帽变换处理，并将所有底帽变换所得结果叠加在一起生成血管主要区域图像；
网络图像处理生成模块，用于对血管主要区域图像进行骨骼化处理，并生成骨骼化网络；
交点判定模块，用于根据骨骼化网络中各个像素点八连通区域内的通路数判断各个像素点是交点还是端点，并形成所有交点集合和所有端点集合；还用于计算各个端点到与之相连通的交点的距离，如果该距离小于设定阈值，则删除该端点到与之相连通的交点所形成的分支，并分别在端点集合和交点集合中去除该端点和该交点；
区域截取模块，用于依次在血管主要区域图像中找到相互连通的两个交点，并截取这两个交点对应的连通域作为候选特征区域，候选特征区域大小为RS＝LD*2*D，
其中，LD表示相互连通的两个交点间的距离，D为连通域内的点到相互连通的两个交点所成直线的最大距离；
区域处理模块，用于对候选特征区域依次进行顺时针旋转θ度，沿X轴方向高斯滤波，以及压缩至候选特征区域的矩形区域变成正方形区域，其中，θ为候选特征区域内对应的相互连通的两个交点所成直线的倾斜角；
特征矩阵生成模块，用于对压缩后的候选特征区域提取特征点，并将各个候选特征区域提取的特征点归分并为特征矩阵，分别形成参考图像特征矩阵和待配准图像特征矩阵；
特征点匹配模块，用于对两个特征矩阵进行特征点匹配，并滤除误匹配点；
配准结果显示模块，用于将参考图像和待配准图像按照特征点匹配后得到的变换矩阵重新采样到同一坐标系，其中参考图像保持不变，待配准图像作半透明处理，并使用棋盘格的方式进行配准结果显示；
所述的图像输入模块、区域图像处理生成模块、网络图像处理生成模块、交点判定模块、区域截取模块、区域处理模块、特征矩阵生成模块、特征点匹配模块以及配准结果显示模块的具体处理过程均是采用权利要求1-9任一项所述的一种具有仿射不变性的视网膜图像配准方法来实现的。
CN201810305736.8A 2018-04-08 2018-04-08 一种具有仿射不变性的视网膜图像配准方法及其装置 Active CN110363738B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810305736.8A CN110363738B (zh) 2018-04-08 2018-04-08 一种具有仿射不变性的视网膜图像配准方法及其装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810305736.8A CN110363738B (zh) 2018-04-08 2018-04-08 一种具有仿射不变性的视网膜图像配准方法及其装置

Publications (2)

Publication Number Publication Date
CN110363738A true CN110363738A (zh) 2019-10-22
CN110363738B CN110363738B (zh) 2021-08-27

Family

ID=68213687

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810305736.8A Active CN110363738B (zh) 2018-04-08 2018-04-08 一种具有仿射不变性的视网膜图像配准方法及其装置

Country Status (1)

Country Link
CN (1) CN110363738B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114115212A (zh) * 2020-08-26 2022-03-01 宁波方太厨具有限公司 一种清洁机器人的定位方法及采用该方法的清洁机器人

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060285755A1 (en) * 2005-06-16 2006-12-21 Strider Labs, Inc. System and method for recognition in 2D images using 3D class models
CN101093539A (zh) * 2007-07-27 2007-12-26 哈尔滨工程大学 手指静脉特征提取与匹配识别方法
US20130064436A1 (en) * 2011-05-10 2013-03-14 Olympus Medical Systems Corp. Medical image processing apparatus and method of operating medical image processing apparatus
CN103400384A (zh) * 2013-07-22 2013-11-20 西安电子科技大学 结合区域匹配和点匹配的大视角图像匹配方法
CN103871061A (zh) * 2014-03-17 2014-06-18 电子科技大学 一种基于双目视觉的眼底图像处理方法
CN107564048A (zh) * 2017-09-25 2018-01-09 南通大学 基于分叉点特征配准方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060285755A1 (en) * 2005-06-16 2006-12-21 Strider Labs, Inc. System and method for recognition in 2D images using 3D class models
CN101093539A (zh) * 2007-07-27 2007-12-26 哈尔滨工程大学 手指静脉特征提取与匹配识别方法
US20130064436A1 (en) * 2011-05-10 2013-03-14 Olympus Medical Systems Corp. Medical image processing apparatus and method of operating medical image processing apparatus
CN103400384A (zh) * 2013-07-22 2013-11-20 西安电子科技大学 结合区域匹配和点匹配的大视角图像匹配方法
CN103871061A (zh) * 2014-03-17 2014-06-18 电子科技大学 一种基于双目视觉的眼底图像处理方法
CN107564048A (zh) * 2017-09-25 2018-01-09 南通大学 基于分叉点特征配准方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
WEI LI等: "A Fully Affine Invariant Feature Detector", 《PROCEEDINGS OF THE 21ST INTERNATIONAL CONFERENCE ON PATTERN RECOGNITION (ICPR2012)》 *
YILIU HANG等: "Retinal Image Registration Based on the Feature of Bifurcation Point", 《2017 10TH INTERNATIONAL CONGRESS ON IMAGE AND SIGNAL PROCESSING, BIOMEDICAL ENGINEERING AND INFORMATICS (CISP-BMEI)》 *
沈奔: "基于局部血管结构的眼底图像配准", 《中国优秀硕士学位论文全文数据库信息科技辑》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114115212A (zh) * 2020-08-26 2022-03-01 宁波方太厨具有限公司 一种清洁机器人的定位方法及采用该方法的清洁机器人

Also Published As

Publication number Publication date
CN110363738B (zh) 2021-08-27

Similar Documents

Publication Publication Date Title
Hermosilla et al. Total denoising: Unsupervised learning of 3D point cloud cleaning
CN106447708A (zh) 一种oct眼底图像数据配准方法
CN106651827B (zh) 一种基于sift特征的眼底图像配准方法
CN106339998B (zh) 基于对比度金字塔变换的多聚焦图像融合方法
CN104835150B (zh) 一种基于学习的眼底血管几何关键点图像处理方法及装置
CN106127849B (zh) 三维精细血管重建方法及其系统
CN106228528B (zh) 一种基于决策图与稀疏表示的多聚焦图像融合方法
CN113239755B (zh) 一种基于空谱融合深度学习的医学高光谱图像分类方法
CN101667289B (zh) 基于nsct特征提取和监督分类的视网膜图像分割方法
CN110033465A (zh) 一种应用于双目内窥镜医学图像的实时三维重建方法
CA2612223A1 (en) System and method for computer aided polyp detection
CN105869166B (zh) 一种基于双目视觉的人体动作识别方法及系统
CN106934761A (zh) 一种三维非刚性光学相干断层扫描图像的配准方法
CN109166104A (zh) 一种病变检测方法、装置及设备
CN103942772A (zh) 一种多模态多维度的血管融合方法及系统
Bartesaghi et al. An energy-based three-dimensional segmentation approach for the quantitative interpretation of electron tomograms
CN110188767A (zh) 基于深度神经网络的角膜病图像序列化特征提取与分类方法及装置
CN109658423A (zh) 一种彩色眼底图的视盘视杯自动分割方法
CN109816617A (zh) 基于导向滤波和图论显著性的多模态医学图像融合方法
Hacihaliloglu et al. Statistical shape model to 3D ultrasound registration for spine interventions using enhanced local phase features
CN110046555A (zh) 内窥镜系统视频稳像方法及装置
Cao et al. An efficient lens structures segmentation method on as-oct images
CN107146211A (zh) 基于线扩散函数和双边滤波的视网膜血管图像降噪方法
Yan et al. Estimating fiber orientation distribution from diffusion MRI with spherical needlets
Wang et al. Integration of global and local features for specular reflection inpainting in colposcopic images

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