基于二維超聲圖像的肝臟運動跟蹤方法研究

龐 杰，周著黃，吳水才

（北京工業(yè)大學(xué)環(huán)境與生命學(xué)部生物醫(yī)學(xué)工程系，北京 100124）

0 引言

目前，原發(fā)性肝腫瘤是我國第4 位常見惡性腫瘤及第2 位腫瘤致死病因，嚴(yán)重威脅我國人民的生命和健康[1]。肝切除、肝移植和影像引導(dǎo)下的熱消融是目前肝腫瘤最有效的治療方法[2]。由于肝腫瘤早期發(fā)現(xiàn)困難，確診時已是晚期，因此常常錯過最佳治療機會，再加上沒有足夠的肝源等原因，只有小部分患者可以進行肝移植和肝切除治療。熱消融術(shù)因其侵入性小、普適性好、經(jīng)濟有效等原因，在肝腫瘤治療領(lǐng)域逐步推廣。

微波/射頻熱消融治療是肝腫瘤的主要治療方法之一。雖然已有不少研究人員開發(fā)了引導(dǎo)手術(shù)、立體定向?qū)Ш?、機器人輔助穿刺[3-5]等熱消融輔助技術(shù)和方法，但其中大部分的模型都是通過術(shù)前CT 圖像建立，并且假設(shè)肝臟處于靜止?fàn)顟B(tài)。雖然在術(shù)中患者處于麻醉狀態(tài)，但人體呼吸運動引起的腫瘤位置改變給熱消融治療帶來了很大困難[6-9]。為了解決呼吸運動對腫瘤定位造成的影響，早期研究人員提出了外部標(biāo)記法[10-11]來應(yīng)對。該方法通過在患者身上放置傳感器，并對患者進行CT 掃描，記錄體外傳感器的位置和對應(yīng)的組織位置，通過建模從體外標(biāo)志點的位置求取組織位置。但該方法有一些缺點，模型需要在術(shù)前建立，而且手術(shù)過程中該模型可能會發(fā)生偏移。而人體內(nèi)部器官運動復(fù)雜，外部傳感器不能完全對應(yīng)，模型產(chǎn)生的輕微誤差可能導(dǎo)致手術(shù)失敗[12]。

許多醫(yī)學(xué)成像技術(shù)被應(yīng)用于臨床以輔助醫(yī)生治療，包括平面X 射線、MRI 和超聲。雖然X 射線可以實時使用，但由于有電離輻射，不能長時間跟蹤胸腹部器官的運動。MRI 雖然可以提供高分辨力和很好的組織對比度，但其設(shè)備占地面積較大，移動不便且很難實時產(chǎn)生圖像，不適合跟蹤胸腹部器官的運動。隨著計算機技術(shù)的發(fā)展，超聲圖像采集質(zhì)量不斷提高，且超聲所具有的實時采集圖像和無電離輻射的優(yōu)勢，使其成為非侵入式治療目標(biāo)跟蹤的一個很好的選擇。Hallack 等[13]使用logDemons 非線性配準(zhǔn)，并通過感興趣點周圍血管大小自適應(yīng)選取感興趣區(qū)域（region of interest，ROI）對感興趣點進行跟蹤，其平均跟蹤誤差為1.21 mm。缺點是感興趣點在相鄰兩幀中位移大小和ROI 邊長大小近似時容易失敗，且其算法運算速度較慢，平均幀率為12 幀/s。Bharadwaj等[14]對孿生神經(jīng)網(wǎng)絡(luò)進行適當(dāng)改進，在超聲圖像中跟蹤目標(biāo)，其平均跟蹤誤差為1.60 mm，在感興趣點接近超聲圖像邊界、被遮擋以及對比度不高的情況下容易跟蹤失敗、產(chǎn)生較大誤差。該研究通過添加模板匹配模塊和線性卡爾曼濾波模塊進行調(diào)整，使得誤差更低。以上方法表明，需要更高的跟蹤精度以及可以實時運行的跟蹤算法來對肝臟運動進行跟蹤，從而減少消融手術(shù)過程中呼吸運動帶來的誤差影響。

為此，本文提出一種中值流[15]、模板匹配、特征點匹配融合的跟蹤方法，該方法通過超聲圖像跟蹤肝臟運動并進行補償，減少呼吸運動帶來的影響。跟蹤結(jié)果表明該方法較為準(zhǔn)確，且運行速度快，能夠滿足實時要求，可以輔助手術(shù)導(dǎo)航系統(tǒng)校正呼吸運動帶來的影響，提高手術(shù)治療效果。

1 方法

1.1 實驗數(shù)據(jù)與平臺

本研究在Visual Studio 2017 上進行，采用CLUST（Challenge on Liver Ultrasound Tracking）2015 2D[16-19]數(shù)據(jù)集進行測試和驗證。該數(shù)據(jù)集目前包含24 個訓(xùn)練集和39 個測試集。其中有4 種類型超聲圖像序列，不同類型的超聲圖像對比度、圖像質(zhì)量、感興趣點位置分布不同，每個超聲圖像類型至少包含4位受試者，每位受試者有1～5 個感興趣點，感興趣點通常為血管中心。將本方法的跟蹤結(jié)果提交至CLUST 2015 官方，根據(jù)官方返回跟蹤誤差大小判斷本方法的跟蹤精度，并進行誤差分析與可行性分析。

1.2 算法流程

在CLUST 2015 競賽中，跟蹤不同類型的超聲任務(wù)中，參數(shù)是固定的或者自動調(diào)節(jié)的，本文設(shè)定所有參數(shù)都是固定的。由于CLUST 2015 所給出的感興趣點的坐標(biāo)為亞像素精度，所以首先選用中值流算法，因為其可以跟蹤亞像素精度的點。中值流作為短時跟蹤算法，在長時間跟蹤過程中會發(fā)生漂移并且積累誤差，所以需要進行校正。在跟蹤過程中，如果中值流算法成功跟蹤到目標(biāo)，那就在目標(biāo)附近進行模板匹配，通過模板匹配來校正中值流算法跟蹤過程的漂移；如果中值流算法未能跟蹤到目標(biāo)，則需使用ORB[20]（oriented fast and rotated brief）特征點進行匹配，判斷ROI 的位移，重新選定ROI 繼續(xù)進行跟蹤，如果連續(xù)跟蹤失敗，則說明中值流算法不適用該情況下的跟蹤，需使用KCF[21]（kernelized correlation filters）算法代替中值流算法進行跟蹤。跟蹤算法流程圖如圖1 所示。

圖1 跟蹤算法流程圖

1.2.1 中值流算法

中值流算法是基于傳統(tǒng)光流法的改進[22-23]。光流法是利用圖像序列中像素的變化尋找前一幀圖像和當(dāng)前幀圖像間的對應(yīng)關(guān)系，進而得到兩幀圖像間物體運動狀態(tài)的一種方法。光流法具有2 個必要的前提假設(shè)：

（1）同一個物體在圖像中對應(yīng)的像素亮度不變。由于光流法是根據(jù)像素亮度尋找兩幀圖像中目標(biāo)的運動關(guān)系，因此如果像素亮度發(fā)生改變，那么將無法在兩幀圖像中實現(xiàn)同一個像素或物體的匹配。

（2）要求兩幀圖像必須具有較小的運動。光流法只在原圖像點附近搜索對應(yīng)的像素點，因此兩幀圖像像素位置不能有較大的移動。

假設(shè)在時間t，點（x，y）處的像素灰度值為I（x，y，t），在時間t+Δt，位置（x，y）移動到位置（x+Δx，y+Δy），此時的像素灰度值為I（x+Δx，y+Δy，t+Δt）。根據(jù)灰度一致假設(shè)，即像素的灰度在瞬時運動中應(yīng)該保持一致，因此有

將公式（1）右側(cè)項泰勒一階展開，則有：

忽略高階無窮小量，兩邊同時除以dt，化簡為

稀疏光流法（Lucas-Kanade，L-K）是假設(shè)領(lǐng)域內(nèi)其他像素的位移矢量和中心像素相同，然后通過加權(quán)最小二乘法對領(lǐng)域內(nèi)的所有像素點進行計算。而中值流算法是基于稀疏光流法來進行改進的，稀疏光流法只計算一幅圖像中ROI 部分特征點的光流強度。

中值流算法是通過度量前向、后向之間的誤差來完成自我跟蹤精度的提升，即在兩幀相鄰圖像it、it+1中，先在圖像it中選取Harris 角點Hi，通過稀疏光流算法計算出在圖像it+1對應(yīng)的角點集Hi+1，再對Hi+1反向使用稀疏光流算法計算出對應(yīng)的角點集Hi'，如果角點集Hi和Hi'中對應(yīng)角點位移差距過大則濾去，如圖2 所示。左圖中選取2 個Harris 角點（點1、點2）進行前向-反向光流計算。點1 在左右2 張圖中均可成功跟蹤；點2 在右圖中被遮擋導(dǎo)致跟蹤錯誤，反向跟蹤并沒有返回到點2，而是反向跟蹤到了另一個點3，點3 與點2 間距離過大，說明跟蹤過程發(fā)生了漂移。因此點2 在中值流算法跟蹤過程中被濾去，從而解決遮擋、漂移的問題。

圖2 前向-反向追蹤圖

1.2.2 模板匹配

模板匹配是一種最原始、最基本的模式識別方法，主要研究某一特定對象物的圖案位于圖像的什么地方，進而識別對象物，這就是一個匹配問題，是圖像處理中最基本、最常用的匹配方法。模板就是一幅已知的小圖像，而模板匹配就是在一幅大圖像中搜尋目標(biāo)，已知該圖中有要找的目標(biāo)，且該目標(biāo)與模板有相同的尺寸、方向和圖像元素，通過一定的算法可以在圖中找到目標(biāo)，確定其坐標(biāo)位置。

模板匹配的度量方法是歸一化相關(guān)系數(shù)匹配算法。其計算公式為

式中，w、h 分別代表模板的長和寬；I（x，y）代表輸入大圖像中（x，y）處的像素值大??；T（x'，y'）代表模板小圖像（x'，y'）處的像素值；R（x，y）代表輸入的大圖像中（x，y）處的相似性度量結(jié)果。R（x，y）的值為[-1，1]，-1 代表負(fù)相關(guān)，0 代表不相關(guān)，1 代表正相關(guān)。對于本實驗而言，計算結(jié)果越接近1 代表與原始模板越匹配，效果越好。

對于在超聲邊界附近的感興趣點，判斷模板匹配的搜索范圍是否全部在超聲圖像內(nèi)部，如果不在則不能使用模板匹配校正，因為要防止超聲圖像外部的信息對跟蹤過程造成干擾。除此之外，在跟蹤過程中，模板不能更新，應(yīng)一直使用初始模板，防止跟蹤過程中發(fā)生漂移。

1.2.3 ORB 特征點的檢測和匹配

ORB 特征點是因已有特征點檢測算法[如SIFT（scale invariant feature transform）[24]]，不能滿足需求而開發(fā)出來的。SIFT[24]特征點在1999 年首次提出，并在2004 年得到進一步完善。SIFT 特征點在光照、噪聲、縮放和旋轉(zhuǎn)等干擾下仍具有良好的穩(wěn)定性，在對象識別、圖像拼接等領(lǐng)域應(yīng)用出色。但是由于其運算速度過慢，無法應(yīng)用到實時系統(tǒng)中，最主要的是它受到專利保護，學(xué)者不能隨意使用。后續(xù)學(xué)者對SIFT算法進行了改進和優(yōu)化，出現(xiàn)了如SURF[25]（speeded up robust features）、ORB 等更優(yōu)的算法，其中ORB 特征點具有和SIFT 相似的匹配性能，受圖像噪聲影響較小，并且能夠?qū)崟r應(yīng)用。ORB 特征點由FAST[26]角點和BRIEF[27]描述子組成，首先通過FAST 角點確定圖像中與周圍像素存在明顯差異的像素點作為角點，之后計算每個角點的BRIEF 描述子，從而確定唯一的ORB 特征點。FAST 角點的核心思想是如果在某個灰度值較小的區(qū)域存在一個灰度值明顯變大的像素點，那么該像素點在這個區(qū)域中具有明顯的特征，可以作為特征點。FAST 角點計算過程如下：

第一步：選擇某個像素點作為中心點p，其像素值為Ip。

第二步：設(shè)置判定FAST 角點的像素閾值，例如Tp=0.2×Ip。

第三步：比較中心點p 的像素值與半徑為3 的圓周上所有像素的像素值，如果存在N 個像素的像素值大于Ip+Tp或者小于Ip-Tp，那么中心點p 為FAST角點。

第四步：遍歷所有像素點，重復(fù)上述步驟，計算圖中的FAST 角點。

FAST 角點的選擇如圖3 所示。

圖3 FAST 角點（綠色）與周圍像素的灰度值比對

BRIEF 描述子用于描述特征點周圍像素灰度值的變化趨勢，如果2 個圖中具有相同的描述子，那么認(rèn)為2 個特征點是同一個特征點。ORB 特征點匹配是通過計算2 個描述子之間的漢明距離來判斷。圖4 為CLUST 2015 數(shù)據(jù)集中相鄰兩幀ORB 特征點的匹配結(jié)果。

圖4 ORB 特征點匹配示意圖

1.2.4 KCF 算法

KCF 跟蹤算法是一種核相關(guān)濾波算法，其核心思想就是擴充負(fù)樣本數(shù)量以增強跟蹤器的性能，而擴充負(fù)樣本的方法就是采用循環(huán)矩陣的方法進行構(gòu)造，只有基樣本為正樣本，其余均為負(fù)樣本。

該算法采用嶺回歸的方法訓(xùn)練跟蹤器。訓(xùn)練的目標(biāo)是找到一個函數(shù)f（z）=wTz 最小化樣本xi和回歸目標(biāo)yi的平方誤差，其中λ 為正則化參數(shù)，w 的計算公式如下：

式中，X 為循環(huán)矩陣；I 為單位矩陣；y 為回歸矩陣。

循環(huán)矩陣X 如公式（6）所示：

然后對循環(huán)矩陣X 進行離散傅里葉變換對角化，可變?yōu)?/p>

1.2.5 算法融合

對于感興趣點的跟蹤是基于其周圍區(qū)域的跟蹤，在初始幀以感興趣點的中點選取一個40×40 像素大小的正方形，將該區(qū)域作為參考模板和ROI，對該區(qū)域使用中值流算法跟蹤。如果跟蹤成功，在后續(xù)圖像所跟蹤到的ROI 周圍±2 像素進行模板匹配，如果兩者跟蹤到的區(qū)域中點的歐氏距離大于3 像素，則以中值流算法作為結(jié)果，否則以模板匹配作為結(jié)果，且以模板匹配結(jié)果重新初始化中值流算法ROI；如果中值流算法跟蹤失敗，則用ORB 特征點匹配的方法校正，如果ORB 特征點匹配計算出對應(yīng)特征點的位移大于7 則去除，通過求取ORB 特征點平均位移計算出ROI 的位移大小，重新對ROI 區(qū)域初始化。在使用ORB 特征點校正中值流算法跟蹤失敗時，如果ROI 周圍沒有超過10 個ORB 特征點，那么只能使用模板匹配來跟蹤。如果連續(xù)超過30 幀圖像都不能使用中值流算法成功跟蹤以及ORB 特征點匹配矯正失敗，這種情況代表該任務(wù)無法使用中值流算法來完成，那么則使用KCF 跟蹤。同樣地，使用模板匹配進行長時目標(biāo)跟蹤誤差校正。由于呼吸運動具有周期性，可以通過歷史運動軌跡校正跟蹤錯誤。記錄每連續(xù)200 幀圖像的跟蹤結(jié)果，對于后續(xù)所跟蹤到的感興趣點判斷其是否在運動范圍±3 像素，如果沒有超出范圍，則認(rèn)為正確，否則在運動范圍內(nèi)進行模板匹配，重新初始化ROI。

2 實驗結(jié)果

對測試集中某一感興趣點進行測試，分別使用沒有校正的中值流算法和校正后的算法進行比對，其結(jié)果如圖5 所示。從圖5 可以看出，校正后的誤差明顯減小，證明校正使算法跟蹤精確度明顯改善。

圖5 算法校正前后對比

通過CLUST 2015 競賽官方的評判，本方法的跟蹤誤差為（1.29±1.43）mm（95%置信區(qū)間：3.59 mm），詳見表1。

表1 CLUST 2015 各系列數(shù)據(jù)跟蹤誤差統(tǒng)計 單位：mm

由于所有任務(wù)的ROI 大小都是40×40 像素，所以跟蹤速率相對穩(wěn)定，平均幀率大于22 幀/s，基本滿足實時要求。在某些情況下，幀率可能會突然降低，可能是由于感興趣點突然出現(xiàn)大位移或者需要額外計算，例如重置位置、對ORB 特征點進行匹配等。

在所有跟蹤任務(wù)中存在部分任務(wù)跟蹤誤差較大的情況，具體見表2。

表2 跟蹤任務(wù)誤差較大的情況統(tǒng)計 單位：mm

造成這些任務(wù)誤差較大的原因可能有很多，具體如下：MED-07-3_1 誤差大的可能原因是該點代表的血管過大，大于中值流算法40×40 像素的ROI，導(dǎo)致中值流算法不能很好地跟蹤該區(qū)域，以及某些情況下會出現(xiàn)大位移也同樣導(dǎo)致中值流算法產(chǎn)生誤差；MED-11_2 誤差大的原因是感興趣點十分靠近超聲圖像邊界，導(dǎo)致在跟蹤過程中中值流算法和模板匹配應(yīng)用效果較差；CIL-03_2 可能是由于部分遮擋導(dǎo)致目標(biāo)被遮擋時跟蹤誤差較大；ETH-13-1_1可能是由于圖像對比度較低以及感興趣點所代表血管過小導(dǎo)致的；MED-06-1_4 也可能同樣是因為感興趣點所代表的血管過小導(dǎo)致；MED-07-4_1 與MED-07-3_1 導(dǎo)致誤差過大的原因也應(yīng)該是相同的，都是感興趣點所代表血管過大導(dǎo)致的；MED-11_1 和MED-12_2 可能因為感興趣點周圍強度受到干擾，導(dǎo)致跟蹤效果較差。根據(jù)所有跟蹤結(jié)果的平均誤差可知，圖像對比度越高，處于超聲圖像中間、無大位移和遮擋的情況下跟蹤精度越高。

3 討論

使用手術(shù)導(dǎo)航系統(tǒng)進行肝腫瘤射頻消融時，呼吸運動會引起肝臟的運動和偏移，進而導(dǎo)致誤差。為了解決這一難題，本文方法主要基于中值流算法、模板匹配、ORB 特征點匹配和KCF 組合對二維超聲圖像中的感興趣點進行跟蹤。CLUST 2015 數(shù)據(jù)集由24個訓(xùn)練數(shù)據(jù)集和39 個測試數(shù)據(jù)集組成，分別具有53 個和85 個帶標(biāo)注的目標(biāo)。本方法平均跟蹤誤差為（1.29±1.43）mm（95%置信區(qū)間：3.59 mm），具有一定的魯棒性，快速準(zhǔn)確，且不需要對每位患者進行訓(xùn)練。由于算法的運行時間主要由ROI 大小確定，且ROI 大小為固定值（40×40 像素），所以本方法的運行時間很穩(wěn)定，平均幀率大于22 幀/s。Bharadwaj 等[14]使用孿生神經(jīng)網(wǎng)絡(luò)和模板匹配在超聲圖像中跟蹤目標(biāo)，添加模板匹配模塊和線性卡爾曼濾波模塊校正跟蹤過程中的漂移，其平均跟蹤誤差為1.60 mm，缺點是在感興趣點接近超聲圖像邊界、遮擋以及對比度不高的情況下容易跟蹤失敗、產(chǎn)生較大誤差。Hallack 等[13]使用logDemons 非線性配準(zhǔn)以及使用SIFT 特征點匹配校正跟蹤過程中的漂移和誤差，其平均跟蹤誤差為1.21 mm，缺點是SIFT 特征點計算過程很慢，更適合離線操作，算法運算速度較慢，平均幀率為12 幀/s，且感興趣點在相鄰兩幀中位移大小與其所選中ROI 大小近似時容易失敗。本文使用中值流算法以及模板匹配和ORB 特征點匹配進行超聲肝臟的跟蹤，跟蹤過程中的漂移和誤差可以較好地跟蹤肝臟運動，且ORB 特征點匹配具有很好的實時性，所以時間復(fù)雜度低、計算速度快，可以滿足實時需求。

對于后續(xù)的研究，一是研究一種可以自動選擇ROI 大小的跟蹤算法，使得ROI 內(nèi)所攜帶的血管信息和血管周圍的信息比例保持在一個適當(dāng)?shù)姆秶?，從而解決血管過大或過小導(dǎo)致跟蹤誤差變大的問題；二是針對遮擋、感興趣點在超聲邊界以及對比度較差的情況，使用supports[28]方法來解決，該方法可以通過關(guān)聯(lián)點計算目標(biāo)點的位置；三是對于時間延遲問題[29-31]的補償，由于數(shù)據(jù)傳輸、計算以及消融針進入人體內(nèi)均需要一定的時間，如果不考慮這段進針時間或不進行時間補償，消融針到達(dá)腫瘤上的目標(biāo)點時會導(dǎo)致消融針進針位置錯誤。因此，需要對該時間延遲進行補償，可以使用支持向量回歸、卡爾曼濾波等方法來解決。

4 結(jié)論

本文設(shè)計了一種基于中值流的方法在超聲圖像中跟蹤肝臟運動，包括跟蹤、漂移校正、跟蹤失敗的ROI 重置。通過提交本文所提出的跟蹤方法的運行結(jié)果至CLUST 2015，驗證了本方法的跟蹤準(zhǔn)確性。本方法跟蹤精度較高、運行速度較快，可以滿足熱消融手術(shù)中對腫瘤跟蹤的臨床需求。