李文駿 丁 輝 王廣志

(清華大學醫(yī)學院生物醫(yī)學工程系，北京 100084)

引言

超聲成像技術(shù)具有實時、廉價、無電離輻射等優(yōu)點，在臨床診斷和引導下的介入手術(shù)等方面得到了廣泛應用[1-2]。在超聲引導下的介入治療、超聲融合成像和徒手三維超聲成像中[1-3]，為了跟蹤探頭的空間位置姿態(tài)，常常在超聲探頭上固定光學或電磁位置傳感器，通過一系列的空間坐標變換，將超聲圖像和術(shù)前影像、器械位姿等統(tǒng)一到同一個坐標系中。筆者在本文中所討論的超聲探頭空間標定，是指確定其中位置傳感器坐標系到超聲圖像坐標系之間的變換關(guān)系的過程，這對保證超聲引導介入、融合成像或三維重建的系統(tǒng)精度有重要的意義，標定過程中的坐標變換關(guān)系如圖1所示。通過模型中的位置確定N形標記線，在超聲成像中切割出對應的標志點，以及在模型與超聲探頭上定位傳感器的相對位置，可以形成閉環(huán)的空間變換關(guān)系，其中Ts←i就是所需要標定的變換關(guān)系。

圖1 標定過程中的坐標變換關(guān)系Fig.1 Coordinate systems and transforms in calibration

國內(nèi)外的研究者已經(jīng)提出了多種超聲探頭空間標定方法[4]，其中基于N線模型的方法是應用最為廣泛的方法之一[5-6]。這種模型在水箱中設(shè)置多組由繃直的細線構(gòu)成的N形目標，每個N形目標經(jīng)超聲成像形成3個亮斑。拾取圖像上各組亮斑的坐標并確定其對應的N線編號后，就可以利用模型的結(jié)構(gòu)參數(shù)解得最終的標定變換。

目前，拾取亮斑的坐標、判別哪些亮斑屬于同一組N線、確定N線組的編號等工作，是靠人工完成的，這種人工識別和點選亮斑標志點的方法存在一些問題。首先，識別和確定N線編號要求操作者了解N線模型的設(shè)計信息，對醫(yī)護人員或普通工人的實際操作有較大難度;而由于模型中N形目標的數(shù)量較多，即使是有經(jīng)驗的實驗人員也難免出現(xiàn)錯誤判斷N線編號的情況。其次，人工點選超聲圖像上的亮點坐標帶有主觀性，會引入隨機誤差。再次，為獲得較好的標定精度，往往需要獲得多幀圖像上的幾十甚至幾百個N形目標點，人工操作耗時較長，這也更增加了發(fā)生編號錯誤的概率。

為解決上述問題，研究者提出了若干基于N線模型的自動或半自動超聲標定方法。Lindseth等[7]和宋等[8]人工在亮點附近指定初步位置，在其鄰域內(nèi)找到亮度脊，并利用幾何約束進行修正。這種方法減小了選點時的主觀誤差，但并未解決標志點自動編號的問題。Hsu等在模具表面加裝特制的薄膜，通過檢測在超聲圖像上形成的直線，實現(xiàn)了自動標定，但增加了模具設(shè)計和制作的復雜度[9]。Chen等則把N線組的數(shù)目減少到2組，降低了N線識別的難度，但需要拍攝更多幀的圖像[10]。

針對以上問題，筆者在不改變N線模型設(shè)計和實驗流程的條件下，提出一種自動化的超聲探頭標定處理方法，提高了標定的效率和精度。

1 方法

1.1 標定原理

1.1.1 N線模型目標點重建

用于超聲探頭標定的N線模型由多層尼龍線安裝在水槽中組成。這些尼龍靶線構(gòu)成一系列被編號的N線組，每個N線組包括兩條平行線和一條斜線。N線組的成像過程可以抽象為超聲圖像平面與靶線相交。如圖2(a)所示，3條靶線與超聲成像平面相交形成3個交點，即圖中的E、G、F。其中，點F在模型坐標系下的三維坐標可以利用超聲圖像和模型結(jié)構(gòu)信息求得，這個過程被稱為N形目標點的“重建”。從超聲圖像中拾取標志點E、F、G的圖像坐標，計算線段EF與EG長度的比值，有

由于三角形△BEF和△CGF相似，得

對已知編號的N線組，式(2)中B點和C點的坐標可由N線模型的設(shè)計信息確定，由此即可重建出F點在模型坐標系下的三維坐標。

圖2 N形目標點重建。(a)坐標重建原理;(b)實際成像光斑Fig.2 Target point reconstruction of N-fiducials.(a)Schematic;(b)Sample image

在實際圖像中，由于超聲成像平面有一定的厚度，且超聲探頭橫向聚焦并非絕對理想，因此靶線所成的標志點呈現(xiàn)亮斑的形式，如圖2(b)所示。采用人工拾取標志點坐標，一方面難以確定亮斑中心的精確位置，另一方面易受背景噪聲等因素影響而出現(xiàn)錯誤。在超聲成像視場中N線組很多的情況下，人工所記錄的編號也可能出現(xiàn)錯誤，這些都可以通過自動化手段加以改進。

1.1.2 空間變換關(guān)系

標定實驗涉及圖1所示的幾個坐標系，包括世界坐標系{w}、模型坐標系{m}、傳感器坐標系{s}和圖像坐標系{i}。在標定中，各坐標系之間均考慮為三維剛性變換，包含3個旋轉(zhuǎn)參數(shù)和3個位移參數(shù)，共6個自由度。利用齊次坐標，剛性變換可以用矩陣乘法表示。用Τa←b表示從坐標系到坐標系{a}的變換，用xa代表點在{a}系下的坐標，則有

超聲探頭標定的最終目標是求解圖像坐標系到傳感器坐標系的變換，即Ts←i。在實際標定過程中，由三維定位系統(tǒng)的輸出數(shù)據(jù)可得Ts←w。在每組實驗開始時，先用三維定位探針點選出模型外壁和上緣的定位孔，并將小孔在世界坐標系中的位置和在模型坐標系中的位置配準，可以求解出變換Τw←m。對每個完整可見的N線組，用本文1.1.1節(jié)的目標點重建算法，求得F點在模型坐標系下的坐標xm，并由式(4)變換到傳感器系，即xs。另外，由預先標定好的超聲像素尺寸，可以直接從超聲圖像中求得標志點在圖像坐標系下的坐標xi。由3個以上的(xs-xi)坐標對，可求解標定變換Ts←i。在最小二乘意義下，即為

至此，標定問題就抽象為一個已知配對關(guān)系的點對剛性配準問題，可利用Horn等提出的解析算法[1-2]來求解。

1.2 自動標定算法

從圖像中，自動提取亮斑中心準確的坐標，并確定其對應的N線組的編號，是實現(xiàn)基于N線模型的超聲探頭自動標定方法的關(guān)鍵。自動提取亮斑的亮度中心，能消除主觀點選坐標帶來的隨機誤差，提高標志點坐標的精確度。另外，進一步通過自動算法確定各標志點對應的靶線編號，能有效防止實驗者在確定N線編號時偶然出錯的情況。

然而，實現(xiàn)N線組的自動編號需要克服一些困難。為了減小在超聲成像平面法線方向上的標定誤差，應使探頭在不同角度下獲得的圖像均包含數(shù)量較多的N線目標，并且其分布應盡量遍布整個超聲圖像范圍，因此在模型中通常要設(shè)置較多的N線組[9]。所設(shè)計的模型包含67條靶線(30個N線組，部分靶線由兩個相鄰的N線組共用)，每幀圖像中有30～50個亮斑標志點，對應8～16個N線組。也就是說，每幀圖像都只能觀察到模型中的一部分靶線，而這部分靶線的編號和位置都不能事先知道。此外，圖像上出現(xiàn)的亮斑粘連、缺失，以及由于超聲在容器壁多次反射等原因造成的偽影，也是干擾N線標志點識別的因素。針對上述N線模型及其超聲圖像的特點，設(shè)計了一套自動標定算法，其流程如圖3所示。

圖3 自動標定算法流程Fig.3 Flow chart of the automatic calibration method

1)預處理。預處理步驟的目的是初步提取圖像中的標志點亮斑。首先，在本實驗的條件下，以96作為灰度閾值，對圖像進行二值化。由于超聲成像平面的厚度效應和超聲橫向聚焦不理想等原因，標志點亮斑常沿著超聲束的切向擴散，呈橫向扁平狀;而容器壁的反射像(如圖4(a)右下部的高亮線段)、超聲波多次反射相干形成的條帶狀偽影(見圖4(a)下部)和一般的圖像噪聲，則沒有上述橫向擴散的特征。因此，遍歷二值化后圖像的各個連通域，將豎直方向上延展長度超過30像素的連通域刪除。用調(diào)整后的二值圖作掩膜，提取出可能的標志點區(qū)域圖像。預處理后的圖像如圖4(b)所示。

2)計算標志點預測位置。這一步是利用已知信息來預測靶線標志點在圖像中的位置，而這些預測值將作為后面搜索步驟的初始位置。該過程用到了N線模型的設(shè)計形狀尺寸信息，由此得到的標志點預測位置與模型靶線是一一對應的，且編號已知，這為后續(xù)步驟中自動確定N線組編號奠定了基礎(chǔ)。

圖4 自動標定算法各步驟處理結(jié)果。(a)原始圖像;(b)預處理結(jié)果;(c)標志點預測位置;(d)鄰域搜索獲得的亮斑中心位置;(e)含編號的N線組識別結(jié)果Fig.4 Results of the automatic calibration algorithm.(a)Original image;(b)Results of preprocessing;(c)Predicted positions of N-fiducials;(d)Results of neighborhood search that represent refined positions of N-fiducials;(e)Results of N-fiducial identification

處理第一幀圖像和后續(xù)其他圖像的方法有所不同。采集第一幀圖像時，讓探頭對準N線模型的中間位置，并保持超聲成像平面與水模壁基本平行。由此所得的圖像標志點的分布樣式能被預先確定，從而能給出標志點的預測位置。處理第二幀及之后的各幀圖像時，之前已經(jīng)由前面的圖像解得標定變換Ts←i的一個初始解;利用這個初始解，可以由三維定位傳感器的讀數(shù)來預測超聲成像平面的位置和姿態(tài)，求出這個虛擬截面與模型靶線的交點，即為標志點的預測位置。為了數(shù)學上的方便，將模型靶線兩端的控制點坐標P1和P2轉(zhuǎn)換到圖像坐標系下，有:

圖5 線段被x-O-y平面所截Fig.5 Line segment intersecting the x-O-y plane

在式(6)中，Τs←w由三維定位系統(tǒng)的實時讀數(shù)直接得到，Τw←m則在每組實驗開始時配準得到。如圖5所示，在圖像坐標系下，超聲成像平面就是x-O-y平面。設(shè)靶線端點坐標分別為P1(x1，y1，z1)和P2(x2，y2，z2)，則模型靶線的直線方程可寫為

即為該靶線與超聲成像平面交點的坐標。注意，交點必須在兩端點之間，其限制條件為

計算每條靶線與虛擬超聲平面的交點，并排除不滿足限制條件的情況，就得到各標志點坐標的預測值，如圖4(c)所示。

3)鄰域搜索和N線編組。在各預測位置的一定大小的鄰域內(nèi)搜索亮度超過閾值的點。如果不存在，說明該點沒有被成像，則該組N線無法用于標定，自動加以排除;如果存在，則找到鄰域內(nèi)面積最大的亮斑，并求出其亮度中心，認為是該標志點的坐標，如圖4(d)所示。

通過以上處理，所搜得的標志點對應的靶線編號就是已知的，能方便地確定哪些標志點屬于同一個N線組。為了避免搜索到圖像中的噪聲或偽影，或標錯靶線編號，處理中進一步剔除N線組三點共線性太差(連線夾角大于3°)以及其中某兩點重合或接近重合(距離小于10像素)的情況。對符合要求的N線標志點，標記其N線編號，并記錄下對應的k值供重建時使用，如圖4(e)所示。

4)重建和解算。至此，獲得了一幀圖像中所有能夠找到的N線組編號及對應的k值。參照文中“標定原理”部分的方法，可以得到N線目標點在傳感器坐標系和圖像坐標系下的坐標，即增加了式(5)中新的(xs-xi)坐標對。利用Horn算法[11-12]，可以更新最終的標定變換Τs←i。

一幀處理完畢后，判斷是否已處理完全部的標定圖像。若否，讀取下一幀圖像，再次更新標定變換，直到所有圖像處理完畢并輸出結(jié)果。

1.3 實驗驗證

為了驗證上述自動標定算法的性能，使用自行設(shè)計的N線模型進行了標定實驗驗證。模型用有機玻璃制成，外尺寸為28 cm×11 cm×24 cm;模型內(nèi)部分布了7層深度為4～16 cm的靶線，由直徑0.14 mm的尼龍線穿成并繃緊，每層有4～5個N線組;在模型外壁和上緣分布有多個定位孔，定位孔及穿線孔的加工精度均達到0.1 mm;在模型內(nèi)壁和底面放置吸聲材料，以減小反射噪聲。

筆者采用的超聲圖像采集系統(tǒng)包括:超聲成像儀(DC-6，邁瑞醫(yī)療)，配臨床常用的3.5 MHz腹部探頭;交流脈沖電磁三維定位系統(tǒng)(Aurora，Northern Digital Inc.);計算機(xw8400，HP)。在計算機端，用視頻采集卡捕捉超聲圖像，并進行數(shù)字化。與圖像同步的超聲探頭三維定位數(shù)據(jù)，通過串口發(fā)送到計算機。實驗中，超聲成像的深度為17.5 cm。標定時的傳聲介質(zhì)為室溫下去除氣泡的純凈水。由于標定的結(jié)果將應用于人體組織的成像和測量，因此使用超聲儀自帶的液性聲速補償功能，可減小掃描線方向距離測量的誤差。

標定實驗共進行6組，每組包含7幀圖像。每組的第1幀圖像均在探頭正直地對準水槽中央位置時采集，并使水模第4層N線組位于圖像的中央位置，其余6幀圖像對應的探頭位置和角度不固定。

為了通過實驗考察標志點自動識別的成功率，驗證自動識別的可靠性，筆者從精密度(precision)和準確度(accuracy)兩方面來驗證新算法的精確度。精密度是指算法在多次實驗中所得結(jié)果的一致程度，即重復性;準確度則是指測量結(jié)果與真值之間的誤差大小。

采用一個直觀評價精密度的方法，就是考察標定變換矩陣各元素在多次實驗中的標準差[7]。但是，由于各研究中實際的標定變換互不相同，這種方法并不適用于標定方法的橫向比較。以往廣泛采用的評價方法是將超聲圖像中的某個固定點(如圖像中心和4個角點)變換到傳感器坐標系下，統(tǒng)計在多次實驗所得的標定變換下點云的散布情況。記某點在超聲圖像坐標系下的坐標為xi，用各次實驗所得的變換將其變換到傳感器坐標系下的坐標均值為，則標定精密度可定義為

常用目標配準誤差(target registration error，TRE)作為標定算法準確度的評價指標[5，7，9]。每組實驗分別抽出一幀圖像，將該幀圖像內(nèi)的N形目標點作為測試用的注冊目標。利用其余各幀圖像(訓練集)得到標定變換矩陣，計算測試幀中各目標點由所得的變換轉(zhuǎn)換到位置傳感器系下的坐標，與由模型結(jié)構(gòu)信息重建所得坐標間的歐氏距離的方均根，即

將自動方法與人工方法進行對比，給出了自動標定算法在識別正確率、標定精度和運行速度方面的結(jié)果。

3 結(jié)果

3.1 標志點自動識別情況

由于實際標定所采集的超聲圖像存在噪聲、模糊和對比度不均勻等圖像退化，因此人工點選和自動搜索都有未被提取的目標。對實驗用的6組共42幀圖像，通過自動識別或人工點選的方法，總共確定了461個N形目標，其中被自動算法成功識別的有414個，占總數(shù)的89.8%。N形目標識別結(jié)果的分類統(tǒng)計如表1所示。對所有自動識別得到的N形目標點進行了檢查，未發(fā)現(xiàn)誤識別(即把噪聲誤認為靶線標志點)或誤編號的情況。自動算法未成功識別的目標點占總數(shù)的10.2%，數(shù)量可以接受，一般是由于標志點附近有較強的反射回聲偽影引起的。比較自動識別和人工點選的目標點坐標，兩者平均距離為1.25像素。以上結(jié)果說明，標志點自動識別是可靠的。

表1 自動和人工選點結(jié)果分類統(tǒng)計Tab.1 Classification of automatically and manually identified N-fiducials

3.2 精確度驗證結(jié)果

如前所述，精確度包含精密度(precision)和準確度(accuracy)兩個方面。精密度是指算法在多次實驗中所得結(jié)果的一致程度，即重復性;準確度則是指測量結(jié)果與真值之間的誤差大小。

根據(jù)式(11)，將采集的6組實驗(每組包含7幀圖像)由自動標定方法所得的精密度結(jié)果與傳統(tǒng)方法(人工點選N形標志點，經(jīng)坐標重建后利用解析算法，求解標定變換[5])的進行比較，如表 2所示。

表2 標定精密度對比Tab.2 Comparison of calibration precision

另外，式(12)所得出的準確度結(jié)果(見表3)表明，新算法與傳統(tǒng)算法相比，目標配準誤差顯著減小，標定準確度有所提高。

表3 目標注冊誤差比較Tab.3 Comparison of target registration error(TRE)

隨著參與計算的圖像數(shù)量的增加，標定精確度逐漸得到改善，并最終收斂到一個較小的值上。以第4組實驗為例，給出目標注冊誤差隨幀數(shù)的關(guān)系曲線，如圖6所示。采用新算法得到的TRE隨幀數(shù)增加呈指數(shù)曲線下降，且收斂速度比傳統(tǒng)算法更快。從第6幀到第7幀，TRE下降的比例為1.5%，表明已經(jīng)基本收斂。

3.3 運行速度

本研究的自動標定算法全部在Matlab R2011a平臺上實現(xiàn)。在搭載Intel?CoreTMi5-2450M CPU@2.50 GHz、4 GB內(nèi)存的計算機上運行，處理每幀圖像所需的時間平均為0.067 s，最大為0.13 s。相比人工逐個點選標志點，大大提高了效率。

4 討論

圖6 目標注冊誤差隨幀數(shù)的關(guān)系曲線Fig.6 Target registration error(TRE)in terms of number of frame

筆者查詢了近年來國內(nèi)外利用N線模型進行超聲探頭標定且被標探頭參數(shù)與本方法相近的文獻，本方法所得標定精確度與相關(guān)文獻比較的情況如表4所示?？梢姡痉椒ǐ@得的精密度和準確度與國內(nèi)文獻比較均有較為明顯的提高;在設(shè)備相對落后的條件下，超過或接近了國外文獻報道的水平。

表4 本方法標定精確度與近年國內(nèi)外文獻的比較Tab.4 Comparison of calibration precision and accuracy to the recent literature

筆者進一步分析了自動識別標志點提高標定精確度的原因。對由自動識別或人工點選得到的所有標志點進行分類統(tǒng)計，如表1所示。其中，被自動方法選出但在人工點選時被略去的標志點約占13.0%，主要分布在超聲圖像的上部和邊緣位置。這部分標志點一般是亮斑亮度偏暗、面積偏小，人工點選時難以辨識，因此只好被放棄。首先，自動方法利用了這些標志點，使N形目標點更均勻地覆蓋到超聲圖像的全部范圍，這是提高最終標定精確度的主要原因之一。其次，自動方法通過計算標志點亮斑的中心位置，消除了人工點選標志點過程中引入的隨機位置誤差。對自動和人工方法均選出的標志點(占總數(shù)的76.8%)，由兩種方法所確定的坐標的平均距離為1.25像素，這反映了自動方法在消除隨機誤差上的貢獻。再次，自動方法所得標志點的共面性比人工點選更優(yōu)。標志點是模型靶線被超聲成像平面截得的交點，理想情況下每幀圖像所得的標志點應當是共面的。但由于誤差的存在，實際重建的標志點并不完全共面;研究表明，共面性的丟失是影響最終標定精度的重要因素[13]。以第4組實驗為例，比較了在各幀圖像中由自動和人工方法分別得到的標志點坐標與其擬合平面間的平均距離和最大距離，如表5所示?？梢姡詣臃椒ㄋ脴酥军c的共面性有顯著改善，從而提高了最終的標定精度。

表5 自動和人工方法所得標志點坐標與其擬合平面間的距離比較Tab.5 Statistics of distance from reconstructed fiducials to the fitting plane

筆者通過本研究提出了自動標定方法，也提高了超聲探頭標定的效率。采用傳統(tǒng)的人工點選方法，實驗者點選處理一幀圖像的時間一般需要幾分鐘。而如本文3.2節(jié)所述，自動標定算法每處理一幀的時間平均只需0.67 s，因此大大節(jié)約了人工點選標志點所需的時間。

5 結(jié)論

超聲探頭標定是確定超聲圖像與三維定位系統(tǒng)傳感器之間的空間變換關(guān)系的過程，是超聲引導介入手術(shù)、超聲融合成像、徒手三維超聲等應用中的基礎(chǔ)技術(shù)環(huán)節(jié)。面向目前廣泛使用的基于N線模型的標定方法，本研究對傳統(tǒng)的人工點選標志點的方法進行了改進，提出了一種新的自動化方法。通過實驗驗證，證實新方法能提高超聲探頭標定的精度和效率。

[1]Xu S，Krucker J，Turkbey B，et al.Real-time MRI-TRUS fusion for guidance of targeted prostate biopsies[J].Computer Aided Surgery，2008，13(5):255-264.

[2]Solberg O，Lango T，Tangen G，et al.Navigated ultrasound in laparoscopic surgery[J].Minimally Invasive Therapy＆ Aided Technologies，2009，18(1):36-53.

[3]Prager R，Rohling R，Gee A，et al.Rapid calibration for 3-D freehand ultrasound systems[J].Ultrasound in Medicme ＆Biology，1998，24(6):855-869.

[4]Mercier L，Lango T，Lindseth F，et al.A review of calibration techniques for freehand 3-D ultrasound systems[J].Ultrasound in Medicine＆ Biology，2005，31(2):143-165.

[5]Pagoulatos N，Haynor D，Kim Y.A fast calibration method for 3-D tracking of ultrasound images using a spatial localizer[J].Ultrasound in Medicine＆ Biology，2001，27(9):1219-1229.

[6]羅楊宇，徐靜，魯通，等.基于磁定位器的手動三維超聲圖像標定[J].中國生物醫(yī)學工程學報，2008，27(2):250-254.

[7]Lindseth F，Tangen G，Lango T，et al.Probe calibration for freehand 3-D ultrasound[J].Ultrasound in Medicine＆ Biology，2003，29:1607-1623.

[8]宋章軍，徐靜，陳懇.一種快速、自動的三維超聲標定方法[J]. 北京生物醫(yī)學工程，200928(6):601-605，613.

[9]Hsu PW，Prager R，Gee A，Treece G，et al.Real-time freehand 3D ultrasound calibration.Ultrasound in Medicine ＆ Biology，2008，34(2):239-251.

[10]Chen TK，Thurston A，Ellis R，et al.A real-time freehand ultrasound calibration system with automatic accuracy feedback and control[J].Ultrasound in Medicine ＆ Biology，2009，35(1):79-93.

[11]Horn BKP.Closed-form solution of absolute orientation using unit quaternion[J].J Opt Soc Am A，1987，4(4):629-642.

[12]Horn BKP，Hilden HM，Negahdaripour S.Closed-form solution of absolute orientation using orthogonal matrices[J].J Opt Soc Am A ，1988，5(7):1127-1135.

[13]朱立人，李文駿，丁輝，等.一種提高N線模型探頭標定精度的解算方法[J].中國生物醫(yī)學工程學報，2012，31(3):337-343.