一道很好的几何题。要求是给出一些圆，按顺序覆盖上去，问哪些圆是可以被看见的。
看刘汝佳的书，开始的时候不明白“每个可见部分都是由一些‘小圆弧’围城的”这样又怎么样，直到我看了他的代码以后才懂。其实意思就是，因为每个可见部分都是由小圆弧围成，所以，如果我们将小圆弧中点（不会是圆与圆的交点）相对于弧的位置向里或向外稍微移动一下，然后找到有哪个圆最后覆盖在在那上面，那么这个圆必定能被看见。
刚接触几何，表示还真想不到可以这样做，所以就觉得这个做法实在太妙了！
代码如下：
View Code 1 #include2 #include 3 #include 4 #include 5 #include 6 #include 7 #include 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 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 pt;108 Poly() {}109 Poly(vector 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 &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 &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 : coincide164 int circleCircleIntersect(Circle C1, Circle C2, vector &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 &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 &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 A209 // ptB : points of tangency on circle B210 int tangent(Circle A, Circle B, vector &ptA, vector &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_back249 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 vis;264 vis.clear();265 REP(i, n) {266 vector 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