计算几何学习笔记 [基础 III]

凸包

这里介绍 Andrew 算法求凸包。

先把所有点按照 \(x\) 坐标为第一关键字，\(y\) 坐标为第二关键字排序。

之后维护一个栈。排序后最小的元素一定在凸包上，所以把它进栈。

之后向后升序扫描每一个点。设栈顶的点是 \(A\)，栈顶下一个点为 \(B\)，当前点为 \(C\)。那么当栈有至少 \(2\) 个元素（也就是 \(A\) 和 \(B\) 存在）时，我们检查：\(\vec{BA} \times \vec{CA}\) 的符号。

• \(\vec{BA} \times \vec{CA} \le 0\)：一直弹栈，并重复检查，直到不满足该条件，或者栈内元素不足 \(2\) 个。
• \(\vec{BA} \times \vec{CA} > 0\)：当前点入栈。

感性理解一下就是从最左下角的点一直加入新点，每次加入往左边转。由于这个凸包我们从左下角点逆时针走一定是不断向左转的，所以这个是对的。

这样其实就求出了凸包的下凸壳。

稍微改动一下，倒序枚举所有点，即可求出上凸壳。

Tips：这里的讲解、实现 不会把在凸包边上的点计算进凸包。如果要把凸包边上的点也算进去，只需要把 \(\le\) 改为 \(<\)，\(>\) 改成 \(\ge\) 即可。

算法的模拟过程将会在代码之后给出。

时间复杂度 \(O(n \log n)\)，瓶颈在排序。

代码：

inline bool operator < (const Point &A, const Point &B) {
	return A.x == B.x ? A.y < B.y : A.x < B.x;
}

int n; Point p[N];
int stk[N], top;
bool used[N];
Point h[N]; int num;

void Andrew()
{
    top = 0;
    sort(p + 1,p + 1 + n);
    stk[++ top] = 1;
    for(int i = 2; i <= n; i ++)
    {
        while(top >= 2 && 
              Cross(p[stk[top]] - p[stk[top - 1]], p[i] - p[stk[top]]) <= 0)
            used[stk[top --]] = 0;
        used[i] = 1;
        stk[++ top] = i;
    }
    int size = top;
    for(int i = n - 1; i > 0; i --)
    {
        if(!used[i])
        {
            while(top > size && 
                  Cross(p[stk[top]] - p[stk[top - 1]], p[i] - p[stk[top]]) <= 0)
                used[stk[top --]] = 0;
            used[i] = 1;
            stk[++ top] = i;
        }
    }
    for(int i = 1; i <= top; i ++)
        h[i] = p[stk[i]];
    num = top - 1;
}

最后的 h[] 数组下标为 1 ~ num 的位置存储了所有在凸包上的点，num + 1 位置上额外存储了凸包上最开始的第一个点（双关键字排序后最小的点）

图解：

旋转卡壳

个人感觉就是一个利用凸包的单调性优化枚举的东西。

求最远点对

经典的求平面内最近点对可以拿分治法做。

当然，也有不少发扬人类智慧的“旋转坐标系跑暴力”之类的做法。

所以最远点对怎么做？

发扬人类智慧乱搞

首先是求一个凸包，这个容易理解。最远点对一定在凸包上。之后我们就变成求凸包直径。

怎么做？这里引入 旋转卡壳 的方法。

逆时针枚举凸包的边，那么离得最远的点必然也在逆时针转，这个过程是具有某种单调性的。这样，就可以把暴力枚举点对的 \(O(n^2)\) 优化到 \(O(n)\)。

具体地，由于求出凸包之后点是逆时针存储的，枚举凸包上的点 \(i\)，点 \(i\) 和 \(i+1\) 就构成边 \((i,i+1)\)，达到了我们枚举凸包边的目的；之后记录一个当前最优点 \(j\)，不断尝试更新它为更优的点 \(j+1\)，每次更新都再尝试更新最远点对距离。这么转一圈就行。

如图，用 \((i,i+1)\) 这个线的平行线去“卡”最优点。

更新最优点的过程也如图。既然只需要判断 \(j\) 和 \(j+1\) 到 \((i,i+1)\) 的距离，那么我们直接把距离的比较转化为同底三角形面积的比较，可以直接用叉积大小判断了。

写的时候要稍微注意细节。

代码：在凸包的基础上加入以下代码：

void work()
{
    int j = 3;
    if(num < 3) 
    {
        ans = dis(h[2], h[1]);
        return;
		// 特判只有两个点
    }
    for(int i = 1; i <= num; i ++)
    {
        while(Cross(h[i + 1] - h[i], h[j] - h[i + 1]) <= 
              Cross(h[i + 1] - h[i], h[j % num + 1] - h[i + 1]))
            j = j % num + 1; // 尝试更新最优点
        ans = max({ans, dis(h[i + 1], h[j]), dis(h[i], h[j])});
		// 尝试更新答案
    }
	
}

需要注意的细节大概就是怎么“转一圈”了。实现方式很多，但是只要能“转一圈”就可以。主要是处理环上最后一个点和第一个点连的边，最优点转一圈从最后一个点回到第一个点之类的东西的时候不要下标越界，其他怎么都好说。

求最小矩形覆盖

题目链接：P3187 [HNOI2007]最小矩形覆盖

大意：\(3 \le n \le 5 \times 10^4\) 个点，求一个面积最小的能覆盖住所有点的矩形，要求 输出面积，以 \(y\) 坐标最小的点（相同则以 \(x\) 坐标最小的点）为第一个点，逆时针顺序输出矩形四个顶点。

先看看怎么求这个矩形：

首先搞一个凸包，之后旋转卡壳维护一个最远点，这个还是可以做到的。但是这样实际上只维护了一对平行线，没有维护矩形。

那么我们还得维护矩形剩下两条边。这两条边都垂直于已经维护的平行线，那么意味着类似于旋转卡壳，这两条平行线也是有一定单调性的。旋转卡壳维护它们经过的凸包上的点即可。

具体地，考虑这个图：

枚举到 \(AB\) 这个边（也就是 \((i,i+1)\)），之后矩形左边过 \(C\)，右边过 \(D\)，\(AB\) 的对边过 \(E\)。

我们肯定是会求 \(AB\) 和 \(E\) 的，关键是怎么把 \(C\) 和 \(D\) 转出来。

因为 \(C\) 在左边，记 \(C\) 对应点在点集中下标为 \(l\)；\(D\) 在右边，记 \(D\) 为 \(r\)。

考虑怎么转 \(D\):

我们主要是要看 \(D\) 这个点在 \(AB\) 方向投影要尽可能靠右。所以我们考虑 使用点积。考虑“试着转一下”，看一下点积的正负即可。

具体地，求 \(\overrightarrow{AB} \cdot \overrightarrow{P_rP_{r+1}}\)。它要是 \(>0\) 那么肯定是往右转，把 \(r\) 更新到 \(r+1\)；否则就不转了，现在的 \(r\) 已经是最右边的 \(D\) 了。

那么怎么求 \(l\)？一样的！但是仔细思考一下，直接变个条件好像不行，它每次好像会转一圈...所以我们倒序转 \(l\) 即可。倒序枚举 \(i\)，每次把 \(l\) 往前转。

所以实现的时候就把旋转拆开了，必须离线记录所有的可能最优的三元组 \((j,l,r)\) 表示对边的点、左边的点、右边的点。不像之前，枚举 \(i\) 的顺序只有一个所以可以直接在线搞。

关键代码：

struct node {
    int j, l, r;
}t[N];

void Rotate1()
{
    int j = 3;
    for(int i = 1; i <= num; i ++)
    {
        while(Cross(h[i + 1] - h[i], h[j] - h[i + 1]) <= 
              Cross(h[i + 1] - h[i], h[j % num + 1] - h[i + 1]))
            j = j % num + 1;  
        t[i].j = j;    
    }
}

void Rotate2()
{
    int r = 2;
    for(int i = 1; i <= num; i ++)
    {
        while(Dot(h[i + 1] - h[i], h[r % num + 1] - h[r]) > 0)
            r = r + 1;
        t[i].r = r;
    }
}

void Rotate3()
{
    int l = num;
    for(int i = num; i; i --)
    {
        int pre = l - 1; if(!pre) pre = num;
        while(Dot(h[i + 1] - h[i], h[pre] - h[l]) < 0)
        {
            l = pre;
            pre = l - 1; if(!pre) pre = num;
        }
        t[i].l = l;
    }
}

之后考虑怎么把矩形的面积和四个端点算出来。

结合上面那个图，各种叉积点积旋转乱搞出各种比例之类的，就能拼出来了。这个说难也不难，说简单的话也要一点点推导和乱搞的技巧才能尽可能简单地实现这一部分。

之后考虑怎么搞这个排序。

我们直接记录下“排序后的第一个点”的下标 \(id\) 作为一个偏移量，每次用 \((i + id) \bmod 4\) 就可以了。注意：不要真的丢进 sort 排序了，那样你的逆时针的顺序就全乱了。

Point Ans[4];
void work()
{
    Rotate1(); Rotate2(); Rotate3();
    for(int i = 1; i <= num; i ++)
    {
        Point A = h[i]; Point B = h[i + 1];
        Vector AB = h[i + 1] - h[i];
        Point C = h[t[i].l];
        Point D = h[t[i].r];
        Point E = h[t[i].j];
        db lenAB = Length(AB);
        db S = fabs(Cross(E - A, B - A));
        db l = fabs(Dot(D - B, B - A) / lenAB) + fabs(Dot(C - A, B - A) / lenAB) + lenAB;
        // r ~ i + 1 + l ~ i + i ~ i + 1
        db h = S / lenAB;
        if(l * h < ans)
        {
            ans = l * h;
            Ans[0] = A + AB * (Dot(C - A, B - A) / lenAB) / lenAB;
            Ans[1] = B + AB * (Dot(D - B, B - A) / lenAB) / lenAB;
            Ans[2] = Ans[0] + Rotate90(Ans[1] - Ans[0]) / l * h;
            Ans[3] = Ans[2] + (Ans[1] - Ans[0]);
            swap(Ans[2], Ans[3]);
        }
    }
    printf("%.5lf\n", ans);
    
    int id = 0;
    for(int i = 1; i < 4; i ++)
    {
        if(dcmp(Ans[i].y - Ans[id].y) < 0 ||
           (dcmp(Ans[i].y - Ans[id].y) == 0 && dcmp(Ans[i].x - Ans[id].x) < 0))
            id = i;
    }
    
    for(int i = 0; i < 4; i ++)
    {
        db x = Ans[(i + id) % 4].x, y = Ans[(i + id) % 4].y;
        if(fabs(x) < 1e-5) x = 0; if(fabs(y) < 1e-5) y = 0;
        printf("%.5lf %.5lf\n", x, y);
    }
}

（后面的特判是因为据说这个题卡精度）

完整代码：

#include 
#define DEBUG puts("QAQ")
#define openFile(a) freopen(a".in","r",stdin),freopen(a".out","w",stdout)
#define NOSYNC ios::sync_with_stdio(false); cin.tie(0); cout.tie(0)
#define FOR(i,j,k)  for(int (i) = (j); (i) <= (k); ++ (i))
#define RFOR(i,j,k) for(int (i) = (j); (i) >= (k); -- (i))
#define For(i,j,k)  for(int (i) = (j); (i) < (k); ++ (i))
#define RFor(i,j,k) for(int (i) = (j); (i) > (k); -- (i))
#define SC(...) scanf(__VA_ARGS__)
#define PR(...) printf(__VA_ARGS__)
#define N 50005
using namespace std;
typedef double db;
const double eps = 1e-10;
const double PI = acos(-1.0);
int dcmp(db x)
{
    if(fabs(x) < eps) return 0;
    return x > 0 ? 1 : -1;    
}

struct Vector {
    db x, y;
    Vector(db x = 0, db y = 0) : x(x), y(y) {}
    Vector operator + (const Vector &A) const {
        return Vector(x + A.x, y + A.y);
    }
    Vector operator - (const Vector &A) const {
        return Vector(x - A.x, y - A.y);
    }
    Vector operator * (const db &A) const {
        return Vector(x * A, y * A);
    }
    Vector operator / (const db &A) const {
        return Vector(x / A, y / A);
    }
    bool operator == (const Vector &A) const {
        return dcmp(x - A.x) == 0 && dcmp(y - A.y) == 0;
    }
};

typedef Vector Point;

db Dot(Vector A, Vector B) {
    return A.x * B.x + A.y * B.y;
}

db Cross(Vector A, Vector B) {
    return A.x * B.y - A.y * B.x;
}

db Length(Vector A) {
    return sqrt(Dot(A, A));
}

Vector Rotate90(Vector A) {
    return A = Vector(-A.y, A.x);
}

inline bool operator < (const Point &A, const Point &B) {
	return A.x == B.x ? A.y < B.y : A.x < B.x;
}

int n; Point p[N];
int stk[N], top;
bool used[N];
Point h[N]; int num;
db ans = 2e10;

void Andrew()
{
    top = 0;
    sort(p + 1,p + 1 + n);
    stk[++ top] = 1;
    for(int i = 2; i <= n; i ++)
    {
        while(top >= 2 && 
              Cross(p[stk[top]] - p[stk[top - 1]], p[i] - p[stk[top]]) <= 0)
            used[stk[top --]] = 0;
        used[i] = 1;
        stk[++ top] = i;
    }
    int size = top;
    for(int i = n - 1; i > 0; i --)
    {
        if(!used[i])
        {
            while(top > size && 
                  Cross(p[stk[top]] - p[stk[top - 1]], p[i] - p[stk[top]]) <= 0)
                used[stk[top --]] = 0;
            used[i] = 1;
            stk[++ top] = i;
        }
    }
    for(int i = 1; i <= top; i ++)
        h[i] = p[stk[i]];
    num = top - 1;
}

struct node {
    int j, l, r;
}t[N];

void Rotate1()
{
    int j = 3;
    for(int i = 1; i <= num; i ++)
    {
        while(Cross(h[i + 1] - h[i], h[j] - h[i + 1]) <= 
              Cross(h[i + 1] - h[i], h[j % num + 1] - h[i + 1]))
            j = j % num + 1;  
        t[i].j = j;    
    }
}

void Rotate2()
{
    int r = 2;
    for(int i = 1; i <= num; i ++)
    {
        while(Dot(h[i + 1] - h[i], h[r % num + 1] - h[r]) > 0)
            r = r + 1;
        t[i].r = r;
    }
}

void Rotate3()
{
    int l = num;
    for(int i = num; i; i --)
    {
        int pre = l - 1; if(!pre) pre = num;
        while(Dot(h[i + 1] - h[i], h[pre] - h[l]) < 0)
        {
            l = pre;
            pre = l - 1; if(!pre) pre = num;
        }
        t[i].l = l;
    }
}

Point Ans[4];
void work()
{
    Rotate1(); Rotate2(); Rotate3();
    for(int i = 1; i <= num; i ++)
    {
        Point A = h[i]; Point B = h[i + 1];
        Vector AB = h[i + 1] - h[i];
        Point C = h[t[i].l];
        Point D = h[t[i].r];
        Point E = h[t[i].j];
        db lenAB = Length(AB);
        db S = fabs(Cross(E - A, B - A));
        db l = fabs(Dot(D - B, B - A) / lenAB) + fabs(Dot(C - A, B - A) / lenAB) + lenAB;
        // r ~ i + 1 + l ~ i + i ~ i + 1
        db h = S / lenAB;
        if(l * h < ans)
        {
            ans = l * h;
            Ans[0] = A + AB * (Dot(C - A, B - A) / lenAB) / lenAB;
            Ans[1] = B + AB * (Dot(D - B, B - A) / lenAB) / lenAB;
            Ans[2] = Ans[0] + Rotate90(Ans[1] - Ans[0]) / l * h;
            Ans[3] = Ans[2] + (Ans[1] - Ans[0]);
            swap(Ans[2], Ans[3]);
        }
    }
    printf("%.5lf\n", ans);
    
    int id = 0;
    for(int i = 1; i < 4; i ++)
    {
        if(dcmp(Ans[i].y - Ans[id].y) < 0 ||
           (dcmp(Ans[i].y - Ans[id].y) == 0 && dcmp(Ans[i].x - Ans[id].x) < 0))
            id = i;
    }
    
    for(int i = 0; i < 4; i ++)
    {
        db x = Ans[(i + id) % 4].x, y = Ans[(i + id) % 4].y;
        if(fabs(x) < 1e-5) x = 0; if(fabs(y) < 1e-5) y = 0;
        printf("%.5lf %.5lf\n", x, y);
    }
}

int main()
{
	SC("%d", &n);
	FOR(i,1,n) SC("%lf%lf", &p[i].x, &p[i].y);
	Andrew();
	work();
	return 0;
}