基于改进Harris角点检测的视网膜图像配准

汪 立，蒋念平

(上海理工大学 光电信息与计算机工程学院，上海 200093)

基于改进Harris角点检测的视网膜图像配准

汪 立，蒋念平

(上海理工大学 光电信息与计算机工程学院，上海 200093)

基于Harris角点的视网膜图像配准中，Harris角点检测算法无法针对不同的图像设定通用的阈值，检测到的角点也可能存在过多过密的问题。针对这些问题，提出一种改进的Harris角点检测算法，利用角点响应函数的均值来得到自适应阈值，并通过自动控制高斯模糊窗口大小，使角点数量在合理的范围。该算法采用Matlab语言来编程实现。实验结果表明，该方法能较好地配准各种视网膜图像，配准速度约为通用双引导就近点搜索算法的51%。

Harris角点；视网膜图像；局部特征； 配准；Matlab

视网膜图像配准是将不同时间、不同成像设备或不同条件下获取的两幅或多幅视网膜图像进行空间对齐的过程。

通过对视网膜图像进行配准，眼科医生能够更好的诊断各种眼科疾病，如老年黄斑变性，青光眼[1]等。视网膜图像配准的方法可以分为基于区域匹配和基于特征匹配的方法[2-3]。基于特征匹配的方法由于操作简单、配准速度较快、精度较高，是目前应用较多的配准方法。

特征检测是基于视网膜图像匹配的第一步，也是一个重要环节。以往大多使用视网膜血管节点作为图像的特征点[4]，但由于各种眼科疾病对视网膜的影响，使得需要配准的视网膜图像血管模糊，很难提取足够多的血管节点，从而使配准不成功。

Harris角点检测算子是一种常用的特征点检测算子，具有很好的旋转不变性，同时对亮度和对比度的变化也不敏感，文中选用其检测的角点来代替血管节点作为特征点，不仅可检测到更多的特征点，且检测到的特征点较为稳定。但Harris角点检测需要手动设定一个阈值，该阈值大小对角点检测结果影响较大。因此，本文通过改进Harris角点检测算子，使其不需要手动设定阈值，同时通过控制检测到的角点数量来提高配准速度或准确度。检测完角点后，本文选用与目前比较成熟的通用双引导就近点搜索算法(GDB-ICP)不同的方法来描述角点局部特征，匹配特征点来实现视网膜图像配准。

1 改进的Harris角点检测算法

1.1 Harris角点检测算法原理

Harris角点检测算法是由C. Harris和M. J. Stephens在1988年提出的。人眼对角点的识别通常是在一个较小的区域来完成的。若在各个方向上移动这个特定的小窗口，窗口内灰度值都发生了较大变化，则认为窗口内存在角点。若在某一个方向移动这个特定的小窗口，窗口内灰度值才有较大变化，则窗口内可能为一条直线的部分线段。如果在各个方向移动这个特定的小窗口，窗口内灰度值都没有变化，则认为窗口内没有角点。

Harris在其文献[5]内给出的角点计算方法并不计算具体的特征值，而是通过计算角点响应R，其计算表达式如下

R=detM-μ(traceM)2

(1)

其中，det和trace分别为M的行列式和迹； 为一个常数，通常取值0.04～0.06。若R大于给定的阈值 ，则这一点就被检测为角点。

Harris给出的角点响应函数含有参数μ，μ的取值大小会对R值大小产生影响，因此，本文选用新的角点响应函数[6]，该角点响应函数避免了参数μ的选择，其表达式如下

(2)

式中， 为一个极小量，用以防止分母为0。

1.2 Harris角点检测算法的改进

在利用Harris算法检测视网膜图像角点时，用来与角点响应函数R比较的阈值T是提前手动设定的。对于不同的图像，为了得到良好的角点检测效果，可能需要设定不同的阈值T。考虑到算法的通用性，不可能在每张图片处理时都进行多次试验来比较得到最佳阈值。考虑到设定的阈值T是用来与角点响应函数R比较大小，从而得到角点。由此想到可利用角点响应函数R来得到阈值，这样便可不用手动去设置。本文选用角点响应函数R的平均值Rave的k倍作为阈值。通过实验知道k取0.1时，检测效果较好。Harris角点检测的步骤如下：

(1)计算图像I(x,y)在X和Y轴方向的梯度Ix，Iy；

(2)计算图像两个方向梯度的乘积

(3)

(3)使用高斯模糊函数对第2步得到的三个值进行加权，生成矩阵M；

(4)计算每个像点的角点响应R，并计算其平均值Rave以及阈值T；

(5)将R中

(4)

(6)使用s×s(一般为3×3或5×5)大小窗口进行非最大值压制，窗口内最大值就是角点；

(7)提取角点坐标。

在图像配准中，特征描述耗时较大。用Harris检测到的角点越多，进行特征描述时耗时越长。本文在实验中注意到，视网膜图像检测的特征点多于一定数量时，再增加角点数对图像配准精度的提升没有明显的帮助，当检测的角点太少时，对图像配准精度影响较大，甚至导致配准失败。因此，可通过控制角点的数量来提高配准的速度和准确度。当检测到的角点多于500个时，可通过逐步增大进行高斯模糊的窗口W的大小，使检测的角点在500以内，这样便可提高配准速度，同时对配准精度影响也不大，当检测到的角点少于200个时，通过逐步减少W的大小，增加检测的角点数量。当然窗口W也不能无限小，否则就没有意义。所以当W的半径<1时，停止减少W，即使此时角点数仍<200个。

2 图像配准

图像配准的第一步特征点检测在上节中已介绍，下文介绍了图像配准的其他步骤，图像配准步骤如图1所示。

图1 配准流程图

2.1 特征点描述、匹配与提纯

在描述特征点之前，需要给特征点分配一个主方向。Sift算法运用方向直方图[7]来计算主方向。在多模图像中，通过方向直方图计算出来的主方向可能指向不相关的方向，并且方向直方图是离散的，导致产生错误的配准点对。本文采用一个连续的方法，平均平方梯度[8-9]来计算每个特征点的主方向，该方法相比方向直方图计算效率更高，结果更准确。

给每个特征点分配一个主方向后，接下来就要提取特征点局部特征。文中运用Jian Chen等人提出的对称描述子[9]来提取局部特征。生成对称描述子的主要步骤如下：

(1)将坐标轴旋转为关键点的方向，以确保旋转不变性；

(2)以关键点为中心取16×16邻域窗口，并分解为16个4×4的子窗口，分别计算每个像素点的梯度幅值和方向，方向限定在 ,然后对每个梯度幅值乘以对应的高斯权重参数；

(3)统计每个4×4的子窗口的方向直方图，每个直方图含有8柱，每22.5°一个柱，生成4×4×8共128维的第一个子描述子H。将H旋转180°得到第二个子描述子Q。最后将两个子描述子按照下面的公式结合起来得到最终的描述子Des

(5)

其中，c1，c2是用来调整对称描述子幅值的比例的参数。

提取局部特征后，下一步进行特征点对匹配，也就是找出参考图像与待配准图像特征点之间的对应关系。本文选用Silpa-Anan等提出的随机kd树算法来搜索对应的特征点，具体算法请参考文献[10～11]。

随机kd树搜索到的匹配点对并不能保证100%的正确，为了提纯匹配点，剔除错误的匹配点，选用RANSAC算法[12]来排除错误的匹配点。

2.3 选择变换方法来配准图像

计算变换参数常用线性正投影变换，仿射变换，二阶多项式变换。优先选用更高阶的变换，但越高阶需要的控制点越多。因此，本文根据控制点的多少来选择变换方法。当控制点为2个时选择线性正投影变换，当控制点数>2个并<6个时，选择仿射变换，当控制点≥6个时选择二阶多项式变换。这样可在保证计算变换参数成功的前提下尽量提高精确度。

3 实验结果及分析

为了验证本算法的配准效果，选取多组不同的视网膜图像进行配准，并与经典的GDB-ICP配准算法进行对比。实验平台硬件环境采用Intel(R) Core(TM) i3-2310M CPU，2.10 GHz，4 GB内存的笔记本电脑，操作系统为Windows 7 32位，算法采用Matlab语言来实现[13-15]。如图2是选用的3组视网膜图片，第1组为视角稍有不同的两张视网膜图片，分辨率为403×384，如图2(a)和图2(b)所示。第2组为相对旋转的两张视网膜图片，分辨率为700×605，如图2(c)和图2(d)所示。第3组为血管模糊的两张视网膜图片，分辨率为239×225，如图2(e)和图2(f)所示。

图2 3组实验图片

首先，对比改进Harris角点检测算法前后，检测到的角点数，匹配点对数及配准时间，除角点检测步骤外，其他配准步骤相同。从表1中数据可看出，改进的Harris角点检测算法在角点少于200时，可增加角点数，提高匹配点数，当角点过多时，可限制角点数，防止角点过于多和密，提高配准速度。

表1 Harris角点检测算法改进前后检测到的角点数和配准时间

其次，用本文实现的改进Harris后的配准算法与GDB-ICP算法来配准图2中的3组图像，本文算法平均用时22.98 s，GDB-ICP平均用时42 s，配准后融合结果如图3所示。图3(a)～图3(c)分别为第1组～第3组视网膜图像通过本文方法配准后融合的结果，图3(d)～图3(f)分别为第1组～第3组视网膜图像通过GDB-ICP方法配准后融合的结果。

图3 3组图片通过本文方法和GDB-ICP算法配准后融合的结果

4 结束语

本文通过改进Harris角点检测算法，使得无需针对不同视网膜图片去设定阈值T，提高了算法的通用性，同时通过灵活控制高斯模糊函数窗口的大小，控制特征点数量在合适的范围，在对精度影响不大的情况下，提高了配准速度。实验表明，该方法配准视网膜图片效果较好，对模糊的视网膜图片也可成功配准，配准速度也比GDB-ICP算法有了提升。但本算法对于旋转角度很大的视网膜图片还不能配准，下一步将对这一问题进行研究。

[1] Sanchez-Galeana C,Bowd C,Blumenthal E Z, et al. Using optical imaging summary data to detect glaucoma[J].Ophthalmology,2001, 108(10):1812 -1818.

[2] Saxena S,Singh R K.A survey of recent and classical image registration methods[J]. International Journal of Signal Processing Image Processing & Pattern Recognition,2014,7(4):167-176.

[3] Zitová B,Flusser J. Image registration methods: a survey[J].Image & Vision Computing,2010,21(11):977-1000.

[4] 赵晓芳,林土胜.视网膜血管图像特征点自动提取和分类[J].计算机工程与应用,2011,47(8):14-17.

[5] Harris C,Stephens M J.A combined corner and edge detector[C].Manchester: Proceedings of the 4th Alvey Vision Conference,1988.

[6] Smith A M,Brady J M.SUSAN: a new approach to low level image processing [J].International Journal of Computer Vision,1997,23(1):45-78.

[7] Lowe D G.Distinctive image features from scale-invariant keypoints[J]. International Journal of Computer Vision,2004,60(2):91-110.

[8] Hong L,Wan Y,Jain A. Fingerprint image enhancement:algorithm and performance evaluation[J].IEEE Transactions on Pattern Analysis & Machine Intelligence,1998,20(8):777-789.

[9] Chen J,Smith R,Tian J,et al.A novel registration method for retinal images based on local features[C].Shanghai：International Conference of the IEEE Engineering in Medicine & Biology Society IEEE Engineering in Medicine & Biology Society Conference，2008.

[10] Silpaanan C,Hartley R.Optimised KD-trees for fast image descriptor matching[C].Alaska:IEEE Conference on Computer Vision & Pattern Recognition,IEEE,2008.

[11] 丁南南,刘艳滢,张叶,等.基于SURF-DAISY算法和随机kd树的快速图像配准[J].光电子.激光,2012(7):1395-1402.

[12] Huang Y Q,Yu F U,Guang-Kun M A. Cylindrical panoramic image stitching method based on RANSAC algorithm[J].Journal of Shenyang University of Technology,2008,30(4):461-465.

[13] Rafael C Gonzalez,Richard E Woods,Steven L Eddins.数字图像处理:Matlab版[M].2版.阮秋琦,译.北京:电子工业出版社,2013.

[14] 陈显毅.图像配准技术及其Matlab编程实现[M].北京:电子工业出版社,2009.

[15] 邓兵华,张明.Matlab环境下图像配准方法[J].电子科技,2014,27(1):139-141.

Retinal Image Registration Based on Improved Harris Corner Detection

WANG Li，JIANG Nianping

(School of Optical-Electrical and Computer Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China)

In the process of retinal image registration based on Harris corner detection, the Harris corner detection algorithm fails to set a universal threshold for each image, and it is a problem that corner detected may be too many. To solve these problems, an improved Harris corner detection algorithm is proposed. This algorithm uses the mean of corner response function to get an auto-adaptive threshold and make the number of corner points within a reasonable range by automatically adjusting the size of the Gaussian blur window. The algorithm is implemented in Matlab language. Experimental results show that this method can register a variety of retinal images at a registration rate of about 51% of that by the dual bootstrap-ICP algorithm.

Harris corner point; retinal image; local feature; registration; Matlab

2016- 04- 07

汪立(1990-)，男，硕士研究生。研究方向：图像处理等。蒋念平(1957-)，男，副教授。研究方向：计算机应用等。

10.16180/j.cnki.issn1007-7820.2017.02.031

TP391.41

A

1007-7820(2017)02-119-04