算法模板（7）：计算几何（2）

zhezhidashi

已于 2023-06-14 13:31:51 修改

阅读量729

点赞数

分类专栏：算法模板文章标签：算法图论 c++

于 2023-06-12 09:42:45 首次发布

本文链接：https://blog.csdn.net/qq_45812711/article/details/131153784

版权

算法模板专栏收录该内容

31 篇文章 0 订阅

订阅专栏

计算几何

旋转卡壳

其实旋转卡壳就是枚举每一条凸包上的边，然后找一找离这条边最远的点是哪个。确定一条边，找那个点用的就是双指针法。因此，旋转卡壳大多是枚举的边，在边上找特征。

2938. 周游世界

题意：给定一个二维平面，平面上有 N 个点。每个点的位置可由一对整数坐标 (x,y) 来表示，不同的点位置不同。请你求出平面上距离最远的点对之间的距离是多少。

#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;

#define x first
#define y second

typedef pair<int, int> P;
const int maxn = 50010;
int N;
P q[maxn];
int stk[maxn], top;
bool used[maxn];

P operator-(P a, P b) {
	return { a.x - b.x, a.y - b.y };
}

int operator*(P a, P b) {
	return a.x * b.y - a.y * b.x;
}

int area(P a, P b, P c) {
	return (b - a) * (c - a);
}

int get_dist(P a, P b) {
	int dx = a.x - b.x;
	int dy = a.y - b.y;
	return dx * dx + dy * dy;
}

void get_convex() {
	sort(q, q + N);
	for (int i = 0; i < N; i++) {
		while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0) {
			if (area(q[stk[top - 2]], q[stk[top - 1]], q[i]) < 0) {
				used[stk[--top]] = false;
			}
			else top--;
		}
		stk[top++] = i;
		used[i] = true;
	}
	used[0] = false;
	for (int i = N - 1; i >= 0; i--) {
		if (used[i]) continue;
		while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0)
			top--;
		stk[top++] = i;
	}
	//因为0号点在起点和终点都插入了一次，因此删掉一个。
	top--;
}

int rotating_calipers() {
	if (top <= 2) return get_dist(q[0], q[N - 1]);

	int res = 0;
	for (int i = 0, j = 2; i < top; i++) {
		//这里 i + 1 不用取模，因为根据上面那个凸包的函数，此时的 stk[top] 就是第一个点的下标
		auto d = q[stk[i]], e = q[stk[i + 1]];
		while (area(d, e, q[stk[j]]) < area(d, e, q[stk[j + 1]])) j = (j + 1) % top;
		res = max(res, max(get_dist(q[stk[j]], d), get_dist(q[stk[j]], e)));
	}
	return res;
}

int main() {
	scanf("%d", &N);
	for (int i = 0; i < N; i++) scanf("%d%d", &q[i].x, &q[i].y);
	get_convex();
	printf("%d\n", rotating_calipers());

	return 0;
}

2142. 最小矩形覆盖

已知平面上不共线的一组点的坐标，求覆盖这组点的面积最小的矩形。输出矩形的面积和四个顶点的坐标。
上面那个证明出来，矩形至少有一条边和凸包的一条边重合。那么，我们其实就是枚举凸包的每一条边，然后找到离这条边最远的边，即投影最大的。然后，找两边的点时，就是找投影最小的。这样矩形的三个点就找到了，最后一个点可以根据矩形的边的关系算出来。

#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
using namespace std;

#define x first
#define y second

typedef pair<double, double> P;
const int maxn = 50010;
const double eps = 1e-12, PI = acos(-1), INF = 1e20;

int N;
P q[maxn], ans[10];
double min_area = INF;

int stk[maxn], top;
bool used[maxn];

int sign(double x) {
	if (fabs(x) < eps) return 0;
	if (x < 0) return -1;
	return 1;
}
int dcmp(double x, double y) {
	if (fabs(x - y) < eps) return 0;
	if (x < y) return -1;
	return 1;
}
P operator+(P a, P b) {
	return { a.x + b.x, a.y + b.y };
}
P operator-(P a, P b) {
	return{ a.x - b.x, a.y - b.y };
}
double operator*(P a, P b) {
	return a.x * b.y - a.y * b.x;
}
P operator*(P a, double t) {
	return { a.x * t, a.y * t };
}
P operator/(P a, double t) {
	return { a.x / t, a.y / t };
}
double operator&(P a, P b) {
	return a.x * b.x + a.y * b.y;
}
double area(P a, P b, P c) {
	return (b - a) * (c - a);
}
double get_len(P a) {
	return sqrt(a & a);
}
double project(P a, P b, P c) {
	//AC 在 AB 上面的投影。
	return ((b - a) & (c - a)) / get_len(b - a);
}

P unit(P a) {
	//单位向量
	return a / get_len(a);
}
P rotate(P a, double b) {
	return { a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b) };
}

void get_convex() {
	sort(q, q + N);
	for (int i = 0; i < N; i++) {
		while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0) {
			if (area(q[stk[top - 2]], q[stk[top - 1]], q[i]) < 0)
				used[stk[--top]] = false;
			else top--;

		}
		stk[top++] = i;
		used[i] = true;
	}
	used[0] = false;
	for (int i = N - 1; i >= 0; i--) {
		if (used[i]) continue;
		while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0)
			top--;
		stk[top++] = i;
	}
	top--;
}

void rotating_calipers() {
	//由题意知不用判断是否共线
	for (int i = 0, a = 2, b = 2, c = 2; i < top; i++) {
		auto d = q[stk[i]], e = q[stk[i + 1]];
		while (dcmp(area(d, e, q[stk[a]]), area(d, e, q[stk[a + 1]])) < 0)
			a = (a + 1) % top;
		while (dcmp(project(d, e, q[stk[b]]), project(d, e, q[stk[b + 1]])) < 0)
			b = (b + 1) % top;
		//c这个点一开始必须从a开始往左走，这样才能保证c是在往极小值点移动。
		if (!i) c = a;
		//注意这里是 > 0 了，因为在左边投影是负的。
		//一定要小心这个 > 0 写在括号外面，不然会陷入死循环。
		while (dcmp(project(d, e, q[stk[c]]), project(d, e, q[stk[c + 1]])) > 0) 
			c = (c + 1) % top;
			
		
		auto x = q[stk[a]], y = q[stk[b]], z = q[stk[c]];
		//矩形的高和宽
		auto h = area(d, e, x) / get_len(e - d);
		auto w = ((y - z) & (e - d)) / get_len(e - d);
		if (h * w < min_area) {
			min_area = h * w;
			ans[0] = d + unit(e - d) * project(d, e, y);
			ans[3] = d + unit(e - d) * project(d, e, z);

			auto u = unit(rotate(e - d, -PI / 2));

			ans[1] = ans[0] + u * h;
			ans[2] = ans[3] + u * h;
		}
	}
}

int main() {
	scanf("%d", &N);
	for (int i = 0; i < N; i++) scanf("%lf%lf", &q[i].x, &q[i].y);
	get_convex();
	rotating_calipers();

	printf("%.5f\n", min_area);

	int k = 0;
	for (int i = 1; i < 4; i++) {
		if (dcmp(ans[i].y, ans[k].y) < 0 ||
			!dcmp(ans[i].y, ans[k].y) && dcmp(ans[i].x, ans[k].x) < 0)
			k = i;
	}
	for (int i = 0; i < 4; i++, k++) {
		auto x = ans[k % 4].x, y = ans[k % 4].y;
		//直接输出的话，会有-0.00000
		if (!sign(x)) x = 0;
		if (!sign(y)) y = 0;
		printf("%.5f %.5f\n", x, y);
	}
	return 0;
}

三角剖分

3034. 望远镜

一个圆，其圆心位于原点，半径为 R。一个 N 个顶点的简单多边形。求这两个图形相交的面积。
方法：就是按照逆时针求多边形有向面积的方法，只不过这次求的是三角形与圆相交的有向面积，把这些加起来就是答案。

讨论三角形与圆的交点

#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>

#define x first
#define y second

using namespace std;

typedef pair<double, double> P;
const int maxn = 55;
const double eps = 1e-8;
//千万别写成 cos(-1)
const double PI = acos(-1);

double R;
int N;
P q[maxn], o;

int sign(double x) {
	if (fabs(x) < eps) return 0;
	if (x < 0) return -1;
	return 1;
}

//小心，这里 < eps 千万别写成 < 0.
int dcmp(double x, double y) {
	if (fabs(x - y) < eps) return 0;
	if (x < y) return -1;
	return 1;
}

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

double area(P a, P b, P c) {
	return (b - a) * (c - a);
}
double get_len(P a) {
	return sqrt(a & a);
}
double get_dist(P a, P b) {
	return get_len(b - a);
}
double project(P a, P b, P c) {
//AC 在 AB 上面的投影
	return ((b - a) & (c - a)) / get_len(b - a);
}
P rotate(P a, double b){
	return { a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b) };
}
P unit(P a) {
	return a / get_len(a);
}
bool on_segment(P p, P a, P b) {
	return !sign((p - a) * (p - b)) && ((p - a) & (p - b)) <= 0;
}
P get_line_intersection(P p, P v, P q, P w) {
	P u = p - q;
	double t = (w * u) / (v * w);
	return p + v * t;
}

double get_sector(P a, P b) {
	//a 和 b 形成的扇形的面积，很容易看出圆心角一定小于 PI.
	auto angle = acos((a & b) / get_len(a) / get_len(b));
	//如果 oab 形成有向面积是负的，那么返回的扇形的面积也要是负的。
	if (sign(a * b) < 0) angle = -angle;
	return R * R * angle / 2;
}

double get_circle_line_intersection(P a, P b, P& pa, P& pb) {
	auto e = get_line_intersection(a, b - a, o, rotate(b - a, PI / 2));
	//mind 表示 o 到 线段 ab 的距离
	auto mind = get_dist(o, e);
	//判断 e 是否在 ab上面
	if (!on_segment(e, a, b)) mind = min(get_dist(o, a), get_dist(o, b));
	//线段ab与圆无交点
	if (dcmp(R, mind) <= 0) return mind;
	//求垂足到线段与圆交点的距离
	auto len = sqrt(R * R - get_dist(o, e) * get_dist(o, e));
	pa = e + unit(a - b) * len;
	pb = e + unit(b - a) * len;
	return mind;
}

double get_circle_triangle_area(P a, P b) {
	//一定要注意，是按照从 a 到 b 逆时针的顺序计算叉积和扇形面积。
	//a 和 b 共线
	if (!sign(a * b)) return 0;

	auto da = get_dist(o, a), db = get_dist(o, b);
	//第一种情况
	if (dcmp(da, R) <= 0 && dcmp(db, R) <= 0) return a * b / 2;
	P pa, pb;  //与圆的交点
	//mind 表示 o 到线段 ab 的距离
	auto mind = get_circle_line_intersection(a, b, pa, pb);
	
	//第二种情况
	if (dcmp(R, mind) <= 0) return get_sector(a, b);

	//第三种情况
	if (dcmp(da, R) <= 0) return a * pb / 2 + get_sector(pb, b);
	//第四种情况
	if (dcmp(db, R) <= 0) return pa * b / 2 + get_sector(a, pa);
	//第五种情况
	return get_sector(a, pa) + get_sector(pb, b) + pa * pb / 2;
}

double work() {
	//题目没说是顺时针还是逆时针。
	//如果给的是逆时针，那么求得面积是正的。
	//如果给的是顺时针，那么求的面积是负的。
	double res = 0;
	for (int i = 0; i < N; i++) {
		res += get_circle_triangle_area(q[i], q[(i + 1) % N]);
	}
	return fabs(res);
}

int main() {
	while (cin >> R >> N) {
		for (int i = 0; i < N; i++) scanf("%lf%lf", &q[i].x, &q[i].y);
		printf("%.2f\n", work());
	}
	return 0;
}

扫描线

主要是用来求面积或者周长

3068. 扫描线

在二维平面中给定 n 个两条边分别与 x 轴和 y 轴平行的矩形，请你求出它们的面积并。
先把举行切割成一片一片的，再把每一片的矩形做一个区间和并（排序），复杂度 $O(n^2logn)$ .

#include<bits/stdc+
using namespace std;

#define x first
#define y second

const int maxn = 1010;
typedef long long ll;
typedef pair<ll, ll> P;
P l[maxn], r[maxn];
vector<ll> xs;
int N;
P q[maxn];  //存储高度的区间

ll range_area(ll a, ll b) {
	int cnt = 0;
	for (int i = 0; i < N; i++) {
		if (l[i].x <= a && b <= r[i].x) {
			q[cnt++] = { l[i].y, r[i].y };
		}
	}
	
	if (!cnt) return 0;
	sort(q, q + cnt);
	ll st = q[0].first, ed = q[0].second;
	ll res = 0;
	for (int i = 1; i < cnt; i++) {
		if (q[i].first > ed) {
			res += ed - st;
			st = q[i].first, ed = q[i].second;
		}
		else ed = max(ed, q[i].second);
	}
	res += ed - st;
	return res * (b - a);
}

int main() {
	scanf("%d", &N);
	
	for (int i = 0; i < N; i++) {
		scanf("%lld%lld%lld%lld", &l[i].x, &l[i].y, &r[i].x, &r[i].y);
		xs.push_back(l[i].x), xs.push_back(r[i].x);
	}
	sort(xs.begin(), xs.end());
	ll ans = 0;
	for (int i = 0; i + 1 < xs.size(); i++) {
		if (xs[i] != xs[i + 1]) {
			ans += range_area(xs[i], xs[i + 1]);
		}
	}
	printf("%lld\n", ans);
	return 0;
}

2801. 三角形面积并

给出 n 个三角形，求它们并的面积。
$O(n^3logn)$ .

#include<iostream>
#include<cstring>
#include<algorithm>
#include<vector>
#include<cmath>
using namespace std;

#define x first
#define y second

typedef pair<double, double> P;
const int maxn = 110;
const double eps = 1e-8, INF = 1e6;
int N;
P tr[maxn][3];
P q[maxn];
vector<double> xs;

int sign(double x) {
	if (fabs(x) < eps) return 0;
	if (x < 0) return -1;
	return 1;
}

int dcmp(double x, double y) {
	if (fabs(x - y) < eps) return 0;
	if (x < y) return -1;
	return 1;
}

P operator+(P a, P b) {
	return { a.x + b.x, a.y + b.y };
}

P operator-(P a, P b) {
	return { a.x - b.x, a.y - b.y };
}

P operator*(P a, double t) {
	return { a.x * t, a.y * t };
}

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

double operator&(P a, P b) {
	return a.x * b.x + a.y * b.y;
}
//这里的 on_segment 需要和下面的 get_line_intersection，才能求线段交点。
bool on_segment(P p, P a, P b) {
	return sign((p - a) & (p - b)) <= 0;
}

P get_line_intersection(P p, P v, P q, P w) {
	//注: p, q 是线段起点，v, w 是方向向量。要是想得到线段终点需要 p + v, q + w.
	//两条线段平行。
	if (!sign(v * w)) return { INF, INF };
	auto u = p - q;
	auto t = w * u / (v * w);
	auto o = p + v * t;
	//需要判断 o 是否在两条线段上。
	if (!on_segment(o, p, p + v) || !on_segment(o, q, q + w)) {
		return { INF, INF };
	}
	return o;
}

double line_area(double a, int side) {
	int cnt = 0;
	for (int i = 0; i < N; i++) {
		auto t = tr[i];

		//只有 a 在 t[0].x 到 t[2].x 之间，三角形才与 x = a 有交点。 
		if (dcmp(t[0].x, a) > 0 || dcmp(t[2].x, a) < 0) continue;
		
		//需要特判三角形的一条边与 x = a 重合
		if (!dcmp(t[0].x, a) && !dcmp(t[1].x, a)) {
			if (side) q[cnt++] = { t[0].y, t[1].y };
		}
		else if (!dcmp(t[1].x, a) && !dcmp(t[2].x, a)) {
			if (!side) q[cnt++] = { t[1].y, t[2].y };
		}
		else {
			//一般情况，三角形与 x = a 要么是两个交点，要么是三个交点（顶点在x = a上）
			double d[3];
			int u = 0;
			for (int j = 0; j < 3; j++) {
				auto o = get_line_intersection(t[j], t[(j + 1) % 3] - t[j], { a, -INF }, { 0, INF * 2 });
				if(dcmp(INF, o.x)) d[u++] = o.y;
			}
			if (u) {
				sort(d, d + u);
				q[cnt++] = { d[0], d[u - 1] };
			}
		}
	}
	if (!cnt) return 0;
	for (int i = 0; i < cnt; i++) {
		if (q[i].first > q[i].second)
			swap(q[i].first, q[i].second);
	}
	sort(q, q + cnt);
	double res = 0, st = q[0].x, ed = q[0].y;
	for (int i = 1; i < cnt; i++) {
		if (q[i].x <= ed) ed = max(q[i].y, ed);
		else {
			res += ed - st;
			st = q[i].x, ed = q[i].y;
		}
	}
	res += ed - st;
	return res;
}

double range_area(double a, double b) {
	return (line_area(a, 1) + line_area(b, 0)) * (b - a) / 2;
}

int main() {
	scanf("%d", &N);
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < 3; j++) {
			scanf("%lf%lf", &tr[i][j].x, &tr[i][j].y);
			xs.push_back(tr[i][j].x);
		}
		sort(tr[i], tr[i] + 3);
	}
	for (int i = 0; i < N; i++) {
		for (int j = i + 1; j < N; j++) {
			for (int x = 0; x < 3; x++) {
				for (int y = 0; y < 3; y++) {
					auto o = get_line_intersection(tr[i][x], tr[i][(x + 1) % 3] - tr[i][x],
						tr[j][y], tr[j][(y + 1) % 3] - tr[j][y]);
					if(dcmp(INF, o.x)) xs.push_back(o.x);
				}
			}
		}
	}
	sort(xs.begin(), xs.end());
	double res = 0;
	for (int i = 0; i + 1 < xs.size(); i++) {
		if (dcmp(xs[i], xs[i + 1])) {
			res += range_area(xs[i], xs[i + 1]);
		}
	}
	printf("%.2f\n", res);
	return 0;
}

自适应辛普森积分

当图形过于复杂时，可以用这个自适应辛普森积分求面积。
辛普森积分法是用二次函数来模拟这个函数，是在数值分析上，根据经验发现效果非常好的一个方法。
那么，可以这样近似

$\int_a^bf(x)dx \approx \frac{b-a}{6}\biggl(f(a) + 4f(\frac{a+b}{2})+f(b)\biggr)$

如果原函数本身的次数不超过2，那么可以求出很精确的值。

3074. 自适应辛普森积分

题意：给定两个整数 a,b，请计算如下积分：

$\int_a^b \frac{Sin(x)}{x}dx$

#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
using namespace std;

const double eps = 1e-12;

double f(double x) {
	return sin(x) / x;
}

double simpson(double l, double r) {
	auto mid = (l + r) / 2;
	return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}

double asr(double l, double r, double s) {
	auto mid = (l + r) / 2;
	auto left = simpson(l, mid), right = simpson(mid, r);
	if (fabs(left + right - s) < eps) return left + right;
	return asr(l, mid, left) + asr(mid, r, right);
}

int main() {
	double l, r;
	scanf("%lf%lf", &l, &r);
	printf("%.6f\n", asr(l, r, simpson(l, r)));
	return 0;
}

3069. 圆的面积并

给出 N 个圆，求其面积并。
如果用扫描线去写，会非常麻烦，因为面积很难算。但是需要指出的是，辛普森积分不一定比扫描线简单。一般是有弧线时用辛普森积分相对简单。大雪菜还说，如果辛普森积分卡精度，那么先把坐标系随意旋转一个角度就行。
这个题，设 $f(x_0)$ 为直线 $x=x_0$ 与圆交集的长度。

#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>

#define x first
#define y second

using namespace std;

const double eps = 1e-8;
typedef pair<double, double> P;
const int maxn = 1010;

int N;
struct Circle {
	P o;
	double R;
}c[maxn];

P q[maxn];

int dcmp(double x, double y) {
	if (fabs(x - y) < eps) return 0;
	if (x < y) return -1;
	return 1;
}

double f(double x) {
	int cnt = 0;
	for (int i = 0; i < N; i++) {
		auto X = fabs(x - c[i].o.x), R = c[i].R;
		if (dcmp(X, R) < 0) {
			auto Y = sqrt(R * R - X * X);
			q[cnt++] = { c[i].o.y - Y, c[i].o.y + Y };
		}
	}
	if (!cnt) return 0;
	sort(q, q + cnt);
	double res = 0, st = q[0].x, ed = q[0].y;
	for (int i = 1; i < cnt; i++) {
		if (q[i].x <= ed) ed = max(q[i].y, ed);
		else {
			res += ed - st;
			st = q[i].x, ed = q[i].y;
		}
	}
	return res + ed - st;
}

double simpson(double l, double r) {
	auto mid = (l + r) / 2;
	return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}

double asr(double l, double r, double s) {
	auto mid = (l + r) / 2;
	auto left = simpson(l, mid), right = simpson(mid, r);
	if (fabs(left + right - s) < eps) return left + right;
	return asr(l, mid, left) + asr(mid, r, right);
}

int main() {
	scanf("%d", &N);
	for (int i = 0; i < N; i++) {
		scanf("%lf%lf%lf", &c[i].o.x, &c[i].o.y, &c[i].R);
	}
	double l = -2000, r = 2000;
	printf("%.3f\n", asr(l, r, simpson(l, r)));
	return 0;
}

zhezhidashi

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
算法模板（7）：计算几何（2）

已知平面上不共线的一组点的坐标，求覆盖这组点的面积最小的矩形。输出矩形的面积和四个顶点的坐标。上面那个证明出来，矩形至少有一条边和凸包的一条边重合。那么，我们其实就是枚举凸包的每一条边，然后找到离这条边最远的边，即投影最大的。然后，找两边的点时，就是找投影最小的。这样矩形的三个点就找到了，最后一个点可以根据矩形的边的关系算出来。
复制链接

扫一扫