convex_hull.cpp

ACM凸包标程,可以快速使用凸包的各种应用

#include <iostream>
#include <cmath>

using namespace std;

#define EPS 1e-10

#define MAXN 1000

struct point
{
	double x, y;
	int id;
	point() {}
	point(double nx, double ny) : x(nx), y(ny) {}
	point initi() const
	{
		double l = hypot(x, y);
		return point(x / l, y / l);
	}
	point rotate(double ang) const
	{
		return point(x * cos(ang) - y * sin(ang), x * sin(ang) + y * cos(ang));
	}
};

point p[MAXN]; int pnum;

point h[MAXN]; int chn;

int bt[MAXN], tp[MAXN];
int bt_top, tp_top;

struct circle
{
	double x, y;
	double r;
	int id;
};

circle c[MAXN]; int cnum;
circle t[MAXN];

void process()
{
	int i, j, k = 0;
	for(i = 0; i < cnum; i++)
	{
		for(j = 0; j < cnum; j++)
		{
			if(i == j) continue;
			double dist = hypot(c[i].x - c[j].x, c[i].y - c[j].y);
			if(c[j].r >= c[i].r && c[j].r - c[i].r >= dist) break;
		}
		if(j == cnum && c[i].r > EPS) t[k++] = c[i];
	}
	//printf("%d\n", k);
	cnum = k;
	
	for(i = 0; i < cnum; i++)
	{
		c[i] = t[i];
		c[i].id = i;
		//printf("%lf %lf %lf\n", c[i].x, c[i].y, c[i].r);
	}
	//putchar('\n');
}

void add_point()
{
	pnum = 0;
	int i, j;
	for(i = 0; i < cnum; i++)
	{
		for(j = i + 1; j < cnum; j++)
		{
			circle c1 = c[i], c2 = c[j];
			if(c1.r < c2.r) swap(c1, c2);
			double d = hypot(c1.x - c2.x, c1.y - c2.y);
			double h = c1.r - c2.r;
			double theta = fabs(acos(h / d));
			
			point p1, p2;

			point vt, vl, vr;
			vt.x = c2.x - c1.x, vt.y = c2.y - c1.y;

			vt = vt.initi();

			vl = vt.rotate(theta), vr = vt.rotate(-theta);

			p1.x = c1.x + c1.r * vl.x;
			p1.y = c1.y + c1.r * vl.y;
			p1.id = c1.id;

			p2.x = c2.x + c2.r * vl.x;
			p2.y = c2.y + c2.r * vl.y;
			p2.id = c2.id;

			p[pnum++] = p1;
			p[pnum++] = p2;

			p1.x = c1.x + c1.r * vr.x;
			p1.y = c1.y + c1.r * vr.y;
			p1.id = c1.id;

			p2.x = c2.x + c2.r * vr.x;
			p2.y = c2.y + c2.r * vr.y;
			p2.id = c2.id;

			p[pnum++] = p1;
			p[pnum++] = p2;
		}
	}
	/*
	for(i = 0; i < pnum; i++)
	{
		printf("%lf %lf %d\n", p[i].x, p[i].y, p[i].id);
	}
	putchar('\n');
	*/
}

bool cmp(const point &p1, const point &p2)
{
	return (p1.x < p2.x) || (fabs(p1.x - p2.x) < EPS && p1.y < p2.y);
}

int dcmp(double d)
{
	if(fabs(d) < EPS) return 0;
	return d > 0 ? 1 : -1;
}

int det(double x1, double y1, double x2, double y2)
{
	return dcmp(x2 * y1 - x1 * y2);
}

int cross(point a, point b, point c)
{
	return det(b.x - a.x, b.y - a.y, c.x - b.x, c.y - b.y);
}

void convex_hull()
{
	int i;
	sort(p, p + pnum, cmp);
	bt[0] = 0, tp[0] = 0, bt[1] = 1, tp[1] = 1;
	bt_top = 1, tp_top = 1;
	for(i = 0; i < pnum; i++)
	{
		if(cross(p[bt[1]], p[0], p[i]) <= 0) bt[1] = i;
		if(cross(p[tp[1]], p[0], p[i]) >= 0) tp[1] = i;
	}
	for(i = bt[bt_top] + 1; i < pnum; i++)
	{
		while(1)
		{
			if(cross(p[bt[bt_top - 1]], p[bt[bt_top]], p[i]) >= 0) bt_top--;
			else
			{
				bt[++bt_top] = i;
				break;
			}
		}
	}
	for(i = tp[tp_top] + 1; i < pnum; i++)
	{
		while(1)
		{
			if(cross(p[tp[tp_top - 1]], p[tp[tp_top]], p[i]) <= 0) tp_top--;
			else
			{
				tp[++tp_top] = i;
				break;
			}
		}
	}
	chn = 0;
	for(i = 0; i <= bt_top; i++) h[chn++] = p[bt[i]];
	for(i = tp_top - 1; i >= 1; i--) h[chn++] = p[tp[i]];
}

void calc_total_area()
{
	h[chn] = h[0];
	double area = 0;
	int i;
	/*
	for(i = 0; i < chn; i++)
	{
		printf("%lf %lf %d\n", h[i].x, h[i].y, h[i].id);
	}
	putchar('\n');
	*/
	for(i = 0; i < chn; i++)
	{
		area += h[i].x * h[i + 1].y - h[i].y * h[i + 1].x;
	}
	area = fabs(area / 2.0);
	//printf("area = %lf\n", area);
	for(i = 0; i < chn; i++)
	{
		if(h[i].id == h[i + 1].id)
		{
			int d = h[i].id;
			point ct = point(c[d].x, c[d].y);
			int f = cross(h[i], h[i + 1], ct);
			
			double dist = hypot(h[i].x - h[i + 1].x, h[i].y - h[i + 1].y);
			double theta = fabs(2 * asin(dist / (2 * c[d].r)));

			//printf("%lf\n", theta);

			double arc_area = c[d].r * c[d].r * (theta - sin(theta)) * .5;

			if(f > 0)
			{
				//printf("YES\n");
				arc_area = M_PI * c[d].r * c[d].r - arc_area;
			}
			area += arc_area;
		}
	}
	printf("%lf\n", area);
}

int main(void)
{
	int i;
	while(scanf("%d", &cnum) != EOF)
	{
		for(i = 0; i < cnum; i++)
		{
			scanf("%lf %lf %lf", &c[i].x, &c[i].y, &c[i].r);
		}
		process();//去掉包含和半径为0的圆
		if(cnum == 1)
		{
			printf("%lf\n", M_PI * c[0].r * c[0].r);
			continue;
		}
		add_point();
		convex_hull();
		calc_total_area();
	}
	return 0;
}