一種改進的最大一致性點云幾何基元擬合算法

劉修國，楊 準，王紅平，梁 棟

（中國地質大學（武漢）信息工程學院，湖北 武漢430074）

平面、二次曲面（球、圓柱、圓錐）等基本幾何基元的高精度擬合是點云數據處理的基礎工作.高精度幾何基元擬合獲取幾何基元可靠參數，可實現標靶高精度定位，完成基于標靶［1］或基于幾何基元［2］的多測站配準；也有助于點云分割［3-5］中過分割、欠分割問題的解決，并有效支撐基于點云數據的表面建模［6-7］.

近年來，幾何基元的高精度擬合仍是點云數據處理的研究熱點之一.文獻［8-9］采用了基于特征值最小二乘的擬合方法，用于解決含有粗差的點云數據平面基元高精度擬合.文獻［10-11］采用整體最小二乘方法擬合幾何基元，該方法顧及點云數據在X，Y，Z3個方向存在的誤差，擬合精度有所提高.文獻［12］提出迭代平方距離函數最小化的曲面擬合方法，用平方距離度量點集與曲面基元的幾何距離，用于擬合曲面基元.對于預先剔除粗差點后的點云數據，文獻［10-12］的擬合方法能較高精度地擬合出平面或簡單曲面幾何基元.目前應用廣泛且較成熟的是 隨 機 抽 樣 一 致 （random sample consensus，RANSAC）幾何基元擬合算法［3-5，13］；該算法支持對平面、二次曲面等多種幾何基元的擬合，能在含有一定粗差的點云中可靠地擬合出幾何基元.但RANSAC算法不能自適應地設定距離閾值以區(qū)分粗差點與非粗差點：閾值過大則不能完全剔除粗差，造成擬合精度不高甚至擬合失?。婚撝颠^小會剔除較多非粗差點，難以保證特定幾何基元點集的完整性［14］，降低擬合精度.為此，文獻［14］提出了 MCMD（maximum consistency with minimum distance and robust Z-score）算法.該算法基于主成分分析（principal component analysis，PCA）［15］選取可靠平面初值，再采用穩(wěn)健Z分數法一次性剔除粗差，有效避免了距離閾值設置以區(qū)分粗差、非粗差點的問題.但MCMD算法基于PCA的初值選取準則僅適于平面基元，不適于二次曲面基元；同時，在點云數據粗差統(tǒng)計分布未知的復雜場景中，穩(wěn)健Z分數法難以一次性剔除全部粗差，限制了該算法的適用性.

1 MCMD平面基元擬合算法思想

在獲取可靠的平面初值后，采用穩(wěn)健Z分數方法來剔除粗差，設點集Q中任意點qi∈Q到初始平面的垂直距離為di，則該點的穩(wěn)健Z分數Zr，j定義為

式中：di，med表示距離中位數；di，mad為距離的中位數中誤差，其定義為

當Zr，j≥k0（2.0≤k0≤2.5）時，可認為qi點是粗差點［16］.剔除掉點集Q中所有Zr≥k0的粗差點后，剩余點認為是屬于同一平面基元的點，對其進行擬合獲取最終的平面基元參數.

2 改進MCMD幾何基元擬合算法

本文算法首先構建點子集的候選集合.設采用隨機抽樣思想從點集Q中產生的候選集合中有T個點子集，設點集Q中某一個點為非粗差點的概率為ε，確定待擬合基元形狀參數的最少點個數為N0，循環(huán)T次獲取的T個點子集樣本中至少有一個樣本不含粗差的概率為Pr，則有1-Pr=（1-εN0）T，即

采用距離和最小初值選取準則從點子集的候選集合中選取最佳點子集，可靠的模型初值為最佳點子集對應的擬合幾何基元參數.

改進的算法流程見圖1.

2.1 距離和最小的初值選取準則

距離和最小的初值選取準則，即點子集到擬合模型的距離和越小則點共表面程度越大.以距離和最小作為最優(yōu)模型初值的選取準則，擴展MCMD算法使之適用二次曲面基元的擬合.

距離和是指選取的點子集中每個點到擬合模型距離的總和.對平面而言，距離和是點到平面模型d=ax+by+cz的垂直距離和，其計算公式為

式中：Dsum表示距離和；a，b，c，d均為平面參數；h為點子集的數目.

對于二次曲面（球體、圓柱、圓錐）而言，距離和是指點子集中每個點到曲面的平方距離和.與垂直距離相比，平方距離能更加精確地表示點與擬合模型的幾何距離［12，17］.平方距離和的計算公式為

圖1 改進的MCMD算法流程Fig.1 Flowchart of improved MCMD algorithm

式中：對于三維表面S，eSD，i為第i個點到模型表面的平方距離；si代表S上與qi最近的點，si處的曲率半徑為ρi，1和ρi，2；ni，1，ni，2為主曲率方向向量，ni，3為法向量；最短距離為.若曲率中心與點qi在表面S的兩側，則di＜0；反之di＞0，當0＜di＜ρi，j或者ρi，j為無窮大時，αi，j=0.

2.2 粗差剔除與模型參數求解

在復雜目標場景中，粗差點統(tǒng)計分布規(guī)律是未知的，此時粗差難以一次性剔除.本文算法采用穩(wěn)健Z分數方法循環(huán)剔除粗差策略逐步剔除粗差，即每次循環(huán)依據距離和最小準則獲取可靠初始模型，重新計算垂直距離中位數di，med、中位數中誤差di.mad，并計算點集Q中每個點的Zr值；將Zr≥k0的點從點集中剔除，當點集中所有點的Zr＜k0時停止循環(huán).

對剔除粗差后的點集采用加權最小二乘法迭代優(yōu)化獲取最終的幾何基元形狀參數.

對平面基元：d=ax+by+cz，待求參數為單位法向量n=（a，b，c），原點到平面的距離d.點集中任意點qi=（xi，yi，zi）到平面的有向距離為

對球面基元：（xi-x0）2+（yi-y0）2+（zi-z0）2=r，待求參數為球中心o=（x0，y0，z0）、球半徑r.則點集中任意點qi=（xi，yi，zi）到球面的有向距離為

對圓柱基元：（xi-x0）2+（yi-y0）2+（zi-z0）2-［nx（xi-x0）+ny（yi-y0）+nz（zi-z0）］2=r2，待求參數為圓柱單位軸向n=（nx，ny，nz）、軸向上一點o=（x0，y0，z0）及圓柱半徑r.則點集中任意點qi=（xi，yi，zi）到圓柱的有向距離為

式中：Pi為權函數；D（k）i為第k次迭代點到幾何基元的有向距離；D（k）為距離的范圍.

設置迭代閾值δ0，當兩次迭代獲取的對應參數差值小于δ0，停止迭代獲取幾何基元形狀參數.

3 實驗結果與分析

3.1 平面擬合實驗

首先采用模擬數據評估算法的粗差識別能力及擬合精度.模擬的平面方程為2=x+y+z，坐標（x，y）在區(qū)間［0，1］上服從均勻分布，坐標z的值由模擬的平面方程計算得到.在模擬點坐標（x，y，z）均加上誤差ν～N（0，0.002）形成非粗差點.模擬兩種粗差點分布（圖2）：①粗差點分布在模型一側（分布A），即對模擬點坐標（x，y，z）加上均值分別為（0.8，0.9，1.0）、方差均為0.5的高斯分布；②粗差點分布在模型兩側（分布B），即一半粗差點與分布A相同，另一半粗差點則是對模擬點坐標（x，y，z）加上均值分別為（-0.8，-0.9，-1.0）、方差均為0.5的高斯分布.

圖2 兩種粗差分布對應的模擬數據（粗差比率為20%）Fig.2 Simulated point cloud corresponding to distribution of two kinds of outliers（outlier rate is 20%）

對于兩種粗差分布，在不同粗差比率（10%，20%，30%，40%，50%）下分別模擬產生1 000組數據集，每組數據集中包含1 000個點.利用本文算法、MCMD和RANSAC算法對模擬數據進行平面基元擬合.針對RANSAC難以設定距離閾值，本文通過設置不同距離閾值進行RANSAC算法實驗，從中選取最優(yōu)結果作為對比；本模擬數據中RANSAC算法距離閾值設為3mm.

采用正確識別比率（correct identificaiton rate，CIR）、錯誤識別比率（swamping rate，SR）兩個指標［14］來評估擬合算法能否正確區(qū)分粗差點、非粗差點.擬合精度采用擬合平面與真實平面法向量夾角θ以及坐標原點到擬合平面距離與真實距離值間的差值ΔD評估.CIR指標為識別的粗差點中真實粗差點個數與總的粗差點個數之比，SR指標為識別的粗差點中非粗差點的個數與總的非粗差點個數之比.3種算法的CIR與SR均值統(tǒng)計見表1，擬合精度指標如圖3與圖4所示.

表1 平面擬合的CIR與SR指標Tab.1 CIR and SR of plane fitting

從表1可以看出，RANSAC算法的CIR指標為100%，即能剔除全部粗差.但該算法SR指標均大于30%，破壞了平面基元點的完整性，也造成平面擬合精度較低.MCMD算法SR指標均小于0.4%，該算法的平面點集完整性與平面擬合精度較

RANSAC顯著提高，但是其在粗差比率為50%的分布B模擬數據時CIR為98.8%，不能剔除所有粗差，這與第1節(jié)分析得出的MCMD算法局限②相吻合.本文算法CIR指標全部為100%，能夠完全剔除粗差.而SR指標均低于1.6%，略高于MCMD算法.這是因為本文算法采用了循環(huán)穩(wěn)健Z分數剔除粗差策略，使得在復雜場景中能完全剔除粗差且盡可能保證屬于同一幾何基元點集的完整性.

圖3 粗差分布A對應的平面擬合精度Fig.3 Plane fitting accuracy corresponding to distribution of outlier A

圖4 粗差分布B對應的平面擬合精度Fig.4 Plane fitting accuracy corresponding to distribution of outlier B

從圖3與圖4可以看出，3種算法中RANSAC算法擬合精度最低，而本文算法與MCMD算法擬合精度相當.值得注意的是，對粗差比率為50%的分布B模擬數據，MCMD算法不能剔除全部粗差，使得擬合失效.

采用真實數據驗證本文算法對目標多樣、存在遮擋、點密度變化顯著的復雜場景的粗差剔除能力.如圖5所示（第1行為正視圖，第2行為頂視圖），點云數據中主墻平面為待擬合平面；場景中存在遮擋主墻平面的電纜線（圖5標號1）、磚塊（圖5標號2）、樹木（圖5標號3）等不同目標；存在位于樹葉之后點密度稀疏的墻面區(qū)域（圖5標號4）.對此復雜場景，分別采用MCMD，RANSAC與本文算法進行粗差剔除，其中RANSAC算法距離閾值設置為3 mm.各算法保留的主墻平面基元點集結果見圖5，法向夾角θ與距離偏差ΔD統(tǒng)計結果見表2.

表2 墻平面擬合結果Tab.2 Results of facade plane fitting

3.2 球面擬合實驗

采用地面三維激光掃描儀獲取的標靶球進行擬合并評估其精度.將獲取的數據場景中5個標靶球（圖6）分割出來，對5個標靶分別進行擬合，其中RANSAC算法距離閾值設為3.5mm，相關的擬合結果見表3.其中識別率為算法自動識別的球面點與人工識別的球面點的比值；球心誤差為算法自動擬合球心與人工識別的點擬合得到的球心歐氏距離.

圖5 復雜場景下剔除粗差后的墻平面Fig.5 Plane data after romving outliers in complex scene

圖6 真實場景標靶球數據Fig.6 Sphere target data from scanning scene

表3 標靶球擬合結果Tab.3 Results of sphere target fitting

從擬合結果來看，RANSAC算法過多地剔除了非粗差點，擬合精度不高；而本文算法擬合結果較好，即使在獲取的標靶球點數目較少的情況下，本文算法的擬合結果也與人工識別的結果相當.

3.3 圓柱與圓錐擬合實驗

對地面三維激光掃描儀獲取的圓柱進行擬合并評估擬合精度.圓柱基元的參數包括軸線上的一個點、軸線方向及其軸半徑.采用RANSAC算法與本文算法對圓柱點云進行擬合，其中RANSAC算法距離閾值設為4mm，相應的擬合結果如圖7與表4所示.其中識別率為算法自動識別的圓柱點與人工識別的圓柱點的比值；軸線距離為算法自動擬合軸線與人工擬合軸線間的最小距離；軸線誤差為算法自動擬合軸線與人工擬合軸線間夾角；半徑誤差為算法自動擬合半徑與人工擬合半徑差的絕對值.

表4 圓柱擬合結果Tab.4 Results of cylinder fitting

圖7 粗差剔除后屬同一圓柱基元的點集Fig.7 Cylinder data after removing outliers

對圓錐擬合實驗則采用頂點為（0，0，0）、軸向為（0，0，1.0）、頂角為90°的圓錐模擬數據.圓錐模擬數據生成方法與平面模擬數據生成類似，模擬驗證數據粗差比率為30%.擬合結果見圖8與表5，其中RANSAC算法距離閾值設為2.5mm.頂點誤差為算法擬合頂點與模擬模型頂點的歐式距離，頂角誤差為算法擬合頂角與模擬模型頂角差值.

表5 圓錐擬合結果Tab.5 Results of cone fitting

圖8 粗差剔除后屬同一圓錐基元的點集Fig.8 Cone data after removing outliers

從圓柱、圓錐基元擬合結果可以看出，本文算法同樣適用，且能有效剔除粗差，保證較高的識別率，獲得高精度的基元參數，各項擬合誤差均比RANSAC算法擬合誤差小.

4 結論

［1］ 陳西江，花向紅，楊榮華，等.分帶K-均值聚類的平面標靶定位［J］.武漢大學學報：信息科學版，2013，28（2）：167.CHEN Xijiang，HUA Xianghong，YANG Ronghua，etal.Planar target location based on the zoning K-means clustering［J］.Geomatics and Information Science of Wuhan University，2013，28（2）：167.

［2］ Theiler P W，Schindler K.Automatic registration of terrestrial laser scanner point clouds using natural planar surfaces［C］∥XXII ISPRS Congress，Technical Commission III.Melbourne：Copernicus Publications，2012：173-178.

［3］ PU Shi，Vosselman G.Automatic extraction of building features from terrestrial laser scanning［J］.International Archives of Photogrammetry，Remote Sensing and Spatial Information Sciences，2006，36（5）：25.

［4］ Schnabel R，Wahl R，Klein R.Efficient RANSAC for pointcloud shape detection［J］.Computer Graphics Forum，2007，26（2）：214.

［5］ Zhang G，Karasev P，Brilakis I，etal.A sparsity-inducing optimization algorithm for the extraction of planar structures in noisy point-cloud data［C］∥ Proceedings of the 2012 ASCE International Conference on Computing in Civil Engineering.Clearwater Beach：ASCE Publications，2012：317-324.

［6］ Kaucic R， Hartley R， Dano N.Plane-based projective reconstruction［C］∥ Eighth IEEE International Conference on Computer Vision.Vancouver：IEEE，2001：420-427.

［7］ NAN Liangliang，Sharf A，ZHANG Hao，etal.Smart boxes for interactive urban reconstruction［J］.ACM Transactions on Graphics，2010，29（4）：1.

［8］ 官云蘭，程效軍，施貴剛.一種穩(wěn)健的點云數據平面擬合方法［J］.同濟大學學報：自然科學版，2008，36（7）：981.GUAN Yunlan，CHENG Xiaojun，SHI Guigang.A robust method for fitting aplane to point clouds［J］.Journal of Tongji University：Natural Science，2008，36（7）：981.

［9］ 王峰，丘廣新，程效軍.改進的魯棒迭代最小二乘平面擬合算法［J］.同濟大學學報：自然科學版，2011，39（9）：1350.WANG Feng，QIU Guangxin，CHENG Xiaojun.An improved robust method for iterating least-squares plane fitting［J］.Journal of Tongji University：Natural Science，2011，39（9）：1350.

［10］ 魯鐵定，周世健，張立亭，等.基于整體最小二乘的地面激光掃描標靶球定位方法［J］.大地測量與地球動力學，2009，29（4）：102.LU Tieding，ZHOU Shijian，ZHANG Liting，etal.Sphere target fixing of point cloud data based on TLS［J］.Journal of Geodesy and Geodynamics，2009，29（4）：102.

［11］ 官云蘭，劉紹堂，周世健，等.基于整體最小二乘的穩(wěn)健點云數據平面擬合［J］.大地測量與地球動力學，2011，31（5）：80.GUAN Yunlan，LIU Shaotang，ZHOU Shijian，etal.Robust plane fitting of point clouds based on TLS［J］.Journal of Geodesy and Geodynamics，2011，31（5）：80.

［12］ WANG Jun，YU Zeyun.Quadratic curve and surface fitting via squared distance minimization［J］.Computers &Graphics，2011，35（6）：1035.

［13］ Fischler M A，Bolles R C.Random sample consensus：A paradigm for model fitting with applications to image analysis and automated cartography［J］.Communications of the ACM，1981，24（6）：381.

［14］ Nurunnabi A，West G，Belton D.Robust outlier detection and saliency features estimation in point cloud data［C］∥ 10th International Conference on Computer and Robot Vision（CRV）.Regina：IEEE，2013：98-105.

［15］ Hoppe H，DeRose T，Duchamp T，etal.Surface reconstruction from unorganized points［C］∥ Proceedings of the 19th Annual Conference on Computer Graphics and Interactive Techniques.New York：ACM Press，1992：71-78.

［16］ 楊元喜.抗差估計理論及其應用［M］.北京：八一出版社，1993.YANG Yuanxi.Robust estimation theory and its applications［M］.Beijing：Bayi Publishing House，1993.

［17］ Pottmann H，Hofer M.Geometry of the squared distance function to curves and surfaces［C］∥3rd International Workshop on Visualization and Mathematics.Berlin：Springer，2002：221-242.