基于单目视觉的视网膜三维重建技术研究

王玉亮　刘　飞　沈建新　刘贤喜

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年生。山东农业大学机械与电子工程学院教授、博士研究生导师。