基于單目視覺的視網膜三維重建技術研究

王玉亮　劉　飛　沈建新　劉賢喜

1.山東農業大學，泰安，271018　　2.南京航空航天大學，南京，210016

?

基于單目視覺的視網膜三維重建技術研究

王玉亮1,2劉飛1沈建新2劉賢喜1

1.山東農業大學，泰安，2710182.南京航空航天大學，南京，210016

為了在視網膜疾病診斷中獲得精度高、可視化好的視網膜，提出一種基于單目視覺的視網膜三維重建算法。該算法通過對兩幅預處理圖像的特征點提取確定匹配點對之間的對應關系，采用RANSAC算法去除誤匹配點對，準確率高，用鏈匹配的方法使得多幅圖像間的特征點匹配從而得到整體最優化結果。根據不同圖像間特征點的對應關系，使用4通光束平差法確定空間中攝像機相對位置關系，使用PMVS算法實現視網膜的三維重建。實驗結果表明，該算法具有很好的魯棒性和穩定性，能較好地實現視網膜的三維重建。

視網膜圖像；尺度不變特征變換；鏈匹配；4通光束平差法；三維重建

0　引言

通過視網膜可以非創傷性地直接觀察到較深層次的微血管，常見一些疾病如高血壓、糖尿病等疾病可以在視網膜上被觀察到[1]，因此視網膜的研究對一些疾病的早期發現、診斷和治療有著重要作用。視網膜三維重建可以獲得精度和可視化效果好的視網膜，使得三維信息的完整性得以保證，從而幫助醫生對視網膜疾病作出準確判斷，改變以往單純依靠經驗的方式。

Choe等[2]提出了基于熒光素圖像序列進行視網膜三維重建和圖像配準的方法，該方法需要注射熒光素進行檢測，對眼睛有損害。Lin等[3]提出一種基于多模態視網膜圖像進行二維和三維配準的方法，該方法也是基于熒光素圖像進行圖像的配準，該方法同樣對眼睛有損害。Lusk等[4]提出一種通過相機標定和圖像矯正進行自動三維重建立體視網膜圖像的方法，在方法中，變形算法已被一個更直接的相機幾何估計算法取代，轉化成理想的立體相機的幾何形狀，從而大大簡化了視差圖生成和現有算法的自動分析階段，但導致了重建精度的下降。

針對目前視網膜三維重建方法的不足，本文提出基于單目視覺的視網膜三維重建技術，主要包括圖像獲取、圖像預處理、圖像特征提取與匹配、稠密重建，其中匹配和稠密重建是視網膜三維重建的關鍵技術。

1　多模態視網膜圖像匹配

本文匹配方法包括三個主要步驟：對圖像的預處理，增強視網膜圖像的邊緣信息對比度；然后本文對預處理的圖像進行SIFT特征提取，采用迭代最近鄰域匹配算法進行圖像對間的匹配，采用隨機抽樣一致性(random sample consensus, RANSAC)算法提純；最后采用鏈匹配的方法進一步提高匹配結果[5]。

1.1視網膜圖像的預處理

為了獲取大量的SIFT特征點，本文采用增強圖像的對比度進行圖像的預處理。采用限制對比度自適應均衡化(CLAHE)的方法[6]來增強圖像邊緣信息對比度，如圖1所示。

(a)原彩圖1　　　　(b)圖1a的CLAHE結果

(c)原彩圖2　　　　(d)圖1c的CLAHE結果圖1　視網膜圖像CLAHE增強前后對比

1.2視網膜圖像間的特征提取及匹配

采用SIFT特征提取點的方法進行特征點的提取[7]。預處理后的視網膜圖像獲得更多的匹配點對，誤匹配點對多，無法獲得圖像間精確的單應性關系。本文采用迭代最近鄰域匹配的算法來解決這個問題[8]。首先對圖像Ii和Ij的單應性關系進行估計，得到單應性矩陣Hi,j，對于在圖像Ii上的特征點x，可以在以Hi,j(x)為中心的一個小的鄰域內采用迭代最近鄰匹配，能得到好的匹配結果。隨著特征匹配對數量的增加，可以得到精確的單應性關系。

本文的方法表述如下：首先，選擇一個較大的搜索區域r=rmax并進行單應性矩陣Hi,j初始估計。每進行一次迭代，屬于圖像Ii的每個特征點x與以Hi,j(x)為搜索中心、r為半徑的區域內的圖像Ij特征點進行最近鄰域匹配。然后，采用RANSAC算法進行精確的單應性估計用于下一次迭代。最后，在搜索區域的最小半徑到達設定值時停止迭代[9]。通過這種方法可以獲得更多可信賴的匹配。

1.3鏈匹配

采用鏈匹配的方法使圖像間的匹配得到整體最優化。如I2不能直接和I3匹配，但I3和I1可以直接匹配，I1和I2可以直接匹配，可以以I1為媒介進行I2和I3的最近鄰匹配，初始化H2,3=H1,3H2,1、r=rmin，在匹配點為中心的粗略小區域內直接執行最近鄰匹配。這樣可以得到足夠的匹配點對，從而獲得精確的單應性矩陣。在進行鏈匹配時，用到了價值函數。選擇好的圖像進行鏈匹配，用最短路徑算法進行計算，使Ii到Ij達到最短路徑的圖像用來進行鏈匹配[10]。價值函數被定義為內點數倒數，如下式所示：

Ci,j=1/inlier(Hi,j)

(1)

其中，函數inlier(·)用來計算內點數量，價值函數Ci,j的數值越大說明匹配過程得到的單應性內點越少，可以通過如上例所描述初始化H2,3=H1,3H2,1、r=rmin進行I2和I3的最短路徑匹配，相應的內點數也會增加，以此類推，當內點不再增加便停止匹配。所有圖像通過最短路徑再次匹配，單應性關系和內點得到優化。采用沿著最短路徑的連鎖單應性作為最佳的單應性估計。最后，選取相對其他所有圖像匹配最小代價的圖像作為參考圖像Iref，表達式如下所示：

(2)

2　視網膜的三維重建

由于視網膜的表面是近似平面的，這種特殊的結構使得視網膜的三維重建工作和其他物體的三維重建方式有所區別[10]。

2.1對極幾何估計

在三維重建過程中，首先要估計對極幾何，經典的7點、8點算法不能很好地對近平面結構進行基礎矩陣估計。由于取樣視場的限制，對應的采樣點通常不能夠提供有效的約束計算正確的基礎矩陣，因此，在進行基礎矩陣估計前驗證這些對應采樣點的有效性就顯得很有必要。

本文提出一種平面+視差的方法進行基礎矩陣F的估計，基礎矩陣通過下式求出：

F=[e′]×H

(3)

(a)平面+視差法估計基礎矩陣

其中，e′是極點，[e′]為極點e′坐標的矩陣形式，H為單應性矩陣，在本方法中引入一個平面π，如圖2a所示。由于本文已對圖像對間的單應性矩陣進行了估計，選擇兩個不在同一平面內的對應點操作，如圖2b所示，確定極點e′，最后可以得到精確的基礎矩陣F。通過這種方法，可以計算出每幅圖像Ii和參考圖像Iref之間的基礎矩陣Fi。

(b)定位核點圖2　平面+視差法估計基礎矩陣和定位核點示意圖

根據下式求出本質矩陣：

E=KTFK

(4)

其中，K為相機內參，由相機標定得出。然后對E進行奇異值分解得到外參矩陣R和t。最后通過下式獲得投影矩陣P[11]：

P=K[R|t]

(5)

通過式(5)可以確定投影矩陣，獲得的基礎矩陣和相應的內點將被下一步的4通光束平差法用來進行相機位置的估計。

2.24通光束平差法估計相機位置

二維匹配以后，根據所有匹配對中最短路徑選擇參考框架Iref。用平面+視差法估計圖像序列Ii相對于參考圖像Iref的基礎矩陣Fi，對應的基礎矩陣內點表示為τi，然后代入到4通光束平差法精確地估計相機的位置[3]。

在光束平差的第一步，首先假定視網膜結構Xi在平面上并平行于參考相機，采用最小重投影誤差的方法估計相機的位置：

(6)

(7)

如圖3c所示。Pi和Xi已經在前兩步求出。最后選擇相對于參考圖像Iref最寬基線的圖像為Ibest估計視網膜稀疏結構，具體步驟如圖3d所示。

(a)估計相機位置 　　　　(b)估計三維結構

(c)精確相機位置和三維結構(d)確定稀疏結構圖3　4通光束平差法示意圖

通過光束平差調整后，稀疏結構和相機位置得到整體的優化。獲得所有攝像機位置以后，需要建立稠密對應，獲得視網膜的稠密點云[12]。

2.3視網膜稠密點云獲取及三維重建

通過建立稠密對應獲得視網膜稠密對應結構，采用基于面片的多視圖立體重構PMVS算法，將獲得的稀疏點云生成視網膜稠密點云。在上一步的光束平差法已經進行攝像機矩陣和稀疏結構的計算，然后采用PMVS算法獲得稠密點云。在此算法中，擴散和濾波是迭代進行的，本文進行3次循環[13]，獲得視網膜圖像的空間點云數據后，采用泊松表面重建的方法進行重建，最后進行紋理映射得到三維重建后的視網膜。

3　實驗結果及分析

系統平臺實現：本文是在硬件平臺為Intel(R) Core(TM) i3 CPU 3.20 GHz、內存為4G，軟件為MATLAB R2012b(64位)條件下實現的。每組采用9幅圖像，像素為1600×1200，迭代最小搜索半徑設為25個像素，匹配耗時125.23 s，獲得匹配點對為7921對，視網膜重建總耗時約42 min。

SIFT算法的閾值選擇對提取特征點的數量、匹配及視網膜三維重建的精度、耗時有直接影響。圖4a所示為匹配點對隨閾值的變化情況，閾值越大匹配點對越多，誤匹配點對越多，選取0.7為閾值，圖4b為不同方法的匹配率對比，本文使用迭代最近鄰提純和鏈匹配方法匹配正確率較高，能獲得好的匹配精度，圖4c為選取的視網膜圖像，圖4d為兩幅視網膜圖像匹配結果，初匹配為284對，提純后為273對，匹配正確率為96.12%。

(a)匹配點對隨閾值的變化

(b)不同方法的匹配率對比

(c)9幅視網膜圖像

(d)兩幅視網膜圖像匹配結果圖4　特征提取與匹配對比

(a)文獻[10]方法重建的病變視網膜圖像

(b)本文方法重建的病變視網膜圖像

(c)文獻[10]方法重建低質量視網膜圖像

(d)本文方法重建的低質量視網膜圖像

(e)文獻[10]方法重建的一般視網膜圖像

(f)本文方法重建的一般視網膜圖像

(g)文獻[10]方法重建的高質量視網膜圖像

(h)本文方法重建的高質量視網膜圖像圖5　視網膜三維重建結果對比

圖6　平均局部深度方差對比

本文采用50組不同的視網膜圖像進行重建，其中1～10組是病變圖像、11～20組是低質量圖像、21～30組是一般圖像、31～50組是高質量圖像。圖5所示分別為病變圖像、低質量圖像、一般圖像和高質量圖像4組圖像，采用文獻[10]重建方法與本文重建方法進行三維重建結果對比，可以得出本文重建方法可以對不同條件下的視網膜得到較好的重建效果。本文對50組圖像通過對每個點云附近100像素范圍內求局部深度方差，最后得到平均局部深度方差如圖6所示，可以看出本文所用方法重建出的視網膜結構更加平滑符合視網膜近平面的結構特性；圖7所示為點云數據與重建表面距離平均偏差對比，偏差越小表示重建誤差越小，由圖7可知，本文重建精度較高。本文重建方法操作簡單，對眼睛沒有損害，較好地實現了視網膜的三維重建。

圖7　距離平均偏差對比

4　結語

針對視網膜的特點，對目前三維重建方法進行了改進，提出采用限制對比度自適應均衡化的方法進行視網膜圖像的預處理，提出迭代最近鄰匹配方法進行特征匹配，采用4通光束平差法進行攝像機位置和視網膜稀疏結構估計，最后通過基于面片的多視圖立體重構PMVS算法的方法進行視網膜的三維重建。結果表明，本文算法具有很好的重建效果。

[1]鄧可欣. 基于超邊圖匹配的視網膜眼底圖像配準算法[J]. 清華大學學報(自然科學版), 2014，54(5):568-574.

DengKexin.RetinalImageRegistrationBasedonHyper-edgeGraphMatching[J].TsinghuaUniv.(Sci. &Technol.), 2014，54(5):568-574.

[2]ChoeTE,MedioniG,CohenI,etal. 2-DRegistrationand3-DShapeInferenceoftheRetinalFundusfromFluoresceinImages[J].MedicalImageAnalysis， 2008, 12(2):174-190.

[3]LinY,MedioniG.RetinalImageRegistrationfrom2Dto3D[C]//IEEEConferenceonComputerVisionandPatternRecognition.Anchorage,AK, 2008:1-8.

[4]LuskJA,NutterB.Automated3-DReconstructionofStereoFundusImagesviaCameraCalibrationandImageRectification[C]//22ndIEEEInternationalSymposiumonComputer-basedMedicalSystems(CBMS).NewYork:IEEE, 2009:1-7.

[5]王玉亮，沈建新，廖文和，等. 基于CUDA的眼底圖像快速自動配準與拼接[J]. 中國機械工程, 2013(13):1749-1754.

WangYuliang,ShenJianxin,LiaoWenhe,etal.FastAutomaticFundusImagesRegistrationandMosaicBasedonCUDA[J].ChinaMechanicalEngineering, 2013(13):1749-1754.

[6]陳萌夢，熊興良，張琰，等. 一種視網膜眼底圖像增強的新方法[J]. 重慶醫科大學學報, 2014, 39(8):1087-1090.

ChenMengmeng，XiongXingliang，ZhangYan,etal.ANewMethodforRetinalFundusImageEnhancement[J].JournalofChongqingMedicalUniversity, 2014, 39(8):1087-1090.

[7]LoweDG.DistinctiveImageFeaturesfromScale-invariantKey-points[J].InternatinalJournalofComputerVision, 2004, 60(2):91-110.

[8]朱利成，姚明海. 基于SIFT算法的目標匹配和識別[J]. 機電工程, 2009, 26(12):74-81.

ZhuLicheng,YaoMinghai.ObjectMatchingandRecognitionBasedonSIFTAlgorithm[J].Mechanical&ElectricalEngineeringMagazine, 2009, 26(12):74-81.

[9]趙燁，蔣建國，洪日昌. 基于RANSAC的SIFT匹配優化[J]. 光電工程, 2014(8):58-65.

ZhaoYe,JiangJianguo,HongRichang.AnOptimizedSIFTMatchingBasedonRANSAC[J].Opto-electronicEngineering, 2014(8):58-65.

[10]TaeEC,MedioniG. 3-DMetricReconstructionandRegistrationofImagesofNear-planarSurfaces[C]//11thIEEEInternationalConferenceonComputerVision.NewYork:IEEE, 2007:1895-1902.

[11]楊敏，沈春林. 基于對極幾何約束的景象匹配研究[J]. 南京航空航天大學學報, 2004, 36(2):235-239.

YangMin,ShenChunlin.StudyonSceneMatchingBasedonEpipolarGeometricConstraint[J].JournalofNanjingUniversityofAeronautics&Astronautics, 2004, 36(2):235-239.

[12]JinH,SoattoS,YezziA.Multi-viewStereoReconstructionofDenseShapeandComplexAppearance[J].InternationalJournalofComputerVision, 2005, 63(3):175-189.

[13]史利民，郭復勝，胡占義. 利用空間幾何信息的改進PMVS算法[J]. 自動化學報,2011,37(5):560-568.

ShiLimin,GuoFusheng,HuZhanyi.AnImprovedPMVSthroughSceneGeometricInformation[J].ActaAutomaticaSinica, 2011,37(5):560-568.

(編輯王艷麗)

Research on 3D Reconstruction of Retina Based on Monocular Vision

Wang Yuliang1,2Liu Fei1Shen Jianxin2Liu Xianxi1

1.Shandong Agricultural University,Taian,Shandong,271018 2.Nanjing University of Aeronautics and Astronautics,Nanjing,210016

In order to obtain high accuracy and good visualization in the diagnosis of retinal diseases, a new algorithm for 3D reconstruction was proposed based on monocular vision. Through two pieces of preprocessing image feature points extraction the algorithm was used to determine the matching relationship between pairs of points and RANSAC algorithm was adopted to remove the false matching points to obtain the feature point to improve the accuracy, the feature points matching of multiple images were obtained by using the method of chain registration. According to the corresponding relationship among the feature points in different images, the relative position of the camera in space was determined by using the method of 4 pass bundle adjustment, PMVS algorithm was used to achieve three-dimensional reconstruction of the retina. Experimental results show that the algorithm has good robustness and stability, may better realize the three-dimensional reconstruction of the retina.

retinal image; scale-invariant feature transform(SIFT); chain registration; 4 pass bundle adjustment; three-dimensional reconstruction

2015-10-08

江蘇省科技支撐計劃資助項目(BE2010652)

TP391.41

10.3969/j.issn.1004-132X.2016.18.011

王玉亮，男，1979年生。山東農業大學機械與電子工程學院講師、博士。主要研究方向為機器視覺、醫學圖像處理。發表論文10余篇。劉飛，男，1988年生。山東農業大學機械與電子工程學院碩士研究生。沈建新，男，1969年生。南京航空航天大學機電學院教授、博士研究生導師。劉賢喜，男，1963年生。山東農業大學機械與電子工程學院教授、博士研究生導師。