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;
else
ans += 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;
else
ans += 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;
else
ans += 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;
else
ans += 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;
else
ans += 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 10
0 0
2
4
4 0
6 0
6 6
4 6
0 0 2000
100 0
100
4
1900 -100
2000 -100
2000 100
1900 100
0 0 2000
100 0
100
4
1900 100
2000 100
2000 -100
1900 -100
5 5 10
0 0
100
4
1 0
2 0
2 1
1 1
5
5
10
0 0
1 4
0 0
100 0
100 100
0 100
*/