基于曲率圖的三維點云數據配準

葛寶臻 ，彭 博 ，田慶國

(1. 天津大學精密儀器與光電子工程學院，天津 300072；2. 光電信息技術教育部重點實驗室，天津 300072)

三維點云數據的配準與融合是計算機視覺領域一個基本而重要的問題．其研究成果在醫學、逆向工程、虛擬現實、計算機視覺等領域得到了廣泛的應用[1]．根據不同視圖點云數據間是否存在畸變，現有的配準算法大致可以分為基于剛性變換的配準和非剛性配準兩類．本文只討論基于剛性變換的配準算法．給定具有重疊部分的三維點云數據，其任務是求取點云數據間的最佳剛性變換，使不同視圖的點云數據準確地融合在一起[2-4]．

迭代最近點(iterative closest point，ICP)算法[5]是當前最流行的剛性配準算法，該算法以相對位置經過粗略配準的兩片點云數據為基礎，迭代求取兩片點云數據間的對應點及相應的剛性變換，直到對應點之間的距離誤差評價函數最?。撍惴ǜ鶕牲c云數據點之間的距離來確定對應點，所以需要待配準的兩點云數據之間基本已經配準．實際上待配準的點云數據間的相對位置往往未知，怎樣利用粗略配準算法將點云數據間的相對位置變換到ICP 算法的收斂域成為問題的關鍵．

目前普遍的解決方法是通過不依賴于三維旋轉和平移變換的點云特征信息來建立點云數據間的對應關系．這些特征信息包括法向[6]、曲率[7]、體積分[8]等簡單的一維特征描述函數，還包括旋轉圖像[9]、形狀內容[10]等高維特征描述函數，有的甚至還包括顏色、強度值等輔助信息．但存在的問題是，低維特征描述函數計算簡單，但所包含的特征信息較少，且對噪聲敏感；高維特征函數能很好地反映特征信息，但計算復雜，時間、空間開銷較大，且怎樣合理地進行特征比較是一個難題．

筆者致力于通過新穎的特征信息和匹配算法快速、準確地得到點云數據間的對應關系．受Gatzke等[11]的啟發，在以點幾何表示的三維點云數據中引入曲率圖的概念，并將其作為三維點云數據的特征描述函數．該描述函數的優點在于通過選擇不同的環數及各環半徑值可得到不同的曲率圖，很容易實現多比例空間的特征保持分析．有了曲率圖的概念以后，本文的主要算法思想是通過一維曲率圖粗選特征點和潛在對應點，然后對特征點和潛在對應點進行高維曲率圖計算，并再次通過高維曲率圖進一步精簡潛在對應點．這樣既保證了豐富的特征信息以獲得正確的對應點，也大大減小了由于對整片點云進行高維曲率圖計算所花費的時間和空間開銷．在特征點的匹配算法上，對已有的Mitra 貪婪算法[12]進行了改進，大大提高了其穩定性．最后通過對真實三維點云數據進行配準，驗證了該算法的有效性．

1 三維點云數據的曲率及曲率圖

在微分幾何中，對于三維歐幾里得空間中可微曲面上的任一點A，給定了兩個主曲率 1k 、2k[13]來衡量該點在不同方向上的彎曲程度．三維點云數據是對真實物體表面進行掃描采樣獲得的離散數據，點云數據中的每個頂點同樣具有主曲率，不但如此，它還是點云數據的一項重要特性．經常提及和主曲率相關的幾個曲率概念還有高斯曲率 K = k1k2、均值曲率

對于三維點云數據，估算其頂點主曲率 1k 、2k 的方法有很多，具代表性的有多項式擬合(polynomial fitting)法[15]、主成分分析(principal component analysis)法[16]以及移動最小二乘(moving least-squares)法[17]等．其中多項式擬合法對噪聲敏感程度最小，穩定性高，但計算相對復雜；主成分分析法求解過程相對簡單，但只適用于采樣點稠密的大型點云；移動最小二乘法具有較好的法向估算精度，但需要求解非線性方程．本文根據不同的點云數據采取與之相適合的主曲率估算方法．但僅利用單一點的均方根曲率C 并不能很好地反映該點周圍點的特征，并且對噪聲比較敏感．曲率圖利用該點及其周圍點的曲率信息分布，能夠準確地反映特征信息的同時有效地抑制噪聲．下面分別介紹以點幾何表示的三維點云數據中一維曲率圖和高維曲率圖的構造方式．這里以降采樣后的點云數據Stanford bunny00[18]為例，對其表面上的任意一頂點A，一維曲率圖的構造比較簡單，如圖1(a)所示，以目標點A 為球心，對半徑為R 的球體內所有頂點的均方根曲率 Ci進行平均，得到的平均值即為該點的一維曲率圖．

圖1 曲率圖的概念示意Fig.1 Sketch of curvature map

式中n 為距離A 點的距離小于R 的頂點個數．對于高維曲率圖，如圖1(b)所示，同樣以目標點A 為球心，利用半徑為 rj(0 ＜ j ≤m)的同心球將鄰域空間劃分為m 個球殼層，然后對每個球殼內的所有頂點的均方根曲率進行平均得到 f(rj)(0＜ j≤m)，高維曲率圖 kmap即為 f(rj)隨 rj的變化關系．用向量表示為

2 基于一維曲率圖提取特征點

有了曲率圖的概念以后，對三維點云數據的每個頂點進行一維曲率圖運算．根據所有頂點的一維曲率圖分布情況進行特征點的提取．這里采用的思想是對所有頂點的一維曲率圖進行統計分析，求取一維曲率圖的平均值μ和方差σ，選取一維曲率圖 f(R)處于μ± kσ之外的點作為特征點．一片點云的數據量很大，使得一維曲率圖往往近似服從正太分布，實驗中發現 k = 1時，獲得的特征點數量大致為點云總數的10%，通過調整系數k 可選擇不同程度的特征點．為了使選取的特征點能更好地反映點云的特征信息，抑制噪聲，對每個頂點選取不同的鄰域范圍做多比例空間的特征分析．這里通過改變一維曲率圖的構造半徑R ，例如對于半徑為 Ri和 Ri+1時的鄰域范圍分別做特征分析得到特征點集 Pi和 Pi+1，選取 Pi和Pi+1的交集作為特征點．最后的特征點集為

式中n 為做多比例空間特征分析選取的半徑個數．對Stanford bunny00 進行了基于一維曲率圖的特征點提取及特征保持分析．算法流程為：采用空間包圍盒法將三維點云數據結構化；在此基礎上按歐式距離進行點云的K 鄰域搜索；采用主成分分析法估算主曲率及均方根曲率(K＝20)；對每個半徑R 分別計算每個頂點的一維曲率圖 f(R)(R ＝0.000 5、0.001、0.002、0.005,mm)；對每個半徑R 分別進行曲率圖統計分析，選取 f(R)處于μ± kσ之外的點作為該曲率圖半徑下的特征點集(k ＝1)；對上述4 個曲率圖半徑R進行多比例空間的特征保持分析，得到特征保持點集；對特征保持點集按一定間隔距離值(0.005,mm)進行采樣．

圖2所示為該算法的特征點提取過程．圖2(a)為R ＝0.000 5,mm 時提取到的特征點，可以看到，通過一維曲率圖提取到的特征點基本上都是處在兔子的耳朵、脖子、腿等輪廓特征比較明顯的部位，驗證了一維曲率圖提取特征點的有效性；圖2(b)顯示了通過多比例空間的特征保持分析后得到的特征保持點集，對比圖2(a)可知，特征保持分析能有效去除一些非特征點或特征不夠顯著的點；圖 2(c)為按0.005,mm 間隔對特征保持點集采樣的結果，對特征保持點集按一定距離采樣是為了進一步減少特征點數量，且讓特征點間有一定的距離，為后面的配準算法打下了基礎．

圖2 特征點的提取過程Fig.2 Process of extracting feature points

3 高維曲率圖精簡潛在對應點

對于要配準的兩片點云數據，通常稱之為源點集和目標點集．理論上，只要在源點集和目標點集之間至少找出3 對對應點，就能夠計算出這兩片點云數據間的剛性變換關系，從而將這兩片點云數據配準到一起．因此，準確找到點集間的對應點成為配準的關鍵．

本文的思想是對源點集進行基于一維曲率圖的特征點提取，然后對目標點集進行相應的一維曲率圖運算，一維曲率圖半徑為R ．對源點集中的每個特征點根據半徑為R 時一維曲率圖的相似性在目標點集中求取相應的對應點，即選取 fp(R) ? fq(R )≤Δ的點作為該特征點的潛在對應點，Δ為設定閾值，設定原則是保證理論上具有對應點的特征點都能在潛在對應點中找到正確的對應點．通過實驗發現，該值設為所有源點集中特征點的一維曲率圖平均值的10%時，基本能夠滿足上述要求．此時源點集中的一點 pi往往對應目標點集中的多個點 qi1，qi2，…，qik．為進一步減少潛在對應點的數量，有效去除錯誤對應點，對pi及其潛在對應點 qi1，qi2，…，qik進行高維曲率圖運算，高維曲率圖的半徑大小及環數根據點云采樣密度決定．然后將 pi的高維曲率圖分別與 qi1，qi2，…，qik的高維曲率圖做比較，比較的方式直接采用高維空間的歐式距離，即E2中選擇與 ie 沒有交集的2 點對 je ，使得 ie 、je 組成的4 點對的距離均方差dRMS 最小．每構成一個新的4 點對，便從E2中刪除包含這4 點對中任意一點的2 點對，遍歷直到E2為空．最后得到4 點對集E4，同樣按照距離均方差的值對E4升序排列．這里距離均方差dRMS 的定義為

(3) 粗略匹配對應點．同理可以對E4進行組合，得到8 點對集E8，乃至16 點對集E16等．在對應點對不是太多的情況下一般得到E4、E8即可．選擇E4或E8中距離均方差最小的4 點對或8 點對作為粗略的對應點．

(4) 精確匹配對應點．在目標點集中以粗略對應點為球心，半徑為 CR 做球，處于球內的點作為潛在對應點，遍歷所有的對應點組合，選擇使dRMS 最小的點對組合作為最終的對應點．CR 的選取值越大得到的結果更精確，但耗時也越大．

對Stanford bunny00 和Stanford bunny45 進行上述匹配算法，點對組合至E4，CR 取0.005,mm，得到4對對應點，如圖3 所示．從圖中可以看到，兩視圖間的對應點對基本上處于兔子的同一個位置，證實了該匹配算法的有效性．

式中：m 為高維曲率圖的維數；fp(rj)、fq(ri)分別為點 pi、qi的高維曲率圖分量．對 pi的所有潛在對應點按d 的大小進行升序排列，選取靠前的M 個對應的 qi作為 pi最終的潛在對應點．M 的選取原則是盡量保證正確的對應點包括在其中，而且值還不能過大，不然會給后面的匹配算法帶來困難．經過實驗驗證M 值選取為10，即一個點具有10 個對應點時，點云重疊部分的特征點基本都能找到正確的對應點，而且對后續的匹配算法不會增加太多負擔．

4 精確匹配對應點

目前主要匹配算法思想是利用形狀的空間變換不變性來有效匹配對應點．在眾多的算法中，Mitra的貪婪算法雖然不能保證找到最佳的對應點，但它能高效快速地找到接近最佳匹配的對應點對．為了較為精確地匹配對應點，筆者提出下面的兩步匹配算法：首先，通過貪婪算法粗略匹配對應點，然后同樣利用形狀空間變換不變性，遍歷這些點某個鄰域范圍內的所有點對組合，找出最佳對應點．具體算法有如下4 個步驟．

(1) 構造點對．對于源特征點集中的每個點對組合 pi和 pj，從其潛在對應點中選擇 qik和 qjl，使得

圖3 精確匹配對應點Fig.3 Accurately matching of corresponding points

5 結果與討論

(2) 組合點對．順序從E2中取出2 點對 ei，再從

針對Stanford bunny00 和Stanford bunny45，配準前如圖4(a)所示，利用匹配算法得到對應點以后，運用奇異值分解(SVD)法[19]直接求取剛性變換，將Stanford bunny45 粗略配準到Stanford bunny00 上，得到圖4(b)所示結果，可以看到兩片點云數據基本上被配準在一起，但仍出現大片單一視圖的點云數據，需要進一步的精確配準．

為此，將粗配準后的兩片點云數據作為輸入，利用經典ICP 算法進一步精配準，ICP 算法的具體算法流程為：

(1) 通過隨機采樣的方式從源點集中選擇一部分點集P；

中日醫院呼吸中心煙草病學與戒煙中心主任肖丹教授指出，無論吸什么樣的煙，只要吸煙就是有害健康的。大量研究表明，煙草煙霧中有超過7000種化學物質，包括69種致癌物，依然含有尼古丁、丙二醇、丙三醇、重金屬等有害物質，會增加患慢性阻塞性肺疾病、肺癌、心血管疾病的風險，可以引起人體中多個器官的損害。

(2) 對點集P 中的任何一點，在目標點集中找尋與之歐式距離最近的點作為對應點，得到點集P的對應點集Q；

(3) 從對應點集P 和Q 中刪除距離值大于平均距離的點，剩下的點作為最后的對應點；

(4) 通過對應點求取坐標變換矩陣(四元素法或SVD)；

(5) 構造并計算誤差評價函數，其實就是第4 節提到的距離均方差；

(6) 判斷誤差評價函數的大小是否滿足精度要求，如果滿足則算法終止，返回結果；如不滿足，迭代進入步驟(1)，此時的目標點集為原始目標點集經過步驟(4)獲得的變換后的結果．

通過ICP 算法的精配準得到最終的配準結果如圖4(c)所示，可以看到兩片點云數據交錯有致，達到了較好的配準效果．

圖4 兔子點云的配準過程Fig.4 Process of registration for Stanford bunny

為了驗證本文配準算法的通用性和適應性，利用本實驗室的掃描系統，對圖5 所示的人頭模型進行掃描實驗．但這次掃描條件比較特殊，只利用其中的一個攝像機，分別從不同角度進行掃描，得到一系列具有重疊部分的三維點云視圖．

圖5 人頭模型Fig.5 Model of head

圖 6 所示為以任意兩角度拍攝的兩點云視圖．運用本文的配準算法模塊，對這兩點云視圖進行配準，得到如圖7 所示的結果．圖7(a)為配準前的兩點云視圖在同一坐標系中的表示，可以看出兩點云的視圖旋轉了一定的角度，且存在少量的平移；圖7(b)所示為本文算法粗配準后的點云視圖，這時兩點云數據已經相當接近配準狀態；圖7(c)為經過ICP 精配準的結果，可以看到視圖間的點云交錯有致，達到了良好的配準效果，驗證了本文算法的有效性．

為說明該算法的抗噪性能，對Stanford bunny00和Stanford bunny45 加入不同程度的白高斯噪聲．噪聲方差δ分別設為原始點云相鄰頂點間平均距離的0.1、0.3 和1.0 倍．為了和同類算法進行比較，在1.0倍噪聲時，分別利用旋轉圖像算法[6]和本文的算法進行粗配準．圖8(a)、8(b)分別表示在0.1 倍和0.3 倍噪聲下本文算法的配準結果，圖8(c)、8(d)分別表示在1.0 倍噪聲下利用本文算法和旋轉圖像算法的粗配準結果．從圖中可以看到，本文的配準算法具有很好的抗噪性，即使在1.0 倍噪聲下也能達到不錯的配準效果，遠遠優越于旋轉圖像算法．

圖6 不同角度的兩點云視圖Fig.6 Two views of head with different angles

圖7 人頭點云的配準過程Fig.7 Process of registration for the head

圖8 不同高斯噪聲下的點云配準結果Fig.8 Results of registration for Stanford bunny with different noise levels

為了進一步說明該算法在尋找對應點上的時間復雜度，以體積分(integral volume descriptors)算法[4]作對比．針對點云數據Bunny、Chicken、Dragon[15]，以Bunny000、Chicken_view1 和Dragonstandright_0作為源點集，Bunny045、Chicken_view2 和Dragonstandright_24 作為目標點集．在相同硬件條件下，在源點集中選取相同的100 個特征點，然后分別利用曲率圖和體積分算子在對應的目標點集中尋找對應點，每個特征點的對應點數量設定為20，實驗結果如表1所示．由于本文的匹配算法是通過兩步法的思想尋找對應點，只對很少的部分點集進行了高維曲率圖計算，所以耗費的時間也大大少于體積分算子．

表1 曲率圖和體積分算子找到對應點的耗時對比Tab.1 Comparison of the point matching time for some models

6 結 語

本文在點幾何表示的三維點云數據中，構造了一維、高維曲率圖．并利用曲率圖分兩步確立了具有重疊部分的三維點云數據間的對應點．先利用一維曲率圖粗略找對應點，然后對找出的點進行高維曲率圖運算，剔除大量錯誤對應點．在此基礎上，通過改進的貪婪算法實現了對應點之間的精確匹配．找到匹配點對后，通過SVD 求解剛性變換，將兩片點云數據粗略配準到ICP 算法的收斂域，最后借助經典ICP算法實現了點云數據的精確配準．通過各種配準實驗，體現了該算法不但具有很強的抗噪性，在計算時間上也大大優于其他基于特征點的配準算法．

致 謝：

特別感謝卡內基美隆大學機器人研究所的Daniel Huber 教授提供了利用旋轉圖像進行三維點云數據配準的可執行代碼．同時感謝斯坦福大學的計算機圖形學實驗室提供的三維點云數據．

2006.

［9］Johnson A，Hebert M. Using spin images for efficient object recognition in cluttered 3D scenes [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence，1999，21(5)：674-686.

［10］Frome A，Huber D，Kolluri R，et al. Recognizing objects in range data using regional point descriptors[J].Computer Vision-ECCV，2004，3023：224-237.

［11］Gatzke T，Grimm C，Garland M，et al. Curvature maps for local shape comparison[C]//IEEE International Conference on Shape Modeling and Applications. 2005：244-253.

［12］Mitra N J，Gelfand N，Pottmann H，et al. Guibas registration of point cloud data from a geometric optimization perspective [C]//Proc Symposium on Geometry Processing. 2004：23-32.

［13］Gray A. Modern Differential Geometry of Curves and Surfaces[M]. Boca Raton，FL，USA：CRC Press，1993.

［14］Huy Tho Ho，Gibbins Danny. Multi-scale feature extraction for 3D models using local surface curvature [C]//Digital Image Computing：Techniques and Applications.2008：16-23.

［15］Cazals F，Pouget M. Estimating differential quantities using polynomial fitting of osculating jets[J]. Computer Aided Geometric Design，2005，22(2)：121-146.

［16］Lange C，Polthier K. Anisotropic smoothing of point sets[J]. Compute Aided Geometric Design ， 2005 ，22(7)：680-692.

［17］Alexa M，Behr J，Cohen-Or D，et al. Point set surfaces[C]// Proceedings of IEEE Visualization. 2001：21-28.

［18］Stanford University. Stanford 3D Scanning Repository[EB/OL]. http：//www-graphics. stanford. edu/data/3Dscanrep/.

［19］Arun K S，Huang T S，Blostein S D. Least-squares fitting of two 3-D point sets[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence ， 1987 ，9(5)：698-700.