凸包算法

BlueHeart0621 · 简书 · · 2020-03-07 00:15

凸包类型的题算法主要有三种： JarvisMarch 算法、 Graham 算法和 Andrew 算法，这三种算法时间性能上递增。

1. JarvisMarch 算法

1.1 思想

纵坐标最小然后横坐标最小的点一定是凸包上的点，我们将其记为 ${p_0}$ ，从 ${p_0}$ 开始，按逆时针的方向，逐个找凸包上的点，每前进一步找到一个点，所以叫作步进法。

选取下一个点的方法：
假设已找到 ${p_0}$ 、 ${p_1}$ ，则利用跟 ${p_0p_1}$ 向量夹角最小的点作为 ${p_2}$ 。（ ${p_1}$ 则利用 ${p_0p_1}$ 向量和水平线的夹角）

1.2 代码

/******************************************************************
                    Jarvis March的步进算法
算法复杂度：O(nH)。（其中 n 是点的总个数，H 是凸包上的点的个数）
******************************************************************/

//#include <bits/stdc++.h>
#include <queue>
#include <cstdio>
#include <cmath>
#include <algorithm>
 
using namespace std;
 
const int MAXN = 10005;
const double MAXD = 1e9;
const double ACCUR = 1e-9;
 
struct node
{
    double x, y;        //点的坐标
    bool operator < (const node n) const
    {
        if (abs(y-n.y) < ACCUR) return x < n.x;
        else    return y < n.y;
    }
    bool operator == (const node n) const
    {
        return (abs(x-n.x) < ACCUR) && (abs(y-n.y) < ACCUR);
    }
    void operator = (const node n)
    {
        x = n.x;
        y = n.y;
    }
};
 
struct vect
{
    double x, y;
    void operator = (const vect v)
    {
        x = v.x;
        y = v.y;
    }
    double operator *(const vect v) const
    {
        return x*v.x + y*v.y;
    }
};
 
bool equal (const double d1, const double d2)
{
    return abs(d1-d2) < ACCUR;
}
vect vform(const node n1, const node n2)
{
    vect tmpv;
    tmpv.x = n2.x - n1.x;
    tmpv.y = n2.y - n1.y;
    return tmpv;
}
double vlen(const vect v)
{
    return sqrt(v.x*v.x+v.y*v.y);
}
double vcos(const vect v1, const vect v2)
{
    return (v1*v2)/(vlen(v1)*vlen(v2));
}
double area (const node n1, const node n2, const node n3)
{
    double b1, b2, b3;
    b1 = vlen(vform(n1,n2));
    b2 = vlen(vform(n2,n3));
    b3 = vlen(vform(n3,n1));
    double b = (b1+b2+b3)/2;
    return sqrt(b*(b-b1)*(b-b2)*(b-b3));
}
 
node p[MAXN];           //点集
queue <node> bq;        //凸包顶点集
 
int main()
{
    int n;
    while(scanf("%d", &n) == 1)
    {
        if(n == 0)
        {
            break;
        }
        /*【注意】第一个点先不标记，作为循环结束条件（即最后找到第一个点）*/
        int f[MAXN] = {0};      //点集标记数组;
        vect v;                 //v表示上两个点形成的向量。
        node p0, p1;             //p0表示第一个点，p1表示上一个点。
        p0.x = p0.y = MAXD;     //初始化
 
        for (int i = 0; i < n; ++i)
        {
            scanf("%lf%lf", &(p[i].x), &(p[i].y));
            if (p[i] < p0)
            {
                p0 = p[i];
            }
        }
 
        p1 = p0;                //初始化上一个点
        //【注意】初始化向量的选取跟第一个点的选取关。
        //如果第一个点是横坐标最小然后纵坐标最小则初始向量为竖直单位向量
        v.x = 1; v.y = 0;       //初始向量为水平单位向量。
        do
        {
            node p2;            //待判定的点。
            vect v1;            //待判定的向量
            int j;              //带判定的点的下标
            double minvcos = -1, minvlen = MAXD;    //初始最大夹角和最小向量长度。
            for (int i = 0; i < n; ++i)
            {
                if (!f[i])      //判断该点是否已经在凸包上
                {
                    vect tmpv;
                    tmpv.x = p[i].x-p1.x;
                    tmpv.y = p[i].y-p1.y;
                    if (vcos(v,tmpv) > minvcos)
                    {
                        p2 = p[i];
                        v1 = tmpv;
                        j = i;
                        minvcos = vcos(v,tmpv);
                        minvlen = vlen(tmpv);
                    }
                    else if (equal(vcos(v,tmpv),minvcos) && vlen(tmpv) < minvlen)
                    {
                        p2 = p[i];
                        v1 = tmpv;
                        j = i;
                        minvcos = vcos(v,tmpv);
                        minvlen = vlen(tmpv);
                    }
                }
            }
            bq.push(p2);
            p1 = p2;
            v = v1;
            f[j] = 1;
            //printf("minvcos=%f,minvlen=%f\n", minvcos, minvlen);
        }while(!(p1==p0));
 
        /*
        while(!bq.empty())
        {
            node tmpp = bq.front();
            printf("(%f,%f)\n", tmpp.x, tmpp.y);
            bq.pop();
        }
        */
        
        //凸包周长
        double ans = 0;
        node fp, ep;
        fp = p0;
        while(!bq.empty())
        {
            ep = bq.front();
            bq.pop();
            ans += vlen(vform(fp, ep));
            fp = ep;
        }
        printf("%.2f\n", ans);
     
        /*
        //凸包面积
        double ans = 0;
        node fp = bq.front();
        bq.pop();
        node np = bq.front();
        bq.pop();
        while(!bq.empty())
        {
            node ep = bq.front();
            bq.pop();
            ans += area(fp,np,ep);
            np = ep;
            //printf("(%f,%f)\n", tmpp.x, tmpp.y);
        }
        printf("%d\n", (int)ans/50);
        */
    }
    return 0;
}

2. Graham 算法

2.1 思想

把所有点放在二维坐标系中，则纵坐标最小的点一定是凸包上的点，记为 ${p_0}$ 。计算各个点相对 ${p_0}$ 的幅角，按从小到大的顺序对各个点排序。（当幅角相同是，距离 ${p_0}$ 比较近的排在前面）则幅角最小的点和最大的点一定在凸包上。取幅角最小的点记为 ${p_1}$ ，将 ${p_0}$ 、 ${p_1}$ 入栈。连接栈顶的点和次栈顶的点，得到直线 $l$ ，看当前点是在直线的右边还是左边，在右边则栈顶元素不是凸包上的点，将其弹出，返回继续执行。如果在左边，则当前点是凸包上的点。一直到幅角最大的那个点为之。

叉积原理
两个向量的叉积 ${P_1 \times P_2 = x_1y_2 - x_2y_1}$ ，其中用结果的正负代表叉乘结果的方向。该公式本质是两个三维向量（ $z$ 轴分量为0）叉乘的结果（原来结果为 ${(x_1y_2 - x_2y_1)\cdot \vec{k}}$ ，其中 ${\vec{k}}$ 是 $z$ 轴单位正向量）。

2.2 代码

/************************************************************************
                        Graham Scan算法
时间复杂度：O(nlogn)。Scan过程为O(n)，预处理排序为O(nlogn)。
预处理排序：极角排序。
************************************************************************/

//#include <bits/stdc++.h>
#include <stack>
#include <cstdio>
#include <cmath>
#include <algorithm>
 
using namespace std;
 
const int MAXN = 10005;
const double MAXD = 1e9;
const double ACCUR = 1e-9;
 
struct node
{
    double x, y;        //点的坐标
    bool operator < (const node n) const
    {
        if (abs(y-n.y) < ACCUR) return x < n.x;
        else    return y < n.y;
    }
    bool operator == (const node n) const
    {
        return (abs(x-n.x) < ACCUR) && (abs(y-n.y) < ACCUR);
    }
    void operator = (const node n)
    {
        x = n.x;
        y = n.y;
    }
};
 
struct vect
{
    double x, y;
    void operator = (const vect v)
    {
        x = v.x;
        y = v.y;
    }
    //叉积
    double operator *(const vect v) const
    {
        return x*v.y - y*v.x;
    }
};
 
node p0;            //纵坐标最小的点
 
bool equal (const double d1, const double d2)
{
    return abs(d1-d2) < ACCUR;
}
vect vform(const node n1, const node n2)
{
    vect tmpv;
    tmpv.x = n2.x - n1.x;
    tmpv.y = n2.y - n1.y;
    return tmpv;
}
double vlen(const vect v)
{
    return sqrt(v.x*v.x+v.y*v.y);
}
double vcos(const vect v1, const vect v2)
{
    return (v1*v2)/(vlen(v1)*vlen(v2));
}
//极角排序
bool cmpp (const node p1, const node p2)
{
    vect v1, v2;
    v1 = vform(p0, p1);
    v2 = vform(p0, p2);
    if (equal(v1*v2,0))
    {
        return vlen(v1) < vlen(v2);
    }
    else
    {
        return v1*v2 > 0;
    }
}
//叉积判断点（v2的终点）在v1左边还是右边
bool cmpv (const vect v1, const vect v2)
{
    return (v1*v2 > 0) || equal(v1*v2,0);
}
double area (const node n1, const node n2, const node n3)
{
    /*
    //海伦公式
    double b1, b2, b3;
    b1 = vlen(vform(n1,n2));
    b2 = vlen(vform(n2,n3));
    b3 = vlen(vform(n3,n1));
    double b = (b1+b2+b3)/2;
    return sqrt(b*(b-b1)*(b-b2)*(b-b3));
    */
 
    //叉积公式（叉积为平行四边形面积）
    vect v1, v2;
    v1 = vform(n1, n2);
    v2 = vform(n1, n3);
    return abs(v1*v2)/2;
}
 
node p[MAXN];           //点集
stack <node> bs;        //凸包顶点集
 
int main()
{
    int n;
    p0.x = p0.y = MAXD; //初始化第一个点
    scanf("%d", &n);
    for(int i = 0; i < n; ++i)
    {
        scanf("%lf%lf", &(p[i].x), &(p[i].y));
        if(p[i] < p0)
        {
            p0 = p[i];
        }
    }
    sort(p,p+n,cmpp);
    bs.push(p[0]);
    bs.push(p[1]);
    int j = 2;
    while(j < n)
    {
        //取出栈顶和次栈顶
        node p1, p2;
        p2 = bs.top();
        bs.pop();
        p1 = bs.top();
        //构造叉乘向量
        vect v1, v2;
        v1 = vform(p1,p2);
        v2 = vform(p2,p[j]);
        if(cmpv(v1,v2))
        {
            bs.push(p2);
            bs.push(p[j]);
            ++j;
        }
    }
    /*
    while(!bs.empty())
    {
        node tmpp = bs.top();
        printf("(%f,%f)\n", tmpp.x, tmpp.y);
        bs.pop();
    }
    */
    /*
    //计算周长
    double ans = 0;
    node fp, ep;
    fp = p[0];
    while(!bs.empty())
    {
        ep = bs.top();
        bs.pop();
        ans += vlen(vform(fp, ep));
        fp = ep;
    }
    printf("%.2f\n", ans);
    */
 
    //计算面积
    double ans = 0;
    node fp, np, ep;
    fp = bs.top();
    bs.pop();
    np = bs.top();
    bs.pop();
    while(!bs.empty())
    {
        ep = bs.top();
        bs.pop();
        ans += area(fp,np,ep);
        np = ep;
    }
    printf("%d\n", (int)ans/50);
    return 0;
}

3. Andrew 算法

3.1 思想

预处理排序改为水平排序，按照横坐标从小到大进行排序，横坐标相同则按纵坐标从小到大排。按照 graham 算法思想从 ${p_0}$ 、 ${p_1}$ 扫描所有点得到下凸包，再从 ${p_{n-1}}$ 、 ${p_{n-2}}$ 扫描所有点得到上凸包，二者结合即为整个凸包。（注意：这里的 ${p_1}$ 不一定在凸包里）

3.2 代码

/************************************************************************
                    Andrew算法（Graham Scan算法变种）
时间复杂度：O(nlogn)。Scan过程为O(n)，预处理排序为O(nlogn)。
预处理排序：水平排序排序。
************************************************************************/
 
//#include <bits/stdc++.h>
#include <stack>
#include <cstdio>
#include <cmath>
#include <algorithm>
 
using namespace std;
 
const int MAXN = 10005;
const double MAXD = 1e9;
const double ACCUR = 1e-9;
 
struct node
{
    double x, y;        //点的坐标
    //水平排序（与极角排序不一样，只能确定p0和pn-1在凸包内）
    bool operator < (const node n) const
    {
        if (abs(x-n.x) < ACCUR) return y < n.y;
        else    return x < n.x;
    }
    bool operator == (const node n) const
    {
        return (abs(x-n.x) < ACCUR) && (abs(y-n.y) < ACCUR);
    }
    void operator = (const node n)
    {
        x = n.x;
        y = n.y;
    }
};
 
struct vect
{
    double x, y;
    void operator = (const vect v)
    {
        x = v.x;
        y = v.y;
    }
    //叉积
    double operator *(const vect v) const
    {
        return x*v.y - y*v.x;
    }
};
 
bool equal (const double d1, const double d2)
{
    return abs(d1-d2) < ACCUR;
}
vect vform(const node n1, const node n2)
{
    vect tmpv;
    tmpv.x = n2.x - n1.x;
    tmpv.y = n2.y - n1.y;
    return tmpv;
}
//计算向量长度
double vlen(const vect v)
{
    return sqrt(v.x*v.x+v.y*v.y);
}
//计算向量夹角余弦值
double vcos(const vect v1, const vect v2)
{
    return (v1*v2)/(vlen(v1)*vlen(v2));
}
/*
//极角排序
bool cmpp (const node p1, const node p2)
{
    vect v1, v2;
    v1 = vform(p0, p1);
    v2 = vform(p0, p2);
    if (equal(v1*v2,0))
    {
        return vlen(v1) < vlen(v2);
    }
    else
    {
        return v1*v2 > 0;
    }
}
*/
//叉积判断点（v2的终点）在v1左边还是右边
bool cmpv (const vect v1, const vect v2)
{
    return (v1*v2 > 0) || equal(v1*v2,0);
}
double area (const node n1, const node n2, const node n3)
{
    /*
    //海伦公式
    double b1, b2, b3;
    b1 = vlen(vform(n1,n2));
    b2 = vlen(vform(n2,n3));
    b3 = vlen(vform(n3,n1));
    double b = (b1+b2+b3)/2;
    return sqrt(b*(b-b1)*(b-b2)*(b-b3));
    */
 
    //叉积公式（叉积为平行四边形面积）
    vect v1, v2;
    v1 = vform(n1, n2);
    v2 = vform(n1, n3);
    return abs(v1*v2)/2;
}
 
node p[MAXN];           //点集
stack <node> bs;        //凸包顶点集
 
int main()
{
    int n;
    scanf("%d", &n);
    for(int i = 0; i < n; ++i)
    {
        scanf("%lf%lf", &(p[i].x), &(p[i].y));
    }
    sort(p,p+n);
    //正向扫描（上凸包）
    bs.push(p[0]);
    bs.push(p[1]);
    int j = 2;
    while(j < n)
    {
        //取出栈顶和次栈顶
        node p1, p2;
        p2 = bs.top();
        bs.pop();
        //p1不一定在凸包中
        if(bs.empty())
        {
            bs.push(p2);
            bs.push(p[j]);
            ++j;
        }
        else
        {
            p1 = bs.top();
            //构造叉乘向量
            vect v1, v2;
            v1 = vform(p1,p2);
            v2 = vform(p2,p[j]);
            if(cmpv(v1,v2))
            {
                bs.push(p2);
                bs.push(p[j]);
                ++j;
            }
        }
    }
    //反向扫描（下凸包）
    int k = n-2;
    while(k >= 0)
    {
        //取出栈顶和次栈顶
        node p1, p2;
        p2 = bs.top();
        bs.pop();
        p1 = bs.top();
        //构造叉乘向量
        vect v1, v2;
        v1 = vform(p1,p2);
        v2 = vform(p2,p[k]);
        if(cmpv(v1,v2))
        {
            bs.push(p2);
            bs.push(p[k]);
            --k;
        }
    }
    bs.pop();   //p0重复进栈一次
 
    /*
    while(!bs.empty())
    {
        node tmpp = bs.top();
        printf("(%f,%f)\n", tmpp.x, tmpp.y);
        bs.pop();
    }
    */
    /*
    //计算周长
    double ans = 0;
    node fp, ep;
    fp = p[0];
    while(!bs.empty())
    {
        ep = bs.top();
        bs.pop();
        ans += vlen(vform(fp, ep));
        fp = ep;
    }
    printf("%.2f\n", ans);
    */
    
    //计算面积
    double ans = 0;
    node fp, np, ep;
    fp = bs.top();
    bs.pop();
    np = bs.top();
    bs.pop();
    while(!bs.empty())
    {
        ep = bs.top();
        bs.pop();
        ans += area(fp,np,ep);
        np = ep;
    }
    printf("%d\n", (int)ans/50);
 
    return 0;
}

原文地址：访问原文地址
快照地址：访问文章快照

分享到微博