计算几何 学习笔记

一直莫名不想碰的一个专题

基础知识

点、向量模板

const ld eps = 1e-10L, PI = acos(-1);

struct P {
	ld x, y;
	inline P() {}
	inline P(ld x, ld y) : x(x), y(y) {}
	inline P &operator += (P o) { return x += o.x, y += o.y, *this; }
	inline P &operator -= (P o) { return x -= o.x, y -= o.y, *this; }
	inline P &operator *= (ld o) { return x *= o, y *= o, *this; }
	inline P &operator /= (ld o) { return x /= o, y /= o, *this; }
	inline friend P operator + (P a, P b) { return a += b; }
	inline friend P operator - (P a, P b) { return a -= b; }
	inline friend P operator * (P a, ld b) { return a *= b; }
	inline friend P operator / (P a, ld b) { return a /= b; }
	inline friend bool operator < (P a, P b) { return fabs(a.x - b.x) < eps ? a.y < b.y : a.x < b.x; }
	inline friend bool operator == (P a, P b) { return !(a < b) && !(b < a); }
	inline friend ld operator * (P a, P b) { return a.x * b.x + a.y * b.y; }
	inline friend ld operator % (P a, P b) { return a.x * b.y - a.y * b.x; }
	inline friend ld operator ^ (P a, P b) { return a * b / a.l() / b.l(); }
	inline ld a() { return atan2(y, x); }
	inline ld l() { return sqrt(*this * *this); }
	inline void r(ld o) { ld s = sin(o), c = cos(o), X = x * c - y * s, Y = x * s + y * c; x = X, y = Y; }
};

向量点积

上面模板中的 * 运算。

两个向量 $p_1,p_2$ 的点积定义为 $p_1,p_2$ 的模长乘它们夹角的 $\cos$，可以解释为 $p_1,p_2$ 做的功的大小。

$\cos$ 是偶函数，因此点积满足交换律。

$$p_1 * p_2 = x_1x_2 + y_1y_2 = p_2 * p_1$$

向量叉积

上面模板中的 % 运算。

两个向量 $p_1,p_2$ 的叉积可以解释为 $p_1,p_2$ 叉出来的平行四边形的有向面积。

叉积没有交换律，但交换后为原来的相反数。

$$p_1 \times p_2 = x_1y_2 - x_2y_1 = -p_2 \times p_1$$

向量夹角

上面模板中的 ^ 运算。

向量极角

上面模板中的 a 运算。

向量模长

上面模板中的 l 运算。

向量旋转

上面模板中的 r 运算。

向量绕起点逆时针旋转角度 $a$，公式为 $X = x \cos(a) - y \sin(a), Y = x \sin(a) + y \cos(a)$，本质上是个矩阵乘法。

$$\begin{bmatrix} X & Y \end{bmatrix} = \begin{bmatrix} x & y \end{bmatrix} \times \begin{bmatrix}\text{cos}(a) & \text{sin}(a)\\\\ -\text{sin}(a) & \text{cos}(a)\end{bmatrix}$$

直线模板

struct L {
	P a, b;
	inline L() {}
	inline L(P a, P b) : a(a), b(b) {}
	inline friend P operator * (L a, L b) { return a.a + a.b * (b.b % (a.a - b.a) / (a.b % b.b)); }
	inline friend ld operator / (L a, P b) { return fabs(a.b % (b - a.a)) / a.b.l(); }
	inline friend bool operator < (L a, L b) { return fabs(a.b.a() - b.b.a()) < eps ? (b.a - a.a) % a.b < 0 : a.b.a() < b.b.a(); }
};

如果被卡常，可以把 ld 换成 double，同时提前存下 b.a() 的值。

求两条直线的交点

上面模板中的 * 运算。

设直线 $a$ 的关键点 $a.a = P$，方向向量 $a.b = V$，直线 $b$ 的关键点 $b.a = Q$，方向向量 $b.b = W$。

所以交点 $(x,y)$ 满足 $x = P_x + t_1V_x = Q_x + t_2W_x, y = P_y + t_1V_y = Q_y + t_2W_y$。

解得 $t_1 = \frac{(P_y - Q_y)W_x - (P_x - Q_x)W_y}{V_xW_y - V_yW_x} = (W \times (P - Q) / (V \times W)) = (b.b \times (a.a - b.a) / (a.b \times b.b))$。

然后...就做完了...

顺便一提，线段交点就把直线交点求出来然后判一下就完了...

点到直线的距离

上面模板中的 / 运算。

用叉积（平行四边形面积）除以低即可。

欧拉定理

平面图满足：点数 $+$ 面数 $=$ 边数 $+ 2$。

判断一个点是否在任意多边形内部

从这个点出发引一条射线，如果这条射线与多边形有奇数个交点则在内部，否则在外部。

为了避免射线与多边形上的点重合，可以将射线的斜率设为无理数（随便乱打一个数）。

三角剖分求多边形面积

把多边形三角剖分，分别计算有向面积，然后求和。

二维凸包

二维上把若干给定点包围在内部的、面积最小的凸多边形。

先排序，然后单调栈维护一下就好了，时间复杂度 $\mathcal O(n \log n)$。

【模板】P2742 【模板】二维凸包 / [USACO5.1]圈奶牛Fencing the Cows

inline vector<P> convex_hull(vector<P> a) {
	sort(a.begin(), a.end());
	int n = unique(a.begin(), a.end()) - a.begin();
	vector<P> p;
	for (int i = 0; i < n; i++) {
		while (p.size() > 1u &&
			(p.back() - p[p.size()-2]) % (a[i] - p[p.size()-2]) < eps)
				p.pop_back();
		p.pb(a[i]);
	}
	ui m = p.size();
	for (int i = n - 2; ~i; i--) {
		while (p.size() > m &&
			(p.back() - p[p.size()-2]) % (a[i] - p[p.size()-2]) < eps)
				p.pop_back();
		p.pb(a[i]);
	}
	p.pop_back();
	return p;
}

int main() {
	int n;
	scanf("%d", &n);
	vector<P> a;
	ld x, y;
	for (int i = 1; i <= n; i++)
		scanf("%Lf %Lf", &x, &y), a.pb(P(x, y));
	a = convex_hull(a);
	a.pb(a[0]);
	ld ans = 0;
	for (ui i = 1; i < a.size(); i++) ans += (a[i] - a[i-1]).l();
	printf("%.2Lf\n", ans);
	return 0;
}

旋转卡壳

用于求平面最远点对。

首先求出凸包，目的是把所有点按逆时针排序。然后枚举逆时针边，定义 $T$ 为当前的对踵点，若 $T$ 逆时针变动能增大到当前边的距离则改变 $T$，直到继续改变会导致距离减小为止。每次求出边（设为 $AB$）的对踵点 $T_i$ 后，分别用 $|AT_i|,|BT_i|$ 来更新答案，即凸多边形的直径。

（对踵点：如果过凸多边形上两点作一对平行线，使得整个多边形都在这两条线之间，那么这两个点被称为一对对踵点。）

时间复杂度 $\mathcal O(n \log n)$。

【模板】P1452 Beauty Contest

const int N = 5e4 + 7;
int n, m;
P a[N], p[N];

inline ld S(P x, P a, P b) { return Abs((a - x) % (b - x) / (a - b).l()); }

int main() {
	ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
	cin >> n;
	for (int i = 1; i <= n; i++) cin >> a[i].x >> a[i].y;
	sort(a + 1, a + n + 1);
	n = unique(a + 1, a + n + 1) - (a + 1);
	for (int i = 1; i <= n; i++) {
		while (m > 1 && (p[m] - p[m-1]) % (a[i] - p[m-1]) <= 0) --m;
		p[++m] = a[i];
	}
	int o = m;
	for (int i = n - 1; i; i--) {
		while (m > o && (p[m] - p[m-1]) % (a[i] - p[m-1]) <= 0) --m;
		p[++m] = a[i];
	}
	if (m < 3) return cout << 0 << endl, 0;
	if (m == 3) {
		ld ans = (p[1] - p[2]).l();
		cout << (ll)(ans * ans + 0.5) << endl;
		return 0;
	}
	int t = 1;
	ld ans = 0;
	for (int i = 1; i < m; i++) {
		while (S(p[t], p[i], p[i+1]) <= S(p[t+1], p[i], p[i+1])) t = t % (m - 1) + 1;
		ans = Max(ans, Max((p[t] - p[i]).l(), (p[t] - p[i+1]).l()));
	}
	cout << (ll)(ans * ans + 0.5) << endl;
	return 0;
}

注意当凸包退化成一个点或一条线时需要特判。

半平面交

给出若干个半平面（有向直线的左侧），求它们的交。

做法与二维凸包比较像。先极角排序，然后单调队列维护即可，时间复杂度 $\mathcal O(n \log n)$。

【模板】P4196 [CQOI2006]凸多边形

半平面交存在三种情况：无解，有解且有界，有解但无界。

下面的 half_plane 中，如果无解，则返回的 deque<P> 大小为 $0$，否则为半平面交。

本题中不可能出现有解但无界的情况，如果可能出现，我们需要在最外面套上一个巨大的正方形转化成有界的情况，最后如果有解，但有边在外层的正方形上，则说明无界。

inline deque<P> half_plane(vector<L> l) {
	sort(l.begin(), l.end());
	deque<P> r;
	deque<L> q;
	for (L o : l) {
		while (r.size() && (r.back() - o.a) % o.b > 0)
			q.pop_back(), r.pop_back();
		while (r.size() && (r.front() - o.a) % o.b > 0)
			q.pop_front(), r.pop_front();
		if (q.size() && fabs(o.b % q.back().b) < eps) {
			if (fabs(o.b.a() - q.back().b.a()) > eps)
				return r.clear(), r;
			q.pop_back();
			if (r.size()) r.pop_back();
		}
		if (q.size()) r.pb(q.back() * o);
		q.pb(o);
	}
	while (r.size() && (r.back() - q.front().a) % q.front().b > 0)
		q.pop_back(), r.pop_back();
	if (q.size() < 3u) r.clear();
	else r.pb(q.front() * q.back());
	return r;
}

int main() {
	int n;
	vector<P> p;
	vector<L> l;
	rd(n);
	for (int i = 1, m, x, y; i <= n; i++) {
		rd(m), p.clear();
		while (m--) rd(x), rd(y), p.pb(P(x, y));
		p.pb(p[0]);
		for (ui i = 1; i < p.size(); i++) l.pb(L(p[i-1], p[i] - p[i-1]));
	}
	deque<P> r = half_plane(l);
	if (!r.size()) return prints("0.000"), 0;
	ld ans = r.back() % r.front();
	while (r.size() > 1u) {
		P o = r.back();
		r.pop_back();
		ans += r.back() % o;
	}
	cout << fixed << setprecision(3) << ans / 2 << endl;
	return 0;
}

扫描线

这一一个跟计算几何没关系却又有关系的算法...

线段树维护，由于扫描线的特殊性质可以不用下传标记。

【模板】P5490 【模板】扫描线

#define ls (p << 1)
#define rs (ls | 1)
#define mid ((t[p].l + t[p].r) >> 1)
#define ll long long

const int N = 1e5 + 6;
int n;
struct P {
	int x, y, z, k;
	inline bool operator < (const P &o) const {
		return x < o.x;
	}
} a[N<<1];
int raw[N<<1];
map<int, int> val;
struct T {
	int l, r, cnt, len;
} t[N<<3];
ll ans;

void build(int p, int l, int r) {
	t[p].l = l, t[p].r = r;
	if (l == r) return;
	build(ls, l, mid), build(rs, mid + 1, r);
}

void change(int p, int l, int r, int k) {
	if (l <= t[p].l && r >= t[p].r) return t[p].len = ((t[p].cnt += k) ? raw[t[p].r+1] - raw[t[p].l] : (t[p].l == t[p].r ? 0 : t[ls].len + t[rs].len)), void();
	if (l <= mid) change(ls, l, r, k);
	if (r > mid) change(rs, l, r, k);
	t[p].len = (t[p].cnt ? raw[t[p].r+1] - raw[t[p].l] : t[ls].len + t[rs].len);
}

int main() {
	rd(n);
	for (int i = 1; i <= n; i++) {
		int k = i << 1, y, z;
		rd(a[k-1].x), rd(y), rd(a[k].x), rd(z);
		raw[k-1] = a[k-1].y = a[k].y = y;
		raw[k] = a[k-1].z = a[k].z = z;
		a[k-1].k = 1, a[k].k = -1;
	}
	n <<= 1;
	sort(raw + 1, raw + n + 1);
	int m = unique(raw + 1, raw + n + 1) - (raw + 1);
	for (int i = 1; i <= m; i++) val[raw[i]] = i;
	sort(a + 1, a + n + 1);
	build(1, 1, m - 1);
	for (int i = 1; i < n; i++) {
		int y = val[a[i].y], z = val[a[i].z] - 1;
		change(1, y, z, a[i].k);
		ans += 1ll * t[1].len * (a[i+1].x - a[i].x);
	}
	print(ans);
	return 0;
}

曼哈顿距离和切比雪夫距离

设 $A(x_1, y_1), B(x_2, y_2)$。

• 曼哈顿距离：$|x_{1}-x_{2}|+|y_{1}-y_{2}|$。
• 切比雪夫距离：$\max(|x_{1}-x_{2}|,|y_{1}-y_{2}|)$。

相互转化

• 曼哈顿距离 $\to$ 切比雪夫距离：$(x, y) \to (x + y, x - y)$。
• 切比雪夫距离 $\to$ 曼哈顿距离：$(x, y) \to (\frac{x+y}{2},\frac{x-y}{2})$。

【例题】P3964 [TJOI2013]松鼠聚会

const int N = 1e5 + 7;
int n, x[N], y[N], p[N], q[N];
ll sx[N], sy[N], ans = 4e18;

int main() {
	rd(n);
	for (int i = 1, a, b; i <= n; i++) rd(a), rd(b), x[i] = p[i] = a + b, y[i] = q[i] = a - b;
	sort(p + 1, p + n + 1), sort(q + 1, q + n + 1);
	for (int i = 1; i <= n; i++) sx[i] = sx[i-1] + p[i], sy[i] = sy[i-1] + q[i];
	for (int i = 1; i <= n; i++) {
		int px = lower_bound(p + 1, p + n + 1, x[i]) - p, py = lower_bound(q + 1, q + n + 1, y[i]) - q;
		ans = min(ans, 1ll * px * x[i] - sx[px] + sx[n] - sx[px] - 1ll * (n - px) * x[i] + 1ll * py * y[i] - sy[py] + sy[n] - sy[py] - 1ll * (n - py) * y[i]);
	}
	print(ans >> 1);
	return 0;
}

参考资料

• 刘汝佳 & 陈锋 《算法竞赛进阶指南·入门经典》
• 李煜东 《算法竞赛进阶指南》
• wjyyy 《计算几何 学习笔记合集【凸包】【旋转卡壳】【半平面交】》
• 自为风月马前卒 曼哈顿距离与切比雪夫距离及其相互转化
• mcfx 《计算几何》课件
• 吴清月 《计算几何》课件