前言

        最近闲暇之余,无事喜欢研究一些经典的算法,这两天看了Delauny三角网的一些算法,以其最经典的一种为例进行分享,今天又是反对内卷的一天。

什么是Delauny三角网?

        首先说三角网,三角网是由一系列连续三角形构成的网状的平面控制图形,用于描述二、三维物体表面起伏特征的网格状图形,如下图所示:

        所谓的Delauny(德劳内)三角网,就是由Delauny三角形所组成的三角网,称为Delauny三角网,关于Delauny三角网的定义,可以参照连接:Delaunay三角网是什么? - 知乎 (zhihu.com)

        还是那句话,把定义搞懂了,代码都是小事,不想看定义的可以直接跳到代码部分。

Bowyer-Watson

Delauny三角网和Bowyer-Watson的关系

        了解了Delauny三角网是什么,那么什么又是Bowyer-Watson呢?

        新手容易搞混哈,要想弄明白也简单,Delauny三角网只是一种理想模型,一种理论目标,并不是一种算法,并不指向具体过程,就比如:领导给你交给你一项任务,最终要形成什么样子的效果,达到什么目的,中间的过程是什么就靠自己去实现了,最终完成目标即可。Delauny三角网就是一个目标,要构建Delauny三角网,算法有很多种,这里笔者用的是Bowyer-Watson算法。

核心过程

假定已生成了连接若干个顶点的 Delaunay 三角网格加入一个新的节点,找出所有外接圆包含新加入节点的三角形,并将这些三角形删除,形成一个空腔空腔的节点与新加入的节点连接,形成新的 Delaunay 三角形网格不断循环直到遍历完所有点、

这里笔者就只写一个核心过程,具体关于Bowyer-Watson构建Delauny三角网的过程可以参照博客:Bowyer-Watson算法-CSDN博客

效果演示       

网上关于三角网的算法有很多,现成的库也有,但是笔者为啥非得自己去研究代码呢,这可能是闲得慌吧,哈哈,笔者对于任何东西都不只是想停留在表面,学习原理是进步的一种方式,研究算法的原理过程中所获得的心得是直接看书的两倍甚至以上,每次笔者在研究的过程中都会产生新的灵感与思路,只有学习了原理才能在原来的基础之上加入自己的东西,而不是单靠着一些现成功能的调用就觉得事情很简单,调用谁不会呢。

        当然了,我也不提倡什么都研究,如果对自己学习和研究方向有用的东西,可以适当的去研究一下,说不定你们也能获得新的东西。

        话说远了哈,就此打住。

代码

主代码:

///

///

/// 德劳内三角剖分(效率较低,待优化)

///

public class DelaunyTriangulations

{

///

/// 外部点集

///

private List Coordinates { get; set; }

///

/// 超级三角

///

public Triangle SuperTriangle { get; private set; }

public DelaunyTriangulations(List Coordinates)

{

this.Coordinates = Coordinates;

this.SuperTriangle = GetSuperTriangle();

}

///

/// 获取超级三角(这里我生成了一个等边三角)

///

///

private Triangle GetSuperTriangle()

{

//先获取点集的最大外包矩形

IExtent Envelope = new IExtent();//先创建一个空矩形对象

Coordinates = Coordinates.AsEnumerable().OrderBy(pt => pt.Y).ToList();//(Y升序排列)

double Ymin = Coordinates[0].Y, Ymax = Coordinates[Coordinates.Count - 1].Y;//获取最小Y值和最大Y值

Coordinates = Coordinates.AsEnumerable().OrderBy(pt => pt.X).ToList();//(X升序排列)

double Xmin = Coordinates[0].X, Xmax = Coordinates[Coordinates.Count - 1].X;//获取最小X值和最大X值

Envelope.MinPoint = new IPoint(Xmin, Ymin, 0);//赋值最小点(左下角)

Envelope.MaxPoint = new IPoint(Xmax, Ymax, 0);//赋值最大点(右上角)

Envelope.Scale(1.5, 1.5);//防止后续有点集在外部大三角边上的情况

return ISpatialMath.GetEquilateralTriangle(Envelope);//返回矩形的外接大三角(等边三角形)

}

///

/// 构建三角网

///

///

public List Build()

{

if (Coordinates.Count < 3) return new List();//传入的点集少于3个直接返回空三角网

List DelaunyTriangles = new List();//德劳内三角(最终的结果)

List ReSplitTriangles = new List();//临时列表,用于存放需要重新分割的三角形以及保存分割后的三角形

for (int i = 0; i < Coordinates.Count; i++)

{

IPoint curPoint = Coordinates[i];//当前顶点

if (i == 0)

{

ReSplitTriangles.Add(SuperTriangle);//将外部最大三角放入分割列表

ReSplitTriangles = SplitTriangleWith(ReSplitTriangles, curPoint);

DelaunyTriangles.AddRange(ReSplitTriangles);//将分割后的三角形存放至德劳内三角列表

continue;//第一个点先进行分割得到基础的德劳内三角网

}

//判断当前点位于哪些德劳内三角形的外接圆内

//删除所在外接圆的三角形,并用当前点重新分割删除的三角形

//这一步可以借助空间索引筛选

for (int c = 0; c < DelaunyTriangles.Count; c++)

{

Triangle triangle = DelaunyTriangles[c];

if (triangle.Area <= 0)//移除空三角

{

DelaunyTriangles.RemoveAt(c);

c--;

continue;

}

Circle circle = triangle.Tag as Circle;

//当前点到三角形的外接圆圆心距离如果小于等于圆的半径,则认为点在圆内部

if (circle.Center.DistanceTo(curPoint, true) <= circle.Radius)

{

ReSplitTriangles.Add(triangle);//将当前三角纳入重新分割的列表

DelaunyTriangles.RemoveAt(c);//将当前三角从德劳内三角列表移除

c--;

}

}

ReSplitTriangles = SplitTriangleWith(ReSplitTriangles, curPoint);

DelaunyTriangles.AddRange(ReSplitTriangles);//将分割后的三角形存放至德劳内三角列表

if (ReSplitTriangles.Count != 0) ReSplitTriangles.Clear();//清理临时列表

}

FilterDelauny(DelaunyTriangles);//过滤三角网(将包含外部大三角顶点的德劳内三角删除)

return DelaunyTriangles;

}

///

/// 通过点将三角形拆分成多个三角形

///

/// 拆分点

///

private List SplitTriangleWith(List Triangles, IPoint Point)

{

List triangles = new List();

if (Triangles.Count == 0) return triangles;//如果没有待分割的三角形,则返回空值

List Vertex = new List();//先获取三角形的所有顶点

for (int i = 0; i < Triangles.Count; i++)

{

IPoint top = Triangles[i].Top;

IPoint left = Triangles[i].LeftBottom;

IPoint right = Triangles[i].RightBottom;

//如果点集内不存在当前点坐标则加入

var q1 = from pt in Vertex where pt.Equals(top) select pt;

if (q1.Count() == 0)

{

Vertex.Add(top);

}

var q2 = from pt in Vertex where pt.Equals(left) select pt;

if (q2.Count() == 0)

{

Vertex.Add(left);

}

var q3 = from pt in Vertex where pt.Equals(right) select pt;

if (q3.Count() == 0)

{

Vertex.Add(right);

}

}

//将点集按照与拆分点连线的水平方向夹角进行排序(逆时针)

Vertex = Vertex.AsEnumerable().OrderBy(pt => Point.VectorTo(pt).Degree).ToList();

for (int i = 0; i < Vertex.Count; i++)

{

int next = i == Vertex.Count - 1 ? 0 : i + 1;

IPoint left = Vertex[i];

IPoint right = Vertex[next];

Triangle tri = new Triangle(Point, left, right);//创建三角形

CircumCircle circumCircle = new CircumCircle();//调用外接圆工具

tri.Tag = circumCircle.BuildFromTriangle(tri);//绑定外接圆形至三角

triangles.Add(tri);//将分割后的三角形添加到结果

}

Vertex.Clear();//清理临时列表

return triangles;

}

///

/// 过滤德劳内三角(将包含外部大三角顶点的德劳内三角删除)

///

///

private void FilterDelauny(List Delauny)

{

List outPoints = new List();//记录外部大三角顶点

outPoints.Add(SuperTriangle.Top);

outPoints.Add(SuperTriangle.LeftBottom);

outPoints.Add(SuperTriangle.RightBottom);

for (int i = 0; i < Delauny.Count; i++)

{

IPoint top = Delauny[i].Top;

IPoint left = Delauny[i].LeftBottom;

IPoint right = Delauny[i].RightBottom;

var q = from pt in outPoints where pt.Equals(top) || pt.Equals(left) || pt.Equals(right) select pt;

if (q.Count() > 0)

{

Delauny.RemoveAt(i);

i--;

}

}

outPoints.Clear();

}

}

辅助类

Vector

///

/// 向量

///

public class Vector

{

///

/// 向量X方向增量

///

public double X { get; set; }

///

/// 向量Y方向的增量

///

public double Y { get; set; }

///

/// 向量Z方向的增量

///

public double Z { get; set; }

///

/// 向量的模长

///

public double Norm { get { return GetNormOfVector(); } }

public Vector(double X, double Y, double Z = 0)

{

this.X = X;

this.Y = Y;

this.Z = Z;

}

///

/// 计算与其他向量的夹角

///

///

///

public double DegreeTo(Vector OtherVector)

{

// 计算向量长度

double vector1Length = Math.Sqrt(this.X * this.X + this.Y * this.Y + this.Z * this.Z);

double vector2Length = Math.Sqrt(OtherVector.X * OtherVector.X + OtherVector.Y * OtherVector.Y + OtherVector.Z * OtherVector.Z);

// 计算向量点积

double dotProduct = this.X * OtherVector.X + this.Y * OtherVector.Y + this.Z * OtherVector.Z;

// 计算夹角(单位为弧度)

double angle = Math.Acos(dotProduct / (vector1Length * vector2Length));

return angle.RadiansToDegree();

}

///

/// 求向量的模长

///

///

private double GetNormOfVector()

{

return Math.Sqrt(this.X * this.X + this.Y * this.Y + this.Z * this.Z);

}

}

IMath

public static class IMath

{

///

/// 弧度转角度

///

/// 弧度值

///

public static double RadiansToDegree(this double Radians)

{

return (180.00 / Math.PI) * Radians;

}

///

/// 角度转弧度

///

/// 角度值

///

public static double DegreeToRadians(this double Degree)

{

return Degree * Math.PI / 180.00;

}

}

Triangle

public class Triangle

{

///

/// 顶点

///

public IPoint Top { get { return _Top; } set { _Top = value; } }

private IPoint _Top;

///

/// 左底点

///

public IPoint LeftBottom { get { return _LeftBottom; } set { _LeftBottom = value; } }

private IPoint _LeftBottom;

///

/// 右底点

///

public IPoint RightBottom { get { return _RightBottom; } set { _RightBottom = value; } }

private IPoint _RightBottom;

///

/// 左边

///

public ILine LeftSide { get { return new ILine(Top, LeftBottom); } }

///

/// 右边

///

public ILine RightSide { get { return new ILine(Top, RightBottom); } }

///

/// 底边

///

public ILine Bottom { get { return new ILine(LeftBottom, RightBottom); ; } }

///

/// 周长

///

public double Length { get { return LeftSide.Length + RightSide.Length + Bottom.Length; } }

///

/// 面积

///

public double Area { get { return Math.Sqrt((this.Length / 2) * ((this.Length / 2 - RightSide.Length) * (this.Length / 2 - LeftSide.Length) * (this.Length / 2 - Bottom.Length))); } }

///

/// 判断是否有效三角

///

public bool IsValid { get { return !(this.Area == double.NaN || this.Area <= 0); } }

///

/// 外部挂接信息(可自定义对象)

///

public object Tag { get; set; }

///

/// 实例化

///

///

///

///

public Triangle(IPoint Top,IPoint LeftBottom,IPoint RightBottom)

{

this.Top = Top;

this.RightBottom = RightBottom;

this.LeftBottom = LeftBottom;

}

}

IPoint

///

/// 点对象

///

public class IPoint

{

public double X

{

get

{

//无精度限定

if (LeaderConfig.GeoDights==-1)

{

return _x;

}

else

{

return Math.Round(_x, LeaderConfig.GeoDights);

}

}

set

{

_x = value;

}

}

private double _x;

public double Y

{

get

{

if (LeaderConfig.GeoDights==-1)

{

return _y;

}

else

{

return Math.Round(_y, LeaderConfig.GeoDights);

}

}

set

{

_y = value;

}

}

private double _y;

public double Z

{

get

{

if (LeaderConfig.GeoDights==-1)

{

return _z;

}

else

{

return Math.Round(_z, LeaderConfig.GeoDights);

}

}

set

{

_z = value;

}

}

private double _z;

///

/// 实例化

///

///

///

///

public IPoint(double X, double Y, double Z)

{

this.X = X;

this.Y = Y;

this.Z = Z;

}

///

/// 两点距离

///

/// 第二点

/// 是否忽略Z值

///

public double DistanceTo(IPoint SecondPoint,bool IgnoreZvalue=false)

{

if (IgnoreZvalue)

{

return Math.Sqrt(Math.Pow(this.X - SecondPoint.X, 2) + Math.Pow(this.Y - SecondPoint.Y, 2));

}

else

{

return this.VectorTo(SecondPoint).Norm;

}

}

///

/// 获取到另一点的向量

///

///

public Vector VectorTo(IPoint SecondPoint)

{

return new Vector(SecondPoint.X - this.X, SecondPoint.Y - this.Y, SecondPoint.Z - this.Z);

}

///

/// 判断两点是否相等

///

///

///

public bool Equals(IPoint iPoint)

{

return this.X == iPoint.X && this.Y == iPoint.Y && this.Z == iPoint.Z;

}

}

ILine

///

/// 线段

///

public class ILine

{

///

/// 起点

///

public IPoint StartPoint { get { return _startPoint; } set { _startPoint = value; } }

private IPoint _startPoint;

///

/// 终点

///

public IPoint EndPoint { get { return _endPoint; } set { _endPoint = value; } }

private IPoint _endPoint;

public double Length { get { return StartPoint.DistanceTo(EndPoint); } }

///

/// 实例化

///

///

///

public ILine(IPoint start, IPoint end)

{

_startPoint = start;

_endPoint = end;

}

}

CicumCircle

public class CircumCircle

{

///

/// 生成三角形的外接圆

///

///

///

public Circle BuildFromTriangle(Triangle triangle)

{

if (!triangle.IsValid) return null;

IPoint Center = GetCircleCenterWith(triangle.Top, triangle.LeftBottom, triangle.RightBottom);

return new Circle(Center, triangle.Top.DistanceTo(Center,true));

}

///

/// 已知圆上3点获取圆心

///

///

///

///

///

private IPoint GetCircleCenterWith(IPoint P1, IPoint P2, IPoint P3)

{

double a13 = P1.X - P3.X;

double a13_ = P1.X + P3.X;

double b13 = P1.Y - P3.Y;

double b13_ = P1.Y + P3.Y;

double a12 = P1.X - P2.X;

double a12_ = P1.X + P2.X;

double b12 = P1.Y - P2.Y;

double b12_ = P1.Y + P2.Y;

double a12b12_2 = a12 * a12_ + b12 * b12_;

double a13b13_2 = a13 * a13_ + b13 * b13_;

double a13b12 = 2 * a13 * b12;

double a12b13 = 2 * a12 * b13;

if (a12b13 - a13b12 == 0) return new IPoint((P2.X + P1.X) / 2.0, (P2.Y + P1.Y) / 2.0, 0);

double af = a12b13 - a13b12;

double bf = a13b12 - a12b13;

double az = b13 * a12b12_2 - b12 * a13b13_2;

double bz = a13 * a12b12_2 - a12 * a13b13_2;

double a = az / af;

double b = bz / bf;

return new IPoint(a, b, 0);

}

}

Circle

public class Circle:IEntity

{

///

/// 圆心

///

public IPoint Center { get; set; }

///

/// 半径

///

public double Radius { get; set; }

///

/// 面积

///

public double Area { get { return Math.PI * Math.Pow(Radius, 2); } }

///

/// 周长

///

public double Length { get { return Math.PI * Radius * 2; } }

public bool IsValid { get { return Radius > 0; } }

///

/// 点径构圆

///

///

///

public Circle(IPoint Center,double Radius)

{

this.Center = Center;

this.Radius = Radius;

}

///

/// 两点构圆

///

///

///

public Circle(IPoint Center, IPoint PointOnCircle)

{

this.Center = Center;

this.Radius = Center.DistanceTo(PointOnCircle);

}

}

结语

        如果要自己改代码,一定得先看懂BowyerWatson算法的过程,重要的事说三遍,不懂的可以私聊我。

参考链接

评论可见,请评论后查看内容,谢谢!!!
 您阅读本篇文章共花了: