一种多发生元Voronoi图的栅格生成方法

江 威,吴 艳 兰,谭 树 东,马 艺 文

(1.武汉市测绘研究院，湖北 武汉 430022；2.安徽大学资源与环境工程学院，安徽 合肥 230601；3.国家海洋信息中心，天津300171；4.武汉市国土资源和规划信息中心，湖北 武汉 430014)



一种多发生元Voronoi图的栅格生成方法

江 威1,吴 艳 兰2*,谭 树 东3,马 艺 文4

(1.武汉市测绘研究院，湖北 武汉 430022；2.安徽大学资源与环境工程学院，安徽 合肥 230601；3.国家海洋信息中心，天津300171；4.武汉市国土资源和规划信息中心，湖北 武汉 430014)

利用距离变换和栅格叠加分析，提出一种实现任意距离定义的2-site Voronoi图生成方法。首先进行距离变换得到距离图，然后通过邻近关系对边界进行划分得到2-site Voronoi图，最后将生成的距离图和2-site Voronoi图叠加。实验表明，该文提出的2-site Voronoi图生成方法可以快速构建多种距离类型和不同邻近关系下的Voronoi图，共生成了21种距离函数下的最远、最邻近和次邻近Voronoi图，解决了Voronoi图的多样性问题。该方法并不局限于点状发生元，可以生成任意形态发生元Voronoi图，并可以扩展生成N-site Voronoi图，生成的广义距离图可用来模拟成组的发生元在诸多约束条件下的区域增长过程。

距离函数；Voronoi图；距离变换；生长模型

0 引言

Voronoi图(简称V图)因其在空间邻接和空间分割中的独特性质，被广泛应用于计算机辅助设计、机器人路径规划、模式识别等领域，也是计算几何的一个重要研究内容。影响V图外观和属性的要素有发生元、距离类型和空间类型,其中，距离类型由空间类型决定。常规V图是由平面上n个发生元通过其欧氏距离最近空间点把平面划分成n个独立区域，其发生元是一组点集，每个Voronoi区域只包含一个点发生元，但这种V图并不适用于发生元具有预先分组信息的场合。本文引用文献[1]中的术语将常规V图称为1-site V图；类似地，将发生元具有预先分组信息的V图称为N-site V图，其主要外观特征是每个Voronoi区域对应着一组点发生元而不是某个发生元。与单发生元V图相比，多发生元V图除在发生元的形态和数量上不同之外，其距离定义也不同。单发生元V图中的距离值是每个像元到其最近单发生元的距离，而多发生元V图的距离含义由距离函数决定，如距离和函数对应于空间每个像元到多个发生元的距离和，此时的距离是由多个发生元共同决定的。

目前对于单发生元V图已有广泛研究[2,3],对于多发生元V图的研究则较少，并且主要研究了2-site V图，即成对发生元V图。点发生元成组之后会产生许多距离类型，因此与常规V图相比，2-site V图有更多的变形体。目前对于2-site V图的研究，主要是从理论上探讨不同距离类型对应的2-site V图的复杂度。Barequet定义的2-site的距离函数包括距离和、距离差、距离乘积、三角形面积等，并利用矢量生成方法中的分治法有效生成了对应的V图。但对于更复杂的距离函数，如这3个点组成的三角形的周长，其2-site V图有待进一步研究。对此，Hanniel证明了距离函数为三角形周长的2-site V图的计算复杂度为O(n2+ε)(ε>0)[4]。Dickerson等提出一个统一的模型形象地描述各种2-site V图，如距离和函数和三角形周长对应的V图，但没有介绍构建每个V图的方法[5]。最近，人们提出一些新的2-site距离函数的变形。例如，文献[6，7]提出“视角”函数，文献[7]还提出一些基于圆的距离函数(如3个点外接圆的半径)。Vyatkina等进一步研究了平面上曼哈顿距离和切比雪夫距离下的2-site V的复杂度，并可以扩展到闵可夫斯基距离度量和更高维的空间[8]。Dickerson等利用2-site V图，根据邻近关系，给地理网络图中的点加了标注[9]。

2-site V图具有很强的趣味性和特殊性。在生成方法上，Barequet等利用矢量方法中的分治法生成2-site V图，但存在算法复杂等矢量方法固有的缺点，此外，其距离函数主要是算术距离和几何距离(即三角形面积和三角形周长)，不具有灵活性。对于任意距离函数的2-site V图，这种基于矢量的生成方法并不实用，目前缺乏适用于任意距离函数的2-site V图的生成方法。本文利用距离变换和栅格叠加分析，提出一种实现任意距离定义的2-site V图生成方法，并将其扩展为适用于任意形态发生元、任意预先分组的多发生元V图。

1 方法原理

生成V图，最重要的就是进行距离变换，也就是模拟发生元的区域增长过程。类似的，生成2-site V图，首先也是进行距离变换得到距离图，由于此时控制距离尺度的距离函数有多种，就会得到多种扩展的距离图；然后，通过邻近关系对边界进行划分，从而得到2-site V图；最后将生成的距离图和2-site V图叠加，就得到了最终用来模拟区域增长过程的V图。本文利用距离变换和栅格叠加，提出一种2-site V图的生成方法，关键步骤为：1)由每个发生元分别生成距离图，称为单发生元距离图(记为1-site距离图)；2)对每对发生元，将其对应的两个1-site 距离图通过栅格叠加生成新的距离图，称为单对发生元距离图(记为1-pair-site距离图)，代表该距离图中只有一对发生元；3)对于包含多对的一组发生元，再次通过栅格叠加运算，将其对应的单对距离图合成为多对发生元距离图(记为N-pair-site距离图)；4)由N-pair-site 距离图得到对应的2-site V图。

值得说明的是：步骤2中栅格叠加对应的叠加算子为距离函数，步骤3中栅格叠加对应的叠加算子与步骤2不同，可定义为最小值(或最大值)函数，由此得到多对发生元的最近(或最远)距离图。为便于描述，将步骤2和步骤3中的栅格叠加分别称为第一次栅格叠加和第二次栅格叠加。本文利用ArcObjects开发了一个构建2-site V图和N-site V图的实验程序(图1),可快速实现欧氏距离变换、栅格叠加(两次)和V图可视化。

图1 2-site V图生成流程

Fig.1 Flow chart of generation of 2-site V diagram

1.1 任意距离类型的单对发生元距离图

单个发生元在平面上的增长过程可用围绕发生元的同心等距线的1-site距离图直观表示。对于2-site V图，由于区域增长过程中有两个中心，因而其等距线图的形状更复杂。Barequet等基于矢量方法，利用膨胀的同心椭圆和同心带等形状的等距线，生成了前文提到的几种距离函数对应的V图[2],但在其后来的文章里提到的更加复杂的距离函数，其等距线图很难用矢量的方法生成[5]。通过对1-site距离图的分析，并结合距离变换和栅格叠加的特点，可在图上直观地表示平面上的一对发生元的区域增长过程。与常规距离图相比，1-pair-site距离图是由一对发生元生成，因此每个像元值为像元到这对发生元的距离，而不是像元到一个发生元的距离。通过第一次栅格叠加可生成1-pair-site距离图。

如图2所示，p和q为空间上任意两点，其分别生成了1-site距离图。因此任意空间点v到点p和q的距离可通过查询这两个1-site距离图中的像元值获取。在此基础上，用不同的距离函数作为叠加算子，将两个1-site距离图中的像元值叠加，可生成一个新的距离图，即1-pair-site距离图。由于叠加算子由距离函数控制，叠加算子不同会直接影响1-pair-site距离图的外观和性质。因此，构建不同距离函数下的2-site V图，在第一次栅格叠加时，要使用对应的距离函数作为叠加算子。例如，加、减、乘算子可用来生成距离函数分别为和、差、乘积的1-pair-site距离图，并最终生成2-site V图。图2a为1-site 距离图的生成说明，叠加算子使用的是3×3窗口，并用不同距离函数指定的叠加算子生成各类1-pair-site距离图；图2b为一对点发生元p和q利用不同的叠加算子生成1-pair-site距离图的实例。

图2 1-pair-site距离图的生成

Fig.2 Generation of 1-pair-site DM

1.2 任意邻近关系的多对发生元距离图

与1-pair-site距离图类似，N-pair-site距离图是一种新的距离图，用来模拟平面上N个点对的区域增长过程。本文中，N-pair-site距离图的生成过程包含两个关键步骤：1)计算每个像元到所有成对发生元的距离值；2)给每个像元分配最小(或最大)距离值，可以生成N-pair-site最近(或最远)距离图。

步骤1可通过查询1-pair-site距离图获得,通过查询图中的像元值，可获得每个像元到这对发生元的距离。因此，从N组发生元生成的一组1-pair-site距离图，在每个1-pair-site距离图中，通过查询在相同位置的每个像元的值，可获得每个像元分别到这N组发生元的距离。

步骤2的详细说明见图3a。假设有4组发生元，ID分别为p1、p2、p3、p4，利用距离乘函数分别生成这4组发生元的1-pair-site距离图，并用ID值代表生成的4个1-pair-site距离图。对于每个像元，都可以查询到其在这4个1-pair-site距离图中的像元值，即为每个像素到这些在特殊的距离函数下生成的1-pair-site距离图的距离值。例如，这4个1-pair-site距离图对应的3×3像素窗口的左上角的值分别为0、2、2、2。执行了第二次叠加之后，每个1-pair-site距离图中的像素值通过一个特定的算子合并成一个值，并生成了4-pair-site距离图和2-site V图。同样以3×3窗口的左上角像素值为例，在使用最小算子的情况下，4-pair-site距离图的左上角像素值为0。图3b为生成4-pair-site距离图的实例，两次叠加所用的叠加算子分别为距离乘函数和最近的邻近关系。步骤2通过图4所示的第二次叠加完成，其与第一次叠加类似，但是添加的数据和叠加算子不同。第二次栅格叠加，将一组1-pair-site距离图作为添加的数据，在合成的新距离图中，1-pair-site距离图中的每个像元会生成一个新的像素值。这个合成图即为N-pair-site距离图。此处使用的叠加算子根据邻近关系的不同而不同。例如，最小(或最大)算子用来寻找最近(或最远)领域，并生成包含任意邻近关系的N-pair-site距离图。

图3 N-pair-site距离图和2-site V图的生成

Fig.3 Generation ofN-pair-site DM and 2-site V diagram

1.3 成对发生元V图的生成和可视化

生成距离图后，对应的V图可通过连接距离图中的最大距离值或提取等距线相遇处的像素生成。使用生成4-pair-site距离图时的相关信息，可生成2-site V图。如图3a所示，N-pair-site距离图中的每个像素的距离值均来源于某个1-pair-site距离图,这些像素值的来源决定了V图的空间像素分配。因此，在生成4-pair-site距离图的过程中，直接把像素值来源的1-pair-site距离图的ID值赋给该像素，即可同时生成2-site V 图。考虑3×3窗口左上角的像素值，由于4-pair-site中该像素的值0来源的1-pair-site距离图的ID值为p1，则2-site V图中该像素的值为p1。实际上，N-pair-site距离图中的距离值来源决定了Voronoi区域每个像元的值。改变叠加算子，如使用最远和次邻近的邻域关系，可生成对应的4-pair-site距离图和2-site V图(图3b)。

V图可将平面划分为若干区域或多边形，将每个点对的ID值对应的颜色渲染给对应的区域，可将V图可视化。渲染之后的Voronoi区域图中，能够迅速查询出每个区域所归属的点对和邻近关系，然而发生元的传播方式和每个发生元到各个区域的距离等信息无法从V图中获取。在一幅图中展现V图的发生元、传播距离和Voronoi区域这3个要素，就不会破坏V图的完整性。尤其是2-site V图，相比1-site V图，其与距离定义和邻近关系间的关系更密切。对于栅格叠加，关键步骤是把发生元和V图叠加在距离图之上，并将V图调整至一定的透明度。

2 实验结果

为验证该方法的灵活性，本文使用文献[2，5]中讨论的算术距离函数和几何距离函数及新定义的一些距离函数，共21种(表1)。利用这些距离函数，执行第一次栅格叠加，可生成对应的1-pair-site距离图(图4)。由表1和图4可见，21种距离函数对应的距离图存在形态相似情况，因此可将表1中的21种距离函数和对应的图6中的1-pair-site距离图分为几组，每组距离函数本质上相同。如图4(1)和图4(10)，其分别对应距离函数(f1)和2-site的三角形周长公式(f10)，因这两个函数的差为一常数，所以图4(1)和图4(10)的等距线形状相同。类似的还有3组：(f2，f20)，(f7，f9)，(f11，f14，f15，f16)。对这4组距离函数，每组只选择一个函数进行后续实验,则后续实验所用距离函数为16种。

表1 21种距离函数的定义

Table 1 The definition of twenty-one kinds of distance functions

距离函数说明f1d(v，p)+d(v，q)距离和f2|d(v，p)-d(v，q)|距离差的绝对值f3d(v，p)×d(v，q)距离乘积f4Max(d(v，p)，d(v，q))Min(d(v．p)，d(v，q))距离相除f5Max(d(v，p)，d(v，q))取较大距离f6Min(d(v，p)，d(v，q))取较小距离f7d(v，lpq)空间点v到直线lpq的距离f8d(v，pq)空间点v到线段pq的距离f9SΔvpqΔvpq的面积f10PcΔvpq带参数的Δvpq周长：PcΔvpq=d(v，p)+d(v，q)+c·d(p，q)且c≥-1。当c=1，PcΔvpq为Δvpq的周长f11R1c(v，p，q)Δvpq的外接圆半径距离函数说明f12R2cv，p，qΔvpq的内切圆半径f13R3Ovpq包含v、p、q的最小圆的半径f14d(o，pq)o到线段pq的距离，o为Δvpq内切圆圆心f15SΔopqΔopq的面积，o为Δvpq内切圆圆心f16pΔopqΔopq的周长，o为Δvpq内切圆圆心f17∠pvq夹角∠pvqf18w1·d(v，p)+w2·d(v，q)加权距离，w1和w2分别为d(v，p)和d(v，q)的权重f19d2(v，p)+d2(v，q)距离平方和f20Var(d(v，(p，q)))V到点p和q距离的方差：Var(d(v，p)，d(v，q))=[(d(v，p)-d)2+(d(v，q)-d)2]/2且d=(d(v，p)+d(v，q))/2f21Var(d(v，p，q))v，p，q3个点之间的3个距离的方差：Var(d(v，p，q))=[(d(v，p)-d)2+(d(v，q)-d)2+(d(p，q)-d)2]/3且d=(d(v，p)+d(v，q)+d(p，q))/3

图4 21种距离函数下的1-pair-site距离图

Fig.4 1-pair-site DMs under 21 kinds of distance functions

3 讨论

对本文方法的特性作如下讨论：1)由图4和图5可见，本文方法可以生成多种类型的2-site V图，从而实现2-site V图的多样性。通过改变距离函数和邻近关系，可以控制生成的2-site V图的外观和性质。本方法的主要特点是两次栅格叠加过程：第一次栅格叠加的叠加算子是定义的21种距离函数，第二次栅格叠加的叠加算子为邻近关系。通过改变第一次叠加的算子，不同的距离函数会生成不同结构的2-site距离图(图4)。类似地，改变第二次叠加的算子，可以生成不同邻近关系的N-pair-site距离图和2-site V图(图5)。由于对多种距离函数和邻近关系的一般控制，本文的方法具有灵活性和可控性。2)本文方法并不局限于点状发生元，其可以用于生成任意形态发生元的V图。这种基于栅格的处理方式，由于任何发生元(点、弧段和面)都按照相同的方法计算距离，因而极大地简化了不同形态发生元2-site V图的生成过程。3)本文方法可以生成N-site V图，此时的距离函数定义为每个点到多个发生元的距离。例如，3-site V图的一个距离和函数定义的距离为点到3个发生元的距离。只需改变第一次叠加，该方法可生成3-site V图。4)本文方法可以利用生成的距离图模拟多种距离函数和邻近关系下的成组发生元的区域增长过程(图5),且本文方法简单、健全且容易实现，其主要步骤是距离变换和栅格叠加,计算复杂度小，稳定性高。5)作为一种基于栅格的方法，当发生元数量比较多时，该方法比较费时。为了提高效率，可以结合并行计算等方法，从而节省运算时间。

4 结语

本文提出一种基于距离变换和栅格叠加的多发生元生成方法，生成了21种距离函数下的最远、次邻近和最邻近3种邻近关系下的2-site V图，并将其扩展生成为多发生元V图。实验表明：与现有矢量方法相比，本文方法具有普遍性和灵活性，具体体现在：1)可根据任意距离函数和邻近关系，生成多种2-site V图；2)方法并不局限于点发生元，可以应用于任意形态发生元(如弧段和面)；3)由于采用栅格运算方式，发生元由点类型扩展到任意形态类型，并没有增加算法的复杂度；4)可以在成对发生元V图2-site V图的基础上扩展生成N-site V图；5)除V图外，本文方法生成的广义距离图，可用来模拟成组的发生元在诸多约束条件下的区域增长过程。因此，它很可能成为分析生长模型在相关应用领域(如材料科学和植物生态学)的一个有力工具。

[1] BAREQUET G,DICKERSON M T,SCOT DRYSDALE R L.2-Point site Voronoi diagrams[J].Discrete Applied Mathematics,2002,122(1):37-54.

[2] CHEN Z G,XIAO Y Y,CAO J.Approximation by piecewise polynomials on Voronoi tessellation[J].Graphical Models,2014,76(5):522-531.

[3] EMIRIS I Z,MANTZAFLARIS A,MOURRAIN B.Voronoi diagrams of algebraic distance fields[J].Computer-Aided Design,2013,45(2)：511-516.

[4] HANNIEL I,BAREQUET G.On the triangle-perimeter two-site Voronoi diagram[A].6th International Symposium on Voronoi Diagrams in Science and Engineering[C].2009.

[5] DICKERSON M T,EPPSTEIN D.Animating a continuous family of two-site Voronoi diagrams (and a proof of a bound on the number of regions)[A].Proceedings of the 25th Annual Symposium on Computational Geometry[C].2009.

[6] ASANO T,TAMAKI H.Angular Voronoi diagram with applications[A].3rd International Symposium on Voronoi Diagrams in Science and Engineering[C].2006.

[7] BAREQUET G,DICKERSON M T,EPPSTEIN D,et al.On 2-site Voronoi diagrams under geometric distance functions[A].8th International Symposium on Voronoi Diagrams in Science and Engineering[C].2011.

[8] VYATKINA K,BAREQUET G.On 2-site Voronoi diagrams under arithmetic combinations of point-to-point distances[A].7th International Symposium on Voronoi Diagrams in Science and Engineering[C].2010.

[9] DICKERSON M T,GOODRICH M T.Two-site Voronoi diagrams in geographic networks[A].16th ACM SIGSPATIAL International Conference on Advances in Geopraphic Information Systems[C].2008.

[10] VACAVANT A.Fast distance transformation on irregular two-dimensional grids[J].Pattern Recognition,2010,43(10)：3348-3358.

[11] CHEN J.A raster-based method for computing Voronoi diagrams of spatial objects using dynamic distance transformation[J].International Journal of Geographical Information Science,1999,13(3)：209-225.

A Raster-Based Method to Generate Two-Site Voronoi Diagrams

JIANG Wei1,WU Yan-lan2,TAN Shu-dong3,MA Yi-wen4

(1.WuhanGeomaticInstitute,Wuhan430022;2.SchoolofResoucesandEnvironmentEngineering,AnhuiUniversity,Hefei230601;3.NationalMarineDataandInformationService,Tianjin300171;4.WuhanLandResourcesandPlanningBureau,Wuhan430014,China)

In contrast with the regular (1-site) Voronoi Diagram,the 2-site VD has more variants due to the change of the distance function.Generating the 2-site VD with respect to arbitrary distance functions and constraints remains a challenge in the field.This paper proposes a flexible and general approach to generate 2-site VD by combining the distance transform and raster overlay into a unified framework.This framework is characterized by the two procedures of raster overlay：the first overlay procedure whose operator is used to control the distance function,and the second overlay procedure whose operator is specified by the neighbor relationship.By manipulating the two overlay operators,the nearest-(furthest-,etc.) neighbor 2-site VDs with respect to the various distance functions can be obtained.The proposed approach was implemented and tested with 21 kinds of distance functions.The results show improved flexibility and robustness over existing vector-based approaches and emphasize the convenience of extending to general sites andN-site VD.The proposed approach also producesN-pair-site distance map：a new type of distance map,providing an easy and convenient way to simulate the region-growing process of the sites previously grouped into different pairs.

distance function;Voronoi diagram;distance transform;growth model

2014-12-22；

2015-03-11

国家自然科学基金项目(41271443);安徽省自然科学基金项目(1308085MD52)

江威(1988-)，男，硕士研究生，研究方向为地图学与地理信息系统。*通讯作者E-mail：wylmq@sina.com

10.3969/j.issn.1672-0504.2015.05.009

P208

A

1672-0504(2015)05-0039-05