Re:从零开始的动态凸壳
动态凸包概述
动态凸包(Dynamic Convex Hull)问题是指在二维平面上,维护一个点集的凸包结构,并支持在点集中动态地插入或删除点,同时能够高效地查询当前凸包的状态或属性(如周长、面积、顶点集合等)。该问题在计算几何、计算机图形学、机器人路径规划、动态规划等领域具有广泛应用。
前置知识
向量
定义 1
向量向量顾名思义是具有方向(向)和大小(量)的量,通常使用有向线段(带箭头的线段)表示,如下图:
表示
向量有三种表示方法,如下图的向量:
其可表示为:
- \(\vec a\);
- \(\vec{AB}\);
- \((-2.5,1)\)。
其中第三种被称为向量的坐标表示。向量的坐标为终点坐标减起点坐标。
定义 2
如两向量方向相反称这两向量互为相反向量记作 \(\vec v = - \vec u\)。
运算
加减
向量的加法可以看成是一段连续位移的叠加,如图:其中 \(\vec v + \vec u = \vec w\),这种加法法则被称为三角形法则。
由此引申可得平行四边形法则,如图:其中 \(\vec v \mathop{//} f, \vec u \mathop{//} g, \vec v + \vec u = \vec w\)。
向量的减法可以看作加法的逆运算,例如 \(\vec v - \vec u = \vec w\) 表示从 \(\vec v\) 中减去 \(\vec u\) 的效果,等价于 \(\vec v + (-\vec u)\)。
容易得到向量加减的坐标运算为 \((x,y)\pm(a,b)=(x\pm a,y\pm b)\)。
模
向量模为向量的长度,符号与绝对值相同,如图:
其中 \(\left | \vec u \right | = \frac{\sqrt{29}}{2}\)。
容易得到向量的模的坐标运算为 \(\left | (a,b) \right | = \sqrt{a^2+b^2}\)。
定义 3
模为 \(1\) 的向量为单位向量,一般用 \(\vec e\) 表示。
叉乘
在二维空间中,向量的叉乘表示两个向量构成的平行四边形的有向面积,如图:其中 \(\vec w \times \vec v = S_a = -2.25\)。可以得到 \(\vec w \times \vec v = |w||v|\sin<w,v>\)。
可以得到叉乘的坐标运算为 \((a,b)\times(x,y)=ay-bx\)。
证明:
设 \(x\) 轴的单位向量为 \(\vec e_1\),\(y\) 轴的为 \(\vec e_2\)。\((a,b)\times(x,y) = (a\vec e_1 + b\vec e_2)\times(x\vec e_1 + y\vec e_2)\)。
展开得 \(a x \vec e_1 \times \vec e_1 + a y \vec e_1 \times \vec e_2 + b x \vec e_2 \times \vec e_1 + b y \vec e_2 \times \vec e_2\)。
由于 \(\vec e_1 \times \vec e_1 = 0\),\(\vec e_2 \times \vec e_2 = 0\),且 \(\vec e_1 \times \vec e_2 = -\vec e_2 \times \vec e_1 = 1\)(\(\vec e_1\perp\vec e_2\))。
化简得 \(a y - b x\),得证。
叉乘符号则表示两个向量的相对方向关系:
- 若结果为正,表示第二个向量在第一个向量的逆时针方向;
- 若结果为负,表示第二个向量在第一个向量的顺时针方向;
- 若结果为 0,说明两向量共线。
凸包
定义一
在平面上能包含所有给定点的最小凸多边形叫做凸包。
如图:
实际上可以理解为用一个橡皮筋包含住所有给定点的形态。
定义二
凸包中相邻的边若斜率愈来愈大则称下凸壳,相邻的边若斜率愈来愈小称上凸壳。
Andrew's Monotone Chain 算法
以下来自 oi-wiki。
该算法的时间复杂度为 \(O(n \log n)\),其中 \(n\) 为待求凸包点集的大小。复杂度的瓶颈在于对所有点坐标的双关键字排序。
具体过程如下:
- 排序:
将所有点以横坐标为第一关键字、纵坐标为第二关键字进行升序排序。排序后,最小和最大的元素一定在凸包上。 - 构造凸壳:
利用单调栈维护凸壳。由于上下凸壳旋转方向不同,先升序枚举构造下凸壳,再降序枚举构造上凸壳。 - 判断方向:
从栈顶取出两个点 \(S_1\)、\(S_2\)(其中 \(S_1\) 为栈顶),若即将入栈的点 \(P\) 与这两个点构成的方向为右拐(即叉积小于 0):\(\vec{S_2S_1} \times \vec{S_1P} < 0\) 则弹出栈顶 \(S_1\),重复此过程直到满足:\(\vec{S_2S_1} \times \vec{S_1P} \leq 0\) 或栈中只剩一个元素为止。 - 保留边界点(可选):
若需保留凸包边界上的点,可将上述条件中的<
改为≤
,相应地将反向判断条件>
改为≥
。
最终将上下凸壳合并,即可得到完整的凸包。
code
点击查看代码
// stk[] 是整型,存的是下标
// p[] 存储向量或点
tp = 0; // 初始化栈
std::sort(p + 1, p + 1 + n); // 对点进行排序
stk[++tp] = 1;
// 栈内添加第一个元素,且不更新 used,使得 1 在最后封闭凸包时也对单调栈更新
for (int i = 2; i <= n; ++i) {
while (tp >= 2 // 下一行 * 操作符被重载为叉积
&& (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0)
used[stk[tp--]] = 0;
used[i] = 1; // used 表示在凸壳上
stk[++tp] = i;
}
int tmp = tp; // tmp 表示下凸壳大小
for (int i = n - 1; i > 0; --i)
if (!used[i]) {
// ↓求上凸壳时不影响下凸壳
while (tp > tmp && (p[stk[tp]] - p[stk[tp - 1]]) * (p[i] - p[stk[tp]]) <= 0)
used[stk[tp--]] = 0;
used[i] = 1;
stk[++tp] = i;
}
for (int i = 1; i <= tp; ++i) // 复制到新数组中去
h[i] = p[stk[i]];
int ans = tp - 1;
正文 —— Andrew's Monotone Chain 算法变种
现在考虑已经有一个凸壳(为了方便讲解,以上凸壳为例),如下图:
现在考虑插入一个点,分两种情况讨论:
情况一:点在凸壳下或在凸壳上
如图:若 \(G\) 为插入点,\(C\) 为插入点在凸包中的前驱(\(x\) 为第一关键字,\(y\) 为第二关键字排序),\(D\) 为插入点在凸包中的后继,设 \(\vec{OG} = \vec u, \vec{OC} = \vec v, \vec{OD} = \vec w\),可得 \((\vec v - \vec u) \times (\vec w - \vec u) \le 0\),如图:
证明:根据向量叉乘的几何意义,\(\vec{GC} \times \vec{GD}\) 的结果表示以 \(GC\) 和 \(GD\) 为邻边的平行四边形的有向面积。若该值 \(\leq 0\),则说明向量 \(\vec{GD}\) 位于向量 \(\vec{GC}\) 的顺时针方向或共线。
设 \(C\)、\(D\) 为插入点 \(G\) 在凸包中的前驱和后继(按 \(x\) 坐标排序),则凸包上 \(C\) 到 \(D\) 的线段是上凸壳的一部分,其斜率单调递减(上凸壳的性质)。若 \(G\) 在凸包下方或凸包上,则 \(G\) 必须位于线段 \(CD\) 的下方或线段上。此时,线段 \(GC\) 的斜率 \(k_1\) 必然大于等于线段 \(GD\) 的斜率 \(k_2\)(即 \(k_1 \geq k_2\)),否则 \(G\) 会位于凸包外部,导致原凸包不再是“最小凸多边形”。
根据叉乘的坐标运算:
设 \(G(x_0,y_0)\),\(C(x_1,y_1)\),\(D(x_2,y_2)\),则:\(\vec{GC} = (x_1-x_0, y_1-y_0)\),\(\vec{GD} = (x_2-x_0, y_2-y_0)\),\(\vec{GC} \times \vec{GD} = (x_1-x_0)(y_2-y_0) - (x_2-x_0)(y_1-y_0)\)
若 \(k_1 \geq k_2\),即 \(\frac{y_1-y_0}{x_1-x_0} \geq \frac{y_2-y_0}{x_2-x_0}\)(假设 \(x_1 < x_0 < x_2\),分母为正),交叉相乘得:\((y_1-y_0)(x_2-x_0) \geq (y_2-y_0)(x_1-x_0)\)
移项后:\((x_1-x_0)(y_2-y_0) - (x_2-x_0)(y_1-y_0) \leq 0\)
即 \(\vec{GC} \times \vec{GD} \leq 0\)。
此时容易发现,这个点不可能改变凸壳。
情况二:点在凸壳上
此时容易发现,这个点会改变凸壳。
于是我们向左右两边更新凸壳,具体如下。
更新
更新左端时,每次我们不断取出加入点(设为 \(G\))在凸壳中的前驱(设为 \(C\))以及前驱的前驱(设为 \(B\))。若 \(BC\) 相对与 \(BG\) 顺时针旋转,即 \(\vec{BC} \times \vec{BG} \ge 0\),则删除 \(C\)(在 \(G\) 加入之后 \(C\) 若在凸壳上则不满足凸壳的性质),如图:
\(\begin{matrix}\Huge\Longrightarrow\\ \\ \\ \\ \\ \\ \\ \\ \end{matrix}\)
重复此操作直到 \(\vec{BC} \times \vec{BG} < 0\)。
右端如法炮制即可,如图:
\(\begin{matrix}\Huge\Longrightarrow\\ \\ \\ \\ \\ \\ \\ \\ \end{matrix}\)
具体实现见代码。
code
using DB = double;
struct Point { // 定义二维平面点的结构体
DB x, y;
// 返回两点构成的向量
Point operator - (const Point &a) const { return {x - a.x, y - a.y}; }
// 定义以 x 为第一关键字 y 为第二关键字的排序
bool operator < (const Point &a) const {
if (x != a.x) return x < a.x;
return y < a.y;
}
};
using Vector = Point;
// 计算二维向量叉积
DB cross(const Vector &a, const Vector &b) {
return a.x * b.y - b.x * a.y;
}
set<Point> up; // 用有序集合维护上凸壳的点集(按x坐标排序)
// 获取迭代器的前驱(上一个元素)
set<Point>::iterator pre(set<Point>::iterator it) {
return it == up.begin() ? up.end() : --it;
}
// 获取迭代器的后继(下一个元素)
set<Point>::iterator next(set<Point>::iterator it) {
++it;
return it == up.end() ? up.end() : it;
}
// 检查点p是否在上凸壳内部或边界上
bool check(const Point &p) {
auto it = up.lower_bound(p); // 找到第一个x >= p.x的点
if (it == up.end()) return false; // 无后继点,不在凸壳内
if (it == up.begin()) { // 是第一个点
if (it->x == p.x) return it->y >= p.y; // x相同且y >= p.y,在边界上
return false; // 否则不在凸壳内
}
auto prev_it = pre(it); // 获取前驱点
Vector a = *prev_it - p, b = *it - p; // 向量pa和pb
return cross(a, b) <= 0; // 叉积<=0表示p在pa和pb构成的线段下方或共线
}
// 尝试删除凸壳上的点it,返回是否成功删除
bool erase(set<Point>::iterator it) {
if (it == up.end() || it == up.begin()) return false; // 首尾点不删除
auto prev_it = pre(it), next_it = next(it);
if (next_it == up.end()) return false; // 无后继点,无法删除
Vector a = *it - *prev_it, b = *next_it - *it; // 向量prev->it和it->next
if (cross(a, b) >= 0) { // 叉积>=0表示三点共线或右拐,it点冗余
up.erase(it);
return true;
}
return false;
}
// 插入点p到上凸壳中,维护凸壳结构
void insert(const Point &p) {
if (check(p)) return; // 若p在凸壳内或边界上,直接返回
auto it = up.insert(p).first; // 插入p并获取迭代器
auto prev_it = pre(it), next_it = next(it); // 获取前驱和后继
// 从插入点左侧开始删除冗余点(左链优化)
for (auto p = pre(it); p != up.end() && p != up.begin(); p = pre(it))
if (!erase(p)) break; // 无法删除时停止
// 从插入点右侧开始删除冗余点(右链优化)
for (auto a = next(it); a != up.end(); a = next(it))
if (!erase(a)) break; // 无法删除时停止
}
例题 ——P2521 [HAOI2011] 防线修建
思路
定睛一看,这不是动态凸壳板子吗!
考虑题中所给的操作时删除,于是离线,把删点化加点,再在处理时更新长度即可。
code
#include <bits/stdc++.h>
using namespace std;
using DB = double;
struct Point {
DB x, y;
Point operator - (const Point &a) const { return {x - a.x, y - a.y}; }
bool operator < (const Point &a) const {
if (x != a.x) return x < a.x;
return y < a.y;
}
};
using Vector = Point;
DB cross(const Vector &a, const Vector &b) {
return a.x * b.y - b.x * a.y;
}
DB length(const Vector &a) { // 计算向量的模长(两点间距离)
return sqrt(a.x * a.x + a.y * a.y);
}
set<Point> up;
DB current_length = 0; // 维护当前上凸壳的总边长,避免每次查询时重新计算
set<Point>::iterator pre(set<Point>::iterator it) {
return it == up.begin() ? up.end() : --it;
}
set<Point>::iterator next(set<Point>::iterator it) {
++it;
return it == up.end() ? up.end() : it;
}
bool check(const Point &p) {
auto it = up.lower_bound(p);
if (it == up.end()) return false;
if (it == up.begin()) {
if (it->x == p.x) return it->y >= p.y;
return false;
}
auto prev_it = pre(it);
Vector a = *prev_it - p, b = *it - p;
return cross(a, b) <= 0;
}
bool erase(set<Point>::iterator it) {
if (it == up.end() || it == up.begin()) return false;
auto prev_it = pre(it), next_it = next(it);
if (next_it == up.end()) return false;
Vector a = *it - *prev_it, b = *next_it - *it;
if (cross(a, b) >= 0) {
// 更新凸壳总长度:减去被删除点两侧的边,加上直接连接前后点的新边
current_length -= length(*it - *prev_it);
current_length -= length(*next_it - *it);
current_length += length(*next_it - *prev_it);
up.erase(it);
return true;
}
return false;
}
void add_point(const Point &p) {
if (check(p)) return;
auto it = up.insert(p).first;
auto prev_it = pre(it), next_it = next(it);
// 调整总长度:若前后点存在,先减去它们原有的边,再加上与新点连接的两边
if (prev_it != up.end() && next_it != up.end())
current_length -= length(*next_it - *prev_it);
if (prev_it != up.end())
current_length += length(*it - *prev_it);
if (next_it != up.end())
current_length += length(*next_it - *it);
for (auto p = pre(it); p != up.end() && p != up.begin(); p = pre(it))
if (!erase(p)) break;
for (auto a = next(it); a != up.end(); a = next(it))
if (!erase(a)) break;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
int n, x, y;
cin >> n >> x >> y;
Point capital = {DB(x), DB(y)}; // 首都点坐标
// 初始化上凸壳:插入左右端点(0,0)、(n,0)和首都点,构成初始凸壳
up.insert({0, 0});
up.insert({DB(n), 0});
up.insert(capital);
// 计算初始凸壳长度:(0,0)到首都,加上首都到(n,0)的距离
current_length = length(capital - Point{0, 0}) + length(Point{DB(n), 0} - capital);
int m;
cin >> m;
vector<Point> cities(m); // 存储所有城市的坐标
for (int i = 0; i < m; ++i) {
cin >> cities[i].x >> cities[i].y;
}
int q;
cin >> q;
vector<pair<int, int>> queries(q); // 存储查询操作:1表示移除城市,2表示查询长度
vector<bool> active(m, true); // 标记城市是否存在
for (int i = 0; i < q; ++i) {
int op;
cin >> op;
if (op == 1) {
int u;
cin >> u;
queries[i] = {1, u - 1}; // 记录要移除的城市
active[u - 1] = false; // 标记该城市为不存在
} else {
queries[i] = {2, -1}; // 记录查询操作
}
}
// 初始添加所有存在城市(未被提前移除的城市)到凸壳中
for (int i = 0; i < m; ++i)
if (active[i]) {
add_point(cities[i]);
}
vector<DB> ans; // 存储查询结果
// 逆序处理查询:将"移除城市"转化为"添加城市",简化动态维护逻辑
for (int i = q - 1; i >= 0; --i)
if (queries[i].first == 2) {
ans.push_back(current_length); // 记录当前凸壳长度
} else {
int u = queries[i].second;
add_point(cities[u]); // 恢复被移除的城市(逆序处理时变为添加)
}
reverse(ans.begin(), ans.end()); // 反转结果,恢复原查询顺序
cout << fixed << setprecision(2); // 保留两位小数
for (DB val : ans) {
cout << val << '\n';
}
return 0;
}