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

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

1 引 言

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

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

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

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

2 基本原理

2.1 模型的粗配準

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

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

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

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

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

2.2 模型的精配準

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

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

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

(1)

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

(2)

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

(3)

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

(4)

當平均距離dk+1小于一定的閾值τ時，迭代終止，否則重復迭代直到dk+1<τ或迭代次數達到預設的閾值為止。

2.2.2 改進的ICP算法

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

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

圖1 雙向距離搜索示意圖

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

(5)

(6)

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

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

(7)

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

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

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

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

(8)

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

本文以pi值作為權重，用一個新的函數表示加權最小二乘法：

(9)

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

(10)

使用SVD算法求解對應點集，具體步驟如下：

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

(11)

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

(12)

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

Step4 由U和V，計算旋轉矩陣R和平移矢量t。

(13)

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

3 實驗與結果分析

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

3.1 斯坦福點云數據配準實驗

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

圖4 斯坦福bunny數據配準結果Fig.4 Stanford bunny data registration results

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

圖5 bunny數據配準結果分析Fig.5 Bunny data registration result analysis

表1 bunny數據配準的定量分析結果

3.2 心臟點云數據配準實驗

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

圖6 兩個圖譜與目標圖像之間的配準結果Fig.6 Registration result between two maps and the target image

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

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

圖7 心臟數據配準結果分析

Fig.7 Analysis of cardiac data registration results

表2 心臟數據配準的定量分析結果

Tab.2 Quantitative analysis of cardiac data registration

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

4 結 論

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