彭磊 楊秀云 張裕飛 李光耀

摘 要：非剛性點集配準算法中，能否找到正確的對應關系對配準結果起著至關重要的作用，而通常兩個點集中的對應點除了距離比較接近之外還具有相似的鄰域結構，因此提出基于全局與局部相似性測度的非剛性點集配準算法。首先，使用一致性點漂移（CPD）算法作為配準框架，采用高斯混合模型對點集進行建模。然后，對全局局部混合距離進行改進，形成全局與局部相似性測度準則。最后，采用期望最大化（EM）算法迭代地求解對應關系和變換公式：在迭代初期局部相似性所占比重較大，從而能夠盡快地找到正確的對應關系;隨著迭代的進展全局相似性比重逐漸增大，從而確保得到較小的配準誤差。實驗結果表明，與薄板樣條魯棒點匹配（TPS-RPM）算法、高斯混合模型點集配準（GMMREG）算法、基于L2E估計的魯棒點匹配算法（RPM-L2E）、基于全局局部混合距離與薄板樣條的點集配準算法（GLMDTPS）和CPD算法相比，所提算法的均方根誤差（RMSE）分別下降了39.93%、42.45%、32.51%、22.36%和11.76%，說明該算法具有較好的配準效果。

關鍵詞：點集配準;圖像配準;非剛性配準;高斯混合模型;相似性測度

中圖分類號：TP391.41

文獻標志碼：A

Abstract： In the non-rigid point set registration algorithm， whether the correct correspondence can be found plays an important role. Generally the corresponding points in two point sets have similar neighborhood structures besides the close distance. Therefore， a non-rigid point set registration algorithm based on global and local similarity measurement was proposed. Firstly， the Coherent Point Drift （CPD） algorithm was used as the registration framework， and the Gaussian mixture model was used to model the point sets.Secondly， the global and local mixture distancewas improved to form the global and local similarity measurement criterion.Finally， the correspondence and the transformation formulawere solved by the Expectation Maximization （EM） algorithm. In the initial stage of the iteration， the proportion of local similarity was larger so that the correct correspondence was able to be found rapidly; with the progress of the iteration， the proportion of global similarity was increased to ensure the smaller registration error. Experimental? results show that compared with the Thin Plate Spline Robust Point Matching （TPS-RPM） algorithm，the Gaussian Mixture Models point set REGistration （GMMREG） algorithm，the Robust Point Matching algorithm based on L2E estimation （RPM-L2E）， the Global and Local Mixture Distance and Thin Plate Spline based point set registration algorithm （GLMDTPS） and the CPD algorithm， the proposed algorithm has the Root Mean Squared Error （RMSE） decreased by 39.93%， 42.45%， 32.51%， 22.36% and 11.76% respectively， indicating the proposed algorithm has better registration performance.Key words：? point set registration; image registration; non-rigid registration; Gaussian mixture model; similarity measurement

0 引言

點集配準是一個基礎而關鍵的問題，廣泛應用在計算機視覺[1]、醫(yī)學圖像分析[2-4]、遙感圖像處理[5]和模式識別[6]等領域。其目標是尋找兩組給定點集之間的對應關系，通過空間變換將一個點集對齊到另一個點集。根據(jù)變換方式的不同可分為剛性配準和非剛性配準。剛性點集配準的變換方式包括放縮、平移、旋轉(zhuǎn)，僅涉及少量的參數(shù)。非剛性點集配準的變換方式往往是未知的、復雜的、難以建模的，通常包含較多的變換參數(shù)。

迭代最近點算法[7]是著名的點集配準方法之一，主要用于自由曲線和曲面的配準，它假設兩個點集間距離最近的點為對應點，通過迭代地最小化均方距離來達到目標函數(shù)的收斂。Chui等[8]提出了一個非剛性點匹配的通用框架，使用薄板樣條作為空間變換模型，并通過軟分配和確定性退火技術聯(lián)合求解優(yōu)化問題。一致性點漂移（Coherent Point Drift， CPD）算法[9]將配準兩個點集考慮為概率密度估計問題，使用運動一致性理論進行約束，通過迭代方法將高斯混合模型（Gaussian Mixture Model， GMM）[9-10]的質(zhì)心擬合到數(shù)據(jù)點。Jian等[11]使用兩個GMM來表示點集，并通過兩個GMM之間的差異最小化來解決點集配準問題。文獻[12-13]采用高斯域的方法實現(xiàn)了異常點移除和點集配準。這些方法主要使用點之間的全局距離來確定對應關系，基本沒有考慮點的鄰域結構信息。

Belongie等[6]引入了形狀上下文（Shape Context，SC）的概念，用來測量兩個形狀之間的相似度。Ma等[14-15]首先使用SC描述符建立粗糙的對應關系，然后用L2E估計器來估算變換。Yang等[16]采用魯棒的全局和局部混合距離（Global and Local Mixture Distance，GLMD）測度準則進行非剛性點集配準。雖然這些算法使用點集的結構信息作為相似性測量，但是它們在查找對應關系和估計變換模型時是分別處理的。

為了充分利用點之間的距離和點的鄰域結構信息，快速找到正確的對應關系，本文以CPD算法[9]為框架，采用GMM表示點集，對GLMD[16]進行改進，并融入到GMM中，形成了全局與局部相似性測度準則。GMM的高斯分量反映了全局相似性測度準則，而每個高斯分量使用不同的比例系數(shù)，則反映了點的局部相似性測度準則。通過權重系數(shù)調(diào)節(jié)二者的比重。采用期望最大化（Expectation Maximization， EM）算法迭代地求解對應關系和變換公式。在迭代初期，局部相似性所占比重較大，從而盡快找到正確的對應關系;迭代過程中逐漸提高全局相似性比重，最終得到更小的配準誤差。該算法在保證快速找到點的正確對應關系的同時，能夠得到較好的配準效果。

1 一致性點漂移算法

設D維（通常D=2或3）空間的兩個點集X和Y，浮動點集為X={xm|m=1，2，…，M}，參考點集為Y={yn|n=1，2，…，N}。非剛性點集配準算法的目標是根據(jù)相似性測度準則找到X和Y之間點的對應關系，并估算一個非剛性變換公式T，使X通過變換后對齊到Y。

CPD算法將兩個點集的配準看成是概率密度估計問題。使用GMM對問題進行建模[9]，把X中的點看成GMM的質(zhì)心，Y中的點看成是由GMM產(chǎn)生的數(shù)據(jù)點。方法的核心是保持點集的拓撲結構，迫使GMM的質(zhì)心作為一個整體一致地向參考點集移動。在優(yōu)化過程中，浮動點集通過變換T與參考點集逐漸對齊，并可以根據(jù)后驗概率獲得點的對應關系。GMM的概率密度函數(shù)為：

式（1）中的均勻分布p（yn|M+1）=1/N，用來處理點集中存在的噪聲和異常點，所占比重為ω（0≤ω≤1）。其他所有高斯組成部分所占比重為1-ω，其中每個高斯分量使用相同的比例系數(shù)ηmn=1/M。

為了使變換比較平滑，采用Tikhonov正則化框架[9]。定義再生核希爾伯特空間的正則項λφ（T）/2，其中φ（T）為一個平滑函數(shù)，λ為正則項的權重系數(shù)。使用θ表示參數(shù)集合，通過極大似然估計來求得參數(shù)。將負對數(shù)似然函數(shù)作為CPD算法的能量函數(shù)，如式（3）所示：

2 全局和局部混合距離

2.1 全局距離

浮動點集上一點xm與參考點集上一點yn之間的全局距離[16]定義為其歐氏距離的平方，如式（4）所示：

Gmn可以描述兩個點集之間的全局性結構差異。通常兩個點集基本對齊后，距離較近的點很可能就是對應點。

2.2 局部距離

某個點q的鄰域定義為其所在點集中距離q最近的K個點。用φ（q）k來表示第k個最近點。局部距離[16]用于度量兩個點xm與yn的鄰域結構相似性，定義為：

其中ψ是位移函數(shù)，表示點xm及其鄰域點整體向點yn移動直至xm與yn點重合，其定義為：

3 本文算法

3.1 局部相似性的改進

GLMD算法中局部距離是通過計算兩個點的鄰域中具有相同編號k的點對之間距離的和得到的。本文將式（5）中φ（yn）的下標由k改為f（k），即尋找兩個鄰域中兩兩點對之間距離之和最小的值作為局部相似性。式（5）改寫為：

f為兩個鄰域中點的某種一一映射關系，使得兩個鄰域（分別包含K個點）之間K個點對的距離之和達到最小，即能夠使Lmn得到最小值。相當于二分圖的最大匹配問題，可以通過匈牙利算法來求解f。

如圖1所示，兩個點p和q，各自的鄰域內(nèi)分別有編號為1、2、3、4的點，虛線表示兩個鄰域之間點的對應關系。通過計算對應點間距離之和得到局部距離。GLMD按照編號確定對應點，有可能產(chǎn)生錯誤的結果，從而得到錯誤的局部距離，如圖1（a）所示。本文算法能夠找到一個最佳匹配來確定對應點，因此能夠求得正確的局部距離，同時所求局部距離也是最小的，如圖1（b）所示。

3.2 全局局部相似性測度

CPD算法利用了點之間的全局距離，即認為全局距離較近的點很可能為對應點。因此，式（2）可以寫為：

CPD算法中每個高斯分量使用相同的比例系數(shù)ηmn=1/M，沒有考慮點的鄰域結構。然而兩個點的鄰域結構越相似，即局部相似性越高，它們成為對應點的可能性也越大。為了充分利用點的局部相似性，對GMM中的每個高斯分量使用不同的比例系數(shù)：

其中：S=∑Mi=1 exp（-βLin）， β為權重系數(shù)。兩個點的局部相似性越大時，對應的ηmn越大，即該高斯分量在GMM中所占的比重越大，其成為對應點的概率就越大。

式（8）和（9）共同體現(xiàn)了全局與局部相似性結合的測度準則。權重系數(shù)β控制全局相似性和局部相似性之間的比例關系。當β非常大時，相當于最小化局部相似性Lmn;當β很小時，趨向于最小化全局相似性Gmn。全局與局部相似性測度提供了一種靈活的方式，通過最小化兩個點集之間的全局相似性和局部相似性上的差異來估算對應關系。

3.3 全局與局部測度引導的配準過程

在非剛性點集配準中，能否找到正確的對應關系對配準效果有很大的影響。當兩個點集形態(tài)近似時，距離較近的點，往往是對應點;而當兩個點集形態(tài)差別較大時，具有相似鄰域結構的點，很有可能是對應點。因此采用全局與局部相似性結合的測度準則引導配準過程。在優(yōu)化算法迭代的初期局部相似性起較大的作用，從而盡快找到正確的對應關系;在迭代后期，兩個點集基本對齊，此時全局相似性測度起的作用變大，從而確保算法收斂到較小的配準誤差。

配準過程首先以浮動點集X和參考點集Y為輸入;其次使用GMM進行建模，把X看成GMM的質(zhì)心，Y看成由GMM產(chǎn)生的數(shù)據(jù)點;再次構造目標函數(shù)并進行參數(shù)設置;然后使用EM優(yōu)化策略對目標函數(shù)中的參數(shù)進行估計，交替地執(zhí)行期望步（Expectation step，E-step）和最大化步（Maximization step， M-step）兩步，直至算法收斂;最后輸出浮動點集的變換公式和兩個點集之間的對應關系。本文算法流程如圖2所示。

3.3.2 參數(shù)設置

優(yōu)化算法最大迭代次數(shù)設為100。鄰域中點的個數(shù)設為K=4。初始時權重系數(shù)β取一個較大的值，一般設β=K2，每一次迭代均乘以一個系數(shù)r=0.95，類似于退火率，即β←rβ。這樣能夠保證優(yōu)化迭代的前期局部相似性起主導作用，后期全局相似性起主導作用，從而在盡快找到正確對應關系的同時得到較小的配準誤差。

3.3.3 E-step階段

在E-step階段，首先計算Y與T（X）之間的局部相似性，更新每個高斯分量的比例系數(shù)ηmn;然后根據(jù)貝葉斯定理，使用當前參數(shù)計算后驗概率分布Pmn。

3.3.4 M-step階段

在M-step階段，計算并更新參數(shù)。對式（10）求導，并使其等于0，可以求得相應參數(shù)。

非剛性變換公式T定義為點的初始位置加上一個位移函數(shù)v，即T（X）=X+v（X）[9]。使用高斯矩陣核定義再生核希爾伯特空間，可以得到v的函數(shù)形式：

其中：wm是系數(shù)矩陣W=[w1 w2 … wM]T中的元素;Γ（x，xm）是核矩陣。正則項φ（T）在這里的具體形式為WΓWT[14]。則變換公式可表示為T=X+ΓW。同時求得參數(shù)ω和σ2的值。

3.3.5 算法輸出

EM算法交替地執(zhí)行E-step和M-step，直至最大迭代次數(shù)或配準誤差小于給定閾值。最終浮動點集的變換公式由T給出，點的對應關系由后驗概率矩陣P給出。

4 實驗與分析

為了評估算法的性能，設計了二維與三維合成點集配準實驗和真實醫(yī)學圖像配準實驗，并對實驗結果進行了比較和定量分析。所有實驗的軟件環(huán)境為Matlab R2013a，硬件環(huán)境為16GB內(nèi)存、i7-4770K Inter處理器的PC。

4.1 合成點集配準實驗

采用二維的Chui-Rangarajan合成數(shù)據(jù)[8]和三維的人臉點集[9]數(shù)據(jù)進行實驗。且實驗結果與薄板樣條魯棒點匹配（Thin Plate Spline Robust Point Matching， TPS-RPM）算法[8]、高斯混合模型點集配準（Gaussian Mixture Models point set REGistration， GMMREG）算法[11]、基于L2E估計的魯棒點匹配算法（Robust Point Matching algorithm based on L2E estimation， RPM-L2E）[14-15]、基于全局局部混合距離與薄板樣條的點集配準算法（Global and Local Mixture Distance and Thin Plate Spline based point set registration algorithm， GLMDTPS）[16]和CPD算法[9]進行比較。性能評價采用與文獻[8，10，12-15]相同的方法。由于合成數(shù)據(jù)中每組測試樣例的真正對應點已經(jīng)預先標識出，即真正的對應關系是已知的，因此使用配準后兩個點集的真正對應點之間的平均歐氏距離，即均方根誤差（Root Mean Squared Error， RMSE）進行定量分析。誤差計算公式為：

4.1.1 二維點集配準實驗

選取Chui-Rangarajan合成數(shù)據(jù)小魚形狀的變形、遮擋、異常點3組數(shù)據(jù)進行實驗。每組包括由低到高5種不同的退化等級，每一等級包含100個測試樣例。圖3為6種算法的配準結果展示。

圖4是配準性能的定量分析與比較，用誤差棒表示每種退化等級樣例的平均誤差和標準差。

圖4中橫坐標表示5種不同的退化等級。

每種算法所對應的誤差棒的中心表示平均誤差，誤差棒的長度表示標準差。從圖4可以看出，當變形、遮擋、異常點等退化程度較小時，各種算法均有不錯的性能，但隨著退化程度的增大，本文算法的優(yōu)勢越來越明顯。表1為6種算法3組實驗的平均誤差統(tǒng)計，可看出本文算法在各種情況下都有一定的優(yōu)勢。

4.1.2 三維點集配準實驗

三維點集配準實驗選取文獻[9]中的人臉點集。變形、遮擋和異常點3組數(shù)據(jù)，每組選取50個測試樣例。本文算法的配準效果如圖5所示。

在二維和三維點集上的實驗結果表明，本文算法與TPS-RPM、GMMREG、RPM-L2E、GLMDTPS和CPD算法相比，平均RMSE分別下降了39.93%、42.45%、32.51%、22.36%和11.76%。

4.2 真實圖像配準實驗

為了比較在CPD配準框架中引入全局與局部相似性測度前、后的配準效果，選用20組醫(yī)學MRI圖像進行實驗。首先提取圖像的輪廓，生成點集，對點集進行配準;然后根據(jù)變換公式生成配準后的圖像。圖6為一組人腦MRI圖像的測試樣例，展示了CPD算法與本文算法的配準結果。圖6第一行為初始的浮動圖像與參考圖像，以及兩種算法配準后浮動圖像與參考圖像的灰度差圖像;第二行為對應的點集。

配準效果的定量分析采用輪廓點集配準誤差和圖像灰度誤差兩種方法。

首先采用類似文獻[17]人工標注點的方法，在兩個輪廓點集上查找具有明顯特征的點對，人工標注50個標準對應點。計算兩個輪廓對應標注點之間的RMSE。兩種算法的平均誤差如表3中人工標注點的平均誤差所示。召回率定義為配準后誤差小于給定閾值的對應點與全部標注的對應點在數(shù)量上的比值。圖7展示了20組測試樣例在各種閾值下CPD與本文算法召回率的變化曲線。實驗結果表明，本文算法優(yōu)于CPD算法。

其中：Y（u，v）與X（u，v）是兩個圖像在坐標點（u，v）上的灰度值;U×V是圖像的分辨率。表3中像素灰度的均方根誤差展示了20組測試樣例使用兩種算法配準后的圖像灰度誤差，同樣可以看出本文算法優(yōu)于CPD算法。

5 結語

在非剛性點集配準算法中，為了充分利用點的鄰域結構信息來查找對應關系，將改進的全局和局部混合距離融入到CPD算法的高斯混合模型配準框架中，提出了基于全局與局部相似性的測度準則。使用權重系數(shù)調(diào)節(jié)全局相似性測度和局部相似性測度二者的比例。采用EM算法迭代地求解對應關系和變換公式，迭代過程中局部相似性所占比重由大變小，全局距離所占比重由小變大，從而在保證快速找到正確對應關系的前提下，提高配準的精確度。通過合成數(shù)據(jù)與真實數(shù)據(jù)的實驗結果比對和分析，可以看出本文算法具有較好的配準性能。然而將點集鄰域結構的計算引入到相似性測度中，加大了算法的計算量、增加了算法的運行時間，可能無法滿足實時性要求較高的應用場合。因此未來的工作在進一步研究如何更好地利用點集結構信息的同時，還需要考慮如何提高算法的執(zhí)行速度。

參考文獻（References）

[1] MA J， ZHAO J， JIANG J， et al. Locality preserving matching[J]. International Journal of Computer Vision， 2019， 127（5）： 512-531.

[2] FERRANTE E， PARAGIOS N. Slice-to-volume medical image registration： a survey[J]. Medical Image Analysis， 2017， 39： 101-123.

[3] GUAN S， WANG T， MENG C， et al. A review of point feature based medical image registration[J]. Chinese Journal of Mechanical Engineering， 2018， 31： No.76.

[4] TANG H， PAN A， YANG Y， et al. Retinal image registration based on robust non-rigid point matching method[J]. Journal of Medical Imaging and Health Informatics， 2018， 8（2）： 240-249.

[5] WU Y， MA W， ZHANG J， et al. Point-matching algorithm based on local neighborhood information for remote sensing image registration[J]. Journal of Applied Remote Sensing， 2018， 12（1）： No.016002.

[6] BELONGIE S， MALIK J， PUZICHA J. Shape matching and object recognition using shape contexts[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2002， 24（4）： 509-522.

[7] BESL P J， McKAY N D. Method for registration of 3-D shapes[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 1992， 14（2）： 239-256.

[8] CHUI H， RANGARAJAN A. A new point matching algorithm for non-rigid registration[J]. Computer Vision and Image Understanding， 2003， 89（2/3）： 114-141.

[9] MYRONENKO A， SONG X. Point set registration： coherent point drift[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2010， 32（12）： 2262-2275.

[10] WANG G， CHEN Y. Fuzzy correspondences guided Gaussian mixture model for point set registration[J]. Knowledge-Based Systems， 2017， 136： 200-209.

[11] JIAN B， VEMURI B C. Robust point set registration using Gaussian mixture models[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2011， 33（8）： 1633-1645.

[12] WANG G， ZHOU Q， CHEN Y. Robust non-rigid point set registration using spatially constrained Gaussian fields[J]. IEEE Transactions on Image Processing， 2017， 26（4）： 1759-1769.

[13] WANG G， CHEN Y， ZHENG X. Gaussian field consensus： a robust nonparametric matching method for outlier rejection[J]. Pattern Recognition， 2018， 74： 305-316.

[14] MA J， QIU W， ZHAO J， et al. Robust L2E estimation of transformation for non-rigid registration[J]. IEEE Transactions on Signal Processing， 2015， 63（5）： 1115-1129.

[15] MA J， ZHAO J， TIAN J， et al. Robust estimation of nonrigid transformation for point set registration[C]// Proceedings of the IEEE 2013 Conference on Computer Vision and Pattern Recognition. Piscataway： IEEE， 2013： 2147-2154.

[16] YANG Y， ONG S H， FOONG K W C. A robust global and local mixture distance based non-rigid point set registration[J]. Pattern Recognition， 2015， 48（1）： 156-173.

[17] 王麗芳， 王雁麗， 藺素珍， 等. 基于改進的Zernike矩的局部描述符與圖割離散優(yōu)化的非剛性多模態(tài)腦部圖像配準[J]. 計算機應用， 2019， 39（2）： 582-588. （WANG L F， WANG Y L， LIN S Z， et al. Non-rigid multi-modal brain image registration based on improved Zernike moments based local descriptor and graph cuts discrete optimization[J]. Journal of Computer Applications， 2019， 39（2）： 582-588.）