基于双目视觉的视网膜血管三维重建

周嘉悦，韩少峰，郑昱，吴壮志，梁庆丰，杨洋

1 北京航空航天大学 机械工程及自动化学院，北京市，100191

2 北京航空航天大学 计算机科学与技术学院，北京市，100191

3 北京同仁医院，北京市，100062

0 引言

眼底是眼球内后节组织，包括视网膜、视乳头、黄斑和视网膜中央动静脉，眼底如图1所示。人眼球的横径约为24 mm，眼球纤维膜的后5/6部分为巩膜。视网膜血管以视乳头为中心，呈网状分布在视网膜上。通过选取样本进行放大测量可以得知视网膜血管的动脉直径约为193.0～203.7 μm，静脉血管直径约为249.5～263.4 μm[1]。普通的眼底影像技术如眼底照相机无法获得眼底参数的三维信息，影响疾病的诊断、治疗和手术定位。

图1 眼底示意图Fig.1 Fundus diagram

手术操作点定位是眼科显微手术机器人研究中的重要研究方向[2]，在视网膜血管搭桥[3]、视网膜血管光凝以及视网膜血管注药等手术中，都需要定位手术器械刺入目标位置（视网膜血管根部、视网膜血管堵塞处、视网膜血管出血点等），并以此为反馈控制机器人运动。目前，手术器械定位依赖医生的肉眼观察。考虑到人类医生的尺度感知能力有限、微尺度的视网膜血管难以精准定位，肉眼观察定位的方式易造成手术误操作，引起器官损伤。因此，需要精确地获取视网膜血管上的刺入点和血管的脉络走向，为医生或机器人提供精准的手术操作点位置信息。双目视觉技术，通过识别目标点的三维坐标等信息，可对目标点进行精确定位，具有应用于眼科手术中的可能。因此，该文中采用双目视觉进行视网膜血管三维重建。

目前，双目立体视觉的研究多侧重于宏观尺度物体的定位和重建，对小尺度物体重建的研究以及在医学方面的应用相对较少。天津大学的高礼圳等[4]使用“距离空间图”匹配算法，在3幅图像上建立匹配关系，误差为1～3 mm。长春理工大学的冯进良等[5]使用双目视觉系统监控焊点位置，精度可达0.68 mm。日本研究人员开发出能在手术过程中透视体内淋巴结、血管等组织的手术“导航”系统，在手术中可借助彩色图像清晰地观察到上述组织的状态[6]。上海理工大学的喻海中等[7]利用双目视觉原理获取了医疗器械的3D数据，与实际测量出的数据相比精度更高，可用于模拟手术和手术导航。

本研究首先搭建双目视觉系统，并对摄像机进行标定，之后获取视网膜血管模型的双目图像并校正为共面平行的两平面图像。对图像进行滤波提取分割出视网膜血管，去除背景，减小误匹配产生的误差，最后利用双目视觉中的视差原理得到视网膜血管的点云图，并进行误差分析。

1 问题描述

为验证双目视觉测量出点的坐标的准确性，采用Solidworks绘制眼球模型作为测量参考，通过对Solidworks中测量出目标点的三维坐标与双目视觉测量出的三维坐标进行对比，验证三维模型的准确性。在Solidworks中绘制的眼球模型，直径38 mm，深度26 mm，内部刻有宽度和深度均为1 mm的凹槽用来模拟视网膜血管。视网膜血管模型，如图2所示。

图2 Solidworks中绘制的视网膜血管模型Fig.2 Retinal vessel model in Solidworks

双目立体成像原理如图3所示，基线距离B表示两摄像机的投影中心连线的距离，摄像机在同一时间观察同一点P，分别获取了左图像和右图像上的点P的坐标Pl=(Xl,Yl)；Pr=(Xr,Yr)。经过立体校正后的图像位于同一平面上，故Yl=Yr=Y。视差D=Xl-Xr，由几何关系可得点P在摄像机坐标系下的三维坐标如式(1)～(3)所示：

图3 双目立体视觉三维重建原理Fig.3 Principle of 3D-reconstruction

2 眼底血管的双目重建

2.1 双目视觉系统搭建

实验中使用的相机为大恒工业USB3.0工业相机MER-500-14U-3C，分辨率为2 592×1 944。双目视觉系统如图4所示。

图4 双目视觉系统Fig.4 Binocular vision system

摄像机镜头选用computar镜头，镜头焦距为12～36 mm，最大成像尺寸8.8 mm×6.6 mm，芯片尺寸H×V=4.8 mm×3.6 mm。通过式(4)可以计算出镜头的实际焦距为13.7 mm。

式(4)中，f表示镜头焦距，W表示物距，F表示视场大小，H×V表示芯片尺寸。对于双目视觉系统，改变参数(如基线长度、物距等)可以减小测量误差，提高系统的精度，双目视觉系统在x、y、z三个方向的理论分辨率如式(5)～(7)所示。

其中，h表示物距，f表示焦距，b表示相机光心距离又称基线距离，△d表示像素精度。双目视觉系统在x、y、z方向上的分辨率与物距和焦距有关，在z方向上的分辨率还与基线的长度有关。由于模型尺寸较小，所以在进行双目立体视觉平台搭建时，需要尽可能保证精度。考虑到工作距离和拍摄视野的要求，尽量减小物距，增大基线距离。

为了满足摄像机公共视野的要求，选择基线距离为68 mm，物距约为500 mm，摄像机精度为0.001 85 mm，根据式(5)～(7)可以求出该双目视觉系统在x、y方向上的理论分辨率为0.06 mm，在z方向的理论分辨率为0.49 mm。

2.2 摄像机标定

本研究采用“张正友标定法”[8]对摄像机进行标定，成像模型如式(8)所示。

式(8)中，k代表比例因子，K为摄像机坐标系的内参矩阵，[X Y Z1]T为模板平面上点的齐次坐标，[u v1]T为模板平面上点投影到图像平面上对应点的齐次坐标，如图5所示。R和T分别是摄像机坐标系相对世界坐标系的旋转矩阵和平移矢量。

图5 平面标定算法原理图Fig.5 Zhang's calibration algorithm

由于摄像机内参矩阵中有5个未知参数，所以当所取得的图像数目大于等于3时，可线性唯一求解出K。

相机标定时采用棋盘格靶标，精度为0.001 mm，棋盘格为边长2.6 mm的正方形。使用双目摄像机对靶标拍摄图片各15张。将拍摄的图片导入Matlab中进行标定，得到标定结果如表1和表2所示。通过标定得到重投影误差为0.5个像素，约为0.03 mm。

2.3 立体校正

图6中P点为空间中某一点的位置，Ol，Or分别为左右摄像机光心，Pl，Pr分别为P点在左右图像上的成像点。gl，gr分别为左右极线。需要通过立体校正将实际中共面行不对准的两幅图像如图6(a)，校正为共面行对准的两幅图像如图6(b)。

表1 相机内部参数Tab.1 Camera internal references

表2 左右相机相对位姿Tab.2 Relative pose of two cameras

图6 立体校正过程Fig.6 Rectify process

理想的双目系统中两摄像机图像平面平行，光轴和图像平面垂直，极点处于无穷远处。要满足这个条件，必须使两幅图像平面和两个摄像机坐标原点的连线平行。计算Rrect矩阵如式(9)所示，使极点处于无穷远处。

式(10)中，T是右摄像机相对于左摄像机的平移向量。e1和e2正交，求出e2如式(11)所示。

e3和e1、e2正交，如式(12)所示。

可以解出矩阵Rrect。将Rrect左乘到R分解后作用于左右相机坐标系的矩阵，即可得到最终的立体校正矩阵。

根据Solidworks绘制的视网膜模型，通过3D打印将其打印为实物，拍摄实物图片，利用上文标定出的数据对其进行立体校正，对该模型拍摄的两张图片如图7所示，能够清楚地看到血管的分布。

图7 视网膜血管模型图片Fig.7 Retinal vessel model images

校正后得到共面平行的两张图片如图8所示。

图8 视网膜血管模型校正后图像Fig.8 Rectified retinal vessel model image

2.4 图像分割

与传统的直接进行立体匹配的方法[9]不同，本研究采用先将血管进行分割，再进行匹配的方法，消除了周围背景区域对匹配的影响，只留下视网膜血管进行后续的立体匹配和三维重建，使重建后的视差图更加清晰完整。

采用对比度受限的自适应直方图均衡算法（CLAHE）[10]对图像进行处理，这种方法采用直方图匹配方法来逐个处理图像中的较小区域，然后使用双线性内插方法将相邻的小片组合起来，从而消除引入的边界。在灰度均匀的区域，通过限制对比度来避免放大噪声。

对CLAHE后的图像进行均值滤波，然后求取均值滤波后的图像与原图的差值，再进行二值化，分割出视网膜血管轮廓。

分割出来的视网膜血管图像是二值的，用其作为掩膜应用于视网膜血管灰度图像中，即可分割出去除背景的视网膜血管灰度分布图，如图9所示。

图9 视网膜血管分割过程Fig.9 Retinal vessel segment

2.5 双目立体匹配

立体匹配的目的是找到左右图像的对应点，需要在一定的范围内进行搜索。由于本研究中只需要对视网膜血管上的点进行匹配，所以采用基于局部的立体匹配算法能够更好地对图像上的感兴趣区域进行匹配。参考图像的坐标(u，v)和视差搜索范围d就构成了视差图。视差图中的每一点(u，v，d)的像素表示的是位于参考图像上的点(u，v)和位于匹配图像上的点(u+d，v)之间的“距离”。SAD算法[11]是一种基于局部的图像匹配算法，常用于图像块匹配，原理如式(13)所示。

计算机实验教学中心建设要考虑安全化和节约化。中心建设要牢固树立“安全第一，预防为主”和“安全无小事”的理念，坚决克服麻痹懈怠思想，认真做好安全防范工作，明确安全责任人，确保安全。实验中心建设要做好五防：防火、防水、防雷、防盗、防病。

SAD算法的基本流程如下：

（1）构造一个类似卷积核的小窗口；

（2）用窗口覆盖左边图像，选出覆盖区域的像素点；

（3）用窗口覆盖右边图像并选出覆盖区域像素点；

（4）求出左右区域像素点差的绝对值之和；

（5）移动右边图像上的窗口，重复(3)、(4)，直到找出差值最小的窗口，即找到了与左图匹配的像素块。

直接对原图像进行匹配，背景区域会对视网膜血管上的点的匹配造成影响，因此本研究先将视网膜血管分割，再进行匹配，两种匹配方式得到的视差图，如图10所示。

图10 视网膜血管视差图Fig.10 Disparity map of retinal vessel

从图10(b)中可以清晰地看出血管的分布和走向，说明先对血管进行分割去除背景后再进行匹配能够达到更好的匹配效果。

3 结果与误差分析

采用第1节中的双目视觉三维重建原理可以求出视网膜血管上每点的三维坐标，在Matlab中绘制出视网膜血管的三维点云图像如图11所示。

图11 三维重建出的点云图像Fig.11 Point cloud of 3D reconstruction

任取视网膜模型图像上的7个血管交叉点来对实验结果准确性进行验证，如图12所示。

图12 7个目标点图像Fig.12 7 target points in image

由于测量目标点实际的三维坐标较为困难，文中将通过每个点的相对位置来对误差进行估计。

提取出每个点的像素坐标，并通过双目视觉三维重建算法即式(1)～(3)计算出每个点的三维坐标，将其称为点集A，其中每个点的三维坐标如表3所示。将模型在Solidworks坐标系中的三维坐标称为点集B，其中每个点的三维坐标如表4所示。寻找点集B相对点集A的旋转矩阵R和平移矩阵T，使得点集B旋转平移后与点集A之间的误差最小。求解出中心点，将点集中心化后计算点集之间的协方差矩阵，通过奇异值分解(SVD)可求解出R、T。

式(14)、(15)中，μA、μB是两个点集的中心点，H是两个点集之间的协方差矩阵。最终得到R、T如式(16)、(17)所示。

将点集B旋转平移前后与点集A的相对位置如图13所示。

计算出每两个对应点之间的误差和7个点的平均误差如表5所示。

表3 所取7个目标点通过双目视觉系统求出的三维坐标Tab.3 The coordinates of seven selected target points obtained by binocular vision system

表4 所取7个目标点通过Solidworks测量出的三维坐标Tab.4 The coordinates of seven selected target points measured with Solidworks

表5 两个点集之间的误差Tab.5 The error between 2 sets of points

从表5中可以看出，平均误差为0.68 mm，产生误差的原因是相机精度有限、相机标定过程中产生的重投影误差和拍照时光线和角度的影响以及3D打印的误差。可以通过选择精度更高，视野更广的相机、增加标定精度和改善图片质量来进一步提高测量精度。

图13 两坐标系下的点集配准Fig.13 Point registration between two coordinates

4 结论

本研究采用视网膜模型以及双目视觉系统对眼底视网膜血管进行图像处理，通过先分割后匹配的方法得到完整的视差图以进行三维重建，得到眼底视网膜血管点云图，对于血管上的任意一点，能较为精确地求出其三维坐标，控制误差在1 mm之内，基本满足对视网膜血管上任意一点进行定位的要求。本方法可对小尺度的视网膜血管模型进行三维重建，并获得定位点的坐标，后续的研究中，将使用真实眼底视网膜血管进行实验，从而实现对真实眼底视网膜血管的三维重建。