首页 > 技术文章 > 三角剖分求多边形面积的交 HDU3060

ITUPC 2016-09-20 23:31 原文

  1 //三角剖分求多边形面积的交 HDU3060
  2 
  3 #include <iostream>
  4 #include <cstdio>
  5 #include <cstring>
  6 #include <stack>
  7 #include <queue>
  8 #include <cmath>
  9 #include <algorithm>
 10 using namespace std;
 11 
 12 const int maxn = 555;
 13 const int maxisn = 10;
 14 const double eps = 1e-8;
 15 const double pi = acos(-1.0);
 16 
 17 int dcmp(double x) {
 18     if(x > eps) return 1;
 19     return x < -eps ? -1 : 0;
 20 }
 21 
 22 struct Point {
 23     double x, y;
 24     Point() {
 25         x = y = 0;
 26     }
 27     Point(double a, double b) {
 28         x = a, y = b;
 29     }
 30     inline Point operator-(const Point &b)const {
 31         return Point(x - b.x, y - b.y);
 32     }
 33     inline Point operator+(const Point &b)const {
 34         return Point(x + b.x, y + b.y);
 35     }
 36     inline double dot(const Point &b)const {
 37         return x * b.x + y * b.y;
 38     }
 39     inline double cross(const Point &b, const Point &c)const {
 40         return (b.x - x) * (c.y - y) - (c.x - x) * (b.y - y);
 41     }
 42 };
 43 
 44 Point LineCross(const Point &a, const Point &b, const Point &c, const Point &d) {
 45     double u = a.cross(b, c), v = b.cross(a, d);
 46     return Point((c.x * v + d.x * u) / (u + v), (c.y * v + d.y * u) / (u + v));
 47 }
 48 
 49 double PolygonArea(Point p[], int n) {
 50     if(n < 3) return 0.0;
 51     double s = p[0].y * (p[n - 1].x - p[1].x);
 52     p[n] = p[0];
 53     for(int i = 1; i < n; ++ i)
 54         s += p[i].y * (p[i - 1].x - p[i + 1].x);
 55     return fabs(s * 0.5);
 56 }
 57 
 58 double CPIA(Point a[], Point b[], int na, int nb) { //ConvexPolygonIntersectArea
 59     Point p[maxisn], tmp[maxisn];
 60     int i, j, tn, sflag, eflag;
 61     a[na] = a[0], b[nb] = b[0];
 62     memcpy(p, b, sizeof(Point) * (nb + 1));
 63     for(i = 0; i < na && nb > 2; ++ i) {
 64         sflag = dcmp(a[i].cross(a[i + 1], p[0]));
 65         for(j = tn = 0; j < nb; ++ j, sflag = eflag) {
 66             if(sflag >= 0) tmp[tn ++] = p[j];
 67             eflag = dcmp(a[i].cross(a[i + 1], p[j + 1]));
 68             if((sflag ^ eflag) == -2)
 69                 tmp[tn ++] = LineCross(a[i], a[i + 1], p[j], p[j + 1]);
 70         }
 71         memcpy(p, tmp, sizeof(Point) * tn);
 72         nb = tn, p[nb] = p[0];
 73     }
 74     if(nb < 3) return 0.0;
 75     return PolygonArea(p, nb);
 76 }
 77 
 78 double SPIA(Point a[], Point b[], int na, int nb) { //SimplePolygonIntersectArea
 79     int i, j;
 80     Point t1[4], t2[4];
 81     double res = 0, if_clock_t1, if_clock_t2;
 82     a[na] = t1[0] = a[0], b[nb] = t2[0] = b[0];
 83     for(i = 2; i < na; ++ i) {
 84         t1[1] = a[i - 1], t1[2] = a[i];
 85         if_clock_t1 = dcmp(t1[0].cross(t1[1], t1[2]));
 86         if(if_clock_t1 < 0) std::swap(t1[1], t1[2]);
 87         for(j = 2; j < nb; ++ j) {
 88             t2[1] = b[j - 1], t2[2] = b[j];
 89             if_clock_t2 = dcmp(t2[0].cross(t2[1], t2[2]));
 90             if(if_clock_t2 < 0) std::swap(t2[1], t2[2]);
 91             res += CPIA(t1, t2, 3, 3) * if_clock_t1 * if_clock_t2;
 92         }
 93     }
 94     return PolygonArea(a, na) + PolygonArea(b, nb) - res;
 95 }
 96 
 97 Point p1[maxn], p2[maxn];
 98 int n1, n2;
 99 
100 int main() {
101     int i;
102     while(scanf("%d%d", &n1, &n2) != EOF) {
103         for(i = 0; i < n1; ++ i) scanf("%lf%lf", &p1[i].x, &p1[i].y);
104         for(i = 0; i < n2; ++ i) scanf("%lf%lf", &p2[i].x, &p2[i].y);
105         printf("%.2f\n", SPIA(p1, p2, n1, n2) + eps);
106     }
107     return 0;
108 }

 

推荐阅读