解题报告-HDU 2892 (计算几何)

http://acm.hdu.edu.cn/showproblem.php?pid=2892

感觉有点恶心的几何，搞了差不多一天才过。

题目说的是给出原始的坐标，还有飞行速度和重力加速度，然后求炮弹最后毁灭的面积。

首先是求炮弹落地的坐标，这个很好求，只要有点物理基础的人都能做出来，主要是怎么求圆和多边形的交面积。

枚举每两个相邻的顶点，判断这两个顶点和两个点构成的交线和圆的关系，主要有下面的4种情况：

1、p[i]和p[i+1]都在圆内，此时只需要求三角形面积就行了。

2、p[i]在圆内，p[i+1]在圆外，此时求出直线和圆的交点inter，然后求三角形面积（红色）加上扇形面积（蓝色）。

3、p[i]在圆外，p[i+1]在院内，和第二种情况类似，用三角形面积加上扇形面积。

4、p[i]和p[i+1]都在圆外，这时候又分为三种情况：

（1）p[i]和p[i+1]构成的直线和圆没有交点 

（2）p[i]和p[i+1]构成的直线和圆又一个交点

（3）p[i]和p[i+1]构成的直线和圆没有交点

很容易发现，第一和第二种情况实际上是一样的，我们可以把他们归在一类种处理，这种情况很好处理，直接求扇形面积就可以了。

对于第三种情况，我们仿照前面的做法，把图形剖分，先求出圆和直线的两个交点inter1和inter2，通过这两个交点对多边形进行划分，把红、蓝黑、三种颜色的面积加起来就是答案了。

但是，题目说的是不规则多边形，并不一定是凸多边形。因为用叉积算出来的面积是有向面积，对于三角形剖析求面积，求出的结果取绝对值就是面积了，但是如果直接用夹角求扇形面积，求出来的扇形面积并不是有向的，即恒为正。为了与有向面积统一，我给扇形面积也加上了一个方向，对于扇形圆弧上的两个端点，做一次叉积运算，把叉积的正负号给扇形面积，这样的扇形就有正负号了。

这种方法只需要遍历一次所有的顶点，时间复杂度为O(n)。

贴上一份自己敲的傻逼代码：

#include <iostream>#include <cmath>#include <cstdio>using namespace std;const int maxn = 100000 + 10;const double eps = 1e-10;typedef struct _point{double x,y;_point (double x = 0, double y = 0):x(x),y(y) {}}Point,Vector;Point p[maxn];typedef struct _circle{Point c;double r;}Circle;typedef struct _line{Point p;Vector v;}Line;Vector operator + (Vector a, Vector b) {return Vector(a.x + b.x, a.y + b.y);}Vector operator - (Vector a, Vector b) {return Vector(a.x - b.x, a.y - b.y);}double Dot (Vector a, Vector b) {return a.x * b.x + a.y * b.y;}double Length (Vector a) {return sqrt (Dot (a, a));}double Cross (Vector a, Vector b) {return a.x * b.y - a.y * b.x;}double Distance (Point a, Point b) {return Length (a - b);}double Angle (Vector a, Vector b) {return acos (Dot (a, b) / Length (a) / Length (b));}int dcmp (double x){if (fabs (x) < eps) return 0;else return x < 0 ? -1 : 1;}int getLineCircleIntersection (Line L, Circle C, double& t1, double& t2){//返回线段和圆的交点个数double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;double e = a * a + c * c, f = 2 * (a * b + c * d), g = b * b + d * d - C.r * C.r;double delta = f * f - 4 * e * g;//判别式if (dcmp (delta) < 0) return 0;if (dcmp (delta) == 0){t1 = t2 = -f / (2 * e);return 1;}t1 = (-f - sqrt (delta)) / (2 * e);t2 = (-f + sqrt (delta)) / (2 * e);return 2;}double slove (Circle r, int n){int i,j;double ans = 0;double x;for (i = 0; i < n; i++){if (dcmp(Distance(p[i], r.c) - r.r) <= 0 && dcmp(Distance(p[(i + 1) % n], r.c) - r.r) <= 0){//两个顶点都在圆内，直接求三角形面积ans += Cross (p[i] - r.c, p[(i + 1) % n] - r.c) / 2;}else if (dcmp(Distance(p[i], r.c) - r.r) <= 0 && dcmp(Distance(p[(i + 1) % n], r.c) - r.r) > 0){//p[i]在圆内，p[i+1]在圆外，求与圆的交点double t1,t2;Line L;L.p = p[i];L.v = p[(i + 1) % n] - p[i];getLineCircleIntersection (L, r, t1, t2);t1 = t1 > 0 ? t1 : t2;Point inter;inter.x = L.p.x + L.v.x * t1; inter.y = L.p.y + L.v.y * t1;ans += Cross (p[i] - r.c, inter - r.c) / 2;//三角形面积x = (Angle (inter - r.c, p[(i + 1) % n] - r.c) / 2.0) * r.r * r.r;//扇形面积if (Cross (inter - r.c, p[(i + 1) % n] - r.c) < 0)ans -= x;elseans += x;}else if (dcmp(Distance(p[i], r.c) - r.r) > 0 && dcmp(Distance(p[(i + 1) % n], r.c) - r.r) <= 0){//p[i]在圆外，p[i+1]在圆内double t1,t2;Line L;L.p = p[i];L.v = p[(i + 1) % n] - p[i];getLineCircleIntersection (L, r, t1, t2);t1 = t1 < t2 ? t1 : t2;Point inter;inter.x = L.p.x + L.v.x * t1; inter.y = L.p.y + L.v.y * t1;ans += Cross (inter - r.c, p[(i + 1) % n] - r.c) / 2;//三角形面积x = (Angle (p[i] - r.c, inter - r.c) / 2.0) * r.r * r.r;//扇形面积if (Cross (p[i] - r.c, inter - r.c) < 0)ans -= x;elseans += x;}else if (dcmp(Distance(p[i], r.c) - r.r) > 0 && dcmp(Distance(p[(i + 1) % n], r.c) - r.r) > 0){//两个点都在圆外double t1,t2;Line L;L.p = p[i];L.v = p[(i + 1) % n] - p[i];if (getLineCircleIntersection (L, r, t1, t2) <= 1){//只有一个或者没有交点x = (Angle (p[(i + 1) % n] - r.c, p[i] - r.c) / 2.0) * r.r * r.r;//扇形面积if (Cross (p[i] - r.c, p[(i + 1) % n] - r.c) < 0)ans -= x;elseans += x;}else{Point inter1,inter2;if (t1 > t2)swap (t1,t2);inter1.x = L.p.x + L.v.x * t1; inter1.y = L.p.y + L.v.y * t1;inter2.x = L.p.x + L.v.x * t2; inter2.y = L.p.y + L.v.y * t2;ans += Cross (inter1 - r.c, inter2 - r.c) / 2;//三角形面积x = (Angle (p[i] - r.c, inter1 - r.c) / 2.0) * r.r * r.r;//扇形面积if (Cross (p[i] - r.c, inter1 - r.c) < 0)ans -= x;elseans += x;x = (Angle (inter2 - r.c, p[(i + 1) % n] - r.c) / 2.0) * r.r * r.r;if (Cross (inter2 - r.c, p[(i + 1) % n] - r.c) < 0)ans -= x;elseans += x;}}}return ans;}int main (){double h;Circle c;while (scanf ("%lf%lf%lf",&c.c.x, &c.c.y, &h) != EOF){double x1,y1;scanf ("%lf%lf",&x1, &y1);//求圆心double t = sqrt (h / 5.0);c.c.x += x1 * t;c.c.y += y1 * t;scanf ("%lf",&c.r);int n;scanf ("%d",&n);int i,j;for (i = 0; i < n; i++){scanf ("%lf%lf",&p[i].x,&p[i].y);}printf ("%.2lf\n",fabs (slove (c, n)));}return 0;}/*5 5 100 0244 06 06 64 60 0 2000100 0100 41900 -1002000 -1002000 1001900 1000 0 2000100 0100 41900 1002000 1002000 -1001900 -1005 5 100 010041 02 02 11 155100 01 40 0100 0100 1000 100*/