海量點(diǎn)云數(shù)據(jù)輪廓特征線的快速生成算法

程效軍，賈東峰，劉燕萍

（1.同濟(jì)大學(xué) 測量與國土信息工程系，上海200092；2.同濟(jì)大學(xué) 現(xiàn)代工程測量國家測繪地理信息局重點(diǎn)實(shí)驗(yàn)室，上海200092；3.同濟(jì)大學(xué) 浙江學(xué)院 土木工程系，浙江 嘉興314000）

隨著點(diǎn)云數(shù)據(jù)規(guī)模的不斷增大，如何有效地結(jié)合逆向工程（RE）和快速成型技術(shù)（RP），從海量的點(diǎn)云數(shù)據(jù)中獲取平面輪廓特征點(diǎn)，并進(jìn)行輪廓線的重建，是實(shí)現(xiàn)基于特征的模型重建的關(guān)鍵.Lee等[1]根據(jù)點(diǎn)的提取率對點(diǎn)云進(jìn)行分層，提出角偏差法來獲取每層數(shù)據(jù)的曲率突變點(diǎn)，并通過同倫算法生成點(diǎn)云的輪廓線模型.由于曲率突變點(diǎn)的獲取容易受到噪聲的影響，所以該方法抗噪能力差.Wu等[2]提出了基于RE和RP的迭代算法，通過定義形狀因子，計(jì)算每層數(shù)據(jù)的形狀誤差，以此來調(diào)節(jié)點(diǎn)云分層的厚度，實(shí)現(xiàn)分層劃分的自適應(yīng)，但是其輪廓特征線的生成效率較低.Liu等[3]提出了一種基于中間點(diǎn)的曲線模型法（IPCM）來集成RE和RP.該方法效率較高，但未解決點(diǎn)云切片中的多環(huán)現(xiàn)象.國內(nèi)的柯映林等[4]提出了基于點(diǎn)集的“緊包圍盒”分割來建立散亂點(diǎn)之間的鄰接信息，構(gòu)建平面切片數(shù)據(jù)，并根據(jù)散亂點(diǎn)間的拓?fù)潢P(guān)系，重建平面曲線.吳杭彬等[5]根據(jù)點(diǎn)云數(shù)據(jù)的特點(diǎn)，提出了點(diǎn)云數(shù)據(jù)的分層模型，實(shí)現(xiàn)了等值線的分層提取，并根據(jù)等值線的細(xì)節(jié)程度不同，對等值線進(jìn)行分類，實(shí)現(xiàn)了多細(xì)節(jié)表達(dá).

以上方法，通過集成RE和RP，避免了曲面擬合和STL文件生成的過程[6－7]，其應(yīng)用領(lǐng)域多為機(jī)械制造，目的是把曲線模型作為RP系統(tǒng)的輸入.本文結(jié)合逆向工程和快速成型兩種技術(shù)，針對海量散亂點(diǎn)云，提出一種基于切片的海量點(diǎn)云數(shù)據(jù)輪廓特征線快速生成算法，并通過實(shí)例分析了算法的效率和可行性.

1 算法描述

1.1 點(diǎn)云切片的生成

設(shè)有散亂點(diǎn)集S＝｛p1，p2，…，pn｝，pi＝｛xi，yi，zi｝∈R3，則點(diǎn)云切片的生成可描述為：采用一組平行平面沿給定方向?qū)θS點(diǎn)云進(jìn)行劃分，并將三維點(diǎn)集轉(zhuǎn)化為二維切片數(shù)據(jù).點(diǎn)集S的坐標(biāo)范圍為（xmin，ymin，zmin）～（xmax，ymax，zmax），有一平面集Τ，平行于xoy平面，法矢n指向Z軸正方向，該平面集Τ由一組Z坐標(biāo)序列表示：Z＝（Z0，Z1，…，Zn）且Z0＜Z1＜…＜Zn，其中：

通過定義參考平面位于分層數(shù)據(jù)的中間位置，將每層數(shù)據(jù)向參考平面投影即可得到切片數(shù)據(jù).

1.2 特征點(diǎn)提取

本文采用基于數(shù)字圖像的方法提取切片數(shù)據(jù)的特征點(diǎn)，該方法由 Chen等[8]和 Zhang等[9]提出.首先建立一個(gè)柵格化的數(shù)字平面（取參考平面），將每層點(diǎn)云數(shù)據(jù)投影到數(shù)字平面上，對于落入點(diǎn)的柵格賦于值“1”（黑塊），對于未落入點(diǎn)的柵格賦于值“0”（白塊）.如圖1a為將點(diǎn)云數(shù)據(jù)投影到柵格化的數(shù)字平面；圖1b為投影到數(shù)字平面的第24層切片.

式中，Zpitch是分層厚度.本文通過獲取對象的高度H，計(jì)算分層厚度.為了控制切片中投影帶的寬度，點(diǎn)云的分層厚度計(jì)算如下：

圖1 將空間點(diǎn)投影到柵格化的數(shù)字平面Fig.1 Projecting spatial points onto digital grid plane

文獻(xiàn)[3，6，10]采用基于細(xì)化模板的特征點(diǎn)提取算法來提取投影后的數(shù)據(jù)點(diǎn).該方法由于未考慮柵格的邊長對特征提取的影響，會造成特征點(diǎn)的刪除，因此使用上述方法提取特征點(diǎn)后，需采用一個(gè)基于最大距離誤差的算法來判斷刪除的冗余點(diǎn)中是否含有特征點(diǎn)，并根據(jù)距離判斷以恢復(fù)刪除的特征點(diǎn)，該過程會影響特征點(diǎn)提取的效率.

為了提高平面輪廓特征點(diǎn)提取的效率，在建立柵格化的參考平面時(shí)，可以通過討論柵格邊長對特征提取的影響，提高輪廓特征點(diǎn)提取的效率，該方法的具體過程如下.

1.2.1 數(shù)字平面柵格邊長的計(jì)算

平面柵格邊長的大小直接影響到特征點(diǎn)的提取和特征線的構(gòu)造精度.為了減少平面柵格邊長對特征提取的影響，提高算法的穩(wěn)健性，本文采用以下方法計(jì)算平面柵格的邊長.

首先，計(jì)算平面點(diǎn)集的最小包圍盒，獲取平面點(diǎn)集的坐標(biāo)范圍（xmin，ymin）～（xmax，ymax），柵格的邊長size計(jì)算公式如下：

式中：α是尺度因子，用于調(diào)節(jié)柵格邊長的大小，n是點(diǎn)的個(gè)數(shù).文獻(xiàn)[11]討論了α的最佳取值范圍是[1，1.5]，本文取α的值為1.2.

1.2.2 特征點(diǎn)的提取

將空間點(diǎn)投影到柵格化的參考平面，以柵格中的點(diǎn)云數(shù)據(jù)重心取代柵格內(nèi)的點(diǎn)作為特征點(diǎn)，則第i個(gè)柵格內(nèi)的特征點(diǎn)的坐標(biāo)計(jì)算公式為：

如圖2a為將空間點(diǎn)投影到按照式（3）計(jì)算的柵格平面中；圖2b為提取每一個(gè)柵格中的點(diǎn)云重心作為輪廓特征點(diǎn).

圖2 特征點(diǎn)的提取Fig.2 The extraction of feature points

1.3 輪廓特征線的生成

如何根據(jù)提出的特征點(diǎn)快速構(gòu)建輪廓特征線是本文研究的重點(diǎn)，由于分層厚度和柵格邊長的設(shè)置可以有效地控制切片投影帶的寬度，根據(jù)平面柵格中投影帶的特點(diǎn)，本文提出了雙向索引連通法來快速構(gòu)建特征線，方法原理介紹如下.

1.3.1 雙向索引連通法

首先按X方向建立柵格的索引，判斷柵格內(nèi)是否含有特征點(diǎn)，分層構(gòu)造特征線段；再按Y方向建立柵格索引，判斷柵格內(nèi)是否含有特征點(diǎn)，分層構(gòu)造特征線段，圖3所示是按雙向索引連通法構(gòu)造特征線的原理圖；圖3a顯示的是柵格平面的特征點(diǎn)；圖3b和圖3c是分別按X和Y方向根據(jù)柵格的索引分層連接特征點(diǎn)；圖3d是輪廓特征線的生成圖.

圖3 雙向索引連通法生成輪廓特征線Fig.3 Contour generation using bidirectional index connecting method（BICM）

1.3.2 對角鄰域、田字形鄰域以及鄰域缺失的連通

基于雙向索引連通法能快速連接特征點(diǎn)的四鄰域，但是對于如圖4所示的對角鄰域、田字形鄰域及鄰域缺失的情況造成輪廓特征線的斷裂和偏轉(zhuǎn)，需要設(shè)計(jì)特定的判斷條件進(jìn)行連通.圖中，P，pi表示點(diǎn)位.

圖4 幾種鄰域類型Fig.4 Several types of neighborhood

首先根據(jù)X方向的索引存儲每層特征線段的首尾特征點(diǎn)（端點(diǎn)），判斷這些端點(diǎn)含有四鄰域點(diǎn)的個(gè)數(shù).文中采用的存儲容器為C＋＋標(biāo)準(zhǔn)模板庫中的vector容器.

對角鄰域的連通：從端點(diǎn)容器中任意取出一點(diǎn)，如果其含有兩個(gè)四鄰域點(diǎn)，則從容器中刪除該點(diǎn)；如果小于等于一個(gè)四鄰域點(diǎn)則保存該點(diǎn)到對角鄰域容器，循環(huán)判斷直到遍歷所有端點(diǎn).判斷完畢，從對角鄰域容器中任意選出一點(diǎn)，根據(jù)索引判斷特征點(diǎn)容器中是否含有該點(diǎn)的對角鄰域點(diǎn)，如果有則連接兩點(diǎn)，然后從容器中剔除這兩個(gè)點(diǎn)；如果沒有對角鄰域點(diǎn)，則保存該點(diǎn)到鄰域缺失容器.

鄰域缺失的連通：經(jīng)過對角鄰域的連通，容器中剩余的點(diǎn)為鄰域缺失的點(diǎn)，任意取出一點(diǎn)，判斷該點(diǎn)到其他端點(diǎn)（可選擇相鄰層的端點(diǎn)）的距離，由于切片數(shù)據(jù)存在不閉合情況，如建筑點(diǎn)云切片的門窗處形成點(diǎn)云的斷裂，要設(shè)定距離閾值進(jìn)行判斷，本算法設(shè)置2size作為距離閾值，如果端點(diǎn)距離大于距離閾值則認(rèn)為是自然斷裂無需連接，如果小于距離閾值則選擇并連接離該點(diǎn)最近的端點(diǎn)，完成鄰域缺失的連通.

田字形鄰域的連通：為了準(zhǔn)確找到田字形鄰域的位置，首先需要判斷特征點(diǎn)的鄰域關(guān)系，如果點(diǎn)P存在如圖5所示的情況（Ynum為Y方向上的序號），則存儲該點(diǎn)和其鄰域值，循環(huán)判斷直到遍歷所有特征點(diǎn)；找到田字形鄰域后，需要對田字形鄰域進(jìn)行合并，本文采用均值合并的方法對其進(jìn)行合并，由于可按上下鄰域取平均或左右鄰域取平均，算法首先分別按兩個(gè)方向取平均，判斷某方向合并后的均值點(diǎn)是否存在兩個(gè)四鄰域點(diǎn)，如果存在則用其替代“田”字鄰域中的4個(gè)點(diǎn)，并刪除容器中的四鄰域點(diǎn)，然后對均值點(diǎn)進(jìn)行連接.

圖5 田字形鄰域的判斷條件Fig.5 The judging condition of“田”shaped neighborhood

2 實(shí)例分析

本文在Visual C＋＋6.0環(huán)境下結(jié)合COIN3D庫測試上述算法，測試用的點(diǎn)云數(shù)據(jù)，如圖6所示，圖6a為標(biāo)準(zhǔn)構(gòu)件的點(diǎn)云數(shù)據(jù)（點(diǎn)云個(gè)數(shù)4102）；圖6b為建筑物的點(diǎn)云數(shù)據(jù)（點(diǎn)云個(gè)數(shù)1 985 617）；圖6c為華佗雕像的點(diǎn)云數(shù)據(jù)（點(diǎn)云個(gè)數(shù)708 987）.為了說明算法的可行性，本文針對以上3種點(diǎn)云數(shù)據(jù)進(jìn)行實(shí)驗(yàn).

第一組實(shí)驗(yàn)：測試算法對對角鄰域、田字形鄰域以及鄰域數(shù)據(jù)缺失的連通情況.如圖7a所示為標(biāo)準(zhǔn)構(gòu)件的第2層點(diǎn)云切片；圖7b是點(diǎn)云切片生成輪廓線過程中的對角鄰域造成了輪廓線的斷裂；圖7c是根據(jù)本算法提出的方法對對角鄰域進(jìn)行連通，從圖中可以看出算法判斷出了對角鄰域點(diǎn)并對其進(jìn)行了連通.

圖6 點(diǎn)云數(shù)據(jù)Fig.6 Point cloud

圖7 對角鄰域的連通Fig.7 The connecting of diagonal neighborhood

圖8a所示為標(biāo)準(zhǔn)構(gòu)件的第5層點(diǎn)云切片；圖8b是輪廓線生成過程中的鄰域缺失造成輪廓線的斷裂；圖8c是根據(jù)本算法提出的方法實(shí)現(xiàn)了鄰域缺失情況的連通.

圖8 鄰域數(shù)據(jù)缺失的連通圖Fig.8 The connecting in case of missing neighbor data

圖9a所示是華佗雕像的第16層點(diǎn)云切片；圖9b是根據(jù)本算法計(jì)算并提取切片數(shù)據(jù)的特征點(diǎn)；圖9c由于田字形鄰域造成輪廓線的局部偏轉(zhuǎn)；圖9d根據(jù)本算法判斷出田字形鄰域的位置；圖9e根據(jù)均值合并的方法對田字形鄰域進(jìn)行合并；圖9f輪廓線的生成，可以看出本算法已經(jīng)消除了田字形鄰域造成的輪廓線偏轉(zhuǎn)，實(shí)現(xiàn)了輪廓線的連通.

第二組實(shí)驗(yàn)：本算法生成的輪廓線和逆向工程軟件Geomagic以及導(dǎo)入計(jì)算機(jī)輔助設(shè)計(jì)軟件（CAD）手繪產(chǎn)生的輪廓線效果進(jìn)行比較.如圖10a所示是某建筑物的第7層點(diǎn)云切片；圖10b是根據(jù)本文提出的算法生成的輪廓線，準(zhǔn)確表達(dá)了點(diǎn)云切片數(shù)據(jù)的輪廓特征；圖10c是逆向工程軟件Geomagic生成的輪廓線，圖中有許多斷裂；圖10d是將點(diǎn)云切片導(dǎo)入CAD后通過手工繪制的輪廓線，雖然可以表達(dá)點(diǎn)云切片的輪廓特征，但是人工判斷以及手工連接需要花費(fèi)大量的時(shí)間.通過比較可以看出本算法比其他商業(yè)軟件可以更準(zhǔn)確地表達(dá)出點(diǎn)云數(shù)據(jù)的輪廓特征.

圖9 田字形鄰域的連通Fig.9 The connecting of“田”shaped neighborhood

圖10 輪廓線的生成效果比較Fig.10 The comparison of contour generation of different methods

第三組實(shí)驗(yàn)，本算法和凸包算法的時(shí)間復(fù)雜度分析比較.Graham[12]提出的平面點(diǎn)集凸包Graham掃描算法，其時(shí)間復(fù)雜度為O（nlgn）.文獻(xiàn)[13]提出了平面點(diǎn)集凸包的最優(yōu)實(shí)時(shí)算法，使其復(fù)雜度達(dá)到了O（n）.本算法根據(jù)X、Y雙向搜索連接，并對局部的連通性進(jìn)行判斷，其時(shí)間復(fù)雜度是線性的，可以達(dá)到O（n）.圖11是采用本算法生成的雕像點(diǎn)云的輪廓線模型.圖11a，b，c分別是華佗點(diǎn)云的正視、前視、俯視輪廓線模型.

圖11 華佗輪廓線模型Fig.11 The contour model of Huatuo statue

3 結(jié)論

在針對海量散亂點(diǎn)云的算法中，基于點(diǎn)云特征的表面重建是近年研究的熱點(diǎn)，如何從海量散亂點(diǎn)云數(shù)據(jù)中提取特征點(diǎn)構(gòu)造特征線，并基于特征線實(shí)現(xiàn)海量點(diǎn)云的表面重建是算法研究的關(guān)鍵.

本文通過結(jié)合逆向工程和快速成型技術(shù)提出了一種快速構(gòu)造點(diǎn)云輪廓特征線的方法.該方法：（1）采用數(shù)字圖像的方法，通過建立數(shù)字柵格平面提取點(diǎn)云切片的特征點(diǎn)；（2）量化了數(shù)字平面的柵格邊長并以落入柵格中的點(diǎn)云重心作為特征點(diǎn)；（3）基于特征點(diǎn)在柵格中的索引值，提出了雙向連通索引法快速構(gòu)造點(diǎn)云特征線.通過和商業(yè)軟件以及凸包算法的分析和比較，本算法可以快速、準(zhǔn)確地生成點(diǎn)云的輪廓曲線模型.進(jìn)一步研究的重點(diǎn)是討論柵格邊長的大小和曲線模型的精度的關(guān)系，并建立二者之間的數(shù)學(xué)模型.

[1] Lee K H，Woo H.Direct integration of reverse engineering and rapid prototyping[J].Computer Industry Engineering，2000，38：21.

[2] WU Y F，WONG Y S，Loh H T，et al.Modeling cloud data using an adaptive slicing approach [J].Computer－Aided Design，2004，36：231.

[3] Liu G H，Wong Y S，Zhang Y F，et al.Error based segmentation of cloud data for direct rapid prototyping [J].Computer－Aided Design，2002，35：633.

[4] 柯映林，王青.反求工程中的點(diǎn)云切片算法研究[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2005，17（8）：1798.KE Yinglin，WANG Qing.Research on point cloud slicing technique in reverse engineering [J].Journal of Computer－Aided Design ＆Computer Graphics，2005，17（8）：1798.

[5] 吳杭彬，劉春.激光掃描數(shù)據(jù)的等值線分層提取和多細(xì)節(jié)表達(dá)[J].同濟(jì)大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，37（2）：267.WU Hangbin，LIU Chun.Point cloud－based stratified contour extraction and its multi－LOD expression with ground laser range scanning [J].Journal of Tongji University：Natural Science，2009，37（2）：267.

[6] Kumbhar V K，Pandey P M，Rao P V M.Improved intermediate point curve model for integrating reverse engineering and rapid prototyping [J].The International Journal of Advanced Manufacturing Technology，2007，37：553.

[7] YIN Zhongwei.Direct integration of reverse engineering and rapid prototyping based on properties of NURBS or B－spline[J].Precision Engineering，2004，28：293.

[8] CHEN Y H，NG C T，WANG Y Z.Data reduction in integrated reverse engineering and rapid prototyping [J].International Journal of Computer Integrated Manufacturing，1999，12：97.

[9] ZHANG Y F，WONG Y S，Loh H T，et al.An adaptive slicing approach to modeling cloud data for rapid prototyping [J].Mater Process Technology，2003，140：105.

[10] Jang B K，Chin R T.Reconstruc
Table parallel thinning[J].International Journal of Pattern Recognition Artificial Intelligence，1993，7（5）：1145.

[11] Piegl L A，Tiller W.Algorithm for finding all k nearest neighbors[J].Computer－Aided Design，2002，34：167.

[12] Graham R L.An efficient algorithm for determine the convex hull of a finite linear set[J].Information Processing Letters，1972（1）：132.

[13] 王志強(qiáng)，洪嘉振，肖立瑾.平面點(diǎn)集凸包的最優(yōu)實(shí)時(shí)算法[J].計(jì)算機(jī)學(xué)報(bào)，1998，8（21）：351.WANG Zhiqiang，HONG Jiazhen，XIAO Lijin.An optimal real time algorithm for determining the convex hull of a set of points in a plane[J].Chinese Journal of Computers，1998，8（21）：351.