Description
Given one triangle and one circle in the plane. Your task is to calculate the common area of these two figures.
Input
The input will contain several test cases. Each line of input describes a test case. Each test case consists of nine floating point numbers, x1, y1, x2, y2, x3, y3, x4, y4 and r, where (x1, y1), (x2, y2) and (x3, y3) are the three vertices of the triangle and (x4, y4) is the center of the circle and r is the radius. We guarantee the triangle and the circle are not degenerate.
Output
For each test case you should output one real number, which is the common area of the triangle and the circle, on a separate line. The result should be rounded to two decimal places.
题目大意：求一个三角形和一个圆形的交的面积。
思路：圆心和三个三角形的三个点连线，把一个三角形划分为3个三角形，利用有向面积来算。然后就变成了求一个三角形和圆的交的面积，其中三角形的一个顶点为圆心。然后各种分情况讨论（但是要分的情况起码会比直接算一个普通三角形和圆形的交要少得多）。我的姿势似乎不是很高级，好像有些什么奇怪的公式……
PS：调了一下发现居然是以前用的模板错了……
代码（1938MS）：
1 #include2 #include 3 #include 4 #include 5 #include 6 using namespace std; 7 8 const int MAXN = 10010; 9 const double EPS = 1e-10; 10 const double PI = acos(-1.0);//3.14159265358979323846 11 const double INF = 1; 12 13 inline int sgn(double x) { 14 return (x > EPS) - (x < -EPS); 15 } 16 17 inline double sqr(double x) { 18 return x * x; 19 } 20 21 struct Point { 22 double x, y, ag; 23 Point() {} 24 Point(double x, double y): x(x), y(y) {} 25 void read() { 26 scanf("%lf%lf", &x, &y); 27 } 28 bool operator == (const Point &rhs) const { 29 return sgn(x - rhs.x) == 0 && sgn(y - rhs.y) == 0; 30 } 31 bool operator < (const Point &rhs) const { 32 if(y != rhs.y) return y < rhs.y; 33 return x < rhs.x; 34 } 35 Point operator + (const Point &rhs) const { 36 return Point(x + rhs.x, y + rhs.y); 37 } 38 Point operator - (const Point &rhs) const { 39 return Point(x - rhs.x, y - rhs.y); 40 } 41 Point operator * (const double &b) const { 42 return Point(x * b, y * b); 43 } 44 Point operator / (const double &b) const { 45 return Point(x / b, y / b); 46 } 47 double operator * (const Point &rhs) const { 48 return x * rhs.x + y * rhs.y; 49 } 50 double length() { 51 return sqrt(x * x + y * y); 52 } 53 double angle() { 54 return atan2(y, x); 55 } 56 Point unit() { 57 return *this / length(); 58 } 59 void makeAg() { 60 ag = atan2(y, x); 61 } 62 void print() { 63 printf("%.10f %.10f\n", x, y); 64 } 65 }; 66 typedef Point Vector; 67 68 double dist(const Point &a, const Point &b) { 69 return (a - b).length(); 70 } 71 72 double cross(const Point &a, const Point &b) { 73 return a.x * b.y - a.y * b.x; 74 } 75 //ret >= 0 means turn right 76 double cross(const Point &sp, const Point &ed, const Point &op) { 77 return cross(sp - op, ed - op); 78 } 79 80 double area(const Point& a, const Point &b, const Point &c) { 81 return fabs(cross(a - c, b - c)) / 2; 82 } 83 //counter-clockwise 84 Point rotate(const Point &p, double angle, const Point &o = Point(0, 0)) { 85 Point t = p - o; 86 double x = t.x * cos(angle) - t.y * sin(angle); 87 double y = t.y * cos(angle) + t.x * sin(angle); 88 return Point(x, y) + o; 89 } 90 91 double includedAngle(const Point &a, const Point &b, const Point &o) { 92 double ret = abs((a - o).angle() - (b - o).angle()); 93 if(sgn(ret - PI) > 0) ret = 2 * PI - ret; 94 return ret; 95 } 96 97 struct Seg { 98 Point st, ed; 99 double ag;100 Seg() {}101 Seg(Point st, Point ed): st(st), ed(ed) {}102 void read() {103 st.read(); ed.read();104 }105 void makeAg() {106 ag = atan2(ed.y - st.y, ed.x - st.x);107 }108 };109 typedef Seg Line;110 111 //ax + by + c > 0112 Line buildLine(double a, double b, double c) {113 if(sgn(a) == 0 && sgn(b) == 0) return Line(Point(sgn(c) > 0 ? -1 : 1, INF), Point(0, INF));114 if(sgn(a) == 0) return Line(Point(sgn(b), -c/b), Point(0, -c/b));115 if(sgn(b) == 0) return Line(Point(-c/a, 0), Point(-c/a, sgn(a)));116 if(b < 0) return Line(Point(0, -c/b), Point(1, -(a + c) / b));117 else return Line(Point(1, -(a + c) / b), Point(0, -c/b));118 }119 120 void moveRight(Line &v, double r) {121 double dx = v.ed.x - v.st.x, dy = v.ed.y - v.st.y;122 dx = dx / dist(v.st, v.ed) * r;123 dy = dy / dist(v.st, v.ed) * r;124 v.st.x += dy; v.ed.x += dy;125 v.st.y -= dx; v.ed.y -= dx;126 }127 128 bool isOnSeg(const Seg &s, const Point &p) {129 return (p == s.st || p == s.ed) ||130 (((p.x - s.st.x) * (p.x - s.ed.x) < 0 ||131 (p.y - s.st.y) * (p.y - s.ed.y) < 0) &&132 sgn(cross(s.ed, p, s.st)) == 0);133 }134 135 bool isIntersected(const Point &s1, const Point &e1, const Point &s2, const Point &e2) {136 return (max(s1.x, e1.x) >= min(s2.x, e2.x)) &&137 (max(s2.x, e2.x) >= min(s1.x, e1.x)) &&138 (max(s1.y, e1.y) >= min(s2.y, e2.y)) &&139 (max(s2.y, e2.y) >= min(s1.y, e1.y)) &&140 (cross(s2, e1, s1) * cross(e1, e2, s1) >= 0) &&141 (cross(s1, e2, s2) * cross(e2, e1, s2) >= 0);142 }143 144 bool isIntersected(const Seg &a, const Seg &b) {145 return isIntersected(a.st, a.ed, b.st, b.ed);146 }147 148 bool isParallel(const Seg &a, const Seg &b) {149 return sgn(cross(a.ed - a.st, b.ed - b.st)) == 0;150 }151 152 //return Ax + By + C =0 's A, B, C153 void Coefficient(const Line &L, double &A, double &B, double &C) {154 A = L.ed.y - L.st.y;155 B = L.st.x - L.ed.x;156 C = L.ed.x * L.st.y - L.st.x * L.ed.y;157 }158 //point of intersection159 Point operator * (const Line &a, const Line &b) {160 double A1, B1, C1;161 double A2, B2, C2;162 Coefficient(a, A1, B1, C1);163 Coefficient(b, A2, B2, C2);164 Point I;165 I.x = - (B2 * C1 - B1 * C2) / (A1 * B2 - A2 * B1);166 I.y = (A2 * C1 - A1 * C2) / (A1 * B2 - A2 * B1);167 return I;168 }169 170 bool isEqual(const Line &a, const Line &b) {171 double A1, B1, C1;172 double A2, B2, C2;173 Coefficient(a, A1, B1, C1);174 Coefficient(b, A2, B2, C2);175 return sgn(A1 * B2 - A2 * B1) == 0 && sgn(A1 * C2 - A2 * C1) == 0 && sgn(B1 * C2 - B2 * C1) == 0;176 }177 178 double Point_to_Line(const Point &p, const Line &L) {179 return fabs(cross(p, L.st, L.ed)/dist(L.st, L.ed));180 }181 182 double Point_to_Seg(const Point &p, const Seg &L) {183 if(sgn((L.ed - L.st) * (p - L.st)) < 0) return dist(p, L.st);184 if(sgn((L.st - L.ed) * (p - L.ed)) < 0) return dist(p, L.ed);185 return Point_to_Line(p, L);186 }187 188 double Seg_to_Seg(const Seg &a, const Seg &b) {189 double ans1 = min(Point_to_Seg(a.st, b), Point_to_Seg(a.ed, b));190 double ans2 = min(Point_to_Seg(b.st, a), Point_to_Seg(b.ed, a));191 return min(ans1, ans2);192 }193 194 struct Circle {195 Point c;196 double r;197 Circle() {}198 Circle(Point c, double r): c(c), r(r) {}199 void read() {200 c.read();201 scanf("%lf", &r);202 }203 double area() const {204 return PI * r * r;205 }206 bool contain(const Circle &rhs) const {207 return sgn(dist(c, rhs.c) + rhs.r - r) <= 0;208 }209 bool contain(const Point &p) const {210 return sgn(dist(c, p) - r) <= 0;211 }212 bool intersect(const Circle &rhs) const {213 return sgn(dist(c, rhs.c) - r - rhs.r) < 0;214 }215 bool tangency(const Circle &rhs) const {216 return sgn(dist(c, rhs.c) - r - rhs.r) == 0;217 }218 Point pos(double angle) const {219 Point p = Point(c.x + r, c.y);220 return rotate(p, angle, c);221 }222 };223 224 double CommonArea(const Circle &A, const Circle &B) {225 double area = 0.0;226 const Circle & M = (A.r > B.r) ? A : B;227 const Circle & N = (A.r > B.r) ? B : A;228 double D = dist(M.c, N.c);229 if((D < M.r + N.r) && (D > M.r - N.r)) {230 double cosM = (M.r * M.r + D * D - N.r * N.r) / (2.0 * M.r * D);231 double cosN = (N.r * N.r + D * D - M.r * M.r) / (2.0 * N.r * D);232 double alpha = 2 * acos(cosM);233 double beta = 2 * acos(cosN);234 double TM = 0.5 * M.r * M.r * (alpha - sin(alpha));235 double TN = 0.5 * N.r * N.r * (beta - sin(beta));236 area = TM + TN;237 }238 else if(D <= M.r - N.r) {239 area = N.area();240 }241 return area;242 }243 244 int intersection(const Seg &s, const Circle &cir, Point &p1, Point &p2) {245 double angle = includedAngle(s.ed, cir.c, s.st);246 double B = dist(cir.c, s.st);247 double a = 1, b = -2 * B * cos(angle), c = sqr(B) - sqr(cir.r);248 double delta = sqr(b) - 4 * a * c;249 if(sgn(delta) < 0) return 0;250 double x1 = (-b - sqrt(delta)) / (2 * a), x2 = (-b + sqrt(delta)) / (2 * a);251 Vector v = (s.ed - s.st).unit();252 p1 = s.st + v * x1;253 p2 = s.st + v * x2;254 return 1 + sgn(delta);255 }256 257 double CommonArea(const Circle &cir, Point p1, Point p2) {258 if(cir.contain(p1) && cir.contain(p2)) {259 return area(cir.c, p1, p2);260 } else if(!cir.contain(p1) && !cir.contain(p2)) {261 Point q1, q2;262 int t = intersection(Line(p1, p2), cir, q1, q2);263 if(t == 0) {264 double angle = includedAngle(p1, p2, cir.c);265 return 0.5 * sqr(cir.r) * angle;266 } else {267 double angle1 = includedAngle(p1, p2, cir.c);268 double angle2 = includedAngle(q1, q2, cir.c);269 if(isOnSeg(Seg(p1, p2), q1))return 0.5 * sqr(cir.r) * (angle1 - angle2 + sin(angle2));270 else return 0.5 * sqr(cir.r) * angle1;271 }272 } else {273 if(cir.contain(p2)) swap(p1, p2);274 Point q1, q2;275 intersection(Line(p1, p2), cir, q1, q2);276 double angle = includedAngle(q2, p2, cir.c);277 double a = area(cir.c, p1, q2);278 double b = 0.5 * sqr(cir.r) * angle;279 return a + b;280 }281 }282 283 struct Triangle {284 Point p[3];285 Triangle() {}286 Triangle(Point *t) {287 for(int i = 0; i < 3; ++i) p[i] = t[i];288 }289 void read() {290 for(int i = 0; i < 3; ++i) p[i].read();291 }292 double area() const {293 return ::area(p[0], p[1], p[2]);294 }295 Point& operator[] (int i) {296 return p[i];297 }298 };299 300 double CommonArea(Triangle tir, const Circle &cir) {301 double ret = 0;302 ret += sgn(cross(tir[0], cir.c, tir[1])) * CommonArea(cir, tir[0], tir[1]);303 ret += sgn(cross(tir[1], cir.c, tir[2])) * CommonArea(cir, tir[1], tir[2]);304 ret += sgn(cross(tir[2], cir.c, tir[0])) * CommonArea(cir, tir[2], tir[0]);305 return abs(ret);306 }307 308 struct Poly {309 int n;310 Point p[MAXN];//p[n] = p[0]311 void init(Point *pp, int nn) {312 n = nn;313 for(int i = 0; i < n; ++i) p[i] = pp[i];314 p[n] = p[0];315 }316 double area() {317 if(n < 3) return 0;318 double s = p[0].y * (p[n - 1].x - p[1].x);319 for(int i = 1; i < n; ++i)320 s += p[i].y * (p[i - 1].x - p[i + 1].x);321 return s / 2;322 }323 };324 //the convex hull is clockwise325 void Graham_scan(Point *p, int n, int *stk, int &top) { //stk[0] = stk[top]326 sort(p, p + n);327 top = 1;328 stk[0] = 0; stk[1] = 1;329 for(int i = 2; i < n; ++i) {330 while(top && cross(p[i], p[stk[top]], p[stk[top - 1]]) <= 0) --top;331 stk[++top] = i;332 }333 int len = top;334 stk[++top] = n - 2;335 for(int i = n - 3; i >= 0; --i) {336 while(top != len && cross(p[i], p[stk[top]], p[stk[top - 1]]) <= 0) --top;337 stk[++top] = i;338 }339 }340 //use for half_planes_cross341 bool cmpAg(const Line &a, const Line &b) {342 if(sgn(a.ag - b.ag) == 0)343 return sgn(cross(b.ed, a.st, b.st)) < 0;344 return a.ag < b.ag;345 }346 //clockwise, plane is on the right347 bool half_planes_cross(Line *v, int vn, Poly &res, Line *deq) {348 int i, n;349 sort(v, v + vn, cmpAg);350 for(i = n = 1; i < vn; ++i) {351 if(sgn(v[i].ag - v[i-1].ag) == 0) continue;352 v[n++] = v[i];353 }354 int head = 0, tail = 1;355 deq[0] = v[0], deq[1] = v[1];356 for(i = 2; i < n; ++i) {357 if(isParallel(deq[tail - 1], deq[tail]) || isParallel(deq[head], deq[head + 1]))358 return false;359 while(head < tail && sgn(cross(v[i].ed, deq[tail - 1] * deq[tail], v[i].st)) > 0)360 --tail;361 while(head < tail && sgn(cross(v[i].ed, deq[head] * deq[head + 1], v[i].st)) > 0)362 ++head;363 deq[++tail] = v[i];364 }365 while(head < tail && sgn(cross(deq[head].ed, deq[tail - 1] * deq[tail], deq[head].st)) > 0)366 --tail;367 while(head < tail && sgn(cross(deq[tail].ed, deq[head] * deq[head + 1], deq[tail].st)) > 0)368 ++head;369 if(tail <= head + 1) return false;370 res.n = 0;371 for(i = head; i < tail; ++i)372 res.p[res.n++] = deq[i] * deq[i + 1];373 res.p[res.n++] = deq[head] * deq[tail];374 res.n = unique(res.p, res.p + res.n) - res.p;375 res.p[res.n] = res.p[0];376 return true;377 }378 379 //ix and jx is the points whose distance is return, res.p[n - 1] = res.p[0], res must be clockwise380 double dia_rotating_calipers(Poly &res, int &ix, int &jx) {381 double dia = 0;382 int q = 1;383 for(int i = 0; i < res.n - 1; ++i) {384 while(sgn(cross(res.p[i], res.p[q + 1], res.p[i + 1]) - cross(res.p[i], res.p[q], res.p[i + 1])) > 0)385 q = (q + 1) % (res.n - 1);386 if(sgn(dist(res.p[i], res.p[q]) - dia) > 0) {387 dia = dist(res.p[i], res.p[q]);388 ix = i; jx = q;389 }390 if(sgn(dist(res.p[i + 1], res.p[q]) - dia) > 0) {391 dia = dist(res.p[i + 1], res.p[q]);392 ix = i + 1; jx = q;393 }394 }395 return dia;396 }397 //a and b must be clockwise, find the minimum distance between two convex hull398 double half_rotating_calipers(Poly &a, Poly &b) {399 int sa = 0, sb = 0;400 for(int i = 0; i < a.n; ++i) if(sgn(a.p[i].y - a.p[sa].y) < 0) sa = i;401 for(int i = 0; i < b.n; ++i) if(sgn(b.p[i].y - b.p[sb].y) < 0) sb = i;402 double tmp, ans = dist(a.p[0], b.p[0]);403 for(int i = 0; i < a.n; ++i) {404 while(sgn(tmp = cross(a.p[sa], a.p[sa + 1], b.p[sb + 1]) - cross(a.p[sa], a.p[sa + 1], b.p[sb])) > 0)405 sb = (sb + 1) % (b.n - 1);406 if(sgn(tmp) < 0) ans = min(ans, Point_to_Seg(b.p[sb], Seg(a.p[sa], a.p[sa + 1])));407 else ans = min(ans, Seg_to_Seg(Seg(a.p[sa], a.p[sa + 1]), Seg(b.p[sb], b.p[sb + 1])));408 sa = (sa + 1) % (a.n - 1);409 }410 return ans;411 }412 413 double rotating_calipers(Poly &a, Poly &b) {414 return min(half_rotating_calipers(a, b), half_rotating_calipers(b, a));415 }416 417 /*******************************************************************************************/418 419 Triangle tir;420 Circle cir;421 422 int main() {423 while(scanf("%lf%lf", &tir[0].x, &tir[0].y) != EOF) {424 tir[1].read();425 tir[2].read();426 cir.read();427 //cout<