解题报告-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;
			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
*/