You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
748 lines
38 KiB
748 lines
38 KiB
using System; |
|
using System.Collections.Generic; |
|
using System.Runtime.InteropServices; |
|
using ESRI.ArcGIS.Geometry; |
|
using GeoAPI.Operation.Buffer; |
|
using KGIS.Framework.AE; |
|
using NetTopologySuite.IO; |
|
using NetTopologySuite.Operation.Buffer; |
|
|
|
namespace Kingo.PluginServiceInterface.Helper |
|
{ |
|
/// <summary> |
|
/// 狭长图形检测工具类（静态类） |
|
/// 功能：检测多边形中是否存在狭长区域（如长条状、细带状图斑） |
|
/// 核心算法思路：通过分析多边形环的顶点分布，识别局部密集点对，构建子多边形并通过紧凑度指标判断是否为狭长区域 |
|
/// </summary> |
|
public static class PolygonNarrowAreaChecker |
|
{ |
|
/// <summary> |
|
/// 检测多边形中的局部狭长区域 |
|
/// </summary> |
|
/// <param name="polygon">待检测的输入多边形（支持带内环的复杂多边形）</param> |
|
/// <param name="distanceThreshold">距离阈值（单位：与坐标系一致，默认1）：用于判断两个顶点是否"过近"，作为狭长区域的初步筛选条件</param> |
|
/// <returns>检测到的狭长区域多边形列表（每个元素为一个狭长子多边形）</returns> |
|
public static List<IPolygon> CheckLocalNarrowAreas(IPolygon polygon, double distanceThreshold = 1) |
|
{ |
|
List<IPolygon> narrowAreas = new List<IPolygon>(); // 存储检测到的狭长区域 |
|
//var xxx = GeometryConvertHelper.ConvertIGeoemtryToWKT(polygon); |
|
string sNarrow = ISNarrow(polygon); |
|
if (!string.IsNullOrEmpty(sNarrow)) |
|
{ |
|
#region 删除图形中的重复点（首位点除外） |
|
IGeometryCollection geomColl = polygon as IGeometryCollection; |
|
// 遍历所有子环（外环和内环） |
|
for (int i = 0; i < geomColl.GeometryCount; i++) |
|
{ |
|
IGeometry subGeom = geomColl.get_Geometry(i); |
|
IRing ring = subGeom as IRing; // 环是闭合的线（IRing 继承自 IPolyline） |
|
if (ring == null) continue; |
|
// 处理当前环的重复点 |
|
RemoveDuplicatePointsInRing(ring); |
|
} |
|
#endregion |
|
IGeometryCollection geometryCollection = polygon as IGeometryCollection;// 将多边形转换为几何集合（用于遍历内部环） |
|
if (geometryCollection == null) return narrowAreas; // 非几何集合类型直接返回空列表 |
|
for (int ringIndex = 0; ringIndex < geometryCollection.GeometryCount; ringIndex++)// 遍历多边形的所有环（外环+内环） |
|
{ |
|
IRing ring = geometryCollection.get_Geometry(ringIndex) as IRing; // 获取当前环 |
|
if (ring == null) continue; // 非环类型跳过 |
|
IPointCollection points = ring as IPointCollection;// 将环转换为点集合（用于获取顶点坐标） |
|
int pointCount = points.PointCount; // 环的总顶点数 |
|
bool isClosed = ring.IsClosed; // 判断环是否闭合（多边形环通常需闭合） |
|
int effectivePointCount = isClosed ? pointCount - 1 : pointCount; // 有效顶点数（闭合环需排除最后一个重复点） |
|
// 【步骤1】预处理环的顶点坐标：转换为结构体数组（提升计算性能）和IPoint对象数组（用于几何构建） |
|
var (ringPoints, ringPointsIPoint) = PreprocessRingPoints(points, effectivePointCount); if (ringPoints == null) continue; // 预处理失败时跳过当前环 |
|
HashSet<(double X, double Y)> usedPoints = new HashSet<(double X, double Y)>(); |
|
// 预计算所有线段的AABB（索引k对应线段k→k+1） |
|
var segmentAABBs = new (double MinX, double MaxX, double MinY, double MaxY)[effectivePointCount - 1]; |
|
for (int idx = 0; idx < effectivePointCount - 1; idx++) |
|
{ |
|
IPoint pStart = ringPointsIPoint[idx]; |
|
IPoint pEnd = ringPointsIPoint[idx + 1]; |
|
segmentAABBs[idx] = ( |
|
MinX: Math.Min(pStart.X, pEnd.X), |
|
MaxX: Math.Max(pStart.X, pEnd.X), |
|
MinY: Math.Min(pStart.Y, pEnd.Y), |
|
MaxY: Math.Max(pStart.Y, pEnd.Y) |
|
); |
|
} |
|
// 遍历所有可能的线段对（需优化遍历范围，避免O(n²)复杂度） |
|
for (int k = 0; k < effectivePointCount - 1; k++) |
|
{ |
|
IPoint p1 = ringPointsIPoint[k]; |
|
IPoint p2 = ringPointsIPoint[k + 1]; |
|
// 检查p1/p2是否已被使用（通过坐标判断） |
|
if (usedPoints.Contains((p1.X, p1.Y)) || usedPoints.Contains((p2.X, p2.Y))) continue; |
|
// 获取当前线段k的AABB |
|
var aabbK = segmentAABBs[k]; |
|
for (int m = k + 2; m < effectivePointCount - 1; m++) // 跳过相邻线段 |
|
{ |
|
IPoint p3 = ringPointsIPoint[m]; |
|
IPoint p4 = ringPointsIPoint[m + 1]; |
|
// 检查p3/p4是否已被使用 |
|
if (usedPoints.Contains((p3.X, p3.Y)) || usedPoints.Contains((p4.X, p4.Y))) continue; |
|
// 获取对比线段m的AABB |
|
var aabbM = segmentAABBs[m]; |
|
// AABB预筛选：判断两线段的包围盒是否可能相交 |
|
// 条件：x轴投影重叠且y轴投影重叠 |
|
double expandDistance = 1;//相当于外扩 |
|
bool xOverlap = (aabbK.MaxX + expandDistance) >= (aabbM.MinX - expandDistance) && |
|
(aabbK.MinX - expandDistance) <= (aabbM.MaxX + expandDistance); |
|
bool yOverlap = (aabbK.MaxY + expandDistance) >= (aabbM.MinY - expandDistance) && |
|
(aabbK.MinY - expandDistance) <= (aabbM.MaxY + expandDistance); |
|
if (!xOverlap || !yOverlap) continue; // AABB不重叠，直接跳过后续复杂计算 |
|
// 判断线段是否共点（如p1=p3, p2=p4等） |
|
if (ArePointsEqual(p1, p3) || ArePointsEqual(p1, p4) || ArePointsEqual(p2, p3) || ArePointsEqual(p2, p4)) continue; |
|
// 计算线段(p1,p2)和(p3,p4)的距离（仅保留AABB重叠的情况） |
|
Vector2D vecAB = CalculateVector(p1, p2); |
|
Vector2D vecCD = CalculateVector(p3, p4); |
|
bool isSameDirection = AreVectorsSameDirection(vecAB, vecCD, 0.01); |
|
if (isSameDirection) continue; |
|
var vangle = CalculateAngle(vecAB, vecCD); |
|
if (vangle < 25 || vangle > 155) continue; |
|
double distance = CalculateSegmentDistance(p1, p2, p3, p4); |
|
// var xxx1 = $@" |
|
//{GeometryConvertHelper.ConvertIGeoemtryToWKT(p1)} |
|
//{GeometryConvertHelper.ConvertIGeoemtryToWKT(p2)} |
|
//{GeometryConvertHelper.ConvertIGeoemtryToWKT(p3)} |
|
//{GeometryConvertHelper.ConvertIGeoemtryToWKT(p4)}"; |
|
if (distance < 0.1) |
|
{ |
|
IPoint iIntersection = GetExtendedLineIntersection(p1, p2, p3, p4); //两条线段的交点 |
|
if (iIntersection != null) |
|
{ |
|
double p1distance = PointToSegmentDistance(p1, p3, p4); |
|
double p2distance = PointToSegmentDistance(p2, p3, p4); |
|
double p3distance = PointToSegmentDistance(p3, p1, p2); |
|
double p4distance = PointToSegmentDistance(p4, p1, p2); |
|
int pindex_k = p1distance > p2distance ? k + 1 : k; |
|
int pindex_m = p3distance > p4distance ? m + 1 : m; |
|
double iIntersection12 = PointToSegmentDistance(iIntersection, p1, p2); |
|
double iIntersection34 = PointToSegmentDistance(iIntersection, p3, p4); |
|
IPoint point = null; |
|
if (iIntersection34 < iIntersection12) |
|
{ |
|
point = p1distance > p2distance ? p2 : p1; |
|
} |
|
else |
|
{ |
|
point = p3distance > p4distance ? p4 : p3; |
|
} |
|
try |
|
{ |
|
IPolyline polyline1 = new PolylineClass(); |
|
polyline1.FromPoint = iIntersection; |
|
polyline1.ToPoint = point; |
|
polyline1.SpatialReference = polygon.SpatialReference; |
|
if (polyline1.Length < 0.1) |
|
{ |
|
polyline1 = getExtendLine(polyline1, 3, 0.01); |
|
} |
|
var CutGeometryresult = FeatureAPI.CutGeometry(polygon, polyline1); |
|
Marshal.ReleaseComObject(polyline1); |
|
if (CutGeometryresult != null && CutGeometryresult.Count == 2) |
|
{ |
|
foreach (var item in CutGeometryresult) |
|
{ |
|
sNarrow = ISNarrow(item); |
|
if (sNarrow == "POLYGON EMPTY") |
|
{ |
|
List<IPoint> pointList = null; |
|
var angle = SharpAngle.GetMinAngle1(item, ref pointList, false); |
|
if (angle < 10) continue; |
|
if (isThreeSides(item)) continue; |
|
narrowAreas.Add(item as IPolygon); |
|
usedPoints.Add((p1.X, p1.Y)); |
|
usedPoints.Add((p2.X, p2.Y)); |
|
usedPoints.Add((p3.X, p3.Y)); |
|
usedPoints.Add((p4.X, p4.Y)); |
|
} |
|
else if (!string.IsNullOrEmpty(sNarrow)) |
|
{ |
|
narrowAreas = CheckLocalNarrowAreas(item as IPolygon); |
|
} |
|
} |
|
} |
|
else if (CutGeometryresult != null && CutGeometryresult.Count == 1) |
|
{ |
|
continue; |
|
//IGeometry geo = CutGeometryresult[0]; |
|
//if (geo != null && !geo.IsEmpty) |
|
//{ |
|
// sNarrow = ISNarrow(geo); |
|
// if (!string.IsNullOrEmpty(sNarrow)) |
|
// { |
|
// List<IPoint> pointList = null; |
|
// var angle = SharpAngle.GetMinAngle1(geo, ref pointList, false); |
|
// if (angle < 10) continue; |
|
// narrowAreas.Add(geo as IPolygon); |
|
// usedPoints.Add((p1.X, p1.Y)); |
|
// usedPoints.Add((p2.X, p2.Y)); |
|
// usedPoints.Add((p3.X, p3.Y)); |
|
// usedPoints.Add((p4.X, p4.Y)); |
|
// continue; |
|
// } |
|
//} |
|
//else continue; |
|
} |
|
} |
|
catch (Exception ex) |
|
{ |
|
continue; |
|
} |
|
} |
|
} |
|
} |
|
} |
|
} |
|
return narrowAreas; |
|
} |
|
return new List<IPolygon>(); |
|
} |
|
|
|
private static IPoint GetExtendedLineIntersection(IPoint start1, IPoint end1, IPoint start2, IPoint end2) |
|
{ |
|
// 创建原始线段 |
|
IPolyline polyline1 = new PolylineClass(); |
|
polyline1.FromPoint = start1; |
|
polyline1.ToPoint = end1; |
|
// 延长线段（假设 getExtendLine 可能返回包含中间点的多段线） |
|
IPolyline extendedPolyline1 = getExtendLine(polyline1, 3, 1); |
|
// 提取延长后的起点和终点，重建为仅含两点的线段 |
|
IPolyline simplifiedLine1 = new PolylineClass(); |
|
simplifiedLine1.FromPoint = extendedPolyline1.FromPoint; // 延长后的起点 |
|
simplifiedLine1.ToPoint = extendedPolyline1.ToPoint; // 延长后的终点 |
|
// 对第二条线段执行相同操作 |
|
IPolyline polyline2 = new PolylineClass(); |
|
polyline2.FromPoint = start2; |
|
polyline2.ToPoint = end2; |
|
IPolyline extendedPolyline2 = getExtendLine(polyline2, 3, 1); |
|
IPolyline simplifiedLine2 = new PolylineClass(); |
|
simplifiedLine2.FromPoint = extendedPolyline2.FromPoint; |
|
simplifiedLine2.ToPoint = extendedPolyline2.ToPoint; |
|
// 计算简化后线段的交点 |
|
ITopologicalOperator topoOp = simplifiedLine1 as ITopologicalOperator; |
|
IGeometry intersection = topoOp.Intersect(simplifiedLine2, esriGeometryDimension.esriGeometry0Dimension); |
|
// 处理交点结果（与原逻辑一致） |
|
if (intersection != null && intersection.GeometryType == esriGeometryType.esriGeometryPoint) |
|
{ |
|
return intersection as IPoint; |
|
} |
|
else if (intersection.GeometryType == esriGeometryType.esriGeometryMultipoint) |
|
{ |
|
IPointCollection points = intersection as IPointCollection; |
|
if (points.PointCount > 0) |
|
return points.get_Point(0); |
|
} |
|
return null; |
|
} |
|
/// <summary> |
|
/// 延长线段 |
|
/// </summary> |
|
/// <param name="passLine">传入去的线</param> |
|
/// <param name="mode">模式，1为从FromPoint处延长，2为从ToPint处延长，3为两端延长</param> |
|
/// <param name="dis">延长的距离</param> |
|
/// <returns></returns> |
|
/// 创建人：懒羊羊 |
|
public static IPolyline getExtendLine(IPolyline passLine, int mode, double dis) |
|
{ |
|
IPointCollection pPointCol = passLine as IPointCollection; |
|
switch (mode) |
|
{ |
|
case 1: |
|
IPoint fromPoint = new PointClass(); |
|
passLine.QueryPoint(esriSegmentExtension.esriExtendAtFrom, -1 * dis, false, fromPoint); |
|
pPointCol.InsertPoints(0, 1, ref fromPoint); |
|
break; |
|
case 2: |
|
IPoint endPoint = new PointClass(); |
|
object missing = Type.Missing; |
|
passLine.QueryPoint(esriSegmentExtension.esriExtendAtTo, dis + passLine.Length, false, endPoint); |
|
pPointCol.AddPoint(endPoint, ref missing, ref missing); |
|
break; |
|
case 3: |
|
IPoint fPoint = new PointClass(); |
|
IPoint ePoint = new PointClass(); |
|
object missing2 = Type.Missing; |
|
passLine.QueryPoint(esriSegmentExtension.esriExtendAtFrom, -1 * dis, false, fPoint); |
|
pPointCol.InsertPoints(0, 1, ref fPoint); |
|
passLine.QueryPoint(esriSegmentExtension.esriExtendAtTo, dis + passLine.Length, false, ePoint); |
|
pPointCol.AddPoint(ePoint, ref missing2, ref missing2); |
|
break; |
|
default: |
|
return pPointCol as IPolyline; |
|
} |
|
return pPointCol as IPolyline; |
|
} |
|
/// <summary> |
|
/// 判断图形是否为三角形 |
|
/// </summary> |
|
/// <param name="geometry"></param> |
|
/// <returns></returns> |
|
private static bool isThreeSides(IGeometry geometry) |
|
{ |
|
bool isThreeSides = false; |
|
try |
|
{ |
|
if (geometry is IPolygon polygon1) |
|
{ |
|
IPointCollection pointCollection = polygon1 as IPointCollection; |
|
int vertexCount = pointCollection.PointCount; |
|
// 排除闭合多边形中重复的首尾顶点 |
|
IPoint startingpoint = pointCollection.get_Point(0);//起点 |
|
IPoint terminalpoint = pointCollection.get_Point(vertexCount - 1);//终点 |
|
if (polygon1.IsClosed && vertexCount > 1 && startingpoint.X == terminalpoint.X && startingpoint.Y == terminalpoint.Y) |
|
{ |
|
vertexCount--; |
|
} |
|
isThreeSides = vertexCount == 3; |
|
} |
|
} |
|
catch (Exception ex) |
|
{ |
|
|
|
} |
|
return isThreeSides; |
|
} |
|
/// <summary> |
|
/// 判断两个点是否共点（考虑浮点精度误差） |
|
/// </summary> |
|
/// <param name="p1">点1</param> |
|
/// <param name="p2">点2</param> |
|
/// <param name="tolerance">容差（默认1e-6，与坐标系单位一致）</param> |
|
/// <returns>true：共点；false：不共点</returns> |
|
private static bool ArePointsEqual(IPoint p1, IPoint p2, double tolerance = 1e-6) |
|
{ |
|
return Math.Abs(p1.X - p2.X) < tolerance && Math.Abs(p1.Y - p2.Y) < tolerance; |
|
} |
|
/// <summary> |
|
/// 计算两条线段的最短距离（支持不相交、相交等场景） |
|
/// </summary> |
|
/// <param name="a1">线段1的起点</param> |
|
/// <param name="a2">线段1的终点</param> |
|
/// <param name="b1">线段2的起点</param> |
|
/// <param name="b2">线段2的终点</param> |
|
/// <returns>两条线段的最短距离</returns> |
|
private static double CalculateSegmentDistance(IPoint a1, IPoint a2, IPoint b1, IPoint b2) |
|
{ |
|
// 向量定义：u为a线段方向，v为b线段方向，w为a1到b1的向量 |
|
double ux = a2.X - a1.X; |
|
double uy = a2.Y - a1.Y; |
|
double vx = b2.X - b1.X; |
|
double vy = b2.Y - b1.Y; |
|
double wx = a1.X - b1.X; |
|
double wy = a1.Y - b1.Y; |
|
double a = ux * ux + uy * uy; // u·u（a线段长度平方） |
|
double b = ux * vx + uy * vy; // u·v（方向向量点积） |
|
double c = vx * vx + vy * vy; // v·v（b线段长度平方） |
|
double d = ux * wx + uy * wy; // u·w |
|
double e = vx * wx + vy * wy; // v·w |
|
double denominator = a * c - b * b; // 分母：(u·u)(v·v)-(u·v)² |
|
double t, s; |
|
// 情况1：线段平行或共线（分母为0） |
|
if (Math.Abs(denominator) < 1e-9) |
|
{ |
|
t = 0; // a线段起点投影到b线段 |
|
s = (b * t + e) / c; // 计算s（若c=0则b线段为点） |
|
} |
|
// 情况2：线段不平行，计算参数t和s |
|
else |
|
{ |
|
t = (b * e - c * d) / denominator; // 投影参数t（a线段上） |
|
s = (a * e - b * d) / denominator; // 投影参数s（b线段上） |
|
} |
|
// 限制t和s在[0,1]范围内（超出则取端点） |
|
t = Math.Max(0, Math.Min(1, t)); |
|
s = Math.Max(0, Math.Min(1, s)); |
|
// 计算投影点距离 |
|
double dx = (wx + ux * t - vx * s); |
|
double dy = (wy + uy * t - vy * s); |
|
return Math.Sqrt(dx * dx + dy * dy); |
|
} |
|
private static double PointToSegmentDistance(IPoint p, IPoint segStart, IPoint segEnd) |
|
{ |
|
double segX = segEnd.X - segStart.X; |
|
double segY = segEnd.Y - segStart.Y; |
|
if (segX == 0 && segY == 0) // 线段为点 |
|
return Math.Sqrt(Math.Pow(p.X - segStart.X, 2) + Math.Pow(p.Y - segStart.Y, 2)); |
|
double t = ((p.X - segStart.X) * segX + (p.Y - segStart.Y) * segY) / (segX * segX + segY * segY); |
|
t = Math.Max(0, Math.Min(1, t)); // 投影参数限制在[0,1] |
|
double projX = segStart.X + t * segX; |
|
double projY = segStart.Y + t * segY; |
|
return Math.Sqrt(Math.Pow(p.X - projX, 2) + Math.Pow(p.Y - projY, 2)); |
|
} |
|
/// <summary> |
|
/// 预处理环的顶点坐标（结构体数组和IPoint数组） |
|
/// </summary> |
|
/// <param name="points">环的顶点集合（IPoint类型）</param> |
|
/// <param name="effectivePointCount">有效顶点数（闭合环需排除最后一个重复点）</param> |
|
/// <returns>元组：(顶点坐标结构体数组, IPoint对象数组)；若异常则返回(null, null)</returns> |
|
private static ((double X, double Y)[], IPoint[]) PreprocessRingPoints(IPointCollection points, int effectivePointCount) |
|
{ |
|
try |
|
{ |
|
var ringPoints = new (double X, double Y)[effectivePointCount]; // 存储顶点坐标（结构体，轻量级） |
|
var ringPointsIPoint = new IPoint[effectivePointCount]; // 存储IPoint对象（用于几何操作） |
|
// 遍历有效顶点，填充两个数组 |
|
for (int k = 0; k < effectivePointCount; k++) |
|
{ |
|
IPoint p = points.get_Point(k); // 获取第k个顶点 |
|
ringPoints[k] = (p.X, p.Y); // 转换为结构体 |
|
ringPointsIPoint[k] = p; // 保留IPoint对象 |
|
} |
|
return (ringPoints, ringPointsIPoint); |
|
} |
|
catch (Exception ex) |
|
{ |
|
// 异常处理：顶点读取失败时记录日志（实际开发中建议添加日志记录） |
|
return (null, null); |
|
} |
|
} |
|
/// <summary> |
|
/// 判断两个向量是否同向 |
|
/// </summary> |
|
/// <param name="vec1">向量1（如AB）</param> |
|
/// <param name="vec2">向量2（如CD）</param> |
|
/// <param name="tolerance">浮点数精度容差（默认1e-6）</param> |
|
/// <returns>是否同向</returns> |
|
public static bool AreVectorsSameDirection(Vector2D vec1, Vector2D vec2, double tolerance = 1e-6) |
|
{ |
|
// 1. 判断共线性：叉积的绝对值小于容差（考虑浮点数精度） |
|
double crossProduct = vec1.X * vec2.Y - vec2.X * vec1.Y; |
|
if (Math.Abs(crossProduct) > tolerance) |
|
return false; // 不共线，直接返回false |
|
// 2. 判断同向性：点积大于0（且向量模长不为0，避免零向量干扰） |
|
double dotProduct = vec1.X * vec2.X + vec1.Y * vec2.Y; |
|
double length1 = CalculateVectorLength(vec1); |
|
double length2 = CalculateVectorLength(vec2); |
|
// 模长为0的向量无方向（实际线段向量不会为零，此处为安全校验） |
|
if (length1 < tolerance || length2 < tolerance) |
|
return false; |
|
return dotProduct > tolerance; // 点积为正，方向相同 |
|
} |
|
/// <summary> |
|
/// 计算向量模长 |
|
/// </summary> |
|
public static double CalculateVectorLength(Vector2D vector) |
|
{ |
|
return Math.Sqrt(Math.Pow(vector.X, 2) + Math.Pow(vector.Y, 2)); |
|
} |
|
/// <summary> |
|
/// 计算线段向量 |
|
/// </summary> |
|
/// <param name="startPoint">线段起点</param> |
|
/// <param name="endPoint">线段终点</param> |
|
/// <returns>向量</returns> |
|
public static Vector2D CalculateVector(IPoint startPoint, IPoint endPoint) |
|
{ |
|
if (startPoint == null || endPoint == null) |
|
throw new ArgumentNullException("点对象不能为空"); |
|
|
|
return new Vector2D( |
|
endPoint.X - startPoint.X, |
|
endPoint.Y - startPoint.Y |
|
); |
|
} |
|
/// <summary> |
|
/// 计算两个向量的夹角（角度制，0°~180°） |
|
/// </summary> |
|
/// <param name="vec1">向量1</param> |
|
/// <param name="vec2">向量2</param> |
|
/// <param name="tolerance">精度容差（默认1e-6）</param> |
|
/// <returns>夹角角度（°）</returns> |
|
public static double CalculateAngle(Vector2D vec1, Vector2D vec2, double tolerance = 1e-6) |
|
{ |
|
// 计算模长 |
|
double len1 = CalculateVectorLength(vec1); |
|
double len2 = CalculateVectorLength(vec2); |
|
// 处理零向量（避免除零异常） |
|
if (len1 < tolerance || len2 < tolerance) |
|
throw new ArgumentException("向量模长不能为零"); |
|
// 计算点积 |
|
double dotProduct = vec1.X * vec2.X + vec1.Y * vec2.Y; |
|
// 计算cosθ（点积/(模长乘积)） |
|
double cosTheta = dotProduct / (len1 * len2); |
|
// 修正数值精度误差（cosθ范围需在[-1, 1]内） |
|
cosTheta = Math.Max(Math.Min(cosTheta, 1.0), -1.0); |
|
// 计算弧度并转换为角度 |
|
double radians = Math.Acos(cosTheta); |
|
return radians * (180 / Math.PI); // 弧度转角度 |
|
} |
|
|
|
// 定义向量结构体 |
|
public struct Vector2D |
|
{ |
|
public double X; |
|
public double Y; |
|
|
|
public Vector2D(double x, double y) |
|
{ |
|
X = x; |
|
Y = y; |
|
} |
|
} |
|
|
|
#region 节点去重 |
|
private static void RemoveDuplicatePointsInRing(IRing ring) |
|
{ |
|
IPointCollection pointColl = ring as IPointCollection; |
|
if (pointColl == null || pointColl.PointCount <= 1) return; |
|
// 步骤1：提取所有点到临时列表 |
|
List<IPoint> allPoints = new List<IPoint>(); |
|
for (int i = 0; i < pointColl.PointCount; i++) |
|
{ |
|
allPoints.Add(pointColl.get_Point(i)); |
|
} |
|
// 步骤2：全局去重（支持X/Y/Z/M全维度比较） |
|
HashSet<string> seenCoordinates = new HashSet<string>(); |
|
List<IPoint> uniquePoints = new List<IPoint>(); |
|
foreach (IPoint point in allPoints) |
|
{ |
|
// 生成坐标唯一键（根据需求调整是否包含Z/M） |
|
string coordKey = GetCoordinateKey(point); |
|
if (!seenCoordinates.Contains(coordKey)) |
|
{ |
|
seenCoordinates.Add(coordKey); |
|
uniquePoints.Add(point); |
|
} |
|
} |
|
// 步骤3：重建点集合（先清空再添加去重点） |
|
pointColl.RemovePoints(0, pointColl.PointCount); // 清空原始点 |
|
foreach (IPoint uniquePoint in uniquePoints) |
|
{ |
|
pointColl.AddPoint(uniquePoint); |
|
} |
|
// 步骤4：确保环闭合（最后一点与第一点一致） |
|
ring.Close(); |
|
} |
|
|
|
// 生成坐标唯一键（可根据需求调整比较维度） |
|
private static string GetCoordinateKey(IPoint point) |
|
{ |
|
// 仅比较X/Y（忽略Z/M） |
|
double x = point.X; |
|
double y = point.Y; |
|
return $"{x:F6},{y:F6}"; // 保留6位小数避免精度问题 |
|
} |
|
#endregion |
|
|
|
/// <summary> |
|
/// 判断图形是否存在狭长 |
|
/// </summary> |
|
/// <param name="geometry"></param> |
|
/// <returns></returns> |
|
private static string ISNarrow(IGeometry geometry) |
|
{ |
|
try |
|
{ |
|
var strgeometry = GeometryConvertHelper.ConvertIGeoemtryToWKT(geometry); |
|
NetTopologySuite.Geometries.Geometry ntsGeometry = (NetTopologySuite.Geometries.Geometry)new WKBReader().Read(GeometryConvertHelper.ConvertGeometryToWKB(geometry)); |
|
// 内缩参数（侵蚀） |
|
var innerBufferParams = new BufferParameters |
|
{ |
|
QuadrantSegments = 30, // 内缩平滑度 |
|
JoinStyle = JoinStyle.Bevel, // 斜角连接避免尖锐残留 |
|
EndCapStyle = EndCapStyle.Square // 平端帽处理开放线段 |
|
}; |
|
var ntsGeometry1 = ntsGeometry.Buffer(-0.05, innerBufferParams); // 内缩0.1 |
|
strgeometry = ntsGeometry1.AsText(); |
|
if (strgeometry == "POLYGON EMPTY") |
|
{ |
|
return strgeometry; |
|
} |
|
// 外扩参数（膨胀） |
|
var outerBufferParams = new BufferParameters |
|
{ |
|
QuadrantSegments = 200, // 高平滑度外扩 |
|
JoinStyle = JoinStyle.Mitre, |
|
MitreLimit = 2 // 限制尖锐拐角延伸长度 |
|
}; |
|
var ntsGeometry2 = ntsGeometry1.Buffer(0.05, outerBufferParams); // 外扩0.1 |
|
strgeometry = ntsGeometry2.AsText(); |
|
double length1 = ntsGeometry.Length; |
|
int pointCount1 = ntsGeometry.NumPoints; |
|
double length2 = ntsGeometry2.Length; |
|
int pointCount2 = ntsGeometry2.NumPoints; |
|
double delta_length = length1 - length2; |
|
double avg_halfangle = 180 * (pointCount1 - 1 - 2) / (pointCount1 - 1) / 2; |
|
double conner_normal_length = 2 * 0.05 / Math.Tan(avg_halfangle * (Math.PI / 180)); |
|
if (delta_length > 8 * conner_normal_length * (pointCount1 - 1)) |
|
{ |
|
return strgeometry; |
|
} |
|
} |
|
catch (Exception ex) |
|
{ |
|
} |
|
return string.Empty; |
|
} |
|
|
|
} |
|
|
|
public class EarClipper |
|
{ |
|
// 三角剖分结果（每个三角形为3个点的集合） |
|
public List<List<IPoint>> Triangles { get; private set; } |
|
|
|
public EarClipper() |
|
{ |
|
Triangles = new List<List<IPoint>>(); |
|
} |
|
// 多边形整体方向（顺时针/逆时针） |
|
private bool _isCounterClockwise; |
|
/// <summary> |
|
/// 对多边形执行耳切法三角剖分 |
|
/// </summary> |
|
/// <param name="polygon">输入多边形（需为简单多边形，无自交）</param> |
|
/// <returns>三角剖分结果</returns> |
|
public bool Clip(IPolygon polygon) |
|
{ |
|
if (polygon == null || polygon.IsEmpty) |
|
return false; |
|
// 提取多边形外环的顶点（简单多边形默认处理外环） |
|
IPointCollection pointCollection = polygon as IPointCollection; |
|
if (pointCollection == null || pointCollection.PointCount < 3) |
|
return false; |
|
// 转换为顶点列表（排除最后一个闭合点） |
|
List<IPoint> vertices = new List<IPoint>(); |
|
for (int i = 0; i < pointCollection.PointCount - 1; i++) |
|
{ |
|
IPoint p = pointCollection.get_Point(i); |
|
vertices.Add(p); |
|
} |
|
// 关键：计算多边形整体方向（顺时针/逆时针） |
|
_isCounterClockwise = IsPolygonCounterClockwise(vertices); |
|
// 执行耳切法 |
|
Triangles.Clear(); |
|
ClipEars(vertices); |
|
return Triangles.Count > 0; |
|
} |
|
|
|
/// <summary> |
|
/// 迭代切割耳 |
|
/// </summary> |
|
private void ClipEars(List<IPoint> vertices) |
|
{ |
|
if (vertices.Count < 3) return; |
|
int index = 0; |
|
while (vertices.Count > 3) |
|
{ |
|
// 循环取顶点（处理索引越界） |
|
int prevIndex = (index - 1 + vertices.Count) % vertices.Count; |
|
int currIndex = index % vertices.Count; |
|
int nextIndex = (index + 1) % vertices.Count; |
|
IPoint prev = vertices[prevIndex]; |
|
IPoint curr = vertices[currIndex]; |
|
IPoint next = vertices[nextIndex]; |
|
// 检查当前三点是否构成"耳" |
|
if (IsConvex(prev, curr, next) && IsEar(prev, curr, next, vertices)) |
|
{ |
|
// 记录三角形 |
|
Triangles.Add(new List<IPoint> { prev, curr, next }); |
|
// 移除当前顶点（切割耳） |
|
vertices.RemoveAt(currIndex); |
|
// 重置索引（从下一个顶点开始检查） |
|
index = 0; |
|
} |
|
else |
|
{ |
|
index++; |
|
// 防止死循环（理论上简单多边形必然有耳） |
|
if (index >= vertices.Count * 2) |
|
break; |
|
} |
|
} |
|
// 添加最后一个三角形 |
|
if (vertices.Count == 3) |
|
{ |
|
Triangles.Add(new List<IPoint> { vertices[0], vertices[1], vertices[2] }); |
|
} |
|
} |
|
|
|
/// <summary> |
|
/// 判断顶点curr是否为凸点（通过叉积） |
|
/// </summary> |
|
private bool IsConvex(IPoint prev, IPoint curr, IPoint next) |
|
{ |
|
// 计算向量：prev->curr 和 curr->next |
|
double v1x = curr.X - prev.X; |
|
double v1y = curr.Y - prev.Y; |
|
double v2x = next.X - curr.X; |
|
double v2y = next.Y - curr.Y; |
|
// 叉积（z分量）：v1 × v2 = v1x*v2y - v1y*v2x |
|
double cross = v1x * v2y - v1y * v2x; |
|
// 叉积≥0：逆时针多边形的凸点（根据多边形方向调整符号） |
|
// 关键：根据多边形整体方向判断凸点 |
|
if (_isCounterClockwise) |
|
{ |
|
// 逆时针多边形：叉积 >= 0 为凸点 |
|
return cross >= 0; |
|
} |
|
else |
|
{ |
|
// 顺时针多边形：叉积 <= 0 为凸点（与逆时针相反） |
|
return cross <= 0; |
|
} |
|
} |
|
/// <summary> |
|
/// 判断多边形顶点是顺时针还是逆时针排列 |
|
/// </summary> |
|
private bool IsPolygonCounterClockwise(List<IPoint> vertices) |
|
{ |
|
double sum = 0; |
|
int n = vertices.Count; |
|
for (int i = 0; i < n; i++) |
|
{ |
|
IPoint p = vertices[i]; |
|
IPoint q = vertices[(i + 1) % n]; // 下一个顶点（首尾相连） |
|
sum += (q.X - p.X) * (q.Y + p.Y); |
|
} |
|
// sum > 0 表示逆时针，sum < 0 表示顺时针 |
|
return sum > 0; |
|
} |
|
|
|
/// <summary> |
|
/// 判断三点组成的三角形是否为"耳"（内部无其他顶点） |
|
/// </summary> |
|
private bool IsEar(IPoint a, IPoint b, IPoint c, List<IPoint> vertices) |
|
{ |
|
// 检查所有其他顶点是否在三角形abc内部 |
|
foreach (IPoint p in vertices) |
|
{ |
|
// 跳过a、b、c本身 |
|
if (p.Equals(a) || p.Equals(b) || p.Equals(c)) |
|
continue; |
|
// 判断点p是否在三角形abc内部 |
|
if (IsPointInTriangle(p, a, b, c)) |
|
return false; |
|
} |
|
return true; |
|
} |
|
|
|
/// <summary> |
|
/// 判断点p是否在三角形abc内部（ barycentric坐标法） |
|
/// </summary> |
|
private bool IsPointInTriangle(IPoint p, IPoint a, IPoint b, IPoint c) |
|
{ |
|
double v0x = c.X - a.X; |
|
double v0y = c.Y - a.Y; |
|
double v1x = b.X - a.X; |
|
double v1y = b.Y - a.Y; |
|
double v2x = p.X - a.X; |
|
double v2y = p.Y - a.Y; |
|
// 计算点积 |
|
double dot00 = v0x * v0x + v0y * v0y; |
|
double dot01 = v0x * v1x + v0y * v1y; |
|
double dot02 = v0x * v2x + v0y * v2y; |
|
double dot11 = v1x * v1x + v1y * v1y; |
|
double dot12 = v1x * v2x + v1y * v2y; |
|
// 计算 barycentric坐标 |
|
double invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01); |
|
double u = (dot11 * dot02 - dot01 * dot12) * invDenom; |
|
double v = (dot00 * dot12 - dot01 * dot02) * invDenom; |
|
// 点在三角形内（包括边界） |
|
return (u >= 0) && (v >= 0) && (u + v <= 1); |
|
} |
|
} |
|
} |