Delaunay三角剖分.doc

Delauney三角网剖分算法原理：分为三步：一、凸包生成：1）求出如下四点：min(x-y)、min(x+y)、max(x-y)、max(x+y)并顺次放入一个数组，组成初始凸包；2）对于凸包上的点I，设它的后续点为J，计算矢量线段IJ右侧的所有点到IJ的距离，求出距离最大的点K；3）将K插入I,J之间，并将K赋给J；4）重复2，3步，直到点集中没有在IJ右侧的点为止；5）将J赋给I，J取其后续点，重复2，3，4步，当遍历了一次凸包后，凸包生成完成。二、环切边界法凸包三角剖分：在凸包数组中，每次寻找一个由相邻两条凸包边组成的三角形，在该三角形的内部和边界上都不包含凸包上的任何其他点，然后去掉该点得到新的凸包链表，重复这个过程，最终对凸包数组中的点进行三角剖分成功。三、离散的内插：1）建立三角形的外接圆，找出外接圆包含待插入点的所有三角形，构成插入区域；2）删除插入区域内的三角形公共边，形成由影响三角形顶点构成的多边形；3）将插入点与多边形所有顶点相连，构成新的Delaunay三角形；4）重复1，2，3，直到所有非凸包上的离散点都插入完为止。功能实现流程：1. 在绘图菜单栏下添加一个子菜单项为Delauney，并且在工具栏上添加一个工具项。设置text为Delaunay三角剖分，name为delaunay等属性，添加单击事件，并为单击事件代码2.为事件函数添加如下代码 Graphics gra = panel1.CreateGraphics(); List pts = new List(); foreach (Geometry_T geo in choosegeos.Geofeatures) if (geo.GetType() = typeof(Point_T) Point_T pt = (Point_T)geo; pts.Add(pt); List deltins = DelauneyTin(pts);/根据多点构建delauney三角网 foreach (Tin tin in deltins) Point ctin = new Point3; for (int i = 0; i 3; i+) cp = new Point(int)tin.Pthreei.X, (int)tin.Pthreei.Y); ctini = cp; gra.DrawPolygon(Pens.Red, ctin); 3三角形TIN的数据结构public class Tin Point_T pthree = new Point_T3; Line_T lthree = new Line_T3; public Line_T Lthree get return lthree; set lthree = value; public Point_T Pthree get return pthree; set pthree = value; public Tin() public Tin(Point_T p1, Point_T p2, Point_T p3) pthree0 = p1; pthree1 = p2; pthree2 = p3; lthree0 = new Line_T(p1, p2); lthree1 = new Line_T(p2, p3); lthree2 = new Line_T(p3, p1); 4.圆的数据结构 public class Circle_T:Geometry_T private Point_T cpt; public Point_T Cpt get return cpt; set cpt = value; double radius; public double Radius get return radius; set radius = value; public Circle_T() public Circle_T(Point_T pt, double r) cpt = pt; radius = r; 5.实现Delaunay三角剖分算法1） public List DelauneyTin(List pts)/根据多点构建delauney三角网；分三步：构建凸包；凸包剖分；离散点内插 Graphics gra = panel1.CreateGraphics(); List deltins = new List(); List envpts = EnvelopeTin(pts);/构建凸包 /for (int i = 0; i envpts.Count - 1; i+) / / gra.DrawLine(Pens.Black, new Point(int)envptsi.X, (int)envptsi.Y), new Point(int)envptsi + 1.X, (int)envptsi + 1.Y); / /gra.DrawLine(Pens.Black, new Point(int)envpts0.X, (int)envpts0.Y), new Point(int)envptsenvpts.Count - 1.X, (int)envptsenvpts.Count - 1.Y); List dispts = new List();/非凸包上的离散点 foreach (Point_T pt in pts) if (!envpts.Contains(pt) dispts.Add(pt); List envtins = EnvelopeDivision(envpts);/凸包剖分 /foreach (Tin tin in envtins) / / Point ctin = new Point3; / for (int i = 0; i 3; i+) / / cp = new Point(int)tin.Pthreei.X, (int)tin.Pthreei.Y); / ctini = cp; / / gra.DrawPolygon(Pens.Blue, ctin); / deltins = TinInsert(envtins, dispts);/离散点内插 return deltins; 2）public List EnvelopeTin(List pts)/构建凸包 List envpts = new List(); List othpts = new List(); foreach (Point_T pt in pts) othpts.Add(pt); /构建以x-y，x+y最大最小值组成的初始矩形框 CompareXaddY comxandy = new CompareXaddY(); CompareXsubY comxsuby = new CompareXsubY(); pts.Sort(comxsuby); envpts.Add(pts0); envpts.Add(ptspts.Count - 1); othpts.Remove(pts0); othpts.Remove(ptspts.Count-1); pts.Sort(comxandy); if(!envpts.Contains(pts0) envpts.Insert(1, pts0); if (!envpts.Contains(ptspts.Count - 1) envpts.Add(ptspts.Count - 1); othpts.Remove(pts0); othpts.Remove(ptspts.Count-1); /构建以x-y，x+y最大最小值组成的初始矩形框 int i = 0; int tag = 0; bool over = true; while(ienvpts.Count) Line_T cline; if (i=envpts.Count-1) cline = new Line_T(envptsi, envpts0); else cline = new Line_T(envptsi, envptsi + 1); double dismax=0; for (int j = 0; j dismax) dismax = distance; tag = j; over = false; if (over) i+; else /envpts.RemoveAt(i); envpts.Insert(i+1, othptstag); over = true; return envpts; public List EnvelopeDivision(List pts)/凸包剖分 List envtins = new List(); List cpts = new List(); foreach (Point_T pt in pts) cpts.Add(pt); while (cpts.Count 2) int tag = 0; double minangle = 120; for (int i = 0; i cpts.Count; i+) double angle; if (i = 0) angle = CalcuAngle(cptscpts.Count - 1, cptsi, cptsi + 1); else if (i = cpts.Count - 1) angle = CalcuAngle(cptsi-1, cptsi, cpts0); else angle = CalcuAngle(cptsi-1, cptsi, cptsi + 1); if (angle - 60) minangle) minangle = angle - 60; tag = i; int btag=tag-1; int atag=tag+1; if (tag = 0) btag = cpts.Count - 1; else if (tag = cpts.Count - 1) atag = 0; Tin ctin = new Tin(cptsbtag, cptstag, cptsatag); envtins.Add(ctin); cpts.RemoveAt(tag); return envtins; public List TinInsert(List tins, List pts)/离散点内插 List deltins = new List(); List ctins = new List();/临时凸包 foreach (Tin tin in tins) ctins.Add(tin); foreach (Point_T pt in pts)/对离散点遍历，内插 List cpts = new List();/临时点集 foreach (Tin tin in ctins)/找到外接圆包含离散点的三角形 Circle_T ccir = DelauneyCicle(tin);/构造外接圆 if (IsPointInCircle(pt, ccir)/点是否包含在圆内 /for (int i = 0; i 3; i+) / / if (!cpts.Contains(tin.Pthreei) / / cpts.Add(tin.Pthreei);/记录当前点 / / deltins.Add(tin); /记录保存当前三角形 /List ecpts = EnvelopeTin(cpts);/求点集（外接圆包含离散的三角形）的凸包？，接下来，插入点，构建新三角网 /for (int j = 0; j ecpts.Count;j+ ) / / Tin tin; / if (j = ecpts.Count-1) / / tin = new Tin(ecptsj, ecpts0, pt); / / else / / tin=new Tin(ecptsj,ecptsj+1,pt); / / ctins.Add(tin); / List eli = BorderTin(deltins); foreach (Line_T line in eli) Tin tin = new Tin(line.Frompt, line.Topt, pt); ctins.Add(tin); foreach (Tin tin in deltins)/改变临时三角网（删除deltins保存的三角网） ctins.Remove(tin); deltins.Clear(); return ctins; 3） public bool IsLeftPoint(Point_T pt, Line_T line)/点在线的左边;叉积大于 bool yes = false; if (pt.X - line.Frompt.X) * line.ParaA + (pt.Y - line.Frompt.Y) * line.ParaB 0) yes = true; return yes; public double CalcuAngle(Point_T fp, Point_T mp, Point_T tp)/首，中，尾三点构成的夹角 double angle = 0; Point_T vector1 = new Point_T(fp.X - mp.X, fp.Y - mp.Y); Point_T vector2 = new Point_T(tp.X - mp.X, tp.Y - mp.Y); angle = Math.Acos(vector1.X * vector2.X + vector1.Y * vector2.Y) / (Math.Sqrt(vector1.X * vector1.X + vector1.Y * vector1.Y) * Math.Sqrt(vector2.X * vector2.X + vector2.Y * vector2.Y); return angle; public Circle_T DelauneyCicle(Tin tin)/构建三角形的外接圆 double x1 = tin.Pthree0.X; double x2 = tin.Pthree1.X; double x3 = tin.Pthree2.X; double y1 = tin.Pthree0.Y; double y2 = tin.Pthree1.Y; double y3 = tin.Pthree2.Y; double x = (y2 - y1) * (y3 * y3 - y1 * y1 + x3 * x3 - x1 * x1) - (y3 - y1) * (y2 * y2 - y1 * y1 + x2 * x2 - x1 * x1) / (2 * (x3 - x1) * (y2 - y1) - 2 * (x2 - x1) * (y3 - y1); double y = (x2 - x1) * (x3 * x3 - x1 * x1 + y3 * y3 - y1 * y1) - (x3 - x1) * (x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1) / (2 * (y3 - y1) * (x2 - x1) - 2 * (y2 - y1) * (x3 - x1); Point_T cpt = new Point_T(x, y); double radius=Math.Sqrt(Math.Pow(x1-x),2)+Math.Pow(y1-y),2); Circle_T cir = new Circle_T(cpt,radius); return cir; public bool IsPointInCircle(Point_T pt, Circle_T cir) if (Math.Sqrt(Math.Pow(pt.X-cir