LA 2572 Viva Confetti (Geometry.Circle)
2013-05-01 20:56
399 查看
https://icpcarchive.ecs.baylor.edu/index.php?option=com_onlinejudge&Itemid=8&page=problem_stats&problemid=573&category=
一道很好的几何题。要求是给出一些圆,按顺序覆盖上去,问哪些圆是可以被看见的。
看刘汝佳的书,开始的时候不明白“每个可见部分都是由一些‘小圆弧’围城的”这样又怎么样,直到我看了他的代码以后才懂。其实意思就是,因为每个可见部分都是由小圆弧围成,所以,如果我们将小圆弧中点(不会是圆与圆的交点)相对于弧的位置向里或向外稍微移动一下,然后找到有哪个圆最后覆盖在在那上面,那么这个圆必定能被看见。
刚接触几何,表示还真想不到可以这样做,所以就觉得这个做法实在太妙了!
代码如下:
View Code
——written by Lyon
一道很好的几何题。要求是给出一些圆,按顺序覆盖上去,问哪些圆是可以被看见的。
看刘汝佳的书,开始的时候不明白“每个可见部分都是由一些‘小圆弧’围城的”这样又怎么样,直到我看了他的代码以后才懂。其实意思就是,因为每个可见部分都是由小圆弧围成,所以,如果我们将小圆弧中点(不会是圆与圆的交点)相对于弧的位置向里或向外稍微移动一下,然后找到有哪个圆最后覆盖在在那上面,那么这个圆必定能被看见。
刚接触几何,表示还真想不到可以这样做,所以就觉得这个做法实在太妙了!
代码如下:
View Code
#include <cstdio> #include <cstring> #include <cmath> #include <set> #include <vector> #include <iostream> #include <algorithm> using namespace std; struct Point { double x, y; Point() {} Point(double x, double y) : x(x), y(y) {} } ; template<class T> T sqr(T x) { return x * x;} inline double ptDis(Point a, Point b) { return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));} // basic calculations typedef Point Vec; Vec operator + (Vec a, Vec b) { return Vec(a.x + b.x, a.y + b.y);} Vec operator - (Vec a, Vec b) { return Vec(a.x - b.x, a.y - b.y);} Vec operator * (Vec a, double p) { return Vec(a.x * p, a.y * p);} Vec operator / (Vec a, double p) { return Vec(a.x / p, a.y / p);} const double EPS = 5e-13; const double PI = acos(-1.0); inline int sgn(double x) { return fabs(x) < EPS ? 0 : (x < 0 ? -1 : 1);} bool operator < (Point a, Point b) { return a.x < b.x || (a.x == b.x && a.y < b.y);} bool operator == (Point a, Point b) { return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0;} inline double dotDet(Vec a, Vec b) { return a.x * b.x + a.y * b.y;} inline double vecLen(Vec x) { return sqrt(sqr(x.x) + sqr(x.y));} inline double angle(Vec v) { return atan2(v.y, v.x);} inline double angle(Vec a, Vec b) { return acos(dotDet(a, b) / vecLen(a) / vecLen(b));} inline double crossDet(Vec a, Vec b) { return a.x * b.y - a.y * b.x;} inline double triArea(Point a, Point b, Point c) { return fabs(crossDet(b - a, c - a));} inline Vec vecUnit(Vec x) { return x / vecLen(x);} inline Vec rotate(Vec x, double rad) { return Vec(x.x * cos(rad) - x.y * sin(rad), x.x * sin(rad) + x.y * cos(rad));} Vec normal(Vec x) { double len = vecLen(x); return Vec(- x.y / len, x.x / len); } struct Line { Point s, t; Line() {} Line(Point s, Point t) : s(s), t(t) {} Point point(double x) { return s + (t - s) * x; } Line move(double x) { Vec nor = normal(t - s); nor = nor / vecLen(nor) * x; return Line(s + nor, t + nor); } Vec vec() { return t - s;} } ; typedef Line Seg; inline bool onSeg(Point x, Point a, Point b) { return sgn(crossDet(a - x, b - x)) == 0 && sgn(dotDet(a - x, b - x)) < 0;} inline bool onSeg(Point x, Seg s) { return onSeg(x, s.s, s.t);} // 0 : not intersect // 1 : proper intersect // 2 : improper intersect int segIntersect(Point a, Point c, Point b, Point d) { Vec v1 = b - a, v2 = c - b, v3 = d - c, v4 = a - d; int a_bc = sgn(crossDet(v1, v2)); int b_cd = sgn(crossDet(v2, v3)); int c_da = sgn(crossDet(v3, v4)); int d_ab = sgn(crossDet(v4, v1)); if (a_bc * c_da > 0 && b_cd * d_ab > 0) return 1; if (onSeg(b, a, c) && c_da) return 2; if (onSeg(c, b, d) && d_ab) return 2; if (onSeg(d, c, a) && a_bc) return 2; if (onSeg(a, d, b) && b_cd) return 2; return 0; } inline int segIntersect(Seg a, Seg b) { return segIntersect(a.s, a.t, b.s, b.t);} // point of the intersection of 2 lines Point lineIntersect(Point P, Vec v, Point Q, Vec w) { Vec u = P - Q; double t = crossDet(w, u) / crossDet(v, w); return P + v * t; } inline Point lineIntersect(Line a, Line b) { return lineIntersect(a.s, a.t - a.s, b.s, b.t - b.s);} // Warning: This is a DIRECTED Distance!!! double pt2Line(Point x, Point a, Point b) { Vec v1 = b - a, v2 = x - a; return crossDet(v1, v2) / vecLen(v1); } inline double pt2Line(Point x, Line L) { return pt2Line(x, L.s, L.t);} double pt2Seg(Point x, Point a, Point b) { if (a == b) return vecLen(x - a); Vec v1 = b - a, v2 = x - a, v3 = x - b; if (sgn(dotDet(v1, v2)) < 0) return vecLen(v2); if (sgn(dotDet(v1, v3)) > 0) return vecLen(v3); return fabs(crossDet(v1, v2)) / vecLen(v1); } inline double pt2Seg(Point x, Seg s) { return pt2Seg(x, s.s, s.t);} struct Poly { vector<Point> pt; Poly() {} Poly(vector<Point> pt) : pt(pt) {} double area() { double ret = 0.0; int sz = pt.size(); for (int i = 1; i < sz; i++) { ret += crossDet(pt[i], pt[i - 1]); } return fabs(ret / 2.0); } } ; struct Circle { Point c; double r; Circle() {} Circle(Point c, double r) : c(c), r(r) {} Point point(double a) { return Point(c.x + cos(a) * r, c.y + sin(a) * r); } } ; int lineCircleIntersect(Line L, Circle C, double &t1, double &t2, vector<Point> &sol) { double a = L.s.x, b = L.t.x - C.c.x, c = L.s.y, d = L.t.y - C.c.y; double e = sqr(a) + sqr(c), f = 2 * (a * b + c * d), g = sqr(b) + sqr(d) - sqr(C.r); double delta = sqr(f) - 4.0 * e * g; if (sgn(delta) < 0) return 0; if (sgn(delta) == 0) { t1 = t2 = -f / (2.0 * e); sol.push_back(L.point(t1)); return 1; } t1 = (-f - sqrt(delta)) / (2.0 * e); sol.push_back(L.point(t1)); t2 = (-f + sqrt(delta)) / (2.0 * e); sol.push_back(L.point(t2)); return 2; } int lineCircleIntersect(Line L, Circle C, vector<Point> &sol) { Vec dir = L.t - L.s, nor = normal(dir); Point mid = lineIntersect(C.c, nor, L.s, dir); double len = sqr(C.r) - sqr(ptDis(C.c, mid)); if (sgn(len) < 0) return 0; if (sgn(len) == 0) { sol.push_back(mid); return 1; } Vec dis = vecUnit(dir); len = sqrt(len); sol.push_back(mid + dis * len); sol.push_back(mid - dis * len); return 2; } // -1 : coincide int circleCircleIntersect(Circle C1, Circle C2, vector<Point> &sol) { double d = vecLen(C1.c - C2.c); if (sgn(d) == 0) { if (sgn(C1.r - C2.r) == 0) { return -1; } return 0; } if (sgn(C1.r + C2.r - d) < 0) return 0; if (sgn(fabs(C1.r - C2.r) - d) > 0) return 0; double a = angle(C2.c - C1.c); double da = acos((sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2.0 * C1.r * d)); Point p1 = C1.point(a - da), p2 = C1.point(a + da); sol.push_back(p1); if (p1 == p2) return 1; sol.push_back(p2); return 2; } void circleCircleIntersect(Circle C1, Circle C2, vector<double> &sol) { double d = vecLen(C1.c - C2.c); if (sgn(d) == 0) return ; if (sgn(C1.r + C2.r - d) < 0) return ; if (sgn(fabs(C1.r - C2.r) - d) > 0) return ; double a = angle(C2.c - C1.c); double da = acos((sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2.0 * C1.r * d)); sol.push_back(a - da); sol.push_back(a + da); } int tangent(Point p, Circle C, vector<Vec> &sol) { Vec u = C.c - p; double dist = vecLen(u); if (dist < C.r) return 0; if (sgn(dist - C.r) == 0) { sol.push_back(rotate(u, PI / 2.0)); return 1; } double ang = asin(C.r / dist); sol.push_back(rotate(u, -ang)); sol.push_back(rotate(u, ang)); return 2; } // ptA : points of tangency on circle A // ptB : points of tangency on circle B int tangent(Circle A, Circle B, vector<Point> &ptA, vector<Point> &ptB) { if (A.r < B.r) { swap(A, B); swap(ptA, ptB); } int d2 = sqr(A.c.x - B.c.x) + sqr(A.c.y - B.c.y); int rdiff = A.r - B.r, rsum = A.r + B.r; if (d2 < sqr(rdiff)) return 0; double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x); if (d2 == 0 && A.r == B.r) return -1; if (d2 == sqr(rdiff)) { ptA.push_back(A.point(base)); ptB.push_back(B.point(base)); return 1; } double ang = acos((A.r - B.r) / sqrt(d2)); ptA.push_back(A.point(base + ang)); ptB.push_back(B.point(base + ang)); ptA.push_back(A.point(base - ang)); ptB.push_back(B.point(base - ang)); if (d2 == sqr(rsum)) { ptA.push_back(A.point(base)); ptB.push_back(B.point(PI + base)); } else if (d2 > sqr(rsum)) { ang = acos((A.r + B.r) / sqrt(d2)); ptA.push_back(A.point(base + ang)); ptB.push_back(B.point(PI + base + ang)); ptA.push_back(A.point(base - ang)); ptB.push_back(B.point(PI + base - ang)); } return (int) ptA.size(); } /****************** template above *******************/ #define DEC(i, a, b) for (int i = (a); i >= (b); i--) #define REP(i, n) for (int i = 0; i < (n); i++) #define _clr(x) memset(x, 0, sizeof(x)) #define PB push_back Circle c[111]; int topCircle(Point x, int n) { DEC(i, n - 1, 0) { if (ptDis(c[i].c, x) < c[i].r) return i; } return -1; } int main() { // freopen("in", "r", stdin); int n; while (cin >> n && n) { REP(i, n) cin >> c[i].c.x >> c[i].c.y >> c[i].r; set<int> vis; vis.clear(); REP(i, n) { vector<double> pos; pos.clear(); pos.PB(0.0); pos.PB(2.0 * PI); REP(j, n) { circleCircleIntersect(c[i], c[j], pos); } sort(pos.begin(), pos.end()); int sz = pos.size() - 1; REP(j, sz) { double mid = (pos[j] + pos[j + 1]) / 2.0; for (int d = -1; d <= 1; d += 2) { double nr = c[i].r + d * EPS; int t = topCircle(Point(c[i].c.x + nr * cos(mid), c[i].c.y + nr * sin(mid)), n); if (~t) vis.insert(t); } } } cout << vis.size() << endl; } return 0; }
——written by Lyon
相关文章推荐
- LA - 2572 - Viva Confetti
- 【LA 2572】Viva Confetti, Kanazawa 2002 (计算几何,圆)
- LA 2572 Viva Confetti Kanazawa - 2002/2003 平面上的圆盘
- LA 2572 Viva Confetti 离散化 *
- LA 2572 - Viva Confetti
- Viva Confetti UVALive - 2572
- UVALive - 2572 Viva Confetti 极角排序
- UVaLive2572 poj1418 UVa1308 Viva Confetti
- 《viva la vida》 歌词
- LA 2572 圆盘的相互覆盖问题,圆弧极角排序,中点代替圆弧,轻微扰动的影响判断
- POJ 1418 Viva Confetti
- poj 1418 Viva Confetti(计算几何,圆)
- POJ 1418 Viva Confetti
- poj 1418 || zoj 1696 Viva Confetti
- LA 3263 That Nice Euler Circuit (2D Geometry)
- POJ 1418 Viva Confetti(Japan 2002 Kanazawa)
- 圆盘覆盖,计算几何(圆盘问题,LA 2572)
- zoj 1696 Viva Confetti
- POJ 1418 Viva Confetti(Japan 2002 Kanazawa)
- three.js 源码注释(六十九)extras/geometries/CircleGeometry.js