C++ 利用分治算法求凸包(Convex Hull using Divide and Conquer Algorithm)

如果您喜欢此文章，请收藏、点赞、评论，谢谢，祝您快乐每一天。

在计算几何中，凸包是包含给定点集的最小凸多边形。它是一个基本概念，在计算机图形学、机器人学和图像处理等各个领域都有应用。

凸包的重要性：

凸包在计算几何中非常重要，原因有以下几点：

• 碰撞检测：凸包可用于有效检测 2D 或 3D 空间中物体之间的碰撞。
• 图像处理：凸包可用于从图像中提取有意义的形状，例如物体的轮廓。
• 数据可视化：凸包可用于可视化散点图中数据点的分布。
• 输入是一个由 x 和 y 坐标指定的点数组。输出是这组点的凸包。

利用分治算法求凸包：

先决条件：两个凸多边形之间的切线

JavaScript 两个凸多边形之间的切线:https://blog.csdn.net/hefeng_aspnet/article/details/160182506
C# 两个凸多边形之间的切线:https://blog.csdn.net/hefeng_aspnet/article/details/160182477
Python 两个凸多边形之间的切线:https://blog.csdn.net/hefeng_aspnet/article/details/160182449
Java 两个凸多边形之间的切线:https://blog.csdn.net/hefeng_aspnet/article/details/160182407
C++ 两个凸多边形之间的切线:https://blog.csdn.net/hefeng_aspnet/article/details/160182260

算法：给定一个点集，我们需要找到它的凸包。假设我们已知左半部分点集和右半部分点集的凸包，那么现在的问题是将这两个凸包合并，从而确定整个点集的凸包。这可以通过找到左右凸包的上切线和下切线来实现。如图所示，两个凸多边形之间的切线。设左凸包为 a，右凸包为 b。则下切线和上切线分别标记为 1 和 2，如图所示。红色轮廓线表示最终的凸包。

现在的问题是如何找到左右两半的凸包。这时递归就派上用场了。我们将点集分割，直到点数非常少，比如只有 5 个点，然后我们可以通过穷举算法找到这些点的凸包。将这两半合并，就能得到整个点集的凸包。

注：我们使用暴力算法来寻找少量点的凸包，其时间复杂度为 O(n³)。但有人提出，由三个或更少点构成的凸包就是全部点集。这没错，但问题在于，当我们尝试合并一个由两个点构成的左凸包和一个由三个点构成的右凸包时，程序在某些特殊情况下会陷入无限循环。因此，为了解决这个问题，我直接求出了由五个或更少点构成的凸包。O(n³)算法虽然略微复杂一些，但并不影响算法的整体复杂度。 

以下是具体实现步骤：

// A divide and conquer program to find convex
// hull of a given set of points.
#include<bits/stdc++.h>
using namespace std;

// stores the centre of polygon (It is made
// global because it is used in compare function)
pair<int, int> mid;

// determines the quadrant of a point
// (used in compare())
int quad(pair<int, int> p)
{
    if (p.first >= 0 && p.second >= 0)
        return 1;
    if (p.first <= 0 && p.second >= 0)
        return 2;
    if (p.first <= 0 && p.second <= 0)
        return 3;
    return 4;
}

// Checks whether the line is crossing the polygon
int orientation(pair<int, int> a, pair<int, int> b,
                pair<int, int> c)
{
    int res = (b.second-a.second)*(c.first-b.first) -
            (c.second-b.second)*(b.first-a.first);

    if (res == 0)
        return 0;
    if (res > 0)
        return 1;
    return -1;
}

// compare function for sorting
bool compare(pair<int, int> p1, pair<int, int> q1)
{
    pair<int, int> p = make_pair(p1.first - mid.first,
                                p1.second - mid.second);
    pair<int, int> q = make_pair(q1.first - mid.first,
                                q1.second - mid.second);

    int one = quad(p);
    int two = quad(q);

    if (one != two)
        return (one < two);
    return (p.second*q.first < q.second*p.first);
}

// Finds upper tangent of two polygons 'a' and 'b'
// represented as two vectors.
vector<pair<int, int>> merger(vector<pair<int, int> > a,
                            vector<pair<int, int> > b)
{
    // n1 -> number of points in polygon a
    // n2 -> number of points in polygon b
    int n1 = a.size(), n2 = b.size();

    int ia = 0, ib = 0;
    for (int i=1; i<n1; i++)
        if (a[i].first > a[ia].first)
            ia = i;

    // ib -> leftmost point of b
    for (int i=1; i<n2; i++)
        if (b[i].first < b[ib].first)
            ib=i;

    // finding the upper tangent
    int inda = ia, indb = ib;
    bool done = 0;
    while (!done)
    {
        done = 1;
        while (orientation(b[indb], a[inda], a[(inda+1)%n1]) >=0)
            inda = (inda + 1) % n1;

        while (orientation(a[inda], b[indb], b[(n2+indb-1)%n2]) <=0)
        {
            indb = (n2+indb-1)%n2;
            done = 0;
        }
    }

    int uppera = inda, upperb = indb;
    inda = ia, indb=ib;
    done = 0;
    int g = 0;
    while (!done)//finding the lower tangent
    {
        done = 1;
        while (orientation(a[inda], b[indb], b[(indb+1)%n2])>=0)
            indb=(indb+1)%n2;

        while (orientation(b[indb], a[inda], a[(n1+inda-1)%n1])<=0)
        {
            inda=(n1+inda-1)%n1;
            done=0;
        }
    }

    int lowera = inda, lowerb = indb;
    vector<pair<int, int>> ret;

    //ret contains the convex hull after merging the two convex hulls
    //with the points sorted in anti-clockwise order
    int ind = uppera;
    ret.push_back(a[uppera]);
    while (ind != lowera)
    {
        ind = (ind+1)%n1;
        ret.push_back(a[ind]);
    }

    ind = lowerb;
    ret.push_back(b[lowerb]);
    while (ind != upperb)
    {
        ind = (ind+1)%n2;
        ret.push_back(b[ind]);
    }
    return ret;

}

// Brute force algorithm to find convex hull for a set
// of less than 6 points
vector<pair<int, int>> bruteHull(vector<pair<int, int>> a)
{
    // Take any pair of points from the set and check
    // whether it is the edge of the convex hull or not.
    // if all the remaining points are on the same side
    // of the line then the line is the edge of convex
    // hull otherwise not
    set<pair<int, int> >s;

    for (int i=0; i<a.size(); i++)
    {
        for (int j=i+1; j<a.size(); j++)
        {
            int x1 = a[i].first, x2 = a[j].first;
            int y1 = a[i].second, y2 = a[j].second;

            int a1 = y1-y2;
            int b1 = x2-x1;
            int c1 = x1*y2-y1*x2;
            int pos = 0, neg = 0;
            for (int k=0; k<a.size(); k++)
            {
                if (a1*a[k].first+b1*a[k].second+c1 <= 0)
                    neg++;
                if (a1*a[k].first+b1*a[k].second+c1 >= 0)
                    pos++;
            }
            if (pos == a.size() || neg == a.size())
            {
                s.insert(a[i]);
                s.insert(a[j]);
            }
        }
    }

    vector<pair<int, int>>ret;
    for (auto e:s)
        ret.push_back(e);

    // Sorting the points in the anti-clockwise order
    mid = {0, 0};
    int n = ret.size();
    for (int i=0; i<n; i++)
    {
        mid.first += ret[i].first;
        mid.second += ret[i].second;
        ret[i].first *= n;
        ret[i].second *= n;
    }
    sort(ret.begin(), ret.end(), compare);
    for (int i=0; i<n; i++)
        ret[i] = make_pair(ret[i].first/n, ret[i].second/n);

    return ret;
}

// Returns the convex hull for the given set of points
vector<pair<int, int>> divide(vector<pair<int, int>> a)
{
    // If the number of points is less than 6 then the
    // function uses the brute algorithm to find the
    // convex hull
    if (a.size() <= 5)
        return bruteHull(a);

    // left contains the left half points
    // right contains the right half points
    vector<pair<int, int>>left, right;
    for (int i=0; i<a.size()/2; i++)
        left.push_back(a[i]);
    for (int i=a.size()/2; i<a.size(); i++)
        right.push_back(a[i]);

    // convex hull for the left and right sets
    vector<pair<int, int>>left_hull = divide(left);
    vector<pair<int, int>>right_hull = divide(right);

    // merging the convex hulls
    return merger(left_hull, right_hull);
}

// Driver code
int main()
{
    vector<pair<int, int> > a;
    a.push_back(make_pair(0, 0));
    a.push_back(make_pair(1, -4));
    a.push_back(make_pair(-1, -5));
    a.push_back(make_pair(-5, -3));
    a.push_back(make_pair(-3, -1));
    a.push_back(make_pair(-1, -3));
    a.push_back(make_pair(-2, -2));
    a.push_back(make_pair(-1, -1));
    a.push_back(make_pair(-2, -1));
    a.push_back(make_pair(-1, 1));

    int n = a.size();

    // sorting the set of points according
    // to the x-coordinate
    sort(a.begin(), a.end());
    vector<pair<int, int> >ans = divide(a);

    cout << "convex hull:\n";
    for (auto e:ans)
    cout << e.first << " "
            << e.second << endl;

    return 0;
}

输出

凸包(convex hull)：
-5 -3 
-1 -5 
1 -4 
0 0 
-1 1

时间复杂度：合并左右凸包需要 O(n) 时间，由于我们将点集分成两个相等的部分，因此上述算法的时间复杂度为 O(n * log n)。 

辅助空间： O(n)

如果您喜欢此文章，请收藏、点赞、评论，谢谢，祝您快乐每一天。