楊凡帆,張勤儉,范巨峰,李海源,晁明輝,褚浩杰

(1.北京信息科技大學 儀器科學與光電工程學院,北京100192;2.北京信息科技大學 機電工程學院,北京100192;3.首都醫(yī)科大學附屬北京朝陽醫(yī)院,北京100020;4.北京郵電大學 自動化學院,北京100876)

0 引言

超聲成像是目前臨床診斷和手術中最流行的醫(yī)療成像方式之一[1],將超聲影像結合醫(yī)學圖像處理技術和三維重建技術可以獲得更豐富的病灶信息。近年來,徒手三維超聲成像因其掃描方式靈活、能提供更大的視野以及便于醫(yī)生操作,而成為超聲影像引導手術的主要研究方向,在圖像引導類以及導航類的微創(chuàng)手術中具有實際意義。

常用的三種獲取超聲影像的方式為機械式掃描[2]、徒手掃描(帶位置傳感器[3-4]或不帶傳感器[5-6])和容積探頭方式[7]。其中,帶位置傳感器的徒手超聲三維重建方式魯棒性更強、更方便、成本更低,目前使用更為廣泛[8]。這種方式通常是由傳統(tǒng)的超聲設備結合某種空間定位裝置來獲取一系列二維超聲圖像及其空間方位信息,再結合標定技術和三維重建技術,實現(xiàn)超聲影像的三維重建。市面上已經有一些較為成熟的空間定位方法,包括機械(臂)定位、聲學定位、電磁定位和光學定位[9]。其中,電磁定位體積小、無遮擋、精度和價格適中,比較適用于術中手術設備需要移動空間且移動范圍不大的穿刺手術。

徒手掃描方式獲取到的是空間中非規(guī)則分布的二維序列圖像,與計算機斷層掃描(computed tomography,CT)、磁共振成像(magnetic resonance imaging,MRI)等獲取的規(guī)則分布的斷層圖像不同,超聲三維重建需要先計算二維圖像到超聲探頭的轉換信息和一系列空間坐標轉換,這一過程即為超聲探頭標定。探頭標定的基本原理是通過掃描一個已知幾何特征及其空間位置的模型,根據(jù)這些特征在超聲圖像中和在實際模型中的位置關系,來求解未知標定矩陣中的參數(shù)。標定流程包括標定模型制作、圖像特征提取和標定方程求解3個步驟[10]。根據(jù)標定模型內特征的形狀可將標定模型分為點式、線式和面式三大類,經典的標定模型有交叉線形、三線形、N線形等[11-14]。N線形模型制作簡單,便于掃描,且應用廣泛,符合本文對重建的要求。

徒手超聲重建是將超聲探頭掃描圖中的二維像素點重建到三維體數(shù)據(jù)中,重建流程包括:構建體數(shù)據(jù)、像素的重新分配和計算體數(shù)據(jù)值。構建體數(shù)據(jù)是預先確定一個三維重建空間,可以通過主成分分析法和包圍盒等方式。分配像素值是根據(jù)求解的坐標變換,將二維平面上的像素插入三維空間中,可以利用不同的插值方法將像素分配到重建區(qū)域。計算體數(shù)據(jù)值是對插值后的三維體進行填充,將空值補充完整,得到一個完整的三維數(shù)據(jù)。如何在圖像質量和計算時間之間取得平衡是目前三維成像重建算法面臨的挑戰(zhàn)[15]。

1 設備和坐標系統(tǒng)

圖1展示了超聲重建過程,先掃描標定模型完成探頭標定,再掃描待重建區(qū)域進行三維重建。

圖1 磁導航徒手超聲三維重建流程Fig.1 Process of freehand ultrasound 3D reconstruction based on magnetic navigation

通過視頻采集卡采集B超探頭畫面,將磁導航位置傳感器固定在B超探頭上,在模型表面模擬醫(yī)生掃體,每一幅超聲圖像都對應一個空間方位數(shù)據(jù)(x,y,z,α,β,γ)。其中,(x,y,z)代表位置傳感器相對于電磁發(fā)射器的坐標,單位為mm;(α,β,γ)為位置傳感器相對于電磁發(fā)射器的方位角、俯仰角和滾轉角,單位為°。利用方向和位置信息構建三維重建空間,再對該區(qū)域進行插值重建,最后將三維數(shù)據(jù)可視化得到超聲的三維影像。

超聲影像重建中涉及的坐標系如圖2所示。空間導航設備的位置傳感器被固定在一個傳統(tǒng)的超聲探頭上。重建需要將二維超聲圖像中的像素點轉化到三維模型空間。具體過程是先將圖像(像素)坐標系(用p來表示)轉化到傳感器坐標系(用s來表示),然后轉化到電磁發(fā)射器坐標系(默認為世界坐標系,用w來表示),再轉換到模型坐標系(用m來表示),由此可將像素點放入三維模型空間。得到掃描序列圖之間的實際物理關系,重建結果才有實際意義。

圖2 坐標轉換過程Fig.2 Coordinate conversion process

將超聲幀中的一個像素(u,v)轉換到模型中的一個點(x,y,z)的過程如式(1)所示(默認模型坐標系為重建坐標系)。

M=Tm←w·Tw←s·Ts←p·P

(1)

式中:Ta←b代表從坐標系b向坐標系a的空間變換矩陣;P為像素點在二維超聲圖中的坐標;M為像素點在模型中的坐標。

2 探頭標定實驗

2.1 標定模型

標定模型選擇使用最廣泛且更適用于臨床的N線形模型,本文基于N線形模型思路提出一種“三層簡易N線”標定模型,以少量的點計算出圖像中多組橫縱方向上的數(shù)據(jù),以此降低標定環(huán)節(jié)的復雜度。標定模型的殼體是使用樹脂材料通過三維打印的方式制成,4條魚線分為3層固定在殼體內部,第1層是兩條平行線,兩線間距|AC|=|BD|=25 mm,其余兩層是方向相反的斜線,每層之間距離相等,間距為10 mm。以模型中心位置為例進行說明,掃描時超聲探頭方向如圖3(a)所示,探頭與AB垂直且與AC平行,E、F、G、H為超聲平面在模型中心位置與魚線的交點,對應的俯視示意圖如圖3(b)。掃描結果如圖3(c),提取標定點后得圖3(d)。為便于后續(xù)的說明,連接GH,點F投影到線段GH上得到點F′,可知|EF|=|FF′|=1/2|EF′|。

圖3 標定模型和標定點Fig.3 Calibration model and calibration point

2.2 標定計算

本文將點的連通域中心作為其坐標,標定的目的是計算不同坐標系之間的變換關系。

變換矩陣Tw←s可由位置傳感器的輸出數(shù)據(jù)(x,y,z,α,β,γ)得到,如式(2)所示:

(2)

將式(1)展開為如下形式:

(3)

式中:sx和sy分別為橫向像素比例和縱向像素比例,兩者均未知;Tm←w和Ts←p未知。未知參數(shù)和矩陣的求解順序為:首先,通過標定點之間的關系求出超聲圖像的像素比例sx和sy;然后,通過模型中的點在傳感器坐標系下的對應關系求出Tm←w;最后,結合之前的求解結果,通過標定點在圖像和模型中的對應關系求出Ts←p。

2.2.1 求解像素比例

超聲圖像的像素比例單位為mm/像素,即圖中的一個像素代表多少毫米。不同種類探頭的像素比例不同,且一幅超聲圖中橫向和縱向的比例也不同,計算式為

(4)

式中:X為超聲圖像中像素的橫向距離,x為模型中與之對應的距離;Y為超聲圖像中像素的縱向距離,y為模型中與之對應的距離。

從超聲圖像中計算出的距離和在標定模型中測量得到的距離均存在一定誤差,在實際求解像素比例時可選取n組距離進行計算,求其平均值作為最后的結果,如式(5)所示:

(5)

2.2.2 求解Tm←w

首先通過測量得到標定模型上的點在模型坐標系下的坐標M和在傳感器坐標系下的坐標S,然后通過轉換矩陣Tw←s將S轉換為電磁發(fā)射器坐標系下的坐標W,最后利用最小二乘法求解式(6)中的坐標轉換矩陣Tm←w。

M=Tm←w·W

(6)

2.2.3 求解Ts←p

對式(1)移項得式(7),可用多幅圖中的F點在圖像和模型中的坐標關系求出Ts←p的最小二乘解。

(7)

式中:PF為F點在超聲圖像中的坐標(sxu,sxv),可直接通過計算得出。MF為F點在模型中的坐標(xF,yF,zF),可利用相似三角形原理得到式(8),則(xF,yF,zF)如式(9)所示,其中z坐標值固定為15.5 mm。

(8)

(9)

至此,標定環(huán)節(jié)結束,得到了圖像像素點到模型坐標系下的變換。由此空間坐標變換可以將像素點分配到三維重建空間。

3 三維重建及可視化

包圍盒(bounding box)技術是構建重建空間的簡易方法,它通過填充一系列二維超聲幀的像素,再找到最大點和最小點來計算體積,構成一個由空體素組成的三維體?；诖?本文提出簡易快速包圍盒方法:在一次掃描得到的序列中直接利用所有坐標中的最小值和最大值快速確定重建區(qū)域范圍。如圖4所示。

圖4 簡易快速包圍盒算法流程Fig.4 Simple and fast bounding box algorithm flow

根據(jù)插值方式的不同可以將超聲影像三維重建分為基于體素的方法(voxel based method,VBM)、基于像素的方法(pixel based method,PBM)和基于函數(shù)的方法(function based method,FBM)。VBM遍歷每個體素,選擇其周圍的像素值填入,常用的最近鄰域體素法(voxel nearest neighbor,VNN)將距當前體素最近的像素的值插入。FBM對二維超聲圖像和其對應的三維重建空間建立函數(shù)關系,將像素值映射到目標體素中去。PBM遍歷圖像中每個像素,將像素值分配給體素。最鄰近像素算法(pixel nearest neighbor,PNN)是一種經典的PBM,也是本文使用的重建方式。圖5展示了PNN的兩個基本步驟:像素重新分配和孔洞填補。步驟①遍歷二維像素值,通常與權值一起存儲,計算最鄰近的一組像素集合的加權平均值賦值給三維體素。步驟②遍歷三維體素,對鄰域內非空體素進行加權平均后賦值給相應的空體素點。

圖5 PBM流程示意圖Fig.5 Diagram of the PBM process

重建的結果是一個三維格式,常見的可視化方法分為兩大類:基于面繪制的方法和基于體繪制的方法。本文采用經典的光線投射(ray casting,RC)算法進行可視化。當光線穿過體數(shù)據(jù)時,RC方法會沿著光線方向對經過的每個體素的顏色和阻光度進行積累,最后在平面上繪制出投影。

4 實驗

實驗采用NDI的3D Guidance trakSTAR 電磁導航定位儀,深圳邁瑞公司DP-7 超聲設備,使用epiphan(DVI2USB 3.0)的視頻采集卡連接電腦采集B超探頭畫面。

本文使用圖3(d)計算像素比例,將相應的數(shù)值代入式(5)計算像素比例得式(10):

(10)

式中:XGF、XFH、XGE、XEH、XGH表示超聲圖像中兩點之間的橫向距離;YEF、YFG、YFH、YEG、YEH表示超聲圖像中兩點之間的縱向距離。由圖3(a)可知,上述距離對應其在標定模型中的實際距離為12.5 mm、12.5 mm、12.5 mm、12.5 mm、25 mm和10 mm、10 mm、10 mm、20 mm、20 mm。

表1記錄了求解坐標轉換矩陣Tm←w時所選取的標定模型上的8個點分別在兩個不同坐標系下的空間位置。在計算時需要統(tǒng)一單位。先利用Tw←s將S下的點轉換到W下,再將表1中的點帶入式(6)中求出Tm←w。

表1 點在模型坐標系下的位置和用傳感器測得的位置

通過式(8)和式(9)計算出F點在模型中的坐標MF,在求解Ts←p時需計算多幅超聲圖像的F點在模型中的坐標并帶入式(7)求最小二乘解。至此得到標定關系Tm←wTw←sTs←p。

標定探頭時可將探頭固定在滑軌上,保持探頭平行于AC和BD且垂直于模型底面,由AC邊掃描至BD邊。選取掃描后較為清晰的31幅超聲圖像,其中10幅圖像用于上述標定計算,剩余21幅圖像用于測試標定誤差。比較F點在標定模型中的實際坐標值與計算坐標值,以此來測試探頭標定誤差。圖6為誤差曲線,x軸的平均誤差為3.28 mm,y軸的平均誤差為0.53 mm,z軸誤差忽略不計(探頭固定在滑軌上平行移動)。以較小的婦科腫瘤為例,腫瘤直徑一般是2～3 cm,穿刺的誤差要求通常在毫米級別內,醫(yī)生通常會努力將穿刺誤差控制在2 mm或更小的范圍內。本文探頭標定誤差(不含z軸)平均為1.90 mm,初步達到了手術需求。

圖6 誤差分析Fig.6 Error analysis

因為不同病灶區(qū)和影像設備的清晰度不一樣,像核磁或CT影像清晰且易分割,但超聲畫面?zhèn)斡按蟆⒃朦c多,往往很模糊。若提前將目標區(qū)域分割出來會減少重建壓力。本文通過點選出目標區(qū)域,再利用樣條曲線擬合輪廓后進行閾值分割。RC方法和最大密度投影(maximum intensity projection,MIP)方法是基于體繪制的兩種常見方法。本文用球體模型進行實驗,用兩種可視化方法對比分割前后的重建效果。結果如圖7所示。

圖7 目標模型的可視化結果Fig.7 Visualization of the target model

根據(jù)結果可以明顯看出:模型內部有干擾介質,未經分割的模型不能很好地展示目標區(qū)域,分割后重建效果更好,且RC方法的可視化效果更好。

5 結束語

圖像引導的微創(chuàng)介入手術逐漸成為當今接受度更高的手術方式之一,也一直是國際前沿的研究熱點之一。本文對徒手超聲三維重建的關鍵技術進行了研究,主要包括超聲標定、超聲重建及三維可視化等內容。研究改進了標定模型,簡化了標定環(huán)節(jié),并優(yōu)化了構建重建區(qū)域的方法,在一定程度上簡化了整個徒手超聲重建過程。

為了提高重建效率,可以考慮從以下幾個方面進行優(yōu)化,如:在數(shù)據(jù)采集完畢后對圖像進行增強或分割來優(yōu)化重建范圍;優(yōu)化標定模型,降低標定環(huán)境的復雜度;對重建和可視化的算法公式進行優(yōu)化。

在未來的工作中,將會采集更多超聲數(shù)據(jù),并把這種簡易的標定思路運用到自動標定的研究中。此外,還會基于超聲影像研究自適應的分割方式,以此來提高重建的效率。