基于改進(jìn)Harris角點(diǎn)檢測(cè)的視網(wǎng)膜圖像配準(zhǔn)

汪 立，蔣念平

(上海理工大學(xué) 光電信息與計(jì)算機(jī)工程學(xué)院，上海 200093)

基于改進(jìn)Harris角點(diǎn)檢測(cè)的視網(wǎng)膜圖像配準(zhǔn)

汪 立，蔣念平

(上海理工大學(xué) 光電信息與計(jì)算機(jī)工程學(xué)院，上海 200093)

基于Harris角點(diǎn)的視網(wǎng)膜圖像配準(zhǔn)中，Harris角點(diǎn)檢測(cè)算法無(wú)法針對(duì)不同的圖像設(shè)定通用的閾值，檢測(cè)到的角點(diǎn)也可能存在過多過密的問題。針對(duì)這些問題，提出一種改進(jìn)的Harris角點(diǎn)檢測(cè)算法，利用角點(diǎn)響應(yīng)函數(shù)的均值來得到自適應(yīng)閾值，并通過自動(dòng)控制高斯模糊窗口大小，使角點(diǎn)數(shù)量在合理的范圍。該算法采用Matlab語(yǔ)言來編程實(shí)現(xiàn)。實(shí)驗(yàn)結(jié)果表明，該方法能較好地配準(zhǔn)各種視網(wǎng)膜圖像，配準(zhǔn)速度約為通用雙引導(dǎo)就近點(diǎn)搜索算法的51%。

Harris角點(diǎn)；視網(wǎng)膜圖像；局部特征； 配準(zhǔn)；Matlab

視網(wǎng)膜圖像配準(zhǔn)是將不同時(shí)間、不同成像設(shè)備或不同條件下獲取的兩幅或多幅視網(wǎng)膜圖像進(jìn)行空間對(duì)齊的過程。

通過對(duì)視網(wǎng)膜圖像進(jìn)行配準(zhǔn)，眼科醫(yī)生能夠更好的診斷各種眼科疾病，如老年黃斑變性，青光眼[1]等。視網(wǎng)膜圖像配準(zhǔn)的方法可以分為基于區(qū)域匹配和基于特征匹配的方法[2-3]?；谔卣髌ヅ涞姆椒ㄓ捎诓僮骱?jiǎn)單、配準(zhǔn)速度較快、精度較高，是目前應(yīng)用較多的配準(zhǔn)方法。

特征檢測(cè)是基于視網(wǎng)膜圖像匹配的第一步，也是一個(gè)重要環(huán)節(jié)。以往大多使用視網(wǎng)膜血管節(jié)點(diǎn)作為圖像的特征點(diǎn)[4]，但由于各種眼科疾病對(duì)視網(wǎng)膜的影響，使得需要配準(zhǔn)的視網(wǎng)膜圖像血管模糊，很難提取足夠多的血管節(jié)點(diǎn)，從而使配準(zhǔn)不成功。

Harris角點(diǎn)檢測(cè)算子是一種常用的特征點(diǎn)檢測(cè)算子，具有很好的旋轉(zhuǎn)不變性，同時(shí)對(duì)亮度和對(duì)比度的變化也不敏感，文中選用其檢測(cè)的角點(diǎn)來代替血管節(jié)點(diǎn)作為特征點(diǎn)，不僅可檢測(cè)到更多的特征點(diǎn)，且檢測(cè)到的特征點(diǎn)較為穩(wěn)定。但Harris角點(diǎn)檢測(cè)需要手動(dòng)設(shè)定一個(gè)閾值，該閾值大小對(duì)角點(diǎn)檢測(cè)結(jié)果影響較大。因此，本文通過改進(jìn)Harris角點(diǎn)檢測(cè)算子，使其不需要手動(dòng)設(shè)定閾值，同時(shí)通過控制檢測(cè)到的角點(diǎn)數(shù)量來提高配準(zhǔn)速度或準(zhǔn)確度。檢測(cè)完角點(diǎn)后，本文選用與目前比較成熟的通用雙引導(dǎo)就近點(diǎn)搜索算法(GDB-ICP)不同的方法來描述角點(diǎn)局部特征，匹配特征點(diǎn)來實(shí)現(xiàn)視網(wǎng)膜圖像配準(zhǔn)。

1 改進(jìn)的Harris角點(diǎn)檢測(cè)算法

1.1 Harris角點(diǎn)檢測(cè)算法原理

Harris角點(diǎn)檢測(cè)算法是由C. Harris和M. J. Stephens在1988年提出的。人眼對(duì)角點(diǎn)的識(shí)別通常是在一個(gè)較小的區(qū)域來完成的。若在各個(gè)方向上移動(dòng)這個(gè)特定的小窗口，窗口內(nèi)灰度值都發(fā)生了較大變化，則認(rèn)為窗口內(nèi)存在角點(diǎn)。若在某一個(gè)方向移動(dòng)這個(gè)特定的小窗口，窗口內(nèi)灰度值才有較大變化，則窗口內(nèi)可能為一條直線的部分線段。如果在各個(gè)方向移動(dòng)這個(gè)特定的小窗口，窗口內(nèi)灰度值都沒有變化，則認(rèn)為窗口內(nèi)沒有角點(diǎn)。

Harris在其文獻(xiàn)[5]內(nèi)給出的角點(diǎn)計(jì)算方法并不計(jì)算具體的特征值，而是通過計(jì)算角點(diǎn)響應(yīng)R，其計(jì)算表達(dá)式如下

R=detM-μ(traceM)2

(1)

其中，det和trace分別為M的行列式和跡； 為一個(gè)常數(shù)，通常取值0.04～0.06。若R大于給定的閾值 ，則這一點(diǎn)就被檢測(cè)為角點(diǎn)。

Harris給出的角點(diǎn)響應(yīng)函數(shù)含有參數(shù)μ，μ的取值大小會(huì)對(duì)R值大小產(chǎn)生影響，因此，本文選用新的角點(diǎn)響應(yīng)函數(shù)[6]，該角點(diǎn)響應(yīng)函數(shù)避免了參數(shù)μ的選擇，其表達(dá)式如下

(2)

式中， 為一個(gè)極小量，用以防止分母為0。

1.2 Harris角點(diǎn)檢測(cè)算法的改進(jìn)

在利用Harris算法檢測(cè)視網(wǎng)膜圖像角點(diǎn)時(shí)，用來與角點(diǎn)響應(yīng)函數(shù)R比較的閾值T是提前手動(dòng)設(shè)定的。對(duì)于不同的圖像，為了得到良好的角點(diǎn)檢測(cè)效果，可能需要設(shè)定不同的閾值T?？紤]到算法的通用性，不可能在每張圖片處理時(shí)都進(jìn)行多次試驗(yàn)來比較得到最佳閾值。考慮到設(shè)定的閾值T是用來與角點(diǎn)響應(yīng)函數(shù)R比較大小，從而得到角點(diǎn)。由此想到可利用角點(diǎn)響應(yīng)函數(shù)R來得到閾值，這樣便可不用手動(dòng)去設(shè)置。本文選用角點(diǎn)響應(yīng)函數(shù)R的平均值Rave的k倍作為閾值。通過實(shí)驗(yàn)知道k取0.1時(shí)，檢測(cè)效果較好。Harris角點(diǎn)檢測(cè)的步驟如下：

(1)計(jì)算圖像I(x,y)在X和Y軸方向的梯度Ix，Iy；

(2)計(jì)算圖像兩個(gè)方向梯度的乘積

(3)

(3)使用高斯模糊函數(shù)對(duì)第2步得到的三個(gè)值進(jìn)行加權(quán)，生成矩陣M；

(4)計(jì)算每個(gè)像點(diǎn)的角點(diǎn)響應(yīng)R，并計(jì)算其平均值Rave以及閾值T；

(5)將R中

(4)

(6)使用s×s(一般為3×3或5×5)大小窗口進(jìn)行非最大值壓制，窗口內(nèi)最大值就是角點(diǎn)；

(7)提取角點(diǎn)坐標(biāo)。

在圖像配準(zhǔn)中，特征描述耗時(shí)較大。用Harris檢測(cè)到的角點(diǎn)越多，進(jìn)行特征描述時(shí)耗時(shí)越長(zhǎng)。本文在實(shí)驗(yàn)中注意到，視網(wǎng)膜圖像檢測(cè)的特征點(diǎn)多于一定數(shù)量時(shí)，再增加角點(diǎn)數(shù)對(duì)圖像配準(zhǔn)精度的提升沒有明顯的幫助，當(dāng)檢測(cè)的角點(diǎn)太少時(shí)，對(duì)圖像配準(zhǔn)精度影響較大，甚至導(dǎo)致配準(zhǔn)失敗。因此，可通過控制角點(diǎn)的數(shù)量來提高配準(zhǔn)的速度和準(zhǔn)確度。當(dāng)檢測(cè)到的角點(diǎn)多于500個(gè)時(shí)，可通過逐步增大進(jìn)行高斯模糊的窗口W的大小，使檢測(cè)的角點(diǎn)在500以內(nèi)，這樣便可提高配準(zhǔn)速度，同時(shí)對(duì)配準(zhǔn)精度影響也不大，當(dāng)檢測(cè)到的角點(diǎn)少于200個(gè)時(shí)，通過逐步減少W的大小，增加檢測(cè)的角點(diǎn)數(shù)量。當(dāng)然窗口W也不能無(wú)限小，否則就沒有意義。所以當(dāng)W的半徑<1時(shí)，停止減少W，即使此時(shí)角點(diǎn)數(shù)仍<200個(gè)。

2 圖像配準(zhǔn)

圖像配準(zhǔn)的第一步特征點(diǎn)檢測(cè)在上節(jié)中已介紹，下文介紹了圖像配準(zhǔn)的其他步驟，圖像配準(zhǔn)步驟如圖1所示。

圖1 配準(zhǔn)流程圖

2.1 特征點(diǎn)描述、匹配與提純

在描述特征點(diǎn)之前，需要給特征點(diǎn)分配一個(gè)主方向。Sift算法運(yùn)用方向直方圖[7]來計(jì)算主方向。在多模圖像中，通過方向直方圖計(jì)算出來的主方向可能指向不相關(guān)的方向，并且方向直方圖是離散的，導(dǎo)致產(chǎn)生錯(cuò)誤的配準(zhǔn)點(diǎn)對(duì)。本文采用一個(gè)連續(xù)的方法，平均平方梯度[8-9]來計(jì)算每個(gè)特征點(diǎn)的主方向，該方法相比方向直方圖計(jì)算效率更高，結(jié)果更準(zhǔn)確。

給每個(gè)特征點(diǎn)分配一個(gè)主方向后，接下來就要提取特征點(diǎn)局部特征。文中運(yùn)用Jian Chen等人提出的對(duì)稱描述子[9]來提取局部特征。生成對(duì)稱描述子的主要步驟如下：

(1)將坐標(biāo)軸旋轉(zhuǎn)為關(guān)鍵點(diǎn)的方向，以確保旋轉(zhuǎn)不變性；

(2)以關(guān)鍵點(diǎn)為中心取16×16鄰域窗口，并分解為16個(gè)4×4的子窗口，分別計(jì)算每個(gè)像素點(diǎn)的梯度幅值和方向，方向限定在 ,然后對(duì)每個(gè)梯度幅值乘以對(duì)應(yīng)的高斯權(quán)重參數(shù)；

(3)統(tǒng)計(jì)每個(gè)4×4的子窗口的方向直方圖，每個(gè)直方圖含有8柱，每22.5°一個(gè)柱，生成4×4×8共128維的第一個(gè)子描述子H。將H旋轉(zhuǎn)180°得到第二個(gè)子描述子Q。最后將兩個(gè)子描述子按照下面的公式結(jié)合起來得到最終的描述子Des

(5)

其中，c1，c2是用來調(diào)整對(duì)稱描述子幅值的比例的參數(shù)。

提取局部特征后，下一步進(jìn)行特征點(diǎn)對(duì)匹配，也就是找出參考圖像與待配準(zhǔn)圖像特征點(diǎn)之間的對(duì)應(yīng)關(guān)系。本文選用Silpa-Anan等提出的隨機(jī)kd樹算法來搜索對(duì)應(yīng)的特征點(diǎn)，具體算法請(qǐng)參考文獻(xiàn)[10～11]。

隨機(jī)kd樹搜索到的匹配點(diǎn)對(duì)并不能保證100%的正確，為了提純匹配點(diǎn)，剔除錯(cuò)誤的匹配點(diǎn)，選用RANSAC算法[12]來排除錯(cuò)誤的匹配點(diǎn)。

2.3 選擇變換方法來配準(zhǔn)圖像

計(jì)算變換參數(shù)常用線性正投影變換，仿射變換，二階多項(xiàng)式變換。優(yōu)先選用更高階的變換，但越高階需要的控制點(diǎn)越多。因此，本文根據(jù)控制點(diǎn)的多少來選擇變換方法。當(dāng)控制點(diǎn)為2個(gè)時(shí)選擇線性正投影變換，當(dāng)控制點(diǎn)數(shù)>2個(gè)并<6個(gè)時(shí)，選擇仿射變換，當(dāng)控制點(diǎn)≥6個(gè)時(shí)選擇二階多項(xiàng)式變換。這樣可在保證計(jì)算變換參數(shù)成功的前提下盡量提高精確度。

3 實(shí)驗(yàn)結(jié)果及分析

為了驗(yàn)證本算法的配準(zhǔn)效果，選取多組不同的視網(wǎng)膜圖像進(jìn)行配準(zhǔn)，并與經(jīng)典的GDB-ICP配準(zhǔn)算法進(jìn)行對(duì)比。實(shí)驗(yàn)平臺(tái)硬件環(huán)境采用Intel(R) Core(TM) i3-2310M CPU，2.10 GHz，4 GB內(nèi)存的筆記本電腦，操作系統(tǒng)為Windows 7 32位，算法采用Matlab語(yǔ)言來實(shí)現(xiàn)[13-15]。如圖2是選用的3組視網(wǎng)膜圖片，第1組為視角稍有不同的兩張視網(wǎng)膜圖片，分辨率為403×384，如圖2(a)和圖2(b)所示。第2組為相對(duì)旋轉(zhuǎn)的兩張視網(wǎng)膜圖片，分辨率為700×605，如圖2(c)和圖2(d)所示。第3組為血管模糊的兩張視網(wǎng)膜圖片，分辨率為239×225，如圖2(e)和圖2(f)所示。

圖2 3組實(shí)驗(yàn)圖片

首先，對(duì)比改進(jìn)Harris角點(diǎn)檢測(cè)算法前后，檢測(cè)到的角點(diǎn)數(shù)，匹配點(diǎn)對(duì)數(shù)及配準(zhǔn)時(shí)間，除角點(diǎn)檢測(cè)步驟外，其他配準(zhǔn)步驟相同。從表1中數(shù)據(jù)可看出，改進(jìn)的Harris角點(diǎn)檢測(cè)算法在角點(diǎn)少于200時(shí)，可增加角點(diǎn)數(shù)，提高匹配點(diǎn)數(shù)，當(dāng)角點(diǎn)過多時(shí)，可限制角點(diǎn)數(shù)，防止角點(diǎn)過于多和密，提高配準(zhǔn)速度。

表1 Harris角點(diǎn)檢測(cè)算法改進(jìn)前后檢測(cè)到的角點(diǎn)數(shù)和配準(zhǔn)時(shí)間

其次，用本文實(shí)現(xiàn)的改進(jìn)Harris后的配準(zhǔn)算法與GDB-ICP算法來配準(zhǔn)圖2中的3組圖像，本文算法平均用時(shí)22.98 s，GDB-ICP平均用時(shí)42 s，配準(zhǔn)后融合結(jié)果如圖3所示。圖3(a)～圖3(c)分別為第1組～第3組視網(wǎng)膜圖像通過本文方法配準(zhǔn)后融合的結(jié)果，圖3(d)～圖3(f)分別為第1組～第3組視網(wǎng)膜圖像通過GDB-ICP方法配準(zhǔn)后融合的結(jié)果。

圖3 3組圖片通過本文方法和GDB-ICP算法配準(zhǔn)后融合的結(jié)果

4 結(jié)束語(yǔ)

本文通過改進(jìn)Harris角點(diǎn)檢測(cè)算法，使得無(wú)需針對(duì)不同視網(wǎng)膜圖片去設(shè)定閾值T，提高了算法的通用性，同時(shí)通過靈活控制高斯模糊函數(shù)窗口的大小，控制特征點(diǎn)數(shù)量在合適的范圍，在對(duì)精度影響不大的情況下，提高了配準(zhǔn)速度。實(shí)驗(yàn)表明，該方法配準(zhǔn)視網(wǎng)膜圖片效果較好，對(duì)模糊的視網(wǎng)膜圖片也可成功配準(zhǔn)，配準(zhǔn)速度也比GDB-ICP算法有了提升。但本算法對(duì)于旋轉(zhuǎn)角度很大的視網(wǎng)膜圖片還不能配準(zhǔn)，下一步將對(duì)這一問題進(jìn)行研究。

[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] 趙曉芳,林土勝.視網(wǎng)膜血管圖像特征點(diǎn)自動(dòng)提取和分類[J].計(jì)算機(jī)工程與應(yīng)用,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算法和隨機(jī)kd樹的快速圖像配準(zhǔn)[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.數(shù)字圖像處理:Matlab版[M].2版.阮秋琦,譯.北京:電子工業(yè)出版社,2013.

[14] 陳顯毅.圖像配準(zhǔn)技術(shù)及其Matlab編程實(shí)現(xiàn)[M].北京:電子工業(yè)出版社,2009.

[15] 鄧兵華,張明.Matlab環(huán)境下圖像配準(zhǔn)方法[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-)，男，副教授。研究方向：計(jì)算機(jī)應(yīng)用等。

10.16180/j.cnki.issn1007-7820.2017.02.031

TP391.41

A

1007-7820(2017)02-119-04