桂志強(qiáng)，姚裕友，張高峰，徐本柱，鄭利平

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

0 引 言

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

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

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

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

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

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

本文的主要貢獻(xiàn)：

（1）將JFA擴(kuò)展應(yīng)用于3D空間power圖的生成。

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

2 相關(guān)工作

2.1 Vor onoi圖與power圖

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

Power圖作為Voronoi圖的擴(kuò)展，為每個(gè)Voronoi站點(diǎn)x i賦予權(quán)重w i，本文用歐氏距離平方重新定義距離d：

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

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

值得注意的是，當(dāng)各站點(diǎn)權(quán)重相等時(shí)，power圖便退化為Voronoi圖［9］。

2.2 質(zhì)心容量限制power圖

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

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

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

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

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

2.3 JFA算法

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

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

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

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

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

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

3.1 問(wèn)題描述

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

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

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

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

（2）進(jìn)行一般JFA迭代，生成3D-power圖。

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

（4）進(jìn)行一次步長(zhǎng)為1的JFA迭代，以平滑邊界。

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

3.2 估值算法

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

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

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

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

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

圖4 估值算法Fig.4 Estimating algorithm

算法1 估值算法

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

輸出 3D-power圖和A ij

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

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

Step 3 由式（3），計(jì)算power距離di和dj

Step 4 如果dj

Step 4.1 count[i][j]++

Step 4.2 更新站點(diǎn)中保存的“信息”

Step 5 結(jié)束

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

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

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

3.3 3D-power圖構(gòu)造算法

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

算法2 基于GPU的3D-power圖構(gòu)造算法

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

權(quán)重集：W={wi，i=1，2，…，n}

輸出 3D-power圖

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

Step 2 重復(fù)

Step 2.1 若當(dāng)前像素點(diǎn)中存在站點(diǎn)，則

Step 2.1.1 將站點(diǎn)ID保存在此像素點(diǎn)中

Step 2.2 且

Step 2.2.1 像素點(diǎn)中保存“-1”

Step 2.3 結(jié)束

Step 3 直到所有像素點(diǎn)完成初始化

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

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

Step 6 對(duì)l=N/2；l≥1；l=l/2

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

Step 7 結(jié)束

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

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

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

Step 11 ///子程序JFA_Kernel（）

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

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

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

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

Step 15 若dj

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

Step 16 結(jié)束

3.4 3D-CCCPD算法

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

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

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

權(quán)重集：W=｛wi，i=1，2，…，n｝

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

容量精度停止條件：tw

質(zhì)心精度停止條件：tx

輸出 3D-CCCPD

Step 1 重復(fù)

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

Step 1.2 重復(fù)

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

Step 1.2.2W←W+δw

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

Step 1.3 直到‖?wF‖≤tw

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

Step 1.5X←b

Step 2 直到‖?xF‖≤tx

Step 3 輸出3D-CCCPD

4 實(shí)驗(yàn)分析

實(shí) 驗(yàn) 環(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并行計(jì)算平臺(tái)為CUDA 8.0。實(shí)驗(yàn)結(jié)果用MATLAB R2019b渲染。

4.1 分辨率設(shè)置

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

4.2 實(shí)驗(yàn)結(jié)果

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

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

4.3 算法性能

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

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

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

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

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

4 結(jié) 論

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

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