一種基于HAMMER算法的軟組織形變檢測(cè)新方法*

馬昂昂 廖琪梅 趙 菲② 盧虹冰

近年來，計(jì)算機(jī)輔助外科(CAS)技術(shù)廣泛應(yīng)用于神經(jīng)外科、整形外科和耳鼻喉科手術(shù)。在這些手術(shù)領(lǐng)域，導(dǎo)航技術(shù)可以有效提高手術(shù)精確度，確保患者手術(shù)安全[1-2]。目前，也有研究小組在胸腹部手術(shù)中引入導(dǎo)航技術(shù)[3]。然而，由于受到呼吸、心跳、患者移動(dòng)和手術(shù)操作的影響，在這類手術(shù)中，術(shù)前圖像無法精確反映術(shù)中解剖結(jié)構(gòu)的變化。雖然通過對(duì)術(shù)前模型及其性質(zhì)的快速數(shù)學(xué)計(jì)算可以在一定程度上反映術(shù)中的形變，但這類工具還未能有效地應(yīng)用于臨床[4]。因此，必須獲得術(shù)中影像數(shù)據(jù)以反映器官形變。

超聲成像具有實(shí)時(shí)、無創(chuàng)的優(yōu)點(diǎn)，因而廣泛應(yīng)用于臨床診斷中[5]。但由于超聲圖像的分辨率和掃描區(qū)域因斑點(diǎn)噪聲影響而受到局限，若能將術(shù)前器官組織高分辨率、高對(duì)比度的三維MRI圖像和術(shù)中超聲圖像進(jìn)行融合，可能是實(shí)現(xiàn)術(shù)中形變器官跟蹤和軟組織形變檢測(cè)的有效方法[6]。雖然三維超聲成像已成功應(yīng)用于術(shù)中的實(shí)時(shí)導(dǎo)航，但目前相關(guān)研究主要集中在形變器官的跟蹤和可視化方面，基于多模態(tài)影像的軟組織形變檢測(cè)和校正的研究相對(duì)開展較少[3]。

除了術(shù)中成像，另一個(gè)實(shí)現(xiàn)形變校正的關(guān)鍵問題是如何從多模態(tài)圖像中檢測(cè)到形變的器官，并對(duì)其產(chǎn)生的形變進(jìn)行有效校正[7]。其中，作為一種高分辨率、基于體素的非剛性形態(tài)學(xué)配準(zhǔn)方法，HAMMER算法利用從圖像中派生出的屬性向量，基于層次形變機(jī)制，實(shí)現(xiàn)了不同個(gè)體腦MRI圖像的準(zhǔn)確融合，具有很強(qiáng)的解剖結(jié)構(gòu)檢測(cè)和配準(zhǔn)能力[8]。

本研究中，為了探索融合術(shù)前三維MRI圖像和術(shù)中二維超聲圖像，進(jìn)行腹部軟組織器官形變檢測(cè)和校正的可行性，通過利用微型電磁導(dǎo)航系統(tǒng)，結(jié)合改進(jìn)的HAMMER算法，實(shí)現(xiàn)了一種形變檢測(cè)和校正的新方法。該方法通過將MRI圖像和超聲圖像進(jìn)行配準(zhǔn)和融合，對(duì)腹部軟組織器官的形變進(jìn)行檢測(cè)和校正。為了定量評(píng)價(jià)所提出的方法，建立了超聲圖像計(jì)算機(jī)仿真模塊，通過對(duì)腹部MRI圖像進(jìn)行剛性和非剛性形變，并利用形變后的MRI圖像產(chǎn)生模擬的超聲圖像。在此基礎(chǔ)上，結(jié)合電磁導(dǎo)航定位數(shù)據(jù)和改進(jìn)的HAMMER算法，對(duì)形變的超聲圖像和原始MRI圖像進(jìn)行配準(zhǔn)，同時(shí)采用平均位移和互信息對(duì)形變檢測(cè)和校正的結(jié)果進(jìn)行定量評(píng)價(jià)[9]。

1 方法和實(shí)驗(yàn)原理

1.1 軟組織形變實(shí)時(shí)跟蹤系統(tǒng)

在課題的前期研究中，建立軟組織形變實(shí)時(shí)跟蹤系統(tǒng)的框架和方法[10]。在這套系統(tǒng)中，通過將高分辨率的術(shù)前CT和(或)MRI圖像、實(shí)時(shí)超聲圖像與電磁定位系統(tǒng)有機(jī)結(jié)合，初步實(shí)現(xiàn)了基于超聲的軟組織形變實(shí)時(shí)跟蹤。基于該框架，可以對(duì)感興趣區(qū)域軟組織的三維結(jié)構(gòu)和形變進(jìn)行有效跟蹤，從而提高腹部介入治療、微創(chuàng)手術(shù)以及開放手術(shù)的精度。

1.2 超聲圖像的計(jì)算機(jī)仿真

在上述跟蹤系統(tǒng)中，主要通過將二維超聲探頭與微型電磁定位系統(tǒng)結(jié)合，對(duì)探頭及二維超聲成像的位置進(jìn)行定位，并注冊(cè)到患者所在物理空間，實(shí)現(xiàn)術(shù)中成像的獲取和定位。然而，由于術(shù)中人體的形變多樣，獲取數(shù)據(jù)的條件及環(huán)境因素復(fù)雜，如果直接利用術(shù)中數(shù)據(jù)，將會(huì)因缺乏真實(shí)的標(biāo)準(zhǔn)而很難對(duì)提出的方法進(jìn)行有效評(píng)價(jià)[11]。為此，開發(fā)了超聲圖像計(jì)算機(jī)仿真模塊，其目的是基于形變后的MRI圖像而產(chǎn)生模擬的超聲圖像。

為了提高仿真圖像的真實(shí)性，在超聲圖像計(jì)算機(jī)仿真過程中，需要建立兩個(gè)重要的模型：超聲采樣模型和噪聲模型[12]。其中，超聲采樣模型通過對(duì)超聲束的發(fā)射以及回波采集過程的分析來模擬真實(shí)超聲圖像的獲取過程。斑點(diǎn)噪聲模型則是通過對(duì)圖像施加符合一定統(tǒng)計(jì)規(guī)律的隨機(jī)噪聲來生成最終的超聲圖像。超聲圖像計(jì)算機(jī)仿真的流程如圖1所示。

圖1 超聲圖像計(jì)算機(jī)仿真的流程

圖2 超聲圖像采樣過程

1.2.1 基于目標(biāo)圖像的超聲采樣

從目標(biāo)圖像(本實(shí)驗(yàn)中為MRI圖像)中生成超聲圖像，其過程可通過計(jì)算發(fā)射超聲波束與目標(biāo)圖像的交點(diǎn)，并進(jìn)一步估計(jì)不同深度的超聲回波值來模擬得到。整個(gè)采樣算法的參數(shù)共有7個(gè)：①超聲束數(shù)目n；②沿每條超聲束方向上的采樣點(diǎn)數(shù)目m；③掃描角度Θ；④圖像寬度ω；⑤扇區(qū)原點(diǎn)高度y0；⑥軸向距離dmin；⑦軸向距離dmax。通過合理設(shè)置這些參數(shù)，可以模擬出具有真實(shí)感的超聲圖像。上述參數(shù)的定義和采樣過程如圖2所示。

從圖2可以得出目標(biāo)圖像的笛卡爾坐標(biāo)系(x,y)和超聲采樣點(diǎn)的極坐標(biāo)系(ρ,θ)間的關(guān)系公式1：

1.2.2 插值

采樣完成后需要對(duì)新圖像中的空白點(diǎn)進(jìn)行插值。考慮到超聲圖像的反射特性，本研究采用sin c函數(shù)作為卷積核[13]，如公式2所示，每個(gè)空白點(diǎn)的值由臨近像素的卷積獲得，臨近像素對(duì)采樣點(diǎn)的貢獻(xiàn)由它們距采樣點(diǎn)的距離x決定：

1.2.3 加噪

超聲圖像主要受乘性噪聲影響，在乘性噪聲中，斑點(diǎn)噪聲被認(rèn)為是超聲圖像退化的主要原因[14]。本研究中采用下面的模型公式3反映噪聲對(duì)超聲圖像的影響：

其中，I(x,y)為原始未污染的圖像，在本文中即為經(jīng)過采樣和插值后所得到的圖像。I^(x,y)為噪聲污染后的圖像。ηm為斑點(diǎn)噪聲，服從瑞利分布(見公式4)：

1.3 改進(jìn)的HAMMER算法

Ding[15]提出的HAMMER算法是一種綜合考慮圖像灰度和幾何結(jié)構(gòu)特征的彈性配準(zhǔn)算法。相較于其他配準(zhǔn)算法，HAMMER的創(chuàng)新主要體現(xiàn)在兩個(gè)方面：

(1)HAMMER算法在每個(gè)體素點(diǎn)上定義一個(gè)屬性向量，以反映底層解剖結(jié)構(gòu)的幾何特性。在形變過程中，屬性向量通過外形的相似性建立了解剖結(jié)構(gòu)的一致性。每個(gè)屬性向量a(x)由三部分組成：①邊界類型；②灰度強(qiáng)度；③幾何不變矩(GMIs)。

(2)HAMMER算法的層次形變機(jī)制。部分解剖結(jié)構(gòu)有特定的形狀特征，這些形狀特征在初始化形變階段被稱為"錨點(diǎn)"，剩下的部分緩慢影響形變從而避免錯(cuò)誤的目標(biāo)。此外，在形變子空間時(shí)，形變機(jī)制形變圖像相關(guān)的大部分不是單個(gè)體素，然后在整個(gè)子空間中估計(jì)屬性向量的相似性，使得形變更加真實(shí)可信。有關(guān)于HAMMER算法更多的細(xì)節(jié)描述可以參見文獻(xiàn)[4][8]。

考慮到軟組織形變的非剛性以及多模圖像強(qiáng)度分布的不同，在軟組織的形變跟蹤中，基于邊緣特性的彈性配準(zhǔn)方法應(yīng)該更加適用[16]。有鑒于此，本研究選擇HAMMER算法，用于軟組織器官的形變檢測(cè)和校正。

原HAMMER算法主要用于腦MRI單模態(tài)圖像的配準(zhǔn)，為了使HAMMER算法更適合MRI和腹部超聲圖像的配準(zhǔn)，對(duì)算法做了相應(yīng)的優(yōu)化和改進(jìn)。特別是每個(gè)體素的屬性向量不再局限于GMIs，而是由一系列選定的屬性組成，這些屬性可從鄰近像素計(jì)算得到。本實(shí)驗(yàn)中增加了像素強(qiáng)度、梯度幅值以及梯度幅值的方差作為屬性向量。

2 實(shí)驗(yàn)和結(jié)果

2.1 MRI圖像數(shù)據(jù)的采集

本實(shí)驗(yàn)的MRI圖像源自第四軍醫(yī)大學(xué)附屬西京醫(yī)院放射科，采用3.0 T MR(MAGNETOM Trio, Siemens AG, Erlangen, Germany)掃描獲得。掃描參數(shù)設(shè)置為：TR=250 ms，TE=2.3 ms，圖像像素尺寸為256×256, 視野為220 mm, 層厚為1.00 mm。

2.2 MRI圖像形變和超聲圖像仿真

為了模擬腹部不同形變的效果，對(duì)腹部肝臟MRI圖像進(jìn)行了如下3種幾何變換，即旋轉(zhuǎn)變換(逆時(shí)針旋轉(zhuǎn)15o)，仿射變換(變換矩陣)和投影變換(變換矩陣)，其結(jié)果如圖3所示：

圖3 原始和形變的MRI圖像

基于形變的MRI圖像，通過超聲圖像計(jì)算機(jī)仿真，得到相應(yīng)的模擬超聲圖像，如圖4所示：

圖4 模擬超聲圖像

為了實(shí)現(xiàn)對(duì)軟組織形變的檢測(cè)與校正，使用優(yōu)化的HAMMER算法，分別對(duì)形變后的MRI圖像和原始MRI圖像以及模擬的形變超聲圖像和原始MRI圖像進(jìn)行配準(zhǔn)，其結(jié)果如圖5所示。

圖5 配準(zhǔn)結(jié)果

需要注意的是，HAMMER算法中需要確定一些參數(shù)，但其中最基本的參數(shù)就是整體搜索大小(δ)和當(dāng)前迭代順序(i)，其他參數(shù)則可進(jìn)行相應(yīng)估計(jì)。本實(shí)驗(yàn)中，依據(jù)經(jīng)驗(yàn)，δ設(shè)置為9，i設(shè)置為25。圖5示為經(jīng)過相應(yīng)校正后的MRI圖像和超聲圖像。

(a1)-(a2)對(duì)旋轉(zhuǎn)變換進(jìn)行校正后的MRI和超聲圖像，(b1)-(b2)對(duì)仿射變換進(jìn)行校正后的MRI和超聲圖像，(c1)-(c2)對(duì)投影變換進(jìn)行校正后的MRI和超聲圖像為了定量評(píng)價(jià)形變檢測(cè)和校正的效果，使用文獻(xiàn)[17]中關(guān)于平均位移的計(jì)算方法分別計(jì)算不同變換條件下，校正前及校正后圖像的平均位移，以定量評(píng)價(jià)形變校正的效果。其結(jié)果可以看出，經(jīng)過配準(zhǔn)和校正后，雖然對(duì)MRI單模圖像的校正效果整體優(yōu)于對(duì)MRI和US多模圖像的校正，但是在3種變換下，校正后超聲圖像的平均位移均＜4個(gè)像素，圖像形變得到很好的補(bǔ)償(見表1)。

表1 形變圖像和校正后圖像的整體平均位移

由于系統(tǒng)的檢測(cè)和校正過程涉及兩種成像模態(tài)，為了進(jìn)一步評(píng)價(jià)圖像整體的校正效果，計(jì)算了原始MRI圖像和校正后MRI圖像、超聲圖像的互信息值(見表2)。

表2 原始圖像和校正后圖像的互信息值

3 結(jié)語

目前，在手術(shù)中軟組織形變的檢測(cè)和校正已經(jīng)成為軟組織導(dǎo)航技術(shù)的熱點(diǎn)和難點(diǎn)問題。在研究中，提出一種新的軟組織形變檢測(cè)和校正方法，以改善跟蹤和校正的效果。提出的新方法將術(shù)前MRI高分辨率圖像和術(shù)中實(shí)時(shí)超聲圖像通過微型電磁定位系統(tǒng)有機(jī)結(jié)合，并利用改進(jìn)的HAMMER算法對(duì)形變器官軟組織進(jìn)行檢測(cè)和校正。基于平均位移和互信息的定量評(píng)價(jià)結(jié)果表明，所提出的新方法具有可行性。

在后續(xù)工作中，進(jìn)一步擬將提出的新方法擴(kuò)展至三維圖像，并利用臨床實(shí)時(shí)超聲掃描數(shù)據(jù)和定位信息，對(duì)提出新方法的臨床可行性進(jìn)行驗(yàn)證。此外，提出的新方法還將與導(dǎo)航系統(tǒng)硬件相結(jié)合，將進(jìn)一步開發(fā)軟組織導(dǎo)航系統(tǒng)。

[1]Xing L,Thorndyke B,Schreibmann E,et al. Overview of image-guided radiation therapy[J]. Med Dosim,2006,31(2):91-112.

[2]Pluim JP, Fitzpatrick JM. Image registration[J]. IEEE Trans Med Imaging,2003,22(11):1341-1343.

[3]Wein W,Brunke S,Khamene A,et al.Automatic CT-ultrasound registration for diagnostic imaging and image-guided intervention[J].Med Image Anal,2008,12(5):577-585.

[4]Shevchenko N, Schwaiger J, Markert M, et al. Evaluation of a resectable ultrasound liver phantom for testing of surgical navigation systems[J].Conf Proc IEEE Eng Med Biol Soc,2011,2011:916-919.

[5]Siddiq A, El Sayed T.Ultrasonic-assisted manufacturing processes: variational model and numerical simulations[J]. Ultrasonics,2011,52(4):521-529.

[6]Coenegrachts K. Magnetic resonance imaging of the liver: New imaging strategies for evaluating focal liver lesions[J].World J Radiol,2009,1(1):72-85.

[7]Xu HX,Lu MD,Liu LN,et al.Magnetic n a v i g a t i o n i n u l t r a s o u n d-g u i d e d interventional radiology procedures[J].Clin Radiol,2011,67(5):447-454.

[8]Ding GS,Christos D,HAMMER:hierarchical attribute matching mechanism for elastic registration[J].IEEE Transactions on Medical Imaging,2001,21:1421-1439.

[9]羅建文,范宇,白凈.基于位移量的軟組織應(yīng)變測(cè)量快速算法[J].聲學(xué)學(xué)報(bào)(中文版),2002,27(1):18-22.

[10]馬寶秋.基于實(shí)時(shí)超聲影像的軟組織形變跟蹤技術(shù)[D].西安:第四軍醫(yī)大學(xué),2010.

[11]Mazza E, Barbarino GG. 3D mechanical modeling of facial soft tissue for surgery simulation[J].Facial Plast Surg Clin North Am,2011,19(4):623-637.

[12]張偉志,剛鐵,王軍.超聲波檢測(cè)計(jì)算機(jī)模擬和仿真的研究及應(yīng)用現(xiàn)狀[J].應(yīng)用聲學(xué),2003,22(3):39-44.

[13]Huang QH,Zheng YP, Lu MH, et al. Development of a portable 3D ultrasound imaging system for musculoskeletal tissues[J]. Ultrasonics,2005,43(3):153-163.

[14]Hiremath PS, Prema TA, Sharan B. Speckle Reducing Contour let Transform for Medical Ultrasound Images[J].International Journal of Computer and Information Engineering,2010,4(4):284-291.

[15]Ding GS. Image registration by hierarchical matching of local spatial intensity histograms[C].Medical Image Computing and Computer Assisted Intervention,France, 2004:582-590.

[16]Taylor ZA, Comas O, Cheng M, et al. Modelling anisotropic viscoelasticity for real-time soft tissue simulation[J]. Med Image Comput Comput Assist Interv,2008,11(1):703-710.

[17]Robert H, Zdenko B, Martin H. Methodology for determination of modal parameters by digital image correlation[C].Modelling of mechanical and mechatronic systems 2011,Slovak Republic,2011:193-201.