1 using System;
2 using System.Collections.Generic;
3 using System.Linq;
4 using System.Text;
5 using ESRI.ArcGIS.Geometry;
6 using RGeos.Geometry
7
8 namespace RGeos.Geometry
9 {
10 public class CulmulateDistance
11 {
12 /// <summary>
13 /// 点到折线最短距离处距离折线起点的累积距离
14 /// </summary>
15 /// <param name="P">任意一点</param>
16 /// <param name="polyline">折线</param>
17 /// <returns></returns>
18 public static double CulmulateDist_Point_to_Polyline(IPoint P, IPolyline polyline)
19 {
20 ISegmentCollection segs = polyline as ISegmentCollection;
21 double min = double.MaxValue;
22 int segIndex = -1;//最短距离所在档的索引
23 for (int i = 0; i < segs.SegmentCount; i++)
24 {
25 //点到每条线段的最短距离
26 double dis = Dist_Point_to_Segment(P, segs.get_Segment(i));
27 if (dis < min)
28 {
29 min = dis;
30 segIndex = i;//取出最小的一个
31 }
32 }
33 double culmulateDis = 0;
34 for (int i = 0; i < segs.SegmentCount; i++)
35 {
36 if (segIndex != i)
37 {
38 culmulateDis += segs.get_Segment(i).Length;
39 }
40 else
41 {
42 ISegment current = segs.get_Segment(i);
43 Vector3d v = new Vector3d();
44 v.X = current.ToPoint.X - current.FromPoint.X;
45 v.Y = current.ToPoint.Y - current.FromPoint.Y;
46
47 Vector3d w = new Vector3d();
48 w.X = P.X - current.FromPoint.X;
49 w.Y = P.Y - current.FromPoint.Y;
50
51 double c1 = dot(w, v);//投影长度
52 if (c1 <= 0)//这种情况最短距离在该线段的起点处
53 {
54 break;
55 }
56 double c2 = dot(v, v);
57 if (c2 <= c1)//这种情况最短距离在该线段的终点处
58 {
59 culmulateDis += Math.Sqrt(c2);
60 break;
61 }
62
63 double b = c1 / c2;
64 culmulateDis += Math.Sqrt(c1);
65 IPoint Pb = new PointClass();
66 Pb.X = current.FromPoint.X + b * v.X;
67 Pb.Y = current.FromPoint.Y + b * v.Y;
68 break;
69 }
70 }
71 return culmulateDis;
72 }
73
74 public static IPoint GetCentrePoint(IPolygon polygon)
75 {
76 IArea pArea = polygon as IArea;
77 IPoint pt = new PointClass();
78 pt.X = pArea.Centroid.X;
79 pt.Y = pArea.Centroid.Y;
80 return pt;
81 }
82
83 public static double Dist_Point_to_Segment(IPoint P, ISegment S)
84 {
85 Vector3d v = new Vector3d();
86 v.X = S.ToPoint.X - S.FromPoint.X;
87 v.Y = S.ToPoint.Y - S.FromPoint.Y;
88
89 Vector3d w = new Vector3d();
90 w.X = P.X - S.FromPoint.X;
91 w.Y = P.Y - S.FromPoint.Y;
92
93 double c1 = dot(w, v);
94 if (c1 <= 0)
95 return d(P, S.FromPoint);
96
97 double c2 = dot(v, v);
98 if (c2 <= c1)
99 return d(P, S.ToPoint);
100
101 double b = c1 / c2;
102 IPoint Pb = new PointClass();
103 Pb.X = S.FromPoint.X + b * v.X;
104 Pb.Y = S.FromPoint.Y + b * v.Y;
105 return d(P, Pb);
106 }
107 /// <summary>
108 /// 向量的模
109 /// </summary>
110 /// <param name="v"></param>
111 /// <returns></returns>
112 public static double norm(Vector3d v)
113 {
114 return Math.Sqrt(dot2(v, v)); // norm = length of vector
115 }
116 /// <summary>
117 /// 2D数量积，点乘
118 /// </summary>
119 /// <param name="u"></param>
120 /// <param name="v"></param>
121 /// <returns></returns>
122 public static double dot2(Vector3d u, Vector3d v)
123 {
124 return ((u).X * (v).X + (u).Y * (v).Y);
125 }
126 public static double dot(Vector3d u, Vector3d v)
127 {
128 return ((u).X * (v).X + (u).Y * (v).Y + (u).Z * (v).Z);
129 }
130
131 public static double d(IPoint P, IPoint P1)
132 {
133 return Math.Sqrt((P1.X - P.X) * (P1.X - P.X) + (P1.Y - P.Y) * (P1.Y - P.Y));
134 }
135
136 /// <param name="x">增量X</param>
137 /// <param name="y">增量Y</param>
138 /// <returns>象限角</returns>
139 public static double GetQuadrantAngle(double x, double y)
140 {
141 double theta = Math.Atan(y / x);
142 if (x > 0 && y == 0) return 0;
143 if (x == 0 && y > 0) return Math.PI / 2;
144 if (x < 0 && y == 0) return Math.PI;
145 if (x == 0 && y < 0) return 3 * Math.PI / 2;
146
147 if (x > 0 && y > 0) return theta;
148 if (x > 0 && y < 0) return Math.PI * 2 + theta;
149 if (x < 0 && y > 0) return theta + Math.PI;
150 if (x < 0 && y < 0) return theta + Math.PI;
151 return theta;
152 }
153 }
154 }
时间： 2025-01-15 14:06:03