一種改進的表面重建算法及其并行化研究

王元+曹靜+婁澤坤

摘 要： 對于泊松三維表面重建中估算指示函數(shù)出現(xiàn)誤差的情況，引入屏蔽因子并使用線性插值的方法糾正誤差。該方法在等值面提取時使用的索引表查找法存在二義性的問題，需使用漸近線方法來解決。改進算法需要較大的存儲空間來存儲頂點信息和法向量以及較長的時間來運算，因而引入了CUDA架構(gòu)的并行化計算來提高運行效率。

關(guān)鍵詞： 多視圖三維重建； 泊松表面重建； 二義性； CUDA架構(gòu)

中圖分類號：TP391 文獻標(biāo)志碼：A 文章編號：1006-8228（2017）05-17-03

An improved 3D surface reconstruction method and it's parallelization

Wang Yuan， Cao Jing， Lou Zekun

（School of Computer and Information Engineering， Henan University， Kaifeng， Henan 475000， China）

Abstract： In this thesis the shielding factor is introduced to solve the problem of errors of indication function to estimate the indicator function in the Poisson surface reconstruction algorithm. In the part of isosurface extraction in surface reconstruction， the triangulation process is performed by searching the index table. But this method will lead to ambiguity， so an asymptotic method is used to avoid it. This surface reconstruction algorithm will improve the holes and depressions problems of complex objects. At the same time， it leads to larger storage space to store vertex and normal vectors， and a longer time to process data. Therefore， CUDA architecture and parallelization are introduced to solve the problem of efficiency.

Key words： multi view 3D reconstruction； Poisson surface reconstruction； ambiguity； CUDA architecture

0 引言

2006年，Michael Kazhdan等人提出了泊松表面重建算法[1]，但泊松表面重建本身存在諸多缺陷。如該算法在矯正指示函數(shù)的時候使用的是一次性全局的偏移[2]。若在實際的實現(xiàn)過程中出現(xiàn)計算指示函數(shù)不可避免的誤差，此時再用一次性全局偏移的方法估算指示函數(shù)則會出現(xiàn)重建結(jié)果誤差較大的情況[3]；另外等值面提取部分用移動立方體算法，該算法通過索引表查找各個體素內(nèi)三角面片的連接模式確定等值面的分布方式，但是索引表算法會產(chǎn)生連接二義性問題。

鑒于第一個問題，最主要的突破口是找到一個合適的全局偏移量來解指示函數(shù)的近似解，本文嘗試使用顯性插值方法對樣本點進行插值計算。對于存在連接二義性的情況，改進方法是漸近線算法。同時對改進之后算法的效率提升問題引入支持計算統(tǒng)一設(shè)備架構(gòu)（CUDA），采用GPU并行計算，其存儲器帶寬和并行處理能力有CPU處理無法逾越的優(yōu)勢。

1 引入屏蔽因子和漸近線算法的表面重建

1.1 引入屏蔽因子

在泊松方程的求解過程中[5]，對于目的有向點集求解尺度指示函數(shù)，使尺度函數(shù)的梯度與有向點集形成最佳逼近即尺度函數(shù)最小化問題

在此，本文給定帶權(quán)重的點集，為尺度函數(shù)最小化增加一些約束重新計算樣本點函數(shù)，見公式⑵：

其中Area（S）是要重建的表面區(qū)域，α是一個權(quán)值，用來衡量擬合梯度和擬合值，也可叫做屏蔽因子，w（p）是樣本點權(quán)重，在該算法實現(xiàn)中，根據(jù)局部采樣密度估計得到表面區(qū)域，每個樣本點的權(quán)重w（p）設(shè)置為1。上述公式簡化為公式⑶：

表示單位體素在函數(shù)空間上的表示形式，這個值是由采樣點的函數(shù)加權(quán)和得到的，見公式⑷：

公式⑷將梯度約束和離散點約束值整合進空間域，隨即指示函數(shù)的方程最小化為或者。

此時，若指示函數(shù)滿足，則此等式就是屏蔽泊松方程。

通過引入屏蔽因子解屏蔽泊松方程的方法，使用顯性插值計算，解決一次性全局估算指示函數(shù)誤差較大的情況。

1.2 使用漸進性方法解決二義性

移動立方體算法提供了以查找索引表的方式確定體素內(nèi)等值面的連接方式，但這種方法也過于依賴直觀的構(gòu)造，忽略了構(gòu)造三角面片時可能存在的二義性問題。如圖1，在立方體體素的某一等值面上若出現(xiàn)對角線的兩端同時為+，另一對角線兩端同時為-的情況時，不同的連接方式會構(gòu)造不同的結(jié)果，如圖1。

對于存在連接二義性的情況改進方法是漸近線算法，通過邊界面與等值面相交的理論得知，二者相交的結(jié)果是一對雙曲線，計算雙曲線的漸近線，通過雙曲線、雙曲線的漸近線與體素邊界面三種體元的位置關(guān)系來判斷是否產(chǎn)生二義性，對于有面二義性的情況，做漸近線處理判別方法，作圖2（h）和（i）：

當(dāng)頂點為+和頂點為-分別出現(xiàn)在一對對角線兩端時出現(xiàn)連接二義性，此時若漸近線交點為+如圖2（h），則漸近線交點和B點D點處于同一區(qū)域，此時等值點M1和M2連接，M3和M4連接；若漸近線交點為-如圖2（i），則漸近線交點和A點C點處于同一區(qū)域，此時等值點M1和M3連接，M2和M4連接。

2 關(guān)于改進表面重建算法的GPU并行優(yōu)化

從CUDA架構(gòu)可知，CUDA并行算法的設(shè)計需要程序的數(shù)據(jù)相互獨立，在實現(xiàn)的過程中需要處理若干相互不相干的子問題[4]，這些子問題被分割放置在每一個kernel函數(shù)中同時執(zhí)行，對于八叉樹結(jié)構(gòu)每個深度層次上的數(shù)據(jù)計算可用CUDA架構(gòu)并行化處理，可根據(jù)串行程序中耗時不同設(shè)計并行算法。

根據(jù)程序耗時分析，設(shè)計數(shù)據(jù)離散化八叉樹生成部分進行GPU并行計算[6-7]。

2.1 包圍盒的并行計算

包圍盒，可作為八叉樹的根節(jié)點，計算各個樣本點的相對位置。在我們的并行設(shè)計算法中，處理并行計算包圍盒用的是CUDPP的方法，CUDPP是一種強大的數(shù)據(jù)并行CUDA庫。

CUDPP的具體調(diào)用過程如下。

⑴ 初始化，調(diào)用cudppCreate實例化一個CUDPP，生成一個CUDPPManager對象，將其地址返回給theCudpp。

⑵ 為主機端和設(shè)備端使用函數(shù)分配存儲空間，對主機端所要輸入?yún)?shù)集賦值。

⑶ 將數(shù)據(jù)從主機端復(fù)制到設(shè)備端指定的存儲空間，使用cudaMemcpy2D（）函數(shù)進行主機端與設(shè)備端之間的二維數(shù)據(jù)傳輸。數(shù)據(jù)傳輸方向是主機端到設(shè)備端，cudaMemcpyHostToDevice（）。

⑷ 生成“CUDPP計劃”。CUDPP計劃是CUDPP方法中特有的一部分內(nèi)容，具體做法是定義一個CUDPPHandle，生成一個cudppPlan將句柄賦值給scanPlan。

⑸ 執(zhí)行原語。

在CUDPP計劃配置參數(shù)完成之后即要開始執(zhí)行這一過程的kernel函數(shù)。最后，用cudaThreadSynchronize（）函數(shù)同步處理。

⑹ 將輸出數(shù)據(jù)復(fù)制到主機端。繼續(xù)由主機端執(zhí)行邏輯控制。依然使用函數(shù)cudaMemcpy2D（）進行拷貝，數(shù)據(jù)傳輸方向是設(shè)備端到主機端cudaMemcpyDeviceToHost（）。

⑺ 所有處理完成之后釋放所分配的存儲空間。

2.2 生成鍵值序列及鍵值排序并行算法設(shè)計

鍵值序列在樹形數(shù)據(jù)結(jié)構(gòu)中表節(jié)點的權(quán)值，也代表了多維空間中數(shù)據(jù)的表達方式，計算鍵值序列的基礎(chǔ)是計算莫頓碼（Morton碼）[8]。

在八叉樹的數(shù)據(jù)結(jié)構(gòu)中給定一個采樣點p，在d層（1<=d<=D）計算p的鍵值序列xyz，以及為鍵值序列排序。分別計算鍵值xd、yd、zd，定義Cdx為d-1層包含p點的節(jié)點Od-1的形心，p.x為采樣點p相對于形心的位置，若p.x位于節(jié)點Od-1的x軸平分線的左側(cè)，x鍵值為0，否則為1；y鍵值和z鍵值也以同樣的方法計算。

如此，深度為D的八叉樹鍵值序列表示為x1y1z1x2y2z2…xDyDzD，意義為從根節(jié)點到當(dāng)前節(jié)點的路徑。但由于輸入點云數(shù)據(jù)是亂序的，由此生成的每個節(jié)點的鍵值序列也有可能是亂序的，還需對鍵值序列進行排序，可以使用CUDPP庫sort方法處理所生成的鍵值序列。

2.3 生成節(jié)點數(shù)組并生成完整八叉樹結(jié)構(gòu)

由于八叉樹結(jié)構(gòu)的每一個節(jié)點都有0個孩子節(jié)點或7個孩子節(jié)點，因此需要對生成的惟一節(jié)點數(shù)組補全剩余的7個節(jié)點數(shù)組信息。首先需要計算其余7個節(jié)點數(shù)組的起始索引地址，使用scan算法并行計算前綴求和的方法在上一個節(jié)點數(shù)組地址基礎(chǔ)上計算鄰接的節(jié)點數(shù)組地址，得到完整的八叉樹節(jié)點數(shù)組地址。

3 仿真結(jié)果

實驗所需輸入點云數(shù)據(jù)是所參與項目中使用的古建筑點云數(shù)據(jù)集jyd.ply。

使用改進過的表面重建算法，生成結(jié)果與原始算法結(jié)果對比截圖如圖3。

圖3左為jyd原始表面重建算法的重建結(jié)果，圖3右為本章改進算法的重建結(jié)果，從圖中可看出屋脊邊緣部分的孔洞問題得到解決。

分別使用這兩組數(shù)據(jù)串行計算和并行計算實驗三次，結(jié)果如表1。

根據(jù)表格中數(shù)據(jù)可知jyd三次的加速比分別為6.05、5.98、5.78，平均為5.93。

4 結(jié)論

通過引入屏蔽因子使用顯性插值方法對樣本點進行插值計算，避免了一次性全局偏移的方法估算指示函數(shù)出現(xiàn)重建結(jié)果誤差較大的情況，通過引入漸近線構(gòu)造等值線正確連接方式的方法，避免了移動立方體算法中三角化過程可能出現(xiàn)的面二義性問題。實驗結(jié)果證明，改進的表面重建算法，可改善模型表面空洞問題。

同時并行算法設(shè)計方案實驗結(jié)果證明GPU并行算法提高了原有表面重建算法的執(zhí)行效率。

參考文獻（References）：

[1] Lin Y， Chen C， Song M， et al. Dual-RBF based surface reconstruction[J]. The Visual Computer，2009.25（5-7）：599-607

[2] Li X， Wan W， Cheng X， et al. An improved poisson surface reconstruction algorithm[C]. International Conference on. IEEE，2010：1134-1138

[3] Yin C， Gang D， Zhi-quan C， et al. An algorithm of CUDA-based Poisson surface reconstruction[C]. Audio Language and Image Processing （ICALIP）， 2010 International Conference on，2010：203-207

[4] Zhou K， Gong M， Huang X， et al. Data-parallel octrees for surface reconstruction[J]. Visualization and Computer Graphics， IEEE Transactions on，2011.17（5）：669-681

[5] 陳建華，趙飛，葛永斌.求解3維泊松方程的一種新方法[J].江西師范大學(xué)學(xué)報（自然科學(xué)版），2013.4：411-415

[6] 游安軍編著.電路數(shù)學(xué)[M].北京電子工業(yè)出版社，2014.

[7] 蔣建剛.薛定諤和泊松方程有限元法求解[D].廣東工業(yè)大學(xué)碩士學(xué)位論文，2014.

[8] 楊宇，闞凌雁，于佳，等.基于激光掃描的人臉三維重建方法[J].紅外與激光工程，2014.12：3946-3950