移動最小二乘重采樣在三維重建中的應用

康傳利,時滿星,程 耀,張臨煒,顧峻峰,陳 洋

(桂林理工大學 a.廣西空間信息與測繪重點實驗室;b.測繪地理信息學院,廣西 桂林 541006)

隨著三維激光掃描儀的普及和計算機硬件水平的提高,地面三維激光掃描儀被廣泛應用于數字化城市[1]、古建筑文物保護[2]、虛擬現實[3]等領域。為了真實重建具有不規則特征的實體,通常直接構建其表面三角網,這種曲面重建的方法主要分為兩類,即顯式曲面重建和隱式曲面重建。顯式曲面重建主要是基于Delaunay三角剖分和voronoi圖的計算幾何方法,如Crust算法[4]、tight Cocone算法[5]、3D alpha-shape算法都是在Delaunay三角剖分的基礎上發展而來的,若點云數據非均勻且含有噪聲,使用該種重建方法會丟失曲面重建的精度，并導致幾何模型變形。三維激光掃描儀采集的點云數據按點與點之間的排列方式可以分為:線式點云數據、陣列式點云數據、部分有序點云數據和完全散亂點云數據。地面三維激光掃描儀的點云數據屬于陣列式結構,具有明顯的空間拓撲關系，但為了全方位地獲得實體表面的空間數據信息,需進行多站點間的配準。因每兩個測站至少要有20%以上的重疊區域,經配準后的點云數據拓撲關系將會丟失,數據冗余度變大;且受配準精度的影響,任意兩站點的數據不能完全融合,小尺度噪聲將會增加,不利于曲面重建。Aurich等[6]研究了基于線性濾波器的高斯濾波算法,缺點是無法適用于不存在拓撲關系的散亂點云數據。馬先明等[7]對無人機影像生成的點云數據利用雙邊濾波法進行了去噪實驗和分析,通過建立點云鄰域對小型噪聲點進行了糾正,缺點是若噪聲點影響到鄰域特征的估計,則會減弱該噪聲算法的穩定性。宋陽等[8]針對點云數據中體外孤點不易辨別難以去噪的問題,提出了一種改進的C均值算法,有效去除了大尺度噪聲,但該方法對點云數據中小尺度噪聲的去噪效果不佳。

針對上述散亂點云數據噪聲難以用統計分析的方法消除以及數據存在冗余從而影響曲面重建質量和效率的問題,本文提出了一種基于移動最小二乘(MLS)平滑去噪重采樣的方法,既獲得了隱式曲面高抗噪性的光順效果,又保留了實體表面的細節特征。

1 MLS算法描述

Lancaster等以最小二乘法為理論基礎最早提出了移動最小二乘法,之后，相關學者對其進行了改進,提高了MLS的自適應能力,同時也繼承了最小二乘法高精度的數值計算能力。本文MLS算法的計算流程如圖1所示：首先建立空間數據點Pi的鄰域,有基于規則采樣的歐氏鄰域和基于不規則采樣的k近鄰鄰域;其次構建以Pi鄰域為局部區域的擬合函數,擬合函數中基函數的階數不同，其擬合精度也不同;之后確定權函數,擬合函數能夠很好地繼承權函數的連續性,如果權函數是C1階連續的,則擬合函數也是C1階連續的[9],不同的權函數平滑效果也不同;通過設置體素化網格(VoxelGrid)的大小,將點云數據所在的空間進行體素柵格化處理,用體素網格內所有點的重心近似代替網格中所有的數據點,從而簡化冗余數據點,最后得到重采樣后的點云。目前大多數平滑去噪算法在進行去噪時都沒有考慮到點模型的優化,點云平滑前后數量未發生改變,本文算法將體素化網格下采樣模型引入到MLS平滑去噪中,與單獨使用MLS重采樣算法相比,在點云平滑去噪的同時實現三維點模型的優化,能夠有效去除小尺度噪聲并解決點云冗余問題。

圖1 MLS算法流程圖Fig.1 Algorithm flow chart of MLS

2 MLS算法基本原理

2.1 擬合函數的建立

在擬合區域的一個局部子域U上,其擬合函數f(x)表示為

(1)

式中:P∈U;α(x)=[α1(x),α2(x),…,αm(x)]T為待求系數,它是一個關于坐標x的函數;P(x)=[P1(x),P2(x),…,Pm(x)]T稱為基函數,是一個k階完備多項式;m是基函數的項數。

對于二維區域,通常基函數的形式為

線性基P=[1,x,y]T,m=3;

二次基P=[1,x,y,x2,xy,y2]T,m=6。

在實際運用中,二次基函數使用較多,所以式(1)也可以描述為

f(x)=α0(x)+α1(x)x+α2(x)y+α3(x)x2+

α4(x)xy+α5(x)y2。

(2)

(3)

式中:n為受影響區域內節點的數目;f(x)為擬合函數;yi=y(xi)是x=xi處的節點值;w(x-xi)是節點xi的權函數。

α=(BWBT)-1BWy,

(4)

W為n×n的對角矩陣;y=[y(x1),y(x2),…,y(xn)]T。

由式(4)能夠求出式(1)中的系數, 再將待定點代入式(1), 可得到該點的擬合函數值。

2.2 確定權函數

移動最小二乘法、最小二乘和加權最小二乘的主要區別是在權函數的運用上,如引入緊支(compact support)的概念。認為點x處的值僅受x子域中相鄰節點的影響,該子域稱為點x的影響區域,影響區域一般是半徑為h的圓形,權函數在x的影響區域內不為0,在影響區域外全為0，并且權函數w(x-xi)應該是非負的, 權值隨著‖x-xi‖2的增加而單調遞減。常用的權函數主要有:

①三次樣條權函數

(5)

②高斯權函數

(6)

式中:di(x)表示第i個原始采樣點Pi與重采樣點P之間的距離；u表示影響區域中的特征與重采樣間的影響因子。

在空間點云數據重采樣的過程中,重采樣點的位置相對于原始采樣點會發生變化,因此式(4)中權因子w(x-xi)應隨著待求點位置的不同而進行移動,本文選用高斯權函數,能夠對式(4)進行重新計算,得到擬合函數新的系數矩陣

α(x)=(BW(x)BT)-1BW(x)y,

(7)

式中:W(x)是由式(6)計算得到的對角矩陣。

2.3 體素化(VoxelGrid)網格下采樣模型

體素化網格下采樣的基本原理是[10-11]:在點云所在坐標系下建立該輸入點云的最小包圍盒,包圍盒是與各坐標平面平行且完全包含點云數據的最小六面體,根據用戶輸入采樣距離的大小,將六面體等間距均勻分解成m×n×l個小網格立方體,所有的點云數據被再次劃分到小網格立方體中,這些小網格立方體可看作是子包圍盒,刪除子包圍盒之外的小網格立方體,再分別計算子包圍盒中所有點的重心,以該重心來近似代替這些點。相對于其他點云簡化方法,該種方法簡單易實現,特別是對于散亂點云不需要再次建立點云間的拓撲關系,減小了計算的復雜度。該采樣模型的實現步驟如下:

①最小包圍盒的確定。找出輸入點云坐標系下各個坐標軸方向、點云坐標的最大值和最小值,分別記為Xmax、Xmin、Ymax、Ymin、Zmax、Zmin。

②根據用戶輸入的采樣距離L,確定最小包圍盒在各坐標軸方向被均勻分割出的立方體的個數,即m、n、l的值。

(8)

③建立子包圍盒的索引。點云數據中的任一點Pi(xi,yi,zi)都被點云數據最小包圍盒的子包圍盒所包圍, 所以能夠根據這個點的空間位置確定其子包圍盒的空間索引(mindex,nindex,lindex)。

(9)

④子包圍盒內近似代表點的確定。建立好子包圍盒的索引后,根據該子包圍盒內所有點的個數,求出它們的重心,并用這個重心來代表這些點。若某個子包圍盒內有n個點,則代表點的坐標為

(10)

3 實驗結果與分析

3.1 實驗環境

在VS2010環境下采用C++語言和PCL實現了本文的算法,即完成空間點云數據的平滑去噪和重采樣。實驗的硬件配置為Intel(R) Xeon(R) CPU,2.40 GHz,4.00 GB內存,Windows 7操作系統,對采樣后的點分別在Geomagic Studio 2013和Matlab 2014中進行了顯式曲面重建。

3.2 實驗數據

主要測試數據來源于地面三維激光掃描儀Riegl vz1000對桂林理工大學雁山校區東門前景觀石的掃描, 為了全方位的獲取景觀石表面的空間數據信息, 此次掃描共架設5個站點, 掃描角分辨率為0.13°, 掃描時長30 min, 在掃描儀自帶的聯機控制軟件Pro riscan進行了多站點的ICP(iterative closest point)配準, 整體配準精度為0.003 2 m。圖2a為景觀石前方表面單個站點的掃描數據,點云排列整齊、凹凸有致;圖2b為配準后的掃描數據,點云的空間關系被完全打亂,且受配準精度的影響,具有相同位置的點云并不能完全重合,點云的冗余度較大。

圖2 配準前后的點云Fig.2 Point cloud before and after registration

3.3 基于Delaunay三角剖分的Crust三維重建

設置體素化網格小立方體的邊長為0、2.0、2.5、3.0 cm,分別與圖3和表1中下采樣1(僅MLS重采樣)、下采樣2、下采樣3、下采樣4相對應,原始點云是指景觀石的原始掃描數據。

點云名稱點云個數三角形個數點數簡化率/%Delaunay三角化耗時/s曲面重建總耗時/s原始點云455092905936085.84185.95下采樣1455092910172095.88194.44下采樣211033822066875.822.2846.99下采樣37403314806083.715.0531.61下采樣45301310602088.410.6522.43

從圖3b—d可以看出,經體素化網格下采樣后的點云在空間結構上具有一定的拓撲關系,且下采樣點云密度在具有表面特征的區域較大,如景觀石表面刻字部分的點個數比同等范圍平緩區域要多,即曲率變化大的區域點云分布比較集中,避免了細節特征信息的丟失。從圖4a和表1可以看出,存在冗余數據的點云會嚴重影響曲面重建的時間開銷，且三維重建效果不佳,正常分辨率下因重建網格過于稠密無法顯示出三維立體效果;從圖4d中可以看出景觀石表面“字跡”較為模糊,主要因采樣距離設置過大,劃分進入子包圍盒的點數較多,相應刪除了較多的點,從而造成局部特征不明顯。從表1可以看出,經MLS重采樣的點云,其中滿足Delaunay三角剖分空外接圓和最小內角最大準則條件的點數增加,故下采樣1重建三角形的個數要高于原始點云,曲面重建時間有所增加，通過對Delaunay三角剖分和曲面重建總耗時的比較,采樣間距越大點云簡化率越高，其表面重建時間越短。

對圖3、 圖4和表1在點云分布情況、 曲面重建立體化效果和曲面重建時間開銷等進行綜合性分析可以得出, 體素化網格下采樣3(2.5 cm×2.5 cm×

圖4 Crust三維重建部分截圖Fig.4 Partial screenshot of the Crust 3D reconstruction

2.5 cm)能夠在保證重建質量下同時提高重建效率, 點數簡化率為83.7%, 重建時間縮短了83%。 因此，針對點云數據冗余的情況, 能夠根據點云的平均密度和局部細節特征選取合適的采樣間距進行簡化處理,提高曲面重建效率。

3.4 基于Geomagic Studio的三維重建

為了驗證本文算法的適用性及在實際建模中的魯棒性,對MLS重采樣的點云在逆向工程軟件Geomagic Studio中進行了真三維重建,并截取了景觀石表面帶有明顯細節特征(“理”的字跡)的部分點云作進一步驗證和分析(圖2b)。圖5a—c是原始點云分別在無平滑去噪、平滑去噪強度低、平滑去噪強度高的條件下封裝成的三角網格表面；圖5d為無去噪三角網格圖5a表面經Laplace網格平滑去噪后的效果圖；圖5e是經MLS重采樣后的點云在無平滑去噪條件下封裝成的網格表面。

在無去噪的情況下對配準后的點云數據進行網格化重建,其重建表面布滿了鱗片狀的褶皺,字跡處也出現了漏洞(圖5a),說明因配準生成的小尺度噪聲會嚴重影響曲面重建的質量;比較圖5b、c表面的重建效果,隨著去噪強度的提高,重建表面也愈趨近于光滑,但圖5c表面的字跡因過于光滑已經變得模糊,造成細節特征信息的丟失;比較圖5d和e,因Laplace是采用多次迭代的方法將數據點向該點鄰域中心移動基于網格的平滑去噪算法,故無法修復圖5a中由噪聲形成的漏洞；但從圖5e中可以看出，經本文算法重采樣后的點云其表面重建效果較好,不僅整體表面較為平滑,同時字跡特征也能清晰地顯現出來且無漏洞,視覺效果優于Laplace平滑去噪算法。

基于本文算法,將景觀石重采樣的點云數據導入到Geomagic Studio中直接封裝成網格模型,如圖6a所示,可以看出表面的整體光滑度高,局部細節特征清晰;為了給模型增加真實感的視覺效果,對圖6a進行紋理貼圖,生成景觀石的真三維實體模型如圖6b所示。

圖6 景觀石實體模型Fig.6 Solid model of landscape stone

4 結 論

點云數據的平滑去噪和簡化,對模型三維立體重建的質量和快束存儲、處理、顯示等具有重要的意義。將移動最小二乘重采樣運用到地面三維激光掃描儀對大型構筑物的三維重建上,能夠很好地解決因配準誤差不可避免而造成數據冗余和小尺度噪聲增加的問題。經實例分析表明:

(1)選擇合適的體素化網格小立方體的邊長,能夠有效簡化海量點云數據,且重采樣后的點云在曲率變化大的地方分布多，平緩區域分布少,在提高曲面重建效率的同時也保證了曲面重建的質量。

(2)將移動最小二乘重采樣后的點云直接在逆向工程軟件中封裝成網格模型,其表面細節特征清晰、光滑度好,跳過了后期打磨、松弛和去除特征等修理環節,再對網格模型進行紋理貼圖生成真三維實體模型,逼真程度高、視覺效果好,優于Laplace網格平滑去噪算法。

圖5 平滑去噪效果Fig.5 Results of smoothing and denoising

因此,移動最小二乘重采樣在地面三維激光掃描儀三維重建中具有很高的實用價值。