魏麗芳 潘 林 余 輪

（福州大學計算機圖象圖形研究所，福州 350002）

引言

眼底圖像配準技術(shù)在輔助眼科診斷和治療過程中有著廣泛的應(yīng)用［1-2］。通過配準眼底圖像，醫(yī)生能夠更好地診斷和檢測各種與眼底相關(guān)的疾病，如糖尿病、青光眼、黃斑部變性等，圖1給出兩幅病變眼底圖像示例。

近年來，各種眼底圖像配準技術(shù)已得到廣泛研究，目前眼底圖像配準方法主要分為基于灰度的方法和基于特征的方法［3］?；趫D像灰度的配準方法直接利用圖像的灰度值來確定配準的空間變換，其充分利用圖像中所包含的信息；但基于灰度的方法容易陷入一個局部極小值，而且變換模型的搜索空間是巨大的?；谔卣鞯姆椒ㄊ且源錅蕡D像對之間存在的對應(yīng)特征點（如血管分支和交叉點）為基礎(chǔ)，通過測量對應(yīng)點之間的相似性最值完成。鑒于兩種方法的特點，基于特征的眼底圖像配準算法得到越來越多的應(yīng)用。Can等提出層次匹配方法［2］，基于血管分支和交叉點，對于血管分支和交叉點的提取利用了血管追蹤和特征抽取的模式。Zana等利用形態(tài)學方法提取血管樹，將血管樹分支和交叉點（特征點）用鄰近血管的方向進行標注，實現(xiàn)每對匹配的特征點之間的角度不變性［4］。但這類方法仍然依賴于眼底血管結(jié)構(gòu)的提取。文獻［3］將適合曲面配準的迭代最近點iterative closest point，ICP）算法引入眼底圖像配準，提出雙重引導的迭代最近點算法（dual-bootstrap iterative closest point，DB-ICP）。該算法對于眼底的曲面結(jié)構(gòu)可以獲得精確的配準，在較好初值的前提下，算法具有較好的收斂性。但DB-ICP秉承了ICP算法耗時大的缺點，而且對于低質(zhì)量的眼底圖像，也很難配準成功。文獻［5］則將尺度不變性特征變換（scale invariable feature transformation，SIFT）引入眼底圖像配準，但使用的隨機抽樣一致性算法（random sample consensus，RANSAC）算法進行初始估計，它在去除誤匹配同時也去除了大量正確匹配特征。圖1為糖尿病誘發(fā)的眼底病變示例。

圖1 不同病變的眼底圖像。（a）增殖性糖尿病視網(wǎng)膜五期玻璃體病變；（b）非增殖性糖尿病視網(wǎng)膜三期出血灶病變Fig.1 Different pathological retinal images.（a）proliferative diabetic retinopathy Ⅴ vitreous gel pathologic；（b）nonproliferative diabetic retinopathy Ⅲ hemorrhagic foci pathologic

針對病變眼底圖像血管結(jié)構(gòu)模糊同時也具有復雜的球體解剖學結(jié)構(gòu)和光學圖像系統(tǒng)的特點，提出了一種基于不變特征的眼底圖像配準方法。首先提取SIFT特征作為匹配的特征點，提出了雙邊“或”的best-bin-first（BBF）特征匹配算法對特征點進行初始匹配，通過SIFT特征點的方向特征和空間幾何特性去除誤匹配，精化匹配特征對，并保留足夠的正確匹配對；根據(jù)獲得的匹配特征通過迭代加權(quán)最小二乘法和M-估計進行參數(shù)估計，實現(xiàn)病變眼底圖像的精確配準。

1 特征提取與匹配

1.1 特征點的提取

圖像不變特征（如，角點、SIFT、SURF）和描述符的不斷發(fā)展為圖像精確配準提供了有利的前提［6-9］。在之前的研究中已經(jīng)證明SIFT算法在眼底圖像配準中的有效性［5］。SIFT特征檢測算法詳見文獻［8-9］，是一種基于尺度空間的、對圖像縮放、旋轉(zhuǎn)甚至仿射變換保持不變性的圖像局部特征描述算子，該算法首先在尺度空間進行特征檢測，并確定關(guān)鍵點的位置和關(guān)鍵點所處的尺度，然后使用關(guān)鍵點鄰域梯度的主方向作為該點的方向特征，以實現(xiàn)算子對尺度和方向的無關(guān)性。其對旋轉(zhuǎn)、尺度縮放、亮度變化保持不變性，對視角變化、仿射變換、噪聲也保持一定程度的穩(wěn)定性，并對每一個關(guān)鍵點產(chǎn)生一個128維的特征描述符以作為匹配依據(jù)。

1.2 特征匹配

為了提高特征正確匹配率，同時確保具有足夠的匹配對應(yīng)用于圖像的兩兩配準，提出一種雙邊“或”的Best-Bin-First（BBF）算法建立兩個圖像特征點之間的相應(yīng)性，利用特征點的方向信息和空間幾何信息去除誤匹配，盡可能地保留正確的匹配。

根據(jù)特征點的鄰域信息建立高維特征描述符，通過最近鄰比值方法，由逼近高維空間的最近鄰搜索的BBF算法找到待匹配圖像中對應(yīng)的特征匹配點對［9-10］；BBF算法通過尋找最近鄰和次近鄰，同時對最近鄰與次近鄰的比值設(shè)定閾值T（文中T取0.5）。

設(shè)圖像I1和I2的SIFT描述符分別是DES1和DES2，計算相應(yīng)的距離［10］：

式中，·代表向量的點積。這個集合包含了des與I2中所有描述符之間的距離。dess對應(yīng)于Dis的最大元素，表示des的最近鄰。如果最近鄰距離與次近鄰距離的比值小于T，則說明構(gòu)成最近鄰的兩點是匹配的，否則就不匹配。該算法可有效達到穩(wěn)定的匹配：在匹配對中，正確的匹配對應(yīng)的最短距離應(yīng)該比錯誤匹配對應(yīng)的的最短距離短得多，而在高維特征空間中相似的距離內(nèi)可能還存在著一些其他的錯誤匹配，利用次近鄰匹配作為這部分空間錯誤匹配密度的估計。

以上匹配過程是單邊匹配的，然而在單邊匹配時，存在一些遺漏的正確匹配。對于具有曲面結(jié)構(gòu)的眼底圖像，正確的匹配對越多且分布越廣泛（重疊區(qū)域內(nèi)），獲得的圖像對之間的變換關(guān)系越精確?；诒Ａ舾嗾_匹配對的考慮，建立雙邊的匹配，即找出I1到I2之間匹配和I2到I1匹配，建立I1和I2之間的完全匹配，進一步解決了漏匹配問題。雙邊“或”的BBF算法在獲得更多正確匹配的同時，也獲得了更多的誤匹配。

為進一步去除初步建立的匹配特征對中的誤匹配，首先利用SIFT特征旋轉(zhuǎn)不變性，根據(jù)匹配特征的方向一致性初步去除誤匹配。然后通過建立對應(yīng)特征對之間的空間幾何關(guān)系進行一致性檢測，有效去除誤匹配。匹配點一致性檢測運用類似于文獻［11］中的原理，代替計算對應(yīng)線段比作為它的檢測數(shù)據(jù)。通過對應(yīng)點空間幾何關(guān)系-空間斜率和距離進行一致性檢測，去除誤匹配。它基于這樣的事實：對于所有對應(yīng)特征點對，它們的空間斜率和空間距離應(yīng)該具有一致性。

當所有點對都是物理對應(yīng)點時，對應(yīng)的空間斜率和空間距離應(yīng)是一樣的。一致性檢測首先對所有線段比值數(shù)據(jù)進行分級聚類，找出一致性數(shù)值。然后通過空間斜率和空間距離數(shù)據(jù)處理，迭代刪除偽匹配對?？臻g幾何關(guān)系的一致性檢測算法詳細描述如下：

步驟1：計算所有對應(yīng)特征點對的空間斜率p_theta和空間距離p_dis

式中，（xi，yi）和（xj，yj）分別為I1和I2中相對應(yīng)的匹配特征點對。

步驟2：對得到的兩組數(shù)據(jù)分別進行分級聚類，使聚類結(jié)果的類別數(shù)少且類內(nèi)距離足夠小。聚類結(jié)果若有某一類樣本個數(shù)大于3，且大于其它類的樣本個數(shù)，則該類樣本值稱為一致性數(shù)值。

步驟3：統(tǒng)計其余樣本的非一致性數(shù)據(jù)個數(shù)和一致性數(shù)據(jù)的方差。

步驟4：將非一致性數(shù)據(jù)個數(shù)值最大的樣本集作為候選樣本，然后將方差最大的樣本集刪除。

步驟5：再次統(tǒng)計剩余樣本的非一致性數(shù)據(jù)個數(shù)和一致性數(shù)據(jù)的方差。

步驟6：若非一致性數(shù)據(jù)的樣本集數(shù)仍大于3，且至少有一個樣本集統(tǒng)計的值大于0，則轉(zhuǎn)到步驟4；否則，輸出保留的一致性數(shù)據(jù)樣本集，兩種幾何關(guān)系一致性檢測后保留的樣本集中的數(shù)據(jù)對應(yīng)的點對即為對應(yīng)的匹配對。

圖2為特征匹配結(jié)果。

圖2 特征匹配。（a）初始匹配；（b）最終匹配Fig.2 Feature matching.（a）initial match；（b）final match

2 變換模型及參數(shù)估計

眼底結(jié)構(gòu)通常近似于剛體，而眼底表面是一個空間的曲面，擁有所有曲面的通性，即能以一定數(shù)學表達式來表達，并反映其空間形態(tài)的本質(zhì)特征［12］。文獻［13-14］中認為眼底是一個二次曲面，而二次曲面包括了橢球面、球面、拋物面及雙曲面。眼底的二次曲面模型，即多項式模型是在相機的弱透視模型和對眼底曲面結(jié)構(gòu)的近似，因此以二次曲面的數(shù)學模型近似眼底表面是合理的。點p=（x，y）T和q=（x’，y’）T分別是圖像I1和I2中的相應(yīng)匹配點，兩幅圖像的空間二次曲面模型的變化關(guān)系可表示為［2］

其中θm，n是2×6參數(shù)矩陣

M估計是一種魯棒的參數(shù)估計方法，也可以消除誤匹配帶來的外點，從而進一步提高估算精度。對θ的M估計算法具體如式（7）和式（8）所示。

式中，k表示待配準的兩圖像具有k對匹配的特征點對，ρ是具有魯棒性損失函數(shù)。進行M估計的一個重要技術(shù)就是根據(jù)實際情況合理選擇ρ。常用的三種ρ函數(shù)的性能如圖3所示［16］，圖中橫坐標為歸一化的誤差u，縱坐標為ρ的權(quán)重值。根據(jù)魯棒估計理論［15］，σ可通過式（9）計算：

圖3 不同損失函數(shù)權(quán)重值示意Fig.3 Plots of the weight functions for the robust loss functions

由圖3可見，Beaton-Tukey函數(shù)反應(yīng)最靈敏，收斂最快，可以有效去除最多的外點，故采用這個函數(shù)。

對θ的估計通過迭代加權(quán)最小二乘法（iteratively-reweighted least-squares，IRLS）實現(xiàn)。通常情況下，IRLS可以很快實現(xiàn)收斂。

在M估計中，權(quán)值函數(shù)求解方法為：

由式（10）可知，Beaton-Tukey函數(shù)對應(yīng)的權(quán)值函數(shù)為

IRLS迭代的初始矩陣取自仿射變換矩陣。

3 算法評估

源圖像由福建省附屬第一醫(yī)院眼科提供，圖像大小為3 072像素×2 048像素。目標圖像共28組，需要說明的是，所討論方法中一致性檢測的實現(xiàn)是基于同組圖像尺度變化較小的前提條件。為定量評價本方法的性能，運用正確匹配特征的數(shù)量和圖像配準精度進行評估。

匹配特征對的數(shù)量和分布，是決定參數(shù)估計準確性的重要因素。把所采用的雙邊“或”匹配方法與單邊方法及雙邊“與”方法進行評測。匹配特征的數(shù)量評測準則定義如下：令θ表示圖像I1和I2之間的映射關(guān)系，點p=（x，y）T和q=（x’，y’）T分別是圖像I1和I2中的相應(yīng)匹配點，理想情況下q=θX（p）。實際上，當p和q在空間位置上足夠接近時認為匹配關(guān)系存在，即‖q-θX（p）‖≤t，則認為p和q是一對正確匹配。實驗中選擇距離門限t=1.5。

圖像配準精度是衡量圖像配準算法性能的有效方法。這里采用均方根誤差（root of mean square error，RMSE）進行配準精度的測量：

式中，N表示最終得到的正確匹配對的數(shù)量。

在檢驗所提出方法在特征點匹配及配準精度等方面的優(yōu)勢時，對病變眼底圖像進行實驗。首先通過病變眼底圖像的配準結(jié)果，說明本方法的有效性；然后利用匹配特征的數(shù)量評估指標，對采用雙邊“或”匹配方法與單邊方法的文獻［5］和文獻［9］及雙邊“與”方法的文獻［10］進行評測；同時計算配準成功圖像對的均方根誤差評估本方法，并與文獻［5］方法所獲得的配準精度進行對比，定量描述本方法的優(yōu)勢。

4 結(jié)果

對28組實驗數(shù)據(jù)進行實驗，其中的27組實驗數(shù)據(jù)配準成功。圖4給出了兩組不同病變程度眼底圖像的配準結(jié)果示例，圖4（a）的圖像受玻璃體病變影響使圖像模糊、渾濁，尤其病灶區(qū)域血管結(jié)構(gòu)已無法提取，圖4（b）的圖像受出血灶影響也很難提取完整的血管結(jié)構(gòu)。從兩組圖像的配準結(jié)果中，可以看出配準后的圖像保留了有效信息，病變區(qū)域及模糊的血管結(jié)構(gòu)均實現(xiàn)了良好的細節(jié)對齊，無明顯的錯誤對齊現(xiàn)象，能夠滿足眼底的醫(yī)學診斷要求。

在分析正確匹配特征數(shù)量的實驗中，分析本方法對保留正確匹配對的影響。表1給出了圖4中兩組示例圖像組產(chǎn)生的匹配特征情況，同時比較文獻［5］、［9］采用單邊方法和文獻［10］采用雙邊“與”方法取得的正確匹配特征的數(shù)量。從表1可以看出，本文方法保留了足夠的正確匹配對，與雙邊“與”方法相比能夠獲得更多的正確匹配對，與單邊方法相比也有一定的優(yōu)勢。

圖4 病變眼底圖像配準結(jié)果。（a）示例1；（b）示例2Fig.4 Pathological retinal image registration results.（a）example 1；（b）example 2

表1 特征匹配方法比較Tab.1 Comparison of feature matching method

利用均方根誤差RMSE分析本研究配準成功的27組圖像的配準精度。同時與文獻［5］采用一般的特征點單邊匹配方法結(jié)合常用的RANSAC算法（RANdom SAmple Consensus，即隨機抽樣一致性算法），去除誤匹配點的基于不變特征配準方法行進對比實驗，其中，經(jīng)實驗統(tǒng)計在28組圖像序列中，本方法配準成功27組病變眼底圖像，文獻［5］配準成功25組。從圖5中可看出本方法對每組病變眼底圖像的RMSE值均小于1，并浮動于0.5左右，遠達到了亞像素級的配準精度。從兩種方法的RMSE值的對比可以看出，本方法中每組圖像數(shù)據(jù)的RMSE值均小于文獻［5］方法的RMSE值，說明本方法優(yōu)于文獻［5］方法的配準精度。

圖5 配準精度對比Fig.5 Comparison of registration accuracy

5 討論和結(jié)論

從實驗結(jié)果可見，采用雙邊“或”的特征匹配和利用空間幾何特性的一致性檢測去除誤匹配的配準方法的優(yōu)越性。同時說明在基于特征的配準方法中，準確與豐富的匹配特征在改善和提高配準性能、實現(xiàn)精確配準方面的重要性。因此，采用合適的去除誤匹配和保留足夠的匹配特征對的方法，對于基于特征的配準方法提高配準精度是有效可行的。本研究的創(chuàng)新點在于合理地運用了特征的雙邊匹配策略，在基于方向特征去除誤匹配的前提下加入了空間幾何特性的一致性檢測，有助于保留足夠的正確匹配特征對，為提高參數(shù)估計的精度奠定了基礎(chǔ)。

病變眼底圖像配準是目前醫(yī)學眼底圖像處理和分析中的一項挑戰(zhàn)性課題，其研究工作的開展和研究成果的獲得將直接影響到該學科領(lǐng)域的發(fā)展和臨床眼科診斷水平的提高。在前期研究工作的基礎(chǔ)上，提出了一種基于特征的病變眼底圖像配準方法。通過28組病變眼底圖像進行實驗仿真和分析表明，對于不同類型的病變眼底圖像該方法都具有較好的配準精度，具有一定的臨床參考價值和應(yīng)用價值。在對匹配特征進行一致性檢測時，是基于同組眼底圖像尺度變化不大的前提下進行的，如何實現(xiàn)病變眼底圖像大尺度變化下的一致性檢測及多幅眼底圖像配準是下一步研究的方向。

［1］王翠平，李新光，季亞成，等.高血壓和糖尿病患者視網(wǎng)膜病變與腦梗死的關(guān)系［J］.實用心腦肺血管病雜志，2007，15（2），152-156.

［2］Can A，Stewart C，Roysam B，et al.A feature-based，robust，hierarchical algorithm for registering pairs of images of the curved human retina［J］.IEEE Tran Pattern Ana.Mach Intel，2002，24：347-364.

［3］Stewart C，Tsai CL，Roysam B.The dual-bootstrap iterative closestpoint algorithm with application to retinal image registration［J］.IEEE Trans Med Imag，2003，22（11）：1379-1394.

［4］Zana F，Klein J.A multimodal registration algorithm of eye fundus images using vessels detection and hough transform［J］.IEEE Trans Biomed Eng，1999，18（5）：419-428.

［5］Wei LiFang，Huang LinLin，Pan，Lin，et al.The retinal image mosaic based on invariant feature and hierarchial transformation models［C］//Proceedings Image and Signal Processing，CISP.Tianjin：IEEE，2009：3463-3467.

［6］Harris C，.Stephens MJ.A combined corner and edge detector［C］//Proceedings of the 4th Alvey Vision Conference.Manchester：IEEE，1988：147-152.

［7］Bay H，Tuytelaars T，Gool LV.Surf：Speeded up robust features［C］//Proceedings of the 9th European Conference on Computer Vision.Graz：Springer，2006：404-417.

［8］Lowe D.Object Recognition from local scale-invariant features［C］//Proceedings of the International Conference on Computer Vision.Washington，DC：IEEE Computer Society，1999：1150-1157.

［9］Lowe D.Distinctive Image Features from Scale-Invariant Keypoints［J］.International Journal of Computer Vision，2004，60（2）：91-110.

［10］Chen Jian，Smith R，Tian Jie，et al.A novel registration method for retinal images based on local features［C］//Proceedings of IEEE Eng Med Biol Soc.Vancouver：IEEE，2008：2242-2245.

［11］LI H，Manjunath BS，Mitra SK.A contour-based approach to multisensor image registration［J］.IEEE Trans on Image Processing，1995，4（3）：320-334.

［12］邵婷婷，鄭穗聯(lián)，王波，等.正常個體人眼角膜空間形態(tài)數(shù)學建模路線研究［J］.眼視光學雜志，2005，7（4）：253-256.

［13］Burek H，Douthwaite WA.Mathematical models of the general corneal Surface［J］.Ophthal Physiol Opt，1993，13（1）：68-72.

［14］Richard L，George S，et al.Descriptors of corneal shape［J］.Op2tom Vis Sci，1998，75（2）：156-158.

［15］Rousseeuw PJ，Leroy AM.Robustregression and outlier detection［M］.New York：John Wiley ＆ Sons，1987：216-247.

［16］Stewart C.Robust parameter estimation in computer vision［J］.Society for Industrial and Applied Mathematics.1999，41（3）：513-537.