張習(xí)民, 余小清， 萬(wàn)旺根， 張 娟

1.上海大學(xué)通信與信息工程學(xué)院，上海200444

2.上海大學(xué)智慧城市研究院，上海200444

3.河南教育學(xué)院物理與電子工程學(xué)院，鄭州450046

三維激光掃描技術(shù)以其非接觸性、高精度等特點(diǎn)廣泛應(yīng)用于空間信息的獲取，但掃描得到的海量點(diǎn)云數(shù)據(jù)不僅占用了大量的存儲(chǔ)空間，而且不利于傳輸及后期運(yùn)算處理.因此，在保證物體的形狀精度和數(shù)據(jù)處理效率的前提下有必要減少點(diǎn)云的數(shù)據(jù)量.目前，專家學(xué)者對(duì)點(diǎn)云精簡(jiǎn)及壓縮方面的研究主要集中于網(wǎng)格壓縮法和基于空間點(diǎn)模型的壓縮法.網(wǎng)格壓縮法是一種針對(duì)網(wǎng)格拓?fù)溥B接信息和幾何信息的編碼方法，可分為單一位率壓縮法（single-rate compression）[1]和漸進(jìn)壓縮法（progressive compression）[2].在單一位率壓縮中，以網(wǎng)格的拓?fù)溥B接信息和幾何數(shù)據(jù)作為一個(gè)整體進(jìn)行編碼壓縮、傳輸和解碼；漸進(jìn)式壓縮允許多層次細(xì)節(jié)（level of details,LODs）的傳輸和重建.網(wǎng)格編碼壓縮從早期以拓?fù)溥B接信息為主的連接信息驅(qū)動(dòng)壓縮算法到后來(lái)的幾何信息驅(qū)動(dòng)型的漸進(jìn)網(wǎng)格編碼算法[3-4]，需要額外對(duì)掃描點(diǎn)進(jìn)行網(wǎng)格化和拓?fù)溥B接信息的處理，運(yùn)算效率不高.基于點(diǎn)模型的壓縮方法節(jié)省了網(wǎng)格化的開(kāi)銷，處理和渲染更加方便，更易于體現(xiàn)復(fù)雜的拓?fù)浣Y(jié)構(gòu)，加之3D激光掃描設(shè)備的大量應(yīng)用，引起了廣泛關(guān)注.其中，文獻(xiàn)[5]針對(duì)體和面的點(diǎn)云數(shù)據(jù)，利用八叉樹(shù)進(jìn)行空間分割，并結(jié)合行程編碼（run length coding,RLC）和霍夫曼編碼（Huffman coding,HC）對(duì)點(diǎn)云模型進(jìn)行壓縮編碼.文獻(xiàn)[6]給出了一種單一位率的壓縮編碼方法，通過(guò)建立預(yù)測(cè)樹(shù)對(duì)點(diǎn)云模型進(jìn)行預(yù)測(cè)和熵編碼.文獻(xiàn)[7]在點(diǎn)云數(shù)據(jù)八叉樹(shù)結(jié)構(gòu)基礎(chǔ)上引入雙緩存處理技術(shù)，實(shí)現(xiàn)對(duì)點(diǎn)云的差值編碼，獲得了較高的壓縮率，但基于空間點(diǎn)的壓縮方法若要保留其顯著的特征點(diǎn)，仍需要對(duì)點(diǎn)的法矢量進(jìn)行隱式求解，算法依然復(fù)雜.Daribo等[8-11]提出了一種曲線基的點(diǎn)云壓縮方法，將靜態(tài)網(wǎng)格投影到物體的表面，通過(guò)相機(jī)捕獲投影網(wǎng)格點(diǎn)信息，將空間點(diǎn)以曲線基的方式進(jìn)行表達(dá)，應(yīng)用競(jìng)爭(zhēng)的預(yù)測(cè)編碼機(jī)制對(duì)單位率的點(diǎn)云進(jìn)行壓縮.該方法提高了壓縮中的率失真性能，但通用性不大.

本文根據(jù)激光掃描的特點(diǎn)，借鑒Daribo曲線基的思路將掃描點(diǎn)換為二維矩陣，創(chuàng)新性地利用壓縮感知原理直接對(duì)點(diǎn)云的坐標(biāo)位置信息進(jìn)行壓縮,既不必考慮點(diǎn)云的拓?fù)潢P(guān)系，又不必重建隱式特征信息，從而實(shí)現(xiàn)了對(duì)點(diǎn)云的保真壓縮.

1 基于壓縮感知（CS）的點(diǎn)云數(shù)據(jù)壓縮

1.1 壓縮感知的原理介紹

壓縮感知理論認(rèn)為：信號(hào)只要在某一正交空間具有稀疏性，就能以較低的頻率采樣，而以較高概率重構(gòu)該信號(hào)[12-14].運(yùn)用壓縮感知原理對(duì)信號(hào)進(jìn)行壓縮采樣的前提是信號(hào)必須稀疏，但并非所有的信號(hào)都滿足該條件，不過(guò)大多信號(hào)經(jīng)過(guò)適當(dāng)?shù)淖儞Q是可以滿足條件的.因此，壓縮感知的第一步工作就是對(duì)原始信號(hào)進(jìn)行某種正交稀疏變換.首先，考慮一維空間RN中一個(gè)有限長(zhǎng)度數(shù)字信號(hào)x，可看成一個(gè)N×1的列向量，于是用一組正交基進(jìn)行稀疏變換

式中，ψ為一個(gè)正交基矩陣[ψ1,ψ2,···,ψN].x為空間域，s表示在基ψ下的變換域表示，是N×1的列向量，顯然s和x都能代表信號(hào).xψ如果變換后的s含有k個(gè)較大的非零值，則可以說(shuō)信號(hào)x是k稀疏的.

其次，用一個(gè)觀測(cè)矩陣φ對(duì)信號(hào)x進(jìn)行觀測(cè)，得到觀測(cè)值y

式中，φ為一個(gè)M×N矩陣.變換后的信號(hào)由RN空間壓縮到RM空間（M＜N），實(shí)現(xiàn)了對(duì)原始信號(hào)的壓縮處理.需要強(qiáng)調(diào)的是：采樣和壓縮在這個(gè)過(guò)程中是同時(shí)進(jìn)行的.

觀測(cè)矩陣的選擇也非常重要，應(yīng)滿足約束等容條件（restricted isometric property,RIP）

約束等容的等價(jià)描述如下：測(cè)量矩陣φ和稀疏正交化矩陣ψ的基不相關(guān)，能夠滿足約束等容性條件的測(cè)量矩陣可以是獨(dú)立同分布的高斯矩陣，也可以是伯努利分布的矩陣，且滿足

式中，c是一個(gè)較小的常量[15].

1.2 信號(hào)的重構(gòu)

若已知測(cè)量矢量y、觀測(cè)矩陣φ、稀疏化正交矩陣ψ，希望重構(gòu)出k稀疏信號(hào)s，即求解如下欠定方程組：

該方程可轉(zhuǎn)化為l0范數(shù)問(wèn)題

式（6）的求解是l0范數(shù)下的優(yōu)化問(wèn)題，是一個(gè)NP問(wèn)題，可以轉(zhuǎn)換為l1范數(shù)的求解問(wèn)題.

該優(yōu)化算法是一個(gè)凸優(yōu)化的過(guò)程，可以轉(zhuǎn)化為線性規(guī)劃的問(wèn)題，可以通過(guò)基追蹤（basic pursuit,BP）算法、匹配追蹤（match pursuit,MP）等算法來(lái)完成，本文采用正交匹配追蹤（orthogonal matching pursuit,OMP）算法.

2 激光掃描點(diǎn)云的降維處理

圖1 點(diǎn)云降維框圖Figure 1 Dimensions reduction principle of the point cloud data

分析FARO激光掃描儀獲取點(diǎn)云數(shù)據(jù)發(fā)現(xiàn)：除了在掃描跳變點(diǎn)外，連續(xù)相鄰數(shù)據(jù)間的差異很小，有很強(qiáng)的相關(guān)性，而這種相關(guān)性正是進(jìn)行數(shù)據(jù)壓縮的前提.如果能根據(jù)數(shù)據(jù)特點(diǎn)對(duì)點(diǎn)云進(jìn)行重新組合，相關(guān)性就會(huì)更強(qiáng).由此，可對(duì)激光掃描儀所獲取的點(diǎn)云進(jìn)行降維處理，處理的過(guò)程見(jiàn)圖1.令原始掃描點(diǎn)云的點(diǎn)集為P={P1,P2,···,Pi,Pi+1,···,Pn}，其中Pi(x,y,z,r,g,b)為其中一個(gè)掃描點(diǎn)，包含浮點(diǎn)型的位置坐標(biāo)分量（x,y,z）和整型的色彩信息分量（r,g,b）.為統(tǒng)一數(shù)據(jù)類型，可對(duì)點(diǎn)集中的數(shù)據(jù)進(jìn)行截短處理，截掉色彩分量，只保留浮點(diǎn)型的位置坐標(biāo)信息分量（x,y,z）.這樣，截短后一個(gè)點(diǎn)可表示為Pi(x,y,z).需要注意的是在截短的數(shù)據(jù)中，同一點(diǎn)的3個(gè)不同分量相關(guān)性較小.考慮相鄰點(diǎn)的相同坐標(biāo)分量有更強(qiáng)相關(guān)性，為獲得更好的壓縮效果，可對(duì)截短的點(diǎn)云數(shù)據(jù)重排序，抽取點(diǎn)云中相同坐標(biāo)分量x（y或z）排為一列，然后截取連續(xù)2n個(gè)坐標(biāo)分量數(shù)據(jù)作為一個(gè)列矢量，最后形成大小2nN的矩陣X來(lái)存貯新生成的點(diǎn)云數(shù)據(jù)，本文取2n=128.

3 激光掃描點(diǎn)云的CS壓縮

降維處理后點(diǎn)云矩陣X（128×N）的相關(guān)性加強(qiáng)，但進(jìn)行CS壓縮還需稀疏化處理.本文選用簡(jiǎn)單且行之有效的Harr小波基，并用MATLAB語(yǔ)言編寫(xiě)生成三級(jí)變換矩陣T（128×128），實(shí)現(xiàn)對(duì)矩陣X列向量的稀疏處理，變換后的矩陣可表示為

變換后矩陣列向量行號(hào)較小的部分相當(dāng)于低頻分量，行號(hào)較大的部分相當(dāng)于高頻分量.

接著進(jìn)行采樣，矩陣X經(jīng)過(guò)三級(jí)Harr變換后較大系數(shù)值（即k值）約在16～20之間，根據(jù)式（4）將M取3～4倍的k值（本文中取M為80），生成該矩陣的MATLAB語(yǔ)句為：“φ=rand(M,N)”，M和N即為采樣矩陣的行和列.該壓縮采樣可用公式表示為

式中，Y為壓縮后的數(shù)據(jù).

4 壓縮點(diǎn)云數(shù)據(jù)的重建與統(tǒng)計(jì)濾波

4.1 壓縮點(diǎn)云數(shù)據(jù)的重建

若已知采樣值Y和測(cè)量矩陣φ，數(shù)據(jù)恢復(fù)采用OMP算法.恢復(fù)X1中第n列的步驟如下：

步驟1 令殘差向量r=Yn，Yn代表測(cè)量值中的第n列；

步驟2 用測(cè)量矩陣的各列與r進(jìn)行內(nèi)積運(yùn)算〈r,φi〉，i=1,···,128，找出內(nèi)積最大的一列φl(shuí),l為列號(hào)；

步驟3 用最小二乘法獲得X1中第n列第l位置上元素值：

步驟4 修正殘差值

步驟5 在觀測(cè)矩陣φ中剔除掉φl(shuí)列，返回步驟2，確定X1第n列下一個(gè)元素值及其位置，直到待定出全部的k個(gè)元素，該列的其他位置置0.

這樣就完成了對(duì)X1矩陣一列數(shù)據(jù)的恢復(fù)，同理可恢復(fù)其他各列.整個(gè)X1矩陣恢復(fù)后，可以通過(guò)Harr逆變換得到點(diǎn)云矩陣X.最后通過(guò)矩陣的逆排序進(jìn)行原始點(diǎn)云的恢復(fù).

4.2 點(diǎn)云重建統(tǒng)計(jì)濾波

在CS進(jìn)行稀疏化處理過(guò)程中，如果相鄰點(diǎn)位置坐標(biāo)跳變，或在不同掃描線過(guò)渡中存在較大坐標(biāo)差異，都會(huì)引起重建后點(diǎn)云數(shù)據(jù)的稀疏外點(diǎn)噪聲干擾，因此有必要對(duì)重建后的點(diǎn)云數(shù)據(jù)進(jìn)行統(tǒng)計(jì)濾波.統(tǒng)計(jì)濾波的原理如下：首先計(jì)算點(diǎn)云中每個(gè)點(diǎn)的k鄰域均值距離μ和標(biāo)準(zhǔn)偏差σ，并認(rèn)為其全局統(tǒng)計(jì)結(jié)果服從Gauss分布.濾波時(shí)將那些落入μ±λσ之外的點(diǎn)作為噪聲外點(diǎn)予以濾除，其中λ為倍乘因子.該統(tǒng)計(jì)濾波用c語(yǔ)言編程實(shí)現(xiàn).

5 實(shí)驗(yàn)結(jié)果及分析

實(shí)驗(yàn)中數(shù)據(jù)源為Faro激光掃描儀所獲得的點(diǎn)云，采用MATLAB和VC混合編程的方法，PC配置為Intel?CoreTMi5-2400 CPU＠3.10 GHz,4 G內(nèi)存，Win7操作系統(tǒng).

首先將一個(gè)模型猴子選作壓縮處理對(duì)象，預(yù)處理后獲得的原始點(diǎn)云數(shù)據(jù)見(jiàn)圖2（a），包含700 672個(gè)點(diǎn)，壓縮復(fù)原后的點(diǎn)云見(jiàn)圖2（b），點(diǎn)數(shù)為437 920.在壓縮算法中，高斯矩陣采樣的比率為80：128，采樣壓縮比為1.6.為獲得壓縮的客觀評(píng)價(jià)，本文求出壓縮前后各對(duì)應(yīng)點(diǎn)的距離差值平方，求和后再求平均.計(jì)算壓縮前后的誤差公式為

式中，Xi為原始點(diǎn)云數(shù)據(jù)中一個(gè)點(diǎn)的列向量，Ri為重構(gòu)后點(diǎn)云數(shù)據(jù)中一個(gè)點(diǎn)的列向量，N為點(diǎn)云數(shù)據(jù)點(diǎn)的總數(shù)量.運(yùn)算后得到的誤差為3.4398×10-7m2，可見(jiàn)誤差較小.

圖2 點(diǎn)云壓縮比較Figure 2 Comparison of point cloud data before and after compression

為了驗(yàn)證算法的通用性，換用大場(chǎng)景數(shù)據(jù)和不同類型的點(diǎn)云進(jìn)行測(cè)試，壓縮所對(duì)應(yīng)的數(shù)據(jù)見(jiàn)表1，實(shí)驗(yàn)結(jié)果見(jiàn)圖3～5.

表1 不同點(diǎn)云壓縮對(duì)比Table 1 Compression comparison of different point cloud data

從表1中可以看出，本算法壓縮比都在2倍以上，壓縮后較好地保持了原始點(diǎn)云數(shù)據(jù)的特征信息，且整體的誤差較小，如圖2、圖4、圖5所示.從圖2和5中還可以得到：本算法不僅能較好地實(shí)現(xiàn)對(duì)較平坦的點(diǎn)云的壓縮（見(jiàn)圖2），而且能夠?qū)崿F(xiàn)對(duì)曲率變化較大的點(diǎn)云數(shù)據(jù)的壓縮（見(jiàn)圖5），說(shuō)明算法具有一定的通用性.不足之處是：雖然對(duì)密度較大（對(duì)應(yīng)的文件也較大）的點(diǎn)云的壓縮效果較好（見(jiàn)圖4和5）；密度較小點(diǎn)云則存在部分細(xì)節(jié)的丟失，見(jiàn)圖3復(fù)原圖方框部分，且誤差也較大（0.039 6 m2），這主要是由于激光掃描獲得的原始點(diǎn)云密度較小的緣故.

本文中的壓縮率可以通過(guò)調(diào)節(jié)采樣矩陣的大小進(jìn)行調(diào)節(jié)，但一定要滿足約束等容條件，否則會(huì)使復(fù)原后的點(diǎn)云質(zhì)量下降，文中的采樣矩陣為80×128.

圖3 上大圖書(shū)館點(diǎn)云壓縮比較Figure 3 Comparison of SHU library point cloud data before and after compression

圖4 書(shū)本點(diǎn)云數(shù)據(jù)壓縮比較Figure 4 Comparison of a book point cloud data before and after compression

圖5 紙杯點(diǎn)云數(shù)據(jù)壓縮比較Figure 5 Comparison of a paper bottle point cloud data before and after compression

當(dāng)然，為了去除壓縮恢復(fù)后點(diǎn)云外部噪聲點(diǎn)，可用前面提及的統(tǒng)計(jì)濾波對(duì)猴子點(diǎn)云數(shù)據(jù)進(jìn)行處理，其中濾波均值k=30，標(biāo)準(zhǔn)偏差σ=1，對(duì)比結(jié)果見(jiàn)圖6，左半部分濾波前的點(diǎn)云數(shù)據(jù)為700 762個(gè)，右半部分濾波后的點(diǎn)云數(shù)據(jù)為619 292，濾除的點(diǎn)數(shù)為81 470.由圖6中（a）和（b）所標(biāo)記的綠色方框顯然可以看到，濾波前后噪聲外點(diǎn)得到明顯的控制，從而使恢復(fù)后的點(diǎn)云邊緣部分更加清晰.

圖6 統(tǒng)計(jì)濾波前后對(duì)比圖Figure 6 Comparison of point cloud data after and before statistical f iltering

6 結(jié)語(yǔ)

本文結(jié)合激光掃描點(diǎn)云數(shù)據(jù)的特點(diǎn)，創(chuàng)新性地將CS理論用于點(diǎn)云壓縮，該算法不需要點(diǎn)云數(shù)據(jù)的拓?fù)潢P(guān)系，只需根據(jù)掃描數(shù)據(jù)的先后順序，直接對(duì)數(shù)據(jù)坐標(biāo)點(diǎn)進(jìn)行采樣壓縮.通過(guò)實(shí)驗(yàn)對(duì)比，可以發(fā)現(xiàn)壓縮效果良好.但也存在一些需要改進(jìn)的地方：隨著壓縮比的增大，噪聲點(diǎn)明顯增多，于是可以考慮采用更優(yōu)的采樣矩陣及恢復(fù)算法加以克服；壓縮算法中未考慮彩色信息，這也是下一步需要解決的問(wèn)題.

[1]ALLIEZ P,DESBRUNM.Valence-driven connectivity encoding for 3D meshes[C]//European Association for Computer Graphics,2001:480-489.

[2]ALLIEZ P,DESBRUN M.Progressive encoding for lossless transmission of triangle meshes[C]//Association for Computing Machinery's Special Interest Group on Computer Graphics and Interactive Techniques,2001:198-205.

[3]PENG J, KUO C C J.Progressive geometry encoder using octree-based space partitioning[C]//Proceedings of the 2004 IEEE International Conference on Multimedia and Exposition,2004:1-4.

[4]PENG J,KUO C C J.Geometry-guided progressive lossless 3D mesh coding with octree(OT)decomposition[C]//Association for Computing Machinery's Special Interest Group on Computer Graphics and Interactive Techniques,2005:609-616.

[5]SIDDIQUI R A,CELASUN I,OCTREE U B.Based compression of volumetric and surface 3D point cloud data[C]//Proceedings of 13th Intl Conference on Virtual Systems and Multimedia,Brisbane,Australia,2007:278-282.

[6]GUMHOLD S,KAMI Z,ISENBURG M,SEIDEL H P.Predictive point-cloud compression[C]//Proceedings of the sixth Israel-Korea Binational Conference,2005:125-129.

[7]KAMMERL J,BLODOw N,RUSU R B,GEDIKLI S,BEETZ M,STEINBACH E.Real-time compression of point cloud streams[C]//Robotics and Automation(ICRA),Saint Pual,2012:778-785.

[8]DARIBO I,FURUKAwA R,SAGAwA R,KAwASAKI H,HIURA S,ASADA N.Point cloud compression for grid-pattern-based 3D Scanning System[C]//Visual Communications and Image Processing(VCIP),2011:1-4.

[9]DARIBO I,FURUKAwA R,SAGAwA R,KAwASAKI H.Dynamic compression of curve-based point cloud[C]//SIVT'11 Proceedings of the 5th Pacif ic Rim conference on Advances in Image and Video Technology-Volume Part II,2011:323-334.

[10]DARIBO I,FURUKAwA R,SAGAwA R,KAwASA K H,HIURAS,ASADAN.Efficient rate-distortion compression of dynamic point cloud for grid-pattern-based 3D scanning systems[C]//3D Research,2012,3(1):1-9.

[11]DARIBOI,FURUKAwA R,SAGAwA R,KAwASAKIH,HIURA S,ASADA N.Curve-based representation of point cloud for efficient compression[C]//The 14th Meeting on Image Recognition and Understanding,2011:1385-1390.

[12]DONOHO D．Compressed sensing[J]．IEEE Transactions on Information Theory,2006,52(4):1289-1306．

[13]CANDE S E,EMAMNUEL J C.Compressive sampling[C]//Proceedings of International Congress of Mathematicians Switzerland:European Mathematical Society Publishing House,2006:1433-1452．

[14]CANDE'S E,ROMBERG J,TAO T．Robust uncertainty principles:exact signal reconstruction from highly incomplete frequency information[J]．IEEE Transactions on Information Theory,2006,52(2):489-509．

[15]BARANIUK R．A lecture on compressive sensing[J].Signal Processing Magazine,2007,24(4):118-121．