基于交叉累計(jì)剩余熵的圖像配準(zhǔn)中插值方法的改進(jìn)

相 艷，賀建峰，張?jiān)拼海?莉

（1.云南省計(jì)算機(jī)技術(shù)應(yīng)用重點(diǎn)實(shí)驗(yàn)室，昆明 650500；2.昆明理工大學(xué)信息工程與自動化學(xué)院，昆明 650500；3.云南大學(xué)軟件學(xué)院，昆明 650091）

基于交叉累計(jì)剩余熵的圖像配準(zhǔn)中插值方法的改進(jìn)

相 艷1，2，賀建峰2，張?jiān)拼?，蔡 莉3

（1.云南省計(jì)算機(jī)技術(shù)應(yīng)用重點(diǎn)實(shí)驗(yàn)室，昆明 650500；2.昆明理工大學(xué)信息工程與自動化學(xué)院，昆明 650500；3.云南大學(xué)軟件學(xué)院，昆明 650091）

交叉累計(jì)剩余熵（CCRE）比傳統(tǒng)互信息在配準(zhǔn)強(qiáng)噪聲圖像時(shí)更具優(yōu)勢，但采用部分體積（PV）插值的CCRE在網(wǎng)格點(diǎn)容易產(chǎn)生局部極值，不利于變換參數(shù)的優(yōu)化。針對該問題，研究基于3階B樣條函數(shù)的PV插值（BPV）、哈寧窗sinc函數(shù)的PV插值（HPV）和Blackman-Harris窗sinc函數(shù)的PV插值（BHPV）方法在CCRE中的應(yīng)用，提出一種新的插值方法。該方法采用靈活的鄰域中心，將插值點(diǎn)對聯(lián)合直方圖貢獻(xiàn)的權(quán)重分散到臨近的9個(gè)點(diǎn)上，并使用高斯函數(shù)作為PV插值的核函數(shù)，避免權(quán)重突變。實(shí)驗(yàn)結(jié)果表明，與BPV，HPV和BHPV插值方法相比，該方法對噪聲圖像的配準(zhǔn)率較高，配準(zhǔn)速度較快，更適合應(yīng)用于CCRE的計(jì)算。

圖像配準(zhǔn)；交叉累積剩余熵；部分體積插值；高斯函數(shù)；3階B樣條函數(shù)；sinc函數(shù)

DO I：10.3969/j.issn.1000-3428.2015.10.037

1 概述

基于互信息（Mutual Information，M I）的圖像配準(zhǔn)方法具有自動化程度高、配準(zhǔn)精度高等優(yōu)點(diǎn)。但當(dāng)圖像采集過程中受噪聲影響嚴(yán)重，或2幅待匹配圖像的重疊部分較少或分辨率較低時(shí)，互信息曲線會出現(xiàn)局部極值，變得極不光滑。 累積剩余熵（Cumulative Residual Entropy，CRE）是由文獻(xiàn)［1］提出的一種新的信息度量方式，該方法將香農(nóng)熵定義中的概率密度函數(shù)用概率分布函數(shù)來代替。Wang等學(xué)者將累計(jì)剩余熵用在圖像配準(zhǔn)中，提出將2幅圖像的交叉累計(jì)剩余熵（Cross Cumulative Residual Entropy，CCRE）作為配準(zhǔn)的相似性測度［2］。CCRE是統(tǒng)計(jì)大于某一灰度值的所有直方圖信息，保持了局部區(qū)域灰度信息的連續(xù)性，從而克服噪聲對局部極值的影響，能得到較 M I更光滑的配準(zhǔn)目標(biāo)函數(shù)［3-4］。

盡管如此，CCRE仍然會出現(xiàn)局部極值，其原因一方面與待配準(zhǔn)的2幅圖像存在一定的局部匹配有關(guān)，另一方面與浮動圖像經(jīng)空間變換后的插值方法有關(guān)。CCRE與M I一樣需要計(jì)算2幅圖像的聯(lián)合直方圖，因此，MI的插值方法可供借鑒。傳統(tǒng)的插值方法有最近鄰法、雙線性法等，但兩者精度都不高。目前被廣泛采用的是部分體積插值法（Partial Volume，PV）。該方法無需進(jìn)行浮動圖像的灰度插值，而是直接更新聯(lián)合直方圖。然而，PV插值的一個(gè)明顯缺陷是在網(wǎng)格點(diǎn)產(chǎn)生聯(lián)合熵的突變，從而引起局部極值［5-6］。 針對這一問題，文獻(xiàn)［6］利用均勻B樣條基函數(shù)拓展了PV插值法的鄰域，文獻(xiàn)［7-8］提出一種新的基于 Hanning窗的PV插值方法，文獻(xiàn)［9］用Blackman-Harris窗sinc函數(shù)取代傳統(tǒng) PV插值法的核函數(shù)。上述方法都是將某點(diǎn)對聯(lián)合直方圖的影響從4個(gè)鄰域點(diǎn)擴(kuò)展到16個(gè)，有效抑制了MI函數(shù)曲線中的局部極值現(xiàn)象。

本文將PV插值方法以及上述3種改進(jìn)的PV插值方法應(yīng)用到CCRE的計(jì)算中，驗(yàn)證3種改進(jìn)方法同樣能抑制PV插值的CCRE曲線中網(wǎng)格點(diǎn)的局部極值，但運(yùn)行時(shí)間較長。同時(shí)提出一種新的插值方法，使用高斯函數(shù)作為PV插值的核函數(shù)，將插值點(diǎn)對聯(lián)合直方圖的貢獻(xiàn)分散到其臨近的9個(gè)點(diǎn)上，以緩解網(wǎng)格點(diǎn)權(quán)重的突變，得到平滑的CCRE曲線，并且提高運(yùn)行速度。

2 交叉累計(jì)剩余熵的計(jì)算原理

CRE是一種新的對隨機(jī)變量的信息測度。設(shè)χ是R中的一個(gè)隨機(jī)變量，并且F（λ）是累計(jì)剩余分布函數(shù)。把香農(nóng)熵定義中的概率密度函數(shù)用其自身的累計(jì)剩余分布函數(shù)代替，則χ的CRE定義為：

其中，χ∈R+。

累積剩余熵比香農(nóng)熵更具有一般性，其定義在連續(xù)和離散區(qū)域都是合理的［10］。圖像 X和圖像 Y的CCRE則定義為：

具體用IR（χ）表示參考圖像，IT（χ）表示浮動圖像，變換矩陣為 g（χ；μ），μ為一系列的變換參數(shù)。配準(zhǔn)問題即轉(zhuǎn)化為尋找最佳變換參數(shù) μ，使得經(jīng)過變換矩陣變換后的浮動圖像 IT°g（χ；μ）與參考圖像IR（χ）的CCRE達(dá)到最大值，即：

定義P（l，k；μ）為（IT°g（χ；μ），IR）的聯(lián)合概率密度，PT（l；μ）和PR（k）分別代表浮動圖像和參考圖像的邊緣概率密度，LT和LR分別代表浮動和參考圖像的離散密度點(diǎn)集，則2幅圖像的CCRE可具體表示為：

3 不同插值方法及其在CCRE計(jì)算中的應(yīng)用

3.1 PV插值法及其改進(jìn)

計(jì)算式（4）中的聯(lián)合概率密度時(shí)，浮動圖像的像素經(jīng)過空間變換映射到參考圖像中時(shí)通常為非整數(shù)點(diǎn)，需要進(jìn)行插值處理。常用的PV插值法是根據(jù)權(quán)重分配原則，將每對像素對聯(lián)合直方圖的貢獻(xiàn)按權(quán)重分配給相鄰的4個(gè)像素對。然而，當(dāng)插值點(diǎn)正好落在網(wǎng)格點(diǎn)處時(shí)，則該網(wǎng)格點(diǎn)處權(quán)重為1，而其他相鄰的3個(gè)網(wǎng)格點(diǎn)權(quán)重為0。因此，非網(wǎng)格點(diǎn)和網(wǎng)格點(diǎn)的PV插值所更新的像素對比例為4：1，聯(lián)合直方圖的雜散程度的差異明顯，導(dǎo)致CCRE容易在整數(shù)倍位移處出現(xiàn)局部極值。文獻(xiàn)［6-9］的3種插值方法都是將插值點(diǎn)對聯(lián)合直方圖的貢獻(xiàn)從4個(gè)鄰域點(diǎn)擴(kuò)展到16個(gè)，采用不同的PV插值核函數(shù)，從而使得插

值點(diǎn)的權(quán)重分布更合理。PV插值改進(jìn)方法示意圖如圖1所示。

圖1 PV插值改進(jìn)方法示意圖

如圖1所示，設(shè)浮動圖像 IT（χ）上某點(diǎn) P經(jīng)空間變換后與參考圖像 IR（χ）上的點(diǎn) q對應(yīng)，點(diǎn) q的鄰域點(diǎn)為qi（i=1，2，…，16），其聯(lián)合直方圖的計(jì)算如下：

其中，ωi為每個(gè)領(lǐng)域點(diǎn)qi的權(quán)重；

在式（6）中，diχ和 diy為點(diǎn)qi與q在 χ和 y方向上的距離，而PV插值核函數(shù)f有以下3種形式：

（1）f為3階 B樣條函數(shù)［6，11］：

（2）f為Hanning窗sinc函數(shù)的近似［7-8］：

（3）f為Blackman-Harris窗sinc函數(shù)［9］：

這3種PV插值的改進(jìn)方法的共同點(diǎn)是：當(dāng)點(diǎn)q落在非網(wǎng)格點(diǎn)上時(shí)，參與聯(lián)合直方圖統(tǒng)計(jì)的鄰域點(diǎn)數(shù)目有16個(gè)；而當(dāng)點(diǎn)q落在網(wǎng)格點(diǎn)q6上時(shí)，參與聯(lián)合直方圖統(tǒng)計(jì)的鄰域點(diǎn)數(shù)目是最少的，有9個(gè)。能分配到權(quán)重的最多鄰域點(diǎn)數(shù)與最少鄰域點(diǎn)數(shù)的比例為16：9，比傳統(tǒng) PV插值法的比例4：1要小。這個(gè)比例越小，說明聯(lián)合直方圖的概率分布越分散［9］，CCRE出現(xiàn)局部極值的機(jī)率降低。

3.2 本文插值方法

上述插值方法存在一定的不足：（1）能分配到權(quán)重的最多鄰域點(diǎn)數(shù)與最少鄰域點(diǎn)數(shù)的比例沒有達(dá)到1:1。（2）插值點(diǎn)權(quán)重需要分散到鄰域16個(gè)點(diǎn)，應(yīng)用到CCRE配準(zhǔn)中速度較慢。（3）利用插值點(diǎn)與鄰域點(diǎn)的空間距離來決定權(quán)重大小，但是固定的鄰域范圍影響了權(quán)重分配的公平性。這類插值方法都以 q的整數(shù)點(diǎn) q6為中心形成16個(gè)點(diǎn)的鄰域。當(dāng)q較靠近q6時(shí)，這種鄰域范圍是比較合理的。但是當(dāng) q遠(yuǎn)離q6時(shí)，則這一鄰域范圍對于權(quán)重分配不太合理。例如當(dāng)q與q6在 χ和 y方向的距離都為0.9時(shí)，q與q1的歐氏距離為2.687 0，與q2和q5的歐氏距離為2.102 4，與q4和q13的歐式距離為2.195 4。而q與鄰域范圍外的q19和 q23的歐式距離為2.102 4，與q18和q22的歐氏距離為2.284 7，與q20和q24歐式距離為2.230 7。若按空間距離的遠(yuǎn)近來決定權(quán)重分配的大小，則這些鄰域范圍外的點(diǎn)也理應(yīng)被分配到權(quán)重，但在這3種插值方法中因?yàn)椴辉?6個(gè)鄰域點(diǎn)之中而被忽略。

根據(jù)以上分析，本文提出一種基于9個(gè)鄰域點(diǎn)的插值方法。

如圖2所示，其中，虛線網(wǎng)格表示整數(shù)點(diǎn)之間的平分線。設(shè)浮動圖像上某點(diǎn) P經(jīng)空間變換后與參考圖像IR（χ）上的點(diǎn) q對應(yīng)，q的整數(shù)點(diǎn)為 q5。選取 q5，q6，q8，q9中與點(diǎn) q距離最近的點(diǎn)作為鄰域中心，將權(quán)重分配給包括該中心在內(nèi)的 9個(gè)鄰域點(diǎn)。例如圖 2中 q離 q5較近，則權(quán)重分配給 qi（i=1，2，…，9）；而另一插值點(diǎn) q′的整數(shù)點(diǎn)同樣為q5，但其離 q9較近，則其權(quán)重分配給以 q9為中心的9個(gè)鄰域點(diǎn)，即

圖2 本文插值方法示意圖

確定了9個(gè)鄰域點(diǎn)后，各點(diǎn)的聯(lián)合直方圖也采

用諸如式（5）和式（6）進(jìn)行計(jì)算，只是式（5）中的16個(gè)點(diǎn)變?yōu)?個(gè)點(diǎn)。同時(shí)，本文提出一種新的核函數(shù)f來計(jì)算式（6）。

根據(jù)權(quán)重分配的原理，核函數(shù) f應(yīng)當(dāng)滿足［6］：（1）qi與 q距離越小，式（6）中的 ωi越大；（2）ωi非負(fù)經(jīng)過合理設(shè)計(jì)，形如式（10）的高斯函數(shù)也能滿足以上要求：

式（10）中標(biāo)準(zhǔn)差σ決定了權(quán)重ωi的分布，歸一化系數(shù)α的作用是將9個(gè)領(lǐng)域點(diǎn)的權(quán)重和調(diào)整到約等于1。經(jīng)分析，當(dāng)時(shí)，權(quán)重和最接近1。式（10）可簡化為：

式（11）的函數(shù)曲線如圖 3所示。可以看出，f（χ）在χ=1.5時(shí)并非為0。意味著與q相距最遠(yuǎn)的網(wǎng)格點(diǎn)（兩者χ和y方向的距離為1.5）也被分配到權(quán)重，即此時(shí)參與聯(lián)合直方圖統(tǒng)計(jì)的鄰域點(diǎn)數(shù)目仍有9個(gè)。因此，無論q落在何處，能分配到權(quán)重的最多鄰域點(diǎn)數(shù)與最少鄰域點(diǎn)數(shù)的比例為9：9，比BPV，HPV，BHPV的比例還小。

圖3 本文方法采用的高斯函數(shù)

將Vanderbilt圖像配準(zhǔn)數(shù)據(jù)集［12］中training-001的PET第11斷層圖像作為參考圖像，將其沿χ軸平移-1.5～1.5像素作為浮動圖像，同時(shí)為參考和浮動圖像添加密度為0.1的椒鹽噪聲。采用基于PV，BPV，HPV，BHPV以及本文的插值方法計(jì)算參考和浮動圖像的CCRE，平移步長為0.1像素，得到如圖4所示的結(jié)果。可以看出，PV插值的CCRE在網(wǎng)格點(diǎn)位置出現(xiàn)了擾動，而其他4種改進(jìn)方法的CCRE曲線都較光滑。隨著噪聲的增加，PV插值的CCRE在網(wǎng)格點(diǎn)位置將有更為明顯的波動，這些局部極值會造成優(yōu)化算法尋優(yōu)的困難，甚至陷入局部極值造成配準(zhǔn)失敗。

圖4 浮動圖像沿χ軸平移的交叉累計(jì)剩余熵曲線

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

為進(jìn)一步驗(yàn)證插值方法的性能，本文設(shè)計(jì)了2組實(shí)驗(yàn)。第1組實(shí)驗(yàn)的參考和浮動圖像為training-001的PET第11斷層圖像，并添加了密度為0.3的椒鹽噪聲。設(shè)剛體變換參數(shù)T=［tχ，ty，θ］，其中，tχ，ty分別為浮動圖像沿χ和y軸的平移，θ為繞圖像中心點(diǎn)的旋轉(zhuǎn)角度。兩者真實(shí)的配準(zhǔn)參數(shù)為T=［5，5，5］。使用步長為2的Pow ell算法對不同插值方法得到的CCRE進(jìn)行尋優(yōu)，結(jié)果如表1所示。

表1 不同插值方法對同幅PET圖像的配準(zhǔn)結(jié)果

從表1的配準(zhǔn)結(jié)果可知，PV插值法的CCRE使得Powell難以收斂，迭代次數(shù)達(dá)到16次，最終得到錯誤的配準(zhǔn)參數(shù)。BPV，HPV，BHPV和本文方法的CCRE都能得到正確的配準(zhǔn)結(jié)果，且本文方法的χ軸平移誤差最小。另外，本文方法的運(yùn)行時(shí)間較BPV縮短了87.83%，較HPV縮短了92.99%，較BHPV縮短了98.25%，配準(zhǔn)效率顯著提高。

第2組實(shí)驗(yàn)的參考圖像為BrainWeb數(shù)據(jù)集［13］中正常腦部磁共振 T1加權(quán)第90斷層圖像，浮動圖像為相應(yīng)斷層的T2加權(quán)像，并根據(jù)文獻(xiàn)［14］的方法添加標(biāo)準(zhǔn)差為22的萊斯噪聲。兩者真實(shí)的配準(zhǔn)參數(shù)為T=［6，-8，5］。使用步長為5的Powell算法對5種不同的插值方法得到的CCRE進(jìn)行尋優(yōu)，結(jié)果如表2所示。

表2 不同插值方法對T1和T2圖像的配準(zhǔn)結(jié)果

同樣，PV插值法得到錯誤的配準(zhǔn)參數(shù)，而BPV，BHPV出現(xiàn)了配準(zhǔn)失誤。HPV和本文方法的CCRE都能得到正確的配準(zhǔn)結(jié)果，且本文方法的精度最高，運(yùn)行時(shí)間較HPV縮短了95.12%。

5 結(jié)束語

本文針對使用PV插值計(jì)算CCRE時(shí)在網(wǎng)格點(diǎn)容易產(chǎn)生局部極值的問題，提出一種基于高斯函數(shù)的改進(jìn)PV插值方法。通過對高斯函數(shù)的方差及歸一化系數(shù)的合理設(shè)計(jì)，可以滿足權(quán)重分布對插值核函數(shù)的要求。該方法的優(yōu)勢為：（1）采用靈活的鄰域中心來決定鄰域范圍，保證了權(quán)重分配的公平性。（2）分配到權(quán)重的最多鄰域點(diǎn)數(shù)與最少鄰域點(diǎn)數(shù)的比例達(dá)到1:1。（3）插值點(diǎn)權(quán)重分散到鄰域9個(gè)點(diǎn)，相比16個(gè)點(diǎn)BPV，HPV，BHPV，大幅提高了配準(zhǔn)效率。實(shí)驗(yàn)結(jié)果表明，將本文插值方法應(yīng)用到CCRE中時(shí)，即使待配準(zhǔn)圖像受到強(qiáng)噪聲影響，也能優(yōu)化得到正確的參數(shù)，且運(yùn)行速度較快。下一步將研究該插值方法在其他相似性測度中的應(yīng)用。

［1］ Rao M，Chen Y，Vemuri B C，et al.Cumulative Residual Entropy：A New Measure of Information［J］.IEEE Transactions on Information Theory，2004，50（6）：1220-1228.

［2］ Wang Fei，Vemuri B C.Non-rigid Multi-modal Image Registration Using Cross-cumulative Residual Entropy［J］. International Journal of Computer Vision，2007，74（2）：201-215.

［3］ Hasan M，Pickering M R，Jia Xiuping.Robust Automatic Registration of Multimodal Satellite Images Using CCRE with Partial Volume Interpolation［J］.IEEE Transactions on Geoscience and Remote Sensing，2012，50（10）：4050-4061.

［4］ 江萬壽，彭芳媛，岳春宇，等.利用Ratio梯度和交叉累積剩余熵進(jìn)行多源遙感影像匹配［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2009，34（9）：1047-1050.

［5］ 馮 林，嚴(yán) 亮，賀明峰，等.圖像配準(zhǔn)中確定性擾動PV插值算法［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2005，17（5）：908-914.

［6］ 彭景林，章 兢，李樹濤.基于改進(jìn)PV插值和混合優(yōu)化算法的醫(yī)學(xué)圖像配準(zhǔn)［J］.電子學(xué)報(bào)，2006，34（5）：963-965.

［7］ Lu Xuesong，Zhang Su，Su He，et al.Mutual Information-based Multimodal Image Registration Using a Novel Joint Histogram Estimation［J］.Computerized Medical Imaging and Graphics，2008，32（3）：202-209.

［8］ Lu Xuesong，Zhang Su，Su He，et al.Non-rigid Medical Image Registration with Joint Histogram Estimation Based on Mutual Information［J］.Transactions of Tianjin University，2007，13（6）：452-455.

［9］ 陳偉卿，華順剛，歐宗瑛.互信息醫(yī)學(xué)圖像配準(zhǔn)中 PV插值算法的改進(jìn)［J］.計(jì)算機(jī)工程與應(yīng)用，2010，46（20）：113-115.

［10］ 李 超，陳 錢，錢惟賢.基于交叉累計(jì)剩余熵的多光譜圖像配準(zhǔn)方法［J］.紅外與激光工程，2013，42（7）：1866-1870.

［11］ Lehmann T M，Gonner C，Spitzer K.Survey：Interpolation Methods in Medical Image Processing［J］.IEEE Transactions on Medical Imaging，1999，18（11）：1049-1075.

［12］ Fitzpatrick M.The Retrospective Image Registration Evaluation Project Dataset［DB/OL］.［2014-08-01］.http：// www.insight-journal.org/rire/download-training-data.php.

［13］ Evans A C.BrainWeb：Simulated MRI Volumes for Normal Brain［DB/OL］.［2014-08-01］.http：//brainweb. bic.mni.mcgill.ca/brainweb/selection-normal.html.

［14］ Aja-Fernández S，Tristán-Vega A.Influence of Noise Correlation in Multiple-coil Statistical Models with Sum of Squares Reconstruction［J］.Magnetic Resonance in Medicine，2012，67（2）：580-585.

編輯 金胡考

Improvement of Interpolation Method in Cross Cumulative Residual Entropy-based Image Registration

XIANG Yan1，2，HE Jianfeng2，ZHANG Yunchun3，CAI Li3
（1.Key Lab of Computer Technologies Application of Yunnan Province，Kunming 650500，China；2.Institute of Information Engineering and Automation，Kunming University of Science and Technology，Kunming 650500，China；3.School of Software，Yunnan University，Kunming 650091，China）

The key strength of the Cross Cumulative Residual Entropy（CCRE）over the popular Mutual Information（M I）method is that the former has significantly larger noise immunity.But CCRE using conventional Partial Volume（PV）interpolation will result in the emergency of the local extremes on grid points，which may ham per the optimization algorithm from getting transformation parameters.In order to solve this problem，three improved PV interpolation methods are studied，including 3-order B-spline PV interpolation（BPV），Hanning windowed sinc PV interpolation（HPV）and Blackman-Harris windowed sinc PV interpolation（BHPV）.Meanwhile，a new interpolation method is proposed which uses flexible neighborhood center and makes the interpolation point to distribute the weight of the joint histogram to its adjacent 9 points.Moreover，it uses a Gaussian function as the PV interpolation kernel function to overcome the mutation of weight.Experimental result show s that the registration accuracy and speed of the proposed method is higher and faster，compared with BPV，HPV and BHPV method.So it is more suitable for CCRE computing.

image registration；Cross Cumulative Residual Entropy（CCRE）；Partial Volume（PV）interpolation；Gaussian function；3-order B-spline function；sinc function

1000-3428（2015）10-0199-05

A

TP391

相 艷，賀建峰，張?jiān)拼海?基于交叉累計(jì)剩余熵的圖像配準(zhǔn)中插值方法的改進(jìn)［J］.計(jì)算機(jī)工程，2015，41（10）：199-203.

英文引用格式：Xiang Yan，He Jianfeng，Zhang Yunchun，et al.Improvement of Interpolation Method in Cross Cumulative Residual Entropy-based Image Registration［J］.Computer Engineering，2015，41（10）：199-203.

國家自然科學(xué)基金資助項(xiàng)目（11265007）；云南省計(jì)算機(jī)技術(shù)應(yīng)用重點(diǎn)實(shí)驗(yàn)室開放基金資助項(xiàng)目（2403660120）。

相 艷（1979-），女，講師、碩士，主研方向：醫(yī)學(xué)圖像處理；賀建峰（通訊作者），教授、博士；張?jiān)拼?，講師、博士；蔡 莉，副教授、碩士。

2014-10-23

2014-11-24E-m ail：sharonxiang@126.com