凸包
算法模型
想象在一个平面上钉下了 \(n\) 个钉子。现在有一根橡皮筋,我们把它撑开,期望在松手之后橡皮筋可以收缩,包住所有的 \(n\) 个钉子。事实上,这正是一个凸包。如下图:
凸包定义为:周长最小的包含点集中所有点的 凸多边形 。即使存在某个凹多边形的周长与凸包相等且可以包含所有点,这个凹多边形也一定不是凸包。如下图,这个凹多边形不是该点集的凸包:
这里只有用法 Andrew 求凸包的方法。
前置知识
向量。
向量叉积
两个向量 \(\overrightarrow{x1, y1}\) 和 \(\overrightarrow{x2, y2}\) 的外积称之为叉积,它的结果是一个向量。叉积的计算方法是 \((x1y2) - (x2y1)\) ,记为 \(\overrightarrow{x1, y1} \times \overrightarrow{x2, y2}\)。
对于两个向量 \(\overrightarrow{ab}, \overrightarrow{bc}\) ,如果它们的叉积 \(> 0\) ,说明 \(\overrightarrow{ab}\) 在 \(\overrightarrow{bc}\) 的顺时针方向;如果它们的叉积 \(= 0\),说明 \(\overrightarrow{ab}\) 和 \(\overrightarrow{bc}\) 共线;如果它们的叉积 \(< 0\),说明 \(\overrightarrow{ab}\) 在 \(\overrightarrow{bc}\) 的逆时针方向。
算法思想
凸包
首先把点集按照以 \(x\) 坐标为第一关键字,以 \(y\) 坐标为第二关键字的方式进行双关键字从小到大排序,排序后的第一个点就是我们选出的 极点 。两个关键字的顺序可以调换。如下图,点 \(1\) 就是该点集的极点。
接着,我们从极点开始 逆时针 考虑将每一个点都加入凸包。显然我们排序后的第一个点和最后一个点一定在凸包上。从第二个点开始,我们假设当前点可以加入凸包。
设凸包上此时有 \(m\) 个点,第 \(m - 1\) 个点和第 \(m\) 个点分别是 \(a, b\),当前要加入凸包的点为 \(c\) 。如果向量 \(\overrightarrow{bc}\) 在向量 \(\overrightarrow{ab}\) 的逆时针方向,说明我们此时可以将点 \(c\) 加入凸包。用叉积来表示为 \(\overrightarrow{ab} \times \overrightarrow{bc} > 0\) 。
反之,因为如果存在两点共线的情况,显然我们选择与凸包上的点距离较长的点会更优,所以如果它们的叉积 \(\leq 0\) ,说明此时不合法,把凸包上的第 \(m - 1\) 个点删除,再判断第 \(m - 2\) 个点是否合法……直到凸包上只剩下第 \(1\) 个点为止。最后把当前点加入凸包。
从极点开始做,最终只能求出凸包的一半,称之为下包。所以我们还需要从第 \(n\) 个点开始,再次逆时针求一遍凸包,求出凸包的另一半,称之为上包。
假设下包上点的数量为 \(k\) ,那么我们在求上包的时候不应该弹出倒数第 \(k\) 个元素及之后的点,否则凸包的正确性就会被破坏。上包和下包实质上是可以互换的。最终存储凸包的序列应该是首尾相连的,所以如果点集中点的数量 \(> 1\),我们还需要把凸包上点的个数 \(- 1\) 。最终得到的凸包序列应该是从极点开始逆时针排序的。
值得注意的是,如果存在多点共线的情况,也就是两个向量的叉积 \(= 0\) ,这时我们也应该弹出凸包中的点。否则,因为最坏情况下所有点共线,它们都会入序列 \(2\) 次,所以我们的凸包序列应该开够 \(2\) 倍的空间。
这里提供一组 \(hack\) 数据, \(n = 4, x_i = \{1, 1, 4, 5, 1, 4\}, y_i = \{0, 0, 0, 0, 0, 0\}\) 。
总时间复杂度为 \(O(nlogn)\) 。
参考代码
两点共线不弹出:
#include <cstdio>
#include <cmath>
#include <stack>
#include <algorithm>
using namespace std;
const int maxn = 1e5 + 5;
struct vec {
double x, y;
vec(): x(), y() {}
vec(double _x, double _y): x(_x), y(_y) {}
bool operator < (const vec& rhs) const {
if (x != rhs.x) {
return x < rhs.x;
}
return y < rhs.y;
}
} a[maxn], st[maxn * 2];
int n;
int top;
double cross(vec a, vec b) {
return (a.x * b.y) - (b.x * a.y);
}
double calc(vec a, vec b, vec c) {
vec ab = vec(b.x - a.x, b.y - a.y);
vec ac = vec(c.x - a.x, c.y - a.y);
return cross(ab, ac);
}
double dis(vec a, vec b) {
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
double Andrew() {
sort(a + 1, a + n + 1);
st[0] = a[1];
st[1] = a[2];
top = 1;
for (int i = 3; i <= n; i++) {
while (top && calc(st[top - 1], st[top], a[i]) < 0) {
top--;
}
st[++top] = a[i];
}
st[++top] = a[n - 1];
for (int i = n - 2; i >= 1; i--) {
while (top && calc(st[top - 1], st[top], a[i]) < 0) {
top--;
}
st[++top] = a[i];
}
return top;
}
int main() {
// freopen("P2742_1.in", "r", stdin);
double ans = 0;
scanf("%d", &n);
for (int i = 1; i <= n; i++) {
scanf("%lf%lf", &a[i].x, &a[i].y);
}
if (n >= 2) {
int m = Andrew();
for (int i = 1; i < m; i++) {
ans += dis(st[i], st[i + 1]);
}
ans += dis(st[1], st[m]);
printf("%.2lf\n", ans);
} else {
puts("0");
}
return 0;
}
两点共线则弹出:
#include <cstdio>
#include <cmath>
#include <stack>
#include <algorithm>
using namespace std;
const int maxn = 1e5 + 5;
struct vec {
double x, y;
vec(): x(), y() {}
vec(double _x, double _y): x(_x), y(_y) {}
bool operator < (const vec& rhs) const {
if (x != rhs.x) {
return x < rhs.x;
}
return y < rhs.y;
}
} a[maxn], st[maxn];
int n;
double cross(vec a, vec b) {
return (a.x * b.y) - (b.x * a.y);
}
double calc(vec a, vec b, vec c) {
vec ab = vec(b.x - a.x, b.y - a.y);
vec ac = vec(c.x - a.x, c.y - a.y);
return cross(ab, ac);
}
double dis(vec a, vec b) {
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
double Andrew() {
int m = 0;
sort(a, a + n);
for (int i = 0; i < n; i++) {
while (m > 1 && calc(st[m - 1], st[m - 2], a[i]) <= 0) {
m--;
}
st[m++] = a[i];
}
int k = m;
for (int i = n - 2; i >= 0; i--) {
while (m > k && calc(st[m - 1], st[m - 2], a[i]) <= 0) {
m--;
}
st[m++] = a[i];
}
if (n > 1) {
m--;
}
return m;
}
int main() {
// freopen("P2742_1.in", "r", stdin);
double ans = 0;
scanf("%d", &n);
for (int i = 0; i < n; i++) {
scanf("%lf%lf", &a[i].x, &a[i].y);
}
if (n >= 2) {
int m = Andrew();
for (int i = 0; i < m - 1; i++) {
ans += dis(st[i], st[i + 1]);
}
ans += dis(st[0], st[m - 1]);
printf("%.2lf\n", ans);
} else {
puts("0");
}
return 0;
}
旋转卡壳
求点集中最远的点对距离,即求凸包直径。
这个问题同样可以用 \(Andrew\) 算法解决。显然距离最远的两点只有可能在凸包上。假设我们已经求出了该点集的凸包。我们可以想象有两条 平行线 将凸包严丝合缝地卡在中间,而这两条线则不断按逆时针顺序转动。在某个时刻这两条线的位置应该有两种情况:
- 某一条线与凸包上的一条线段部分重合,另一条线卡在了凸包中的点上。
- 两条线都卡在了凸包中的点上。
因为第二种情况不好考虑,并且当两点 \(a, b\) 之间的距离最长时,我们在考虑 \(b\) 与其相邻点 \(c\) 的边 \(bc\) 到 \(a\) 的距离时统计了这种情况,所以我们只需要考虑第一种情况即可。
假设现在一条平行线与凸包上的线段 \(ab\) 部分重合。这时我们考虑另一条平行线卡住的点 \(c\) ,我们将 \(c\) 称为 对踵点。我们需要做的就是枚举每一条线段,然后找出与其两端点距离最长的对踵点。
我们可以发现,当线段的位置被固定时,它最多有一个对踵点,并且凸包上的点到这条线段的距离应该是先递增到达峰值点,转而再递减的。也就是说,凸包上的点到该线段的距离是成单峰函数的。
因为对踵点到达对应线段的距离比其他点要大,也就是说,以对踵点和线段的两个端点为顶点形成的三角形面积应该是最大的。我们可以利用叉积的几何意义来求解对踵点:
因为叉积可以用来求解以向量 \(\overrightarrow{ab}\) 和向量 \(\overrightarrow{bc}\) 为邻边组成的平行四边形的面积或是三角形的面积。我们已知三点 \(a, b, c\),有结论 $|\overrightarrow{ab} \times \overrightarrow{bc}| \div 2 = $ 以 \(a, b, c\) 为顶点的三角形面积。
因为三角形的两个顶点已经固定,第三个顶点(对踵点)决定三角形的高。对踵点越远,高越长,三角形面积越大。想要求出最远的对踵点,等价于求出令该三角形面积最大的点。因此我们可以枚举对踵点,叉积最大的点就是对应线段的对踵点。
因为对踵点到线段的距离成单峰函数,因此当我们逆时针枚举下一条线段时,只需要从当前的对踵点开始继续枚举,找到第一个当前叉积大于下一叉积的点即可。
最终时间复杂度为 \(O(nlogn)\) 。
参考代码
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const int maxn = 5e4 + 5;
struct vec {
double x, y;
vec(): x(), y() {}
vec(double _x, double _y): x(_x), y(_y) {}
bool operator < (const vec& rhs) const {
if (x != rhs.x) {
return x < rhs.x;
}
return y < rhs.y;
}
} a[maxn], st[maxn];
int n, m;
double dis(vec a, vec b) {
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
double cross(vec a, vec b) {
return (a.x * b.y) - (b.x * a.y);
}
double side(vec a, vec b, vec c) {
vec ab = vec(b.x - a.x, b.y - a.y);
vec ac = vec(c.x - a.x, c.y - a.y);
return cross(ab, ac);
}
int andrew() {
sort(a, a + n);
for (int i = 0; i < n; i++) {
while (m > 1 && side(st[m - 2], st[m - 1], a[i]) <= 0) {
m--;
}
st[m++] = a[i];
}
int k = m;
for (int i = n - 2; i >= 0; i--) {
while (m > k && side(st[m - 2], st[m - 1], a[i]) <= 0) {
m--;
}
st[m++] = a[i];
}
if (n > 1) {
m--;
}
return m;
}
double calc() {
int j = 1;
double ans = 0;
for (int i = 0; i < m; i++) {
while (side(st[i], st[i + 1], st[j]) < side(st[i], st[i + 1], st[j + 1])) {
j = (j + 1) % m;
}
ans = max(ans, max(dis(st[i], st[j]), dis(st[i + 1], st[j])));
}
return ans;
}
int main() {
scanf("%d", &n);
for (int i = 0; i < n; i++) {
scanf("%lf%lf", &a[i].x, &a[i].y);
}
andrew();
printf("%.0lf\n", calc());
return 0;
}