3D-power圖的快速生成方法

桂志強，姚裕友，張高峰，徐本柱，鄭利平

（合肥工業(yè)大學 計算機與信息學院，安徽合肥230009）

0 引 言

Voronoi圖在日常生活和自然界中普遍存在，因其優(yōu)良的幾何特性，被廣泛應用于可視化、幾何建模、建造設計等領域。power圖是Voronoi圖的擴展，即對Voronoi圖站點引入權重概念，并重新定義距離。對power圖的區(qū)域施加容量限制，得到容量限制power圖（capacity constrained power diagram，CCPD）。進一步對CCPD站點施加質心限制，得到質心容量限制 power圖（centroidal capacity constrained power diagram，CCCPD）。

Power圖因其優(yōu)良的特性被廣泛應用于2D平面中的采樣和點畫［1］、藍噪聲［2-4］、計算機動畫［5］、3D空間流體仿真［6-8］等。

AURENHAMMER［9-10］提出power圖的概念，并對power圖的性質、應用進行了總結，提出用分治法構建power圖；BALZER等提出用“試位法”優(yōu)化 容 量 約 束，給 出 有 限 區(qū) 域［11］和 連 續(xù) 區(qū) 域［12］的CCPD算法，但時間復雜度較高，且收斂性差；GOES等［2］提 出 用 牛 頓 法 優(yōu) 化 權 重，并 結 合MULLEN等［13］提出的自適應步長梯度下降法優(yōu)化質心位置，二者交替迭代，生成CCCPD，但優(yōu)化過程中存在相互干擾，收斂速度較慢；XIN等［14］提出一種超線性收斂算法，將L-BFGS與牛頓法相結合，生成CCCPD。

相對2D領域，3D-power圖的生成復雜性大幅提升，目前大部分研究聚焦于3D-Voronoi圖［15-16］，并未將其推廣至power圖領域，其可能性和適應性也未得到充分證明；YAN等［17］提出了一種以計算幾何算法庫（computational geometry algorithms library，CGAL）［18］為基礎，用邊界剪裁方法計算3DVoronoi圖 的 算 法；RAY等［19］提 出 了 無 網(wǎng) 格3DVoronoi圖的圖形處理器（graphics processing unit，GPU）算法；LIU等［20］在YAN等［17］的基礎上提出了基于GPU的3D-Voronoi圖計算算法。

GOES等［6］實現(xiàn)了基于VORO++［21］的線程安全的并行power圖構造算法，計算效率得到較大提高；ZHAI等［8］根 據(jù)VORO++提 出 了 一 種 基 于GPU的裁切并行算法。這些算法適用于站點分布較為均勻的情況，如流體仿真的power圖，對隨機分布的站點表現(xiàn)不佳。GOLIN等［22］證明了在N個隨機站點中3D-Voronoi圖的復雜度（即頂點、邊和面的數(shù)量）達到O（n2），而如果站點均勻獨立分布，則復雜度可降至O（n）。對于隨機分布的站點，目前比較穩(wěn)定和有效的是基于CGAL和VORO++的power圖構造方法。

RONG等［23-24］提出了基于GPU的跳轉泛洪算法（jump flooding algorithm，JFA）構造Voronoi圖；ZHENG等［25］將其擴展至2D平面power圖，用連通域算法提取power圖的頂點和邊以優(yōu)化生成CCCPD，算法的時間效率較高，可擴展至3D空間，但因3D-power圖遠較2D平面圖復雜，連通域算法并不適合計算3D-power圖中的頂點、邊和面。因此，本文以ZHENG等［25］提出的基于GPU的構造算法為基礎，將其擴展至3D空間，并提出了一種新的估值算法，以優(yōu)化生成3DCCCPD。

本文的主要貢獻：

（1）將JFA擴展應用于3D空間power圖的生成。

（2）提出一種估值算法，計算離散3D-power圖中各區(qū)域的面積，與Lloyd算法和牛頓法結合優(yōu)化生成3D-CCCPD。

2 相關工作

2.1 Vor onoi圖與power圖

定義Voronoi圖在域Ω∈E d內(nèi)，基于直線距離d(x，x i)=‖x-x i‖對站點集X={x i|x i∈Ω，i=1，2，…，n}進行區(qū)域劃分。每個Voronoi區(qū)域vi定義為

Power圖作為Voronoi圖的擴展，為每個Voronoi站點x i賦予權重w i，本文用歐氏距離平方重新定義距離d：

對 于 新 站 點 集 (X，W)={(x，x i)|i=1，2，…，n}，將在域Ω∈E d內(nèi)的power劃分定義為

圖1 3D-power區(qū)域Fig.1 3D-power cell

值得注意的是，當各站點權重相等時，power圖便退化為Voronoi圖［9］。

2.2 質心容量限制power圖

引入站點權重，使得power圖的容量可控。對其施加精確容量限制，即可得容量限制power圖。設在問題域Ω∈E d，d=3中，存在連續(xù)密度函數(shù)ρ(x)，則定義power區(qū)域Pi w i的容量為

各power區(qū)域容量的和為

在3D空間，當密度函數(shù)為常數(shù)時，各power區(qū)域的容量即為其體積。

給 定 預 設 容 量 集 合C={ci|i=1，2，…，n}，中的權重，使得每個區(qū)域的容量滿足mi=ci，即得到嚴格滿足容量限制的power圖。

在某些情況下，可能會出現(xiàn)站點不在所屬power區(qū)域的情況，達不到預期的應用效果。為此，進一步對容量限制power圖施加質心約束，使每個站點x i均在其power區(qū)域Pi w i的質心bi內(nèi)，即3DCCCPD。

2.3 JFA算法

JFA算法最初由RONG等［23］提出，隨后陸續(xù)得到了JFA的性質和變體（1+JFA，JFA+1，halving resolution等）［24］。因受當時硬件條件限制，僅給出了3D-Voronoi圖的簡單擴展。對N×N×N的三維離散空間，首先考慮將其轉化為N個N×N的二維平面，然后將各個站點映射至N個平面，逐層調(diào)用JFA。JFA無法保證每次都能生成正確的Voronoi圖。RONG等［23］解釋了JFA產(chǎn)生誤差的原因，提出了1+JFA算法，并用實驗驗證了1+JFA算法可有效降低JFA算法的錯誤率。RONG等［24］提出了另一種JFA變體，即halving resolution算法，在相同分辨率下，其較一般JFA速度更快。

用JFA算法構造3D-power圖的過程如圖2所示。將每個像素P（x，y，z）所存儲的“信息”（站點坐標和權重）以步長l擴散至像素點P（x+k1×l，y+k2×l，z+k3×l），其中，k1，k2，k3∈{-1，0，1}且l∈{n/2，n/4，…，1}。每個像素點在接收到“信息”后按式（3）計算power距離，并保存與power距離最小的站點“信息”，繼續(xù)迭代，后輪迭代步長取前輪迭代步長的1/2，……，直至l=1，最終生成3Dpower圖。

圖2中，（a）為初始狀態(tài)，在8×8×8的立方體中，左前下角有一個種子點。首先，以步長為4擴散為狀態(tài)（b），然后，狀態(tài)（b）以步長為2擴散為狀態(tài)（c），最后，狀態(tài)（c）以步長為1擴散為狀態(tài)（d）。即N×N×N的立方體只需經(jīng)過logN次迭代便可將種子點擴散至整個立方體。

圖2 JFA算法在3D空間的迭代過程Fig.2 The iterative process of JFA algorithm in 3D space

3 基于GPU的3D-CCCPD生成算法

隨著power圖的應用領域越來越廣，問題域的規(guī)模越來越大、限制條件越來越復雜、生成速度要求越來越快，尤其是3D空間，以CGAL為基礎基于CPU的power圖構造方法已無法滿足實際應用需要。為此，本文引入GPU，以提高3D-power圖的構造效率，并結合CPU中的Lloyd算法和牛頓法優(yōu)化質心和容量，以大幅提升3D-power圖的生成效率，使其能適應大規(guī)模站點等三維應用場景。

3.1 問題描述

從能量優(yōu)化角度，對質心容量限制power圖，定義能量函數(shù)F(X，W)：

文獻［2-6］證明了能量函數(shù)F(X，W)：

本文將JFA變體halving resolution算法擴展至3Dpower圖的計算。步驟如下：

（1）將N×N×N分辨率下的初始像素圖映射至N/2×N/2×N/2分辨率下（即將8個像素映射在一個像素中，如果此8個像素中有多個站點值，則只選擇其中1個）。

（2）進行一般JFA迭代，生成3D-power圖。

（3）將N/2×N/2×N/2分辨率下生成的3Dpower圖再次映射至N×N×N分辨率下。

（4）進行一次步長為1的JFA迭代，以平滑邊界。

再通過Lloyd算法和牛頓法迭代求解能量函數(shù)F(X，W)的最小值，加速生成3D-CCCPD，提高3D-power圖的適用性。

3.2 估值算法

ZHENG等［25］提出的連通域算法，將由JFA生成的2D平面power圖紋理中的像素點劃分為4種情況，據(jù)此提取power圖的頂點，連接得到具有幾何結構數(shù)據(jù)的power圖。在3D-power圖紋理中，很難對像素點的分布進行具體劃分。即使用能夠提取3D-power圖的頂點構造三維凸包，也會嚴重影響時間效率。因此，需提出一種適合3D空間的計算power圖幾何數(shù)據(jù)的方法。

考慮3D-power圖生成的最小計算單位為像素點，由式（8）知，在用牛頓法優(yōu)化容量時較難計算相鄰power區(qū)域相交的面積，而power圖中邊和面均由像素點構成，若能直接統(tǒng)計相交面上像素點的個數(shù)，并以此作為其面積的估值，則能達到優(yōu)化3Dpower圖容量的目的。

圖3中，（c）為（a）中紅色標記區(qū)域平滑邊界前的放大圖，（d）為（b）中紅色標記區(qū)域平滑邊界后的放大圖。觀察到在JFA變體halving resolution算法中，最后一步以步長為1迭代平滑邊界，只計算了power圖“邊界”（2D-power圖中為邊，3D-power圖中為面）像素點。所以，對平滑邊界，只需對這部分像素點進行標記，統(tǒng)計其數(shù)量并作為當前power區(qū)域的面積估值。在平滑邊界的同時計算power圖的幾何數(shù)據(jù)，以提高時間性能。

圖3 JFA平滑邊界Fig.3 Smooth boundary of JFA

以圖4的估值算法為例，在問題域內(nèi)有紅色（x r），綠色（x g）和藍色（x b）3個站點，經(jīng)JFA和halving resolution算法迭代，得到如圖4所示的結果。在像素點P1內(nèi)存儲有x b的信息，在像素點P2內(nèi)存儲有x r的信息。當進行步長l=1的迭代時，經(jīng)計算，d b1

圖4 估值算法Fig.4 Estimating algorithm

算法1 估值算法

輸入 3D-power圖和站點數(shù)n

輸出 3D-power圖和A ij

Step 1初始化統(tǒng)計數(shù)組count[n][n]

Step 2 將保存在像素點P(x，y，z)中的站點ID擴散至P(x+k1，y+k2，z+k3)，k1，k2，k3∈{-1，0，1}

Step 3 由式（3），計算power距離di和dj

Step 4 如果dj

Step 4.1 count[i][j]++

Step 4.2 更新站點中保存的“信息”

Step 5 結束

Step 6 //統(tǒng)計，取同一個面估值的均值

Aij=(count[i][j]+count[j][i])/2

Step 7 輸出：3D-power圖及Aij

3.3 3D-power圖構造算法

用1+JFA和halving resolution算法構造3Dpower圖。在構造過程中調(diào)用估值算法，計算每個power區(qū)域的面積Aij。

算法2 基于GPU的3D-power圖構造算法

輸入 站點集：X={xi，i=1，2，…，n}，

權重集：W={wi，i=1，2，…，n}

輸出 3D-power圖

Step 1 //初始化GPU紋理，為每個站點分配唯一ID∈{1，2，…，n}

Step 2 重復

Step 2.1 若當前像素點中存在站點，則

Step 2.1.1 將站點ID保存在此像素點中

Step 2.2 且

Step 2.2.1 像素點中保存“-1”

Step 2.3 結束

Step 3 直到所有像素點完成初始化

Step 4 將3D紋理分辨率減半至N=N/2

Step 5 調(diào)用JFA_Kernel（），l=1

Step 6 對l=N/2；l≥1；l=l/2

Step 6.1 調(diào)用JFA_Kernel（），l

Step 7 結束

Step 8 將得到的紋理映射到N=2N的3D紋理中

Step 9 調(diào)用算法1//估值算法

Step 10 輸出渲染后的3D-power圖

Step 11 ///子程序JFA_Kernel（）

Step 12 //GPU并行地將“信息”傳播至每個像素，為每個站點分配唯一ID

Step 13 //26個傳播方向，將保存在像素點P(x，y，z)中的站點xi的ID擴散到P(x+k1×l，y+

k2×l，z+k3×l)，k1，k2，k3∈{-1，0，1}

Step 14 由式（3），計算power距離di與dj

Step 15 若dj

Step 15.1 將P(x+k1×l，y+k2×l，z+k3×l)中的“信息”更新為站點xj的ID

Step 16 結束

3.4 3D-CCCPD算法

3D-CCCPD的容量優(yōu)化算法——牛頓法，其海森矩陣的數(shù)據(jù)來源于GPU端的估值算法，但牛頓法迭代步長的計算運行于CPU端。由于涉及矩陣求逆和矩陣相乘，計算較為耗時。

算法3 3D-power圖優(yōu)化生成算法

輸入 站點集：X=｛xi，i=1，2，…，n｝

權重集：W=｛wi，i=1，2，…，n｝

容量集：C=｛ci，i=1，2，…，n｝

容量精度停止條件：tw

質心精度停止條件：tx

輸出 3D-CCCPD

Step 1 重復

Step 1.1 GPU：調(diào)用算法2構建3D-power圖

Step 1.2 重復

Step 1.2.1 //牛頓法優(yōu)化容量CPU：計算牛頓法迭代步長δw

Step 1.2.2W←W+δw

Step 1.2.3 GPU：調(diào)用算法2構造3D-power圖

Step 1.3 直到‖?wF‖≤tw

Step 1.4 //GPU-Lloyd算法優(yōu)化質心GPU：并行計算每個區(qū)域的質心bi

Step 1.5X←b

Step 2 直到‖?xF‖≤tx

Step 3 輸出3D-CCCPD

4 實驗分析

實 驗 環(huán) 境 為Microsoft Windows 10，Intel（R）Core（TM）i7-4720 CPU@2。60 GHz，16.0 GB RAM。編譯環(huán)境為Microsoft Visual C++2012，GPU為NVIDIA GeForce GTX 960M，GPU并行計算平臺為CUDA 8.0。實驗結果用MATLAB R2019b渲染。

4.1 分辨率設置

圖5為不同分辨率及不同站點數(shù)下用算法2構造3D-power圖的時間性能對比?？梢?，在相同分辨率、不同站點數(shù)下，算法2的時間性能相對穩(wěn)定；在相同站點數(shù)下，構造3D-power圖的時間隨分辨率的增加而增長。因此需選取一個合適的分辨率。由圖5，綜合時間性能和準確率，選取分辨率N=256的正方體作為問題域。

4.2 實驗結果

設質心精度tx=2，容量精度tw=0.000 3。圖6為生成的3D-CCCPD立體圖，第1列展示的為站點數(shù)分別為1 000，2 000，5 000時的3D-CCCPD立體圖效果，第2列展示的為站點數(shù)分別為1 000，2 000，5 000時的3D-CCCPD透視圖，第3～5列為Z坐標軸逐層拆解得到的X-Y平面圖，從左到右分別為第20，80和160層（共N=256層平面圖）?？芍?，本文算法生成的3D-CCCPD其結構較均勻緊湊。

圖5 其于GPU構造3D-power圖的算法性能Fig.5 The computational performance of the 3D-power diagram constructing algorithm based on GPU

4.3 算法性能

表1為不同站點數(shù)下，本文算法2生成power圖的時間性能對比?？芍?，在不同數(shù)量級的站點數(shù)下，算法2生成power圖的時間性能相對穩(wěn)定，適用于較多站點情況下3D-power圖的生成。

圖6 3D-CCCPD立體圖及截面圖Fig.6 The stereogram and cross-sectional views of 3D-CCCPD

表2給出了以CGAL為基礎的算法與本文算法3及生成3D-CCCPD的時間性能對比。其中，TCGAL為以CGAL為基礎的算法生成3D-CCCPD的時間性能，Tours為本文算法3生成3D-CCCPD的時間性能，容量精度停止條件為tw=0.000 3，質心精度停止條件為tx=2，即2個像素點的距離。定義CGAL與本文算法3的加速比（sp）為當站點數(shù)由2 000增至5 000時，加速比下降，這是因為本文算法3在站點數(shù)為5 000時，時間消耗發(fā)生了突增，其中，在基于CPU平臺的牛頓法優(yōu)化過程中，迭代步長的計算時間突然增加，其受站點數(shù)影響較大。另外，本文算法3的時間消耗隨站點數(shù)增加而增加，較以CGAL為基礎的算法優(yōu)勢明顯，且隨著站點數(shù)的增加，加速效果更明顯。

表1 不同站點數(shù)下本文算法2生成3D-power圖的時間性能Table 1 Computational performance of 3D-power diagrams construction of algorithm 2 with different number of sites

表2 不同站點數(shù)下以CGAL為基礎的算法與本文算法3的時間性能對比Table 2 Comparison results of algorithm based on CGAL and the proposed algorithm 3 with different number of sites

4 結 論

將JFA應用于3D-power圖的構造，提出了一種牛頓法優(yōu)化3D-CCCPD的容量估值算法，通過與JFA變體halving resolution算法結合，極大提高了生成3D-power圖的時間性能，結合最優(yōu)化算法，能夠快速生成滿足約束條件的3D-CCCPD。

實驗中發(fā)現(xiàn)，算法3的牛頓法迭代步長計算對最終的時間性能影響較大。下一步將繼續(xù)研究純GPU平臺下生成3D-CCCPD的算法，以及在保證時間性能的同時生成更高精度的3D-CCCPD。