计算几何 学习笔记
Visits: 635
一直莫名不想碰的一个专题
基础知识
点、向量模板
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 《计算几何》课件
- 吴清月 《计算几何》课件