LA 2572 Viva Confetti (Geometry.Circle)

最新推荐文章于 2024-03-09 12:25:24 发布

转载最新推荐文章于 2024-03-09 12:25:24 发布 · 116 阅读

0 ·

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/LyonLys/archive/2013/05/01/LA_2572_Lyon.html

本文详细介绍了如何解决一道几何题，通过将每个可见部分视为由小圆弧围成，并利用刘汝佳的代码进行深入理解。文章解释了这种解题方法的原理，包括如何通过移动小圆弧中点来确定哪些圆可以被看见，提供了一种直观且高效的方法来解决此类问题。

https://icpcarchive.ecs.baylor.edu/index.php?option=com_onlinejudge&Itemid=8&page=problem_stats&problemid=573&category=

　　一道很好的几何题。要求是给出一些圆，按顺序覆盖上去，问哪些圆是可以被看见的。

　　看刘汝佳的书，开始的时候不明白“每个可见部分都是由一些‘小圆弧’围城的”这样又怎么样，直到我看了他的代码以后才懂。其实意思就是，因为每个可见部分都是由小圆弧围成，所以，如果我们将小圆弧中点（不会是圆与圆的交点）相对于弧的位置向里或向外稍微移动一下，然后找到有哪个圆最后覆盖在在那上面，那么这个圆必定能被看见。

　　刚接触几何，表示还真想不到可以这样做，所以就觉得这个做法实在太妙了！

代码如下：

View Code

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <cmath>
  4 #include <set>
  5 #include <vector>
  6 #include <iostream>
  7 #include <algorithm>
  8 
  9 using namespace std;
 10 
 11 struct Point {
 12     double x, y;
 13     Point() {}
 14     Point(double x, double y) : x(x), y(y) {}
 15 } ;
 16 template<class T> T sqr(T x) { return x * x;}
 17 inline double ptDis(Point a, Point b) { return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));}
 18 
 19 // basic calculations
 20 typedef Point Vec;
 21 Vec operator + (Vec a, Vec b) { return Vec(a.x + b.x, a.y + b.y);}
 22 Vec operator - (Vec a, Vec b) { return Vec(a.x - b.x, a.y - b.y);}
 23 Vec operator * (Vec a, double p) { return Vec(a.x * p, a.y * p);}
 24 Vec operator / (Vec a, double p) { return Vec(a.x / p, a.y / p);}
 25 
 26 const double EPS = 5e-13;
 27 const double PI = acos(-1.0);
 28 inline int sgn(double x) { return fabs(x) < EPS ? 0 : (x < 0 ? -1 : 1);}
 29 bool operator < (Point a, Point b) { return a.x < b.x || (a.x == b.x && a.y < b.y);}
 30 bool operator == (Point a, Point b) { return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0;}
 31 
 32 inline double dotDet(Vec a, Vec b) { return a.x * b.x + a.y * b.y;}
 33 inline double vecLen(Vec x) { return sqrt(sqr(x.x) + sqr(x.y));}
 34 inline double angle(Vec v) { return atan2(v.y, v.x);}
 35 inline double angle(Vec a, Vec b) { return acos(dotDet(a, b) / vecLen(a) / vecLen(b));}
 36 inline double crossDet(Vec a, Vec b) { return a.x * b.y - a.y * b.x;}
 37 inline double triArea(Point a, Point b, Point c) { return fabs(crossDet(b - a, c - a));}
 38 inline Vec vecUnit(Vec x) { return x / vecLen(x);}
 39 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));}
 40 Vec normal(Vec x) {
 41     double len = vecLen(x);
 42     return Vec(- x.y / len, x.x / len);
 43 }
 44 
 45 struct Line {
 46     Point s, t;
 47     Line() {}
 48     Line(Point s, Point t) : s(s), t(t) {}
 49     Point point(double x) {
 50         return s + (t - s) * x;
 51     }
 52     Line move(double x) {
 53         Vec nor = normal(t - s);
 54         nor = nor / vecLen(nor) * x;
 55         return Line(s + nor, t + nor);
 56     }
 57     Vec vec() { return t - s;}
 58 } ;
 59 typedef Line Seg;
 60 
 61 inline bool onSeg(Point x, Point a, Point b) { return sgn(crossDet(a - x, b - x)) == 0 && sgn(dotDet(a - x, b - x)) < 0;}
 62 inline bool onSeg(Point x, Seg s) { return onSeg(x, s.s, s.t);}
 63 
 64 // 0 : not intersect
 65 // 1 : proper intersect
 66 // 2 : improper intersect
 67 int segIntersect(Point a, Point c, Point b, Point d) {
 68     Vec v1 = b - a, v2 = c - b, v3 = d - c, v4 = a - d;
 69     int a_bc = sgn(crossDet(v1, v2));
 70     int b_cd = sgn(crossDet(v2, v3));
 71     int c_da = sgn(crossDet(v3, v4));
 72     int d_ab = sgn(crossDet(v4, v1));
 73     if (a_bc * c_da > 0 && b_cd * d_ab > 0) return 1;
 74     if (onSeg(b, a, c) && c_da) return 2;
 75     if (onSeg(c, b, d) && d_ab) return 2;
 76     if (onSeg(d, c, a) && a_bc) return 2;
 77     if (onSeg(a, d, b) && b_cd) return 2;
 78     return 0;
 79 }
 80 inline int segIntersect(Seg a, Seg b) { return segIntersect(a.s, a.t, b.s, b.t);}
 81 
 82 // point of the intersection of 2 lines
 83 Point lineIntersect(Point P, Vec v, Point Q, Vec w) {
 84     Vec u = P - Q;
 85     double t = crossDet(w, u) / crossDet(v, w);
 86     return P + v * t;
 87 }
 88 inline Point lineIntersect(Line a, Line b) { return lineIntersect(a.s, a.t - a.s, b.s, b.t - b.s);}
 89 
 90 // Warning: This is a DIRECTED Distance!!!
 91 double pt2Line(Point x, Point a, Point b) {
 92     Vec v1 = b - a, v2 = x - a;
 93     return crossDet(v1, v2) / vecLen(v1);
 94 }
 95 inline double pt2Line(Point x, Line L) { return pt2Line(x, L.s, L.t);}
 96 
 97 double pt2Seg(Point x, Point a, Point b) {
 98     if (a == b) return vecLen(x - a);
 99     Vec v1 = b - a, v2 = x - a, v3 = x - b;
100     if (sgn(dotDet(v1, v2)) < 0) return vecLen(v2);
101     if (sgn(dotDet(v1, v3)) > 0) return vecLen(v3);
102     return fabs(crossDet(v1, v2)) / vecLen(v1);
103 }
104 inline double pt2Seg(Point x, Seg s) { return pt2Seg(x, s.s, s.t);}
105 
106 struct Poly {
107     vector<Point> pt;
108     Poly() {}
109     Poly(vector<Point> pt) : pt(pt) {}
110     double area() {
111         double ret = 0.0;
112         int sz = pt.size();
113         for (int i = 1; i < sz; i++) {
114             ret += crossDet(pt[i], pt[i - 1]);
115         }
116         return fabs(ret / 2.0);
117     }
118 } ;
119 
120 struct Circle {
121     Point c;
122     double r;
123     Circle() {}
124     Circle(Point c, double r) : c(c), r(r) {}
125     Point point(double a) {
126         return Point(c.x + cos(a) * r, c.y + sin(a) * r);
127     }
128 } ;
129 
130 int lineCircleIntersect(Line L, Circle C, double &t1, double &t2, vector<Point> &sol) {
131     double a = L.s.x, b = L.t.x - C.c.x, c = L.s.y, d = L.t.y - C.c.y;
132     double e = sqr(a) + sqr(c), f = 2 * (a * b + c * d), g = sqr(b) + sqr(d) - sqr(C.r);
133     double delta = sqr(f) - 4.0 * e * g;
134     if (sgn(delta) < 0) return 0;
135     if (sgn(delta) == 0) {
136         t1 = t2 = -f / (2.0 * e);
137         sol.push_back(L.point(t1));
138         return 1;
139     }
140     t1 = (-f - sqrt(delta)) / (2.0 * e);
141     sol.push_back(L.point(t1));
142     t2 = (-f + sqrt(delta)) / (2.0 * e);
143     sol.push_back(L.point(t2));
144     return 2;
145 }
146 
147 int lineCircleIntersect(Line L, Circle C, vector<Point> &sol) {
148     Vec dir = L.t - L.s, nor = normal(dir);
149     Point mid = lineIntersect(C.c, nor, L.s, dir);
150     double len = sqr(C.r) - sqr(ptDis(C.c, mid));
151     if (sgn(len) < 0) return 0;
152     if (sgn(len) == 0) {
153         sol.push_back(mid);
154         return 1;
155     }
156     Vec dis = vecUnit(dir);
157     len = sqrt(len);
158     sol.push_back(mid + dis * len);
159     sol.push_back(mid - dis * len);
160     return 2;
161 }
162 
163 // -1 : coincide
164 int circleCircleIntersect(Circle C1, Circle C2, vector<Point> &sol) {
165     double d = vecLen(C1.c - C2.c);
166     if (sgn(d) == 0) {
167         if (sgn(C1.r - C2.r) == 0) {
168             return -1;
169         }
170         return 0;
171     }
172     if (sgn(C1.r + C2.r - d) < 0) return 0;
173     if (sgn(fabs(C1.r - C2.r) - d) > 0) return 0;
174     double a = angle(C2.c - C1.c);
175     double da = acos((sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2.0 * C1.r * d));
176     Point p1 = C1.point(a - da), p2 = C1.point(a + da);
177     sol.push_back(p1);
178     if (p1 == p2) return 1;
179     sol.push_back(p2);
180     return 2;
181 }
182 
183 void circleCircleIntersect(Circle C1, Circle C2, vector<double> &sol) {
184     double d = vecLen(C1.c - C2.c);
185     if (sgn(d) == 0) return ;
186     if (sgn(C1.r + C2.r - d) < 0) return ;
187     if (sgn(fabs(C1.r - C2.r) - d) > 0) return ;
188     double a = angle(C2.c - C1.c);
189     double da = acos((sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2.0 * C1.r * d));
190     sol.push_back(a - da);
191     sol.push_back(a + da);
192 }
193 
194 int tangent(Point p, Circle C, vector<Vec> &sol) {
195     Vec u = C.c - p;
196     double dist = vecLen(u);
197     if (dist < C.r) return 0;
198     if (sgn(dist - C.r) == 0) {
199         sol.push_back(rotate(u, PI / 2.0));
200         return 1;
201     }
202     double ang = asin(C.r / dist);
203     sol.push_back(rotate(u, -ang));
204     sol.push_back(rotate(u, ang));
205     return 2;
206 }
207 
208 // ptA : points of tangency on circle A
209 // ptB : points of tangency on circle B
210 int tangent(Circle A, Circle B, vector<Point> &ptA, vector<Point> &ptB) {
211     if (A.r < B.r) {
212         swap(A, B);
213         swap(ptA, ptB);
214     }
215     int d2 = sqr(A.c.x - B.c.x) + sqr(A.c.y - B.c.y);
216     int rdiff = A.r - B.r, rsum = A.r + B.r;
217     if (d2 < sqr(rdiff)) return 0;
218     double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
219     if (d2 == 0 && A.r == B.r) return -1;
220     if (d2 == sqr(rdiff)) {
221         ptA.push_back(A.point(base));
222         ptB.push_back(B.point(base));
223         return 1;
224     }
225     double ang = acos((A.r - B.r) / sqrt(d2));
226     ptA.push_back(A.point(base + ang));
227     ptB.push_back(B.point(base + ang));
228     ptA.push_back(A.point(base - ang));
229     ptB.push_back(B.point(base - ang));
230     if (d2 == sqr(rsum)) {
231         ptA.push_back(A.point(base));
232         ptB.push_back(B.point(PI + base));
233     } else if (d2 > sqr(rsum)) {
234         ang = acos((A.r + B.r) / sqrt(d2));
235         ptA.push_back(A.point(base + ang));
236         ptB.push_back(B.point(PI + base + ang));
237         ptA.push_back(A.point(base - ang));
238         ptB.push_back(B.point(PI + base - ang));
239     }
240     return (int) ptA.size();
241 }
242 
243 /****************** template above *******************/
244 
245 #define DEC(i, a, b) for (int i = (a); i >= (b); i--)
246 #define REP(i, n) for (int i = 0; i < (n); i++)
247 #define _clr(x) memset(x, 0, sizeof(x))
248 #define PB push_back
249 Circle c[111];
250 
251 int topCircle(Point x, int n) {
252     DEC(i, n - 1, 0) {
253         if (ptDis(c[i].c, x) < c[i].r) return i;
254     }
255     return -1;
256 }
257 
258 int main() {
259 //    freopen("in", "r", stdin);
260     int n;
261     while (cin >> n && n) {
262         REP(i, n) cin >> c[i].c.x >> c[i].c.y >> c[i].r;
263         set<int> vis;
264         vis.clear();
265         REP(i, n) {
266             vector<double> pos;
267             pos.clear();
268             pos.PB(0.0);
269             pos.PB(2.0 * PI);
270             REP(j, n) {
271                 circleCircleIntersect(c[i], c[j], pos);
272             }
273             sort(pos.begin(), pos.end());
274             int sz = pos.size() - 1;
275             REP(j, sz) {
276                 double mid = (pos[j] + pos[j + 1]) / 2.0;
277                 for (int d = -1; d <= 1; d += 2) {
278                     double nr = c[i].r + d * EPS;
279                     int t = topCircle(Point(c[i].c.x + nr * cos(mid), c[i].c.y + nr * sin(mid)), n);
280                     if (~t) vis.insert(t);
281                 }
282             }
283         }
284         cout << vis.size() << endl;
285     }
286     return 0;
287 }

——written by Lyon

转载于:https://www.cnblogs.com/LyonLys/archive/2013/05/01/LA_2572_Lyon.html