三角网格上基于测地距离的近似Voronoi Diagram

2012311303             濮云飞
2012310636             许春旭

1  引言
测地线的一个定义可以简单表述为：曲面上的一条曲线，其主曲率方向与曲面的法向相同。直观上讲，这一定义说明这样的曲线，其弯曲程度是完全受限于曲面的，可以看作是经过曲面上两点间的"最直"的曲线（例如，一个球面上两点间的测地线就是经过这两点间的大圆弧）。测地距离就是两点测地线的长度。
三角网格是在计算机中离散化表示三维模型的常用方式。测地线在三角网格上有如下几个特点：(1)在三角形内部，测地线为直线段；(2)如果测地线不经过任何顶点，则将三角面片沿测地线经过的边展开，测地线仍为直线段。这种特殊情况的测地线称 为离散测地线（Discrete geodesic path）。需要注意的是测地线与沿曲面最短路径的不同：两点间可能有多条测地线，但并不是每一条都是最短路径。但在本文中，我们所讲的测地线（测地距离）可简单地认为就是最短路径（最短距离）。
Voronoi Diagram是对空间的一个划分，具体的定义随空间和度量的不同而不同。在三维模型上，基于测地距离的 Voronoi Diagram是一种对模型表面的自然划分形式：其对表面的划分表示了表面上各点到一组site中最近的关系。
本实验主要就三角网格上基于测地距离的近似Voronoi Diagram作一简要的探讨，对一已有算法进行了实现，并对其效果 和性能做了测试和评价。

2  相关工作
离散测地线作为几何中的重要概念，其相关算法在过去的很长一段时间内得到了比较充分的研究。
在理论上，[1]提出了一个由平面的Dijkstra算法演化而来的单源测地线算法（MMP算法），其运行方式与Dijkstra 算法十分类似，时间和空间复杂度为O(n2logn)和O(n2)，这里n为三维模型的复杂度（由于三维模型的顶点数、边数和面边相互成线性关系，所以这里没有刻意区分n的具体意义）。而[2]提出了另一个算法（CH算法），其运行过程类似于对一棵树 的breadth-first search，其时间和空间复杂度为O(n2)和O(n)，相对于MMP算法有了渐近意义上的改进。这两个算法是目前计算离散测地线的最主要的两个算法。
在实际上[3]首次实现了MMP算法，并和CH算法的实现进行比较，其结果说明在实际运行时，MMP算法的时间上界可能是相当不紧的，其时间性能要优于一个普通的CH算法的实现。而[4]中对CH算法提出了改进，通过一定的filter policy在实际计算中减少了计算量；另一方面，类似于MMP算法，通过引入Priority Queue的结构，在实际运行过程中进一步加快了改进后的CH算法的运算 速度。

虽然对于不同维度、不同度量下的Voronoi Diagram有很多关于算法、结构、复杂度等方面的工作，但我们没有找到一个可行的基于测地 距离的Voronoi Diagram的构造算法。这里我们采用的是[5]中提出的称为"ϵ-Voronoi Diagram"算法，该算法的基本框架是在MMP算法的基础上，通过不断局部细分三角网格形成拓扑一致、误差可控的近似Voronoi Diagram。

3  MMP算法简述
MMP算法是[5]中所提出的算法的基础，所以这一部分将对MMP算法作一全面而简要的说明。

MMP算法中的一个重要结构被称为"窗口"（Window）或"区间"（Interval）。[1]中的一个重要结论是：每条边是被这种窗口完全覆盖且窗口间没有重叠。同一窗口中的每个点到源点的最短路径会经过相同的三角面序列，如图1所示：

Figure 1: 同一窗口内的点到源点的最短路径经过相同的三角面序列

在说明MMP算法的运行过程之前，我们需要阐述"凹点"（convex vertex）和"伪源点"（pseudosource）的概念。
凹点是指所有相邻角度之和大于2π的顶点。如果一个三角网格不包含凹点，则任意两个点之间的测地线均不会经过任何顶点。而考虑一下包含凹点的三角网格，由于沿其周围角将相邻三角面展开，所有三角面无法展开到同一平面上而不重叠，可知存在这样的两个顶点，两者之间的最短路径必然经过凹点。这样，对于部分边上的点，其到源点的最短路径需经过该凹点。我们将这凹点称为伪源点，即可将其视为一个新的源点，从而使得部分点需经过该点回到源点。伪源点与源点的不同之处在于，从某点到达对应伪源点后尚未回到源点，离到源点的距离还有σ（σ为从伪源点到源点的最短距离）。这样，我们可以将所有的（伪）源点视为一个有σ域的结构，只不过 源点是σ值为0的特殊伪源点。

类比于平面图的Dijkstra算法，MMP算法的运行过程如下：
算法的整个运行过程中维护一个Priority Queue，其中的元素为Window，而两个Window之间的比较是通过对应的σ值决定的。初始时无窗口，只有一个源点。在与源点对立的所有边上建立一个恰好覆盖该边的窗口，并将其全部加入Priority Queue。
每次从Priority Queue中取出σ值最小的Window，按其对应的（伪）源点，将其投射到相邻三角形的边上生成新的窗口并压入到Priority Queue中；具体的情况见图2。其中黄色的窗口与原窗口有相同的（伪）源点，而红色的窗口以该坐标系的源点代表的顶点作为伪源点。

Figure 2: 将窗口投射到相邻三角形的边上

在投射窗口的过程中，新生成的窗口可能会与对应边上的原有窗口相交。这时需要对窗口进行划分，以保持该边上的窗口互不相交的性质，见图3 。划分的过程实际上是解一个一元二次方程。新生成的窗口加入到Priority Queue中。处理完窗口的划分后，再从Priority Queue中取出下一次待投射的窗口进行处理。这一迭代过程直至Priority Queue为空为止。

Figure 3: 窗口的划分

4  ϵ-Voronoi Diagram算法

4.1  算法基础
基于测地距离的Voronoi Diagram比一般意义下Rd中基于Lp度量的Voronoi Diagram要复杂，其中最重要的一点是：由于伪源点的存在，基于测地距离的VoronoiDiagram可能包含直线段、抛物线段、双曲线段等多种情况。另外，由于度量的计算不再是常数时间，无论是精确还是近似Voronoi Diagram的计算都在很大程度上依赖测地距离算法。

ϵ-Voronoi Diagram算法的基础是MMP算法。与CH算法相比，之所以以MMP算法为基础，除了在实际测试中后者的时间性能较好外，最重要的原因如下：
(1)虽然MMP算法是单源点的测地线算法，但通过前部分对MMP算法运行过程的分析，可知很容易将该算法扩展为多源点，即设定多个源点后，在迭代开始时在所有与源点相对的边上建立覆盖整个边的窗口并插入Priority Queue中；后面的过程与单源点算法完全相同。这一特点与Voronoi Diagram的多个site是相符的；
(2) MMP算法中，窗口的结构是连续的，即：同一条边上的相邻窗口在相邻点处，经过两个窗口中的任何一个到源点，其距离均等于最短距离。如果在窗口的结构中加一个域，表示其对应的源点的ID（注意此处是源点而非伪源点），这样，如果相邻两个窗口的该域不同，则说明两个窗口相 交处必为Voronoi Diagram边上的一点。

定义Voronoi Diagram与三角网格的边相交的点为关键点。对于一个三角形，其边上的某点是关键点当且仅当相邻窗口的交点，且两个相邻窗口对应的源点不同。因而原三角网格中，每个三角形周围的关键点是可以确切找出的。根据三角形周围关键点的个数，可分为如下两种情况：

如果对于某三角形，其周边仅有两个关键点，则可知必然有一条Voronoi Diagram的边（可能是有多段曲线连接成的"边"）通过这个三角形且经过这两个关键点。但不知道其在三角形内部的具体情况，如图4所示。

Figure 4: 三角形周围有两个关键点可能的情况

如果对于某三角形，其周边有三个或以上的关键点，则其内部的情况会更复杂：可能在三角形内部没有Voronoi Diagram边的交点，如图 5所示；可能在三角形内部Voronoi Diagram边有交点，但具体的相交情况也无法得知，如图6所示。

Figure 5: 三角形周围有多于两个关键点，且内部无交点

Figure 6: 三角形周围有多于两个关键点，且内部有交点

假设我们按如下规则直接连接现有三角形周围的关键点：若三角形周围有两个关键点，则将其直接连接；若有多于两个关键点，则选取三角面内部的某点（例如所有关键点的平均值），将其与所有关键点相连。这样的处理方式有如下两个问题：(1)拓扑连接关系与实际情况相差太大；(2) 误差无法控制。[5]中通过细分三角形处理来解决这两个问题。
4.2  局部细分
对三角形的一次细分过程如下：
(1) 对于待细分三角形，选取其最长边，将其等分为两条边，并将切分处作为一个新的顶点，与最长边相对的两个顶点相连构 成两条新边；
(2) 更新所有新生成的边、被更改的边上的距离场（即Window信息）。
具体细分的实现将在"实现细节"一部分中描述。

显然，对所有三角形都进行细分是没有必要的。一个容易想到的区分三角形是否需要进一步划分的判断是：如果三角形周围无关键点，则该三角形必然不需要进一步划分。对于周围有关键点的三角形，同样不是都需要进一步划分的。
对于周围仅有两个关键点的三角形，其内部的拓扑关系是可知的，我们仅需要考虑控制误差。这在种情况下，停止的准则可以简单地确定为当三角形的大小小于一个阈值（误差限决定）；对于周围有两个以上的关键点的三角形，细分准则要考虑到拓扑连接关系和误差这两个问题。 [5]中对于多于两个关键点的情况提出了一个解决这两个问题的停止细分准则，它是基于如下结论的
结论面片f内存在k条Voronoi分支的汇聚，当且仅当f的边上依逆时针序存在k个关键点p1,...,pk，以及{i1, i2, ..., ik } ∈ { 1, ..., m}，使得pj距离vij和vij+1最近，ik+1 = i1。
在这一结论下，我们很容易得知是否需要对某一周围有多于两个关键点的三角形进一步细分。综上，对三角形的细分准则如下：
(1) 三角形的大小小于某一值；
(2) 三角形周围仅有两个关键点，或有多于两个关键点且满足上面的结论。
当然，没有关键点时细分同样需要停止。

5  实现细节
前面我们描述了ϵ-Voronoi Diagram算法的运行过程。这一部分对实现细节做一介绍。
5.1  数据结构
如前所述，ϵ-Voronoi Diagram是建立在MMP算法基础上的。我们已有的MMP算法的实现使用了双向链接边表（DCEL）的结构，如下所示：

class CVertex{
public:
Vector3D m_vPosition;               //顶点的空间位置
unsigned int* m_piEdge;             //由顶点出发的所有半边(ptr)
short m_nValence;                   //顶点的度数
};

class CEdge{
public:
unsigned int m_iVertex[2];          //半边的起点(0)和终点(1)
unsigned int m_iTwinEdge;           //twin边(ptr)
unsigned int m_iNextEdge;           //next边(ptr)
unsigned int m_iFace;               //沿逆时针方向围绕的边(ptr)
};

class CFace{
public:
short m_nType;                      //面的边数，这里始终为3
unsigned int* m_piVertex;           //构成面的顶点列表，逆时针方向存储(ptr)
unsigned int* m_piEdge;             //构成面的半边列表，逆时针方向存储(ptr)
};

一个三角网格模型由三个数组组成，分别用于存放顶点（CVertex）、半边（CEdge）和面（CFace）。所有的DCEL中的链接关系均是由index而非pointer 实现的。

在这样的存储结构下，对一个三角形按某一条边的细分可以很容易实现。详细地实现可参考代码，其基本思路即是增加1个顶点、6条半边、2个面。新增的顶点位于被分割的边的中点；6条半边中，有2条是从原有的2条边中分离出来的，另外4条是新生成的边；2个面均是从与被分割的边相关的两个面分离出来的。另外，有2条原有边和2个原有面的相关信息需要修改，修改包括next关系、twin关系等等。

5.2  Window结构
参考[3]的实现，一个Window的构成见图7。程序中Window结构的定义如下：

struct Window{
unsigned int Edge_idx;
double b1, b2;
double d1, d2;
double sigma;
bitset<MAXSOURCENUM> srcId;
};

Figure 7: 窗口的构成

大部分域的意义与其名称一致。需要注意的是我们用C++ STL中的bitset结构来记录：对应的位为1标示该Window来自对应的源点。

在MMP算法中，一对互为twin的半边上的Window结构是不对称的，即两个半边上的窗口数量、信息并不一定相等，相同位置上的点通过两个半边上对应窗口回到源点的路径、长度也是不一定相等的。因而，为了更方便地寻找一个三角形周围的关键点，在算法开始的部分，以及后面更新Window结构后，需要对每一对twin的半边上的窗口进行合并。这一过程与下一部分所说明的内容是相同的。

5.3  更新新生成的边上的Window结构
在每次细分后，对于新增的6条和修改的2条边上的距离场（Window信息）进行更新。其更新方法如下：
对于修改的2条边，其Window结构只是从原有结构的基础上去除一部分Window而得到；新增的6条边中，由已有的2条边上剪裁下的2条边上的Window结构只需要将对应的部分简单拷贝过来即可；对于剩余的4条边，找到其参与组成的三角形的另外两条边上的窗口投射而形成的新窗口，这一过程中与投射窗口的顺序无关，最重要的一步是前面 提到的窗口的划分。

两个相交窗口的划分等价于求解一个一元二次方程（参见[3]）。这里我们需要解决的问题是：一个窗口可能和多个窗口同时相交。我们实现的思路是：对对应边上已存在的窗口遍历，将每一个窗口与新生成、待裁剪的窗口进行裁剪；裁剪的结果中，来自已有窗口的裁剪结果可以直接按顺序插入到对应边的对应位置上，而来自新窗口的裁剪结果作为待比较的窗口进入遍历的下一轮迭代中。

5.4  误差控制
关于误差控制有多种不同的理解和实现。程序中对误差限采取如下定义：对于某一绝对的error bound，精确Voronoi Diagram上的任一点 与ϵ-Voronoi Diagram上的某点相差不超过error bound；而这一error bound的定义是通过交互输入一个ϵ值并按 式error bound = ϵAvgEdgeLen计算的，其中AvgEdgeLen为所有边的平均值。

5.5  精度控制
实现中的精度控制采用了比较简单的策略，即设定一个小阈值，任何比较都通过这个阈值加以限定（如大于的判定至少需要大于一个阈值才能认为是大于）。但通过后面的分析可以看到，对于浮点误差的控制和处理还有待于进一步完善。

6  演示程序说明
程序使用C++实现，算法部分未使用任何库。界面采用Qt设计，图形显示采用OpenGL实现。程序的界面见图8。

Figure 8: 演示程序界面

演示程序可以读取一个.obj格式的三维模型文件并显示。与三维模型的交互通过鼠标完成：左键单击拖动以实现旋转；中键单击拖动以实现平移；滚轮滚动以缩放；右键打开菜单以选取顶点（其中，非顶点的目地点通过Alt键和鼠标右键点击来选取）。Voronoi Diagram的site通过鼠标右键选取，见图9。

Figure 9: 选取Voronoi Diagram的site

在选取了所有site后，单击菜单栏的"Geodesic"中的"Build Distance Field"运行MMP算法建立距离场；而后通过Lighting标签下 的文本框设置ϵ-Voronoi Diagram的误差限，再单击"Geodesic"中的"Build Voronoi Diagram"运行ϵ-Voronoi Diagram算法建立并显示Voronoi Diagram。此后，还可以过通过菜单项Show Splited Mesh来显示局部细分三角网格的结果，如图10 所示。

Figure 10: 三角网格局部细分的结果

由于程序中嵌入了MMP算法，可以通过构建测地线在一定意义下验证构造的Voronoi Diagram的正确性：选取某一点，右键使之成为目的 顶点；然后通过"Geodesic"菜单中的"Build Geodesic Path"菜单项建立从它到最近的site的测地线，可以验证该site是否是该目 的地所属的Voronoi Diagram cell对应的site.

演示程序的其余功能包括：模型颜色的选择；点、线框、面等不同显示模式；光照效果及环境光、散射光、镜面光成分的调节等。此类功能基本通过OpenGL实现，操作较简单，故此处不再赘述。

7  运行结果及算法分析
我们测试了一些不同的模型。部分结果如图11所示。对部分细节放大，可以看出Voronoi Diagram在交点处的连接关系，如图 12所示。但在部分区域，可以看出部分Voronoi Diagram边上有断点，如图所示，可能的原因是对浮点数误差的处理不够好， 使得部分周围关键点数目不为0 的三角形在误差存在的情况下没能正确找到关键点，从而在最终连接时出现问题。

Figure 11: 部分运行结果

Figure 12: 细节放大

图13是在ϵ = 0.5下程序性能测试，在此参数下，所给数据均可在可接受的时间内给出相应结果。测试机器的参数为：CPU: Intel(R) Core(TM) i7-2600; OS: Windows 7 64bit; Memory: 8G；图14是运行时间与O(n) 与O(nlogn) 的比较，其中程序运行的时间放大至10000倍。由图 大体可以看出，无论是原有的MMP算法还是ϵ-Voronoi Diagram算法，在测试的例子中均没有达到时间上界。由[3]中的测试可以看到， MMP算法的实际运行时间大约是低于2次的。

Figure 13: 运行时间测试1

Figure 14: 运行时间测试2

对于ϵ-Vornoi Diagram算法，一个简单的时间复杂度分析如下：
设最长边的长度为D，则易知对三角形的划分轮数不会超过Titer = 3log(D/ϵ)。设第i轮细分后三角形的数目为Ti，由[1]中知每条边上的窗口最坏情况下可达O(n)，易知更新三角形及其Window结构的时间为O(nTi)；另一方面，搜索一条边上的关键点的时间也是O(nTi)。 考虑最坏的情况，即Voronoi Diagram的site数目多至Voronoi Diagram的边经过的三角面的数量为O(n)，并且每一次划分后Voronoi Diagram的边仍经过所有的三角形。则三角形的数目呈指数增长。我们有运行时间：
T(n) = ∑i=0s O(nTi) = n∑i=0s O(Ti) = n∑i=0s (n ×2i) = n2∑i=0s (2i) = n2(2s−1) = (ϵ/D)−3n2
其中s = Titer。可以看出，实际测试数据上，这一界也远没有达到。一个主要原因是一般不会选取使Voronoi Diagram边经过 O(n)的三角面。

由于MMP算法的空间复杂度为O(n2)，而且经过实际实验，MMP算法本身对空间的消耗是非常大的。在对三角网格进行细分后，易知空间复杂度有类似于时间复杂度的表达式，因而在实际测试中，一方面较大的三角网格会占用相当大的空间，另一方面随着ϵ的减小，程序的空间占用会急剧上升，这一算法实现的空间消耗是一个不足之处。

8  实验总结及未来工作
在本次实验中，我们在MMP算法的基础上实现了[5]中的ϵ-Voronoi Diagram算法，该算法是保证拓扑连接关系并能控制误差 的基于测地距离的Voronoi Diagram构造算法。通过测试，我们对算法的理论性能和实际性能做了一定的分析和比较。我们实现了一个交互 友好的展示程序，不但可以展示了实现的ϵ-Voronoi Diagram算法，还可以通过嵌入的MMP算法在一定程度上证明ϵ-Voronoi Diagram算法的正确性。

可以看到，我们的实现仍有待改进：程序中对浮点误差的控制不够完善是一个待改进的重点；而由于算法本身的空间复杂度，其实际空间消耗也是一个值得进一步探讨的问题。由于这里构造的Voronoi Diagram 没有类似于DCEL的连接关系，所以未来可以在组织Voronoi Diagram结构方面展 开进一步的工作；另外，Voronoi Diagram与Medial Axis的关系也是[5]中的一个重要内容，所以今后在构建三角网格上基于测地距离的 Medial Axis方面也有很多工作可以开展。

References

[1]

Mitchell, Joseph SB, David M. Mount, and Christos H. Papadimitriou. "The discrete geodesic problem." SIAM Journal on Computing 16.4 (1987): 647-668.

[2]

Chen, Jindong, and Yijie Han. "Shortest paths on a polyhedron." Proceedings of the sixth annual symposium on Computational geometry. ACM, 1990.

[3]

Surazhsky, Vitaly, et al. "Fast exact and approximate geodesics on meshes." ACM Transactions on Graphics (TOG). Vol. 24. No. 3. ACM, 2005.

[4]

Xin, Shi-Qing, and Guo-Jin Wang. Ïmproving Chen and Han's algorithm on the discrete geodesic problem." ACM Transactions on Graphics (TOG) 28.4 (2009): 104.

[5]

Zhang, Wenqi. "Topology-consistent Voronoi Diagram on Triangulated Mesh and Application" Thesis

File translated from TEX by TTH, version 4.03.
On 15 Jun 2013, 22:45.