基于膠囊內(nèi)窺鏡圖像的腸胃道三維重建技術(shù)*

胡 超,蔡振宇,尤曉赫,張智煥

(1.浙江大學(xué)寧波理工學(xué)院信息科學(xué)與工程學(xué)院,浙江 寧波 315100;2.浙江大學(xué)控制學(xué)院,杭州 310007)

?

基于膠囊內(nèi)窺鏡圖像的腸胃道三維重建技術(shù)*

胡 超1*,蔡振宇1,尤曉赫2,張智煥1

(1.浙江大學(xué)寧波理工學(xué)院信息科學(xué)與工程學(xué)院,浙江 寧波 315100;2.浙江大學(xué)控制學(xué)院,杭州 310007)

為了改進(jìn)膠囊內(nèi)窺鏡觀測的準(zhǔn)確性和真實(shí)性,提出了基于膠囊內(nèi)窺鏡序列圖像的胃腸道三維重建的方法。首先利用SIFT算法提取前后兩幅序列圖像中盡可能多的對應(yīng)特征點(diǎn);計(jì)算獲取各特征點(diǎn)在成像面上的二維坐標(biāo);進(jìn)一步利用8點(diǎn)算法計(jì)算膠囊內(nèi)鏡運(yùn)動(dòng)變化的旋轉(zhuǎn)矩陣和平移矢量。進(jìn)而計(jì)算得到每個(gè)特征點(diǎn)的相對三維坐標(biāo)和世界三維坐標(biāo);然后,采用Delaunay三角剖分算法對各三維點(diǎn)進(jìn)行網(wǎng)格化,并完成場景的三維重建。實(shí)驗(yàn)表明相機(jī)與被測點(diǎn)距離在100 mm之內(nèi)時(shí),得到的深度誤差小于1 mm;距離250 mm內(nèi)時(shí),相對誤差在3%之內(nèi)。說明所提出的算法是可行的。

膠囊內(nèi)窺鏡;三維重建;胃腸道圖像;8點(diǎn)位姿算法

無線膠囊內(nèi)窺鏡(以下稱膠囊內(nèi)鏡)M2A(圖1)于2000年在以色列Given Image公司首先開發(fā)成功[1-2],次年被美國FDA通過并應(yīng)用于臨床。后來,改名為PillCam。通過吞服,它能拍攝食道、胃、小腸和大腸圖像,完成全消化道系統(tǒng)的觀測。由于能實(shí)現(xiàn)完整小腸(5 m～7 m)的檢查,且避免了給病人帶來大的痛苦,其發(fā)展?jié)摿γ黠@,是內(nèi)窺鏡檢查術(shù)的一個(gè)重大突破。但是,從臨床應(yīng)用中反饋來看,現(xiàn)有膠囊內(nèi)鏡有很多問題需要解決。其中有病灶組織的尺寸確定問題:膠囊內(nèi)窺鏡圖像是二維的,觀測時(shí)不夠直接,對病變組織的判斷不夠準(zhǔn)確。為此,有必要對膠囊內(nèi)鏡圖像進(jìn)行病灶組織的三維測量。

圖1 M2A膠囊內(nèi)窺鏡(現(xiàn)名PillCam)

這種三維重建可利用膠囊內(nèi)窺鏡拍攝的序列圖像和相應(yīng)的膠囊位置和姿態(tài)定位信息來完成。常用的定位方法有磁定位技術(shù)[3-4],能夠提供膠囊的3維位置和2維方向信息,但是不能給出膠囊繞主軸旋轉(zhuǎn)的信息。為此,我們提出用膠囊內(nèi)窺鏡圖像序列進(jìn)行定位,獲取膠囊內(nèi)窺鏡的3維旋轉(zhuǎn)和位移方向參數(shù)。進(jìn)一步結(jié)合磁定位的結(jié)果實(shí)現(xiàn)3維位置和3維姿態(tài)定位和胃腸道病灶組織的三維重建。顯然,這種基于膠囊內(nèi)窺鏡圖像定位是單目相機(jī)定位,前期已有很多的研究成果[5-7];但是,這些成果在在膠囊內(nèi)窺鏡中應(yīng)用,來進(jìn)行胃腸道和病灶組織三維重建尚未有報(bào)道。

圖2 三維重建流程圖

我們提出的腸胃道病灶組織三維重建流程如圖2所示。首先對膠囊內(nèi)鏡的前后圖像進(jìn)行分析,并SIFT算法提取兩幅圖像中相對應(yīng)的多個(gè)特征點(diǎn)[8]。基于特征點(diǎn)在成像面上的坐標(biāo),能用8點(diǎn)算法計(jì)算分解旋轉(zhuǎn)矩陣和平移矢量[9],即得到相機(jī)運(yùn)動(dòng)的位移方向和旋轉(zhuǎn)角度的變化。進(jìn)而計(jì)算得到每個(gè)特征點(diǎn)的相對三維坐標(biāo)和世界三維坐標(biāo)。然后,對各三維點(diǎn)進(jìn)行網(wǎng)格化場景三維重建[10-11]。本文將對這些步驟和方法進(jìn)行闡述。

1 特征點(diǎn)提取

特征點(diǎn)是指膠囊內(nèi)鏡在不同位姿下形成的成對圖像中,對同一目標(biāo)點(diǎn)產(chǎn)生的成像點(diǎn)。目前,對于特征點(diǎn)提取,常采用SIFT算法對匹配得到。SIFT算法采用不同尺度的高斯函數(shù)對原始二維圖像濾波,形成高斯金字塔,然后通過相鄰兩幅圖像的差值來形成查分金字塔圖像。通過極值點(diǎn)周圍的像素灰度信息生成特征向量。如圖3所示,兩幅配對圖像中,利用SIFT,在左圖找到481個(gè)特征點(diǎn),右圖443個(gè)特征點(diǎn),而匹配得到197個(gè)特征點(diǎn)。從實(shí)驗(yàn)知,當(dāng)圖像內(nèi)容較豐富時(shí),可以找到比較多的匹配點(diǎn),但是當(dāng)圖像分辨率低,圖像比較平滑時(shí),只有很少的成功匹配的特征點(diǎn)。

圖3 SIFT算法匹配結(jié)果

2 求取膠囊內(nèi)鏡運(yùn)動(dòng)參數(shù)(八點(diǎn)算法)

在得到前后圖像的匹配點(diǎn)后,可求取膠囊的前后位置和方向的變化(即運(yùn)動(dòng)參數(shù))。如圖4所示,膠囊內(nèi)窺鏡的成像光心從O轉(zhuǎn)移到Or。P=(x,y,z)T為運(yùn)動(dòng)前在XYZ坐標(biāo)系上的空間目標(biāo)點(diǎn)坐標(biāo);P′=(x′,y′,z′)T為膠囊運(yùn)動(dòng)后,該點(diǎn)在新的成像坐標(biāo)系X′Y′Z′上的坐標(biāo)。

圖4 膠囊內(nèi)鏡(單相機(jī))立體成像原理

設(shè)(u,v)和 (u′,v′) 分別代表P和P′的在成像面上坐標(biāo)。當(dāng)z=f=1時(shí),兩個(gè)坐標(biāo)點(diǎn)為 (u=x/z,v=y/z)和(u′=x′/z′,v′=y′/z′),即有

P′=R0·P+T0

(1)

式中:R0為3×3的旋轉(zhuǎn)矩陣,T0=[tx0ty0tz0]T是平移矢量。

應(yīng)用八點(diǎn)算法,取任意非零并同T0共線的矢量T(T×T0=0),在式(1)兩側(cè)做向量積,

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

(2)

即有

(3)

由于T×(u′v′ 1)T和(u′v′ 1)是正交的,在式(3)兩側(cè)點(diǎn)乘以(u′v′ 1)T,得到

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

(4)

式中:T×R0=[T×r1T×r2T×r3],r1,r2,r3為R0的列向量。運(yùn)動(dòng)參數(shù)矩陣E定義為

E=T×R0

(5)

這個(gè)E矩陣定義為基本矩陣。式(4)表明前后兩幅圖像上的對應(yīng)點(diǎn)(u,v)和(u′,v′),與矩陣E滿足齊次線性的:

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

(6)

設(shè)

(7)

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

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

(8)

由式(8)得到

UΛ=0

(9)

||UΛ||2=min

(10)

下面,把基本矩陣E分解為平移向量T和旋轉(zhuǎn)矩陣R0。

(11)

假定

(12)

這樣向量T的解是EET最小特征值對應(yīng)的特征向量。從式(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以后,那么旋轉(zhuǎn)矩陣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)

問題就轉(zhuǎn)變?yōu)?/p>

(18)

這個(gè)方程如同‖R0C-D‖=min,其中C=[C1C2C3],D=[D1D2D3]。這樣就可以通過四元數(shù)組來得到旋轉(zhuǎn)矩陣R0。

定義一個(gè)4×4的矩陣F

(19)

式中:

(20)

定義一個(gè)四元數(shù)組q=(q0q1q2q3)T為矩陣F最小特征值對應(yīng)的特征向量,就得到‖R0C-D‖2=min中的結(jié)果R0為

(21)

3 三角視覺確定目標(biāo)三維坐標(biāo)

如圖4所示,取左邊相機(jī)的光心為世界坐標(biāo)系的原點(diǎn),那么世界坐標(biāo)系為o-xyz,左圖坐標(biāo)為Ol-XlYl,右圖坐標(biāo)為Or-XrYr,三維空間中的一點(diǎn)P(x,y,z)在左成像平面中的投影為Ol(Xl,Yl),在右成像平面中的投影為Or(Xr,Yr),根據(jù)透視模型得到:

(22)

(23)

(24)

(25)

三維世界坐標(biāo)系o-xyz和or-xryrzr可以由下式表示:

(26)

式中:M為空間轉(zhuǎn)換矩陣

(27)

R與T為兩個(gè)測量點(diǎn)之間的旋轉(zhuǎn)矩陣與平移矢量。通過式(22)～(27)的聯(lián)立計(jì)算得到每個(gè)特征點(diǎn)在真實(shí)場景中的三維坐標(biāo),進(jìn)而也可以進(jìn)行三維重建。

(28)

4 離散空間點(diǎn)的三角網(wǎng)格構(gòu)建

由于特征點(diǎn)是三維空間中離散的點(diǎn),要采用三角剖分算法把三維點(diǎn)集剖分為若干個(gè)三角網(wǎng)格,得到三角模型的三維幾何外形。在眾多算法中,Delaunay三角剖分算法和模型數(shù)據(jù)結(jié)構(gòu)簡單,適應(yīng)多種狀況,得到了大量應(yīng)用。Voronoi圖是Delaunay三角剖分的基礎(chǔ)。Voronoi圖可以直觀的定義為從每個(gè)點(diǎn)開始,以一樣的速率向外擴(kuò)展直到遇到另一個(gè)點(diǎn)的擴(kuò)展范圍,如圖5中虛線形成的就是Voronoi圖?？梢钥吹絍oronoi圖為凸的多邊形。Delaunay三角網(wǎng)格由Voronoi圖得到,Voronoi圖中每條邊表示Delaunay三角網(wǎng)中每個(gè)點(diǎn)之間的對應(yīng)的連線。圖4中實(shí)線形成的就是Delaunay三角網(wǎng)格。

圖5 Voronoi圖與Delaunay三角網(wǎng)

在得到圖像特征點(diǎn)三維坐標(biāo),并采用Delaunay三角剖分算法對所有數(shù)據(jù)點(diǎn)結(jié)成三角網(wǎng)格后,就可以在三維空間中繪制出目標(biāo)的三維拓?fù)浣Y(jié)構(gòu)。為了增強(qiáng)三維顯示效果,將圖像中的紋理使用OpenGL紋理貼圖技術(shù)貼到物體模型的表面,這樣就可以使重建出的三維模型具有較強(qiáng)的真實(shí)感。

5 實(shí)驗(yàn)結(jié)果與分析

利用典型的膠囊內(nèi)鏡的圖像,設(shè)計(jì)了仿制的胃腸道。利用單目相機(jī)對實(shí)際構(gòu)成的場景進(jìn)行拍攝,得到圖6所示的兩幅模擬膠囊內(nèi)窺鏡在不同位姿下拍攝的圖像。

圖6 待三維重建的兩幅圖像

三維重建的步驟如下:

①特征點(diǎn)的提取。得到52對特征點(diǎn),如圖7所示。

圖7 特征匹配結(jié)果

②相機(jī)相對運(yùn)動(dòng)的計(jì)算。利用特征點(diǎn)匹配的結(jié)果,通過八點(diǎn)算法計(jì)算得到相機(jī)的運(yùn)動(dòng)參數(shù)平移矢量T和旋轉(zhuǎn)矩陣R0。

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

由于八點(diǎn)算法得到的相機(jī)平移矢量T不具有長度信息,本實(shí)驗(yàn)中結(jié)合磁定位的結(jié)果[4]得到相機(jī)具有長度的平移參數(shù)為:

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

③三維坐標(biāo)的計(jì)算。采用式(28)對52個(gè)特征點(diǎn)計(jì)算得到其三維坐標(biāo),部分坐標(biāo)點(diǎn)見表1,其中Xl、Yl、Xr、Yr分別是該特征點(diǎn)在左右兩幅圖中的坐標(biāo)值,X、Y、Z分別為計(jì)算得到的該點(diǎn)在三維空間中的坐標(biāo)點(diǎn)。

④三維網(wǎng)格的剖分。采用Delaunay三角剖分算法對三維點(diǎn)集進(jìn)行網(wǎng)格化。得到的網(wǎng)格結(jié)果如圖8上圖所示,對網(wǎng)格進(jìn)行著色后如圖8下圖所示,可以從圖中清楚的看出場景的相關(guān)結(jié)構(gòu)。

表1 三維坐標(biāo)結(jié)果

(5)OpenGL紋理映射實(shí)現(xiàn)三維重建。

利用C++和OpenGL通過紋理映射技術(shù)顯示在屏幕上。圖9展示了三維重建結(jié)果在4個(gè)視角下的重建結(jié)果,可以看到,三維重建的模型比較直觀,與原場景較為相似。

為了進(jìn)一步驗(yàn)證三維重建的準(zhǔn)確度,在實(shí)驗(yàn)中對若干個(gè)特征點(diǎn)進(jìn)行指定,之后測量這些特征點(diǎn)與相機(jī)平面之間的距離,再與計(jì)算出的結(jié)果進(jìn)行比對。本實(shí)驗(yàn)中,使用內(nèi)窺鏡相機(jī)在(0,0,0)mm和(30,0,5)mm,相機(jī)主軸角度相差20°的前后兩點(diǎn)上進(jìn)行圖像的拍攝;并選定了6個(gè)特殊點(diǎn)進(jìn)行匹配與測量,實(shí)驗(yàn)結(jié)果如表2所示。

圖8 三維坐標(biāo)網(wǎng)格

圖9 三維重建結(jié)果

特征點(diǎn)目標(biāo)實(shí)際深度/mm計(jì)算得到深度/mm誤差/mm相對誤差/%14747.68210.68211.4526464.86480.86481.3538080.89890.89891.124114116.12902.12901.865168171.42863.42862.046241248.27597.27593.01

由表2可以看到,當(dāng)特征點(diǎn)離相機(jī)鏡頭平面的距離在100 mm之內(nèi)時(shí),算法計(jì)算得到的深度與目標(biāo)實(shí)際深度之間的誤差較小,小于1 mm,系統(tǒng)精度較高。相機(jī)與被測特征點(diǎn)的距離大于100 mm時(shí),由于深度距離變遠(yuǎn),誤差也相應(yīng)增大?？傮w來說,在深度250 mm內(nèi),相對誤差控制在3%之內(nèi)。由于膠囊內(nèi)窺鏡圖像的拍攝距離較近,這樣深度已經(jīng)足夠。且對于膠囊內(nèi)窺鏡臨床應(yīng)用,醫(yī)生建議的目標(biāo)測量精度為10%內(nèi),本方法的精度指標(biāo)能夠達(dá)到這一目標(biāo)。

6 結(jié)論

針對膠囊內(nèi)窺鏡的序列圖像,提出了胃腸道三維重建的方法。首先提取前后兩幅序列圖像中相對應(yīng)的多個(gè)特征點(diǎn);基于特征點(diǎn)計(jì)算各特征點(diǎn)在成像面上的坐標(biāo);進(jìn)一步利用8點(diǎn)算法計(jì)算膠囊內(nèi)鏡運(yùn)動(dòng)變化的旋轉(zhuǎn)矩陣和平移矢量,即得到其運(yùn)動(dòng)的位移方向和旋轉(zhuǎn)角度變化。進(jìn)而計(jì)算得到每個(gè)特征點(diǎn)的相對三維坐標(biāo)和世界三維坐標(biāo);然后,采用Delaunay三角剖分算法對各三維點(diǎn)進(jìn)行網(wǎng)格化,并完成場景的三維重建。實(shí)驗(yàn)結(jié)果表明所提出的算法是可行的。

[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] 任宇鵬,胡超,項(xiàng)圣,等. 一種人體移動(dòng)對膠囊內(nèi)窺鏡磁定位干擾的補(bǔ)償方法[J]. 傳感技術(shù)學(xué)報(bào),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] 佟帥,徐曉剛,易成儔,等. 基于視覺的三維重建技術(shù)綜述。計(jì)算機(jī)應(yīng)用研究[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] 羅肖. 序列內(nèi)窺鏡圖像的三維結(jié)構(gòu)重建[D]. 上海:上海交通大學(xué),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-),男,加拿大阿爾伯特大學(xué)博士,浙江大學(xué)寧波理工學(xué)院三江學(xué)者特聘教授;浙江大學(xué)控制學(xué)院博士生導(dǎo)師。從事自動(dòng)化、機(jī)器人控制和定位跟蹤技術(shù)研究。主持國家863計(jì)劃、支撐計(jì)劃、國家自然基金等科研項(xiàng)目30多項(xiàng);發(fā)表論文200多篇。

項(xiàng)目來源:國家自然科學(xué)基金項(xiàng)目(61273332);與寧波市科技項(xiàng)目(創(chuàng)新團(tuán)隊(duì));計(jì)劃項(xiàng)目(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