From: 已由 Internet Explorer 11 保存
Date: Wed, 29 Jan 2014 21:45:28 +0800

球面Delaunay三角剖分与聚类应用

――计算几何大作业报告

王运涛
2011310600

程  宽
2011210871

赵宇航
2011210800

王  欢
2010310535

下载链接：可执行文件下载 实验报告下载
源文件下载 输入数据预处理模块源码下载（根据夜视图图像像素提取亮点）

以下为简化版的文档说明：

一．研究概述... 3

二．Delaunay三角剖分... 5

2.1 二维Delaunay三角剖分... 6

2.2 球面Delaunay三角剖分... 6

2.2.1
ToLeft测试... 7

2.2.2
InCircle测试... 7

2.2.3
InTriangle测试... 7

2.2.4
构造初始三角形... 8

2.2.5
EdgeFlip过程... 8

2.2.6
DCEL数据结构辅助实现... 8

2.2.7
交互部分... 10

2.2.8球面Delaunay三角剖分复杂度分析... 15

三．基于球面Delaunay三角剖分的聚类... 16

3.1 核心思想... 16

3.2 改进的球面AutoClust算法... 16

3.3 相关代码分析... 21

3.4 复杂度分析... 21

四．球面上聚类的最小覆盖圆计算... 21

4.1 算法主框架... 21

4.2 球面最小覆盖圆算法细节补充... 23

4.3 相关代码分析... 23

4.4 复杂度分析... 23

五．球面上的凸包计算... 24

5.1 算法步骤... 24

5.2 球面凸包算法细节补充... 24

5.3 相关代码分析... 25

5.4 复杂度分析... 25

六．聚类的评价函数... 25

七．应用实例分析... 25

7.1 数据获取... 26

7.1.1
数据来源... 26

7.1.2
数据处理... 26

7.2 夜景图数据展示效果... 27

7.3 地震数据展示效果... 28

八．参考资料... 29

一．研究概述

首先用一个列表显示我们小组所做的工作：

（1）
实现了球面Delaunay三角剖分算法。可以单步执行，观察算法的运行过程；也可以自动执行、一次性执行；

（2）
实现了基于球面Delaunay三角剖分的聚类算法――AutoClust算法，并可以单步执行，观察算法的运行过程；也可以自动执行、一次性执行；

（3）
实现了球面最小覆盖圆算法；

（4）
实现了球面凸包算法――Jarvis
marching算法；

（5）
通过给点分配权值，实现了多分辨率Delaunay三角剖分和多分辨率聚类展示。考虑到功能的相互集成，多分辨率模式需要在启动自动运行模式下，且自动模式时间间隔足够小时，才能流畅的操纵多分辨率。此外，对分辨率级数的设定，也会影响到多分辨率的观看效果。

（6）
使用OpenGL绘图进行展示，用OpenGL实现了画球、画球面三角形、画球面弧、画球面圆、画三维管子，以及光照显示等。

（7）
使用Qt进行界面控制。

（8）
实现了对全球5级以上地震数据的Delaunay三角剖分、聚类实例展示

（9）
实现了对全球夜景图数据的提取，并可以进行Delaunay三角剖分、聚类实例展示。夜景图数据的提取通过分析不同像素点间关联关系的方式实现，用到了OpenCV库。

使用我们的系统需要注意：

（1）我们为读者提供了一定测试实例，输入的三维数据要求的分布于单位球面上。

（2）由于要画大量的管子、球面，所以对电脑显卡要求相对较高。在i5-2300
[email-masked],3.49GB内存，NVIDIA
GEForce GTS 450显卡配置下，操作2000个点在几秒内即可完成，操作20000个点需要几分钟时间。当需要多分辨率实时显示时，操作2000个点大约是其上限。

下面进行详细阐述。

网格划分是有限元分析前处理中的一个关键步骤，其密度、质量等对有限元计算的精度、效率以及收敛性等都有着重要的影响。研究者们已经提出了各种各样的网格生成的方法，其中Delaunay三角剖分(DT)则是最通用的全自动网格生成的方法之一。

目前，Delaunay三角剖分已经不仅仅局限在二维空间中，三维空间中的Delaunay三角剖分在很多实际应用中起到了关键的作用。本文在平面上的Delaunay三角剖分方法的基础上，对三维空间中的球面三角剖分方法进行了探究和实现。同时，我们还提出了一些基于球面Delaunay三角剖分的应用，例如地球夜视图分析、世界地震区域分析等。为了实现这些应用，我们在球面Delaunay三角剖分的基础上，对边进行聚类，以达到对全球经济发达程度(通过对夜视图的分析)和地震多发程度的侦测。

为了证明我们所使用方法的合理性，我们设计了评价函数，对聚类分析的结果进行评价。在聚类效果展示中，我们使用了球面上的最小外接圆算法和凸包算法完成聚类效果展示。

本实验中，我们在Visual
Studio 2008的环境下，用C++语言实现。在实现过程中还用到了OpenGL，OpenCV，Qt。

对地震数据的Delaunay三角剖分和聚类效果如下图所示：

随机生成2000个点后的三角剖分和聚类效果如下：

在下文中，我们将对具体的算法和实现进行详细的介绍。第二部分为Delaunay三角剖分从二维扩展到三维球面的算法；第三部介绍基于球面Delaunay三角剖分的聚类算法；第四、五、六部分分别介绍最小覆盖圆算法、球面凸包算法、聚类评价函数的定义。针对这些算法的实际应用在第七部分进行展示；最后进行未来工作的探究，并对整个工作进行总结。

二．Delaunay三角剖分

本实验的主要内容是将Delaunay三角剖分从二维平面扩展到三维球面。因此，在进行球面Delaunay三角剖分算法的介绍前，我们首先对平面上的Delaunay三角剖分进行简要的分析和实现。

2.1
二维Delaunay三角剖分

该部分详细内容参见文档。

2.2
球面Delaunay三角剖分

球面三角剖分和平面二维三角剖分不同，在球面三角剖分之前，需要定义球面的几个概念。包括球面顶点、 球面边、球面三角形、ToLeft测试以及InCircle测试。

球面顶点定义为在空间中距离原点位置为1的点，即：(x , y , z)，其中x2+y2+z2=1。

球面边定义如下，对于任意两个不重合的球面顶点A、B，如果他们不与原点O共线，则AOB可以确定唯一的一个平面。AOB平面与球面相交，构成一个空间的圆圈C，则圆圈C上以A、B之间的劣弧 为球面边AB。

球面三角形定义如下，球面上不共边(球面边)的三个顶点A、B、C，如果球面边AB、AC和BC都有定义，则将球面变AB、BC、CA依次连接组成的封闭图形即是球面三角形。

ToLeft测试，对于球面边AB和任意一个不在AB边上的点X，规定向量 为球面边的法向量方向。定义球面边AB的左侧为法向量指向的方向。对于点X，判断ToLeft测试成功的方法是 ，如果在线段上则有 ，在右侧有 。

InCircle测试，球面上的不共边(球面边)的三个顶点A、B、C，依次连接球面边可以得到一个球面三角形。而Circle的定义为ABC所在平面与球面的交点组成的圆圈。对于InCircle的判断，需要使用语ToLeft测试中相同的方法，使用定义的法向量来确定球面上一点X是否在Circle中。定义 为面ABC的法向量，InCircle的点的定义为 ，同时，在ABC组成的Circle外部的条件是 ，X在平面ABC上的条件是 。

2.2.1
ToLeft测试

对于球面边AB和球面上任意一个点V，坐标如下：

, ,

定义AOB平面的法向量为 ，则ToLeft函数定义如下：

计算法向量 =，得

如果R > 0，ToLeft测试成功，R = 0，表示V在球面边AB上，R < 0，表示V在AB右侧。

2.2.2
InCircle测试

对于球面边任意不共边(球面边)的三个点A、B、C和球面上任意一个点V，坐标如下：

, ,
，

定义向量如下：

令平面ABC法向量 = ，

得 ，如果R > 0，ToLeft测试成功，R = 0，表示V在球面边AB上，R < 0，表示V在AB右侧。

2.2.3
InTriangle测试

对于球面边任意不共边(球面边)的三个点A、B、C和球面上任意一个点V。如果对球面边AB、AC、BC分别使用ToLeft测试对V进行测试，如果全部测试成功则，表示V点在ABC三角形内部。

2.2.4
构造初始三角形

通过构造4个顶点，12条边，4个三角形来完成初始的Delaunay三角剖分。其中一个顶点的坐标是D(0,0,1)，另外三个顶点ABC和点D构成正四面体，该正四面体的外接球恰好是以0为球心的单位球。

执行外点、边、面构造后，以后没加入一个顶点，就要进行一个rebucket操作，将这些顶点放入对应的三角形内部。

2.2.5
EdgeFlip过程

假设p点对ab边产生怀疑，ab边对应的顶点x落在三角形pab外接圆的内部，那么此时需要进行EdgeFlip操作。EdgeFlip操作的主要步骤如下：

（1）
增加两条边px和xp，并更新这两条边对应的成员信息；

（2）
增加两个三角形pax和pxb，并更新这两个三角形对应的成员信息；

（3）
删除两条边ba和ba；

（4）
修改顶点a和顶点b的成员信息。

（5）
对原三角形pab和abx内部的所有点进行rebucket操作，将其放入三角形pax和三角形pxb内部；

（6）
把边ax和bx放入怀疑边队列内。

2.2.6
DCEL数据结构辅助实现

我们使用DCEL数据结构辅助实现Delaunay三角剖分。SphereVertex表示顶点类，除去记录顶点坐标和顶点权重外，还记录从它出发的任意一条半边。

class
SphereVertex

{

double x,y,z,weight;/*单位球上的顶点及点的权重*/

SphereHalfEdge
*IncidentEdge;/*从该点出发的任意边*/

}

SphereHalfEdge类表示半边，它需要记录该边的起点，它的兄弟边，该边所在的三角形，以及该边的前一条和后一条边。

class
SphereHalfEdge

{

SphereVertex
*Origin;/*边起点*/

SphereHalfEdge
*Twin;/*兄弟边*/

SphereTriangle
*IncidentTriangle;/*边所在三角形*/

SphereHalfEdge
*Next;/*下一条边*/

SphereHalfEdge
*Prev;/*前一条边*/

}

SphereTriangle类表示三角形，它记录该三角形内部的任意一条边，三角形外部边作为保留变量使用，它在程序中不发挥任何作用。此外，还需要有一个与该三角形对应的一个桶，记录三角形内部的剩余点数；若三角形内部没有点，则该桶为空。

class
SphereTriangle

{

SphereHalfEdge
*OuterComponent;/*保留变量，暂时没用*/

SphereHalfEdge
*InnerComponent;/*三角形内部的任意边*/

SphereBucket
*CorrBucket;/*与三角形对应的桶。若等于，表示三角形内部没有剩余点*/

}

SphereBucket类表示桶类，它需要记录一个与之对应的三角形（该三角形不能为空），同时包含桶内部的所有点。

class
SphereBucket

{

public:

std::vector<SphereVertex*>
IncludedVertexVecor;/*桶内部的点，用vector存储*/

SphereTriangle
*CorrTriangle;/*桶对应的三角形，不可以是*/

}

我们用SphereGraph类来组织所有的顶点、半边、三角形。所有的顶点、半边、三角形都用Vector向量组织。同时有一个状态变量SphereGraphState表示当前Delaunay三角剖分状态。

class
SphereGraph

{

public:

vector<SphereVertex*> VertexVector;/*Delaunay顶点*/

vector<SphereHalfEdge*> HalfEdgeVector;/*Delaunay边*/

vector<SphereTriangle*> TriangleVector;/*Delaunay三角形*/

vector<SphereBucket*> BucketVector;/*桶，每个Delaunay三角形对应一个*/

SphereGraphState
CurrentState;/*当前剖分状态*/

SphereGraphState
NextState;/*下一个剖分状态*/

}

2.2.7
交互部分

Ø
预处理

与二维平面上的情形相比较，三维球面的情形较为复杂：任意给定两个顶点
A、B，不一定存在合理的球面边
AB（A、B可能与原点
O共线
）。同样的，任意给定三个顶点
A、B、C，也不一定存在球面三角形，也不一定存在球面三角形
ABC。但可以通过在初始化时，在球体上插入一个“
正四面体
”作为初始三角形集合，来避免这些情况的出现。对于每个球面三角形，建立一个空的桶，用来存储点云数据。

Ø
数据导入

这部分可以实现数据的导入，数据格式为（x
, y , z , area）最后的area可以缺省。数据导入后，首先会将这些数据按照预处理部分生成的三个球面三角形，将所有的数据分布到每个三角形指向的桶中。

Ø
三角剖分状态转移

以下是单步执行的状态机，如下按照每个状态进行说明。整个状态的流程如下图所示：

Begin：最开始进入的状态是Begin，此状态下，完成了数据导入和预处理。

Ready：如果事件包括全部执行，单步执行下一步等驱动了状态转移，判断桶中是否有点，如果所有的桶中都没有未处理的点，则Delaunay三角剖分结束。否则，进入NewVertexSelected状态，在每一次Delaunay三角剖分子过程中，Ready状态是起始状态。界面上的输出结果和输出信息：

Begin：表示只插入了点，还未进行任何剖分，如下图所示。

Ready：表示已完成部分点的剖分，还有一部分点存在桶中，未进行剖分，如下图所示。

NewVertexSelected：表示已经随机的在一个桶中选择了一个顶点，是此次Delaunay三角剖分子过程的处理点，如下图所示。

TriangleFound：表示已经找到了新选择的顶点所在的三角形,
如下图所示。

VertexInserted：表示已经将新选择的顶点插入到已有的三角剖分中，生成了三条新的边，将这三条边加入到DCEL结构中,
如下图所示。

SwapReady：表示准备进行下一步SwapTest，是一个内部循环的起始状态，多次回到这个状态。对于新加入的点，选择一条需要InCircle测试的边，类似于在2.1章节中介绍的方法,
如下图所示。

InCircleTested：选择一条边进行InCircle测试，在这一个状态下，需要判断需不需要边翻转，如果需要进入EdgeFlipped状态，不需要则返回到SwapReady状态进行下一次判断,
如上图第（2,3）个子图所示。

EdgeFlipped：如果InCircle测试成功，需要进行内部的边翻转，是一个对于DCEL结构的操作,
如下图所示,红色表示待删除的边，蓝色表示即将增加的边。

SwapDone：表示所有的SwapTest结束，一个顶点的插入完成，一次Delaunay剖分的子过程完成。在这个状态下，如果所有桶中的点全部处理完，进入DelaunayDone状态，表示所有的剖分完成，否则进入Ready状态，进行下一次子过程。同时维护桶结构，将点重新分配到新生成的三角形对应的桶中。该状态下的绘图如下图所示：

2.2.8球面Delaunay三角剖分复杂度分析

对于DCEL结构的维护上，添加、修改和删除操作都是O(1)的时间复杂度，对于桶的搜索操作为O(h)，其中，每次h是当前存在的三角形的数量。对于rebucketing操作，每次rebucketing时间复杂度是O(logn)，类比于二维平面的Delaunay三角剖分，rebucketing总共需要的时间为3nlogn+O(1)，所以在数据结构层面上，整个算法的时间复杂度为O(nlogn)。

对于Delaunay三角剖分，新加入一个点，需要三次插入操作，如果在InCircle测试中成功，则需要边翻转，需要删除一条边同时增加一条边。类比于二维Delaunay三角剖分，平均需要6次修改操作。所以在这方面，时间复杂度为O(n)。

综上所述，整个球面Delaunay三角剖分的时间复杂度为O(nlogn)。

三．基于球面Delaunay三角剖分的聚类

在三维球面上进行三角剖分后，为了进一步在实际中进行应用，我们希望基于剖分的结果进行聚类，下面对聚类算法进行介绍。

3.1
核心思想

Delaunay三角剖分算法可以用于聚类。聚类算法的思想是把距离较近的点聚成一类，而距离较远的点视为属于不同的类。于是，在Delaunay三角剖分结果的基础上，可以认为较长的边连接的是不同类的点，而较短的边连接的是同一个类中的点。所以，可以把图中较长的边删除，之后求得所有的连通分支，每个连通分支就是一个聚类。我们在传统的AutoClust算法上进行改进，加入调节因子，从而可以人为地控制聚类。传统的AutoClust是针对平面点集的聚类算法，这里将改进为球面算法。

3.2
改进的球面AutoClust算法

Ø
第一步

在已经取得Delaunay三角剖分DD(P)的基础上进行。这里DD(P)表示点集P的三角剖分。

计算与DD(P)中顶点和边相关的若干变量。包括：各个顶点度数，各边的弧长。

计算各顶点的Local_Mean。

.

上式计算
的Local_Mean，
表示
的邻接边集合。
表示
的度数。
含义是
邻接边的平均长度。
表示边
的弧长。

之后计算各顶点的
，代表
邻接边长度的标准差。

.

随后计算
，表示各顶点
的均值。

.

我们根据上面这些量来判别长边和短边。

这里我们引入了调节因子scaleratio，范围在-1到1之间，从而使我们能够调节聚类的精细程度。算法开始时，需要事先输入调节因子scaleratio。最初的AutoClust聚类算法不含参数。在我们的算法中，只需将scaleratio设置为1，就是原本的AutoClust聚类算法。

通过计算，我们将所有点对应的边分为三类：长边、短边、其它边，如下图所示，红色表示长边，绿色表示其它边，蓝色表示短边。注意到一条边对于一个顶点可能是长边，而对于另一个顶点则可能是短边。

下面删除所有的长边。长边的定义如下。其中
。

然后再删除所有的短边。短边定义如下。

如下图所示，删除短边和长边后的结果如下：

Ø
第二步

删除长边和短边之后，剩下的边称之为其它边。此时，求出所有连通分支。下面的操作基于这样一个观察。如果一个顶点拥有其它边，那么所有其它边把该顶点与同一个连通分支相连接。下面是算法的详细步骤。

1.
对于每个顶点
，我们选择一个连通分支C，选择的规则如下：

a)
如果与
邻接的两条短边分别连接两个不同的连通分支，则取较大（顶点数较多）的那个连通分支。

b)
如果有多个连通分支含有的顶点数相同，每个连通分支都由一条短边与
连接，那么我们取对应的短边最短的那个连通分支。

2.
如果
的其它边连接到一个连通分支，该连通分支不是C，那么删除这条其
它边。

3.
将与
邻接的短边都加入到C中。

本步骤可以删除如上图所示两聚类之间连通的情况。就是可以把这种“桥”断开。恢复短边，删除部分中间边后的状态如下：

Ø
第三步

顶点
的k邻接边集合：从
出发的所有小于等于k条边的路径，这些路径上的边构成顶点
的k邻接边集合。
表示顶点
的k邻接边集合。

.

对于每个顶点
，删除所有满足下面条件的边。

such that .

下图中展示了以上三步操作的大致效果。b、c、d分别是经过第一步、第二步、第三步之后的结果。

以上三步是原始的AutoClust算法。但是，经我们分析，该算法存在不完善之处。第二步中，对于只有短边的顶点的处理存在问题。假设有m个顶点，它们之间只有短边相连，并且构成一个连通分支。如下图所示，删除二度距离意义上的长边后，如下图所示：

3.3
相关代码分析

该部分详细内容参见文档。

3.4
复杂度分析

本部分算法在已经获得Delaunay三角剖分的基础上完成，对于Delaunay三角剖分获得的边，边数与顶点数成线性关系，所以边数为O(n)。

第一步，在算法上仅仅是对于点集合和边集合的遍历，所以复杂度为O(n)。

第二步，求解连通分支，复杂度O(n)。之后，对所有顶点的操作都是O(1)，所以总时间O(n)。

第三步，对每个顶点进行深度为2的优先搜索，时间复杂度也是O(n)。

所以在已经获得Delaunay三角剖分的基础上，该算法的时间复杂度是O(n)。

四．球面上聚类的最小覆盖圆计算

在已经求得聚类的基础上，为了能够显示聚类效果，并且为了计算每个聚类的评价函数，需要计算每个聚类的最小覆盖圆。因为算法是在球面点集的基础上进行，所以需要对平面的最小覆盖圆算法进行大幅调整和改进。

4.1
算法主框架

算法中所述的圆，都是指球面上的圆，“圆心”在球面上，“半径”以“圆心”到圆边界的弧长计算。所以最小覆盖圆实际是最小覆盖球冠。

1.
设顶点集合是U，共有n个元素。

2.
将U中的点逐一加入集合V。V在开始时是空集。当加入第i个点，即
时，集合设为
，此时最小覆盖圆设为
。

3.
如果
在
内部或者边界上，说明
；否则
在
外部，那么
就必然在
的边界上，此时需要在已知
上一个顶点
的情况下计算
。

4.

即最小覆盖圆。

在已知
上一个顶点
的情况下计算
，当前点集是
。

1.
将当前最小覆盖圆是
，圆心设为
，半径为0。

2.
遍历
中除了
之外的所有顶点，遍历到第j个顶点，设为
，当前最小覆盖圆是
，强制要求
必须过
。

3.
对于顶点
，如果
在
内部或者边界上，那么
；否则，
在
外部，那么
必然在
边界上，因为
也在
上，所以需要在已知
上两个顶点的情况下计算
。

4.

。

在已知
上两个顶点
和
的情况下计算
，当前点集是
。

1.
将当前最小覆盖圆是
，圆心设为
和
之间球面弧的中点，半径为
和
之间球面弧长度的一半。

2.
遍历
中除了
之外的所有顶点，遍历到第k个顶点，设为
，当前最小覆盖圆是
，强制要求
必须过
和
。

3.
对于顶点
，如果
在
内部或者边界上，那么
；否则，
在
外部，那么
必然在
边界上，因为
和
也在
上，三个顶点即可确定外接圆，故可得
。

4.

。

4.2
球面最小覆盖圆算法细节补充

由于是在球面上求最小覆盖圆。所以一些具体操作与在平面上求最小覆盖圆有很大的不同。

算法中所有的“圆”的圆心都在球面上，而“半径”这里实际是弧。

判断点与圆的位置关系时，计算点到圆心的距离，实际上是点到圆心的弧长。之后再与“半径”作比较。

已知三点求外接圆时，具体方法是求这三点构成平面的法向量，将法向量归一化，就得到“圆心”的坐标。

4.3
相关代码分析

该部分详细内容参见文档。

4.4
复杂度分析

平面最小覆盖圆的时间复杂度是O(n)。球面最小覆盖圆的时间复杂度也是O(n)。分析最小覆盖圆的时间复杂度用到后向分析。

该部分详细内容参见文档。

五．球面上的凸包计算

为了计算聚类的边界。我们改造了凸包算法，使之能够计算球面凸包。球面是个特殊的平面，所以我们可以对平面的凸包算法进行调整，使之能够计算球面凸包。

5.1
算法步骤

1.
找到点集中最下端的点
。即z坐标最小的点。

2.
从顶点
开始，采用与平面凸包算法Jarvis
marching相似的方法，逐一找到边界点。每一个边界点v与它的下一个边界点
以及球心构成一个平面P，剩下的所有点都在P的同侧（左侧）。

3.
转一周之后就可以找到凸包。

5.2
球面凸包算法细节补充

第2步中，如果有顶点落在P上，那么v的下一个边界点
应当是距离v最远的顶点，这里的距离指的是两点之间的弧长。如图，如果A是v，那么
应该是B，而不能是C。

判断顶点在P的哪一侧，需要先计算P的法向量。位于P同侧的顶点与P的法向量的点积同号。

第3步最后，需要判断
是否在
与
所在的圆周上，
是第2
步找到的最后一个顶点。如果
在
与
所在的圆周上，那么删除
。

5.3
相关代码分析

该部分详细内容参见文档。

5.4
复杂度分析

Jarvis
marching算法时间复杂度O(nh)，其中h是凸包包含的顶点数。球面凸包算法对平面凸包算法的一些具体操作进行改动，如上所述。但是这些操作在改动前和改动之后的时间复杂度都是O(1)。所以球面凸包算法的总时间复杂度仍然是O(nh)。

六．聚类的评价函数

该部分详细内容参见文档。

七．应用实例分析

基于我们实现的球面的Delaunay三角剖分和聚类，我们希望将这些算法应用到各种实际问题中去。在本次大作业中，我们选择了两个应用场景：

1.
根据地球上的地震发生情况，对地震发生区域进行划分。

2.
根据地球的夜视图，对地球上经济发达区域进行划分(默认夜景越繁华的地区经济越发达)。

基于以上应用场景，我们针对地球上显示的数据点进行球面Delaunay三角剖分和聚类，从而将地球上的区域按照地震高发情况和经济情况进行划分，再通过最小覆盖圆和球面凸包的算法，将划分出的分类表示出来。在这两个应用场景的实现中，我们将球面上的Delaunay三角剖分、聚类、最小覆盖圆、凸包计算等方法很好的结合，使它们更具有实际应用价值。

7.1
数据获取

7.1.1
数据来源

数据的来源有两部分，地震数据和Google
Earth夜视图数据，其中，地震数据记录了从1556年至2009年震级在5级以上的地震数据。总计187个数据。Google
Earth夜视图数据，使用了1992年和2010年两年的夜视图信息，1992年的数据总计有4393个，2010年的数据总计有3913个。

7.1.2
数据处理

²地震数据

地震作为一种不可以展示的数据，可以直接获得每个数据的经度、纬度和震级等详细的数据，仅仅需要将经度、纬度数据转换为球面空间坐标。

对于经纬度信息 ，映射到球面坐标
为：

；
；

其中 是针对OpenGL纹理映射中的偏移量。

²夜视图数据

Google
Earth夜视图是由美国航天局提供的图片，需要将图片上的数据转为经度、纬度和面积信息，进而使用到球面上。

夜视图数据为一张43201X16801像素的tif格式文件，直接读入会导致内存溢出，所以需要利用PhotoShop将图片分为四份，每份的大小在10000X16801像素，利用金字塔方法，将图片缩小4倍，分辨率变为2500X4200。对获得的图像使用连通域检测，求得每个连通域的中心位置坐标和面积信息，转为经度、纬度信息，进而转为单位球面上的三维点信息。

该数据的详细处理过程参见文档。

在数据处理结束后，将寻找到的连通域映射到背景图上，得到以下效果：

为了使夜景图更好的作为纹理显示，我们对其所有像素取反后的结果如下，这样看起来会更清晰。

7.2
夜景图数据展示效果

如下图所示，夜景图数据点很好的和实际地理位置匹配。由于每个数据点都给了权重，所以可以很方便的进行多分辨率展示。

对非洲大陆夜景图的三角剖分及聚类效果图如下：

7.3
地震数据展示效果

数据显示、三角剖分及聚类效果分别如下图所示。

八．参考资料

[1]Mark de Berg, etc. 邓俊辉译.
计算几何算法与应用，清华大学出版社。

[2]邓俊辉.
计算几何课程讲义。

[3].Estivill-Castro, V. and I. Lee. AUTOCLUST: Automatic
clustering via boundary extraction for mining massive point-data sets. 2000:
Citeseer.

[4].Hipp, J.R., Clustering of Large Event Sets Using Delaunay
Tessellations and K-Medoid Optimization. 2007, DTIC
Document.

[5].2010.袁文.Delaunay三角剖分随机增量算法的演示性实现及球面扩展.

[6].
http://earthquake.usgs.gov/earthquakes/world/historical_country.php

[7].
http://www.ngdc.noaa.gov/dmsp/downloadV4composites.html

[8].

[9]http://visibleearth.nasa.gov/view.php?id=55167

[10].
http://wenku.baidu.com/view/883001868762caaedd33d4bc.html