二维计算几何模板整理

分类： ACMICPC 2013-10-06 18:55 130人阅读 评论(0) 收藏 举报

[cpp] view plain copy print ?

1. /***********三分法求函数极值*************/
2. void solve()
3. {
4. double L, R, m, mm, mv, mmv;
5. while (L + eps < R)
6. {
7. m = (L + R) / 2;
8. mm = (m + R) / 2;
9. mv = calc(m);
10. mmv = calc(mm);
11. if (mv <= mmv) R = mm; //三分法求最大值时改为mv>=mmv
12. else L = m;
13. }
14. }
15. /*************基础***********/
16. int dcmp(double x) {
17. if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;
18. }
19. struct Point {
20. double x, y;
21. Point(double x=0, double y=0):x(x),y(y) { }
22. };
24. typedef Point Vector;
26. Vector operator + (Vector A, Vector B) { return Vector(A.x+B.x, A.y+B.y); }
27. Vector operator - (Point A, Point B) { return Vector(A.x-B.x, A.y-B.y); }
28. Vector operator * (Vector A, double p) { return Vector(A.x*p, A.y*p); }
29. Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }
31. bool operator < (const Point& a, const Point& b) {
32. return a.x < b.x || (a.x == b.x && a.y < b.y);
33. }
35. bool operator == (const Point& a, const Point &b) {
36. return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;
37. }
39. double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; }
40. double Length(Vector A) { return sqrt(Dot(A, A)); }
41. double Angle(Vector A, Vector B) { return acos(Dot(A, B) / Length(A) / Length(B)); }
42. double angle(Vector v) { return atan2(v.y, v.x); }
43. double Cross(Vector A, Vector B) { return A.x*B.y - A.y*B.x; }
44. Vector vecunit(Vector x){ return x / Length(x);} //单位向量
45. Vector Normal(Vector x) { return Point(-x.y, x.x) / Length(x);} //垂直法向量
46. Vector Rotate(Vector A, double rad) {
47. return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
48. }
49. double Area2(const Point A, const Point B, const Point C) { return Length(Cross(B-A, C-A)); }
51. /****************直线与线段**************/
53. //求直线p+tv和q+tw的交点 Cross(v, w) == 0无交点
54. Point GetLineIntersection(Point p, Vector v, Point q, Vector w)
55. {
56. Vector u = p-q;
57. double t = Cross(w, u) / Cross(v, w);
58. return p + v*t;
59. }
61. //点p在直线ab的投影
62. Point GetLineProjection(Point P, Point A, Point B) {
63. Vector v = B-A;
64. return A+v*(Dot(v, P-A) / Dot(v, v));
65. }
67. //点到直线距离
68. double DistanceToLine(Point P, Point A, Point B) {
69. Vector v1 = B - A, v2 = P - A;
70. return fabs(Cross(v1, v2)) / Length(v1); // 如果不取绝对值，得到的是有向距离
71. }
72. //点在p线段上
73. bool OnSegment(Point p, Point a1, Point a2) {
74. return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p)) < 0;
75. }
76. // 过两点p1, p2的直线一般方程ax+by+c=0
77. // (x2-x1)(y-y1) = (y2-y1)(x-x1)
78. void getLineGeneralEquation(const Point& p1, const Point& p2, double& a, double& b, double &c) {
79. a = p2.y-p1.y;
80. b = p1.x-p2.x;
81. c = -a*p1.x - b*p1.y;
82. }
83. //点到线段距离
84. double DistanceToSegment(Point p, Point a, Point b)
85. {
86. if(a == b) return Length(p-a);
87. Vector v1 = b-a, v2 = p-a, v3 = p-b;
88. if(dcmp(Dot(v1, v2)) < 0) return Length(v2);
89. else if(dcmp(Dot(v1, v3)) > 0) return Length(v3);
90. else return fabs(Cross(v1, v2)) / Length(v1);
91. }
92. //两线段最近距离
93. double dis_pair_seg(Point p1, Point p2, Point p3, Point p4)
94. {
95. return min(min(DistanceToSegment(p1, p3, p4), DistanceToSegment(p2, p3, p4)),
96. min(DistanceToSegment(p3, p1, p2), DistanceToSegment(p4, p1, p2)));
97. }
98. //线段相交判定
99. bool SegmentItersection(Point a1, Point a2, Point b1, Point b2)
100. {
101. double c1 = Cross(a2-a1, b1-a1), c2 = Cross(a2-a1, b2-a1),
102. c3 = Cross(b2-b1, a1-b1), c4 = Cross(b2-b1, a2-b1);
103. return dcmp(c1)*dcmp(c2) < 0 && dcmp(c3)*dcmp(c4) < 0;
104. }
105. // 有向直线。它的左边就是对应的半平面
106. struct Line {
107. Point P; // 直线上任意一点
108. Vector v; // 方向向量
109. double ang; // 极角，即从x正半轴旋转到向量v所需要的角（弧度）
110. Line() {}
111. Line(Point P, Vector v):P(P),v(v){ ang = atan2(v.y, v.x); }
112. bool operator < (const Line& L) const {
113. return ang < L.ang;
114. }
115. };
116. //两直线交点
117. Point GetLineIntersection(Line a, Line b) {
118. return GetLineIntersection(a.p, a.v, b.p, b.v);
119. }
121. // 点p在有向直线L的左边（线上不算）
122. bool OnLeft(const Line& L, const Point& p) {
123. return Cross(L.v, p-L.P) > 0;
124. }
126. // 二直线交点，假定交点惟一存在
127. Point GetLineIntersection(const Line& a, const Line& b) {
128. Vector u = a.P-b.P;
129. double t = Cross(b.v, u) / Cross(a.v, b.v);
130. return a.P+a.v*t;
131. }
133. // 半平面交主过程
134. vector<Point> HalfplaneIntersection(vector<Line> L) {
135. int n = L.size();
136. sort(L.begin(), L.end()); // 按极角排序
138. int first, last; // 双端队列的第一个元素和最后一个元素的下标
139. vector<Point> p(n); // p[i]为q[i]和q[i+1]的交点
140. vector<Line> q(n); // 双端队列
141. vector<Point> ans; // 结果
143. q[first=last=0] = L[0]; // 双端队列初始化为只有一个半平面L[0]
144. for(int i = 1; i < n; i++) {
145. while(first < last && !OnLeft(L[i], p[last-1])) last--;
146. while(first < last && !OnLeft(L[i], p[first])) first++;
147. q[++last] = L[i];
148. if(fabs(Cross(q[last].v, q[last-1].v)) < eps) { // 两向量平行且同向，取内侧的一个
149. last--;
150. if(OnLeft(q[last], L[i].P)) q[last] = L[i];
151. }
152. if(first < last) p[last-1] = GetLineIntersection(q[last-1], q[last]);
153. }
154. while(first < last && !OnLeft(q[first], p[last-1])) last--; // 删除无用平面
155. if(last - first <= 1) return ans; // 空集
156. p[last] = GetLineIntersection(q[last], q[first]); // 计算首尾两个半平面的交点
158. // 从deque复制到输出中
159. for(int i = first; i <= last; i++) ans.push_back(p[i]);
160. return ans;
161. }
163. /***********多边形**************/
164. //多边形有向面积
165. double PolygonArea(vector<Point> p) {
166. int n = p.size();
167. double area = 0;
168. for(int i = 1; i < n-1; i++)
169. area += Cross(p[i]-p[0], p[i+1]-p[0]);
170. return area/2;
171. }
172. // 点集凸包
173. // 如果不希望在凸包的边上有输入点，把两个 <= 改成 <
174. // 注意：输入点集会被修改
175. vector<Point> ConvexHull(vector<Point>& p) {
176. // 预处理，删除重复点
177. sort(p.begin(), p.end());
178. p.erase(unique(p.begin(), p.end()), p.end());
180. int n = p.size();
181. int m = 0;
182. vector<Point> ch(n+1);
183. for(int i = 0; i < n; i++) {
184. while(m > 1 && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) <= 0) m--;
185. ch[m++] = p[i];
186. }
187. int k = m;
188. for(int i = n-2; i >= 0; i--) {
189. while(m > k && Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) <= 0) m--;
190. ch[m++] = p[i];
191. }
192. if(n > 1) m--;
193. ch.resize(m);
194. return ch;
195. }
197. // 凸包直径返回 点集直径的平方
198. int diameter2(vector<Point>& points) {
199. vector<Point> p = ConvexHull(points);
200. int n = p.size();
201. if(n == 1) return 0;
202. if(n == 2) return Dist2(p[0], p[1]);
203. p.push_back(p[0]); // 免得取模
204. int ans = 0;
205. for(int u = 0, v = 1; u < n; u++) {
206. // 一条直线贴住边p[u]-p[u+1]
207. for(;;) {
208. int diff = Cross(p[u+1]-p[u], p[v+1]-p[v]);
209. if(diff <= 0) {
210. ans = max(ans, Dist2(p[u], p[v])); // u和v是对踵点
211. if(diff == 0) ans = max(ans, Dist2(p[u], p[v+1])); // diff == 0时u和v+1也是对踵点
212. break;
213. }
214. v = (v + 1) % n;
215. }
216. }
217. return ans;
218. }
219. //两凸包最近距离
220. double RC_Distance(Point *ch1, Point *ch2, int n, int m)
221. {
222. int q=0, p=0;
223. REP(i, n) if(ch1[i].y-ch1[p].y < -eps) p=i;
224. REP(i, m) if(ch2[i].y-ch2[q].y > eps) q=i;
225. ch1[n]=ch1[0]; ch2[m]=ch2[0];
227. double tmp, ans=1e100;
228. REP(i, n)
229. {
230. while((tmp = Cross(ch1[p+1]-ch1[p], ch2[q+1]-ch1[p]) - Cross(ch1[p+1]-ch1[p], ch2[q]- ch1[p])) > eps)
231. q=(q+1)%m;
232. if(tmp < -eps) ans = min(ans,DistanceToSegment(ch2[q],ch1[p],ch1[p+1]));
233. else ans = min(ans,dis_pair_seg(ch1[p],ch1[p+1],ch2[q],ch2[q+1]));
234. p=(p+1)%n;
235. }
236. return ans;
237. }
238. //凸包最大内接三角形
239. double RC_Triangle(Point* res,int n)// 凸包最大内接三角形
240. {
241. if(n<3) return 0;
242. double ans=0, tmp;
243. res[n] = res[0];
244. int j, k;
245. REP(i, n)
246. {
247. j = (i+1)%n;
248. k = (j+1)%n;
249. while((j != k) && (k != i))
250. {
251. while(Cross(res[j] - res[i], res[k+1] - res[i]) > Cross(res[j] - res[i], res[k] - res[i])) k= (k+1)%n;
252. tmp = Cross(res[j] - res[i], res[k] - res[i]);if(tmp > ans) ans = tmp;
253. j = (j+1)%n;
254. }
255. }
256. return ans;
257. }
258. //模拟退火求费马点 保存在ptres中
259. double fermat_point(Point *pt, int n, Point& ptres)
260. {
261. Point u, v;
262. double step = 0.0, curlen, explen, minlen;
263. int i, j, k, idx;
264. bool flag;
265. u.x = u.y = v.x = v.y = 0.0;
266. REP(i, n)
267. {
268. step += fabs(pt[i].x) + fabs(pt[i].y);
269. u.x += pt[i].x;
270. u.y += pt[i].y;
271. }
272. u.x /= n;
273. u.y /= n;
274. flag = 0;
275. while(step > eps)
276. {
277. for(k = 0; k < 10; step /= 2, ++k)
278. for(i = -1; i <= 1; ++i)
279. for(j = -1; j <= 1; ++j)
280. {
281. v.x = u.x + step*i;
282. v.y = u.y + step*j;
283. curlen = explen = 0.0;
284. REP(idx, n)
285. {
286. curlen += dist(u, pt[idx]);
287. explen += dist(v, pt[idx]);
288. }
289. if(curlen > explen)
290. {
291. u = v;
292. minlen = explen;
293. flag = 1;
294. }
295. }
296. }
297. ptres = u;
298. return flag ? minlen : curlen;
299. }
300. //最近点对
301. bool cmpxy(const Point& a, const Point& b)
302. {
303. if(a.x != b.x)
304. return a.x < b.x;
305. return a.y < b.y;
306. }
307. bool cmpy(const int& a, const int& b)
308. {
309. return point[a].y < point[b].y;
310. }
311. double Closest_Pair(int left, int right)
312. {
313. double d = INF;
314. if(left==right)
315. return d;
316. if(left + 1 == right)
317. return dis(left, right);
318. int mid = (left+right)>>1;
319. double d1 = Closest_Pair(left,mid);
320. double d2 = Closest_Pair(mid+1,right);
321. d = min(d1,d2);
322. int i,j,k=0;
323. //分离出宽度为d的区间
324. for(i = left; i <= right; i++)
325. {
326. if(fabs(point[mid].x-point[i].x) <= d)
327. tmpt[k++] = i;
328. }
329. sort(tmpt,tmpt+k,cmpy);
330. //线性扫描
331. for(i = 0; i < k; i++)
332. {
333. for(j = i+1; j < k && point[tmpt[j]].y-point[tmpt[i]].y<d; j++)
334. {
335. double d3 = dis(tmpt[i],tmpt[j]);
336. if(d > d3)
337. d = d3;
338. }
339. }
340. return d;
341. }
343. /************圆************/
344. struct Circle
345. {
346. Point c;
347. double r;
348. Circle(){}
349. Circle(Point c, double r):c(c), r(r){}
350. Point point(double a) //根据圆心角求点坐标
351. {
352. return Point(c.x+cos(a)*r, c.y+sin(a)*r);
353. }
354. };
355. //直线与圆交点 返回个数
356. int getLineCircleIntersection(Line L, Circle C, double& t1, double& t2, vector<Point>& sol){
357. double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;
358. double e = a*a + c*c, f = 2*(a*b + c*d), g = b*b + d*d - C.r*C.r;
359. double delta = f*f - 4*e*g; // 判别式
360. if(dcmp(delta) < 0) return 0; // 相离
361. if(dcmp(delta) == 0) { // 相切
362. t1 = t2 = -f / (2 * e); sol.push_back(L.point(t1));
363. return 1;
364. }
365. // 相交
366. t1 = (-f - sqrt(delta)) / (2 * e); sol.push_back(L.point(t1));
367. t2 = (-f + sqrt(delta)) / (2 * e); sol.push_back(L.point(t2));
368. return 2;
369. }
370. //两圆交点 返回个数
371. int getCircleCircleIntersection(Circle C1, Circle C2, vector<Point>& sol) {
372. double d = Length(C1.c - C2.c);
373. if(dcmp(d) == 0) {
374. if(dcmp(C1.r - C2.r) == 0) return -1; // 重合，无穷多交点
375. return 0;
376. }
377. if(dcmp(C1.r + C2.r - d) < 0) return 0;
378. if(dcmp(fabs(C1.r-C2.r) - d) > 0) return 0;
380. double a = angle(C2.c - C1.c);
381. double da = acos((C1.r*C1.r + d*d - C2.r*C2.r) / (2*C1.r*d));
382. Point p1 = C1.point(a-da), p2 = C1.point(a+da);
384. sol.push_back(p1);
385. if(p1 == p2) return 1;
386. sol.push_back(p2);
387. return 2;
388. }
390. //三角形外接圆
391. Circle CircumscribedCircle(Point p1, Point p2, Point p3) {
392. double Bx = p2.x-p1.x, By = p2.y-p1.y;
393. double Cx = p3.x-p1.x, Cy = p3.y-p1.y;
394. double D = 2*(Bx*Cy-By*Cx);
395. double cx = (Cy*(Bx*Bx+By*By) - By*(Cx*Cx+Cy*Cy))/D + p1.x;
396. double cy = (Bx*(Cx*Cx+Cy*Cy) - Cx*(Bx*Bx+By*By))/D + p1.y;
397. Point p = Point(cx, cy);
398. return Circle(p, Length(p1-p));
399. }
401. //三角形内切圆
402. Circle InscribedCircle(Point p1, Point p2, Point p3) {
403. double a = Length(p2-p3);
404. double b = Length(p3-p1);
405. double c = Length(p1-p2);
406. Point p = (p1*a+p2*b+p3*c)/(a+b+c);
407. return Circle(p, DistanceToLine(p, p1, p2));
408. }
410. // 过点p到圆C的切线。v[i]是第i条切线的向量。返回切线条数
411. int getTangents(Point p, Circle C, Vector* v) {
412. Vector u = C.c - p;
413. double dist = Length(u);
414. if(dist < C.r) return 0;
415. else if(dcmp(dist - C.r) == 0) { // p在圆上，只有一条切线
416. v[0] = Rotate(u, PI/2);
417. return 1;
418. } else {
419. double ang = asin(C.r / dist);
420. v[0] = Rotate(u, -ang);
421. v[1] = Rotate(u, +ang);
422. return 2;
423. }
424. }
426. //所有经过点p 半径为r 且与直线L相切的圆心
427. vector<Point> CircleThroughPointTangentToLineGivenRadius(Point p, Line L, double r) {
428. vector<Point> ans;
429. double t1, t2;
430. getLineCircleIntersection(L.move(-r), Circle(p, r), t1, t2, ans);
431. getLineCircleIntersection(L.move(r), Circle(p, r), t1, t2, ans);
432. return ans;
433. }
435. //半径为r 与a b两直线相切的圆心
436. vector<Point> CircleTangentToLinesGivenRadius(Line a, Line b, double r) {
437. vector<Point> ans;
438. Line L1 = a.move(-r), L2 = a.move(r);
439. Line L3 = b.move(-r), L4 = b.move(r);
440. ans.push_back(GetLineIntersection(L1, L3));
441. ans.push_back(GetLineIntersection(L1, L4));
442. ans.push_back(GetLineIntersection(L2, L3));
443. ans.push_back(GetLineIntersection(L2, L4));
444. return ans;
445. }
447. //与两圆相切 半径为r的所有圆心
448. vector<Point> CircleTangentToTwoDisjointCirclesWithRadius(Circle c1, Circle c2, double r) {
449. vector<Point> ans;
450. Vector v = c2.c - c1.c;
451. double dist = Length(v);
452. int d = dcmp(dist - c1.r -c2.r - r*2);
453. if(d > 0) return ans;
454. getCircleCircleIntersection(Circle(c1.c, c1.r+r), Circle(c2.c, c2.r+r), ans);
455. return ans;
456. }
458. //多边形与圆相交面积
459. Point GetIntersection(Line a, Line b) //线段交点
460. {
461. Vector u = a.p-b.p;
462. double t = Cross(b.v, u) / Cross(a.v, b.v);
463. return a.p + a.v*t;
464. }
465. bool OnSegment(Point p, Point a1, Point a2)
466. {
467. return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p)) < 0;
468. }
469. bool InCircle(Point x, Circle c) { return dcmp(sqr(c.r) - sqr(Length(c.c - x))) >= 0;}
470. bool OnCircle(Point x, Circle c) { return dcmp(sqr(c.r) - sqr(Length(c.c - x))) == 0;}
471. //线段与圆的交点
472. int getSegCircleIntersection(Line L, Circle C, Point* sol)
473. {
474. Vector nor = normal(L.v);
475. Line pl = Line(C.c, nor);
476. Point ip = GetIntersection(pl, L);
477. double dis = Length(ip - C.c);
478. if (dcmp(dis - C.r) > 0) return 0;
479. Point dxy = vecunit(L.v) * sqrt(sqr(C.r) - sqr(dis));
480. int ret = 0;
481. sol[ret] = ip + dxy;
482. if (OnSegment(sol[ret], L.p, L.point(1))) ret++;
483. sol[ret] = ip - dxy;
484. if (OnSegment(sol[ret], L.p, L.point(1))) ret++;
485. return ret;
486. }
487. double SegCircleArea(Circle C, Point a, Point b) //线段切割圆
488. {
489. double a1 = angle(a - C.c);
490. double a2 = angle(b - C.c);
491. double da = fabs(a1 - a2);
492. if (da > PI) da = PI * 2.0 - da;
493. return dcmp(Cross(b - C.c, a - C.c)) * da * sqr(C.r) / 2.0;
494. }
496. double PolyCiclrArea(Circle C, Point *p, int n)//多边形与圆相交面积
497. {
498. double ret = 0.0;
499. Point sol[2];
500. p[n] = p[0];
501. REP(i, n)
502. {
503. double t1, t2;
504. int cnt = getSegCircleIntersection(Line(p[i], p[i+1]-p[i]), C, sol);
505. if (cnt == 0)
506. {
507. if (!InCircle(p[i], C) || !InCircle(p[i+1], C)) ret += SegCircleArea(C, p[i], p[i+1]);
508. else ret += Cross(p[i+1] - C.c, p[i] - C.c) / 2.0;
509. }
510. if (cnt == 1)
511. {
512. if (InCircle(p[i], C) && !InCircle(p[i+1], C)) ret += Cross(sol[0] - C.c, p[i] - C.c) / 2.0, ret += SegCircleArea(C, sol[0], p[i+1]);
513. else ret += SegCircleArea(C, p[i], sol[0]), ret += Cross(p[i+1] - C.c, sol[0] - C.c) / 2.0;
514. }
515. if (cnt == 2)
516. {
517. if ((p[i] < p[i + 1]) ^ (sol[0] < sol[1])) swap(sol[0], sol[1]);
518. ret += SegCircleArea(C, p[i], sol[0]);
519. ret += Cross(sol[1] - C.c, sol[0] - C.c) / 2.0;
520. ret += SegCircleArea(C, sol[1], p[i+1]);
521. }
522. }
523. return fabs(ret);
524. }