有限元彈性配準中的驅(qū)動外力及其網(wǎng)格細化

張靜亞，李周雁

（1.常熟理工學(xué)院 物理與電子工程學(xué)院，江蘇 常熟 215500；2.中國工商銀行吳中支行，江蘇 蘇州 215006）

1 引言

醫(yī)學(xué)影像技術(shù)在現(xiàn)代醫(yī)療技術(shù)中發(fā)揮著不可替代的作用.對于醫(yī)學(xué)影像的分析是將醫(yī)學(xué)影像應(yīng)用于臨床的關(guān)鍵一步.圖像配準是醫(yī)學(xué)影像分析技術(shù)之一，其主要任務(wù)是獲取同一主體由于病情發(fā)展、攝像條件改變或者不同主體差異而產(chǎn)生的形變位移場［1-2］，廣泛應(yīng)用于疾病診斷、手術(shù)引導(dǎo)、治療效果跟蹤、放射計劃制定等環(huán)節(jié)中.

目前，圖像配準模型大致分為基于圖像特征和基于圖像灰度兩大類［3-4］.圖像特征通常包括圖像的邊緣輪廓和圖像特征點.配準時，首先提取圖像特征，然后通過特征匹配來獲取目標在空間上的對齊.該方法具有一定的局限性：首先，由于特征的不確定性使得特征選取成為一個非常耗時又困難的過程；其次，特征數(shù)量越多，算法的計算復(fù)雜度就越高.因此，對于不同配準目標選擇合理的特征種類和特征數(shù)量是一件困難的事情.基于灰度信息的配準技術(shù)是通過優(yōu)化由圖像灰度值定義的相似度函數(shù)來獲取圖像間的空間變換關(guān)系.此類技術(shù)無須進行特征選取，可靈活選擇相似度函數(shù)，并能結(jié)合金字塔多分辨率優(yōu)化方式，體現(xiàn)了更大的優(yōu)越性.基于連續(xù)介質(zhì)力學(xué)理論的彈性體形變模型［5-6］就是一種基于圖像灰度的配準模型，該模型把輸入圖像間的差異看作一個由連續(xù)介質(zhì)的物理運動而產(chǎn)生的形變位移場.與計算機圖形學(xué)所關(guān)注的圖形對象相比，醫(yī)學(xué)圖像一方面具有更為復(fù)雜的細節(jié)，另一方面對圖像內(nèi)容又具有物理學(xué)先驗知識，因此在醫(yī)學(xué)圖像的非剛性配準中，彈性體形變模型具有現(xiàn)實的物理意義.

本文首先采用彈性體形變理論建立能量函數(shù)，并采用有限元法求解當能量函數(shù)達到極值時的位移場.在能量函數(shù)中，驅(qū)動外力是由相似度梯度流推導(dǎo)得到的，因此相似度表達式對配準精度具有至關(guān)重要的影響.常用的相似度準則包括互信息（Mutual Information，MI）和灰度差平方和（Sum of squared differences，SSD）［7］.前者與圖像的灰度統(tǒng)計特性有關(guān)，具有對圖像之間亮度差異不敏感的優(yōu)點.但是計算過程中需要估計聯(lián)合概率密度，計算復(fù)雜度高，并且窗口尺寸對計算結(jié)果影響很大.SSD準則與圖像灰度差值有關(guān)，具有運算效率高的特點，但存在以下3方面缺點：第一，當圖像存在局部病灶異常、亮度不均勻等問題時，驅(qū)動點處的灰度差不能代表解剖結(jié)構(gòu)的差異，因而導(dǎo)致誤配準；第二，該相似度函數(shù)沒有結(jié)合區(qū)域信息，可能出現(xiàn)當全局相似度最優(yōu)時，圖像局部溢出的問題；第三，平方函數(shù)處理方式增加了溢出點處的估計誤差.針對以上問題，本文通過結(jié)合驅(qū)動點鄰域的空間信息和灰度信息，并引入結(jié)構(gòu)張量［8］和魯棒函數(shù)［9］，提出了一個新的驅(qū)動外力表達式，從而減少局部誤配準.

另外，網(wǎng)格剖分是有限元求解必不可少的步驟之一.均勻的三角化剖分是最常見的剖分方法之一，該方法不考慮圖像中由于目標解剖結(jié)構(gòu)的不同而存在的灰度差異，影響了配準的效果.針對該不足，本文采用區(qū)域細化的網(wǎng)格剖分法，該算法以圖像區(qū)域方差作為衡量圖像細節(jié)的標準，在圖像細節(jié)區(qū)域進行網(wǎng)格細化，從而減少了全局最優(yōu)與局部溢出之間的矛盾.本文提出的改進算法分別用于肺部CT圖像組和腦部MR圖像組的彈性配準實驗，實驗結(jié)果表明改進后的模型相對于傳統(tǒng)的有限元彈性模型具有更好的配準效果.

2 有限元驅(qū)動外力的改進

根據(jù)連續(xù)介質(zhì)彈性力學(xué)得到的彈性勢能和由驅(qū)動外力得到的動能建立起如式（1）所示的能量函數(shù)W（u），當能量函數(shù)達到極值時可獲得圖像配準所需的位移場

驅(qū)動外力F是由圖像相似度梯度流導(dǎo)出的，因此驅(qū)動外力的改進可通過對傳統(tǒng)相似度函數(shù)的改進來實現(xiàn).傳統(tǒng)的SSD相似度表達式

其中，r和m分別為參考圖像和浮動圖像灰度值.如前所述，（3）式所示的SSD準則具有運算效率高的優(yōu)點，但是缺乏對驅(qū)動點鄰域信息和圖像結(jié)構(gòu)信息的應(yīng)用，并且平方函數(shù)的處理方式加大了溢出點的估計誤差，導(dǎo)致配準結(jié)果對圖像局部亮度變化敏感，易產(chǎn)生局部溢出.本文針對以上問題，從3方面對（3）式進行了改進.

2.1 結(jié)構(gòu)張量項

對于結(jié)構(gòu)性圖像來說，結(jié)構(gòu)張量的跡能表達圖像中目標的解剖結(jié)構(gòu)，由文獻［8］可知，具有較小跡值的區(qū)域為圖像的平坦區(qū)域，較大跡值的區(qū)域為邊緣區(qū)域或角點區(qū)域.顯然，跡的大小與圖像中各點的灰度梯度密切相關(guān)，而與圖像的亮度無關(guān).根據(jù)以上分析，本文通過在相似度函數(shù)中加入結(jié)構(gòu)張量的跡的差值平方和來解決圖像亮度不均勻問題.x點處的結(jié)構(gòu)張量通常可表示為.結(jié)構(gòu)張量的跡表示為.其中，是x點處的灰度值，是空間梯度算子.將結(jié)構(gòu)張量的跡加入到（3）式中，得到相似度函數(shù),.其中，α為權(quán)重系數(shù)，實驗統(tǒng)計顯示，在處理亮度不均勻圖像時，α取值可稍大α=0，提高結(jié)構(gòu)張量項的比重，其余可取1.

2.2 驅(qū)動點鄰域信息

為了減少局部溢出，本文將式（3）所示的相似度計算從驅(qū)動點的灰度差值拓展到驅(qū)動點鄰域的灰度差值.另外，本文利用一個均值為驅(qū)動點坐標x，標準差為γ的高斯函數(shù)Gr對鄰域內(nèi)各點的灰度差值平方進行加權(quán)，實現(xiàn)由鄰域內(nèi)各點與驅(qū)動點相對位置來調(diào)節(jié)相似度的目的，用Sl（u）表示：

通過實驗統(tǒng)計分析，鄰域半徑過大并不能帶來更好的圖像配準結(jié)果，取值在1～2之間較為合理，對于噪聲小的CT圖像可取1，處理噪聲相對大的MR圖像時可取2.γ與鄰域各點的灰度權(quán)重有關(guān)，本文取γ=3.

2.3 魯棒函數(shù)

由于式（3）中的平方函數(shù)處理方式增加了溢出點處的估計誤差，受文獻［9］的啟發(fā)，本文采用了既能減弱溢出點的不利影響又能保證相似度函數(shù)可微性和凸函數(shù)性的魯棒函數(shù):

3 模型的有限元法求解和區(qū)域網(wǎng)格細化

3.1 彈性力學(xué)模型的有限元求解

將式（1）所示的能量函數(shù)W（u）對u求導(dǎo)并令其為0，可求得能量函數(shù)達到極小值時的u值.

對（6）式所表達的非線性方程可采用有限元法求解.圖像各點的位移值由有限元結(jié)點位移與形狀函數(shù)插值而得到，其中Nd為元素結(jié)點數(shù)，插值表達式為

將式（7）和式（1）代入到式（6），經(jīng)整理可得到如式（8）所示的矩陣方程，該方程的解即為結(jié)點位移向量.

剛度矩陣K：剛度矩陣K由各元素剛度矩陣Kel組合得到，元素剛度矩陣Kel如下：

3.2 區(qū)域網(wǎng)格細化

由式（5）可以看出，驅(qū)動外力F（u，x）與驅(qū)動點處的灰度梯度密切相關(guān).另外，由（7）式看出，圖像中任意點的位移是由結(jié)點位移經(jīng)過有限元插值得到的.因此，在圖像邊緣和細節(jié)較多區(qū)域的有限元結(jié)點位移值對提高配準準確度起到了非常重要的作用.本文通過統(tǒng)計圖像灰度的區(qū)域方差并設(shè)定閾值來標注圖像的邊緣和局部細節(jié)較多的區(qū)域，并在這些區(qū)域插入細化結(jié)點，進行網(wǎng)格細化.

如圖1（a）所示，{Ω1，Ω2，Ω3，Ω4}為4個標志值為1的細化區(qū)域.通常，對區(qū)域進行細化剖分需要兩個步驟：首先，在各細化區(qū)域插入如圖1（a）所示的細化點（由實點表示），然后對所有結(jié)點進行Delaunay三角形網(wǎng)格生成，得到如圖1（c）所示的細化網(wǎng)格.

圖2 給出了大腦MR圖像的剖分結(jié)果，（a ）圖是大腦MR圖像，（b）圖是細化標志矩陣A的圖像顯示，（c）圖是均勻三角化剖分結(jié)果，（d）圖是區(qū)域細化剖分結(jié)果.

圖2 腦部MRI圖像的網(wǎng)格剖分圖示

4 結(jié)果和分析

為驗證本文提出算法的有效性，以大量醫(yī)學(xué)圖像為實驗對象進行了測試與分析.為了定量評估配準算法的性能，本文選用了均方誤差MSE、峰值信噪比PSNR、歸一化互信息NMI、相關(guān)系數(shù)NC和尺度不變特征變換（SIFT）特征點平均距離FDist這5種度量標準.假設(shè)經(jīng)過SIFT特征點匹配及去誤［11-12］后，所得特征點分別為P=（p1，p2，…，pL）和Q=（q1，q2，…，qL），用x（pi）表示特征點pi處的SIFT特征點坐標，則匹配特征點歐氏距離平均值為

均方誤差MSE、峰值信噪比PSNR、歸一化互信息NMI、歸一化相關(guān)系數(shù)NC分別定義為

本文首先對經(jīng)過樣條函數(shù)形變的20組肺部CT圖進行了配準實驗.圖3給出了其中一組實驗結(jié)果，圖3（a）為參考圖像.圖3（b）是形變的肺部CT圖，作為浮動圖像.圖3（c）～3（e）分別給出了傳統(tǒng)彈性力學(xué)模型、驅(qū)動外力改進算法及加入?yún)^(qū)域細化法的配準結(jié)果.圖3（c）的結(jié)果顯示，利用傳統(tǒng)彈性力學(xué)模型可對部分形變進行校準，但是局部邊緣處存在誤差，可從圖3（c）中方形框區(qū)域看出.利用驅(qū)動外力的改進算法和進一步的區(qū)域細化剖分算法后，可取得更好的配準結(jié)果，如圖3（d）、圖3（e）所示.從圖3（d）的框區(qū)域看出，沒經(jīng)過區(qū)域細化的配準結(jié)果在局部細節(jié)處仍存在小形變，而加入細化剖分后在細節(jié)處具有更好的對齊效果，如圖3（e）所示.從3者與參考圖像的差值圖上也可以得到這些結(jié)論.其中黑色代表相同，白色代表差別.表1給出了以上圖像配準結(jié)果的定量評價.

圖3 CT圖的配準結(jié)果

表1 CT圖配準結(jié)果的相似度度量值比較

對20組CT圖像配準結(jié)果進行定量評價，得到了MSE，PSNR，NMI，NC和FDist的平均值，如表2所示.從表中可看出，本文提出的改進驅(qū)動外力配合局部細化剖分算法相對傳統(tǒng)算法取得了顯著的提高，t檢驗得到的P值遠小于0.05.

MR圖表現(xiàn)出清晰的軟組織細節(jié)，對醫(yī)學(xué)診斷具有很好的指導(dǎo)意義.因此，圖像配準必須能較好地保留這些細節(jié)信息.由于參考圖像和浮動圖像之間存在明顯的亮度差異，由圖4（c）可看出，基于傳統(tǒng)的彈性力學(xué)模型的配準結(jié)果中存在明顯的內(nèi)部細節(jié)沒有對齊，并且產(chǎn)生了局部溢出.圖4（d）是采用本文改進的驅(qū)動外力和區(qū)域細化算法后得到的配準結(jié)果.由圖可看出，圖像匹配的效果得到了明顯改善，圖像內(nèi)部細節(jié)都得到了較好的配準.

表2 CT圖配準結(jié)果的相似度度量平均值

圖4 腦部MR圖的配準

表3給出了對以上配準結(jié)果進行定量評價的度量值.為了通過實驗結(jié)果的統(tǒng)計數(shù)值來更客觀地評價算法的性能，我們對20對不同人腦部MR圖像進行了配準實驗，對實驗結(jié)果統(tǒng)計分析得到的度量平均值如表4所示.由表4可以看出本文提出的算法在處理MR圖像配準上有了顯著提高，t檢驗P值小于0.05.實驗結(jié)果充分說明本文對驅(qū)動外力的改進以及局部細化區(qū)域的改進算法對醫(yī)學(xué)圖像配準起到了積極的作用.

表3 腦部MR圖配準結(jié)果相似度度量值比較

表4 20對腦部MR圖配準結(jié)果的相似度平均值比較

5 結(jié)論

本文基于有限元彈性配準模型估計了彈性物體在外力作用下的形變場.在以傳統(tǒng)SSD梯度流推導(dǎo)得到的驅(qū)動外力方法中，本文采用與鄰域空間相關(guān)的高斯函數(shù)對驅(qū)動點鄰域的灰度差平方加權(quán)，并應(yīng)用結(jié)構(gòu)張量和魯棒函數(shù)，從而得到了一個改進的驅(qū)動外力，減少了局部誤配準.在算法的有限元求解中，本文根據(jù)區(qū)域方差概率分布函數(shù)得到一個由相對數(shù)閾值確定的細化標志矩陣，進一步生成了區(qū)域細化網(wǎng)格，從而提高了結(jié)點位移的計算準確度.通過大量CT和MR圖像的配準實驗，分析比較了本文算法的有效性.實驗結(jié)果顯示，本文改進算法的配準準確度相對于傳統(tǒng)彈性配準算法有了顯著提高，t檢驗得到的P值小于0.05.