本帖最后由 始于丨初心 于 2020-5-4 13:52 编辑
不规则三角网(TIN, Triangulated Irregular Network)模型采用一系列相连接的三角形拟合地表或其他不规则表面,常用来构造数字地面模型,特别是数字高程模型。
using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
namespace Create_TIN
{
//Tin_Point类
class Tin_Point
{
//X、Y、Z坐标
public double X { get; set; }
public double Y { get; set; }
public double Z { get; set; }
}
//Tin_Line类
class Tin_Line
{
//两端点在离散点集的下标
public int P0 { get; set; }
public int P1 { get; set; }
//线段的使用次数
public int Count { get; set; }
}
//Tin_Tri类
class Tin_Tri
{
//三条线段在基线数据集的下标
public int L0 { get; set; }
public int L1 { get; set; }
public int L2 { get; set; }
}
//TIN_Tin类
class Tin_Tin
{
private List<Tin_Point> points = new List<Tin_Point>();//存储离散点数据
private List<Tin_Line> lines = new List<Tin_Line>();//存储基线数据
private List<Tin_Tri> tris = new List<Tin_Tri>();//存储三角形数据
//获取离散点数据
public bool create_points(string path)
{
StreamReader sr = new StreamReader(path, Encoding.UTF8);
string data = sr.ReadLine();
while (data != null)
{
string[] info = data.Split(',');
try
{
Tin_Point item = new Tin_Point();
item.X = Convert.ToDouble(info[0]);
item.Y = Convert.ToDouble(info[1]);
item.Z = Convert.ToDouble(info[2]);
points.Add(item);
}
catch
{
return false;
}
data = sr.ReadLine();
}
return true;
}
//获取两点之间的平距
private double Distance(Tin_Point a, Tin_Point b)
{
double deltaX = b.X - a.X;
double deltaY = b.Y - a.Y;
double temp = Math.Sqrt(deltaX * deltaX + deltaY * deltaY);
return temp;
}
//获取三角形角A的大小,单位弧度
private double Angle(Tin_Line line, int p)
{
double c = Distance(points[line.P0], points[line.P1]);
double a = Distance(points[line.P0], points[p]);
double b = Distance(points[line.P1], points[p]);
double temp = Math.Acos((a * a + b * b - c * c) / (2 * a * b));
return temp;
}
//根据最大角原则获取第三点
private int Max(Tin_Line line, List<int> usablepoints)
{
double max = Angle(line, usablepoints[0]);
int index = usablepoints[0];
for (int i = 1; i < usablepoints.Count; i++)
{
double temp = Angle(line, usablepoints[i]);
if (max <= temp)
{
max = temp;
index = usablepoints[i];
}
}
return index;
}
//验证找到的第三点是否满足空外接圆准则
private bool Circumcircle(Tin_Line line, int p)
{
Tin_Point a = points[line.P0];
Tin_Point b = points[line.P1];
Tin_Point c = points[p];
double A1 = 2 * (b.X - a.X);
double A2 = 2 * (c.X - b.X);
double B1 = 2 * (b.Y - a.Y);
double B2 = 2 * (c.Y - b.Y);
double C1 = b.X * b.X + b.Y * b.Y - a.X * a.X - a.Y * a.Y;
double C2 = c.X * c.X + c.Y * c.Y - b.X * b.X - b.Y * b.Y;
double X = (B2 * C1 - B1 * C2) / (A1 * B2 - A2 * B1);
double Y = (-A2 * C1 + A1 * C2) / (B2 * A1 - B1 * A2);
Tin_Point centre = new Tin_Point();
centre.X = X;
centre.Y = Y;
double R = Distance(a, centre);
for (int i = 0; i < points.Count; i++)
{
if (line.P0 != i &&
line.P1 != i &&
p != i)
{
double temp = Distance(centre, points[i]);
if (temp < R)
{
return false;
}
}
}
return true;
}
//获取某点在竖直方向到直线的距离
private double Sign(Tin_Line line, int p)
{
Tin_Point a = points[line.P0];
Tin_Point b = points[line.P1];
Tin_Point c = points[p];
double A = (b.Y - a.Y) / (b.X - a.X);
double B = (b.X * a.Y - a.X * b.Y) / (b.X - a.X);
double Z = c.Y - A * c.X - B;
return Z;
}
//判断两条线段是否为同一条
private bool Equaltion(Tin_Line line1, Tin_Line line2)
{
if (line1.P0 == line2.P0 && line1.P1 == line2.P1)
{
return true;
}
if (line1.P0 == line2.P1 && line1.P1 == line2.P0)
{
return true;
}
return false;
}
//创建不规则三角网
public void create_tin()
{
/*---生成第0个三角形---*/
Tin_Tri tri = new Tin_Tri();
/*寻找与第0点距离最近点*/
double mindis = Distance(points[0], points[1]);
int index = 1;
for (int i = 1; i < points.Count; i++)
{
double temp = Distance(points[0], points[i]);
if (temp < mindis)
{
mindis = temp;
index = i;
}
}
/*初始基线L0*/
Tin_Line line0 = new Tin_Line();
line0.P0 = 0;
line0.P1 = index;
line0.Count = 1;
lines.Add(line0);
tri.L0 = lines.Count - 1;
/*生成可用点集*/
List<int> usablepoints = new List<int>();
for (int i = 1; i < points.Count; i++)
{
if (i != index )
{
usablepoints.Add(i);
}
}
/*根据张角最大原则获取第三点*/
index = Max(line0, usablepoints);
/*生成基线L1*/
Tin_Line line1 = new Tin_Line();
line1.P0 = line0.P0;
line1.P1 = index;
line1.Count = 1;
lines.Add(line1);
tri.L1 = lines.Count - 1;
/*生成基线L2*/
Tin_Line line2 = new Tin_Line();
line2.P0 = line0.P1;
line2.P1 = index;
line2.Count = 1;
lines.Add(line2);
tri.L2 = lines.Count - 1;
tris.Add(tri);
/*第0个三角形生成结束*/
/*对于第K个三角形为拓展三角形时*/
int k = 0;
while (true)
{
/*处理第K个三角形的L0边*/
if (lines[tris[k].L0].Count < 2)
{
/*获取第K个三角形,不是L0边上的端点*/
int test_point = -1;
if (lines[tris[k].L0].P0 == lines[tris[k].L1].P0 ||
lines[tris[k].L0].P1 == lines[tris[k].L1].P0)
{
test_point = lines[tris[k].L1].P1;
}
if (lines[tris[k].L0].P0 == lines[tris[k].L1].P1 ||
lines[tris[k].L0].P1 == lines[tris[k].L1].P1)
{
test_point = lines[tris[k].L1].P0;
}
/*获取该端点在竖直方向到L0线段的距离*/
double sign1 = Sign(lines[tris[k].L0], test_point);
/*清空可用点集*/
usablepoints.Clear();
/*获取可用点集*/
for (int i = 0; i < points.Count; i++)
{
if (lines[tris[k].L0].P0 != i &&
lines[tris[k].L0].P1 != i &&
lines[tris[k].L1].P0 != i &&
lines[tris[k].L1].P1 != i &&
lines[tris[k].L2].P0 != i &&
lines[tris[k].L2].P1 != i)
{
double sign2 = Sign(lines[tris[k].L0], i);
if (sign1 * sign2 < 0)
{
usablepoints.Add(i);
}
}
}
/*生成下一个三角形*/
if (usablepoints.Count > 0)
{
/*根据张角最大原则获取第三点*/
index = Max(lines[tris[k].L0], usablepoints);
/*验证第三点是否满足空外接圆准则*/
if (Circumcircle(lines[tris[k].L0], index))
{
/*生成临时L1、L2基线*/
line1 = new Tin_Line();
line1.P0 = lines[tris[k].L0].P0;
line1.P1 = index;
line2 = new Tin_Line();
line2.P0 = lines[tris[k].L0].P1;
line2.P1 = index;
/*标记临时基线是否存在的标志*/
int lineIdx1 = -1;
int lineIdx2 = -1;
int count1 = 0;
int count2 = 0;
/*判断基线数据集中是否存在这两条边*/
for (int i = 0; i < lines.Count; i++)
{
if (Equaltion(line1, lines[i]))
{
lineIdx1 = i;
count1 = lines[i].Count;
}
else if (Equaltion(line2, lines[i]))
{
lineIdx2 = i;
count2 = lines[i].Count;
}
}
/*生成最终L1、L2基线*/
if (count1 < 2 && count2 < 2)
{
/*新三角形的L0基线*/
tri = new Tin_Tri();
tri.L0 = tris[k].L0;
/*新三角形的L1基线*/
if (count1 == 0)
{
line1.Count = 1;
lines.Add(line1);
tri.L1 = lines.Count - 1;
}
else
{
lines[lineIdx1].Count = 2;
tri.L1 = lineIdx1;
}
/*新三角形的L2基线*/
if (count2 == 0)
{
line2.Count = 1;
lines.Add(line2);
tri.L2 = lines.Count - 1;
}
else
{
lines[lineIdx2].Count = 2;
tri.L2 = lineIdx2;
}
tris.Add(tri);
/*第K个三角形的L0基线拓展三角形生成结束*/
}
}
}
/*不论是否生成新三角形,最终都要将 第K个三角形的L0边 使用次数改为2*/
lines[tris[k].L0].Count = 2;
}
/*处理第K个三角形的L1边*/
if (lines[tris[k].L1].Count < 2)
{
/*获取第K个三角形,不是L1边上的端点*/
int test_point = -1;
if (lines[tris[k].L1].P0 == lines[tris[k].L0].P0 ||
lines[tris[k].L1].P1 == lines[tris[k].L0].P0)
{
test_point = lines[tris[k].L0].P1;
}
if (lines[tris[k].L1].P0 == lines[tris[k].L0].P1 ||
lines[tris[k].L1].P1 == lines[tris[k].L0].P1)
{
test_point = lines[tris[k].L0].P0;
}
/*获取该端点在竖直方向到L1线段的距离*/
double sign1 = Sign(lines[tris[k].L1], test_point);
/*清空可用点集*/
usablepoints.Clear();
/*获取可用点集*/
for (int i = 0; i < points.Count; i++)
{
if (lines[tris[k].L0].P0 != i &&
lines[tris[k].L0].P1 != i &&
lines[tris[k].L1].P0 != i &&
lines[tris[k].L1].P1 != i &&
lines[tris[k].L2].P0 != i &&
lines[tris[k].L2].P1 != i)
{
double sign2 = Sign(lines[tris[k].L1], i);
if (sign1 * sign2 < 0)
{
usablepoints.Add(i);
}
}
}
/*生成下一个三角形*/
if (usablepoints.Count > 0)
{
/*根据张角最大原则获取第三点*/
index = Max(lines[tris[k].L1], usablepoints);
/*验证第三点是否满足空外接圆准则*/
if (Circumcircle(lines[tris[k].L1], index))
{
/*生成临时L1、L2基线*/
line1 = new Tin_Line();
line1.P0 = lines[tris[k].L1].P0;
line1.P1 = index;
line2 = new Tin_Line();
line2.P0 = lines[tris[k].L1].P1;
line2.P1 = index;
/*标记临时基线是否存在的标志*/
int lineIdx1 = -1;
int lineIdx2 = -1;
int count1 = 0;
int count2 = 0;
/*判断基线数据集中是否存在这两条边*/
for (int i = 0; i < lines.Count; i++)
{
if (Equaltion(line1, lines[i]))
{
lineIdx1 = i;
count1 = lines[i].Count;
}
else if (Equaltion(line2, lines[i]))
{
lineIdx2 = i;
count2 = lines[i].Count;
}
}
/*生成最终L1、L2基线*/
if (count1 < 2 && count2 < 2)
{
/*新三角形的L0基线*/
tri = new Tin_Tri();
tri.L0 = tris[k].L1;
/*新三角形的L1基线*/
if (count1 == 0)
{
line1.Count = 1;
lines.Add(line1);
tri.L1 = lines.Count - 1;
}
else
{
lines[lineIdx1].Count = 2;
tri.L1 = lineIdx1;
}
/*新三角形的L2基线*/
if (count2 == 0)
{
line2.Count = 1;
lines.Add(line2);
tri.L2 = lines.Count - 1;
}
else
{
lines[lineIdx2].Count = 2;
tri.L2 = lineIdx2;
}
tris.Add(tri);
/*第K个三角形的L1基线拓展三角形生成结束*/
}
}
}
/*不论是否生成新三角形,最终都要将 第K个三角形的L1边 使用次数改为2*/
lines[tris[k].L1].Count = 2;
}
/*处理第K个三角形的L2边*/
if (lines[tris[k].L2].Count < 2)
{
/*获取第K个三角形,不是L2边上的端点*/
int test_point = -1;
if (lines[tris[k].L2].P0 == lines[tris[k].L0].P0 ||
lines[tris[k].L2].P1 == lines[tris[k].L0].P0)
{
test_point = lines[tris[k].L0].P1;
}
if (lines[tris[k].L2].P0 == lines[tris[k].L0].P1 ||
lines[tris[k].L2].P1 == lines[tris[k].L0].P1)
{
test_point = lines[tris[k].L0].P0;
}
/*获取该端点在竖直方向到L2线段的距离*/
double sign1 = Sign(lines[tris[k].L2], test_point);
/*清空可用点集*/
usablepoints.Clear();
/*获取可用点集*/
for (int i = 0; i < points.Count; i++)
{
if (lines[tris[k].L0].P0 != i &&
lines[tris[k].L0].P1 != i &&
lines[tris[k].L1].P0 != i &&
lines[tris[k].L1].P1 != i &&
lines[tris[k].L2].P0 != i &&
lines[tris[k].L2].P1 != i)
{
double sign2 = Sign(lines[tris[k].L2], i);
if (sign1 * sign2 < 0)
{
usablepoints.Add(i);
}
}
}
/*生成下一个三角形*/
if (usablepoints.Count > 0)
{
/*根据张角最大原则获取第三点*/
index = Max(lines[tris[k].L2], usablepoints);
/*验证第三点是否满足空外接圆准则*/
if (Circumcircle(lines[tris[k].L2], index))
{
/*生成临时L1、L2基线*/
line1 = new Tin_Line();
line1.P0 = lines[tris[k].L2].P0;
line1.P1 = index;
line2 = new Tin_Line();
line2.P0 = lines[tris[k].L2].P1;
line2.P1 = index;
/*标记临时基线是否存在的标志*/
int lineIdx1 = -1;
int lineIdx2 = -1;
int count1 = 0;
int count2 = 0;
/*判断基线数据集中是否存在这两条边*/
for (int i = 0; i < lines.Count; i++)
{
if (Equaltion(line1, lines[i]))
{
lineIdx1 = i;
count1 = lines[i].Count;
}
else if (Equaltion(line2, lines[i]))
{
lineIdx2 = i;
count2 = lines[i].Count;
}
}
/*生成最终L1、L2基线*/
if (count1 < 2 && count2 < 2)
{
/*新三角形的L0基线*/
tri = new Tin_Tri();
tri.L0 = tris[k].L2;
/*新三角形的L1基线*/
if (count1 == 0)
{
line1.Count = 1;
lines.Add(line1);
tri.L1 = lines.Count - 1;
}
else
{
lines[lineIdx1].Count = 2;
tri.L1 = lineIdx1;
}
/*新三角形的L2基线*/
if (count2 == 0)
{
line2.Count = 1;
lines.Add(line2);
tri.L2 = lines.Count - 1;
}
else
{
lines[lineIdx2].Count = 2;
tri.L2 = lineIdx2;
}
tris.Add(tri);
/*第K个三角形的L2基线拓展三角形生成结束*/
}
}
}
/*不论是否生成新三角形,最终都要将 第K个三角形的L1边 使用次数改为2*/
lines[tris[k].L2].Count = 2;
}
k++;
if (k == tris.Count)
{
break;
}
}
}
}
}
|