Mean Shift的漸進無偏變換圖像配準

楊 烜

(深圳大學計算機與軟件學院 深圳 518060)

1 引言

圖像配準的目的是尋找兩幅圖像之間的變換函數，使得圖像中相應目標的位置對準[1]。基于標志點對應關系的圖像配準方法已得到廣泛使用[2-6]，這類方法需要根據對應點對構造變換函數，變換函數具有明確的解析式，易于分析，計算簡單。但這類變換函數往往是不可逆的，只能確保標志點之間的對應關系，不能保證圖像間其他位置的對應關系，即變換的無偏性較差。

如果變換函數及其逆變換的Jacobia在1附近對稱分布，稱為無偏變換[7]。無偏變換對于圖像配準非常重要，它意味著無論圖像A向圖像B配準，或是圖像B向圖像A配準都可以使同一目標得到配準。更重要的是，無偏變換是拓撲保持的，可以使圖像在形變前后的拓撲關系不發生變化，這是大多數圖像配準方法努力實現的目標[8-10]。目前關于無偏變換的研究比較少，文獻[11]定義的對稱代價函數，使源圖像和目標圖像互換后的變換函數互逆，但該方法采用FFD形變模型，計算量較大；文獻[4]提出了標志點的一致性配準方法，該方法可以使標志點在正、反變換達到一致，但會導致標志點對應關系發生較大誤差，還需要引入一致性強度配準進一步調整；文獻[7]定義了變換函數與理想變換之間的對稱KL距離，并定義了對稱目標函數尋找無偏變換函數，但文獻[7]采用的是粘性流體配準模型，計算量仍然很大。

本文針對基于標志點對應關系的圖像配準問題，提出了一種基于無偏變換的漸進式圖像配準方法，該方法只需要很少的幾個初始標志點，通過不斷加入新的標志點逐漸提高配準精度。本文從理論上分析了選擇標志點的原則，利用Mean Shift迭代快速尋找對應標志點，構造無偏的變換函數。本文方法所構造的圖像變換函數既具有較好的無偏特性，同時又保證了標志點之間的對應關系。

2 標志點分布對sKL距離的影響

無偏變換要求變換函數及其逆變換使圖像被壓縮和被膨脹程度是一致的，我們采用文獻[7]提出的對稱Kullback-Leibler(sKL)距離來描述變換的無偏性。假設變換函數為h，逆函數為h-1,h和h-1的Jacobi行列式分別為定義3種概 率 密 度 函 數 pdfh(ξ)=|Dh(ξ)|,pdfh-1(ξ)=|Dh-1(ξ)|,pdfid(ξ)=1，其中pdfid是單位映射(id(x)的概率密度函數。對稱Kullback-Leibler(sKL)距離定義為

其中KL是Kullback-Leibler距離。尋找無偏的變換函數，就是使變換函數的sKL距離盡可能小。對于基于標志點對應關系的配準問題而言，就是尋找標志點對，在配準精度提高的同時，使變換函數的sKL距離盡可能小。

假設源、目標點集分別為P和Q，有N對標志點(pi,qi),i=1,2,…,N。pi∈P,qi∈Q，利用徑向基函數構造變換函數，

其中x是變換空間中的點，R(r)是徑向基函數，ri=||x-pi||，系數wi,i=1,…,N滿足方程KW=Q-P,K是N×N矩陣，Kij=R(||pi-pj||)。徑向基函數構造的h不能保證可逆，即使其逆函數h-1存在，h-1與h的壓縮和膨脹程度也有較大差異，即無偏性較差。

下面我們將在一個簡單2維模型上，分析添加標志點對變換函數sKL的影響。假設當前只有一個源標志點p向右水平移動到目標點q，水平位移量為d，變換函數h(x)及其Jacobi行列式Dh(x)為

當基函數是單調減的緊支撐徑向基函數時，滿足 ?R(r)/ ?r＜ 0 。不失一般性地以p為原點，則x＜ 0 的點(p左側區域)Jacobi值大于1，而x＞0(p右側區域)Jacobi值小于1。在點p的垂直上方L處添加點pa，該點水平位移至qa，假設位移量為dR( | |pa-p| | )+ε，其中dR( ||pa-p| |)是根據原變換函數h(x)對pa的位移量，ε是該位移的小擾動。則新的變換函數ha(x)為

當L相對基函數的支撐集較小時，式(5)可以近似為

我們以緊支撐薄板樣條(CSTPS)[5]作為徑向基函數進行分析，可以推出

式(7)表明，新增點可以使x＞0的點Jacobi值增加，x＜0的 Jacobi值減小，這意味著通過添加點，在點p兩側的Jacobi值更接近1，即變換函數的無偏性改善了。需要說明的是，擾動ε不能過大，否則調整量會過大，可能會使變換函數的無偏性變差。上述分析中對于新增標志點pa的位置并沒有嚴格的限制，一般地講，我們需要在原標志點附近選擇新增點，并且新增點的位移與原變換函數產生的位移是一致，即新加點的對應點選擇在原變換函數的映射位置附近，這樣該點的位移趨勢與原變換函數近似，以滿足ε較小的條件。另一方面，L較小也意味著將在較小的范圍內搜索對應標志點，有利于尋找到準確的對應點。

對于圖像配準問題，上述的要求可以這樣理解：當原始變換函數的無偏性需要改善時，可以在原標志點的附近添加新增標志點，新增標志點的選擇原則是以圖像配準精度得到提高為依據，這樣新增標志點可以對變換函數的無偏性起到改善作用。描述得更形象一些，當形變曲面被某些位移撕扯時，由于這種撕扯的受力不夠均勻，會使圖像在這些撕扯位置附近出現過度擠壓或膨脹的情況。通過增加標志點，可以使撕扯效果得到緩解，變得平滑和均勻，變換函數的無偏性得到改善。

3 Mean Shift確定新增標志點對

假設S,R分別是源圖像和參考圖像，當前源標志點集為P，目標標志點集為Q，標志點對集合為?pi∈P，選擇新增源標志點滿足

其中 L oc(I,x)是圖像I中以x為中心的局部區域，E(S)是源圖像S的邊緣點集合。新增源標志點處于邊緣上，其所處的區域與參考圖像相差較大。同時該點與當前標志點的距離不能太遠，這樣可以利用當前標志點提供一個較理想的初始位置，有利于Mean Shift迭代收斂；同時也不能距離過近，以免局部形變受力不勻。

其中k(x)是核函數，分別是以為中心的鄰域點，是歸一距離，C和Cq是歸一系數，u是灰度量化值，b(x)是量化函數，本文采用直方圖模糊約束聚類方法進行量化。根據Mean Shift迭代，更新為

由于Mean Shift收斂是基于概率密度函數變化最大方向進行的，容易收斂于局部最優。需要對Mean Shift進行幾點改進。首先進行多尺度配準，將配準圖像下采樣到大尺度圖像上進行Mean Shift收斂，可以減小搜索范圍。將大尺度圖像的標志點對映射到小尺度圖像中，得到粗配準結果。然后在小尺度圖像中重復這個配準過程，得到精度更好的配準結果。

其次，由于Mean Shift統計的概率密度是不反映圖像空間位置信息的，Mean Shift收斂后對應的點，可能并不是使局部圖像空間對應關系最優的點。本文在Mean Shift迭代的路徑中計算各點的局部配準度量，選擇局部配準度量最大的點作為收斂結果，以提高對應點精度。另外，Mean Shift迭代路徑要進行限制，以避免沿某個錯誤方向迭代過遠。

在一個尺度下的完整配準過程如下：

步驟 1 輸入源圖像和參考圖像，初始標志點對集合Ω，設置參數 d istThmin,Th，支撐集c。

步驟 2 利用式(2)計算Ω的變換函數h(x)。

步驟 3 如果 d istTmin＞T h，則算法退出，否則

步驟 4 ?pi∈P，確定滿足式(8)的新增源標志點。

步驟 7 計算Ω∪Ωa的變換函數ha(x)和，以及代價函數。

步驟 9 如果集合Ω沒有增大，則 d istTmin=distTmin+ 1 。轉至步驟3。

4 實驗結果

分別針對大形變和小形變配準進行實驗驗證。圖1中第1、第2行是Rose和Clown的大形變配準，分別選擇了4個和7個初始標志點對，采用CSTPS作為徑向基函數，支撐集c=200,d istTmin= 8,Th=16。圖像大小均為 256×256，將圖像下采樣為128×128，進行大尺度配準，再進行小尺度配準。其中Rose圖、Clown圖中各添加了142, 92對標志點，可以看出， Rose的邊界部分和Clown的帽沿這些形變劇烈的地方都得到了較好的配準。

為了比較配準效果，本文同時通過SIFT特征[13]確定對應標志點，將位移大于32的嚴重誤配標志點對剔除。SIFT標志點對應關系的閾值設為使變換函數Jacobi值大于0的最大值，以盡可能多地保留標志點對。SIFT在Rose圖和Clown圖中分別添加了42和31對標志點，SIFT在Clown中尋找的標志點對存在誤差，無法得到Jacobi值大于0的變換結果。這表明 SIFT不適用于形變比較劇烈的彈性配準問題，而本文方法具有較好的魯棒性。

圖1 將源圖像配準到參考圖像的實驗結果(從左到右：待配圖像，參考圖像，本文方法配準結果與參考圖像的差，SIFT配準結果與參考圖像的差)

小形變配準采用了人體腦部掃描切片，待配準圖像源于BrainWeb的subject04個體的腦部體數據切片，分別是第0和第50切片，經過人工小形變，得到參考圖像。第0和第50切片中各選擇了5對初始標志點。這兩幅配準圖像的細節比較豐富，特別是腦切片中部區域對配準要求較高，本文方法在第0，第50切片分別增加了110, 123對標志點，并且在腦中部區域確定了多個精度較高的標志點，提高了這些區域的配準精度，而SIFT在第0，第50切片中各新增65和308對標志點，但這些標志點中有許多點幾乎重合，而在腦切片中部區域選擇的標志點較少，分布也不均勻，導致這些區域配準精度不高。表1列出了本文方法和SIFT配準結果的互信息、均方差(MSD)、相關系數(CC)和sKL距離，加粗的數值表示最優結果。從表1可以看出，對于小形變配準問題，本文配準結果同樣優于 SIFT方法的結果，同時sKL距離表明，本文方法得到的變換函數的無偏性優于SIFT方法。

5 結論

本文針對標志點對應關系的圖像配準問題，提出了一種漸進式無偏變換的配準方法。首先人為選擇數量極少的幾個初始標志點，在初始標志點附近添加適當的標志點，利用Mean Shift迭代尋找相應的對應點；根據新標志點對判定配準度量的變化，確定需要添加到標志點集合的新標志點對；在擴展的標志點集合中繼續利用Mean Shift迭代尋找對應標志點，逐步達到配準的效果。該方法中尋找對應標志點的過程充分考慮了變換無偏性的影響，同時為Mean Shift收斂提供了較準確的初始位置，多尺度的配準策略有利于魯棒地尋找到準確的對應點。實驗比較結果表明本文方法不僅適用于大形變配準問題，也可以用于解決小形變問題，是一種魯棒的配準方法。

表1 配準結果比較(N/A表示變換函數拓撲不能保持)

[1]Rueckert D and Aljabar P. Nonrigid registration of medical images: theory, methods, and applications.IEEE Signal Processing Magazine, 2010, 27(4): 113-119.

[2]Myronenko A and Song X B.Intensity-based image registration by minimizing residual complexity.IEEE Transactions on Medical Imaging, 2010, 29(11): 1882-1891.

[3]Chen J, Tian J, Lee N,et al.. A partial intensity invariant feature descriptor for multimodal retinal image registration.IEEE Transactions on Biomedical Engineering, 2010, 57(7):1707-1718.

[4]Johnson H J and Christensen G E. Consistent landmark and intensity-based image registration.IEEE Transactions on Medical Imaging, 2002, 21(5): 450-461.

[5]Zhang Z X and Yang X. Elastic image warping using a new radial basic function with compact support. Proceedings of International Conference on Bio-inspired Systems and Signal Processing 2008, Funchal, Madeira, Portugal, January 28-31,2008: 216-219.

[6]Shen D and Davatzikos C. HAMMER: hierarchical attribute matching mechanism for elastic registration.IEEE Transactions on Medical Imaging, 2002, 21(11): 1421-1439.

[7]Leow A D, Yanovsky I, Chiang M C,et al.. Statistical properties of Jacobian maps and the realization of unbiased large-deformation nonlinear image registration, 2007, 26(6):822-832.

[8]Lin X B, Qiu T S, Morain-Nicolier F,et al.. A topology preserving non-rigid registration algorithm with integration shape knowledge to segment brain subcortical structures from MRI images.Pattern Recognition, 2010, 43(7):2418-2427.

[9]Sdika M. A fast nonrigid image registration with constraints on the Jacobian using large scale constrained optimization.IEEE Transactions on Medical Imaging, 2008, 27(2): 271-281.

[10]Noblet V, Heinrich C, Heitz F,et al.. 3-D deformable image registration: a topology preservation scheme based on hierarchical deformation models and interval analysis optimization.IEEE Transactions on Image Processing, 2005,14(5): 553-566.

[11]Gholipour A, Kehtarnavaz N, Yousefi S,et al.. Symmetric deformable image registration via optimization of information theoretic measures.Image and Vision Computing,2010, 28(6): 965-975.

[12]Comaniciu D, Ramesh V, and Meer P. Real-time tracking of non-rigid objects using mean shift. Proceedings of IEEE Conference Computer Vision and Pattern Recognition, Los Alamitos: IEEE Computer Society Press, 2000: 142-149.

[13]Lowe D G. Distinctive image features from scale-invariant keypoints.International Journal of Computer Vision, 2004,60(2): 91-110.