【Luogu P4406】「CQOI2005」三角形面积并

Description

给出 (n) 个三角形，求这 (n) 个三角形并的面积。

数据范围：(1 leq n leq 100)，(-10^6 leq x, y leq 10^6)。

Solution

扫描线的做法大家都会，这篇题解的做法是自适应辛普森积分。

记 (f(t)) 表示：给出的 (n) 个三角形与直线 (x = t) 的交集长度。

那么答案即为：

[int_a^b f(x) ext{dx} ]

其中 (a, b) 分别表示输入的 (x) 中的最小值和最大值。

那么可以使用自适应辛普森积分做。现在的关键是要如何求出任意一个 (f(t))。
对于每一个三角形，考虑求出这个三角形与直线 (x = t) 的相交部分是哪一段区间，把这些区间提出来，区间合并即可求出函数值。计算一个函数值的时间复杂度是 (mathcal{O}(n log n)) 的，可以接受。

然后需要提高精度：

• 可以在递归过程中设置一个强制执行的最少迭代层数。
• 可以将每个三角形的最小 (x) 值和最大 (x) 值提出来后排序，对所有相邻点构成的区间求积分，将这些积分相加，就可以得到最后大区间 ([a, b]) 的积分。

Code

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <vector>

using namespace std;

const double eps = 1e-9;
const double INF = 1e6;

int sgn(double n) {
	if (fabs(n) < eps) return 0;
	return n > 0 ? 1 : -1; 
}

int dcmp(double x, double y) {
	return sgn(x - y);
}

struct point {
	double x, y;

	point() { x = y = 0; }
	point(double A, double B) : x(A), y(B) {}
};

typedef point vec;

vec operator + (vec a, vec b) { return vec(a.x + b.x, a.y + b.y); }
vec operator - (vec a, vec b) { return vec(a.x - b.x, a.y - b.y); }
vec operator * (vec a, double b) { return vec(a.x * b, a.y * b); }
vec operator / (vec a, double b) { return vec(a.x / b, a.y / b); }

bool operator < (point a, point b) {
	if (dcmp(a.x, b.x)) return a.x < b.x;
	return a.y < b.y;
}

double cross(vec a, vec b) {
	return a.x * b.y - a.y * b.x;
}

point line_intersection(point A, point B, point C, point D) {
	point p = A, q = C;
	vec x = B - A, y = D - C;

	vec u = q - p;
	double t = cross(u, y) / cross(x, y);
	return p + x * t;
}

const int N = 110;

int n;
vector<point> a[N];

int t;
double pos[N * 2];

int m;
struct range {
	double l, r;

	range() {}
	range(double A, double B) : l(A), r(B) {}
} b[N];

bool ruler(range a, range b) {
	if (dcmp(a.l, b.l)) return a.l < b.l;
	return a.r < b.r;
}

double f(double x) {
	point Ga = point(x, 0), Gb = point(x, 1);

	m = 0;
	for (int i = 1; i <= n; i ++) {
		vector<point> u = a[i];

		if (dcmp(x, u[0].x) < 0) continue;
		if (dcmp(x, u[2].x) > 0) continue;

		if (!dcmp(u[0].x, u[1].x) && !dcmp(u[0].x, x)) {
			b[++ m] = range(u[0].y, u[1].y);
		} else if (!dcmp(u[1].x, u[2].x) && !dcmp(u[1].x, x)) {
			b[++ m] = range(u[1].y, u[2].y);
		} else {
			vector<double> seq;
			seq.clear();

			for (int j = 0; j < 3; j ++) {
				point A = u[j], B = u[(j + 1) % 3];
				if (B < A) swap(A, B);

				if (dcmp(x, A.x) < 0) continue;
				if (dcmp(x, B.x) > 0) continue;

				seq.push_back(line_intersection(A, B, Ga, Gb).y);
			}

			sort(seq.begin(), seq.end());
			b[++ m] = range(seq[0], seq[seq.size() - 1]);
		}
	}

	if (!m) return 0;

	sort(b + 1, b + 1 + m, ruler);

	double ans = 0;
	double st = b[1].l, ed = b[1].r;

	for (int i = 2; i <= m; i ++) {
		if (dcmp(b[i].l, ed) > 0)
			ans += ed - st,
			st = b[i].l, ed = b[i].r;
		else
			ed = max(ed, b[i].r);
	}

	ans += ed - st;

	return ans;
}

double simpson(double Lv, double Mv, double Rv, double len) {
	return (Lv + Rv + 4 * Mv) * len / 6; 
}

double asr(double l, double r, double Lv, double Mv, double Rv, int dep) {
	double mid = (l + r) / 2;
	double A = f((l + mid) / 2), B = f((mid + r) / 2);

	double s = simpson(Lv, Mv, Rv, r - l);
	double Ls = simpson(Lv, A, Mv, mid - l);
	double Rs = simpson(Mv, B, Rv, r - mid);

	if (!dcmp(s, Ls + Rs) && dep <= 0) return Ls + Rs;
	return asr(l, mid, Lv, A, Mv, dep - 1) + asr(mid, r, Mv, B, Rv, dep - 1);
}

int main() {
	scanf("%d", &n);

	for (int i = 1; i <= n; i ++) {
		for (int j = 0; j < 3; j ++) {
			double x, y;
			scanf("%lf%lf", &x, &y);

			a[i].push_back(point(x, y));
		}

		sort(a[i].begin(), a[i].end());

		pos[++ t] = a[i][0].x;
		pos[++ t] = a[i][2].x;
	}

	sort(pos + 1, pos + 1 + t);

	double ans = 0;

	for (int i = 2; i <= t; i ++) {
		if (!dcmp(pos[i - 1], pos[i])) continue;

		double l = pos[i - 1], r = pos[i];
		ans += asr(l, r, f(l), f((l + r) / 2), f(r), 5);
	}

	printf("%.2lf
", ans);

	return 0;
}