分割线算图形面积 圆和多边形面积并

冬令营听了hwd的课，感觉好像计算几何很简单的样子。

最近做了一道计算圆和多边形面积并的面积，用了扫面线的方法，感觉不太好（可能我没写习惯），写了(7k)的代码才A。QAQ

以后记得找月下柠檬树做做。

代码

#include <cstdio>
#include <iostream>
#include <vector>
#include <iomanip>
#include <cmath>
#include <assert.h>
using namespace std;

const int MAXN = 103;
const long double PI = acos(-1);
const long double EPS = 1e-15;

template <class T>
T sqr(const T &x) {
	return x * x;
}

int n, m;

struct Point {
	long double x, y;

	Point () {}

	Point (long double a, long double b):x(a), y(b) {}

	Point operator - (const Point &o) const {
		return Point(x - o.x, y - o.y);
	}

	Point operator + (const Point &o) const {
		return Point(x + o.x, y + o.y);
	}
	
	long double operator * (const Point &o) const {
		return x * o.y - y * o.x;
	}

	long double operator % (const Point &o) const {
		return x * o.x + y * o.y;
	}

	Point operator * (const long double &o) const {
		return Point(x * o, y * o);
	}

	Point operator / (const long double &o) const {
		return Point(x / o, y / o);
	}
};

struct Segment {
	Point first, second;

	Segment() {}

	Segment(Point x, Point y):first(x), second(y) {}
};

struct Circle {
	Point center;
	long double r;
};

Circle A[MAXN];
Point B[MAXN];

Point ret0, ret1;

long double dis(const Point &a, const Point &b) {
	return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));
}

Point rotate(const Point &p, const long double &c, const long double &s) {
	return Point(p.x * c - p.y * s, p.x * s + p.y * c);
}

bool cir_cir(const Circle &a, const Circle &b) {
	long double len = dis(a.center, b.center);
	if (a.r + b.r < len) return false;
	if (len < fabs(a.r - b.r)) return false;
	long double COS = (sqr(a.r) + sqr(len) - sqr(b.r)) / (len * a.r * 2);
	long double SIN = sqrt(1. - sqr(COS));
	Point tmp = b.center - a.center;
	long double k = a.r / len;
	ret0 = a.center + rotate(tmp, COS, SIN) * k;
	ret1 = a.center + rotate(tmp, COS, -SIN) * k;
	return true;
}

Point line_line(const Segment &a, const Segment &b) {
	long double k = (a.second - a.first) * (b.first - a.first);
	k = k / (k - (a.second - a.first) * (b.second - a.first));
	return b.first + (b.second - b.first) * k;
}

bool cir_line(const Circle &a, const Segment &b) {
	Point v = b.second - b.first;
	swap(v.x, v.y);
	v.x = -v.x;
	v = a.center + v;
	Point A = line_line(Segment(v, a.center), b);
	long double d = dis(A - a.center, Point(0, 0));
	if (d > a.r) return false;
	long double len = sqrt(sqr(a.r) - sqr(d));
	long double norm = dis(b.second - b.first, Point(0, 0));
	ret0 = A + (b.second - b.first) * len / norm;
	ret1 = A - (b.second - b.first) * len / norm;
	return true;
}

long double cirSubArea(const Circle &a, const Point &p, const Point &q) {
	long double X = dis(p, q);
	long double angle = acos(1. - sqr(X) * .5 / sqr(a.r));
		//acos((sqr(a.r) * 2 - sqr(X)) * .5 / sqr(a.r));
	long double ret = angle * sqr(a.r) * .5;
	ret -= fabs((p - a.center) * (q - a.center) * .5);
	return ret;
}

struct Query {
	int belong;
	bool enter;
	Point left, right;
	long double key;

	bool operator < (const Query &o) const {
		return key < o.key;
	}
};

int main() {
	freopen("polygon.in", "r", stdin);
	freopen("polygon.out", "w", stdout);

	cin >> n >> m;
	vector<long double> key;

	for (int i = 0; i < n; i ++) {
		cin >> A[i].center.x >> A[i].center.y >> A[i].r;
		key.push_back(A[i].center.x - A[i].r);
		key.push_back(A[i].center.x + A[i].r);
		for (int j = 0; j < i; j ++) {
			if (cir_cir(A[i], A[j])) {
				key.push_back(ret0.x);
				key.push_back(ret1.x);
			}
		}
	}
	for (int i = 0; i < m; i ++) {
		cin >> B[i].x >> B[i].y;
		key.push_back(B[i].x);
	}
	B[m] = B[0];
	for (int i = 0; i < n; i ++) {
		for (int j = 0; j < m; j ++) {
			if (cir_line(A[i], Segment(B[j], B[j + 1]))) {
				if (ret0.x >= min(B[j].x, B[j + 1].x) - EPS &&
						ret0.x <= max(B[j].x, B[j + 1].x) + EPS)
					key.push_back(ret0.x);
				if (ret1.x >= min(B[j].x, B[j + 1].x) - EPS &&
						ret1.x <= max(B[j].x, B[j + 1].x) + EPS)
					key.push_back(ret1.x);
			}
		}
	}

	sort(key.begin(), key.end());
	vector<Query> que;

	long double ans = 0;

	for (int i = 0, _i = key.size(); i + 1 < _i; i ++) {
		if (i < _i && key[i] == key[i + 1]) continue;
		long double mid = (key[i] + key[i + 1]) * .5;
		Segment vecl = Segment(Point(key[i], 0), Point(key[i], 1));
		Segment vecr = Segment(Point(key[i + 1], 0), Point(key[i + 1], 1));
		Segment vecm = Segment(Point(mid, 0), Point(mid, 1));
		que.clear();

		for (int j = 0; j < n; j ++) {
			if (cir_line(A[j], vecm)) {
				if (ret0.y > ret1.y) swap(ret0, ret1);
				long double key0 = ret0.y, key1 = ret1.y;

				cir_line(A[j], vecl);
				Point a = ret0, b = ret1;
				if (a.y > b.y) swap(a, b);
				cir_line(A[j], vecr);
				if (ret0.y > ret1.y) swap(ret0, ret1);
				Query tmp;
				tmp.belong = j;

				tmp.left = a;
				tmp.right = ret0;
				tmp.enter = true;
				tmp.key = key0;
				que.push_back(tmp);
				tmp.left = b;
				tmp.right = ret1;
				tmp.enter = false;
				tmp.key = key1;
				que.push_back(tmp);
			}
		}
		bool in = false;
		for (int j = 0; j < m; j ++) {
			Point t = line_line(Segment(B[j], B[j + 1]), vecm);
			if (((B[j + 1] - B[j]) % (t - B[j]) > 0) == 
					((B[j + 1] - B[j]) % (t - B[j + 1]) > 0)) continue;

			in = true;
			Query tmp;
			tmp.belong = -1;
			tmp.left = line_line(vecl, Segment(B[j], B[j + 1]));
			tmp.right = line_line(vecr, Segment(B[j], B[j + 1]));
			tmp.enter = false;
			tmp.key = t.y;
			que.push_back(tmp);
		}
		if (in) {
			int t = que.size();
			if (que[t - 2].left.y >= que[t - 1].left.y &&
					que[t - 2].right.y >= que[t - 1].right.y)
				que[t - 1].enter = true;
			else
				que[t - 2].enter = true;
		}

		sort(que.begin(), que.end());

		int cover = 0;
		Query lower;
		for (int j = 0, _j = que.size(); j < _j; j ++) {
			Query &cur = que[j];
			if (cur.enter) {
				if (cover == 0) {
					lower = cur;
					if (cur.belong != -1)
						ans += cirSubArea(A[cur.belong], 
								cur.left, cur.right);
				}
				cover ++;
			}
			else {
				cover --;
				if (cover == 0) {
					ans += (cur.left.y + cur.right.y - lower.left.y - lower.right.y) *
						(key[i + 1] - key[i]) * .5;
					if (cur.belong != -1) 
						ans += cirSubArea(A[cur.belong], cur.left, cur.right);
				}
			}
		}
	}

	//cout << ans << endl;
	printf("%.15lf
", (double) ans);

	return 0;
}