基于散亂點(diǎn)云的快速體積計(jì)算法

胡曉彤，陶森柏

(天津科技大學(xué)計(jì)算機(jī)科學(xué)與信息工程學(xué)院，天津 300222)

基于散亂點(diǎn)云的快速體積計(jì)算法

胡曉彤，陶森柏

(天津科技大學(xué)計(jì)算機(jī)科學(xué)與信息工程學(xué)院，天津 300222)

三維可視化體積計(jì)算基本上都是先由散亂點(diǎn)云構(gòu)建出表面網(wǎng)格模型，然后基于網(wǎng)格模型計(jì)算體積，存在計(jì)算量大、速度慢的缺點(diǎn)．針對此問題提出一種快速體積計(jì)算法，首先使用改進(jìn)的增量式 Delaunay三角剖分對散亂點(diǎn)云進(jìn)行四面體剖分；然后利用 K近鄰計(jì)算散亂點(diǎn)的擬合曲面和最小生成樹，得到各點(diǎn)的法向量；由各點(diǎn)法向量剔除體外四面體；最后計(jì)算各四面體體積之和從而得到總體積．實(shí)驗(yàn)表明，該算法不僅保證了計(jì)算準(zhǔn)確度，而且較傳統(tǒng)算法大大提高了效率．

散亂點(diǎn)云；四面體剖分；Delaunay三角剖分；法向量；K近鄰

Abstract：Visual volume calculation in 3D space basically is based on mesh model nowadays which is constructed from scattered point cloud. It exposes inferiority like huge calculation and low speed when volume is needed only. According to that，an algorithm of rapid volume calculation was proposed. First，the convex hull of point cloud was subdivided into tetrahedron with improved incremental Delaunay triangulation. Second，fitting quadric surface and MST of points with KNN were calculated to get normal vectors. Third，those tetrahedron in vitro were removed by using normal vectors. Finally，all tetrahedron’s volume were added up to get object’s volume. Experimental results show that this method can improve efficiency greatly compared with common methods and keep precision of volume calculation.

Keywords：scattered point cloud；tetrahedron subdivision；Delaunay triangulation；normal vector；KNN

科學(xué)計(jì)算可視化(visualization in scientific computing)是 20世紀(jì) 80年代后期提出并發(fā)展起來的一個新領(lǐng)域．它指的是運(yùn)用計(jì)算機(jī)圖形及圖像技術(shù)，將科學(xué)計(jì)算過程及結(jié)果的數(shù)據(jù)轉(zhuǎn)化為屏幕上的圖形及圖像并進(jìn)行交互處理的理論、方法和技術(shù)．對于可視化出來的圖形進(jìn)行交互測量在醫(yī)學(xué)、地質(zhì)勘探、天體物理等領(lǐng)域都有重要意義．

要實(shí)現(xiàn)可視化交互測量，先要獲得目標(biāo)物的三維點(diǎn)云，這通?；跈C(jī)器視覺或設(shè)備采集實(shí)現(xiàn)．基于機(jī)器視覺即由多個攝像機(jī)獲得圖像序列，通過圖像匹配，計(jì)算目標(biāo)物的三維坐標(biāo)；基于設(shè)備采集指由三坐標(biāo)儀或激光掃描儀等設(shè)備采集散亂點(diǎn)云．然后進(jìn)行測量，現(xiàn)行方法基本上都是基于不規(guī)則網(wǎng)格模型的計(jì)算，其中三角網(wǎng)格模型使用最為廣泛．任意網(wǎng)格模型體積計(jì)算的常見方法有：坐標(biāo)法[1]，即由封閉空間的表面基元的坐標(biāo)行列式計(jì)算體積；切片法[2]，用一組互相平行的平面剖切模型，由剖面面積和平面間距計(jì)算體積；投影法[3]，計(jì)算三角網(wǎng)格模型的三角面片及其在投影平面上的投影所圍成的凸五面體的帶符號體積，所有帶符號體積的代數(shù)和即為總體積．

這些方法都需要先重建出表面網(wǎng)格模型，然后在此基礎(chǔ)上計(jì)算體積．雖然體積計(jì)算算法本身的精度和效率都較高，但是整個過程大部分時(shí)間耗費(fèi)在了重建表面網(wǎng)格模型上．在單純需要得到目標(biāo)體積的情況下，把時(shí)間浪費(fèi)在重建表面網(wǎng)格模型上是不需要的．因此，本文提出了一種基于不附帶任何拓?fù)湫畔⒌纳y點(diǎn)云的四面體剖分體積計(jì)算法，該算法不需要重建表面網(wǎng)格模型，可直接由點(diǎn)云得到體積，因而提升了體積計(jì)算的效率．

1 四面體剖分體積算法

1.1 基于凸包分割的四面體剖分

四面體剖分本質(zhì)上就是將點(diǎn)集S?R3的凸包分割成若干以 S中的點(diǎn)為頂點(diǎn)的四面體．傳統(tǒng)的四面體剖分算法在浮點(diǎn)數(shù)系統(tǒng)中由于退化現(xiàn)象而不夠健壯，所以本文采用基于增量式 Delaunay三角剖分的改進(jìn)算法．Edelsbrunner和Seidel于1986年證明：任何維度為d的有限點(diǎn)集S?Rd，其 Delaunay三角剖分都可以從的凸包中獲得[4]．

1.1.1 基本概念

(1)凸包：是指包含點(diǎn)集S?R3中任意兩點(diǎn)連成的線段的最小凸多面體．

(2)Delaunay規(guī)則：三維空間中的四面體的外接球內(nèi)不包含其他的散亂點(diǎn)．

(3)退化現(xiàn)象：在浮點(diǎn)數(shù)系統(tǒng)中，因計(jì)算機(jī)數(shù)值精度有限而導(dǎo)致拓?fù)溴e誤．

1.1.2 算法

(1)三維 Delaunay三角剖分是由二維 Delaunay三角剖分?jǐn)U展而來．首先構(gòu)建一個初始四面體，形成初始化四面體網(wǎng)格．

(2)將散亂點(diǎn)插入當(dāng)前四面體網(wǎng)格中，對于輸入點(diǎn)P，使用隨機(jī)行走方法來尋找包含P的四面體．先指定一個四面體T，如果 P位于該四面體內(nèi)，則完成行走．如果不在四面體內(nèi)，則隨機(jī)指定一個三角面E，如果E所在的平面將T和P分割開(即T和P在平面的兩邊)，下一個訪問的四面體就是共享E的鄰近四面體；否則，就按預(yù)定的順序遍歷其他的面，直到找到分割開T和P的面．

(3)找到包含P的四面體，則分割該四面體為4個小四面體．

(4)如果P位于當(dāng)前四面體網(wǎng)格外，則選擇網(wǎng)格的一個可見面(即P在面的一側(cè))，連接P與該三角面的3個頂點(diǎn)構(gòu)成新的四面體加入到四面體網(wǎng)格中，選擇可見面時(shí)，盡量避免使新生成的四面體是狹長的.

(5)重復(fù)步驟(2)—(4)，直到所有散亂點(diǎn)都被插入四面體網(wǎng)格．

(6)驗(yàn)證Delaunay三角剖分的有效性．首先檢查Delaunay三角剖分?jǐn)?shù)據(jù)結(jié)構(gòu)的連貫性，即四面體的鄰接關(guān)系．然后驗(yàn)證各四面體的方向和由 Delaunay三角剖分獲得的凸包的正確性．

該算法的改進(jìn)主要在于步驟(2)的行走方法．傳統(tǒng)算法使用的是線性行走方法．對于輸入點(diǎn) P，指定一個四面體，如果 P位于該四面體內(nèi)，則完成行走；如果不在該四面體內(nèi)，就構(gòu)建一條射線，原點(diǎn)是當(dāng)前四面體的體內(nèi)一點(diǎn)，標(biāo)記為 C，它的方向?yàn)?C→P，定位到與該射線相交的一個四面體的面，與當(dāng)前四面體相鄰并共用這個面的四面體就是包含 P的下一個候選四面體．這個方法雖然快捷，但有個缺陷，即當(dāng)射線穿過四面體的頂點(diǎn)或邊時(shí)與射線相交的下一個四面體不是相鄰的四面體，導(dǎo)致 P的定位錯誤，而隨機(jī)行走方法有效地避免了該現(xiàn)象[5]．

1.1.3 四面體剖分實(shí)例

測試用的部分散亂點(diǎn)云實(shí)例見圖 1，點(diǎn)云數(shù)據(jù)來自 CGAL標(biāo)準(zhǔn)圖庫[6]．圖 2為由圖 1(b)散亂點(diǎn)云生成的物體凸包以及四面體網(wǎng)格．計(jì)算所有四面體的體積總和就可以得到物體凸包的體積，但是，凸包中還存在大量體外四面體，需要對體外四面體進(jìn)行剔除，才能夠?qū)崿F(xiàn)準(zhǔn)確的體積計(jì)算．

圖1 測試用散亂點(diǎn)云Fig.1 Scattered point cloud for test

圖2 四面體剖分結(jié)果Fig.2 Results of tetrahedron subdivision

上述算法在最壞情況下的時(shí)間復(fù)雜度為O(n2)，在輸入點(diǎn)集為一般物體表面點(diǎn)集的情況下，算法復(fù)雜度接近O(n lg n)．使用 C++語言實(shí)現(xiàn)以上算法進(jìn)行四面體剖分的耗時(shí)數(shù)據(jù)見表 1，運(yùn)行環(huán)境是 Windows XP,32位系統(tǒng)，編譯環(huán)境是 Visual Studio 2008，編譯器為VC++9.0，硬件為Core2 Duo E7400 2.80 GHz，2.0,GB DDR2 800 MHz RAM．

表1 四面體剖分耗時(shí)Tab.1 Time-consuming of tetrahedron subdivision

1.2 法向量的計(jì)算

由于四面體剖分實(shí)質(zhì)上是對凸包的分割，因此部分四面體是體外的冗余數(shù)據(jù)，需要剔除．由于散亂點(diǎn)云本身不包含任何拓?fù)湫畔ⅲ苯哟_定體外四面體比較困難，因此本文借由計(jì)算散亂點(diǎn)—四面體的頂點(diǎn)的法向量來描述四面體所處位置，從而剔除體外四面體．上述計(jì)算主要分為兩步：基于擬合曲面計(jì)算法向量和法向量一致化．

1.2.1 基于擬合曲面計(jì)算法向量

為了計(jì)算各點(diǎn)的法向量，用二次曲面來擬合點(diǎn)和它的K近鄰(指在歐幾里德距離下，包含n個點(diǎn)的數(shù)據(jù)集S?Rd中，找到一個點(diǎn)的 K個最近點(diǎn))．本文使用最小二乘法計(jì)算二次曲面參數(shù)，二次曲面方程可表示為

對于給定散亂點(diǎn)(xi，yi，zi)i＝1，2，…，N，使總誤差 Q最?。?/p>

根據(jù)式(3)可得到方程的參數(shù)，得到二次曲面方程后就可以通過偏微分來計(jì)算該點(diǎn)的法向量[7]．

1.2.2 法向量一致化

因?yàn)槌跏挤ㄏ蛄康闹赶騼?nèi)外不一，得到初始法向量后需要一致化調(diào)整，這個問題可以模型化為圖的優(yōu)化．這個圖的節(jié)點(diǎn)由各散亂點(diǎn)構(gòu)成，如果兩個節(jié)點(diǎn)足夠近的話則構(gòu)建一條邊，邊的權(quán)值用 ni·nj表示，這樣法向量一致化就可歸結(jié)為圖的總權(quán)值最大化問題，這里的主要問題是選取連接圖中的哪對節(jié)點(diǎn)．由于物體表面可以認(rèn)為是單連通元，所以圖應(yīng)該是連通的．最小生成樹就是一個連接鄰近點(diǎn)的簡單連通圖，因此基于散亂點(diǎn)的歐幾里得距離建立最小生成樹．由于初始最小生成樹邊的密度不能滿足本文要求，因此需要添加邊，如果兩個節(jié)點(diǎn)中的一個點(diǎn)是另一個點(diǎn)的 K近鄰則連接兩點(diǎn)，由此生成的圖稱為黎曼圖(Riemannian graph)；然后，再基于該圖生成最小生成樹，選擇一個種子法向量，按深度優(yōu)先準(zhǔn)則在最小生成樹中傳播調(diào)整法向量方向[8]．調(diào)整規(guī)則就是：對于最小生成樹中的兩個近鄰點(diǎn)，如果相應(yīng)法向量的點(diǎn)積 ni·nj＜0，則 ni或 nj應(yīng)當(dāng)被反向．

這一步的時(shí)間復(fù)雜度為O(n2)，近鄰個數(shù) K的選擇是影響法向量精度的一個重要因素，由于對象三維拓?fù)涞膹?fù)雜度不同，可將 K設(shè)置成可調(diào)整的變量．根據(jù)經(jīng)驗(yàn)，在大部分表面比較平滑的情況下 K可以取偏小的值；如果表面紋理復(fù)雜，K可以取偏大的值．

圖 3是奶牛點(diǎn)云的法向量效果圖(經(jīng)過多次測試，這里將 K設(shè)置為 18，效果較好)．可以看到，在散亂點(diǎn)比較密集的頭部和局部平滑的軀干部獲得的法向量準(zhǔn)確度比較高，而散亂點(diǎn)比較稀疏的腿部和尾部得到的法向量則準(zhǔn)確度相對低一些，部分點(diǎn)很難獲得指向正確的法向量．

圖3 奶牛點(diǎn)云的最終法向量Fig.3 Ultimate normal vectors of milk cow point cloud

1.3 體外四面體的剔除

獲得點(diǎn)云的法向量以后，體外和體內(nèi)四面體可通過法向量的指向判斷：體內(nèi)四面體頂點(diǎn)的法向量指向都背離四面體，而體外四面體頂點(diǎn)的法向量則指向四面體內(nèi)部或附近區(qū)域．實(shí)驗(yàn)表明，在獲得正確法向量的前提下，利用頂點(diǎn)的法向量是否和四面體外接球相交來剔除體外四面體可以取得理想的結(jié)果．具體算法是：遍歷四面體，計(jì)算每個四面體的外接球，然后從四面體的每個頂點(diǎn)沿著法向量方向引出一條射線，如果4個頂點(diǎn)的射線都與外接球有交點(diǎn)(不包含射線原點(diǎn))，可以認(rèn)為該四面體為體外四面體，不予考慮；否則，計(jì)算該四面體的體積并加入總體積．如圖 4(b)所示為圖 4(a)中四面體的外接球和 4個頂點(diǎn)的沿法向量指向的射線，射線與外接球有 4個交點(diǎn)，就可以判斷其為體外四面體．

圖4 四面體剔除原理圖Fig.4 Schematic of tetrahedron removing

這一步的時(shí)間復(fù)雜度為 O(n)，經(jīng)大量程序驗(yàn)證，在能得到基本正確的法向量的前提下，本算法能剔除絕大部分體外四面體，證明了上述方法的可行性．

1.4 體積的計(jì)算

經(jīng)過以上步驟，最后計(jì)算物體體積就變得非常簡單，只需要根據(jù)式(4)計(jì)算形體內(nèi)各四面體的體積之和，就可以得到物體的體積

2 實(shí) 驗(yàn)

基于上述方法，對圖 1中散亂點(diǎn)云進(jìn)行了測試，并與傳統(tǒng)算法比較．由于傳統(tǒng)算法的用時(shí)主要在構(gòu)建表面網(wǎng)格模型過程，所以比較時(shí)忽略了傳統(tǒng)算法的后續(xù)體積計(jì)算用時(shí)．本文選擇被廣泛使用的 Poisson重建法[9]重建表面網(wǎng)格模型，其將表面重建問題表示為 Poisson方程的解，通過估計(jì)模型的指示函數(shù)和提取等值面，對表面重建一個無縫的三角逼近．這樣做有諸多優(yōu)點(diǎn)，很多對隱式表面進(jìn)行擬合的方法把數(shù)據(jù)分割到不同區(qū)域進(jìn)行局部擬合，然后使用合成函數(shù)進(jìn)一步合成這些局部擬合結(jié)果，相比而言，Poisson重建是一個全局方法，不借助于啟發(fā)式的分割或合并，這樣它能創(chuàng)建較光滑的表面，穩(wěn)健地近似含有噪聲的數(shù)據(jù).

先使用本文算法計(jì)算物體體積，再使用 Poisson方法重建網(wǎng)格模型，比較兩者耗時(shí)，然后比較散亂點(diǎn)云提供者給出的體積信息與本文計(jì)算結(jié)果，如表2所示．對于 3個實(shí)例的散亂點(diǎn)云，本算法的總計(jì)算時(shí)間僅為Poisson重建時(shí)間的20.9%～40.2%，準(zhǔn)確率達(dá)到94.93%～98.07%．

表2 算法驗(yàn)證與比較Tab.2 Algorithm verifying and comparison

3 結(jié) 語

針對傳統(tǒng)算法先構(gòu)建表面網(wǎng)格模型再計(jì)算物體體積的流程，本文提出一種基于散亂點(diǎn)云直接計(jì)算物體體積的方法．本算法通過基于凸包分割的四面體剖分、法向量的計(jì)算、體外四面體的剔除和體積的計(jì)算等步驟，在保證準(zhǔn)確度的前提下，大幅度提高了效率，使實(shí)時(shí)的快速體積計(jì)算成為可能．

算法的準(zhǔn)確率和效率還有待提高，考慮在復(fù)雜點(diǎn)云的法向量計(jì)算及體外四面體剔除等方面作深入研究，進(jìn)一步提高算法的準(zhǔn)確率和效率．

［1］ 韋進(jìn). 三維空間任意多面體體積的一種坐標(biāo)計(jì)算法[J]. 湖州師專學(xué)報(bào)：自然科學(xué)，1997，19(5)：67–73.

［2］ 龐明勇，戴文俊. 基于體積分布特征匹配的三維實(shí)體網(wǎng)格模型檢索[J]. 系統(tǒng)仿真學(xué)報(bào)，2007，19(1)：30–34.

［3］ 王泉德. 任意三角網(wǎng)格模型體積的快速精確計(jì)算方法[J]. 計(jì)算機(jī)工程與應(yīng)用，2009，45(18)：32–34，58.

［4］ Edelsbrunner H，Seidel R. Voronoi diagrams and arrangements[J]. Discrete & Computational Geometry，1986(1)：25–44.

［5］ Devillers O，Pion S，Teillaud M. Walking in a triangulartion[J]. International Journal of Foundations on Computer Science，2002(13)：181–199.

［6］ European Research Programs. CGAL[EB/OL]. [2010–03–20]. http：//www. cgal. org/download. html.

［7］ 魏永超，蘇顯渝. 基于方向角的散亂點(diǎn)云三角剖分算法[J]. 四川大學(xué)學(xué)報(bào)：工程科學(xué)，2009，41(4)：202–207.

［8］ Hoppe H，DeRose T，Duchamp T，et al. Surface reconstruction from unorganized points[J]. Computer Graphics，1992，71(8)：71–77.

［9］ Kazhdan M，Bolitho M，Hoppe H. Poisson surface reconstruction[C]// Polthier K，Sheffer A. Symposium on Geometry Processing. Switzerland：The Eurographics Association，2006：61–70.

Algorithm of Rapid Volume Calculation Based on Scattered Point Cloud

HU Xiao-tong，TAO Sen-bai
(College of Computer Science and Information Engineering，Tianjin University of Science & Technology，Tianjin 300222，China)

TP391

A

1672-6510(2011)01-0067-05

2010–10–18；

2010–12–03

胡曉彤（1971—），男，北京人，副教授，huxt@tust.edu.cn.