1.向量Vector3d
1 using System;
2
3 namespace RGeos.Geometry
4 {
5 /// <summary>
6 /// 3D向量类
7 /// </summary>
8 public class Vector3d
9 {
10 public double[] vector;
11 private const double E = 0.0000001f;
12 /// <summary>
13 ///
14 /// </summary>
15 /// <param name="x"></param>
16 /// <param name="y"></param>
17 /// <param name="z"></param>
18 ///
19 public Vector3d()
20 {
21 vector = new double[3];
22 }
23 public Vector3d(double x, double y, double z)
24 {
25 vector = new double[3] { x, y, z };
26 }
27 public Vector3d(Vector3d vct)
28 {
29 vector = new double[3];
30 vector[0] = vct.X;
31 vector[1] = vct.Y;
32 vector[2] = vct.Z;
33 }
34 #region 属性
35 /// <summary>
36 /// X向量
37 /// </summary>
38 public double X
39 {
40 get { return vector[0]; }
41 set { vector[0] = value; }
42 }
43 /// <summary>
44 /// Y向量
45 /// </summary>
46 public double Y
47 {
48 get { return vector[1]; }
49 set { vector[1] = value; }
50 }
51 /// <summary>
52 /// Z向量
53 /// </summary>
54 public double Z
55 {
56 get { return vector[2]; }
57 set { vector[2] = value; }
58 }
59 #endregion
60
61 #region 向量操作
62 /// <summary>
63 /// /// <summary>
64 /// 向量加法+
65 /// </summary>
66 /// <param name="lhs"></param>
67 /// <param name="rhs"></param>
68 /// <returns></returns>
69 public static Vector3d operator +(Vector3d lhs, Vector3d rhs)//向量加
70 {
71 Vector3d result = new Vector3d(lhs);
72 result.X += rhs.X;
73 result.Y += rhs.Y;
74 result.Z += rhs.Z;
75 return result;
76 }
77 /// <summary>
78 /// 向量减-
79 /// </summary>
80 /// <param name="lhs"></param>
81 /// <param name="rhs"></param>
82 /// <returns></returns>
83 public static Vector3d operator -(Vector3d lhs, Vector3d rhs)//向量减法
84 {
85 Vector3d result = new Vector3d(lhs);
86 result.X -= rhs.X;
87 result.Y -= rhs.Y;
88 result.Z -= rhs.Z;
89 return result;
90 }
91 /// <summary>
92 /// 向量除
93 /// </summary>
94 /// <param name="lhs"></param>
95 /// <param name="rhs"></param>
96 /// <returns></returns>
97 public static Vector3d operator /(Vector3d lhs, double rhs)//向量除以数量
98 {
99 if (rhs != 0)
100 return new Vector3d(lhs.X / rhs, lhs.Y / rhs, lhs.Z / rhs);
101 else
102 return new Vector3d(0, 0, 0);
103 }
104 /// <summary>
105 /// 向量数乘*
106 /// </summary>
107 /// <param name="lhs"></param>
108 /// <param name="rhs"></param>
109 /// <returns></returns>
110 public static Vector3d operator *(double lhs, Vector3d rhs)//左乘数量
111 {
112 return new Vector3d(lhs * rhs.X, lhs * rhs.Y, lhs * rhs.Z);
113 }
114 /// <summary>
115 /// 向量数乘
116 /// </summary>
117 /// <param name="lhs"></param>
118 /// <param name="rhs"></param>
119 /// <returns></returns>
120 public static Vector3d operator *(Vector3d lhs, double rhs)//右乘数量
121 {
122 return new Vector3d(lhs.X * rhs, lhs.Y * rhs, lhs.Z * rhs);
123 }
124
125 /// <summary>
126 /// 判断量向量是否相等
127 /// </summary>
128 /// <param name="lhs"></param>
129 /// <param name="rhs"></param>
130 /// <returns>True 或False</returns>
131 public static bool operator ==(Vector3d lhs, Vector3d rhs)
132 {
133 if (Math.Abs(lhs.X - rhs.X) < E && Math.Abs(lhs.Y - rhs.Y) < E && Math.Abs(lhs.Z - rhs.Z) < E)
134 return true;
135 else
136 return false;
137 }
138 public static bool operator !=(Vector3d lhs, Vector3d rhs)
139 {
140 return !(lhs == rhs);
141 }
142 public override bool Equals(object obj)
143 {
144 return base.Equals(obj);
145 }
146 public override int GetHashCode()
147 {
148 return base.GetHashCode();
149 }
150 public override string ToString()
151 {
152 return "(" + X + "," + Y + "," + Z + ")";
153 }
154 /// <summary>
155 /// 向量叉积，求与两向量垂直的向量
156 /// </summary>
157 public static Vector3d Cross(Vector3d v1, Vector3d v2)
158 {
159 Vector3d r = new Vector3d(0, 0, 0);
160 r.X = (v1.Y * v2.Z) - (v1.Z * v2.Y);
161 r.Y = (v1.Z * v2.X) - (v1.X * v2.Z);
162 r.Z = (v1.X * v2.Y) - (v1.Y * v2.X);
163 return r;
164 }
165 /// <summary>
166 /// 向量数量积
167 /// </summary>
168 /// <param name="lhs"></param>
169 /// <param name="rhs"></param>
170 /// <returns></returns>
171 public static double operator *(Vector3d lhs, Vector3d rhs)//
172 {
173 return lhs.X * rhs.X + lhs.Y * rhs.Y + lhs.Z * rhs.Z;
174 }
175 /// <summary>
176 /// 内积
177 /// </summary>
178 /// <param name="v1"></param>
179 /// <param name="v2"></param>
180 /// <returns></returns>
181 public static double InnerMultiply(Vector3d v1, Vector3d v2)
182 {
183 double inner = 0.0;
184 inner = v1.X * v2.X + v1.Y * v2.Y + v1.Z * v2.Z;
185 return inner;
186 }
187 /// <summary>
188 /// 求向量长度，向量的模
189 /// </summary>
190 public static double Magnitude(Vector3d v1)
191 {
192 return (double)Math.Sqrt((v1.X * v1.X) + (v1.Y * v1.Y) + (v1.Z * v1.Z));
193 }
194 /// <summary>
195 /// 单位化向量
196 /// </summary>
197 public static Vector3d Normalize(Vector3d v1)
198 {
199 double magnitude = Magnitude(v1);
200 v1 = v1 / magnitude;
201 return v1;
202 }
203 #endregion
204 }
205 }
2. 计算基础
1 using System; 2 using RGeos.Geometry; 3 4 namespace RGeos.Basic 5 { 6 public class RMath 7 { 8 public static double SMALL_NUM = 0.0000000001; // anything that avoids division overflow 9 /// <summary> 10 /// dot product (3D) which allows vector operations in arguments 11 /// </summary> 12 /// <param name="u"></param> 13 /// <param name="v"></param> 14 /// <returns></returns> 15 public static double dot(Vector3d u, Vector3d v) 16 { 17 return ((u).X * (v).X + (u).Y * (v).Y + (u).Z * (v).Z); 18 } 19 /// <summary> 20 /// 2D数量积，点乘 21 /// </summary> 22 /// <param name="u"></param> 23 /// <param name="v"></param> 24 /// <returns></returns> 25 public static double dot2(Vector3d u, Vector3d v) 26 { 27 return ((u).X * (v).X + (u).Y * (v).Y); 28 } 29 /// <summary> 30 /// 2D矢量叉积，定义为(0,0),P1,P2和P1P2包围四边形的带符号面积 31 /// </summary> 32 /// <param name="u"></param> 33 /// <param name="v"></param> 34 /// <returns></returns> 35 public static double perp(Vector3d u, Vector3d v) 36 { 37 return ((u).X * (v).Y - (u).Y * (v).X); 38 } 39 /// <summary> 40 /// 向量的模 41 /// </summary> 42 /// <param name="v"></param> 43 /// <returns></returns> 44 public static double norm(Vector3d v) 45 { 46 return Math.Sqrt(dot(v, v)); // norm = length of vector 47 } 48 /// <summary> 49 /// 两点间距离 50 /// </summary> 51 /// <param name="u"></param> 52 /// <param name="v"></param> 53 /// <returns></returns> 54 public static double d(Vector3d u, Vector3d v) 55 { 56 return norm(u - v); // distance = norm of difference 57 } 58 59 public static double d(RPoint P1, RPoint P2) 60 { 61 return GetDistance(P1, P2); // distance = norm of difference 62 } 63 64 // 判断点P2在直线P0P1的左边还是在右边，还是在直线上 65 //isLeft(): tests if a point is Left|On|Right of an infinite line. 66 // Input: three points P0, P1, and P2 67 // Return: >0 for P2 left of the line through P0 and P1 68 // =0 for P2 on the line 69 // <0 for P2 right of the line 70 public static int isLeft(RPoint P0, RPoint P1, RPoint P2) 71 { 72 double l = ((P1.X - P0.X) * (P2.Y - P0.Y) - (P2.X - P0.X) * (P1.Y - P0.Y)); 73 return (int)l; 74 } 75 /// <summary> 76 /// 获取由两个点所形成的向量的象限角度 77 /// </summary> 78 /// <param name="preCoord">第一个点的坐标</param> 79 /// <param name="nextCoord">第二个点的坐标</param> 80 /// <returns></returns> 81 public static double GetQuadrantAngle(RPoint preCoord, RPoint nextCoord) 82 { 83 return GetQuadrantAngle(nextCoord.X - preCoord.X, nextCoord.Y - preCoord.Y); 84 } 85 /// <summary> 86 /// 由增量X和增量Y所形成的向量的象限角度 87 /// 区别方位角：方位角以正北方向顺时针 88 /// | | 89 /// | / |b / 90 /// | / a |^/ 91 /// |/_)____象限角 |/______方位角 92 /// </summary> 93 /// <param name="x">增量X</param> 94 /// <param name="y">增量Y</param> 95 /// <returns>象限角</returns> 96 public static double GetQuadrantAngle(double x, double y) 97 { 98 double theta = Math.Atan(y / x); 99 if (x > 0 && y == 0) return 0; 100 if (x == 0 && y > 0) return Math.PI / 2; 101 if (x < 0 && y == 0) return Math.PI; 102 if (x == 0 && y < 0) return 3 * Math.PI / 2; 103 104 if (x > 0 && y > 0) return theta; 105 if (x > 0 && y < 0) return Math.PI * 2 + theta; 106 if (x < 0 && y > 0) return theta + Math.PI; 107 if (x < 0 && y < 0) return theta + Math.PI; 108 return theta; 109 } 110 /// <summary> 111 /// 获取由相邻的三个点A-B-C所形成的两个向量之间的夹角 112 /// 向量AB，BC形成的夹角 113 /// </summary> 114 /// <param name="preCoord">第一个点</param> 115 /// <param name="midCoord">中间点</param> 116 /// <param name="nextCoord">第三个点</param> 117 /// <returns></returns> 118 public static double GetIncludedAngle(RPoint preCoord, RPoint midCoord, RPoint nextCoord) 119 { 120 double innerProduct = (midCoord.X - preCoord.X) * (nextCoord.X - midCoord.X) + (midCoord.Y - preCoord.Y) * (nextCoord.Y - midCoord.Y); 121 double mode1 = Math.Sqrt(Math.Pow((midCoord.X - preCoord.X), 2.0) + Math.Pow((midCoord.Y - preCoord.Y), 2.0)); 122 double mode2 = Math.Sqrt(Math.Pow((nextCoord.X - midCoord.X), 2.0) + Math.Pow((nextCoord.Y - midCoord.Y), 2.0)); 123 return Math.Acos(innerProduct / (mode1 * mode2)); 124 } 125 /// <summary> 126 /// 获取由两个点所形成的向量的模(长度) 127 /// </summary> 128 /// <param name="preCoord">第一个点</param> 129 /// <param name="nextCoord">第二个点</param> 130 /// <returns>由两个点所形成的向量的模(长度)</returns> 131 public static double GetDistance(RPoint preCoord, RPoint nextCoord) 132 { 133 return Math.Sqrt(Math.Pow((nextCoord.X - preCoord.X), 2) + Math.Pow((nextCoord.Y - preCoord.Y), 2)); 134 } 135 } 136 }