張 冉，王 雷，夏 威，高 欣

（1.中國(guó)科學(xué)院長(zhǎng)春光學(xué)精密機(jī)械與物理研究所，吉林長(zhǎng)春130033；2.中國(guó)科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所醫(yī)學(xué)影像室，江蘇蘇州215163；3.中國(guó)科學(xué)院大學(xué)，北京100049）

2D／3D圖像配準(zhǔn)中的相似性測(cè)度和優(yōu)化算法

張 冉1，2，3，王 雷1，2，3，夏 威1，2，3，高 欣2

（1.中國(guó)科學(xué)院長(zhǎng)春光學(xué)精密機(jī)械與物理研究所，吉林長(zhǎng)春130033；2.中國(guó)科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所醫(yī)學(xué)影像室，江蘇蘇州215163；3.中國(guó)科學(xué)院大學(xué)，北京100049）

在手術(shù)引導(dǎo)治療中，2D／3D圖像配準(zhǔn)能輔助醫(yī)生準(zhǔn)確定位病人病灶，而準(zhǔn)確的配準(zhǔn)涉及相似性測(cè)度和優(yōu)化算法等眾多方面。為了研究相似性測(cè)度和優(yōu)化算法對(duì)2D／3D圖像剛性配準(zhǔn)的影響，本文結(jié)合6種相似性測(cè)度和4種優(yōu)化方法在配準(zhǔn)“金標(biāo)準(zhǔn)”數(shù)據(jù)上進(jìn)行了2D／3D圖像配準(zhǔn)實(shí)驗(yàn)，并從配準(zhǔn)成功率、平均迭代次數(shù)和平均配準(zhǔn)時(shí)間三個(gè)方面對(duì)配準(zhǔn)結(jié)果進(jìn)行了對(duì)比研究。實(shí)驗(yàn)結(jié)果表明，以模式強(qiáng)度為相似性測(cè)度，用Powell方法進(jìn)行優(yōu)化搜索是最佳配準(zhǔn)組合。并且，在不改變相似性測(cè)度條件下，Powell方法是所用優(yōu)化方法中配準(zhǔn)效果最好的優(yōu)化方法。

2D／3D圖像配準(zhǔn)；剛性配準(zhǔn)；相似性測(cè)度；優(yōu)化方法；配準(zhǔn)評(píng)估

1 引 言

在手術(shù)引導(dǎo)治療中，2D／3D圖像配準(zhǔn)在輔助醫(yī)生準(zhǔn)確定位病人病灶上是一種重要技術(shù)［1－3］，被廣泛的應(yīng)用在圖像引導(dǎo)的微創(chuàng)手術(shù)、放射治療的計(jì)劃制定、術(shù)后治療的效果檢驗(yàn)等介入手術(shù)方面。2D／3D圖像配準(zhǔn)的引入可降低上述醫(yī)療技術(shù)的實(shí)現(xiàn)難度，減少介入性治療的創(chuàng)傷并提高手術(shù)精度。為了實(shí)現(xiàn)準(zhǔn)確配準(zhǔn)，科研工作者提出了大量的配準(zhǔn)算法，基本可以分為基于特征［4－5］、基于灰度［4－7］和基于梯度［3］三類(lèi)方法。其中，基于灰度的配準(zhǔn)算法直接利用圖像的灰度值進(jìn)行配準(zhǔn)，無(wú)需進(jìn)行分割處理和人工介入，可實(shí)現(xiàn)全自動(dòng)配準(zhǔn)且精度較高，因而成為2D／3D醫(yī)學(xué)圖像剛性配準(zhǔn)的主要研究方向。該算法基于數(shù)字重建影像（Digitally Reconstructed Radiographs，DRR）技術(shù)，首先采用Ray casting算法［4］將變換后的3D圖像投影為2D圖像即DRR圖像，選取適當(dāng)?shù)南嗨菩詼y(cè)度作為目標(biāo)函數(shù)來(lái)衡量DRR圖像與X-ray圖像之間的相似程度，然后對(duì)目標(biāo)函數(shù)進(jìn)行優(yōu)化，當(dāng)函數(shù)達(dá)到最優(yōu)時(shí)即代表獲得最佳配準(zhǔn)。

近年來(lái)，在基于灰度的配準(zhǔn)算法中，科研人員對(duì)相似性測(cè)度和優(yōu)化方法做了一定的對(duì)比研究［4－5，8］，以評(píng)定相似性測(cè)度和優(yōu)化算法對(duì)配準(zhǔn)精度和魯棒性的影響。Penney等人［4］對(duì)歸一化相關(guān)（Normalized Cross Correlation，NCC）、互信息（Mutual Information，MI）、差值圖像的熵（Entropy of the Difference Image，EDI）、梯度相關(guān)（Gradient Correlation，GC）、梯度差分（Gradient Difference，GD）和模式強(qiáng)度（Pattern Intensity，PI）6種相似性測(cè)度進(jìn)行了對(duì)比分析；Maes等人［5］研究了Powell方法、單純形法（Simplex，SMP）、梯度下降法（Steepest Gradient Descent，STD）、共軛梯度法（Conjugate Gradient，CJG）、擬牛頓法（Quasi Newton，QSN）和列文伯格－馬夸爾特法（Levenberg Marquardt，LVM）6種優(yōu)化算法對(duì)基于互信息的2D／3D圖像配準(zhǔn)的影響，然而在2D／3D圖像配準(zhǔn)中，僅僅考慮相似性測(cè)度或者優(yōu)化方法對(duì)配準(zhǔn)結(jié)果的影響是不夠的，因此，本文結(jié)合6種相似性測(cè)度和4種優(yōu)化方法在公開(kāi)發(fā)表的“金標(biāo)準(zhǔn)”數(shù)據(jù)［2］上進(jìn)行2D／3D圖像剛性配準(zhǔn)實(shí)驗(yàn)，并利用Kraats等人［6］提出的方法從配準(zhǔn)成功率、平均迭代次數(shù)和平均配準(zhǔn)時(shí)間三個(gè)方面對(duì)配準(zhǔn)結(jié)果進(jìn)行評(píng)估，分析不同的相似性測(cè)度和優(yōu)化方法對(duì)圖像配準(zhǔn)的影響，并給出通用性最好的組合，為2D／3D圖像剛性配準(zhǔn)在臨床中的應(yīng)用提供理論基礎(chǔ)。

2 相似性測(cè)度

在2D／3D圖像剛性配準(zhǔn)中，相似性測(cè)度是用來(lái)度量?jī)煞鶊D像配準(zhǔn)程度的一個(gè)重要指標(biāo)，可以衡量經(jīng)過(guò)空間變換后浮動(dòng)圖像（DRR）與參考圖像（X-ray）在空間上的一致性程度。因此，相似性測(cè)度的選擇是否合理，直接影響到圖像配準(zhǔn)的準(zhǔn)確性和穩(wěn)定性。本文選取比較常用的6種相似性測(cè)度進(jìn)行對(duì)比研究，它們是歸一化相關(guān)、互信息、歸一化互信息、梯度相關(guān)、梯度差分以及模式強(qiáng)度［4，7－10］。為了便于統(tǒng)一標(biāo)識(shí)，定義IX( i，j)和IDRR( i，j)分別為X-ray圖像和DRR圖像中像素(i，j)處的灰度值。

2.1 歸一化互相關(guān)（NCC）

在2D／3D圖像剛性配準(zhǔn)中，歸一化互相關(guān)測(cè)度是通過(guò)求兩幅圖像的相關(guān)系數(shù)來(lái)表征二者間的相關(guān)程度，用公式可表示為：

其中，IX，IDRR分別表示X－ray圖像和DRR圖像像素灰度的平均值。兩幅圖像越相似，該測(cè)度的值越接近1。

2.2 互信息（MI）

互信息是信息論里一種有用的信息度量，是指兩個(gè)事件集合之間的相關(guān)性，一般用熵來(lái)表示。圖像X的熵和圖像X，Y的聯(lián)合熵分別定義為：

其中，聯(lián)合概率分布p( x，y)用歸一化的聯(lián)合直方圖表示，p( x)、p( y)是p( x，y)的邊緣分布概率。

因此，兩幅圖像之間的互信息可定義為：

其中，()H X，()H Y分別是圖像X，Y的熵；H X，()Y是它們的聯(lián)合熵。在2D／3D圖像剛性配準(zhǔn)中，互信息的值越大，兩幅圖像的相似性程度越高；當(dāng)互信息達(dá)到最大值時(shí)，表示兩幅圖像達(dá)到最佳配準(zhǔn)。

2.3 歸一化互信息（NMI）

由于互信息測(cè)度對(duì)配準(zhǔn)圖像間的重疊程度較為敏感，研究人員提出了歸一化互信息的概念，其中Sthdholme等人［10］給出了其中一種形式，可表示為：

2.4 梯度相關(guān)（GC）

基于梯度的相似性測(cè)度首先要計(jì)算出兩幅圖像的梯度圖像（本文采用Sobel算子），然后根據(jù)公式（1）分別計(jì)算出梯度圖像在水平和豎直兩個(gè)方向的歸一化互相關(guān)，它們的和就是兩幅圖像的梯度相關(guān)，可表示為：

2.5 梯度差分（GD）

這種測(cè)度也是基于梯度圖像，在得到各向梯度圖像的基礎(chǔ)上，比較兩幅梯度圖像間的差異，實(shí)際上是比較兩幅圖像的邊緣方向信息。其計(jì)算公式如下：

其中，Aν和Ah是X-ray圖像分別在豎直和水平方向上的方差；s是衰減因子。IdiffV和IdiffH分別表示待配準(zhǔn)圖像在豎直方向和水平方向上的梯度差值圖像。梯度差分測(cè)度通過(guò)使用1／（1＋x2）的形式增強(qiáng)了對(duì)細(xì)線結(jié)構(gòu)的魯棒性［9］，通過(guò)調(diào)整s，Aν和Ah的值可以使該測(cè)度更加平滑。

2.6 模式強(qiáng)度（PI）

該測(cè)度通過(guò)測(cè)量差值圖像（兩幅圖像的灰度差）中所存在的模式是否降到最低來(lái)判斷配準(zhǔn)是否成功，其數(shù)學(xué)表達(dá)式為：

其中，σ，r是常量；Idiff( v，w)表示差值圖像中(v，w)處的灰度值。模式強(qiáng)度認(rèn)為，當(dāng)一個(gè)像素與其臨近的像素的值的差異顯著時(shí)，該像素便屬于一個(gè)模式，圖像配準(zhǔn)的過(guò)程就是要盡量消除這種差異［7，9］。圖像達(dá)到最佳配準(zhǔn)時(shí)，Idiff中待配準(zhǔn)的模式會(huì)消失，其模式強(qiáng)度會(huì)最小，該測(cè)度的值最大。

3 優(yōu)化方法

除了相似性測(cè)度以外，優(yōu)化算法也會(huì)影響2D／ 3D圖像配準(zhǔn)的精度。醫(yī)學(xué)圖像剛性配準(zhǔn)在本質(zhì)上是一個(gè)多參數(shù)的最優(yōu)化問(wèn)題，即尋找相似性測(cè)度達(dá)到最大時(shí)所對(duì)應(yīng)的6個(gè)空間變換參數(shù)的值，可由下式表示：

式中，C是由相似性測(cè)度構(gòu)成的目標(biāo)函數(shù)；IX和I3D分別是2D X-ray圖像和3D體數(shù)據(jù)；μ是作用于體數(shù)據(jù)上的變換參數(shù)，μ′是圖像配準(zhǔn)時(shí)算法搜索到的最優(yōu)參數(shù)。

最優(yōu)化是一個(gè)迭代尋優(yōu)的過(guò)程，可以采用如下迭代格式［11］表示：

其中，dk表示第k次迭代時(shí)的搜索方向；αk是迭代步長(zhǎng)。不同的優(yōu)化算法對(duì)應(yīng)著不同的搜索方向和搜索步長(zhǎng)。文中選取4種優(yōu)化方法用于搜索圖像配準(zhǔn)的最優(yōu)空間變換參數(shù)，優(yōu)化方法具體如下［5，8，11－13］。

3.1 梯度下降法（GRadient Descent，GRD）

梯度下降法又稱(chēng)最速下降法，其搜索方向沿著負(fù)梯度方向，目標(biāo)函數(shù)可以最快的達(dá)到極小值，迭代格式可表示為：

3.2 非線性共軛梯度（Nonlinear Conjugate Gradient，NCG）

對(duì)于非線性共軛梯度，公式（13）中搜索方向dk可以表示為：

其中，βk公式有多種，Y H DAI，Y YUAN［12］給出βDYk的公式，Hestenes和Stiefel［8，12］給出βHSk的公式。本文采用的βk公式是二者的混合，即βk＝分別由下面兩個(gè)公式求得：

3.3 Powell算法

Powell算法［11］是直接搜索法中一種比較有效的方法，它不需要目標(biāo)函數(shù)的導(dǎo)數(shù)信息，只要求目標(biāo)函數(shù)連續(xù)即可。由于原始Powell算法中無(wú)法保證搜索方向是線性無(wú)關(guān)的，因此本文中采用的是改進(jìn)后的Powell算法［11］。

3.4 擬牛頓法（Quasi Newton，QN）

擬牛頓法是牛頓法的直接推廣，基本思想是構(gòu)造海森矩陣的近似矩陣Bk或逆海森矩陣的近似矩陣Hk，通過(guò)在試探點(diǎn)附近的二次逼近引導(dǎo)擬牛頓方程來(lái)確定搜索方向。擬牛頓方程主要有兩種形式：Hk＋1yk＝sk和Bk＋1sk＝y(tǒng)k，本文采用后者，其公式為：

其中，sk＝xk＋1－xk，yk＝gk＋1－gk，gk和gk＋1分別表示x在第k和第（k＋1）次的導(dǎo)數(shù)。

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

本文實(shí)驗(yàn)采用MATLAB軟件平臺(tái)編程實(shí)現(xiàn)配準(zhǔn)過(guò)程，同時(shí)結(jié)合CUDA技術(shù)，電腦配置為Intel Xeon E5645，2.4 GHz處理器，內(nèi)存為48 GB，NVIDIA Tesla C2050顯卡，顯存為4 GB。

4.1 實(shí)驗(yàn)數(shù)據(jù)

在配準(zhǔn)實(shí)驗(yàn)中，我們采用公開(kāi)發(fā)表的2D／3D配準(zhǔn)“金標(biāo)準(zhǔn)”數(shù)據(jù)［2］，該數(shù)據(jù)是對(duì)豬的顱骨進(jìn)行成像得到的，包含3D的CT體數(shù)據(jù)和2D的X-ray圖像。高斯濾波及各向同性采樣到1mm后，CT數(shù)據(jù)的大小為326×326×330，X-ray圖像尺度為410×410，如圖1（a）所示。由于配準(zhǔn)數(shù)據(jù)同時(shí)包含大量的軟組織和硬組織，為提高配準(zhǔn)精度，避免軟組織對(duì)剛性配準(zhǔn)的影響，同時(shí)縮短配準(zhǔn)時(shí)間，實(shí)驗(yàn)選取包含軟組織較少的圓形區(qū)域作為目標(biāo)區(qū)域，其直徑為200個(gè)像素。圖1（b）是去除標(biāo)記點(diǎn)的X-ray圖像的興趣區(qū)。

圖1 豬顱骨X－ray圖像

4.2 配準(zhǔn)評(píng)估

根據(jù)Kraats等人提出的2D／3D配準(zhǔn)評(píng)估準(zhǔn)則［6］，本文從配準(zhǔn)成功率（Success Rate，SR）、平均迭代次數(shù)以及平均配準(zhǔn)時(shí)間三個(gè)方面對(duì)配準(zhǔn)結(jié)果進(jìn)行評(píng)估。配準(zhǔn)成功率可以用平均目標(biāo)配準(zhǔn)誤差（mean target registration error，mTRE）來(lái)計(jì)算，mTRE公式如下：

其中，T和Tgold分別表示配準(zhǔn)后計(jì)算得到的變換矩陣T（由配準(zhǔn)算法決定）和金標(biāo)準(zhǔn)變換矩陣Tgold，pn為在CT體數(shù)據(jù)中隨機(jī)選取的目標(biāo)點(diǎn)集。對(duì)“金標(biāo)準(zhǔn)”變換參數(shù)進(jìn)行隨機(jī)擾動(dòng)后得到初始變換參數(shù)，使得初始mTRE分布在0～10 mm之間，并在每1 mm的間隔內(nèi)再隨機(jī)選取10個(gè)初始位置，這樣就產(chǎn)生了總共100次初始變換。

4.3 參數(shù)設(shè)置

本文中所有參數(shù)均是參照以往的文獻(xiàn)設(shè)置，對(duì)于梯度差分（GD），公式（7）中參數(shù)Aν和Ah均設(shè)為10，公式（8）和公式（9）中參數(shù)s均設(shè)為0.2。模式強(qiáng)度（PI）中參數(shù)σ設(shè)為10，r設(shè)為3，s設(shè)為0.2。

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

文中對(duì)6種相似性測(cè)度和4種優(yōu)化算法的不同組合進(jìn)行配準(zhǔn)對(duì)比，結(jié)果如表1和表2所示。由于“金標(biāo)準(zhǔn)”數(shù)據(jù)含有大量的軟組織并且在使用Ray Casting時(shí)沒(méi)進(jìn)行優(yōu)化處理來(lái)生成DRR圖像，因此本文選取的配準(zhǔn)成功閾值為8 mm，即最終目標(biāo)配準(zhǔn)誤差（mTRE）小于8 mm視為配準(zhǔn)成功。

表1 所有相似性測(cè)度和優(yōu)化算法所進(jìn)行配準(zhǔn)組合的成功配準(zhǔn)率

由表1可以看出，配準(zhǔn)成功率范圍為0.61～0.81，最高成功率（0.81）的配準(zhǔn)組合是PI和Powell組合，而最低成功率（0.61）的配準(zhǔn)組合是MI和NCG組合。從結(jié)果中還可以看出，相同的相似性測(cè)度，Powell是成功率最高的優(yōu)化方法，而NCG的成功率最低。

除了配準(zhǔn)成功率，本文還從平均配準(zhǔn)時(shí)間和平均迭代次數(shù)上對(duì)2D／3D剛性配準(zhǔn)進(jìn)行了評(píng)估，其中，平均配準(zhǔn)時(shí)間是指一次完整配準(zhǔn)所需的平均時(shí)間，平均迭代次數(shù)是指一次完整配準(zhǔn)所需的平均迭代次數(shù)。

表2 所有配準(zhǔn)組合的平均配準(zhǔn)時(shí)間（t）和平均迭代次數(shù)（N）

從表2可以看出，后三種相似性測(cè)度（GC、GD以及PI）無(wú)論用什么優(yōu)化策略進(jìn)行配準(zhǔn)所需時(shí)間都大于前面三種（NCC、MI以及NMI）。在所有的優(yōu)化方法中，Powell所需的平均配準(zhǔn)時(shí)間和平均迭代次數(shù)都是最少的。因此，在2D／3D圖像剛性配準(zhǔn)中，Powell成了最好的優(yōu)化方法。

5 結(jié) 論

本文對(duì)6種相似性測(cè)度和4種優(yōu)化方法的不同組合進(jìn)行了2D／3D剛性配準(zhǔn)對(duì)比研究，并從配準(zhǔn)成功率、平均配準(zhǔn)時(shí)間和平均迭代次數(shù)上進(jìn)行了配準(zhǔn)評(píng)估。結(jié)果表明，除了相似性測(cè)度會(huì)對(duì)2D／3D剛性配準(zhǔn)產(chǎn)生影響外，優(yōu)化方法的選擇同樣會(huì)影響配準(zhǔn)效果。就本文而言，以PI（模式強(qiáng)度）作為相似性測(cè)度，用Powell進(jìn)行優(yōu)化搜索是最好的配準(zhǔn)組合。此外，結(jié)果還表明Powell方法是2D／3D圖像剛性配準(zhǔn)最好的優(yōu)化方法。

［1］ PMarkelj，D Toma?evi?，B Likar，etal.A review of3D-2D registration methods for image-guided interventions［J］. Medical Image Analysis，2012，16（3）：642－661.

［2］ PMarkelj，F(xiàn) Pernu?，SA Pawiro，et al.Validation for 2D／3D registration I：a new gold standard data set［J］.Medical physics，2011，38（3）：1481－1490.

［3］ Barbara R?per，Nassir Navab，Wolfgang Wein.2D／3D registration based on volume gradients［C］／／Medical Imaging International Society for Optics and Photonics，2005：144－150.

［4］ Jürgen Weese，Graeme P Penney，John A Little，et al.A comparison of similaritymeasures for use in 2D3D medical image registration［J］.IEEE Transactions on Medical Imaging，1998，17（4）：586－595.

［5］ FMaes，Dirk Vandermeulen，Paul Suetens.Comparative evaluation ofmultiresolution optimization strategies formultimodality image registration bymaximization ofmutual information［J］.Medical Image Analysis，1999，3（4）：373－386.

［6］ Graeme P Penney，Everine B van de Kraats，Dejan Toma?evi?.Standardized evaluationmethodology for 2D－3D registration［J］.IEEE Transactions on Medical Imaging，2005，24（9）：1177－1189.

［7］ Liang Wei.2D－3D registration ofmedical image［D］. Nanjing：Southeast University，2004.（in Chinese）梁瑋.2D－3D醫(yī)學(xué)圖像配準(zhǔn)研究［D］.南京：東南大學(xué)，2004.

［8］ IM Jvan der Bom，SKlein，M Staring.Evaluation of optimizationmethods for intensity-based 2D－3D registration in x-ray guided interventions［C］／／SPIEMedical Imaging International Society for Optics and Photonics，2001，7962（23）：1－15.

［9］ Wang Kai.The study ofmultimodality medical image regitration［D］.Changsha：Central South University，2008.（in Chinese）王凱.多模態(tài)醫(yī)學(xué)圖像配準(zhǔn)研究［D］.長(zhǎng)沙：中南大學(xué)，2008.

［10］D L G Hill，C Studholme，D JHawkes.An overlap invariant entropymeasure of 3D medical image alignment［J］. Pattern Recognition，1999.

［11］Sun Jingming，Ling Yingchun.Optimal design ofmachine［M］.Beijing：China Machine Press，2006.（in Chinese）孫靖明，梁迎春.機(jī)械優(yōu)化設(shè)計(jì)［M］.北京：機(jī)械工業(yè)出版社，2006.

［12］Y H DAI，Y YUAN.An efficienthybrid conjugate gradient method for unconstrained optimization［J］.Annals of Operations Research，2001，103（1）：33－47.

［13］Ma Changfeng.Optimizationmethods and the design of its matlab program［M］.Beijing：Science Press，2010.（in Chinese）馬昌鳳.最優(yōu)化方法及其Matlab程序設(shè)計(jì)［M］.北京：科學(xué)出版社，2010.

Comparison of sim ilarity measurement and optim ization methods in 2D／3D image registration

ZHANG Ran1，2，3，WANG Lei1，2，3，XIAWei1，2，3，GAO Xin2
（1.Changchun Institute of Optics，F(xiàn)ine Mechanics and Physics，Chinese Academy of Sciences，Changchun 130033，China；2.Medical Imaging Department，Suzhou Institute of Biomedical Engineering and Technology，Chinese Academy of Sciences，Suzhou 215163，China；3.University of the Chinese Academy of Sciences，Beijing 100049，China）

In surgical guide treatment，2D／3Dmedical image registration can provide the precise position of patient for surgeon.Accurate registration involvesmany aspects，such as similaritymeasurements and optimizationmethods.In order to investigate the influence of similarity measurements and optimization methods on 2D／3D image registration，a comparison of six similaritymeasurements in combination with four optimizationmethods is performed using the public and available porcine skull phantom datasets from Medical University Vienna.Comparison is performed for the registration results based on success rate，the number of iterations and execution time.The results show that themost accuracy registration is obtained by pattern intensity combined with Powell.Furthermore，the best 2D／3D registration results are obtained by Powell search strategy with fixed similaritymeasurement.

2D／3D image registration；rigid registration；similaritymeasures；optimization；registration evaluation

TP391

A

10.3969／j.issn.1001-5078.2014.01.022

1001-5078（2014）01-0098-05

國(guó)家自然科學(xué)基金項(xiàng)目（No.81000651）；蘇州市科技計(jì)劃項(xiàng)目（No.SH201210）；江蘇省科技計(jì)劃項(xiàng)目（No.BL2012049）資助。

張 冉（1988－），女，碩士研究生，研究方向?yàn)獒t(yī)學(xué)圖像處理。

2013-05-20