凸包

实验报告

姓名：张晖

班级：地信1002

学号：201001136 指导：刘少华老师

2013.4.17

基于凸包的三角网算法

一，实验背景

数字地形模型是针对地形地貌的一种数字建模,这种建模的结

果通常就是一个数字高程模型（DEM）。不规则三角网（TIN）

模型是DEM中存储和表示非规则数据的理想模型，它既减少

规则网格方法造成的数据冗余，同时在计算效率方面又优于

纯粹基于等高线的方法，所以寻求一种好的TIN算法更能快

速逼真的显示与模拟出地貌三维信息。在所有可能的三角网

中，Delaunay三角网在地形拟合方面表现最为出色，因此常

常用于TIN的生成。

二，实验目的

1，了解TIN的算法种类及各算法特点，尤其是对三角网生长

算法的了解与掌握

2，根据基于凸包的三角网算法，应用C#编写一个自动生成三

角网的程序。

3，学生熟练掌握C#的使用。

三，实验内容

1. 查找关于三角网生成算法的论文，了解与掌握其要义。

2. 分析与研究三角网生成算法

3. 选择一种算法实现三角网的生成（三角网生长算法）

4. 选择一种语言编写三角网的生成代码（C#）

四，算法简介

（1）凸包生成。①求出点集中满足min（x-y）、min（x+y）、max（x-y）、max（x+y）的四个点，并按照逆时针方向组成一个链表。这四个点是离

散点中与包含离散点的外接矩形的四个角点最接近的点，这四个点构成

的多边形作为初始凸包。

②对于每个凸包上的点I，设它的后续点为J，计算矢量线段IJ右侧的所有

点到IJ的距离，求出距离最大的点K。

③将K插入I、J之间，并将K赋给J。

④重复②、③步，直到点集中没有在线段IJ右侧的点为止。

⑤将J赋给I，J取其后续点，重复②、③、④步。

⑥当凸包中任意相邻两点连线的右侧不存在离散点时，结束点集凸包求取

过程。

完成这一步后，形成了包含所有离散点的多边形（凸包）。

（2）环切边界法凸包三角剖分。在凸包链表中每次寻找一个由相邻两条凸包边组成的三角形，在该三角形的内部和边界上都不包含凸包上的任何

其他点。将这个点去掉后得到新的凸包链表。重复这个过程，直到凸包

链表中只剩下三个离散点为止。将凸包链表中的最后三个离散点构成一

个三角形，结束凸包三角剖分过程。

完成这一步后，将凸包中的点构成了若干Delaunay三角形。

（3）离散点内插。在对凸包进行三角剖分之后，不在凸包上的其余离散点，可采用逐点内插的方法进行剖分。其基本过程为：

①找出外接圆包含待插入点的所有三角形，构成插入区域。

②删除插入区域内的三角形公共边，形成由三角形顶点构成的多边形。

③将插入点与多边形所有顶点相连，构成新的Delaunay三角形。

④重复①、②、③，直到所有非凸壳离散点都插入为止。

完成这一步，就完成了Delaunay三角网的构建。

五，源代码

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using System.Drawing;

namespace TINProject

{

class CPoint

{

public int id;

public float x, y;

public CPoint()

{ id = 0; }

}

class CEdge

{

public int id;

public int[] pid;

public int[] tid;

public CEdge()

{

id = 0;

pid = new int[2];

tid = new int[2];

}

}

class CTriangle

{

public int id;

public int[] pid;

public int[] eid;

public CTriangle()

{

id = 0;

pid = new int[3];

eid = new int[2];

}

}

class CTin

{

public List _PntsAry;

public List _EdgesAry;

public List _TrianlgesAry;

public List _TuBaoAry;

public CTin()

{

_PntsAry = new List();

_EdgesAry = new List();

_TrianlgesAry = new List();

_TuBaoAry = new List();

}

public bool IsRightInLine(CPoint p1, CPoint p2, CPoint p)

{

double dx, dy;

dx = p2.x - p1.x;

dy = p2.y - p1.y;

double d = dy * (p.x - p1.x) + dx * (p1.y - p.y);

if (d < 0)

return true;

return false;

}

public double distance(CPoint a, CPoint b

{

double dis;

dis = Math.Sqrt((b.y - a.y) * (b.y - a.y) + (b.x - a.x) * (b.x - a.x));

return dis;

}

public double GetThreePointsCosAngle(CPoint a, CPoint b, CPoint c)

{

double distance1, distance2, distance3;

double cos;

distance1 = Math.Pow((b.x - a.x), 2) + Math.Pow((b.y - a.y), 2);

distance2 = Math.Pow((c.x - a.x), 2) + Math.Pow((c.y - a.y), 2);

distance3 = Math.Pow((b.x - c.x), 2) + Math.Pow((b.y - c.y), 2);

cos = (double)(distance3 + distance2 - distance1) * 1.000000 / (2 * Math.Sqrt(distance3) * Math.Sqrt(distance2));

return cos;

}

public int GetMaxCosAngle(CPoint p1, CPoint p2)

{

double minCos = 1;

int targetIndex = -1;

foreach (CPoint p in_PntsAry)

{

if (IsRightInLine(p1, p2, p))

{

double cosangle = GetThreePointsCosAngle(p1, p2, p);

if (cosangle < minCos)

{

minCos = cosangle;

targetIndex = p.id;

}

}

}

return targetIndex;

}

public int get_leftpoint_cos_index(CPoint p1, CPoint p2)

{

double minCos = 1;

int Index = -1;

foreach (CPoint p in_PntsAry)

{

if (!IsRightInLine(p1, p2, p))

{

double cosangle = GetThreePointsCosAngle(p1, p2, p);

if (cosangle < minCos)

{

minCos = cosangle;

Index = p.id;

}

}

}

return Index;

}

private Pen pen = new Pen(Color.Blue, 5);

private Brush brush = Brushes.Blue;

float r = 10;

public void draw(Graphics g)

{

foreach (CPoint p in_PntsAry)

{

g.FillEllipse(brush, p.x - r / 2, p.y - r / 2, r, r);

}

for (int i = 0; i <= _EdgesAry.Count - 1; i++)

{

g.DrawLine(new Pen(Color.Black), _PntsAry[_EdgesAry[i].pid[0]].x, _PntsAry[_EdgesAry[i].pid[0]].y, _PntsAry[_EdgesAry[i].pid[1]].x,

_PntsAry[_EdgesAry[i].pid[1]].y);

}

for(int j=0;j<_TuBaoAry.Count-1 ;j++)

{

g.DrawLine(new Pen(Color.Red), _TuBaoAry[j].x, _TuBaoAry[j].y, _TuBaoAry[j + 1].x, _TuBaoAry[j + 1].y);

}

}

public void tubao_construct() //建立凸包

{

_TuBaoAry.Clear();

int i, min_index = -1, index = -1;

double min_y, maxangle = -1, angle = -1;

CPoint z=new CPoint();

CPoint p1 = new CPoint();

z = _PntsAry[0];

min_y = z.y;

for (i = 0; i < _PntsAry.Count; i++)

{

if (min_y >= _PntsAry[i].y)

{

min_y = _PntsAry[i].y;

min_index = i;

}

}

p1 = _PntsAry[min_index];

_TuBaoAry.Add(p1);

CPoint p = new CPoint();

p.id = -1;

p.x = p1.x + 200;

for (int j = 0; j < _PntsAry.Count; j++)

{

if (p1.id != _PntsAry[j].id)

{

angle = GetThreePointsCosAngle( p, _PntsAry[j],p1);

if (angle > maxangle)

{

maxangle = angle;

index = j;

}

}

}

if (index >= 0)

{

_TuBaoAry.Add(_PntsAry[index]);

}

int count = _TuBaoAry.Count - 1;

CPoint p0 = _TuBaoAry[count];

while (p0.id != p1.id)

{

count = _TuBaoAry.Count - 1;

CPoint p2 = _TuBaoAry[count - 1];

CPoint temp_point = new CPoint();

temp_point.id =p2.id;

temp_point.x = 2 *p0.x -p2.x;

temp_point.y = 2 *p0.y - p2.y;

maxangle = -1;

index = -1;

for (int j = 0; j < _PntsAry.Count; j++)

{

if (p0.id != _PntsAry[j].id)

{

angle = GetThreePointsCosAngle( temp_point, _PntsAry[j],p0);

if (angle > maxangle)

{

maxangle = angle;

index = j;

}

}

}

if (index >=0)

{

_TuBaoAry.Add(_PntsAry[index]);

count++;

p0 = _TuBaoAry[count];

}

}

}

public void construct_tin()

{

_EdgesAry.Clear();

if (_PntsAry.Count < 3)

int use_point_cout = 0;

int i, count = -1;

CEdge edg = new CEdge();

CEdge one = new CEdge();

CEdge other = new CEdge();

int firstpoint = 0;

double mixdistance = 10000000, dis1;

for (i = 1; i < _PntsAry.Count; i++)

{

dis1 = distance(_PntsAry.ElementAt(firstpoint), _PntsAry.ElementAt(i));

if (mixdistance > dis1)

{

mixdistance = dis1;

count = i;

}

}

if (count < 0) return;

edg.id = 0;

edg.pid[0] = firstpoint;

edg.pid[1] = count;

_EdgesAry.Add(edg);

CPoint p2 = _PntsAry[count];

int k = GetMaxCosAngle(_PntsAry[firstpoint], p2);

if (k >= 0)

{

one.id = 1;

one.pid[0] = firstpoint;

one.pid[1] = k;

_EdgesAry.Add(one);

other.id = 2;

other.pid[0] = k;

other.pid[1] = count;

_EdgesAry.Add(other);

}

int f;

int currentid = 1;

while (currentid <= _EdgesAry.Count - 1)

{

int c, b;

f =

GetMaxCosAngle(_PntsAry.ElementAt(_EdgesAry.ElementAt(currentid).pid[0]),

_PntsAry.ElementAt(_EdgesAry.ElementAt(currentid).pid[1]));

if (f > 0)

{

c = _PntsAry.ElementAt(_EdgesAry.ElementAt(currentid).pid[0]).id;

b = f;

if (!IsExistCedg(c, b))

{

CEdge edg2 = new CEdge();

edg2.id = _EdgesAry.Count;

edg2.pid[0] = c;

edg2.pid[1] = b;

_EdgesAry.Add(edg2);

}

c = f;

b = _PntsAry.ElementAt(_EdgesAry.ElementAt(currentid).pid[1]).id;

if (!IsExistCedg(c, b))

{

CEdge edg1 = new CEdge();

edg1.id = _EdgesAry.Count;

edg1.pid[0] = c;

edg1.pid[1] = b;

_EdgesAry.Add(edg1);

}

}

use_point_cout = (_EdgesAry.Count + 3) / 2;

currentid++;

}

}

public bool IsExistCedg(int i, int j)

{

foreach (CEdge t in_EdgesAry)

{

if (t.pid[0] == i && t.pid[1] == j)

return true;

else if (t.pid[1] == i && t.pid[0] == j)

return true;

}

return false;

}

public void drawline(Graphics g)

{

CPoint p1 = new CPoint();

CPoint p2 = new CPoint();

foreach (CEdge edg in_EdgesAry)

{

for (int i = 0; i < _EdgesAry.Count - 1; i++)

{

for (int j = 0; j < _PntsAry.Count - 1; j++)

{

if (_PntsAry[j].id == _EdgesAry.ElementAt(i).pid[0] &&

_PntsAry[j].id == _EdgesAry.ElementAt(i).pid[1])

{

p1.x = _PntsAry[_EdgesAry.ElementAt(i).pid[0]].x;

p1.y = _PntsAry[_EdgesAry.ElementAt(i).pid[0]].y;

p2.x = _PntsAry[_EdgesAry.ElementAt(i).pid[1]].x;

p2.y = _PntsAry[_EdgesAry.ElementAt(i).pid[1]].y;

g.DrawLine(new Pen(Color.Black), p1.x, p1.y, p2.x, p2.y);

}

}

}

} }}

（1）自己取点

（2）随机取点

在刘老师的指导下，通过一段时间的学习，我对数字高程模型有了一定的了解。TIN（Triangulated Irregular Network）即不规则三角网，Delaunay三角网是其中的一种表现形式，也是一种主要的DTM 的表示法。老师针对几种算法中的生成算法进行了重点讲解和实际操作指导。清楚明白的了解到TIN即不规则三角网的生成原理和过程Delaunay三角网生成算法经过多年的发展已经相当成熟，它在生活、军事、工业等领域都得到了广泛地应用，Delaunay三角网生成算法具有非常重要的意义。