CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~
二维
点和向量化为坐标 Coord 进行运算，使用 stl 中的 complex 实现。
复数相乘的几何意义为长度相乘，极角相加。
用直线上的一点 p 和方向向量 v 表示一条经过 p 的直线，直线上的所有点 q 满足 q=p+t*v，其中 t 是参数；当限制 t≥0 时，该参数方程表示射线；限制 0≤t≤1 时，该参数方程表示线段。
此外，如果已知线段端点 a1 和 a2，可以通过 Line(a1,a2-a1)来得到对应的参数形式。
Morley 定理：三角形每个内角的三等分线相交成等边三角形。
欧拉定理：平面图的点数 V、边数 E 和面数 F 满足 V+F-E=2。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
typedef double lf;
typedef complex<lf> Coord;
const lf EPS = 1e-9;
#define X real()
#define Y imag()
struct Line
{
Coord p, v;
Line(Coord p = Coord(), Coord v = Coord()) : p(p), v(v) {}
Coord point(lf t) const { return p + v * t; }
};
struct Circle
{
Coord c;
lf r;
Circle(Coord c = Coord(), lf r = 0) : c(c), r(r) {}
Coord point(lf t) const { return c + polar(r, t); } //t为参数，幅角
};
/*
Coord(lf x=0,lf y=0);//构造函数
lf real(Coord a);//a的实部（复平面的横坐标）,也可写作a.real()
lf imag(Coord a);//a的虚部（复平面的纵坐标）,也可写作a.imag()
lf abs(Coord a);//向量a的模长，或是点a到原点的距离
lf norm(Coord a);//abs的平方，比abs快，但是要注意浮点数精度溢出
lf arg(Coord a);//a的幅角，与atan2(a.real(),a.imag())等价
Coord polar(lf r,lf t);//极坐标生成方式，r为幅值，t为幅角
//运算符重载+、-、*、/（以及对应的赋值运算，但是赋值运算不能写在表达式中，详见参考地址）、<<、>>（输出括号形式的坐标）
*/
int sgn(lf d) { return (d > EPS) - (d < -EPS); }
bool operator!=(const Coord &A, const Coord &B) { return sgn(A.X - B.X) || sgn(A.Y - B.Y); } //不等运算符，涉及到浮点数比较要重写
bool operator==(const Coord &A, const Coord &B) { return !(A != B); }
bool cmpCoord(const Coord &A, const Coord &B) { return sgn(A.X - B.X) ? sgn(A.X - B.X) < 0 : sgn(A.Y - B.Y) < 0; } //复数没有小于运算，只能这样定义一个比较函数
bool cmpLine(const Line &A, const Line &B) { return sgn(arg(A.v) - arg(B.v)) < 0; } //按极角排序，求凸包中使用
lf Dot(const Coord &A, const Coord &B) { return A.X * B.X + A.Y * B.Y; }
lf Cross(const Coord &A, const Coord &B) { return A.X * B.Y - B.X * A.Y; }
lf Angle(const Coord &A, const Coord &B) { return acos(Dot(A, B) / abs(A) / abs(B)); }
lf Area2(const Coord &A, const Coord &B, const Coord &C) { return Cross(B - A, C - A); } //三角形ABC有向面积的两倍
Coord Rotate(const Coord &A, lf rad) { return A * polar(1.0, rad);} //向量A逆时针旋转rad弧度
Coord Normal(const Coord &A) //A的法向量，把A逆时针旋转九十度并长度化为1
{
lf L = abs(A);
return Coord(-A.Y / L, A.X / L);
}
bool onLeft(const Coord &P, const Line &L) { return sgn(Cross(L.v, P - L.p)) > 0; } //p是否在有向直线L左侧，不含线上
lf DistanceToLine(const Coord &P, const Line &L) { return Cross(L.v, P - L.p) / abs(L.v); } //点到直线距离（有向）
lf DistanceToLine(const Coord &P, const Coord &A, const Coord &B) { return DistanceToLine(P, Line(A, B - A)); }
lf DistanceToSegment(const Coord &P, const Coord &A, const Coord &B) //点到线段的距离（无向）
{
if (A == B)
return abs(P - A);
Coord v1 = B - A, v2 = P - A, v3 = P - B;
if (sgn(Dot(v1, v2)) < 0)
return abs(v2);
if (sgn(Dot(v1, v3)) > 0)
return abs(v3);
return fabs(DistanceToLine(P, Line(A, B - A)));
}
Coord getLineProjection(const Coord &P, const Line &L) { return L.point(Dot(L.v, P - L.p) / norm(L.v)); } //点在直线上的投影
Coord getLineProjection(const Coord &P, const Coord &A, const Coord &B) { return getLineProjection(P, Line(A, B - A)); }
Coord getSymmetry(const Coord &P, const Coord &O) { return O + O - P; } //P关于O的对称点
Coord getSymmetry(const Coord &P, const Line &L) { return getSymmetry(P, getLineProjection(P, L)); } //P关于L的对称点
Coord getLineIntersection(const Line &L1, const Line &L2) { return L1.point(Cross(L2.v, L1.p - L2.p) / Cross(L1.v, L2.v)); } //直线交点,须确保两直线相交
Coord getLineIntersection(const Coord &A1, const Coord &A2, const Coord &B1, const Coord &B2) { return getLineIntersection(Line(A1, A2 - A1), Line(B1, B2 - B1)); }
bool SegmentProperIntersection(const Coord &A1, const Coord &A2, const Coord &B1, const Coord &B2) //线段相交判定，交点不在一条线段的端点
{
lf C1 = Cross(A2 - A1, B1 - A1), C2 = Cross(A2 - A1, B2 - A1),
C3 = Cross(B2 - B1, A1 - B1), C4 = Cross(B2 - B1, A2 - B1);
return sgn(C1) * sgn(C2) < 0 && sgn(C3) * sgn(C4) < 0;
}
bool onSegment(const Coord &P, const Coord &A1, const Coord &A2) { return sgn(Dot(A1 - P, A2 - P)) < 0 && !sgn(Cross(A1 - P, A2 - P)); } //判断点是否在线段上，不包含端点
lf PolygonArea(const vector<Coord> &p) //计算多边形的有向面积，凸多边形即为面积
{
lf s = 0;
for (int i = 2; i < p.size(); ++i)
s += Area2(p[0], p[i - 1], p[i]);
return s / 2;
}
int inPolygon(const Coord &p, const vector<Coord> &poly) //点在多边形内的判定，转角法，正值为内部，0为外部，-1在边界上
{
int ans = 0;
for (int i = 0, k, d1, d2, n = poly.size(); i != n; ++i)
{
if (onSegment(p, poly[i], poly[(i + 1) % n]))
return -1; //在边界上
k = sgn(Cross(poly[(i + 1) % n] - poly[i], p - poly[i]));
d1 = sgn(poly[i].Y - p.Y);
d2 = sgn(poly[(i + 1) % n].Y - p.Y);
if (k > 0 && d1 <= 0 && d2 > 0)
++ans;
if (k < 0 && d2 <= 0 && d1 > 0)
--ans;
}
return ans;
}
vector<Coord> ConvexHull(vector<Coord> p, int collineation = 1) //获得凸包，不希望凸包的边上有输入点第二个参数传0
{
vector<Coord> ans;
sort(p.begin(), p.end(), cmpCoord); //先比横坐标再比纵坐标
for (int i = 0; i < p.size(); ++i) //求出下凸包
{
while (ans.size() > 1 && sgn(Area2(ans[ans.size() - 2], ans[ans.size() - 1], p[i])) < collineation)
ans.pop_back();
ans.push_back(p[i]);
}
for (int i = p.size() - 2, k = ans.size(); i >= 0; --i) //求出上凸包
{
while (ans.size() > k && -sgn(Area2(ans[ans.size() - 1], ans[ans.size() - 2], p[i])) < collineation)
ans.pop_back();
ans.push_back(p[i]);
}
if (p.size() > 1)
ans.pop_back();
return ans;
}
vector<Coord> cutPolygon(const vector<Coord> &poly, const Coord &A, const Coord &B) //用有向直线A->B切割多边形poly， 返回「左侧」。 如果退化，可能会返回一个单点或者线段，复杂度O(n^2)
{
vector<Coord> newpoly;
for (int i = 0, n = poly.size(); i != n; ++i)
{
Coord C = poly[i], D = poly[(i + 1) % n];
if (sgn(Cross(B - A, C - A)) >= 0)
newpoly.push_back(C);
if (!sgn(Cross(B - A, C - D)))
{
Coord ip = getLineIntersection(Line(A, B - A), Line(C, D - C));
if (onSegment(ip, C, D))
newpoly.push_back(ip);
}
}
return newpoly;
}
vector<Coord> getHalfPlaneIntersection(vector<Line> L) //半平面交
{
sort(L.begin(), L.end(), cmpLine); //按极角排序
vector<Coord> p(L.size(), Coord()); //p[i]为q[i]和q[i+1]的交点
int first = 0, last = 0; //双端队列的第一个元素和最后一个元素
vector<Line> q(L.size(), Line()); //双端队列
q[0] = L[0]; //队列初始化为只有一个半平面L[0]
for (int i = 0, n = L.size(); i != n; ++i)
{
while (first < last && !onLeft(p[last - 1], L[i]))
--last;
while (first < last && !onLeft(p[first], L[i]))
++first;
q[++last] = L[i];
if (!sgn(Cross(q[last].v, q[last - 1].v)))
{
--last;
if (onLeft(L[i].p, q[last]))
q[last] = L[i];
}
if (first < last)
p[last - 1] = getLineIntersection(q[last - 1], q[last]);
}
while (first < last && !onLeft(p[last - 1], q[first]))
--last; //删除无用平面
if (last - first <= 1)
return vector<Coord>(); //空集
p[last] = getLineIntersection(q[last], q[first]);
return vector<Coord>(p.begin() + first, p.begin() + last + 1); //从deque复制到输出中
}
int getLineCircleIntersection(const Line &L, const Circle &C, vector<Coord> &sol)
{
lf a = L.v.X,
b = L.p.X - C.c.X,
c = L.v.Y,
d = L.p.Y - C.c.Y,
e = a * a + c * c,
f = 2 * (a * b + c * d),
g = b * b + d * d - C.r * C.r,
delta = f * f - 4 * e * g;
if (sgn(delta) < 0)
return 0;
if (!sgn(delta))
return sol.push_back(L.point(-f / (2 * e))), 1;
sol.push_back(L.point((-f - sqrt(delta)) / (2 * e)));
sol.push_back(L.point((-f + sqrt(delta)) / (2 * e)));
return 2;
}
int getCircleIntersection(const Circle &C1, const Circle &C2, vector<Coord> &sol)
{
lf d = abs(C1.c - C2.c);
if (!sgn(d))
return sgn(C1.r - C2.r) ? 0 : -1; //重合返回-1
if (sgn(C1.r + C2.r - d) < 0 || sgn(fabs(C1.r - C2.r) - d) > 0) //外离或内含
return 0;
lf a = arg(C2.c - C1.c),
da = acos((C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d));
Coord p1 = C1.point(a - da), p2 = C1.point(a + da);
sol.push_back(p1);
if (p1 == p2)
return 1; //相切
return sol.push_back(p2), 2;
}
Line getTangent(const Coord &C, const Coord &P) { return Line(P, Normal(C - P)); } //圆心C，圆上一点P处切线
int getTangents(const Coord &p, const Circle &C, vector<Coord> &sol) //点到圆的切点，返回个数
{
Coord u = p - C.c;
lf d = abs(u);
if (d < C.r)
return 0; //点在圆内
if (!sgn(d - C.r)) //点在圆上
return sol.push_back(p), 1;
lf base = arg(u), ang = acos(C.r / d);
sol.push_back(C.point(base + ang));
sol.push_back(C.point(base - ang));
return 2;
}
int getTangents(Circle A, Circle &B, vector<Coord> &a, vector<Coord> &b) //公共切线的切点
{
int cnt = 0;
if (A.r < B.r)
swap(A, B), swap(a, b); //有时需标记交换
lf d = abs(A.c - B.c),
rdiff = A.r - B.r,
rsum = A.r + B.r;
if (sgn(d - rdiff) < 0)
return 0; //内含
lf base = arg(B.c - A.c);
if (!sgn(d) && !sgn(rdiff))
return -1; //重合，无穷多条切线
if (!sgn(d - rdiff)) //内切，外公切线
{
a.push_back(A.point(base));
b.push_back(B.point(base));
return 1;
}
//有外公切线的情形
lf ang = acos(rdiff / d);
a.push_back(A.point(base + ang));
b.push_back(B.point(base + ang));
a.push_back(A.point(base - ang));
b.push_back(B.point(base - ang));
cnt += 2;
if (!sgn(d - rsum))
{
a.push_back(A.point(base));
b.push_back(B.point(base + M_PI));
++cnt;
}
else if (sgn(d - rsum) > 0)
{
lf ang_in = acos(rsum / d);
a.push_back(A.point(base + ang_in));
b.push_back(B.point(base + ang_in + M_PI));
a.push_back(A.point(base - ang_in));
b.push_back(B.point(base - ang_in + M_PI));
cnt += 2;
}
return cnt;
}
lf AreaCircleWithTriangle(const Circle &C, Coord A, Coord B) //C和三角形OAB的相交面积，如果三角形顶点不在O上则把圆和三角形同时平移，直到有一个顶点在O上
{
int sg = sgn(Cross(A, B));
if (!sg || A == C.c || B == C.c)
return 0;
lf OA = abs(A - C.c), OB = abs(B - C.c), angle = Angle(A, B),
d = DistanceToLine(Coord(), A, B);
if (sgn(OA - C.r) <= 0 && sgn(OB - C.r) <= 0)
return Cross(A, B) / 2;
if (sgn(OA - C.r) >= 0 && sgn(OB - C.r) >= 0 && sgn(d - C.r) >= 0)
return sg * C.r * C.r * angle / 2;
if (sgn(OA - C.r) >= 0 && sgn(OB - C.r) >= 0 && sgn(d - C.r) < 0)
{
Coord prj = getLineProjection(Coord(), A, B);
if (!onSegment(prj, A, B))
return sg * C.r * C.r * angle / 2;
vector<Coord> p;
Line L = Line(A, B - A);
getLineCircleIntersection(L, C, p);
lf s1 = C.r * C.r * angle / 2,
s2 = C.r * C.r * Angle(p[0], p[1]) / 2;
s2 -= fabs(Cross(p[0], p[1]) / 2);
s1 = s1 - s2;
return sg * s1;
}
if (sgn(OB - C.r) < 0)
swap(A, B);
Line L = Line(A, B - A);
vector<Coord> inter;
getLineCircleIntersection(L, C, inter);
Coord inter_point = inter[!onSegment(inter[0], A, B)];
lf s = fabs(Cross(inter_point, A) / 2);
s += C.r * C.r * Angle(inter_point, B) / 2;
return s * sg;
}
lf AreaCircleWithPolygon(const Circle &C, const vector<Coord> &p)
{
lf ans = 0;
for (int i = 0; i < p.size(); ++i)
ans += AreaCircleWithTriangle(C, p[i], p[(i + 1) % p.size()]);
return fabs(ans);
}
Coord getGravityCenter(const vector<Coord> &p) //多边形重心
{
Coord a(0, 0);
lf am = 0, mj;
for (int i = 0; i < p.size(); ++i)
{
mj = Cross(p[i], p[(i + 1) % p.size()]);
a += mj * (p[i] + p[(i + 1) % p.size()]);
am += mj;
}
return a / am / 3.0;
}
三维
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
typedef double lf;
const lf EPS = 1e-9, INF = 1e9;
int sgn(lf d) { return (d > EPS) - (d < -EPS); }
struct Coord3
{
lf X, Y, Z;
friend bool operator!=(const Coord3 &a, const Coord3 &b) { return sgn(a.X - b.X) || sgn(a.Y - b.Y) || sgn(a.Z - b.Z); }
friend bool operator==(const Coord3 &a, const Coord3 &b) { return !(a != b); }
Coord3 &operator+=(const Coord3 &b) { return X += b.X, Y += b.Y, Z += b.Z, *this; }
friend Coord3 operator+(Coord3 a, const Coord3 &b) { return a += b; }
Coord3 &operator-=(const Coord3 &b) { return X -= b.X, Y -= b.Y, Z -= b.Z, *this; }
friend Coord3 operator-(Coord3 a, const Coord3 &b) { return a -= b; }
Coord3 &operator*=(lf d) { return X *= d, Y *= d, Z *= d, *this; }
friend Coord3 operator*(Coord3 a, lf d) { return a *= d; }
friend Coord3 operator*(lf d, Coord3 a) { return a *= d; }
Coord3 &operator/=(lf d) { return X /= d, Y /= d, Z /= d, *this; }
friend Coord3 operator/(Coord3 a, lf d) { return a /= d; }
friend lf Dot(const Coord3 &A, const Coord3 &B) { return A.X * B.X + A.Y * B.Y + A.Z * B.Z; }
friend Coord3 Cross(const Coord3 &A, const Coord3 &B) { return {A.Y * B.Z - A.Z * B.Y, A.Z * B.X - A.X * B.Z, A.X * B.Y - A.Y * B.X}; }
friend lf norm(const Coord3 &A) { return Dot(A, A); }
friend lf abs(const Coord3 &A) { return sqrt(norm(A)); }
friend lf Angle(const Coord3 &A, const Coord3 &B) { return acos(Dot(A, B) / abs(A) / abs(B)); }
friend lf Area2(Coord3 A, Coord3 B, Coord3 C) { return abs(Cross(B - A, C - A)); }
friend lf Volume6(Coord3 A, Coord3 B, Coord3 C, Coord3 D) { return Dot(D - A, Cross(B - A, C - A)); } //四面体体积
friend Coord3 Centroid(Coord3 A, Coord3 B, Coord3 C, Coord3 D) { return (A + B + C + D) / 4.0; } //四面体的重心
friend lf DistanceToPlane(Coord3 p, Coord3 p0, const Coord3 &n) { return Dot(p - p0, n) / abs(n); } //点p到平面p0-n的有向距离
friend Coord3 getPlaneProjection(Coord3 p, Coord3 p0, const Coord3 &n) { return p - n * Dot(p - p0, n); } //点p在平面p0-n上的投影。n必须为单位向量
friend Coord3 LinePlaneIntersection(Coord3 p1, Coord3 p2, Coord3 p0, Coord3 n) //直线p1-p2 与平面p0-n的交点，假设交点唯一存在
{
Coord3 v = p2 - p1;
lf t = Dot(n, p0 - p1) / Dot(n, p2 - p1); //分母为0，直线与平面平行或在平面上
return p1 + v * t; //如果是线段 判断t是否在0~1之间
}
friend lf DistanceToLine(Coord3 P, Coord3 A, Coord3 B) //点P到直线AB的距离
{
Coord3 v1 = B - A, v2 = P - A;
return abs(Cross(v1, v2)) / abs(v1);
}
friend lf DistanceToSeg(Coord3 P, Coord3 A, Coord3 B) //点到线段的距离
{
if (A == B)
return abs(P - A);
Coord3 v1 = B - A, v2 = P - A, v3 = P - B;
if (sgn(Dot(v1, v2)) < 0)
return abs(v2);
if (sgn(Dot(v1, v3)) > 0)
return abs(v3);
return fabs(DistanceToLine(P, A, B));
}
friend bool LineDistance3D(Coord3 p1, Coord3 u, Coord3 p2, Coord3 v, lf &s) //求异面直线 p1+s*u与p2+t*v的公垂线对应的s，如果平行|重合，返回0
{
lf b = Dot(u, u) * Dot(v, v) - Dot(u, v) * Dot(u, v);
if (!sgn(b))
return 0;
lf a = Dot(u, v) * Dot(v, p1 - p2) - Dot(v, v) * Dot(u, p1 - p2);
return s = a / b, 1;
}
friend bool SameSide(Coord3 p1, Coord3 p2, Coord3 a, Coord3 b) { return sgn(Dot(Cross(b - a, p1 - a), Cross(b - a, p2 - a))) >= 0; } //p1和p2是否在线段a-b的同侧
friend bool PointInTri(Coord3 PP, Coord3 P[3]) //点P在三角形P0,P1,p中
{
return SameSide(PP, P[0], P[1], P[2]) &&
SameSide(PP, P[1], P[0], P[2]) &&
SameSide(PP, P[2], P[0], P[1]);
}
friend bool TriSegIntersection(Coord3 P[3], Coord3 A, Coord3 B, Coord3 &PP) //三角形P0P1p是否和线段AB相交，如有则为PP
{
Coord3 n = Cross(P[1] - P[0], P[2] - P[0]);
if (sgn(Dot(n, B - A)) == 0)
return false; //线段A-B和平面P0P1p平行或共面
lf t = Dot(n, P[0] - A) / Dot(n, B - A); //平面A和直线P1-p有惟一交点
if (sgn(t) < 0 || sgn(t - 1) > 0)
return false; //不在线段AB上
return PointInTri(PP = A + (B - A) * t, P);
}
friend bool TriTriIntersection(Coord3 T1[3], Coord3 T2[3]) //空间两三角形是否相交
{
Coord3 P;
for (int i = 0; i < 3; ++i)
if (TriSegIntersection(T1, T2[i], T2[(i + 1) % 3], P) ||
TriSegIntersection(T2, T1[i], T1[(i + 1) % 3], P))
return 1;
return 0;
}
friend lf SegSegDistance(Coord3 a1, Coord3 b1, Coord3 a2, Coord3 b2, Coord3 &ans1, Coord3 &ans2) //空间两直线上最近点对 返回最近距离 点对保存在ans1 ans2中
{
Coord3 v1 = (a1 - b1), v2 = (a2 - b2);
Coord3 N = Cross(v1, v2);
Coord3 ab = (a1 - a2);
lf ans = Dot(N, ab) / abs(N);
Coord3 d1 = b1 - a1, d2 = b2 - a2, cd = Cross(d1, d2);
lf nd = norm(cd), t1 = Dot(Cross(a2 - a1, d2), cd) / nd, t2 = Dot(Cross(a2 - a1, d1), cd) / nd;
return ans1 = a1 + (b1 - a1) * t1, ans2 = a2 + (b2 - a2) * t2, fabs(ans);
}
friend bool InsideWithMinDistance(Coord3 PP, Coord3 *P, lf dist) //判断PP是否在三角形P中，并且到三条边的距离都至少为dist。保证P,A,B,C共面
{
return PointInTri(PP, P) &&
DistanceToLine(PP, P[0], P[1]) >= dist ||
DistanceToLine(PP, P[1], P[2]) >= dist ||
DistanceToLine(PP, P[2], P[0]) >= dist;
}
friend lf minBall(const vector<Coord3> &data, const lf eps = EPS * 1e-3) //模拟退火求最小球覆盖，EPS玄学调整；返回半径
{
lf step = 1, ans = INF;
for (Coord3 z{0, 0, 0}; step > eps; step *= 0.99)
{
int s = 0;
for (int i = 0; i < data.size(); ++i)
if (abs(z - data[s]) < abs(z - data[i]))
s = i;
ans = min(ans, abs(z - data[s]));
z -= (z - data[s]) * step;
}
return ans;
}
};
struct Sphere
{
Coord3 o;
lf r;
Coord3 point(lf lat, lf lng) const //经纬度确定球面上一点
{
lat *= M_PI / 180;
lng *= M_PI / 180;
return o + r * Coord3{cos(lat) * cos(lng), cos(lat) * sin(lng), sin(lat)};
}
};
struct ConvexPolyhedron //空间多边形和凸包问题
{
struct Face
{
int v[3];
Face(int a, int b, int c) { v[0] = a, v[1] = b, v[2] = c; }
Coord3 Normal(const vector<Coord3> &P) const { return Cross(P[v[1]] - P[v[0]], P[v[2]] - P[v[0]]); }
bool CanSee(const vector<Coord3> &P, int i) const { return Dot(P[i] - P[v[0]], Normal(P)) > 0; } //f是否能看见P[i]
};
vector<Face> faces;
vector<Coord3> p;
ConvexPolyhedron(vector<Coord3> P) : p(P)
{
for (int i = 0; i < p.size(); ++i)
P[i] += Coord3{randEPS(), randEPS(), randEPS()};
vector<vector<int>> vis(P.size(), vector<int>(P.size()));
faces.push_back(Face(0, 1, 2)); //由于已经进行扰动，前三个点不共线
faces.push_back(Face(2, 1, 0));
for (int i = 3; i < P.size(); ++i)
{
vector<Face> next;
for (int j = 0; j < faces.size(); ++j) //计算每条边的「左面」的可见性
{
Face &f = faces[j];
int res = f.CanSee(P, i);
if (!res)
next.push_back(f);
for (int k = 0; k < 3; ++k)
vis[f.v[k]][f.v[(k + 1) % 3]] = res;
}
for (int j = 0; j < faces.size(); ++j)
for (int k = 0; k < 3; ++k)
{
int a = faces[j].v[k], b = faces[j].v[(k + 1) % 3];
if (vis[a][b] != vis[b][a] && vis[a][b]) //(a,b)是分界线，左边对P[i]可见
next.push_back(Face(a, b, i));
}
swap(faces, next);
}
}
lf randEPS() { return (rand() / lf(RAND_MAX) - 0.5) * EPS; }
Coord3 centroid() //三维凸包重心
{
Coord3 C = p[0], tot{0, 0, 0};
lf totv = 0;
for (int i = 0; i < faces.size(); ++i)
{
Coord3 p1 = p[faces[i].v[0]], p2 = p[faces[i].v[1]], p3 = p[faces[i].v[2]];
lf v = -Volume6(p1, p2, p3, C);
totv += v;
tot += v * Centroid(p1, p2, p3, C);
}
return tot / totv;
}
lf dist(Coord3 C) //凸包内一点到表面最近距离
{
lf ans = INF;
for (int i = 0; i < faces.size(); ++i)
{
Coord3 p1 = p[faces[i].v[0]], p2 = p[faces[i].v[1]], p3 = p[faces[i].v[2]];
ans = min(ans, fabs(-Volume6(p1, p2, p3, C) / Area2(p1, p2, p3)));
}
return ans;
}
};