HDU2892 area 简单多边形与圆面积交

以圆心为中心将简单多边形划分为n个矢量三角形,对每个三角形与圆求交,根据有向边判断相交面积正负,最后相加取绝对值。

一个顶点在圆心的三角形与圆的交需要讨论的情况比较少,容易计算。

  1 #include<stdio.h>
  2 #include<string.h>
  3 #include<stdlib.h>
  4 #include<math.h>
  5 #include<algorithm>
  6 const int maxn = 111111;
  7 const int maxisn = 21;
  8 const double eps = 1e-8;
  9 const double pi = acos(-1.0);
 10 int dcmp(double x)
 11 {
 12     if(x > eps) return 1;
 13     return x < -eps ? -1 : 0;
 14 }
 15 struct Point
 16 {
 17     double x, y;
 18     Point(){x = y = 0;}
 19     Point(double a, double b)
 20     {x = a, y = b;}
 21     inline Point operator-(const Point &b)const
 22     {return Point(x - b.x, y - b.y);}
 23     inline Point operator+(const Point &b)const
 24     {return Point(x + b.x, y + b.y);}
 25     inline Point operator*(const double &b)const
 26     {return Point(x * b, y * b);}
 27     inline double dot(const Point &b)const
 28     {return x * b.x + y * b.y;}
 29     inline double cross(const Point &b, const Point &c)const
 30     {return (b.x - x) * (c.y - y) - (c.x - x) * (b.y - y);}
 31     inline double Dis(const Point &b)const
 32     {return sqrt((*this - b).dot(*this - b));}
 33     inline bool InLine(const Point &b, const Point &c)const//三点共线
 34     {return !dcmp(cross(b, c));}
 35     inline bool OnSeg(const Point &b, const Point &c)const//点在线段上,包括端点
 36     {return InLine(b, c) && (*this - c).dot(*this - b) < eps;}
 37 };
 38 inline double min(double a, double b)
 39 {return a < b ? a : b;}
 40 inline double max(double a, double b)
 41 {return a > b ? a : b;}
 42 inline double Sqr(double x)
 43 {return x * x;}
 44 inline double Sqr(const Point &p)
 45 {return p.dot(p);}
 46 Point LineCross(const Point &a, const Point &b, const Point &c, const Point &d)
 47 {
 48     double u = a.cross(b, c), v = b.cross(a, d);
 49     return Point((c.x * v + d.x * u) / (u + v), (c.y * v + d.y * u) / (u + v));
 50 }
 51 double LineCrossCircle(const Point &a, const Point &b, const Point &r, 
 52             double R, Point &p1, Point &p2)
 53 {
 54     Point fp = LineCross(r, Point(r.x + a.y - b.y, r.y + b.x - a.x), a, b);
 55     double rtol = r.Dis(fp);
 56     double rtos = fp.OnSeg(a, b) ? rtol : min(r.Dis(a), r.Dis(b));
 57     double atob = a.Dis(b);
 58     double fptoe = sqrt(R * R - rtol * rtol) / atob;
 59     if(rtos > R - eps) return rtos;
 60     p1 = fp + (a - b) * fptoe;
 61     p2 = fp + (b - a) * fptoe;
 62     return rtos;
 63 }
 64 double SectorArea(const Point &r, const Point &a, const Point &b, double R)
 65 //不大于180度扇形面积,r->a->b逆时针
 66 {
 67     double A2 = Sqr(r - a), B2 = Sqr(r - b), C2 = Sqr(a - b);
 68     return R * R * acos((A2 + B2 - C2) * 0.5 / sqrt(A2) / sqrt(B2)) * 0.5;
 69 }
 70 double TACIA(const Point &r, const Point &a, const Point &b, double R)
 71 //TriangleAndCircleIntersectArea,逆时针,r为圆心
 72 {
 73     double adis = r.Dis(a), bdis = r.Dis(b);
 74     if(adis < R + eps && bdis < R + eps) return r.cross(a, b) * 0.5;
 75     Point ta, tb;
 76     if(r.InLine(a, b)) return 0.0;
 77     double rtos = LineCrossCircle(a, b, r, R, ta, tb);
 78     if(rtos > R - eps) return SectorArea(r, a, b, R);
 79     if(adis < R + eps) return r.cross(a, tb) * 0.5 + SectorArea(r, tb, b, R);
 80     if(bdis < R + eps) return r.cross(ta, b) * 0.5 + SectorArea(r, a, ta, R);
 81     return r.cross(ta, tb) * 0.5 + 
 82         SectorArea(r, a, ta, R) + SectorArea(r, tb, b, R);
 83 }
 84 double SPICA(int n, Point r, double R)//SimplePolygonIntersectCircleArea
 85 {
 86     int i;
 87     Point ori, p[2];
 88     scanf("%lf%lf", &ori.x, &ori.y);
 89     p[0] = ori;
 90     double res = 0, if_clock_t;
 91     for(i = 1; i <= n; ++ i)
 92     {
 93         if(i == n) p[i & 1] = ori;
 94         else scanf("%lf%lf", &p[i & 1].x, &p[i & 1].y);
 95         if_clock_t = dcmp(r.cross(p[~i & 1], p[i & 1]));
 96         if(if_clock_t < 0) res -= TACIA(r, p[i & 1], p[~i & 1], R);
 97         else res += TACIA(r, p[~i & 1], p[i & 1], R);
 98     }
 99     return fabs(res);
100 }
101 Point boom;
102 int n;
103 double R;
104 int main()
105 {
106     double sx, sy, h, vx, vy;
107     while(scanf("%lf%lf%lf", &sx, &sy, &h) != EOF)
108     {
109         scanf("%lf%lf%lf", &vx, &vy, &R);
110         h = sqrt(2 * h / 10);
111         boom = Point(h * vx + sx, h * vy + sy);
112         scanf("%d", &n);
113         printf("%.2f\n", SPICA(n, boom, R));
114     }
115     return 0;
116 }
posted @ 2012-09-08 23:35  CSGrandeur  阅读(1091)  评论(0编辑  收藏  举报