高钦泉，黄伟萍，杜 民，韦孟宇,3，柯栋忠

(1.福州大学 物理与信息工程学院，福州 350116； 2.福州大学 福建省医疗器械与医药技术重点实验室，福州 350116；3.澳门大学 模拟与混合信号超大规模集成电路国家重点实验室，澳门 999078)

0 引言

微创手术具有创伤小、疼痛轻、恢复快的优越性，在临床中受到广泛的欢迎[1]。通过内窥镜等监控系统，医生可以通过小孔在人体内施行手术，对人体伤害小，达到理想的手术效果。在传统的盆腔微创手术导航中，医生通过从术前影像(如CT(Computed Tomography)、MRI(Magnetic Resonance Imaging)等)来确定病灶的位置并规划手术实施方案，通过内窥镜从微小的孔洞获取即时的手术视觉信息，引导手术的推进。医生观察解剖结构和手术器械的操作过程都是通过内窥镜摄像头传输至显示器来进行，然而，传统的内窥镜通常只能为医生提供二维的脏器表面的视觉信息，医生在手术过程中很难能够准确地掌握器官内或空腔脏器内部解剖结构的三维深度信息，尤其对于早期肿瘤的定位就更加困难，这对于准确地操作器械提出了更大的挑战和要求。在手术过程中，医生通常能看到的只是暴露在内窥镜视觉下的器官表面，要开展精确的手术，医生往往需要回顾并参考复杂的人体解剖学结构，结合大脑中记忆的术前拟定的手术方案路径，并借助内窥镜捕获的实时画面来寻找定位手术的目标，从而引导手术的开展。在盆腔内实施手术过程中，内窥镜只提供狭窄的二维视觉信息，由于人体内部器官组织错综复杂，很难根据当前二维视觉信息获取当前手术刀在人体当中准确的位置，有可能造成手术刀“迷路”。如果不能准确找到病灶部位，容易造成癌变组织切除过少或者正常组织切除过多，造成肿瘤残留或者对器官功能的过多破坏，极大降低手术效果。因此，临床医生对术中借助先验知识的视觉导航有着急切的需求。目前，已经不同形式的手术导航技术应用于手术当中。大体可以分成两类：1)基于标记点的手术导航，如基于生物标记点的光学跟踪[2]、基于人工标记的手术导航仪[3]，但这类技术过度依赖于标记点；2)基于无标记点的手术导航，如电磁导航跟踪[4]，但其设备较为复杂及昂贵，术中X射线透视[5]、CT成像辅助导航[6]等但实时影像会给病人甚至医生带来高强度的辐射。目前还没有一套较简便的设备能应用于复杂的盆腔微创手术,实现精准视觉导航的系统。为此，本文提出并设计了一种基于立体视觉的盆腔微创手术增强现实(Augmented Reality, AR)导航仿真系统，充分利用术前CT影像重建病人的三维组织信息，并通过增强现实(AR)技术渲染显示在手术屏幕上，为医生提供“透视”的视觉导航，为临床的手术应用提供技术支持。

在仿真手术开始进行时，首先手动调整术前重建的模型使之与拍摄的患者仿真手术视频中相对应的骨盆目标区域大致对齐，然后利用基于可视点颜色一致性的2D/3D配准的配准技术进行微调，实现立体内窥镜下2D手术视野与3D骨盆模型的初始化配准；再利用双目立体视觉跟踪算法对立体内窥镜进行实时跟踪，获取内窥镜的运动信息；最后，根据内窥镜多自由度的变换矩阵来实现术前3D模型与手术视野的实时融合和增强现实显示，将术前3D模型与手术视野画面准确融合并成像于手术显示屏，从而给医生视觉辅助，以定位当前手术刀的位置以及到手术目标的方向和距离，为手术的实施提供良好的视觉导航作用。

1 AR视觉导航仿真系统设计

1.1 AR导航仿真系统框架设计

整个盆腔微创手术AR导航系统的具体实施方法如下：

1)首先进行双目内窥镜的标定[7]，以获得内窥镜的内参数矩阵K，畸变矩阵Q，左右内窥镜的相对姿态。

2)其次，根据病人术前CT影像，利用ITK-SNAP系统分割出骨盆模型，依据第一步标定后双目内窥镜的参数，使用开放图形库(Open Graphics Library, OpenGL)渲染在本导航系统显示界面当中。

3)然后,先手动调整模型，使之与目标区域大致对齐，再利用本文配准算法进行微调，达到更精确的配准。待配准对齐后，将骨盆模型的表面顶点覆盖在内窥镜2D图像的目标区域，达到增强现实的视觉融合。

4)最后，利用双目视觉的左右视图信息来实现双目内窥镜的移动轨迹的实时跟踪，根据内窥镜前后帧移动的多自由度变换矩阵，对术前3D模型进行位置更新，与内窥镜视野下的仿真手术画面进行实时融合显示在导航系统界面，达到增强现实的视觉导航效果。

1.2 盆腔手术仿真视频的制作

由于真实手术场景下的环境错综复杂，且内窥镜真实轨迹很难获取，无法直接在手术中进行算法验证，本文提出一种盆腔手术虚拟双目内窥镜下仿真视频的制备方法，可以获取具有真实纹理的手术场景的仿真视频和内窥镜的真实轨迹。

本文截取骨盆肿瘤的微创手术中真实场景画面作为骨盆3D模型的纹理，使用纹理坐标轴UV展开对3D模型进行纹理贴图，形成具有真实手术纹理信息的骨盆3D模型。模型依据虚拟内窥镜的内参数进行投影变换到2D平面，形成视频帧。先预设虚拟内窥镜的轨迹，作为内窥镜的真实轨迹；再根据模型的模型矩阵与视图矩阵的关系，对骨盆模型作旋转平移变换，与此同时，录制当前渲染画面。通过仿真虚拟内窥镜的移动轨迹，利用该模型制作一套具有真实轨迹的仿真骨盆手术视频。

1.3 基于颜色一致性的双目视觉2D/3D初始化配准

实现术前3D模型与真实手术视觉场景中的目标物体完全融合在一起，是实现增强现实至关重要的步骤。而真实手术视觉通常都是2D图像，需要将病人术前重建的3D模型与术中手术的2D视野进行配准。本文提出了一种在双目内窥镜视觉下的3D模型颜色一致性的配准方法。通过优化双目内窥镜的姿态，使得双目视觉左右视图下3D模型的可视点的颜色距离达到最小，来实现3D模型到2D手术视野配准的目的。

1.3.1 基于颜色一致性的内窥镜姿态初始化配准

在理想状态下，同一个顶点从不同角度或不同时刻观测的色彩空间不变，即具有颜色一致性[8]。因此，如图1所示，在双目内窥镜的左右相机姿态下各自姿态观测骨盆3D模型，当顶点经投影映射到二维坐标系统中所对应的左右视图二维手术画面具有一致的色彩空间时，3D模型能够与目标物体配准对齐。基于上述的假设，本文提出的3D模型与真实物体2D视觉初始化融合配准的目标能量方程可定义为:

(1)

图1 基于立体视觉的颜色一致性

然而，在实际场景中，光照及物体表面的反射率等是影响颜色一致性的重要因素。因此，本文采用自然环境光为光源，以降低光照条件对颜色一致的影响，同时，以左右视图颜色RGB的平均值作为参考值。在为使术前3D模型和真实2D手术视觉能够配准对齐，设3D模型在左相机坐标系下的视图矩阵为Tcm，式(1)的能量方程可进一步表示为：

(2)

当初始化配准对齐后，3D模型与真实物体2D视觉完全融合，此时内窥镜的初始位置可根据当前模型视图矩阵获得，3D模型姿态Tcm与真实相机姿态Tcamera的关系为：Tcamera=(Tcm)-1。

在真实场景中，由于光照角度的影响，导致同一个顶点从不同角度观察可能会不一致，但是，由于本文使用的是双目内窥镜，其左右镜头离得较近且内参数较为相近，受光照角度影响较小，同时观察一个顶点的颜色信息相差不大，而且式(2)旨在寻找全部可视点的颜色误差全局最优。

1.3.2 3D模型的可视点快速检测法

双目立体视觉的颜色一致性是基于3D模型的可视点进行计算和优化的。常见的可视点寻找法有Z缓冲器算法等[9]，但时间复杂度较高。本文提出一种基于双视图颜色模型的可视点快速寻找法。术前重建的3D模型是由多个三角片组成，对每个三角片进行编号，按序号将三角片渲染成不同颜色，形成一个带颜色的3D模型，通过OpenGL的3D图形渲染引擎投影几何关系，生成2D投影图像Iproject，如图2所示，Iproject是由W×H个像素点组成的，并且每个像素点的颜色RGB是已知的。

图2 颜色模型

根据OpenGL的投影变换原理，求出3D模型的顶点Pi在屏幕上面的坐标值(u,v)T：

(3)

如图3(a)所示，三角片的三个顶点分别按投影变换关系投影到二维平面，若投影区域包含一个或多个整型像素点，输出第一个最先被搜寻到的整型像素点坐标值(u,v)T。图3(b)展示搜寻算法，步骤如下。

输入：二维平面投影点p1:(u1,v1),p2:(u2,v2),p3:(u3,v3)。

输出：p:(u,v)。

步骤1 分别计算L1，L2，L3直线的一般方程aix+biy+ci=0(i=1,2,3)。其中:a1=v2-v1，b1=u1-u2，c1=-(a1u1+b1v1)；a2=v2-v3，b2=u3-u2，c2=-(a2u2+b2v2)；a3=v3-v1，b3=u1-u3，c3=-(a3u3+b3v3)。

步骤2 计算z1=a2u1+b2v1+c2；z2=a3u2+b3v2+c3；z3=a1u3+b1v3+c1。

步骤3 计算xmin=min(ui)，xmax=max(ui)，ymin=min(vi)，ymax=max(vi)，i=1,2,3。

步骤4i=ceil(xmin)，j=ceil(ymin)。判断a2i+b2j+c2与z1,a3i+b3j+c3与z2，a1i+b1j+c1与z3是否同正负号。若是，则(i,j)在投影区域里，输出(u=i,v=j)T，算法结束;否则, 继续下一步。

步骤5i++，j++，如果i≤xmax,j≤ymax转步骤4；否则(u,v)T由最近邻插值得到。

搜寻算法得出(u,v)T，判断在Iproject对应颜色RGB与其颜色模型上的三角片颜色是否一致。若是一致的，则认为三角片的三个顶点在当前姿态为可视的；反之不是。

模型表面的三角片投影区域受视点远近影响，当视点较远时，此时三角片投影区域没有包含一个整型像素点，此时，用传统的插值方法(如双线性插值(Bilinear Interpolation, BI)、最近邻插值(Nearest Neighbor, NN)、部分体积插值(Partial Volume, PV))得到投影点的RGB(本文使用的是BI)，再作判断该三角片是否可视。

图3 三角片可视点判断

基于颜色一致性的配准是基于双目立体内窥镜共同的可视点，所以可视点需在左右内窥镜都能被观测到。如图4所示，图4(a)可视点能够同时被左右内窥镜观测，图4(b)中不可视点是由于只被左内窥镜观测。

图4 模型可视点检测

1.4 基于双目视觉的内窥镜运动轨迹跟踪算法

目前有很多基于视觉的跟踪算法被提出，主要分成两类: 基于特征点的跟踪，如PTAM(Parallel Tracking And Mapping)[10]、ORB-SLAM(Oriented FAST and Rotated BRIEF Simultaneous Localization And Mapping)[11-12]等; 基于像素点的跟踪，如DTAM(Dense Tracking And Mapping)[13]、SVO(Semi-direct monocular Visual Odometry)[14-15]等。无论是基于特征点还是像素点，提出相对应的”点”都是为了相机定位准备。相机定位问题属于状态估计问题，主流使用的方法有两种: 基于滤波的方法，如MonoSLAM(Mono Simultaneous Localization And Mapping)[16]，但定位精度不高；基于优化的方法，这也是近几年常用的方法。基于优化的方法分成两个阶段: 构图和跟踪[13],常用集束调整(Bundle Adjustment, BA)[17]。由于有向的Fast算子和可旋转的BIREF描述子(Oriented FAST and Rotated BRIEF, ORB)匹配速度快，且具有较稳定的旋转不变性，因此本文选用ORB作为特征点提取，使用立体BA来优化相机位姿。

(4)

(5)

其中，e(k,Xj)=ρ(xj-π(RkXj+tk))，Ωk为滑动窗口中全部的三维点。

图5 局部地图跟踪

2 实验验证与结果分析

实验环境如下，操作系统为ubuntu 16.04的台式机，中央处理器(Central Processing Unit, CPU)为Intel Core i7-4790 CPU 3.60 GHz×8,随机存取存储器(Random Access Memory, RAM)为8 GB，图形处理器(Graphics Processing Unit, GPU)为NVIDIA GeForce GTX 960，显存为2 GB。开发环境为Qt 5.9.1，并使用C++以及OpenGL实现。

2.1 骨盆的3D模型及纹理模型制备

仿真视频采用的双目相机参数来自伦敦帝国理工学院Hamlyn Centre的双目内窥镜参数：1)http://hamlyn.doc.ic.ac.uk/vision/data/Dataset9/Left_Camera_Calibration_Intrinsic.txt; 2)http://hamlyn.doc.ic.ac.uk/vision/data/Dataset9/Right_Camera_Calibration_Intrinsic.txt; 3)http://hamlyn.doc.ic.ac.uk/vision/data/Dataset9/camera_extrinsic.txt。如图6所示，采用CT影像在ITK-SNAP重建的3D模型，使用UV展开对空白的骨盆模型进行纹理贴图。纹理图选自盆腔手术视频(http://pan.baidu.com/s/1eSBy2X0)，图7的盆腔手术视频截图与图6(b)的表面纹理具有较高的相似度，说明纹理模型能够为仿真手术提供具有类似真实纹理的信息。

图6 骨盆3D模型

图7 手术视频截图

2.2 基于颜色一致性的双目视觉2D/3D初始化配准

在初始配准前，需对模型进行可视点检测，如图8所示，模型顶点总共有340 074个，设置成不同颜色，然后按左右内窥镜的相机参数渲染模型，检测到左右视图共同的可视点有21 474个，耗时62 ms。

将配准融合完全的模型分别绕X、Y、Z轴从偏离-30°旋转到偏离30°，每次旋转0.1°；沿X、Y、Z轴从-20 mm平移到20 mm，每次平移0.1 mm。然后按式(2)计算出可视点在左右视图投影点颜色的均方差，如图9(a)和图9(b)所示。两幅图可以说明，存在模型视图矩阵Tcm使得式(2)达到全局最优的状态。

图8 可视点检测效果

图9 顶点在左右视图的投影点均方差

式(2)的优化是非线性优化。梯度下降法如Newton、Gauss-Newton和Levenberg-Marquardt等非线性优化算法虽然在非线性优化有良好的效果，但是经实验测试，式(2)在使用梯度下降法优化的过程中很容易陷入局部最优状态，故不适合使用梯度下降法进行优化。本文使用无导数优化(derivative-free optimization)中的二次逼近约束优化(Bound Optimization BY Quadratic Approximation, BOBYQA)[18]进行优化, 而且这是一种有约束的无导数优化。将优化变量变成6维的李代数形式：Tcm=eξ∧。Tcm为4×4矩阵，但ξ为6维向量，前三维为旋转向量，后三维为平移向量，优化的变量转化为6个维度。BOBYQA优化将置信区间作为步长，然后沿着每一维度变化，更新二次模型，最后根据二次模型求出最优解。由图9可以得出在0值(ground truth)投影点均方差并不为零，这是由于部分顶点投影颜色是由插值得到，导致左右视图顶点投影颜色无法精确相等，观察图9定义可视点的全部可视点的颜色均方差小于5，则认为3D模型与手术2D画面配准融合完全。配准效果如图10(b)。

图10 配准效果

2.3 基于双目视觉的内窥镜运动轨迹跟踪算法

本文使用g2o框架[19]的Levenberg-Marquardt梯度下降法对相机姿态进行优化。采用此算法在仿真视频上进行跟踪，共700帧，跟踪轨迹与真实轨迹显示如图11所示，实线表示真实轨迹，虚线是本文跟踪算法的轨迹，两条轨迹非常接近，均方根误差为2.393 3 mm，说明跟踪轨迹能够提供精度较高的位置信息。

图12、图13、图14分别表示轨迹图的Z-X视图、Y-X视图、Z-Y视图，从中可以看出，Z-X视图的两条轨迹基本吻合；Y-X视图、Z-Y视图的最大偏差约为5 mm。进一步验证了该算法估计的位置信息精度高。

图11 轨迹三维显示

图12 Z-X 视图

图13 Y-X视图

为了进一步验证跟踪算法的效果，将重建的3D模型利用3D打印机打印出实体模型，并涂上血肉颜色，形成仿真的打印模型。模拟手术的过程，拍摄一段手术视频，在该系统上进行测试。跟踪一段时间后，效果显示如图15所示。从视觉观测，重建的模型与打印模型的融合程度较高，说明具有良好的AR显示效果。

图14 Z-Y视图

图15 AR显示

3 结语

本文设计开发了一套基于立体视觉的盆腔微创手术增强现实(AR)导航仿真系统，在本系统中提出了一种基于立体视觉的2D/3D初始化配准融合，并且实现了基于ORB的立体视觉跟踪算法来估计相机的轨迹。立体内窥镜观测骨盆时，利用左右内窥镜的相对位姿来实现颜色一致性，以完成骨盆模型的初始化配准。通常术前图像(CT或MRI)等模态影像有较高的分辨率和视觉效果，待重建的3D模型与2D画面完全融合，能够很好地为医生提供AR显示，能更快、更准确地确定病灶的位置。基于立体视觉的跟踪算法，在不借助其他硬件设备时，能够为医生提供较准确的位置信息。仿真实验结果显示，AR盆腔手术导航系统具有良好的初始化配准、融合，能够较准确地跟踪内窥镜，为临床医生实施手术提供视觉辅助。

本文的仿真视频制备能够为不易获取真实盆腔手术轨迹的场景提供一种验证算法准确性的途径，在其他手术场景也同样适用。在未来的工作中，将结合实际临床手术来开展深入验证工作，并进一步改进真实手术场景中双目内窥镜的视觉跟踪算法，根据实际临床数据设计AR手术导航的显示方法，为基于双目视觉的手术导航技术研究奠定基础。

[19] KÜMMERLE R, GRISETTI G, STRASDAT H, et al. g2o: a general framework for graph optimization [EB/OL]. [2017- 12- 10]. http://www.willowgarage.com/sites/default/files/g2o-icra11.pdf.