以圆心为中心将简单多边形划分为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 }