王 賓，劉 林，侯榆青，賀小偉

(西北大學(xué) 信息科學(xué)與技術(shù)學(xué)院，西安市影像組學(xué)與智能感知重點實驗室 , 陜西 西安 710127)

1 引 言

三維點云配準(zhǔn)技術(shù)就是求解不同點云之間的旋轉(zhuǎn)平移變換，將源點云調(diào)整到和目標(biāo)點云相同的坐標(biāo)位置，實現(xiàn)點云表達(dá)的實物信息的融合,在三維輪廓測量[1],醫(yī)療診斷和圖像處理[2]，虛擬現(xiàn)實[3]和逆向工程[4]等領(lǐng)域有著重要作用。

在目前流行的醫(yī)學(xué)圖像分割方法中，大多數(shù)研究者將圖譜配準(zhǔn)或者關(guān)鍵特征點配準(zhǔn)與機器學(xué)習(xí)的方法結(jié)合起來，首先通過圖譜配準(zhǔn)或關(guān)鍵點配準(zhǔn)得到圖像分割的初始位置或者統(tǒng)計概率圖，然后通過各種機器學(xué)習(xí)的方法對分割結(jié)果進一步精確微調(diào)。大多數(shù)情況下單個圖譜因不能體現(xiàn)個體間的差異而難以取得滿意的結(jié)果，所以許多研究者通過整合多個圖譜單獨配準(zhǔn)的結(jié)果來彌補這些差異。在基于多圖譜配準(zhǔn)的醫(yī)學(xué)圖像分割方法中[5-6]，由于生物體自身的差異以及器官的形變、蠕動等原因，導(dǎo)致不同樣本中組織器官的形狀、位置等有一些偏差，為了確立各圖譜與目標(biāo)圖像之間的對應(yīng)點關(guān)系，需要選擇目標(biāo)圖像作為參考來配準(zhǔn)所有圖譜。由于心臟是所有器官中形狀最為復(fù)雜的，因此三維點云數(shù)據(jù)模型復(fù)雜，不同圖譜和目標(biāo)圖像之間具有一定差異；同時，在進行CT心臟數(shù)據(jù)分割時，由于輪廓不清晰以及人為分割誤差等因素可能造成數(shù)據(jù)部分殘缺，這幾種情況均加大了配準(zhǔn)的難度。

為了提高三維心臟點云數(shù)據(jù)配準(zhǔn)的精度和速度，本文采用先粗配再細(xì)配的配準(zhǔn)方法。目前常用的粗配準(zhǔn)方法主要包括重心重合法[7],標(biāo)簽法[8]和特征提取法[9]。重心重合法通過計算兩個模型的重心坐標(biāo)將其重心重合，這種方法的優(yōu)點是可以縮小平移差異，缺點是無法縮小旋轉(zhuǎn)差異；標(biāo)簽法通過檢測時人為貼上的特征點進行定位，從而將兩個模型的位置大致調(diào)整到同一個坐標(biāo)系下，這種方法的缺點是過于依賴新加入的特征點，人為因素對于結(jié)果的影響比較大；特征提取法一般根據(jù)提取到的平面特征、輪廓特征和邊角特征等將兩個模型表面特征相同部分重合使其初始位置差異進行縮小，這種方法的缺點是要求模型具有足夠多的明顯特征。迭代最近點算法(Iterative Closet Point, ICP)是Besl和Mckay[10]在1992年提出的一種精確配準(zhǔn)算法。由于其思想簡單，容易理解，方便操作和良好的配準(zhǔn)效果而成為了主流的配準(zhǔn)算法，在各種領(lǐng)域的配準(zhǔn)問題上得到了廣泛地應(yīng)用[11-12]。但是，傳統(tǒng)的ICP算法對于初始值的敏感性比較高，并且在每次迭代中都對點云數(shù)據(jù)進行全局搜索以查找對應(yīng)點，所以計算量大，而且，對于較為復(fù)雜的兩個模型不能得到良好的配準(zhǔn)結(jié)果。針對傳統(tǒng)ICP算法的缺點，中外許多學(xué)者提出了各種改進算法以提高算法的速度和精度。Rusinkiewicz等[13]提出基于法線空間的均勻采樣方法，在具有較少特征點的點云進行配準(zhǔn)時，該方法可以提高配準(zhǔn)速度。Chen等[14]采用待配準(zhǔn)點云的點法線與參考點點云的交點確定對應(yīng)點，該方法能減少迭代次數(shù)并加快算法收斂速度，但在某些情況下魯棒性較差。王欣等[15]提出基于點云邊界特征點的改進ICP算法，提高了逆向工程中點云數(shù)據(jù)配準(zhǔn)的效率和精度。Chetverikov等[16]提出的TrICP算法基于穩(wěn)健統(tǒng)計分析中的最小截斷二乘法的一致性對點云進行配準(zhǔn)，適用于點云初始值比較差的配準(zhǔn)，該方法可以提高配準(zhǔn)精度，但是仍會有異常數(shù)據(jù)影響結(jié)果準(zhǔn)確性，并且該方法在效率上改善不大。

為了解決多圖譜配準(zhǔn)中心臟點云數(shù)據(jù)在具有復(fù)雜形狀和局部殘缺情況下的配準(zhǔn)問題，考慮到主成分分析法(Principal Component Analysis, PCA)在特征降維和圖像去噪的廣泛應(yīng)用[17-18]，本文采用主成分分析法進行粗配準(zhǔn)[19-20]，用KD-tree進行最近點搜索提高搜索速度，在精配準(zhǔn)方面提出了基于雙向距離比例的ICP算法提高配準(zhǔn)的精度，最后使用加權(quán)最小二乘法迭代求解坐標(biāo)轉(zhuǎn)換[21]。

2 基本原理

2.1 模型的粗配準(zhǔn)

在精確配準(zhǔn)之前，由于兩個模型的點集處在同一坐標(biāo)系的不同位置，沒有足夠多的重復(fù)率，初始位置差異較大，因此點集迭代初始值不能滿足精確配準(zhǔn)的要求。在點云數(shù)目較大的情況下，進行最近點搜索時，算法的復(fù)雜度和計算量會劇增，為了提高算法執(zhí)行效率和配準(zhǔn)精度，在進行精確配準(zhǔn)之前，先對兩個點集進行全局的粗配準(zhǔn)，以便為精確配準(zhǔn)提供良好的初始值。

綜合考慮各粗配準(zhǔn)方法的優(yōu)缺點，本文使用的粗配準(zhǔn)方法是主成分分析法(PCA)，算法實現(xiàn)步驟為：

(1)設(shè)n維的三維點數(shù)據(jù)集P={P1,P2,P3,...,Pn}，計算其均值和協(xié)方差矩陣cov；

(2)對矩陣cov進行特征向量分解，得到的正交特征向量作為點集P的三坐標(biāo)軸x,y,z，以其均值作為坐標(biāo)系原點O；

(3)按照上述方法分別對兩組點云數(shù)據(jù)集計算坐標(biāo)系，并求取兩組點云數(shù)據(jù)集的參考坐標(biāo)系間的剛體轉(zhuǎn)換參數(shù)，以實現(xiàn)點云數(shù)據(jù)集間的粗配準(zhǔn)。

2.2 模型的精配準(zhǔn)

精確配準(zhǔn)即局部配準(zhǔn)，在粗配準(zhǔn)的基礎(chǔ)上，如何精確求解兩點云數(shù)據(jù)集之間剛體轉(zhuǎn)換參數(shù)是整個配準(zhǔn)過程中最為重要和核心的部分。本文提出基于雙向距離比例的ICP算法提高配準(zhǔn)的精度，利用迭代算法求解基于加權(quán)最小二乘理論的最優(yōu)解，得到具有較高精度的全局最優(yōu)化剛體轉(zhuǎn)換參數(shù)，完成兩個點云數(shù)據(jù)集的配準(zhǔn)。針對點云配準(zhǔn)的本質(zhì)問題，即求解兩個點集之間的旋轉(zhuǎn)矩陣和平移矢量，本文采用奇異值分解法(Singular Value Decomposition, SVD)進行求解。

2.2.1 傳統(tǒng)的ICP算法

假設(shè)源數(shù)據(jù)點集P：{Pi,i=1,2,...,Np}及目標(biāo)數(shù)據(jù)點集Q：{Qi,i=1,2,...,NQ}，每一次迭代中，執(zhí)行如下步驟：

(1)

(2)通過最小二乘法計算最優(yōu)剛性變換：

(2)

(3)對Pk+1做旋轉(zhuǎn)平移變換：

(3)

(4)計算平均距離dk+1：

(4)

當(dāng)平均距離dk+1小于一定的閾值τ時，迭代終止，否則重復(fù)迭代直到dk+1<τ或迭代次數(shù)達(dá)到預(yù)設(shè)的閾值為止。

2.2.2 改進的ICP算法

針對目前常用的ICP算法存在的問題本文提出了基于雙向距離比例改進的ICP算法，使其更加方便的應(yīng)用在基于多圖譜配準(zhǔn)的醫(yī)學(xué)圖像分割問題上，同時保證配準(zhǔn)精度和速度。

在經(jīng)典的ICP算法中，計算過程中的每一組匹配點對具有相同的權(quán)重，對于錯誤匹配點沒有進行任何處理，這樣會使得一些錯誤匹配點對影響配準(zhǔn)結(jié)果的精確度。目前常用的TrICP方法是基于最小截斷二乘法的思想，即將所有匹配的點對之間的距離計算并排序，把殘差較大的直接去除，該方法可以減小離群點對配準(zhǔn)結(jié)果的影響，但是去除匹配點對的比重不好掌握，如果去除過多，有可能去掉正確匹配點對，如果去除較少，有可能仍舊會殘留錯誤匹配點對，匹配精度無法保證。因此，受加權(quán)最小二乘法思想的啟發(fā)，給每組匹配點對一個權(quán)值，本文提出基于雙向距離比例的ICP算法來提高配準(zhǔn)的精度，既能減少誤匹配對配準(zhǔn)精度的影響，也能保證搜索到更多的點對參與最終配準(zhǔn)結(jié)果的求解。算法的主要思想是:對于目標(biāo)數(shù)據(jù)中的點Qi在源數(shù)據(jù)點集中尋找最近的點Pi，如圖1(a)所示，然后對點Pi反向在目標(biāo)數(shù)據(jù)中尋找最近的點Ql，搜索到的點Ql有可能是點Qi本身，也有可能是其相鄰點，如圖1(b)所示。

圖1 雙向距離搜索示意圖

分別建立雙向匹配，如式(5)和式(6)所示：

(5)

(6)

其中，Rk-1，tk-1分別為每次迭代過程中的旋轉(zhuǎn)矩陣和平移矢量。

采用式(7)基于雙向距離求比值。

(7)

圖2 本文所用指數(shù)函數(shù)特性示意圖Fig.2 Exponential function characteristic diagram used in this paper

在此算法中，基于加權(quán)最小二乘求解最優(yōu)變換，權(quán)重的賦值即概率值是需要解決的重要問題。根據(jù)指數(shù)函數(shù)的特性，引入一個指數(shù)函數(shù)來計算概率值pi，有效反映每一點對雙向距離比值與應(yīng)該賦予權(quán)重之間的關(guān)系，本文所采用指數(shù)函數(shù)特性如圖2所示。其中，pi是與ρi相關(guān)的一個概率值，此概率值可以反應(yīng)點對是否屬于內(nèi)點的概率，即正確匹配的概率。

概率值的計算如式(8)所示：

pi,k=e-λ(ρi,k-1),

(8)

其中：λ是預(yù)設(shè)參數(shù)，依據(jù)指數(shù)函數(shù)的特性,可以反映出隨著雙向距離的比值ρi增大，概率值pi下降的特性。如果雙向距離相等，那么最大概率的被視為內(nèi)點，如果ρi比值增加，概率值會下降。

本文以pi值作為權(quán)重，用一個新的函數(shù)表示加權(quán)最小二乘法：

(9)

通過加權(quán)最小二乘計算最優(yōu)變換：

(10)

使用SVD算法求解對應(yīng)點集，具體步驟如下：

Step2 以重心為三維坐標(biāo)系的原點重建源數(shù)據(jù)點集和目標(biāo)數(shù)據(jù)點集坐標(biāo)系，并以Pi′和Qi′作為下一步的運算數(shù)據(jù)，以簡化運算。

(11)

Step3 以Pi′和Qi′為基礎(chǔ)構(gòu)造矩陣H，利用奇異值分解法分解矩陣H，得到矩陣U和V。

(12)

圖3 改進的ICP算法流程圖Fig.3 Flow chart of improved ICP algorithm

Step4 由U和V，計算旋轉(zhuǎn)矩陣R和平移矢量t。

(13)

改進的ICP算法總體描述為：首先輸入源數(shù)據(jù)和目標(biāo)數(shù)據(jù)，對其進行粗配準(zhǔn)；然后使用KD-tree進行雙向距離搜索，并且對每一點對計算雙向距離比值和概率值；最后利用加權(quán)最小二乘計算最優(yōu)變換，運用奇異值分解法得到最終的變換矩陣。判斷是否收斂，如果不收斂，重復(fù)迭代，如果收斂，源數(shù)據(jù)和目標(biāo)數(shù)據(jù)配準(zhǔn)并輸出最終結(jié)果，算法流程圖如圖3所示。

3 實驗與結(jié)果分析

在CT心臟點云數(shù)據(jù)的配準(zhǔn)中，主要面臨具有較大初始位置差異、較大形狀差異和局部殘缺的心臟點云數(shù)據(jù)配準(zhǔn)的精度以及速度問題，因此，本文設(shè)計兩組實驗驗證本文算法在不同情況下的點云數(shù)據(jù)配準(zhǔn)精度以及速度。一組是斯坦福點云數(shù)據(jù)配準(zhǔn)實驗，一組是心臟點云數(shù)據(jù)配準(zhǔn)實驗。實驗中，設(shè)置式(8)中的λ=6，最大迭代次數(shù)設(shè)為150。本部分實驗均在配置為3.00 GHz Intel Core i5 CPU 和8.00 GB RAM的個人計算機上完成，計算工具為MATLAB R2015b。

3.1 斯坦福點云數(shù)據(jù)配準(zhǔn)實驗

為了驗證本文算法對初始值敏感性和初始姿態(tài)重疊度具有較好魯棒性，首先采用公用的斯坦福點云數(shù)據(jù)設(shè)置一組實驗，源數(shù)據(jù)為斯坦福bunny點云數(shù)據(jù)，目標(biāo)數(shù)據(jù)為源數(shù)據(jù)在三維空間做(Re,te)的較大幅度旋轉(zhuǎn)平移變換，并分別在整個點云數(shù)據(jù)上隨機剔除5%，10%，15%，20%的點云。實驗引入?yún)?shù)間相對的誤差來衡量配準(zhǔn)精度，εR=‖Rk-Re‖F(xiàn)，εt=‖tk-te‖2。其中(Rk,tk)是最終求解的最優(yōu)變換，(Re,te)是人為對目標(biāo)數(shù)據(jù)做的旋轉(zhuǎn)平移變換，‖·‖F(xiàn)是F范數(shù)。為了研究本文方法在配準(zhǔn)精度和配準(zhǔn)速度方面的性能，將本文算法與經(jīng)典ICP算法和文獻(xiàn)18提出的TrICP算法進行了對比。其中，剔除20%比例點云的配準(zhǔn)結(jié)果如圖4所示。

圖4 斯坦福bunny數(shù)據(jù)配準(zhǔn)結(jié)果Fig.4 Stanford bunny data registration results

通過表1和圖5的數(shù)據(jù)分析可以得出在剔除不同比例點云數(shù)目的情況下，TrICP算法通過剔除異常點，較經(jīng)典ICP算法減小了配準(zhǔn)誤差，但精度仍不及本文算法，本文算法的相對旋轉(zhuǎn)誤差較經(jīng)典ICP算法減少了48%～49%，較TrICP算法減少了30%～32%，相對平移誤差較經(jīng)典ICP算法減少了47%～52%，較TrICP算法減少了36%～45%。同時，TrICP算法的效率稍優(yōu)于經(jīng)典ICP算法，而本文算法采用KD-tree進行雙向距離搜索，同時由于配準(zhǔn)精度的提高，算法收斂速度快，因此效率大幅優(yōu)于以上兩種算法，并且以上兩種算法在殘缺程度(剔除點云比例)較大的情況下，相對誤差會增大，本文算法在殘缺程度較大的情況下依舊可以得到精確的配準(zhǔn)結(jié)果，同時在減少點云數(shù)目的情況下旋轉(zhuǎn)誤差和平移誤差基本不受影響，趨于穩(wěn)定，而且由于需要搜索點云數(shù)目的減少，配準(zhǔn)時間減少，4 500～6 000的點云數(shù)目根據(jù)模型殘缺程度以及模型復(fù)雜程度的不同一般只需要0.6～1.3 s的配準(zhǔn)時間。因此使用本文方法可以通過減少目標(biāo)數(shù)據(jù)點云數(shù)目來提高搜索效率，降低時間。通過對斯坦福點云數(shù)據(jù)配準(zhǔn)結(jié)果的分析可以初步確定本文改進的ICP算法可以精確高效穩(wěn)定的解決心臟點云數(shù)據(jù)的配準(zhǔn)問題。

圖5 bunny數(shù)據(jù)配準(zhǔn)結(jié)果分析Fig.5 Bunny data registration result analysis

表1 bunny數(shù)據(jù)配準(zhǔn)的定量分析結(jié)果

3.2 心臟點云數(shù)據(jù)配準(zhǔn)實驗

通過觀察圖6的配準(zhǔn)結(jié)果可以發(fā)現(xiàn)在具有較大形狀差異和殘缺差異的心臟點云數(shù)據(jù)配準(zhǔn)中，本文改進的ICP算法的配準(zhǔn)結(jié)果吻合度更高，特別是圖譜的邊緣部分以及局部形狀差異大的部分。通過表2和圖7的配準(zhǔn)定量分析結(jié)果和對比結(jié)果也可得到雖然由于模型差異較大導(dǎo)致配準(zhǔn)誤差較大，但是本文算法的配準(zhǔn)精度仍舊高于經(jīng)典ICP算法和TrICP算法，配準(zhǔn)平均誤差較經(jīng)典ICP算法減少了約21%，較TrICP算法減少了約13%。通過對比第一組和第二組的誤差可以得到本文方法在第二組心臟配準(zhǔn)的誤差改善更大，通過觀察待配準(zhǔn)圖譜可以發(fā)現(xiàn)，第二個待配準(zhǔn)圖譜較目標(biāo)圖像之間的殘缺更大，特別是心臟上部較目標(biāo)圖像有著明顯的形狀差異和殘缺差異，而第一個待配準(zhǔn)圖譜較目標(biāo)圖像之間差異較小，由此可以體現(xiàn)出本文方法在具有較大殘缺差異和形狀

圖6 兩個圖譜與目標(biāo)圖像之間的配準(zhǔn)結(jié)果Fig.6 Registration result between two maps and the target image

差異時對于精度的提高較經(jīng)典ICP和TrICP的優(yōu)勢更加明顯。在配準(zhǔn)速度方面，TrICP算法稍優(yōu)于經(jīng)典ICP算法，但仍不及本文算法的配準(zhǔn)速度，本文算法由于采用KD-tree進行雙向距離搜索，同時由于配準(zhǔn)精度的提高，算法收斂速度快，4 500左右的點云數(shù)目根據(jù)模型復(fù)雜程度以及殘缺程度只需要1.6～1.8 s左右的配準(zhǔn)時間，因此速度大幅優(yōu)于以上兩種算法。

因此，本文改進的ICP算法可以滿足基于多圖譜配準(zhǔn)的CT心臟分割中具有較大初始位置差異、較大形狀差異和局部殘缺的心臟點云數(shù)據(jù)之間精確高效地配準(zhǔn)。

圖7 心臟數(shù)據(jù)配準(zhǔn)結(jié)果分析

Fig.7 Analysis of cardiac data registration results

表2 心臟數(shù)據(jù)配準(zhǔn)的定量分析結(jié)果

Tab.2 Quantitative analysis of cardiac data registration

方法誤差時間第一組經(jīng)典ICP11.427 018.926 2TrICP10.582 615.524 1本文ICP9.495 71.771 3第二組經(jīng)典ICP14.722 419.173 3TrICP12.985 715.957 6本文ICP10.946 41.622 9

4 結(jié) 論

為了提高具有較大初始位置差異、較大形狀和局部殘缺差異的心臟點云數(shù)據(jù)配準(zhǔn)的精度以及速度問題，本文采用一種先粗配再精配的配準(zhǔn)方法。首先采用主成分分析法進行粗配準(zhǔn)，然后在精配準(zhǔn)方面提出了基于雙向距離比例的ICP算法提高配準(zhǔn)的精度。實驗結(jié)果表明，本文方法具有良好的穩(wěn)定性，較高的精度和效率，配準(zhǔn)平均誤差較經(jīng)典ICP算法減少了約21%，較TrICP算法減少了約13%。改進算法的優(yōu)勢主要體現(xiàn)在兩方面，一是在模型殘缺程度較大的情況下配準(zhǔn)精度基本不受影響，但由于點云數(shù)目減少，配準(zhǔn)時間減少，因此可以通過減少目標(biāo)數(shù)據(jù)點云數(shù)目來提高搜索效率，降低時間；二是在具有較大殘缺差異和形狀差異時本文算法具有很強的魯棒性，在配準(zhǔn)精度方面的優(yōu)勢更加明顯。因此本文算法可以應(yīng)用于基于多圖譜配準(zhǔn)的醫(yī)學(xué)圖像分割中。進一步優(yōu)化配準(zhǔn)方法，以及將本文配準(zhǔn)方法適用于其它應(yīng)用環(huán)境中，比如生物醫(yī)學(xué)的中的不同模態(tài)，不同結(jié)構(gòu)的對象配準(zhǔn)問題，是下一步研究的方向。