改進的Powell醫學圖像配準算法的應用

邢正偉，李海瑛

昆明市延安醫院 設備科，云南 昆明 650051

引言

醫學圖像配準在醫學圖像處理領域中是一個重要的分支。尤其是當前，醫學成像技術的發展速度非常快，不管是在結構還是功能方面，都為臨床診斷提供了大量重要實用的影像數據。在臨床醫學中，單一模態的圖像往往不能提供醫生所需的足夠信息，通常需要配準并融合多模態圖像，幫助醫生了解病變組織或器官的情況，為臨床診斷和手術治療提供更全面和更準確的信息[1]。醫學圖像配準[2-3]就是通過尋求一種空間變換，使得兩幅醫學圖像結構上的對應點達成空間上的一致。在實際的臨床應用中，計算效率的高低往往是影響醫學圖像配準發揮作用的關鍵因素。

醫學圖像配準可以簡單的分為兩類：一類是基于圖像特征信息的配準，通過提取圖像的點、線、面來求取變換參數，達到配準的目的，該方法比較直接、計算量小、速度快，但是精度和魯棒性不好；另一類方法是基于圖像灰度信息的配準，該方法直接利用圖像的灰度進行配準，精度高、魯棒性強、不需要預處理，但是計算量大，速度慢，容易陷入局部極值。

將信息論的原理和方法應用在基于圖像灰度信息的配準是一個新的方向。而且，以互信息為相似性測度的配準已經成為一種主流的方法。醫學圖像配準的實質是配準函數的優化問題，即尋找互信息達到最大時，相關的幾個空間變換參數值，最優化方法在配準過程中具有非常重要的地位；但是配準函數經常不是光滑的，存在許多局部極值，這給求解帶來了很大難度，從而對配準結果有較大的影響[4]。目前使用基于互信息的醫學圖像配準的優化算法[5-8]有Powell算法、粒子群優化算法（Particle Swarm Optimization，PSO）及其改進算法、遺傳算法、單純形法等。其中使用較多的優化算法是Powell算法。Powell算法不需要計算導數，在每一維中使用Brent算法迭代搜索，搜索速度比較快，局部尋優能力極強，在局部搜索中精度要高于其他的優化算法[9]。另外，Powell算法在實現過程中一定程度上需要人工的干預。針對不同規格的醫學圖像，該算法需要專業人員來設定合理的搜索步長，以便獲得較好的配準結果，從而在效率、準確性和穩定性之間獲得良好的平衡效果。

本文在基于圖像灰度信息的配準問題研究上，開展了兩個方面的工作。第一，改進了歸一化互信息計算中聯合直方圖的統計方法，突出了感興趣區域在配準中的決定作用，使得互信息的大小能更準確的反映配準程度。第二，針對Powell優化算法對初始點過度依賴的問題，設計了一種改進算法來尋找配準的初始點，加快了尋優過程中的收斂速度，同時降低了誤配準發生的概率。

1 互信息配準的一般理論和方法

1.1 互信息的基本理論

對于在不同時間和不同條件下獲取的兩幅圖像R(x)和F(x)配準，定義一個相似性測度S，并尋找一個空間變換關系T，使得經過該空間變換后，兩幅圖像間的相似性測度達到最大，也就是圖像R上的每一個點在圖像F上都有唯一的點與之相對應，并且這兩點對應同一解剖位置。

互信息是信息理論中的一個基本概念，通常用于描述兩個系統間的統計相關性，或者是一個系統中所包含的另一個系統中信息的多少，它可以用熵來描述，熵表達的是一個系統的復雜性或者是不確定性[9]。

系統A的熵表示為：

系統B的熵表示為：

系統A，B的聯合熵表示為：

其中a∈A，b∈B。兩個系統的互信息可以表示為：

互信息可以作為相似性測度來完成醫學圖像的配準操作。一般用聯合概率分布pAB(a,b)和完全獨立時的概率分布p(a)A×pB(b)間的廣義距離來估計互信息[10]。

其中聯合概率分布pAB(a,b)用兩幅圖像重疊區域的歸一化聯合直方圖來計算：

其中邊緣概率分布pA(a)、pB(b)的計算方法為：

傳統互信息對兩幅圖像重疊程度的變化非常敏感，使得配準過程的互信息變化曲線不光滑，因此本文采用歸一化互信息（Normalized Mutual Information，NMI）作為配準的相似性測度，提高配準的穩定性。由式(9)可知，NMI的取值范圍是[0, 2]。NMI值的大小反映了配準程度的好壞，當兩幅圖像完全配準時，NMI達到最大（NMI=2）。NMI(A,B)=[H(A)+H(B)]/H(A,B) (9)

1.2 Powell算法的原理

在醫學圖像配準領域中最常用的優化算法是Powell算法，該算法是一種局部多參數局部最優化算法，相對于全局優化算法來說有較快的收斂速度，但對于大規格的圖像來說，仍有提升的空間；該算法容易受到局部極值的影響陷入誤配準，其改進的關鍵在于初始點的選擇。在每一次迭代中，都要從初始點開始進行一維搜索，搜索的方向是Powell算法中的另一個重要參數，用一個多方向集C來表示。當參數的導數不易計算時，應用多方向集算法[11-13]這個概念，就可以把多元函數的求解問題轉化為求取一維極值的問題。

Powell算法具體的優化過程分為若干次迭代，每次迭代都由N+1次一維搜索組成，其中N是搜索空間的維度，例如二維圖像剛性變換配準的搜索空間維度是3，分別是水平方向、垂直方向和旋轉方向。算法每次迭代過程的步驟是，從初始點開始，依次沿方向集的N個方向進行一維搜索，得到一個最佳值及最佳值對應的點。然后沿該點和初始點的連線的方向進行搜索，求得本次迭代的最佳值及最佳值對應的點，然后把連線的方向替換前N個方向中效果不好的方向，形成新的方向集，進行下一次迭代。

具體構造過程如下：① 將方向集ci初始化：ci=e（i=1,2,…,N，e為坐標軸的單位向量）；② 記錄初始點位置P0；③ 從Pi出發，依次沿方向集的各個方向尋找該方向歸一化互信息的極大值點，記該點為Pi+1，以該點為新的出發點，進行下一次迭代搜索；④ 重復步驟（3），直到新找到的點和前一點之間的距離＜ε(ε為設定精度值)，則跳出循環，把該點作為算法的解。

2 算法的改進

2.1 歸一化互信息的改進

使用改進算法進行初始點尋優的時候，對于大規格的圖像會存在精度不準確的情況，因為圖像的背景區域太大，而有效信息區域在整幅圖像中的比例過小。針對這種情況需要對聯合直方圖的計算方法進行改進，以提高改進算法的穩定性。

在計算聯合直方圖的時候需要用到插值算法（圖1），因為當浮動圖像經過變換到參考圖像的坐標網格中時，像素點所在的位置不一定是和參考圖像的像素點重合的。這就需要我們使用插值算法來計算變換后點的像素值。現代配準算法一般采用的是部分體積插值算法（Partial Volume Interpolation，PV），我們也采用PV算法來計算像素值。

如圖1所示，在每一次計算互信息值大小的時候，根據變換參數，變換T都會把浮動圖像F上的點Pf（i,j）變換到參考圖像R上的點Pr，坐標為（x+Δx,y+Δy）。如果要用一般的PV算法來計算點Pr的灰度值，引進了圖像原本不存在的灰度值，計算量太過龐大，在實際的編程中，聯合直方圖采用如下公式的計算方法。

圖1 PV插值算法示意圖

如果要放大感興趣部分在整幅圖像中的權重，則需要把hist中背景部分的計算量縮小，具體的計算公式如下：

其中A，B分別是參考圖像的長和寬，a，b分別是剔除背景后的圖像在圖像坐標系中水平方向和垂直方向的最大和最小坐標。這樣就降低了背景區域在聯合直方圖中的作用，避免了因為背景重合的局部極值造成的誤配準，也增加了尋找初始點算法的精確度和穩定性。

歸一化互信息的各部分分量在配準過程中的變化曲線，見圖2。

圖2 互信息中各分量的變換曲線

如圖2所示，在計算聯合直方圖的時候，a代表參考圖像，b代表浮動圖像。浮動圖像b中的每個像素點都被計算在內，對直方圖的貢獻權值為1，故圖像b的信息熵Hb是一常量，其變化曲線是水平的。而參考圖像a的部分像素點沒有被計入聯合直方圖的計算，其信息熵Ha為一變量。因此在計算歸一化互信息的時候可以保留參考圖像的信息熵，舍棄浮動圖像的信息熵，來簡化歸一化互信息的計算，同時保留與原來相似性測度相同的收斂特性。改進的歸一化互信息計算公式如下所示：

2.2 尋優算法的改進

要想節省優化過程的時間，根本在于減小迭代的次數。以下以二維圖像剛性變換配準為例，分別對3個搜索方向的互信息變換曲線進行研究。

選取一幅圖像作為參考圖像，以參考圖像水平方向向左平移5個像素，垂直方向向上平移13個像素，順時針旋轉10°得到的圖像為浮動圖像，以參考圖像和浮動圖像對齊時默認位置為初始點。把浮動圖像以圖像中心為原點，向左平移20個像素至向右平移20個像素，求取平移后圖像與原參考圖像的互信息，得到在該方向上的互信息測度變化曲線，見圖3a。在垂直方向上下移動20個像素和旋轉方向逆時針至順時針旋轉10°得到的互信息測度變換曲線，見圖3b和圖3c。

圖3 待配準的兩幅圖像相對移動得到的歸一化互信息曲線

由圖3可見該曲線在預設的參數值附近達到最大值，說明完全可以設計一種簡單的算法來求取該值的大小作為Powell算法的初始點，避免Powell算法前期漫長的尋優過程和大量的計算，加快算法的收斂速度，提高圖像配準的速度。

平移變換初始點的計算方法如下：先計算參與配準的兩幅圖像在初始位置V0（0,0,0）的互信息值F0。因為平移變化較旋轉變化的互信息曲線更加光滑，故先從水平方向進行計算。計算初始位置向左和向右平移一個像素單位的互信息值，選取互信息值增加的一側為尋優方向，不斷搜索，直到互信息值開始減小。如果前進一個像素單位互信息值增加，則繼續搜索，如果互信息值繼續減小，則跳出循環，返回前一個點對應的坐標x0作為新的初始點V0（x0,0,0）。從新的初始點V0（x0,0,0）開始，在垂直方向按照同樣的方法尋找極值點y0。由于在旋轉方向互信息曲線的局部極值較多，故使用上面的方法來求取極值點容易陷入局部極值。旋轉角度都有一定的限度，故對于旋轉方向的極值點求取要采取新的方法。計算[-180, 180]范圍內每個整數單位處的互信息值，求其中的最大值，作為旋轉方向的極值點θ0。最終得到點V0（x0,y0,z0）作為下一步Powell算法的初始點。

具體的實現步驟如下：

（1）默認參與配準的兩幅圖像初始位置的點V0（0,0,0）為初始點，令x0=0，y0=0，z0=0，計算該點的互信息值F(0,0,0)。

（2）計算浮動圖像在水平方向平移+1，-1個像素的互信 息 值 F(+1,0,0)和 F(-1,0,0)。 如 果 F(+1,0,0)-F(0,0,0)＞0，則令e=1，否則，令e=-1，此時x0=e。

（3）若 F(x0+e,0,0)-F(x0,0,0)＞0，則 x0=x0+e，重復步驟（3）。

（4） 當 F(x0+e,0,0)-F(x0,0,0)＜0時， 如 果 F(x0+2e,0,0)-F(x0,0,0)＞0，則x0=x0+e，重復步驟（3）；如果F(x0+2e,0,0)-F(x0,0,0)＜0，則跳出循環，返回水平方向的初始值x0。

（5）以點V0（x0,0,0）為初始點，計算該點的互信息值F(x0,0,0)。計算浮動圖像在垂直方向平移+1，-1個單位后的互信息值F(x0,+1,0)和F(x0,-1,0)。如果F(x0,+1,0)-F(x0,0,0)＞0，則令e=1，否則，令e=-1，此時y0=e。

（6）若 F(x0,y0+e,0)-F(x0,y0,0)＞0，則 y0=y0+e，重復步驟（6）。

（7）當F(x0,y0+e,0)-F(x0,y0,0)＜0時，如果F(x0,y0+2e,0)-F(x0,y0,0)＞0，則 y0=y0+e，重復步驟（6）；如果 F(x0,y0+2e,0)-F(x0,y0,0)＜0，則跳出循環，返回垂直方向的初始值y0。

（8）以點 V0（x0,y0,0）為初始點，計算 z0在 [-180，180]范圍內每個整數點處的互信息值Fi，求得其中的最大值Fmax以及最大值對應的點θ0。

（9）以點 V0（x0,y0,θ0）作為Powell算法的初始點。

該算法不僅能減少優化算法的迭代時間，而且對于角度旋轉方向來說，可以越過大部分的局部極值，見圖4，可以看出互信息曲線在角度不斷變化的情況下有許多個局部極值，而本算法在全局范圍內的搜索，跳出了這些局部極值的干擾，能有效的計算出最大值在10°左右。

3 實驗數據及結果

為驗證算法的有效性和穩定性，本文采用了不同的醫學圖像進行驗證。

圖4 旋轉一周的互信息曲線

實驗一采用了Brain Web的腦部磁共振T1加權圖像[14-15]，圖像大小為181×217×181。任意選取其中的一幀作為參考圖像，見圖5a。參考圖像在水平方向、垂直方向和旋轉方向做不同的參數變化，得到的圖像作為浮動圖像，見圖5b。圖5b是圖5a向右平移10個像素，向下平移12個像素，逆時針旋轉10°得到的。

圖5 實驗一所用圖像

制作了5組浮動圖像，每組浮動圖像都是參考圖像在水平、垂直方向平移一定的距離，并旋轉一定的角度形成的。實驗分別采用經典Powell算法，PSO算法（以傳統歸一化互信息為相似性測度）和本文的改進算法3種算法對實驗圖像進行配準操作，每次配準分別獨立運行10次取平均值。

具體的配準結果，見表1。其中x、y、ang分別是圖像在水平，垂直和旋轉方向的變換參數，x、y的單位是像素，ang的單位是度（°）。變換參數是（20,20,1）時，經典Powell算法出現了配準失敗的情況，而改進算法得到的變換參數和預設參數的誤差在亞像素的精度范圍內，并未出現配準失敗的情況，比經典Powell算法更穩定。在速度方面，改進的算法也比傳統的算法快很多。

實驗二采用的圖像是Vanderbilt大學提供的病人的CT圖像，圖像大小為512×512×29[15]。任意選取其中的一幀作為參考圖像[16]，見圖6a。參考圖像在水平方向，垂直方向和旋轉方向做不同的參數變化，得到的圖像作為浮動圖像，見圖6b。圖6b是圖6a向右平移10個像素，向下平移12個像素，順時針旋轉10°得到的。

制作了5組浮動圖像，每組浮動圖像都是參考圖像在水平、垂直方向平移一定的距離，并旋轉一定的角度形成的。實驗分別采用經典Powell算法、PSO算法（以傳統歸一化互信息為相似性測度）和本文的改進算法3種算法對實驗圖像進行配準操作，每次配準分別獨立運行10次取平均值。具體的配準結果，見表2。

表1 實驗一的配準結果

圖6 實驗二所用圖像

表2 實驗二的配準結果

由表2可以看出，本文的改進算法在精度方面，得到的變換參數和預設參數的誤差在亞像素的精度范圍內，滿足配準的理論要求和臨床需要；在速度方面，改進的算法比PSO算法至少快一倍以上，比傳統的Powell算法也要快的多。

4 結論

本文針對傳統互信息圖像配準中Powell算法的自動化程度不高、配準時間長和易受局部極值影響而出現誤配準的問題從兩個方面進行了改進。一是改進了Powell算法中初始點的尋找方式，節省了優化初期過多的迭代量，減少了配準所需要的時間，同時提高了配準的自動化程度。二是改進了歸一化互信息的計算公式及其聯合直方圖的計算方法，突出了感興趣區域在相似性測度中的作用，使得互信息能更準確的反應配準的程度，減小誤配準的幾率。實驗證明了本方法的有效性。

[參考文獻]

[1] 鄭瑩,李光耀.區域和局部信息結合的雙向醫學圖像配準[J].中國圖象圖形學報,2011,16(1):90-96.

[2] 張汗靈,楊帆.基于互信息和混合優化算法的多模醫學圖像配準[J].湖南大學學報:自然科學版,2006,33(1):117-120.

[3] 楊日容.基于遺傳算法及最大互信息的醫學圖像配準研究[D].昆明:昆明理工大學,2005.

[4] 方偉,孫俊,丁彥蕊,等.醫學圖像配準的混合量子粒子群優化算法研究[J].計算機工程與應用,2011,47(3):166-169.

[5] 劉斌,彭嘉雄.圖像配準的小波分解方法[J].計算機輔助設計與圖形學學報,2003,15(9):1070-1073.

[6] Pardalos PM,Shalloway D,Xue G.Optimization methods for computing global minima of nonconvex potential energy functions[J].J Global Optim,1994,4(2):117-133.

[7] Nelder JA,Mead R.A simplex method for function minimization[J].Comput J,1965,7(4):308-313.

[8] Wilkinson B,Allen M.Parallel Programming: Techniques and Applications Using Networked Workstations and Parallel Computers[M].北京:高等教育出版社,2002.

[9] 方偉,孫俊,丁彥蕊,等.醫學圖像配準的混合量子粒子群優化算法研究[J].計算機工程與應用,2011,47(3):166-169.

[10] Pluim JPW,Maintz JBA,Viergever MA.Mutual-informationbased registration of medical images: a survey[J].IEEE T Med Imaging,2003,22(8):986-1004.

[11] 袁偉,時公濤,蔣詠梅.一種基于ROEWA算子和GA-Powell算法的SAR圖像配準方法[A].第十四屆全國信號處理學術年會(CCSP-2009)論文集[C].2009.

[12] Ding H,Bian ZH.A sub-pixel registration approach based on powell algorithm[J].Acta Photonica Sinica,2009,38(12):46-49.

[13] 劉葉青,劉三陽,谷明濤.Powell算法在線性支持向量機中的應用[J].計算機工程,2011,37(12):161-163.

[14] 馮林,管慧娟,滕弘飛.基于互信息的醫學圖像配準技術研究進展[J].生物醫學工程學雜志,2005,22(5):1078-1081.

[15] 彭景林,章兢,李樹濤.基于改進PV插值和混合優化算法的醫學圖像配準[J].電子學報,2006,34(5):962-965.

[16] 馮林,嚴亮,黃德根,等.PSO和Powell混合算法在醫學圖像配準中的應用研究[J].北京生物醫學工程,2005,24(1):8-12.