三角网生长算法构建TIN
程序员文章站
2022-05-20 19:45:45
...
三角网生长算法构建TIN
本文内容实现均用C#
1、构建三角形和点存储结构
public string Filename;
public System.Drawing.Bitmap curBitmap;
public Form1()
{
InitializeComponent();
}
int[] t1 = new int[1000];
int[] t2 = new int[1000];
int[] t3 = new int[1000];
//创建高程点的结构，存储高程点的名称,X、Y坐标,高程H值
public struct Point
{
public int Number;
public string Name; //存储点的名称
public double x; //存储点的X坐标
public double y; //存储点的Y坐标
public double h; //存储点的高程值H
}
Point[] pt = new Point[1000]; //定义初始的点数组大小为1000
int Lines; //记录文件的行数，即点的个数
double xmax, xmin, ymax, ymin; //记录所有点中的x，y坐标最大最小值
int K;
2、打开高程点数据文件
//打开高程点数据文件
private void 打开OToolStripMenuItem_Click(object sender, EventArgs e)
{
OpenFileDialog filename = new OpenFileDialog();
filename.Filter = "All files(*.*)|*.*|txt files(*.txt)|*.txt|dat files(*.dat)|*.dat";
filename.FilterIndex = 2;
//filename.RestoreDirectory = true;
if (filename.ShowDialog() == DialogResult.OK)
{
Filename = filename.FileName.ToString();
string[] lines = File.ReadAllLines(Filename);
Lines = lines.Length;
for (int i = 1; i <= Lines; i++)
{
string[] sArray = lines[i - 1].Split(','); //按","将每一行分割成四个字符串
pt[i].Number = i;
pt[i].Name = sArray[0];
pt[i].x = Convert.ToDouble(sArray[1]);
pt[i].y = Convert.ToDouble(sArray[2]);
pt[i].h = Convert.ToDouble(sArray[3]);
}
}
}
3、确定所有点的范围
//确定所有点的范围
private void Area()
{
xmax = xmin = pt[1].x;
ymax = ymin = pt[1].y;
for (int i = 2; i <= Lines; i++)
{
if (xmax < pt[i].x) xmax = pt[i].x;
if (xmin > pt[i].x) xmin = pt[i].x;
if (ymax < pt[i].y) ymax = pt[i].y;
if (ymin > pt[i].y) ymin = pt[i].y;
}
}
4、计算坐标转换比例因子
//计算坐标转换比例因子
public double CalcScale()
{
Area();
Rectangle m_rect = pictureBox1.ClientRectangle;
double ds = 1.0;
double dsx, dsy;
if ((xmax - xmin != 0) && (ymax - ymin != 0))
{
dsx = Math.Abs((xmax - xmin) / m_rect.Height);
dsy = Math.Abs((ymax - ymin) / m_rect.Width);
ds = Math.Max(dsx, dsy);
}
else
{
if (xmax - xmin != 0)
{
ds = Math.Abs((xmax - xmin) / m_rect.Height);
}
else
{
if (ymax - ymin != 0)
{
ds = Math.Abs((ymax - ymin) / m_rect.Width);
}
else { ds = 1; }
}
}
return ds;
}
5、找到最近的两个高程点
//找到两个最近的高程点
public void MinDistance(Point[] pt, out int pt1, out int pt2)
{
int i, j;
double[,] Distance = new double[Lines, Lines];
//将任意两点间的距离存储到矩阵Distance中
for (i = 1; i <= Lines; i++)
for (j = i + 1; j < Lines; j++)
if (i != j)
Distance[i, j] = Math.Sqrt(Math.Pow(pt[i].x - pt[j].x, 2) + Math.Pow(pt[i].y - pt[j].y, 2));
double[] Mindistance = { 10000, 0, 0 };
//找到矩阵Distance中的最小值，并记录行列号
for (i = 1; i <= Lines; i++)
for (j = i + 1; j < Lines; j++)
if (Mindistance[0] > Distance[i, j])
{
Mindistance[0] = Distance[i, j];
Mindistance[1] = i;
Mindistance[2] = j;
}
pt1 = (int)Mindistance[1];
pt2 = (int)Mindistance[2];
}
6、找到离中点最近的点
//找到离中点最近的点
public void Find(int pt1, int pt2, out int pt3)
{
int i;
double meanx = (pt[pt1].x + pt[pt2].x) / 2;
double meany = (pt[pt1].y + pt[pt2].y) / 2;
double Min = 10000000000;
pt3 = 0;
for (i = 1; i <= Lines; i++)
{
if (i != pt1 && i != pt2)
{
double temp = Math.Sqrt(Math.Pow(pt[i].x - meanx, 2) + Math.Pow(pt[i].y - meany, 2));
if (Min > temp)
{
Min = temp;
pt3 = i;
}
}
}
}
7、判断三角形扩展是否在同一侧
//判断三角形扩展点是否在同一侧
public bool Direction(int point1, int point2, int point3, int point4)
{
//计算直线方程的系数a，b
double a = (pt[point2].y - pt[point1].y) / (pt[point2].x - pt[point1].x);
double b = (pt[point1].x * pt[point2].y - pt[point2].x * pt[point1].y) / (pt[point2].x - pt[point1].x);
double fxy1 = pt[point3].y - (a * pt[point3].x - b);
double fxy2 = pt[point4].y - (a * pt[point4].x - b);
//当位于非同一侧时
if (fxy1 < 0 && fxy2 > 0 || fxy1 > 0 && fxy2 < 0)
return true;
//当位于同一侧时
else return false;
}
8、计算扩展边的角度余弦值
//计算扩展边的角度余弦值
public double Angle(int pt1, int pt2, int pt3)
{
double angle;
double L1 = Math.Sqrt((pt[pt2].x - pt[pt3].x) * (pt[pt2].x - pt[pt3].x) + (pt[pt2].y - pt[pt3].y) * (pt[pt2].y - pt[pt3].y));
double L2 = Math.Sqrt((pt[pt1].x - pt[pt3].x) * (pt[pt1].x - pt[pt3].x) + (pt[pt1].y - pt[pt3].y) * (pt[pt1].y - pt[pt3].y));
double L3 = Math.Sqrt((pt[pt2].x - pt[pt1].x) * (pt[pt2].x - pt[pt1].x) + (pt[pt2].y - pt[pt1].y) * (pt[pt2].y - pt[pt1].y));
angle = (L1 * L1 + L2 * L2 - L3 * L3) / (2 * L1 * L2);
return angle;
}
9、找到扩展边形成张角最大的点
//找到扩展边形成张角最大的点
private int MaxAngle(int[] x, int A, int B, int n)
{
double C = 0, temp, s = 0;
int max = 0;
for (int i = 1; i <= n; i++)
{
if (x[i] != A && x[i] != B)
{
s = Angle(A, B, x[i]);
if (s < 1)
C = Math.Acos(s);
else C = 0;
max = x[i];
break;
}
}
for (int i = 1; i <= n; i++)
{
if (i != A && i != B)
{
s = Angle(A, B, x[i]);
if (s < 1)
temp = Math.Acos(s);
else temp = 0;
if (temp > C)
{
C = temp;
max = x[i];
}
}
}
return max;
}
10、判断三角形的一条边是否已经出现过两次
//判断三角形的一条边是否已经出现过两次
public bool Repeat(int point1, int point2, int L)
{
int sum = 0;
for (int i = 1; i <= L; i++)
{
if (point1 == t1[i] && point2 == t2[i] || point1 == t2[i] && point2 == t1[i] ||
point1 == t2[i] && point2 == t3[i] || point1 == t3[i] && point2 == t2[i] ||
point1 == t3[i] && point2 == t1[i] || point1 == t1[i] && point2 == t3[i])
{
sum++;
if (sum == 2)
return false;
}
}
return true;
}
12、构建并绘制TIN
private void 构建TINCToolStripMenuItem_Click(object sender, EventArgs e)
{
//找到所有点中距离最小的两个点，作为第一个三角形的第一个点和第二个点
MinDistance(pt, out int point1, out int point2);
t1[1] = point1;
t2[1] = point2;
//寻找第一个三角形的第三个点：离第一条边距离最短的点
Find(point1, point2, out int point3);
t3[1] = point3;
//设置计数变量K记录扩展的三角形数
K = 0;
//设置计数变量L记录已经形成的三角形数
int L = 1;
//设置数组存储可能的扩展点
int[] x = new int[Lines];
//扩展三角形
while (K != L)
{
K++;
point1 = t1[K];
point2 = t2[K];
point3 = t3[K];
//第一条扩展边不重复，没有被两个三角形共用
if (Repeat(point1, point2, L))
{
//判断新扩展的边
int t = 0;
x[t++] = 0;
//寻找可能的扩展点
for (int i = 1; i <= Lines; i++)
if (i != point1 && i != point2 && i != point3 && Direction(point1, point2,point3, i))
{
x[t++] = i;
}
//存在扩展点
if (t > 1)
{
int max = MaxAngle(x, point1, point2, t - 1);
L = L + 1;
t1[L] = point1;
t2[L] = point2;
t3[L] = max;
}
}
//第二条扩展边不重复，没有被两个三角形共用
if (Repeat(point1, point3, L))
{
int t = 0;
x[t++] = 0;
for (int i = 1; i <= Lines; i++)
if (i != point1 && i != point3 && i != point2 && Direction(point1, point3, point2, i))
{
x[t++] = i;
}
if (t > 1)
{
int max = MaxAngle(x, point1, point3, t - 1);
L = L + 1;
t1[L] = point1;
t2[L] = point3;
t3[L] = max;
}
}
//第三条扩展边不重复，没有被两个三角形共用
if (Repeat(point2, point3, L))
{
int t = 0;
x[t++] = 0;
for (int i = 1; i <= Lines; i++)
if (i != point2 && i != point3 && i != point1 && Direction(point2, point3, point1, i))
{
x[t++] = i;
}
if (t > 1)
{
int max = MaxAngle(x, point2, point3, t - 1);
L = L + 1;
t1[L] = point2;
t2[L] = point3;
t3[L] = max;
}
}
}
//绘制TIN
Graphics g = pictureBox1.CreateGraphics();
double m_scale = CalcScale();
Pen mypen = new Pen(Color.Red, 1); //创建画笔
Rectangle m_rect = pictureBox1.ClientRectangle; //获得画布大小
g.SmoothingMode = System.Drawing.Drawing2D.SmoothingMode.HighQuality; //消除锯齿
for (int i = 1; i <= L; i++)
{
//由测量坐标计算屏幕坐标
double ix1 = (pt[t1[i]].y - ymin) / m_scale;
double iy1 = m_rect.Height - (pt[t1[i]].x - xmin) / m_scale - 20;
double ix2 = (pt[t2[i]].y - ymin ) / m_scale;
double iy2 = m_rect.Height - (pt[t2[i]].x - xmin) / m_scale - 20;
double ix3 = (pt[t3[i]].y - ymin) / m_scale;
double iy3 = m_rect.Height - (pt[t3[i]].x - xmin) / m_scale - 20;
g.DrawLine(mypen, (float)ix1, (float)iy1, (float)ix2, (float)iy2);
g.DrawLine(mypen, (float)ix1, (float)iy1, (float)ix3, (float)iy3);
g.DrawLine(mypen, (float)ix3, (float)iy3, (float)ix2, (float)iy2);
}
//g.Dispose();
}
展点结果
设计界面
功能菜单设置-1
绘制TIN结果
高程点文件内容:()
1,512851.699,1121834.867,210.310
2,512850.236,1121835.836,210.430
3,512848.802,1121836.824,210.310
4,512847.151,1121814.595,210.290
5,512849.057,1121831.194,210.440
6,512846.319,1121834.680,210.360
7,512847.872,1121828.652,210.430
8,512844.406,1121835.692,210.410
9,512849.974,1121822.530,210.350
10,512847.420,1121826.715,210.330
11,512841.893,1121816.777,210.380
12,512842.204,1121833.716,210.390
13,512845.495,1121825.847,210.410
14,512847.488,1121821.435,210.470
15,512839.877,1121837.289,210.520
16,512842.601,1121829.914,210.470
17,512847.693,1121830.249,210.360
18,512843.665,1121825.828,210.440
19,512841.285,1121830.542,210.350
20,512850.042,1121810.872,210.320
21,512850.820,1121807.228,210.350
22,512839.262,1121831.614,210.430
23,512840.900,1121826.505,210.360
24,512847.948,1121806.836,210.330
25,512839.078,1121822.817,210.360
26,512851.561,1121787.907,210.460
27,512836.039,1121822.946,210.330
28,512843.538,1121805.103,210.410
29,512845.537,1121792.055,211.540
30,512848.374,1121791.586,211.930
31,512849.211,1121789.581,211.860
32,512831.485,1121831.862,210.350
33,512841.485,1121810.862,210.350
34,512835.537,1121810.055,211.540
35,512838.540,1121805.496,210.528
36,512830.158,1121800.693,211.359
37,512840.088,1121790.363,210.358
推荐阅读