[学习笔记]二维凸包

对于一个集合D，D中任意有限个点的凸组合的全体称为D的凸包。

比如下图中组成的图形就是一个凸包

 

 可以发现，二维凸包就像一个橡皮圈把所有点包了起来

我们可以简单地定义：平面凸包为平面上覆盖 $n$ 个点的最小凸多边形

 那么现在我们讨论求凸包的方法

1. 斜率逼近法

这个比较容易想到，不妨设斜率为 $k$

那么首先在所有点中找到最底部的（$y$ 最小的）点 $P_0$

接着从 $P_0$ 出发，逆时针方向寻找 $k$ 最小的点

如果有斜率相同的点，选取最远的加入凸包

依次处理，直到 $P_0$ 再次被扫描到为止

时间复杂度最坏是 $\Theta(n^2)$ 的，但当  $k \to \infty$ 也就是直线平行于 $y$ 轴时无法得到解，所以不常用

1.5 前置知识----叉积

叉积貌似是三位向量特有的运算？

对于两个向量 $\vec{a}=(a_x,a_y,a_z),\vec{b}=(b_x,b_y,b_z)$ ，我们有

$$a\times b=\begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ a_x & a_y & a_z \\ b_x & b_y & b_z \end{vmatrix}=\vec{i}\begin{vmatrix}a_y &a_z \\ b_y & b_z \end{vmatrix}+\vec{j}\begin{vmatrix}a_z &a_x \\ b_z & b_x \end{vmatrix}+\vec{k}\begin{vmatrix}a_x &a_y \\ b_x & b_y \end{vmatrix}=(a_yb_z-a_zb_y,a_zb_x-a_xb_z,a_xb_y-a_yb_x)$$

其中 $\vec{i},\vec{j},\vec{k}$ 分别为空间中的三个单位正交向量

叉积的方向是：伸出右手，从起始向量旋转至结束向量（这里是从 $\vec{a}$ 转至 $\vec{b}$），保持大拇指垂直于这两个向量所在平面，大拇指的方向就是向量的方向

其实这玩意立体几何用的也不少

2. Jarvis步进法

思路和斜率逼近法比较类似，总体可分为 $3$ 步

$\qquad\mathfrak{1:}$ 先找到 $y$ 最小的点 $P_0$ 入栈，遍历点集 $S$，找到使直线 $P_0P_1$ 与水平方向夹角最小的点 $P1$ 入栈

$\qquad\mathfrak{2:}$ 遍历所有点，找到与栈顶两个点的连线夹角最小的点 $P_i$ 入栈，继续操作直到 $P_i=P_1$

$\qquad\mathfrak{3:}$ 栈内的点就是构成凸包的点集

考虑如何找到最小夹角，我们有一个叫做叉积的神必玩意

不妨设栈顶点为 $P$ ，当前最优点为 $A$ 那么对于考虑到的点 $B$，可以用 $\overrightarrow{PA}\times\overrightarrow{PB}$ 进行判断

由于两个向量第三维是 $0$，因此它们的叉积只在第三维上有长度，即为 $a_xb_y-a_yb_x$，如果得出叉积大于 $0$ 那么这是逆时针旋转，肯定选取 $B$ 更优

如果叉积等于 $0$ 两向量共线，选取最远的点

 

#include<cstdio>
#include<cstring>
#include<string>
#define WR WinterRain
#define int long long
using namespace std;
const int WR=1001000,INF=1099588621776;
const double eps=1e-10;
struct Point{
    double x,y;
    Point make_point(double _x,double _y){
        Point res;
        res.x=_x;res.y=_y;
        return res;
    }
    Point operator-(Point &b){
        return make_point(x-b.x,y-b.y);
    }
}a[WR];
int n,stk[WR];
int read(){
    int s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        s=(s<<1)+(s<<3)+ch-48;
        ch=getchar();
    }
    return s*w;
}
double cross(Point a,Point b){
    return a.x*b.y-a.y*b.x;
}
double dis2(Point a,Point b){
    return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}
bool judge(Point p0,Point a,Point b){
    double m=cross(a-p0,b-p0);
    return (m<0||(m==0&&dis2(a,p0)>=dis2(b,p0)));
}
signed main(){
    n=read();
    int frst=1;
    for(int i=1;i<=n;i++){
        scanf("%lf%lf",&a[i].x,&a[i].y);
        if(a[i].y<a[frst].y||(a[i].y==a[frst].y&&a[i].x<a[frst].x)) frst=i;
    }
    int top=0;
    stk[++top]=frst;
    do{
        int tmp=1;
        for(int i=2;i<=n;i++){
            if(judge(a[stk[top]],a[i],a[tmp])) tmp=i;
        }
        stk[++top]=tmp;
    }while(stk[top]!=frst);
    top--;
    printf("%lld\n",top);
    for(int i=1;i<=top;i++){
        printf("%.5lf %.5lf\n",a[stk[i]].x,a[stk[i]].y);
    }
    return 0;
}

 

然后设凸包点集为 $\Re$ 这玩意的时间复杂度是 $\Theta(n\left\vert\Re\right\vert)$ 的显然不可行，因此我们需要更优秀的做法

 

2.5 前置知识----极角排序

首先引入极坐标系的概念：

在平面上取定一点 $O$ ，称为极点。从 $O$ 出发引一条射线 $Ox$ ，称为极轴。再取定一个单位长度，通常规定角度取逆时针方向为正。这样，平面上任一点 $P$ 的位置就可以用线段 $OP$ 的长度 $\rho$ 以及从 $Ox$ 到 $OP$ 的角度 $\theta$ 来确定，有序数对 $(\rho,\theta)$ 就称为 $P$ 点的极坐标，记为 $P(\rho,\theta)$ ；$\rho$ 称为 $P$ 点的极径，$\theta$ 称为 $P$ 点的极角。

引自百度百科

可以用叉积的方法，类似于 $\operatorname{Jarvis}$ 步进法

首先找到 $y$ 值最小的点作为原点，然后把所有点拉到第 $1,2$ 象限进行计算

 

#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#include<algorithm>
#define WR WinterRain
#define int long long
using namespace std;
const int WR=1001000,INF=1099588621776;
const double eps=1e-10;
struct Point{
    double x,y;
    Point make_point(double _x,double _y){
        Point res;
        res.x=_x,res.y=_y;
        return res;
    }
    Point operator-(Point&b){
        return make_point(x-b.x,y-b.y);
    }
}a[WR],p0;
int n;
int frst;
int read(){
    int s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        s=(s<<1)+(s<<3)+ch-48;
        ch=getchar();
    }
    return s*w;
}
double cross(Point a,Point b){
    return a.x*b.y-a.y*b.x;
}
double dis2(Point a,Point b){
    return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}
bool judge(Point a,Point b){
    double m=cross(a-p0,b-p0);
    return (m>0||(m==0&&dis2(a,p0)<=dis2(b,p0)));
}
int Quadrant(Point a){
    if(a.x>0&&a.y>=0) return 1;
    if(a.x<=0&&a.y>0) return 2;
    if(a.x<0&&a.y<=0) return 3;
    if(a.x>=0&&a.y<0) return 4;
    return -1;
}
bool cmp(Point a,Point b){
    if(Quadrant(a)==Quadrant(b)) return judge(a,b);
    return Quadrant(a)<Quadrant(b);
}
signed main(){
    n=read();
    frst=1;
    for(int i=1;i<=n;i++){
        scanf("%lf%lf",&a[i].x,&a[i].y);
        if(a[i].y<a[frst].y||(a[i].y==a[frst].y&&a[i].x<a[frst].x)) frst=i;
    }
    Point tmp=a[frst];
    for(int i=1;i<=n;i++) a[i]=a[i]-tmp;
    p0=a[frst].make_point(0,0);
    sort(a+1,a+1+n,cmp);
    for(int i=1;i<=n;i++) printf("%.3lf %.3lf\n",a[i].x,a[i].y);
    return 0;
}

 

还有一个黑科技：$\operatorname{atan2}$ 函数

它有两个参数：第一个代表已知点的 $y$ 坐标，第二个是已知点的 $x$ 坐标，返回的是这个点与原点的连线和 $x$ 轴的夹角

 

#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#include<algorithm>
#define WR WinterRain
#define int long long
using namespace std;
const int WR=1001000,INF=1099588621776;
const double eps=1e-10;
struct Point{
    double x,y;
    double angle;
    Point make_point(double _x,double _y){
        Point res;
        res.x=_x,res.y=_y;
        return res;
    }
    Point operator-(Point&b){
        return make_point(x-b.x,y-b.y);
    }
}a[WR],p0;
int n;
int frst;
int read(){
    int s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        s=(s<<1)+(s<<3)+ch-48;
        ch=getchar();
    }
    return s*w;
}
bool cmp(Point a,Point b){
    if(a.angle==b.angle) return (a.x*a.x+a.y*a.y)<(b.x*b.x+b.y*b.y);
    return a.angle<b.angle;
}
signed main(){
    n=read();
    frst=1;
    for(int i=1;i<=n;i++){
        scanf("%lf%lf",&a[i].x,&a[i].y);
        if(a[i].y<a[frst].y||(a[i].y==a[frst].y&&a[i].x<a[frst].x)) frst=i;
    }
    Point tmp=a[frst];
    for(int i=1;i<=n;i++){
        a[i]=a[i]-tmp;
        a[i].angle=atan2(a[i].y,a[i].x);
    }
    sort(a+1,a+1+n,cmp);
    for(int i=1;i<=n;i++) printf("%.3lf %.3lf\n",a[i].x,a[i].y);
    return 0;
}

 

3. Graham扫描法

结合下图进行理解

 

 

首先把所有点放在二维坐标系中，则纵坐标最小的点一定是凸包上的点，如图中的 $P_0$

然后把所有点的坐标平移一下，使 $P_0$ 作为原点，如上图

计算各个点相对于 $P_0$ 的幅角 $\alpha$ ，按从小到大的顺序对各个点排序。当 $\alpha$ 相同时，距离 $P_0$ 比较近的排在前面。

例如上图得到的结果为 $P_1$，$P_2$，$P_3$，$P_4$，$P_5$，$P_6$，$P_7$，$P_8$

显然的，第一个点 $P_1$ 和最后一个点 $P_8$ 一定是凸包上的点。

以上，我们已经知道了凸包上的第一个点  $P_0$ 和第二个点 $P_1$，我们把它们放在栈里面

现在把 $P_1$ 后面的那个点拿出来做当前点，即 $P_2$

接下来开始找第三个点：

$\qquad\mathfrak{1:}$ 连接 $P_0$ 和栈顶的那个点，得到直线 $l$

$\qquad\mathfrak{2:}$ 看当前点是在直线 $l$ 的右边还是左边。如果在直线的右边就执行步骤 $3$ ；如果在直线上，或者在直线的左边就执行步骤 $4$。

$\qquad\mathfrak{3:}$ 如果在右边，则栈顶的那个元素不是凸包上的点，把栈顶元素出栈。执行步骤 $2$。

$\qquad\mathfrak{4:}$ 当前点是凸包上的点，把它压入栈，执行步骤5。

$\qquad\mathfrak{5:}$ 检查当前的点 $P_2$ 是不是栈顶。是的话就结束，如果不是的话就把 $P_2$ 后面那个点做当前点，返回步骤 $4$。

可以看看动图

 

 

 

#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#include<algorithm>
#define WR WinterRain
#define int long long
using namespace std;
const int WR=1001000,INF=1099588621776;
const double eps=1e-10;
struct Point{
    double x,y;
    Point make_point(double _x,double _y){
        Point res;
        res.x=_x;res.y=_y;
        return res;
    }
    Point operator-(Point &b){
        return make_point(x-b.x,y-b.y);
    }
}a[WR],stk[WR];
int n,top;
Point p0;
int read(){
    int s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        s=(s<<1)+(s<<3)+ch-48;
        ch=getchar();
    }
    return s*w;
}
double cross(Point a,Point b){
    return a.x*b.y-a.y*b.x;
}
double dis2(Point a,Point b){
    return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}
bool cmp(Point a,Point b){
    double m=cross(a-p0,b-p0);
    return (m>0||(m==0&&dis2(a,p0)<=dis2(b,p0)));
}
signed main(){
    n=read();
    int frst=1;
    for(int i=1;i<=n;i++){
        scanf("%lf%lf",&a[i].x,&a[i].y);
        if(a[i].y<a[frst].y||(a[i].y==a[frst].y&&a[i].x<a[frst].x)) frst=i;
    }
    swap(a[1],a[frst]);p0=a[1];
    sort(a+2,a+1+n,cmp);
    stk[++top]=a[1];
    stk[++top]=a[2];
    for(int i=3;i<=n;i++){
        while(top>1&&cross(stk[top]-stk[top-1],a[i]-stk[top])<0) top--;
        stk[++top]=a[i];
    }
    printf("%lld\n",top);
    for(int i=1;i<=top;i++) printf("%lf %lf\n",stk[i].x,stk[i].y);
    return 0;
}

 

4. Andrew算法求凸包

$\operatorname{Andrew}$ 算法类似于 $\operatorname{Graham}$ 算法的改进，由于省去了 $\operatorname{sort}$ 中的叉积部分复杂度更优

主体思路是一样的，只不过在排序的时候按照 $x,y$ 的顺序排序

但是这样扫描的时候我们会惊讶的发现，

凸包只扫了一半就到头了$\cdots\cdots$

不过不用慌，我们

从 $n$ 向前再扫一遍就行$\cdots\cdots$

 

#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#include<algorithm>
#define WR WinterRain
#define int long long
using namespace std;
const int WR=1001000,INF=1099588621776;
const double eps=1e-10;
struct Point{
    double x,y;
    Point make_point(double _x,double _y){
        Point res;
        res.x=_x;res.y=_y;
        return res;
    }
    Point operator-(Point b){
        return make_point(x-b.x,y-b.y);
    }
    bool operator<(const Point &b)const{
        if(x==b.x) return y<b.y;
        return x<b.x;
    }
}a[WR];
int n,stk[WR],top;
Point p0;
int read(){
    int s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        s=(s<<1)+(s<<3)+ch-48;
        ch=getchar();
    }
    return s*w;
}
double cross(Point a,Point b){
    return a.x*b.y-a.y*b.x;
}
double dis2(Point a,Point b){
    return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}
signed main(){
    n=read();
    int frst=1;
    for(int i=1;i<=n;i++){
        scanf("%lf%lf",&a[i].x,&a[i].y);
    }
    sort(a+1,a+1+n);
    stk[++top]=1;
    stk[++top]=2;
    for(int i=3;i<=n;i++){
        while(top>1&&cross(a[stk[top]]-a[stk[top-1]],a[i]-a[stk[top]])<0) top--;
        stk[++top]=i;
    }
    stk[++top]=n-1;
    for(int i=n-2;i>=1;i--){
        while(top>1&&cross(a[stk[top]]-a[stk[top-1]],a[i]-a[stk[top]])<0) top--;
        stk[++top]=i;
    }
    printf("%lld\n",top);
    for(int i=1;i<=top;i++) printf("%lf %lf\n",a[stk[i]].x,a[stk[i]].y);
    return 0;
}

 

顺便可以求出凸包的周长

 

#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#include<algorithm>
#define WR WinterRain
#define int long long
using namespace std;
const int WR=1001000,INF=1099588621776;
const double eps=1e-10;
struct Point{
    double x,y;
    Point make_point(double _x,double _y){
        Point res;
        res.x=_x;res.y=_y;
        return res;
    }
    Point operator-(Point &b){
        return make_point(x-b.x,y-b.y);
    }
}a[WR],stk[WR];
int n,top;
Point p0;
int read(){
    int s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        s=(s<<1)+(s<<3)+ch-48;
        ch=getchar();
    }
    return s*w;
}
double cross(Point a,Point b){
    return a.x*b.y-a.y*b.x;
}
double dis2(Point a,Point b){
    return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}
bool cmp(Point a,Point b){
    double m=cross(a-p0,b-p0);
    return (m>0||(m==0&&dis2(a,p0)<=dis2(b,p0)));
}
signed main(){
    n=read();
    int frst=1;
    for(int i=1;i<=n;i++){
        scanf("%lf%lf",&a[i].x,&a[i].y);
        if(a[i].y<a[frst].y||(a[i].y==a[frst].y&&a[i].x<a[frst].x)) frst=i;
    }
    swap(a[1],a[frst]);p0=a[1];
    sort(a+2,a+1+n,cmp);
    stk[++top]=a[1];
    stk[++top]=a[2];
    for(int i=3;i<=n;i++){
        while(top>1&&cross(stk[top]-stk[top-1],a[i]-stk[top])<0) top--;
        stk[++top]=a[i];
    }
    // printf("%lld\n",top);
    // for(int i=1;i<=top;i++) printf("%lf %lf\n",stk[i].x,stk[i].y);
    stk[++top]=a[1];
    double ans=0.0;
    for(int i=1;i<top;i++) ans+=sqrt(dis2(stk[i],stk[i+1])); 
    printf("%.2lf\n",ans);
    return 0;
}

 

5. 更优的解法----Melkman算法

据说这个复杂度更好？

算法的策略非常直接，按原始顺序依次处理多段线上的每个点

假定输入多段线为 $S=\{P_0,P_1,\cdots,P_n\}$ 。在每一步，算法将当前为止处理过的所有顶点形成的凸包存储在一个双端队列中

接下来，考虑下一个顶点 $P_k$ 。$P_k$ 有两种情况：

$\qquad\mathfrak{1:}$ $P_k$在当前凸包内；

$\qquad\mathfrak{2:}$ $P_k$在当前凸包外，并且形成一个新的凸包顶点。原来凸包中的点可能变为在新凸包的内部，需要被丢弃，然后再将 $P_k$ 加入队列

首先给双端队列两个标签: $bot$ 和 $top$ ，在这中间的是当前凸包结果

在队列两端，都可以增加或删除元素。在顶部（ $top$ 之上），我们称为 $push / pop$ ；在底部（ $bot$ 之下），我们将增删元素的操作称为 $insert / delete$

不妨将队列记为 $\mathsf{\Phi}=\{D_{bot},\cdots,D_{top}$，$D_i$ 是原多段线中的点。

当 $D_{bot}=D_{top}$ 就形成了一个多边形。在 $\operatorname{Melkman}$ 算法中，处理顶点 $P_k$ 后，$\mathsf{\Phi}_k$ 满足：

$\qquad\mathfrak{1:}$ $$\mathsf{\Phi}_k$ 是多段线 $S=\{P_0,P_1,\cdots,P_k}$ 的逆时针方向的凸包；

$\qquad\mathfrak{2:}$ $D_{bot}=D_{top}$，是最近添加到 $\mathsf{\Phi}_k$ 中的点。

 然后代码没写，不想写了qwq

 

6. 动态凸包

那么结合例题来看

动态凸包就是在可增删直线的情况下，求凸包（或凸包上的点）。

回忆静态二维凸包的构建方法，我们先维护一个上下凸壳，再进行操作

实际上现在就是让这个维护变成动态的

C++ STL里的map会非常好用

开两个map，一个叫top，维护上凸壳；一个叫down，维护下凸壳。

top/down[x]表示横坐标为x的点的纵坐标是什么

这个过程比较类似于 $\operatorname{Melkman}$ 算法

 

#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#include<map>
#include<algorithm>
#define WR WinterRain
#define int long long
using namespace std;
const int WR=100100,INF=1099588621776;
int n,m;
map<int,int>top,bot;
int read(){
    int s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        s=(s<<1)+(s<<3)+ch-48;
        ch=getchar();
    }
    return s*w;
}
bool check_top(int x,int y){
    map<int,int>::iterator i=top.lower_bound(x);
    if(i==top.end()) return false;
    if(i->first==x) return y<=i->second;
    if(i==top.begin()) return 0;
    map<int,int>::iterator j=i;j--;
    return (i->first-j->first)*(y-j->second)-(i->second-j->second)*(x-j->first)<=0;
}
bool check_bot(int x,int y){
    map<int,int>::iterator i=bot.lower_bound(x);
    if(i==bot.end()) return false;
    if(i->first==x) return y>=i->second;
    if(i==bot.begin()) return 0;
    map<int,int>::iterator j=i;j--;
    return (i->first-j->first)*(y-j->second)-(i->second-j->second)*(x-j->first)>=0;
}
bool delete_top(map<int,int>::iterator i){
    if(i==top.begin()) return false;
    map<int,int>::iterator j=i,k=i;
    j--,k++;
    if(k==top.end()) return false;
    if((i->first-j->first)*(k->second-j->second)-(i->second-j->second)*(k->first-j->first)>=0){
        top.erase(i);
        return true;
    }else return false;
}
bool delete_bot(map<int,int>::iterator i){
    if(i==bot.begin()) return false;
    map<int,int>::iterator j=i,k=i;
    j--,k++;
    if(k==bot.end()) return false;
    if((i->first-j->first)*(k->second-j->second)-(i->second-j->second)*(k->first-j->first)<=0){
        bot.erase(i);
        return true;
    }else return false;
}
void insert_top(int x,int y){
    if(check_top(x,y)) return;
    top[x]=y;
    map<int,int>::iterator i=top.find(x),j=i;
    if(i!=top.begin()){
        j--;
        while(delete_top(j++)) j--;
    }
    if(++j!=top.end()) while(delete_top(j--)) j++;
}
void insert_bot(int x,int y){
    if(check_bot(x,y)) return;
    bot[x]=y;
    map<int,int>::iterator i=bot.find(x),j=i;
    if(i!=bot.begin()){
        j--;
        while(delete_bot(j++)) j--;
    }
    if(++j!=bot.end()) while(delete_bot(j--)) j++;
}
signed main(){
    n=read();
    for(int i=1;i<=n;i++){
        int opt=read();
        if(opt==1){
            int x=read(),y=read();
            insert_top(x,y);
            insert_bot(x,y);
        }else{
            int x=read(),y=read();
            if(check_top(x,y)&&check_bot(x,y)) printf("YES\n");
            else printf("NO\n");
        }
    }
    return 0;
}