三角网内插高程点
Imports Autodesk.AutoCAD.ApplicationServicesImports Autodesk.AutoCAD.EditorInput
Imports Autodesk.AutoCAD.DatabaseServices
Imports Autodesk.AutoCAD.Geometry
Module TriangleElevation
'三角形内插入高程点
Sub Main()
Dim docLock As DocumentLock = Core.Application.DocumentManager.MdiActiveDocument.LockDocument() '“非模态窗口，要锁定文档”
NativeMethods.SetFocus(Core.Application.DocumentManager.MdiActiveDocument.Window.Handle) 'CAD获得焦点
' 获取当前文档和数据库
Dim doc As Document = Core.Application.DocumentManager.MdiActiveDocument
Dim db As Database = doc.Database
Dim ed As Editor = doc.Editor
' 提示用户选择三角形（假设三角形是由 3D 面或多段线表示）
Dim peo As New PromptSelectionOptions()
peo.MessageForAdding = "请选择三角形（由 3D 面或多段线表示）: "
Dim psr As PromptSelectionResult = ed.GetSelection(peo)
' 检查选择结果是否有效
If psr.Status <> PromptStatus.OK Then
ed.WriteMessage("未选择任何三角形对象。")
Return
End If
Do
' 遍历选择集中的每个三角形对象
Using tr As Transaction = db.TransactionManager.StartTransaction()
' 提示用户在三角形内拾取一点
Dim ppr As PromptPointResult = ed.GetPoint("请在三角形内拾取一点: ")
' 检查拾取点结果是否有效
If ppr.Status <> PromptStatus.OK Then
ed.WriteMessage("未拾取任何点。")
Exit Do
Else
For Each id As ObjectId In psr.Value.GetObjectIds()
Dim ent As Entity = tr.GetObject(id, OpenMode.ForRead)
' 检查对象是否为 Polyline 或 Polyline3d
If TypeOf ent Is Polyline OrElse TypeOf ent Is Polyline3d Then
' 获取三角形顶点
Dim vertices As Point3dCollection = GetTriangleVertices(ent, tr)
' 检查拾取的点是否在三角形内
Dim pickedPoint As Point3d = ppr.Value
If IsPointInTriangle(vertices, pickedPoint) Then
' 计算拾取点的高程值
Dim elevation As Double = CalculateElevation(vertices, pickedPoint)
' 在图上生成点和文本
CreatePointAndText(db, tr, pickedPoint, elevation)
' 刷新图形视图以显示新创建的点和文本
ed.UpdateScreen()
Exit For
Else
ed.WriteMessage("拾取的点不在三角形内。")
End If
End If
Next
End If
tr.Commit()
End Using
Loop
docLock.Dispose()'解锁文档
End Sub
' 获取三角形的顶点
Function GetTriangleVertices(ent As Entity, tr As Transaction) As Point3dCollection
Dim vertices As New Point3dCollection()
If TypeOf ent Is Polyline Then
Dim pline As Polyline = CType(ent, Polyline)
For i As Integer = 0 To pline.NumberOfVertices - 1
vertices.Add(pline.GetPoint3dAt(i))
Next
ElseIf TypeOf ent Is Polyline3d Then
Dim pline3d As Polyline3d = CType(ent, Polyline3d)
For Each vId As ObjectId In pline3d
Dim v As PolylineVertex3d = CType(tr.GetObject(vId, OpenMode.ForRead), PolylineVertex3d)
vertices.Add(v.Position)
Next
End If
Return vertices
End Function
' 计算拾取点的高程值
Function CalculateElevation(vertices As Point3dCollection, pickedPoint As Point3d) As Double
' 确保三角形是平面的
If vertices.Count <> 3 Then Return 0.0
' 使用重心坐标法计算高程值
Dim x1 As Double = vertices(0).X
Dim y1 As Double = vertices(0).Y
Dim z1 As Double = vertices(0).Z
Dim x2 As Double = vertices(1).X
Dim y2 As Double = vertices(1).Y
Dim z2 As Double = vertices(1).Z
Dim x3 As Double = vertices(2).X
Dim y3 As Double = vertices(2).Y
Dim z3 As Double = vertices(2).Z
Dim x As Double = pickedPoint.X
Dim y As Double = pickedPoint.Y
' 计算重心坐标
Dim detT As Double = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3)
Dim lambda1 As Double = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3)) / detT
Dim lambda2 As Double = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3)) / detT
Dim lambda3 As Double = 1.0 - lambda1 - lambda2
' 计算插值高程
Dim elevation As Double = lambda1 * z1 + lambda2 * z2 + lambda3 * z3
Return elevation
End Function
' 判断点是否在三角形内
Function IsPointInTriangle(vertices As Point3dCollection, point As Point3d) As Boolean
Dim x1 As Double = vertices(0).X
Dim y1 As Double = vertices(0).Y
Dim x2 As Double = vertices(1).X
Dim y2 As Double = vertices(1).Y
Dim x3 As Double = vertices(2).X
Dim y3 As Double = vertices(2).Y
Dim x As Double = point.X
Dim y As Double = point.Y
Dim denominator As Double = ((y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3))
Dim a As Double = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3)) / denominator
Dim b As Double = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3)) / denominator
Dim c As Double = 1 - a - b
' 检查点是否在三角形内
Return 0 <= a AndAlso a <= 1 AndAlso 0 <= b AndAlso b <= 1 AndAlso 0 <= c AndAlso c <= 1
End Function
' 在图上生成点和文本
Sub CreatePointAndText(db As Database, tr As Transaction, point As Point3d, elevation As Double)
Dim bt As BlockTable = CType(tr.GetObject(db.BlockTableId, OpenMode.ForRead), BlockTable)
Dim btr As BlockTableRecord = CType(tr.GetObject(bt(BlockTableRecord.ModelSpace), OpenMode.ForWrite), BlockTableRecord)
' 创建点对象
Dim dbPoint As New DBPoint(point) With {
.ColorIndex = 1 ' 红色
}
AppendEntity(dbPoint)
' 创建文本对象
Dim text As New DBText With {
.Position = New Point3d(point.X, point.Y, point.Z + 0.5), ' 将文本位置稍微偏高
.TextString = elevation.ToString("F2"),
.Height = 1.0,
.ColorIndex = 2 ' 黄色
}
AppendEntity(text)
End Sub
End Module 应该加上四叉树，不然会很影响速度，就没啥实际意义了 判断点是否在三角形内的函数应加上容差判断，当点落于三角形边上时，a、b、c三个变量中会有1个或2个非常接近0，此时该函数有时会判断在外，有时又判断在内，故调整为如下，其中的容差xEps取值为10的负10次方。 ' 判断点在三角形内
Public Function IsPointInTriangle(ByVal vertices As Point3dCollection, ByVal point As Point3d, ByVal xEps As Double) As Boolean
Dim x1 As Double = vertices(0).X
Dim y1 As Double = vertices(0).Y
Dim x2 As Double = vertices(1).X
Dim y2 As Double = vertices(1).Y
Dim x3 As Double = vertices(2).X
Dim y3 As Double = vertices(2).Y
Dim x As Double = point.X
Dim y As Double = point.Y
Dim denominator As Double = ((y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3))
Dim a As Double = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3)) / denominator
Dim b As Double = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3)) / denominator
Dim c As Double = 1 - a - b
If Math.Abs(a) <= xEps Then a = 0.0
If Math.Abs(b) <= xEps Then b = 0.0
If Math.Abs(c) <= xEps Then c = 0.0
' 检查点是否在三角形内
Return 0 <= a AndAlso a <= 1 AndAlso 0 <= b AndAlso b <= 1 AndAlso 0 <= c AndAlso c <= 1
End Function 本帖最后由 czb203 于 2025-2-24 10:50 编辑
using System;
using System.Collections.Generic;
using System.Linq;
public struct Point
{
public double X { get; }
public double Y { get; }
public double Z { get; }
public Point(double x, double y, double z)
{
X = x;
Y = y;
Z = z;
}
}
public class Triangle
{
public Point P1 { get; }
public Point P2 { get; }
public Point P3 { get; }
public double MinX { get; }
public double MaxX { get; }
public double MinY { get; }
public double MaxY { get; }
public Triangle(Point p1, Point p2, Point p3)
{
P1 = p1;
P2 = p2;
P3 = p3;
MinX = Math.Min(p1.X, Math.Min(p2.X, p3.X));
MaxX = Math.Max(p1.X, Math.Max(p2.X, p3.X));
MinY = Math.Min(p1.Y, Math.Min(p2.Y, p3.Y));
MaxY = Math.Max(p1.Y, Math.Max(p2.Y, p3.Y));
}
}
public class QuadTree
{
private const int MaxCapacity = 4;
private readonly List<Triangle> triangles = new List<Triangle>();
private QuadTree[] children;
private readonly double xMin, yMin, xMax, yMax;
public QuadTree(double xMin, double yMin, double xMax, double yMax)
{
this.xMin = xMin;
this.yMin = yMin;
this.xMax = xMax;
this.yMax = yMax;
}
public void Insert(Triangle triangle)
{
if (children != null)
{
foreach (var child in children)
{
if (ChildIntersectsTriangle(child, triangle))
{
child.Insert(triangle);
}
}
return;
}
triangles.Add(triangle);
if (triangles.Count > MaxCapacity)
{
Split();
foreach (var t in triangles)
{
foreach (var child in children)
{
if (ChildIntersectsTriangle(child, t))
{
child.Insert(t);
}
}
}
triangles.Clear();
}
}
private bool ChildIntersectsTriangle(QuadTree child, Triangle triangle)
{
return !(triangle.MinX > child.xMax ||
triangle.MaxX < child.xMin ||
triangle.MinY > child.yMax ||
triangle.MaxY < child.yMin);
}
private void Split()
{
double xMid = (xMin + xMax) / 2;
double yMid = (yMin + yMax) / 2;
children = new QuadTree;
children = new QuadTree(xMin, yMin, xMid, yMid);
children = new QuadTree(xMid, yMin, xMax, yMid);
children = new QuadTree(xMin, yMid, xMid, yMax);
children = new QuadTree(xMid, yMid, xMax, yMax);
}
public List<Triangle> Query(double x, double y)
{
var result = new List<Triangle>();
if (x < xMin || x > xMax || y < yMin || y > yMax)
return result;
if (children != null)
{
foreach (var child in children)
{
result.AddRange(child.Query(x, y));
}
}
else
{
result.AddRange(triangles);
}
return result;
}
}
public static class TriangleHelper
{
public static bool IsPointInTriangle(Point p, Triangle triangle, double epsilon = 1e-6)
{
double denominator = ((triangle.P2.Y - triangle.P3.Y) * (triangle.P1.X - triangle.P3.X)) +
((triangle.P3.X - triangle.P2.X) * (triangle.P1.Y - triangle.P3.Y));
if (Math.Abs(denominator) < epsilon)
return false;
double a = ((triangle.P2.Y - triangle.P3.Y) * (p.X - triangle.P3.X) +
(triangle.P3.X - triangle.P2.X) * (p.Y - triangle.P3.Y)) / denominator;
double b = ((triangle.P3.Y - triangle.P1.Y) * (p.X - triangle.P3.X) +
(triangle.P1.X - triangle.P3.X) * (p.Y - triangle.P3.Y)) / denominator;
double c = 1 - a - b;
return a >= -epsilon && a <= 1 + epsilon &&
b >= -epsilon && b <= 1 + epsilon &&
c >= -epsilon && c <= 1 + epsilon &&
Math.Abs(a + b + c - 1) < 3 * epsilon;
}
public static double InterpolateZ(Point p, Triangle triangle)
{
double denominator = ((triangle.P2.Y - triangle.P3.Y) * (triangle.P1.X - triangle.P3.X)) +
((triangle.P3.X - triangle.P2.X) * (triangle.P1.Y - triangle.P3.Y));
double a = ((triangle.P2.Y - triangle.P3.Y) * (p.X - triangle.P3.X) +
(triangle.P3.X - triangle.P2.X) * (p.Y - triangle.P3.Y)) / denominator;
double b = ((triangle.P3.Y - triangle.P1.Y) * (p.X - triangle.P3.X) +
(triangle.P1.X - triangle.P3.X) * (p.Y - triangle.P3.Y)) / denominator;
double c = 1 - a - b;
return a * triangle.P1.Z + b * triangle.P2.Z + c * triangle.P3.Z;
}
}
public class TerrainInterpolator
{
private readonly QuadTree quadTree;
public TerrainInterpolator(IEnumerable<Triangle> triangles)
{
var triangleList = triangles.ToList();
double xMin = triangleList.Min(t => t.MinX);
double xMax = triangleList.Max(t => t.MaxX);
double yMin = triangleList.Min(t => t.MinY);
double yMax = triangleList.Max(t => t.MaxY);
quadTree = new QuadTree(xMin, yMin, xMax, yMax);
foreach (var triangle in triangleList)
{
quadTree.Insert(triangle);
}
}
public double? Interpolate(Point p)
{
var candidates = quadTree.Query(p.X, p.Y);
foreach (var triangle in candidates)
{
if (TriangleHelper.IsPointInTriangle(p, triangle))
{
return TriangleHelper.InterpolateZ(p, triangle);
}
}
return null;
}
}
// 使用示例
class Program
{
static void Main()
{
// 创建测试三角形
var triangles = new List<Triangle>
{
new Triangle(
new Point(0, 0, 100),
new Point(10, 0, 200),
new Point(5, 10, 300)
)
// 添加更多三角形...
};
var interpolator = new TerrainInterpolator(triangles);
Point testPoint = new Point(5, 5, 0);
double? elevation = interpolator.Interpolate(testPoint);
if (elevation.HasValue)
{
Console.WriteLine($"插值高程值: {elevation.Value}");
}
else
{
Console.WriteLine("点不在任何三角形内");
}
}
}
页:
[1]