超详细的凸包问题的分治算法说明及C++实现

题目要求: 输入数字, 再输入点(x,y); 以顺时针方向输出最外侧的点, 使外侧的点组成的多边形能包含所有的输入点

参考: <算法设计与分析> 第122-123页

准备步骤: 首先, 将输入点按x轴从小到大排序, 如果x相同,则按y排序

这时候, 我们就能得到p[0]和p[n-1]的直线, 将求解凸包问题分解为求解上凸包和下凸包两个问题.

这里以上凸包为例, 我们要找到一个点Pmax, 使三角形P0-Pmax-Pn-1面积最大. 

直观上也很好理解, 假设我们这里只能找三角形, 那么三角形面积越大, 包含的点肯定就越多.

注意, 这里从P0->Pmax->Pn-1是顺时针顺序, 后面还会用到.

那么怎么求三角形面积呢?

我们可能会想到底x高x1/2的方法, 但是, 这里我们不用这个直观上很简单的方法

我自己理解的原因有两个: 

第一,虽然我这里给出的P0-Pn-1是一条平行于x轴的直线, 但如果说, 它是一条斜线, 这个方法就显然很麻烦;

第二, 这个方法求出来的面积是没有方向的, 这不便于我们进一步把凸包分为上凸包和下凸包进行递归调用

看这张图, 我们可以发现, 如果是上凸包, 从P0->Pmax->Pn-1是顺时针顺序;

如果是下凸包, P0->Pmax->Pn-1是逆时针顺序.

换言之, Pn-1->Pmax->P0是顺时针顺序. (注意这种逆向思维, 后面会用作简化代码)

我们就用这种方式区分上凸包和下凸包

(可能又有人会说了, 我直接用y坐标的大小来判断不就好了嘛. 在这里参见原因一, 如果变为斜线以后就不能用y了)

准备工作做到这里, 我们需要一点数学知识了:

如果我们有三个点A(x1,y1),B(x2,y2), C(x3,y3)

我们可以用以下这个行列式求面积 (也就是我们常说的"叉积")

这个行列式对于平面上任意三角形, 求解面积都很方便, 所得结果是三角形ABC面积的两倍 (可以回忆高中知识~)

而且, 当A->B->C为顺时针顺序时, 该值为正, 反之则为负

 (或者按书上的解释, 当C点在AB左侧时, 该值为正; 我个人觉得顺时针/逆时针的方法更好理解~)

我们可以先写好这个函数的代码:

int Djudge(POINT a1, POINT a2, POINT a3)
{
    int calculate = a1.x*a2.y + a3.x*a1.y + a2.x*a3.y - a3.x*a2.y - a2.x*a1.y - a1.x*a3.y;
    return calculate;
}

然后我们再来考虑, 怎么递归: 

我们先来考虑上凸包, 这样下凸包就可以如法炮制了~

假设有个函数DealLeft(int first, int last)

我就来以这张图上的情况, 大致的讲讲是怎么递归的~

第一步, 令first=0, last=n-1, 也就是说, 调用DealLeft(0,n-1),求P0PmaxPn-1最大的点

我们从图上可以看到, 找到的三角形内部的点无需再考虑,但外部的点有可能被遗漏.

我们用个标记数组visit, 令其值为1, 标记一下找到的这个Pmax

第二步, 令first=Pmax, last=Pn-1, 看作一个新的凸包问题进一步进行求解, 把外部的点包含进来.

我们发现这个例子中,新的点Pmax2和Pmax, Pn-1构成一个逆时针的三角形 (Pmax2在PmaxPn-1直线的右边), 

我们这时候就考虑: 要不要再写一个DealRight函数呢?

我看了很多其他的递归程序, 一般都是写了上凸包和下凸包两个函数, 我觉得是没有必要的, 为什么呢?

回忆一下, 我们刚刚讲到:

如果是下凸包, P0->Pmax->Pn-1是逆时针顺序. Pn-1->Pmax->P0是顺时针顺序

那么, 我们在求解下凸包时,只要调用DealLeft(n-1,0),也就是把n-1看作起始点,  

就可以把逆时针的三角形转化为顺时针的三角形, 而不用再另外写一个函数.

这也是我为什么觉得顺时针, 逆时针会比"左右"要更好理解一些.

讲到这里, 基本上这个问题已经可以解决啦~

但是, 我们再讲一个细节问题:

在这个例子中, 我们发现, 圆圈中圈出来的点也在凸包上, 但是也在直线Pmax2-Pn-1上.

写程序的时候很容易忽略这一点.

我最开始写求最大值的函数时就没有考虑到, 而递归到这一步时, 函数会认为没有面积更大的三角形了, 就把这个点忽略掉.

这也导致我写出来的代码一直通不过, 如何保证这个点的输出也耗费了我大量时间...

我们注意到, 如果扫描到圆圈这个点时, 三点共线, 求出来的行列式为0.

这时候, 我们也用visit数组给这一点做出标记.

(关于这一点, 我不是很确定; 因为测试平台上只给出了一组数据, 欢迎大家讨论~)

最后,再讨论一下"顺时针输出"这个要求.

如果充分理解了上面讲到行列式的作用, 就比较好理解.

我们以P0-Pn-1直线为分界线, 先输出上凸包, 然后再输出下凸包即可

有的时候, 开心就是长成这个样子的~~~

==========================

附上完整代码~

#include <iostream>
using namespace std;
#include <algorithm>
#include <stdlib.h>
#define N 10000
int n = 0;
struct POINT
{
    int x, y;
}p[N],ans[N];
int visit[N],mark[N];
int Djudge(POINT a1, POINT a2, POINT a3)
{
    int calculate = a1.x*a2.y + a3.x*a1.y + a2.x*a3.y - a3.x*a2.y - a2.x*a1.y - a1.x*a3.y;
    return calculate;
}
bool cmpxy(const POINT a, const POINT b) //按x轴排序,如果x相同,按y轴排序
{
    if (a.x != b.x)
        return a.x < b.x;
    else
        return a.y < b.y;
}
void DealLeft(int first, int last)
{
    int max = 0, index = -1;
    int i = first;
    if (first < last)
    {
        for (i = first+1; i < last; i++) //注意两端,对于first和last,没必要再进行计算
        {
            int calcu = Djudge(p[first], p[i], p[last]);
            if (calcu == 0) {  visit[i] = 1; } //
            if (calcu > max)
            {
                max = calcu;
                index = i;
            }

        }
    }
    else
    {
        for (i-1; i >last; i--) //如果first>last,重复上述过程,注意这里下界不是0.
        {
            int calcu = Djudge(p[first], p[i], p[last]);
            if (calcu == 0) {visit[i] = 1;} //
            if (calcu >  max)
            {
                max = calcu;
                index = i;
            }
        }
    }
    if (index != -1)
    {
        visit[index] = 1; //对取到的点进行标注  
        DealLeft(first, index);
        DealLeft(index, last);//分治的部分
    }
}

int main()
{
    cin >> n;
    for (int i = 0; i < n; i++)
    {
        cin >> p[i].x >> p[i].y;
        visit[i] = 0;
    }
    visit[0] = 1;
    visit[n - 1] = 1;
    sort(p, p + n, cmpxy);
    DealLeft(0, n - 1); //查找上凸包;
    DealLeft(n - 1, 0); //查找下凸包;
    int t = 0;
    for (int i = 0; i < n; i++)
    {
        if (visit[i] == 1)
        {
            ans[t].x = p[i].x;
            ans[t].y = p[i].y;
            t++;
        }
    }
//顺时针输出
    mark[0] = mark[t - 1] = 1; //数组mark避免重复检查降低效率
    for (int i = 1; i < t - 1; i++)
    {
        mark[i] = 0;
    }
    cout << ans[0].x << " " <<ans[0].y<< endl;
    for (int i =1; i < t-1; i++)
    {
        int d = Djudge(ans[0], ans[t-1], ans[i]);
        if (d >= 0)
        {
            cout << ans[i].x << " " << ans[i].y << endl;
            mark[i] = 1;
        }
    }
    cout << ans[t - 1].x << " " << ans[t - 1].y << endl;
    for (int i = 1; i < t; i++)
    {
        if (mark[i] != 1)
        {
            int d = Djudge(ans[0], ans[t - 1], ans[i]);
            if (d < 0)
            {
                cout << ans[i].x << " " << ans[i].y << endl;
            }
        }
    }

    system("pause");
    return 0;
}