胡 超,蔡振宇,尤晓赫,张智焕

(1.浙江大学宁波理工学院信息科学与工程学院,浙江 宁波 315100;2.浙江大学控制学院,杭州 310007)



基于胶囊内窥镜图像的肠胃道三维重建技术*

胡 超1*,蔡振宇1,尤晓赫2,张智焕1

(1.浙江大学宁波理工学院信息科学与工程学院,浙江 宁波 315100;2.浙江大学控制学院,杭州 310007)

为了改进胶囊内窥镜观测的准确性和真实性,提出了基于胶囊内窥镜序列图像的胃肠道三维重建的方法。首先利用SIFT算法提取前后两幅序列图像中尽可能多的对应特征点;计算获取各特征点在成像面上的二维坐标;进一步利用8点算法计算胶囊内镜运动变化的旋转矩阵和平移矢量。进而计算得到每个特征点的相对三维坐标和世界三维坐标;然后,采用Delaunay三角剖分算法对各三维点进行网格化,并完成场景的三维重建。实验表明相机与被测点距离在100 mm之内时,得到的深度误差小于1 mm;距离250 mm内时,相对误差在3%之内。说明所提出的算法是可行的。

胶囊内窥镜;三维重建;胃肠道图像;8点位姿算法

无线胶囊内窥镜(以下称胶囊内镜)M2A(图1)于2000年在以色列Given Image公司首先开发成功[1-2],次年被美国FDA通过并应用于临床。后来,改名为PillCam。通过吞服,它能拍摄食道、胃、小肠和大肠图像,完成全消化道系统的观测。由于能实现完整小肠(5 m～7 m)的检查,且避免了给病人带来大的痛苦,其发展潜力明显,是内窥镜检查术的一个重大突破。但是,从临床应用中反馈来看,现有胶囊内镜有很多问题需要解决。其中有病灶组织的尺寸确定问题:胶囊内窥镜图像是二维的,观测时不够直接,对病变组织的判断不够准确。为此,有必要对胶囊内镜图像进行病灶组织的三维测量。

图1 M2A胶囊内窥镜(现名PillCam)

这种三维重建可利用胶囊内窥镜拍摄的序列图像和相应的胶囊位置和姿态定位信息来完成。常用的定位方法有磁定位技术[3-4],能够提供胶囊的3维位置和2维方向信息,但是不能给出胶囊绕主轴旋转的信息。为此,我们提出用胶囊内窥镜图像序列进行定位,获取胶囊内窥镜的3维旋转和位移方向参数。进一步结合磁定位的结果实现3维位置和3维姿态定位和胃肠道病灶组织的三维重建。显然,这种基于胶囊内窥镜图像定位是单目相机定位,前期已有很多的研究成果[5-7];但是,这些成果在在胶囊内窥镜中应用,来进行胃肠道和病灶组织三维重建尚未有报道。

图2 三维重建流程图

我们提出的肠胃道病灶组织三维重建流程如图2所示。首先对胶囊内镜的前后图像进行分析,并SIFT算法提取两幅图像中相对应的多个特征点[8]。基于特征点在成像面上的坐标,能用8点算法计算分解旋转矩阵和平移矢量[9],即得到相机运动的位移方向和旋转角度的变化。进而计算得到每个特征点的相对三维坐标和世界三维坐标。然后,对各三维点进行网格化场景三维重建[10-11]。本文将对这些步骤和方法进行阐述。

1 特征点提取

特征点是指胶囊内镜在不同位姿下形成的成对图像中,对同一目标点产生的成像点。目前,对于特征点提取,常采用SIFT算法对匹配得到。SIFT算法采用不同尺度的高斯函数对原始二维图像滤波,形成高斯金字塔,然后通过相邻两幅图像的差值来形成查分金字塔图像。通过极值点周围的像素灰度信息生成特征向量。如图3所示,两幅配对图像中,利用SIFT,在左图找到481个特征点,右图443个特征点,而匹配得到197个特征点。从实验知,当图像内容较丰富时,可以找到比较多的匹配点,但是当图像分辨率低,图像比较平滑时,只有很少的成功匹配的特征点。

图3 SIFT算法匹配结果

2 求取胶囊内镜运动参数(八点算法)

在得到前后图像的匹配点后,可求取胶囊的前后位置和方向的变化(即运动参数)。如图4所示,胶囊内窥镜的成像光心从O转移到Or。P=(x,y,z)T为运动前在XYZ坐标系上的空间目标点坐标;P′=(x′,y′,z′)T为胶囊运动后,该点在新的成像坐标系X′Y′Z′上的坐标。

图4 胶囊内镜(单相机)立体成像原理

设(u,v)和 (u′,v′) 分别代表P和P′的在成像面上坐标。当z=f=1时,两个坐标点为 (u=x/z,v=y/z)和(u′=x′/z′,v′=y′/z′),即有

P′=R0·P+T0

(1)

式中:R0为3×3的旋转矩阵,T0=[tx0ty0tz0]T是平移矢量。

应用八点算法,取任意非零并同T0共线的矢量T(T×T0=0),在式(1)两侧做向量积,

T×(x′y′z′)T=T×[R0(xyz)T]+T×T0

(2)

即有

(3)

由于T×(u′v′ 1)T和(u′v′ 1)是正交的,在式(3)两侧点乘以(u′v′ 1)T,得到

(u′v′ 1)(T×R0)(uv1)T=0

(4)

式中:T×R0=[T×r1T×r2T×r3],r1,r2,r3为R0的列向量。运动参数矩阵E定义为

E=T×R0

(5)

这个E矩阵定义为基本矩阵。式(4)表明前后两幅图像上的对应点(u,v)和(u′,v′),与矩阵E满足齐次线性的:

(u′v′ 1)E(uv1)T=0

(6)

设

(7)

从式(6)和式(7),得到

[u′uu′vu′v′uv′vv′uv1]·

(8)

由式(8)得到

UΛ=0

(9)

||UΛ||2=min

(10)

下面,把基本矩阵E分解为平移向量T和旋转矩阵R0。

(11)

假定

(12)

这样向量T的解是EET最小特征值对应的特征向量。从式(3)有z′T×(u′v′1)T=zT×[R0(uv1)T]=zE(uv1)T

(13)

那么

∑[T×(u′v′ 1)T]T·E(uv1)T>0

(14)

由于z<0,z′<0,T×(u′v′ 1)T和E·(uv1)T有相同的方向,并且都是非零的,所以如果不满足式(14),那么T←(-T)。

得到T以后,那么旋转矩阵R0就可以从E=T×R0确定。如果没有噪声,可以得到

R0=[W2×W3,W3×W1,W1×W2]-T×E

(15)

式中:W1=T×E1,W2=T×E2,W3=T×E3,E1、E2、E3分别为矩阵E的列向量。

如果存在噪声,那么R0满足

‖E-[T]xR0‖2=min

(16)

(17)

问题就转变为

(18)

这个方程如同‖R0C-D‖=min,其中C=[C1C2C3],D=[D1D2D3]。这样就可以通过四元数组来得到旋转矩阵R0。

定义一个4×4的矩阵F

(19)

式中:

(20)

定义一个四元数组q=(q0q1q2q3)T为矩阵F最小特征值对应的特征向量,就得到‖R0C-D‖2=min中的结果R0为

(21)

3 三角视觉确定目标三维坐标

如图4所示,取左边相机的光心为世界坐标系的原点,那么世界坐标系为o-xyz,左图坐标为Ol-XlYl,右图坐标为Or-XrYr,三维空间中的一点P(x,y,z)在左成像平面中的投影为Ol(Xl,Yl),在右成像平面中的投影为Or(Xr,Yr),根据透视模型得到:

(22)

(23)

(24)

(25)

三维世界坐标系o-xyz和or-xryrzr可以由下式表示:

(26)

式中:M为空间转换矩阵

(27)

R与T为两个测量点之间的旋转矩阵与平移矢量。通过式(22)～(27)的联立计算得到每个特征点在真实场景中的三维坐标,进而也可以进行三维重建。

(28)

4 离散空间点的三角网格构建

由于特征点是三维空间中离散的点,要采用三角剖分算法把三维点集剖分为若干个三角网格,得到三角模型的三维几何外形。在众多算法中,Delaunay三角剖分算法和模型数据结构简单,适应多种状况,得到了大量应用。Voronoi图是Delaunay三角剖分的基础。Voronoi图可以直观的定义为从每个点开始,以一样的速率向外扩展直到遇到另一个点的扩展范围,如图5中虚线形成的就是Voronoi图。可以看到Voronoi图为凸的多边形。Delaunay三角网格由Voronoi图得到,Voronoi图中每条边表示Delaunay三角网中每个点之间的对应的连线。图4中实线形成的就是Delaunay三角网格。

图5 Voronoi图与Delaunay三角网

在得到图像特征点三维坐标,并采用Delaunay三角剖分算法对所有数据点结成三角网格后,就可以在三维空间中绘制出目标的三维拓扑结构。为了增强三维显示效果,将图像中的纹理使用OpenGL纹理贴图技术贴到物体模型的表面,这样就可以使重建出的三维模型具有较强的真实感。

5 实验结果与分析

利用典型的胶囊内镜的图像,设计了仿制的胃肠道。利用单目相机对实际构成的场景进行拍摄,得到图6所示的两幅模拟胶囊内窥镜在不同位姿下拍摄的图像。

图6 待三维重建的两幅图像

三维重建的步骤如下:

①特征点的提取。得到52对特征点,如图7所示。

图7 特征匹配结果

②相机相对运动的计算。利用特征点匹配的结果,通过八点算法计算得到相机的运动参数平移矢量T和旋转矩阵R0。

T0=[0.999 3 -0.025 1 -0.027 3]T

由于八点算法得到的相机平移矢量T不具有长度信息,本实验中结合磁定位的结果[4]得到相机具有长度的平移参数为:

T=[28.979 7 -0.727 9 0.791 7]T

③三维坐标的计算。采用式(28)对52个特征点计算得到其三维坐标,部分坐标点见表1,其中Xl、Yl、Xr、Yr分别是该特征点在左右两幅图中的坐标值,X、Y、Z分别为计算得到的该点在三维空间中的坐标点。

④三维网格的剖分。采用Delaunay三角剖分算法对三维点集进行网格化。得到的网格结果如图8上图所示,对网格进行着色后如图8下图所示,可以从图中清楚的看出场景的相关结构。

表1 三维坐标结果

(5)OpenGL纹理映射实现三维重建。

利用C++和OpenGL通过纹理映射技术显示在屏幕上。图9展示了三维重建结果在4个视角下的重建结果,可以看到,三维重建的模型比较直观,与原场景较为相似。

为了进一步验证三维重建的准确度,在实验中对若干个特征点进行指定,之后测量这些特征点与相机平面之间的距离,再与计算出的结果进行比对。本实验中,使用内窥镜相机在(0,0,0)mm和(30,0,5)mm,相机主轴角度相差20°的前后两点上进行图像的拍摄;并选定了6个特殊点进行匹配与测量,实验结果如表2所示。

图8 三维坐标网格

图9 三维重建结果

特征点目标实际深度/mm计算得到深度/mm误差/mm相对误差/%14747.68210.68211.4526464.86480.86481.3538080.89890.89891.124114116.12902.12901.865168171.42863.42862.046241248.27597.27593.01

由表2可以看到,当特征点离相机镜头平面的距离在100 mm之内时,算法计算得到的深度与目标实际深度之间的误差较小,小于1 mm,系统精度较高。相机与被测特征点的距离大于100 mm时,由于深度距离变远,误差也相应增大。总体来说,在深度250 mm内,相对误差控制在3%之内。由于胶囊内窥镜图像的拍摄距离较近,这样深度已经足够。且对于胶囊内窥镜临床应用,医生建议的目标测量精度为10%内,本方法的精度指标能够达到这一目标。

6 结论

针对胶囊内窥镜的序列图像,提出了胃肠道三维重建的方法。首先提取前后两幅序列图像中相对应的多个特征点;基于特征点计算各特征点在成像面上的坐标;进一步利用8点算法计算胶囊内镜运动变化的旋转矩阵和平移矢量,即得到其运动的位移方向和旋转角度变化。进而计算得到每个特征点的相对三维坐标和世界三维坐标;然后,采用Delaunay三角剖分算法对各三维点进行网格化,并完成场景的三维重建。实验结果表明所提出的算法是可行的。

[1] Gavriel I,Gavriel M,Arkady G,et al. Wireless CaTPule Endoscopy[J]. Nature,2000,405(6 785):417.

[2] Max Q H M,Mei T,Pu J X,et al. Wireless Robotic CaTPule Endoscopy:State-of-the-Art and Challenges[C]//Proceedings of the 5th World Congress on Intelligent Control and Automation,2004:5561-5565.

[3] 任宇鹏,胡超,项圣,等. 一种人体移动对胶囊内窥镜磁定位干扰的补偿方法[J]. 传感技术学报,2015,28(11):1640-1646.

[4] Hu C,Li M,Song S,et al. A Cubic 3-Axis Magnetic Sensor Array for Wirelessly Tracking Magnet Position and Orientation[J]. IEEE Sensors Journal,2010,10(5):903-913.

[5] Shen Y,Guturu P,Buckles B. Wireless CaTPule Endoscopy Video Segmentation Using an unsupervised Learning Approach based on Probabilistic Patent Semantic Analysis with Scale Invariant Features[J]. IEEE Transactions on Information Technology in Biomedicine,2012,16(1):98-105.

[6] Hyon L,Sudipta N S,Michael F C,et al. Real-Time Monocular Image-Based 6-DoF Localization[J]. The International Journal of Robotics Research,2015,34(4-5):476-492.

[7] Seo-Yeeon H,Jae-Bok S. Monocular Vision-Based SLAM Indoor Environment Using Corner,Lamp,and Door Features from Upward-Looking Camera[J]. IEEE Transactions on Industrial Electronics,2011,58(10):4904-4812.

[8] 佟帅,徐晓刚,易成俦,等. 基于视觉的三维重建技术综述。计算机应用研究[J]. 2011,28(7):2411-2417.

[9] Hartley R I. In Defense of the Eight-Point Algorithm[J]. IEEE Trans on Pattern Analysis and Machine Intelligence,1997,19(6):580-593.

[10] 罗肖. 序列内窥镜图像的三维结构重建[D]. 上海:上海交通大学,2009.

[11] Sun B,Liu L,Hu C,et al. 3D Reconstruction Based on CaTPule Endoscopy Image Sequences[C]//International Conference on Audio Language and Image Processing,2010:607-612.

3D Recovery of Gastrointestinal Tract Based on CaTPule Endoscopic Images*

HU Chao1*,CAI Zhenyu1,YOU Xiaohe2,ZHANG Zhihuan1

(1.Information Science and Engineering College,Ningbo Institute of Technology,Zhejiang university,Ningbo Zhejiang 315100,China;2.College of Control,Zhejiang university,Hangzhou 310007,China)

In order to improve the accuracy and authenticity for the caTPule endoscopic examination,we present the 3D recovery technique for gastrointestinal tract based on the sequential endoscopic images. First,the corresponding feature points are extracted from the two sequential images by SIFT method,and their 2D coordinates on the imaging plane are acquired. Then the position and orientation shifts are calculated by the 8-point algorithm,which are related to the rotation matrix and transition vector. Further,the 3D relative and global coordinates are computed for all the feature points. Finally,the 3D feature points are gridded by Delaunay triangular subdivision algorithm and the 3D recovery is realized for the real scene. The experimental results show that the depth error is smaller than 1mm in case that the distance is 100 mm from camera to the tested point,and smaller than 3 mm the distance 250 mm. These results validate the feasibility of the proposed method.

caTPule endoscopy;3D recovery;gastrointestinal tract images;8-point algorithm

胡 超(1960-),男,加拿大阿尔伯特大学博士,浙江大学宁波理工学院三江学者特聘教授;浙江大学控制学院博士生导师。从事自动化、机器人控制和定位跟踪技术研究。主持国家863计划、支撑计划、国家自然基金等科研项目30多项;发表论文200多篇。

项目来源:国家自然科学基金项目(61273332);与宁波市科技项目(创新团队);计划项目(2014B82015)

2016-09-13 修改日期:2016-12-26

TP212.9;TP242

A

1004-1699(2017)05-0708-07

C:7230

10.3969/j.issn.1004-1699.2017.05.013