馮兆美，黨 軍，崔崤峣，焦 陽

（1．中國科學(xué)院 蘇州生物醫(yī)學(xué)工程技術(shù)研究所，江蘇 蘇州215163；2．中國科學(xué)院 長春光學(xué)精密機械與物理研究所，吉林 長春130000；3．中國科學(xué)院大學(xué)，北京100049）

醫(yī)學(xué)圖像配準(zhǔn)是指采用適當(dāng)?shù)目臻g變換方法，使一幅醫(yī)學(xué)圖像與另一幅醫(yī)學(xué)圖像上對應(yīng)的解剖結(jié)構(gòu)點、具有診斷意義的點以及手術(shù)感興趣的區(qū)域達(dá)到空間位置上的一致［1］。多模態(tài)醫(yī)學(xué)圖像配準(zhǔn)在醫(yī)學(xué)診斷和治療計劃中有著重要作用。但由于成像方法以及病人器官解剖特征結(jié)構(gòu)的不同，使得準(zhǔn)確、高效的對多模態(tài)醫(yī)學(xué)圖像進(jìn)行非剛性配準(zhǔn)成為一個難題［2］。

例如超聲成像（Ultrasound，US）因具有實時成像、低成本和無輻射等優(yōu)點在臨床領(lǐng)域被廣泛應(yīng)用，特別是用于引導(dǎo)介入和術(shù)中監(jiān)測等方面。CT圖像由于其根據(jù)人體各種組織或器官對X射線吸收不等的特性進(jìn)行掃描，所以對人體各個器官呈現(xiàn)斷層切面且圖像分辨率較高，但CT成像更適合進(jìn)行手術(shù)前的診斷。由于超聲圖像能夠?qū)崟r快速的提供診斷數(shù)據(jù)，CT圖像能夠保證數(shù)據(jù)的準(zhǔn)確詳細(xì)，所以進(jìn)行超聲圖像和CT圖像的配準(zhǔn)研究在臨床診斷效率、實時病情監(jiān)測、提高外科手術(shù)水平等方面具有重要的臨床應(yīng)用價值。近些年，三維圖像非剛性配準(zhǔn)算法逐漸成為專家學(xué)者研究的重心，Wein Wolfgang［3］等人采用剛性與放射變換模型，對超聲圖像與CT進(jìn)行自動化配準(zhǔn)，并對此方法進(jìn)行評估；盧振態(tài)［4］等人提出一種新的基于主成分分析的三維醫(yī)學(xué)圖像快速配準(zhǔn)算法，此方法能夠顯著減少計算量，并且不需要復(fù)雜的優(yōu)化過程；Rueckert Daniel［5］等人提出利用自由形式變換模型對三維乳腺MRI圖像進(jìn)行非剛性配準(zhǔn)；Roche Alexis和Pennec Xavier［6］等提出將基于圖像梯度與灰度的混合配準(zhǔn)模型應(yīng)用于乳腺MRI圖像的非剛性配準(zhǔn)中。不斷有基于新的數(shù)學(xué)模型的算法出現(xiàn)，如 Tanner Christine［7］等人建立了錯切變換與橢球壓縮模型描述乳腺形變；Mark Holden［8］使用微分同胚時空映射來正則化速度場，避免求解Navier－Stokes方程時遇到的奇異點問題。

本文采用加入薄板樣條能量約束項的互信息［9］作為相似性測度對3－D超聲圖像和3－D CT圖像進(jìn)行配準(zhǔn)，聯(lián)合剛性變換和自由形變模型，并使用有限內(nèi)存的擬牛頓優(yōu)化方法（L－BFGS－B）［10］對變形參數(shù)進(jìn)行優(yōu)化搜索。采用基于體素的方法計算相似度時，對于較高分辨率的三維體數(shù)據(jù)集，其包含的數(shù)據(jù)量極大，為減少計算量，計算前先對原圖像進(jìn)行重采樣。同時在非剛性配準(zhǔn)之前，首先對兩幅圖像進(jìn)行剛性配準(zhǔn)，以使得在非剛性配準(zhǔn)前兩幅圖像的初始幾何偏差達(dá)到最小。在配準(zhǔn)過程中，采用多分辨率和多重網(wǎng)格的策略，以提高配準(zhǔn)過程的速度和精度，避免局部極值問題的產(chǎn)生。最后使用互信息值、相對重疊率和運行時間對乳腺仿體的超聲數(shù)據(jù)和CT數(shù)據(jù)配準(zhǔn)的精度和速度進(jìn)行評價。此外，還使用了29組人體肝臟超聲數(shù)據(jù)和3組人體胸腔CT圖像進(jìn)一步驗證算法性能。

1 原理與方法

假設(shè)參考圖像（Reference Image）為fR（X），另一幅待配準(zhǔn)的測試圖像（Testing Image）為fT（X），使用g（X｜μ）描述圖像形變，μ是一系列的變換參數(shù)。將醫(yī)學(xué)圖像配準(zhǔn)看成是一種優(yōu)化問題，那么校準(zhǔn)參考圖像和變換后的測試圖像的主要任務(wù)就是尋找一組最優(yōu)的變換參數(shù)μ使兩幅圖像差異S最小，數(shù)學(xué)表示如下：

1．1 變形模型的建立

基于樣條函數(shù)的配準(zhǔn)方法主要利用控制點和樣條函數(shù)來描述非線性幾何變換域。將位于樣條函數(shù)柵格上的控制點賦值，當(dāng)這些控制點產(chǎn)生空間位移時，樣條函數(shù)被彎曲，那么控制點周圍的圖像空間也會隨著樣條函數(shù)的彎曲而產(chǎn)生不同程度的扭曲。本文使用基于B樣條的自由形變模型（free form deformation，F(xiàn)FD）［11］。在三維空間域中，基于B樣條函數(shù)的自由形變可描述為三個一維三次B樣條函數(shù)的張量積：

其中βi（u）為第i次B樣條的基函數(shù)：

同時用一個旋轉(zhuǎn)矩陣R和一個3元素的平移向量T建立剛性的全局變形模型，這樣應(yīng)用于測試圖像的非線性變換為：

其中Xc是測試圖像體素的中心坐標(biāo)。

基于B樣條的自由形變模型是非剛性配準(zhǔn)的重要手段，它能夠有效地擬合物體的彈性形變。FFD變換的自由度很高，可以實現(xiàn)局部控制并且變換域C2連續(xù)性［11］，能夠確保在優(yōu)化過程中使用基于梯度有約束優(yōu)化方法。但是此變形可能產(chǎn)生不可逆的空間變換函數(shù)，即當(dāng)某些B樣條網(wǎng)格中的節(jié)點位移過大時，會導(dǎo)致圖像折疊［12］。為了得到一個可逆的空間變換函數(shù)以避免折疊現(xiàn)象發(fā)生，需要對節(jié)點的位移大小進(jìn)行約束。根據(jù)3次B樣條的可逆性條件［7］，第（i，j，k）個控制點的位移dx、dy、dz應(yīng)滿足條件dx＜sx／k、dy＜sy／k、dz＜sz／k，k是常數(shù)，在三維空間變換的情況下取值為2．47947。由于限制了節(jié)點位移的大小，一層固定網(wǎng)格間距的B樣條模型就不能完成大規(guī)模形變，需要引入多層次B樣條模型［13］。多層次B樣條模型采用多個由疏到密的網(wǎng)格，分層依次對圖像進(jìn)行配準(zhǔn)。

1．2 互信息相似性測度

互信息［14］是一種基于熵的方法，廣泛應(yīng)用于多種不同模式圖像的配準(zhǔn)，互信息方法同時也是配準(zhǔn)精度和魯棒性最好的可回溯性配準(zhǔn)方法之一。假設(shè)兩幅圖像分別用兩個隨機變量A，B表示，其邊緣概率分布為pA（a）與pB（b），聯(lián)合概率分布為pAB（a，b）。那么用互信息來衡量圖像之間的相似程度，其數(shù)學(xué)表達(dá)式為

非剛性配準(zhǔn)中，變換模型內(nèi)部完全保證幾何變換的平滑可逆是一個病態(tài)問題，配準(zhǔn)過程中會出現(xiàn)一定程度的圖像交叉和重疊問題。本文通過在配準(zhǔn)代價函數(shù)中加入一個幾何變換的約束項解決此問題，即在目標(biāo)函數(shù)中加入薄板樣條能量約束項［17］。

那么配準(zhǔn)的代價函數(shù)為：

其中，λ為常量權(quán)重因子，其為經(jīng)驗值。在λ＝0．01時，使得目標(biāo)函數(shù)中相互競爭的兩部分能夠很好的折衷。通過實驗發(fā)現(xiàn)，當(dāng)控制點網(wǎng)格分辨率較低時，λ對B樣條的平滑性影響不大，但是隨著控制點網(wǎng)格分辨率的提高，約束項的作用愈加明顯。原因是FFD對圖像局部的變形能力隨著控制點網(wǎng)格間距的減少而增強。因此，在致密控制點網(wǎng)格中尤其需要約束項。

1．3 多分辨率優(yōu)化策略

本文采用L－BFGS－B優(yōu)化方法，基于B樣條的自由形變模型的互信息配準(zhǔn)算法運算量比較大，而L－BFGS－B［15］優(yōu)化方法能夠解決求解過程占用大量內(nèi)存而降低運算速度的問題。多分辨率優(yōu)化策略［16］亦是一種有效提高算法的執(zhí)行速度和避免陷入局部極值的方法。因此可以采用多分辨率的分級配準(zhǔn)策略，配準(zhǔn)按照由粗到精的方式執(zhí)行。配準(zhǔn)將分為兩個層次進(jìn)行，第一層只包含剛性變換參數(shù)的粗配準(zhǔn)，第二層只包含自由形變參數(shù)的精細(xì)配準(zhǔn)。將剛性變換得到的參數(shù)作為第二層精細(xì)配準(zhǔn)的初始值。配準(zhǔn)框架如圖1。圖中被細(xì)線框住的部分是算法的核心，包括對目標(biāo)函數(shù)的計算和優(yōu)化，L－BFGS－B為尋找空間變換參數(shù)μ的優(yōu)化器。當(dāng)同時滿足終止1的條件（目標(biāo)函數(shù)值與目標(biāo)函數(shù)的梯度值滿足預(yù)設(shè)的閾值，默認(rèn)值分別為1e－3，1e－4）與終止2的條件（迭代次數(shù)，默認(rèn)值為500），輸出配準(zhǔn)圖像、運行時間以及變形參數(shù)。

2 實驗結(jié)果與分析

為了檢驗本文算法對三維超聲圖像與CT圖像配準(zhǔn)的效果，使用乳腺仿體的超聲圖像和CT圖像進(jìn)行非剛性配準(zhǔn)實驗。此外，使用29組人體超聲圖像對和3組人體胸腔CT圖像對進(jìn)一步驗證算法性能（硬件環(huán)境為Intel（R）Core（TM）i3－2 1 2 0CPU＠3．3 0GHz和4．0 0GB RAM ；軟 件環(huán)境為 Microsoft VS2008，ITK 3．20，VolView 3．4，Matlab R2011b，Windows 7）。

圖1 多分辨率圖像配準(zhǔn)框架

圖2 乳腺仿體的超聲－CT圖像

本文使用的是一種乳腺超聲活檢訓(xùn)練用仿體（052A 模型，CIRS，Tissue Simulation ＆ Phantom Technology，USA）采集超聲與CT數(shù)據(jù)。在使用ATL HDI 5000Ultrasound系統(tǒng)掃描超聲數(shù)據(jù)過程中，將乳腺仿體浸入在水槽中，使用二維線陣探頭進(jìn)行垂直無接觸式掃描，掃描步長為2 mm／幀，將獲得的二維超聲數(shù)據(jù)重建成三維體數(shù)據(jù)（圖2b）。在Philips Gemini TF 64PET／CT系統(tǒng)上收集CT圖像（圖2a）。將采集到的超聲圖像與CT圖像進(jìn)行配準(zhǔn)，圖2是乳腺仿體的超聲體數(shù)據(jù)與CT體數(shù)據(jù)，其中CT圖像的大小為240×145×425，像素間距0．35×0．35×0．33；超聲圖像的大小為482×308×30，像素間距為0．35×0．35×0．33。圖中箭頭標(biāo)識的部分是“腫瘤”部分。圖3是其配準(zhǔn)結(jié)果，圖像的大小為482×308×30，像素間距為0．35×0．35×0．33。我們給出其三維體積圖與橫斷面、冠狀面、矢狀面三個方向上的信息。從圖3可以看出，與原始的超聲圖像相比，配準(zhǔn)后的圖像更加光滑，但是腫塊邊緣仍然有一定程度上的模糊，這是因為：1）超聲圖像有噪聲；2）與CT體積圖425個切片相比，超聲體積圖只有30個切片，使得超聲圖像橫截面的分辨率更低，失去更多的結(jié)構(gòu)信息。

圖3 US－CT配準(zhǔn)結(jié)果

收集29組合適的人體肝臟超聲圖像數(shù)據(jù)，在ATL HDI 5000Ultrasound系統(tǒng)上，使用扇形掃描方式得到超聲體數(shù)據(jù)。超聲圖像的分辨率為512×512×256，像素間距為0．427×0．387×0．823。圖4是同一病人不同時間掃描的肝臟超聲圖像，圖4a為參考圖像，圖4b為測試圖像，其中箭頭標(biāo)志為目標(biāo)區(qū)域。圖5為兩幅超聲圖像的配準(zhǔn)結(jié)果，圖像的大小為512×512×256，像素間距為0．427×0．387×0．823。我們給出其三維體積圖與橫斷面、冠狀面、矢狀面三個方向上的信息。從圖中可以看到，箭頭標(biāo)志部分整體上對齊較好。

圖4 肝部超聲臨床圖像

圖5 肝臟超聲圖像配準(zhǔn)結(jié)果

在Philips Gemini TF 64PET／CT 系統(tǒng)上收集到3組合適的人體胸腔CT圖像，CT圖像的分辨率為512×512×41，像素間距為0．703×0．703×0．817。圖6是同一病人不同時間掃描的胸腔CT圖像，其中圖6a為參考圖像，圖6b為測試圖像圖。圖7為其配準(zhǔn)結(jié)果，圖像的大小為512×512×41，像素間距為0．703×0．703×0．817。我們給出其三維體積圖與橫斷面、冠狀面、矢狀面三個方向上的信息。

圖6 胸腔CT圖像

3 方法評價

圖像配準(zhǔn)是為某一特定目的而在某種特定原則下尋求相對最優(yōu)結(jié)果，針對非剛性多模態(tài)圖像配準(zhǔn)，由于獲取圖像的本質(zhì)上有很大程度的不同，使得像素點之間的對應(yīng)關(guān)系是未知的，所以很難

對這些算法進(jìn)行評估，無法指定一個統(tǒng)一的金標(biāo)準(zhǔn)對任意配準(zhǔn)結(jié)果進(jìn)行客觀判定，但是可以使用多種不同標(biāo)準(zhǔn)提供的信息來評估非剛性配準(zhǔn)算法的性能［17，18］。

感興趣區(qū)域的對齊程度［18］能夠很好的體現(xiàn)出兩幅圖像的配準(zhǔn)效果，即兩個對應(yīng)的分割區(qū)域的重疊率。相對重疊率數(shù)學(xué)表達(dá)如式8，其中P與S是相對應(yīng)的分割區(qū)域。

我們請五位有經(jīng)驗的專家分別從體積圖中手動分割出目標(biāo)區(qū)域。圖8給出配準(zhǔn)前后圖像對應(yīng)目標(biāo)區(qū)域的重疊圖像，用白色表示配準(zhǔn)前的目標(biāo)區(qū)域，彩色表示配準(zhǔn)后的目標(biāo)區(qū)域。其中，圖8a是乳腺仿體中的腫塊區(qū)域重疊圖，圖8b是超聲圖像目標(biāo)區(qū)域的重疊圖像，圖8c是從CT胸腔圖像分割出的肝臟重疊圖。配準(zhǔn)前后的目標(biāo)區(qū)域整體上能夠達(dá)到空間位置上的對齊。

圖7 胸腔CT圖像配準(zhǔn)結(jié)果

圖8 目標(biāo)區(qū)域的重疊圖

在本次配準(zhǔn)實驗中，將本算法與互信息配準(zhǔn)算法［11］、Demons配準(zhǔn)算法［19］進(jìn)行比較，從互信息值、相對重疊率指標(biāo)、程序運行時間等方面評價算法精度與效率。配準(zhǔn)結(jié)果評價見表1，表中 MIp為帶有約束項的互信息相似性函數(shù)，MI配準(zhǔn)前后兩圖像的互信息值，ROI為感興趣區(qū)域的相對重疊率，Time為程序的運行時間。

從表1中可以看出，使用本算法對乳腺仿體的超聲數(shù)據(jù)與CT數(shù)據(jù)進(jìn)行配準(zhǔn)得到了較好的實驗結(jié)果：未改進(jìn)算法的感興趣區(qū)域的相對重疊率為89．03%，帶有薄板樣條約束項的B樣條自由形變非剛性配準(zhǔn)算法得到的感興趣區(qū)域的相對重疊率最高為91．28%，精度提高2．5%；Demons方法運行時間為1137s，本算法運行時間55s，運行時間提高了19．7倍。此外，用29組超聲圖像、3組CT圖像驗證本算法對臨床數(shù)據(jù)配準(zhǔn)的可行性，得到的實驗結(jié)果為：在非剛性配準(zhǔn)中具有薄板樣條約束項的配準(zhǔn)方法得到的感興趣區(qū)域重疊率為89．87%，無約束項的配準(zhǔn)方法其感興趣區(qū)域的相對重疊率為91．38%，精度提高1．7%左右；本算法的平均運行時間為266s，Demons方法運行時間為2276s，運行效率提高了7．6倍。

表1 配準(zhǔn)結(jié)果評價The evaluations of the image registration results

4 結(jié)論

本文將基于互信息的非剛性配準(zhǔn)方法應(yīng)用在超聲圖像和CT圖像配準(zhǔn)中，并在互信息中加入薄板樣條能量約束項，以避免配準(zhǔn)過程中圖像交叉與重疊問題的出現(xiàn)。聯(lián)合剛性變換和B樣條自由形變的非線性變換方法，使用L－BFGS的優(yōu)化方法對變換模型中的參數(shù)進(jìn)行多分辨率優(yōu)化搜索。在此基礎(chǔ)上，使用感興趣區(qū)域的相對重疊率、互信息值、配準(zhǔn)運行時間對本算法進(jìn)行評價，與沒有約束項的互信息方法以及微分同胚的Demons配準(zhǔn)方法比較。得到如下結(jié)論：未改進(jìn)算法得到的感興趣區(qū)域重疊率為89．58%，而帶有薄板樣條約束項的改進(jìn)算法得到的感興趣區(qū)域重疊率為91．35%，精度提高2．0%。Demons方法平均耗時1896s，本文改進(jìn)算法平均耗時195s，運行效率提高8．7倍。實驗結(jié)果表明此算法能夠?qū)Τ暫虲T圖像進(jìn)行配準(zhǔn)并獲得較好的結(jié)果。

［1］ Brown L G．A survey of image registration techniques［J］．ACM Computing Surveys（CSUR），1992，24（4）：325－376．

［2］ 田 捷，戴曉倩，楊 飛．醫(yī)學(xué)成像與醫(yī)學(xué)圖像處理教程［M］．北京：清華大學(xué)出版社，2013．152－229．Tian J，Dai X Q，Yang F．Medical Imaging and Medical Images Processing ［M］．Beijing：Tsinghua University Press，2013．187－189．

［3］ Wein W，Brunke S，Khamene A，et al．Automatic ct－ultrasound registration for diagnostic imaging and image－guided intervention［J］．Medical Image Analysis，2008，12（5）：577－585．

［4］ 盧振泰，陳武凡．基于主成分分析的三維醫(yī)學(xué)圖像快速配準(zhǔn)算法［J］．南方醫(yī)科大學(xué)學(xué)報，2008，28（9）：1591－1593．Lu Z T，Chen W F．A fast 3－d medical image registration algorithm using principal component analysis ［J］．Journal of Southern Medical University，2008，28（9）：1591－1593．

［5］ Rueckert D，Sonoda L I，Hayes C，et al．Nonrigid registration using free－form deformations：application to breast MR images ［J］．IEEE Transactions on Medical Imaging，1999，18（8）：712－721．

［6］ Roche A，Pennec X，Malandain G，et al．Rigid registration of 3－d ultrasound with MR images：a new approach combining intensity and gradient information［J］．IEEE Transactions on Medical Imaging，2001，20（10）：1038－1049．

［7］ Tanner C，Karssemeijer N，Székely G．Deformation models for registering MR and 3dultrasound breast images ［C］．IEEE International Symposium on Biomedical Imaging：From Nano to Macro，2011．582－585．

［8］ Holden M．A review of geometric transformations for nonrigid body registration［J］．IEEE Transactions on Medical Imaging，2008，27（1）：111－128．

［9］ Pluim J P W，Maintz J B A，Viergever M A．Mutual－information－based registration of medical images：a survey ［J］．IEEE Transactions on Medical Imaging，2003，22（8）：986－1004．

［10］ Zou X L，Navon M，Berger M，et al．Numerical experience with limited－memory quasi－newton and truncated newton methods［J］．SIAM Journal on Optimization，1993，3（3）：582－608．

［11］ Mattes D，Haynor D R，Vesselle H，et al．Pet－ct image registration in the chest using free－form deformations［J］．IEEE Transactions on Medical Imaging，2003，22（1）：120－128．

［12］ Rueckert D，Aljabar P，Heckemann R A，et al．Diffeomorphic registration using b－splines［C］．Medical Image Computing and Computer－assisted Intervention Miccai 2006．2006，Springer．p．702－709．

［13］ Xie Z Y，F(xiàn)arin G E．Image registration using hierarchical b－splines［J］．IEEE Transactions on Visualization and Computer Graphics，2004，10（1）：85－94．

［14］ Maes F，Collignon A，Vandermeulen D，et al．Multimodality image registration by maximization of mutual information ［J］．IEEE Transactions on Medical Imaging，1997，16（2）：187－198．

［15］ Nocedal J．Updating quasi－newton matrices with limited storage ［J］．Mathematics of Computation，1980，35（151）：773－782．

［16］ Likar B，Pernu F．A hierarchical approach to elastic registration based on mutual information ［J］．Image and Vision Computing，2001，19（1）：33－44．

［17］ Klein S，Staring M，Pluim J P W．Evaluation of optimization methods for non－rigid medical image registration using mutual information and b－splines［J］．IEEE Transactions on Image Processing，2007，16（12）：2879－2890．

［18］ Christensen G E，Geng X J，Kuhl J G，et al．Introduction to the non－rigid image registration evaluation project ［R］．Biomedical Image Registration，2006．Springer．128－135．

［19］ Vercauteren T，Pennec X，Perchant A，et al．Diffeomorphic demons：efficient non－parametric image registration［J］．NeuroImage，2009，45（1）：S61－S72．