凸包问题的分治算法_凸包问题——Graham算法

本篇文章将介绍二维的Graham扫描算法和它的源码实现。

Graham扫描算法会先将点按照极角的大小顺序排列，接着按顺序遍历每个点，通过夹角的大小判断哪些点在凸包上，算法的伪码如下所示:

求给定二维点集

的凸包

1. 求出最左下角的点，即
   分量的值最小点，若
   分量值最小的点有多个，取
   分量最小的点，设为
   ;
2. 剩下的点集，根据极角大小，逆时针排序，得到
   ;
3. 令栈
   为空，用
   表示栈中第
   个元素，用
   表示栈中的元素个数；若栈中有
   个元素，则
   即为栈顶元素；每一次PUSH}操作，栈元素个数
   加1；每次POP操作，栈元素个数
   减1;
4. PUSH(
   ,
   );
5. PUSH(
   ,
   );
6. for(
   ;
   ;
   )
7.     while (
   && 由
   ,
   和
   的夹角
   )
8.         POP(
   );
9.     end while;
10. end for;
11. if
    , then
12.     点集无法生成一个凸包;
13. else
14.     栈
    中的
    个元素，就是凸包的顶点;
15. end if;

图1. Graham扫描算法过程

以如图1(a)所示为例，给定一个点集，

。第1步，选择最左下角的点

，若

分量值最小的点有多个，取

分量最小的点，特别注意，在实际实现中，可能存在多个相同的

点，必需把它们排除，保证点

是唯一的，否则在第2步中的排序会造成错误。第2步，如图(b)所示，对

按照极角从小到大进行排序，若极角相同，则根据与点

的距离从近到远排序。第3 ~ 5步，如图(c)所示，把

和

两个点压入栈中。逆时针方向处理每个顶点，第6 ~ 10步，如图(d)~(i)，若

且由

,

和

的夹角大于等式180度，则把栈顶元素出栈，否则把新点

压入栈。最终得到的凸包如图(k)所示，点

在线段

上，是否把它作为凸包的顶点，可以根据具体的问题再作决定，上面提供的算法伪码，不允许把点

作为凸包的顶点。

对于对极点的排序，有多种算法实现，比如冒泡排序，快速排序，堆排序等，各种排序的性能对比，如表1的总结。若采用冒泡排序，则Graham扫描算法的复杂度将达到

，采用其它几种排序算法，Graham算法的复杂度都能达到

。采用快速排序，平均复杂度是

，在C++中，提供了qsort()快速排序的算法，但在最坏情况下，Graham算法的复杂度仍然是

。因此，可以选择归并排序或者堆排序。

表1. 有关排序算法的分析结果

基础库源码链接，参见这里，下面是前面所描述的算法的实现。

#include "SrGeometricTools.h"
#include "SrDataType.h"
#include <stdio.h>
#include <list>

namespace {
    /*
brief  找到最左下角的顶点，最左下角的顶点有多个，则只保留一个顶点。
*/
    int GrahamPivot(SrPoint2D* points, int numPoint, int& newNumPoint)
    {
        int i, minIndex = 0, numIndex = 1;

        for (i = 1; i < numPoint; i++)
        {
            if (EQUAL(points[i].x, points[minIndex].x) && EQUAL(points[i].y, points[minIndex].y))
            {
                continue;
            }
            if (LESS(points[i].x, points[minIndex].x))
            {
                minIndex = numIndex;
            }
            else if (EQUAL(points[i].x, points[minIndex].x) && LESS(points[i].y, points[minIndex].y))
            {
                minIndex = numIndex;
            }

            points[numIndex++] = points[i];
        }
        newNumPoint = numIndex;
        return minIndex;
    }

    /*
    brief  对比极角 p0-pivot-p1：
            若极角小于90，return  1；
            若极角大于90，return -1；
            若极角为0且p0到pivot的距离小于p1到pivot的距离，return 1；
            若极角为0且p0到pivot的距离大于p1到pivot的距离，return -1；
            否则，p0与p1重叠，return 0。
    */
    int GrahamCompare(const SrPoint2D& p0, const SrPoint2D& p1, const SrPoint2D& pivot)
    {
        SrReal angle = (p0 - pivot).cross(p1 - pivot);
        if (EQUAL(angle, 0))
        {
            SrReal dif = pivot.distanceSquared(p0) - pivot.distanceSquared(p1);
            if (EQUAL(dif, 0))
                return 0;
            return  GREATER(dif, 0) ? -1 : 1;
        }
        return GREATER(angle, 0) ? 1 : -1;
    }
    /*
    brief  采用堆排序，将顶点按照极角从小到大排序。
    */
    void GrahamHeapSort(SrPoint2D * points, int numPoint, const SrPoint2D& pivot)
    {
        SrPoint2D* heap = new SrPoint2D[numPoint + 1];
        int i, j, n = numPoint, child;
        for (i = 0; i < numPoint; i++)
            heap[i + 1] = points[i];
        //build a minimum heap.
        SrPoint2D y;
        for (i = numPoint / 2; i >= 1; i--)
        {
            y = heap[i];
            child = 2 * i;
            while (child <= numPoint)
            {
                if (child < numPoint&&GrahamCompare(heap[child], heap[child + 1], pivot) < 0)/*heap[child]>heap[child+1]*/
                    child++;
                if (GrahamCompare(y, heap[child], pivot) >= 0)
                    break;
                heap[child / 2] = heap[child];
                child *= 2;
            }
            heap[child / 2] = y;
        }
        //Pop the head of heap and update the minimum heap.
        int indx = 0;
        for (j = 1; j <= numPoint; j++)
        {
            y = heap[n--];
            points[indx++] = heap[1];
            i = 1, child = 2;
            while (child <= n)
            {
                if (child < n&&GrahamCompare(heap[child], heap[child + 1], pivot) < 0)
                    child++;
                if (GrahamCompare(y, heap[child], pivot) >= 0)
                    break;
                heap[i] = heap[child];
                i = child;
                child *= 2;
            }
            heap[i] = y;
        }
        delete[]heap;
    }
    /*
    brief  采用Graham扫描线法计算凸包。
    */
    void    GrahamScanHull(const SrPoint2D* pPoints, int numPoints, std::list<SrPoint2D>& result)
    {
        int i;
        SrPoint2D* buffer = new SrPoint2D[numPoints];
        for (i = 0; i < numPoints; i++)
            buffer[i] = pPoints[i];
        int newNumPoints;
        int minIndex = GrahamPivot(buffer, numPoints, newNumPoints);

        SrPoint2D pivot;
        pivot = buffer[minIndex];
        buffer[minIndex] = buffer[0];
        buffer[0] = pivot;

        GrahamHeapSort(buffer + 1, newNumPoints - 1, pivot);

        result.push_back(buffer[0]);
        result.push_back(buffer[1]);
        std::list<SrPoint2D>::iterator last, current;
        for (i = 2; i < newNumPoints; i++)
        {
            while (result.size() >= 2)
            {
                last = result.end();
                last--;
                current = last;
                last--;
                if (LEQUAL((buffer[i] - *current).cross(*last - *current), 0))
                    result.pop_back();
                else
                    break;
            }
            result.push_back(buffer[i]);
        }

        delete[]buffer;
    }
    /*
    brief  用于判断给定的多边形是否是凸包
    */
    bool InitConvex(const SrPoint2D* point, int numPoint, int &mNumVertex, SrPoint2D*& mVertex)
    {
        if (numPoint <= 2)
            return false;
        int i;
        std::list<SrPoint2D> result;
        GrahamScanHull(point, numPoint, result);
        mNumVertex = result.size();
        if (mNumVertex <= 2)
        {
            mNumVertex = 0;
            return false;
        }
        mVertex = new SrPoint2D[mNumVertex];
        std::list<SrPoint2D>::iterator iter;
        for (iter = result.begin(), i = 0; iter != result.end(); iter++, i++)
            mVertex[i] = *iter;
        return true;
    }
    /*
    brief  Compare the 2d points lexicographically.p<q <==> (px<qx) or ((px=qx)and(py<qy)).
    return If pi<pj, return 1;
            If pi>pj, return -1;
            If pi==pj,return 0;
    */
    int CompareVertex(const SrPoint2D& p0, const SrPoint2D& p1)
    {
        if (LESS(p0.x, p1.x)) return 1;
        if (GREATER(p0.x, p1.x)) return -1;
        if (LESS(p0.y, p1.y)) return 1;
        if (GREATER(p0.y, p1.y)) return -1;
        return 0;
    }

    /*
    brief 判断多边形是否是凸多边形
    */
    bool IsConvex(SrPoint2D* vertex, int numVertex)
    {
        int last, current, next;
        int sign = CompareVertex(vertex[numVertex - 1], vertex[0]), tmpSign;
        int signChange = 0;
        for (last = 0, current = 1; current < numVertex; last = current, current++)
        {
            tmpSign = CompareVertex(vertex[last], vertex[current]);
            if (tmpSign == 0)
                return false;
            if (tmpSign != signChange)
            {
                signChange = tmpSign;
                signChange++;
                if (signChange > 2)
                    return false;
            }
        }
        bool isGreater = false;
        SrReal angle;
        for (last = numVertex - 2, current = numVertex - 1, next = 0; next < numVertex; last = current, current = next, next++)
        {
            angle = (vertex[next] - vertex[current]).cross(vertex[last] - vertex[current]);
            if (LESS(angle, 0))
                return false;
            else if (!isGreater && GREATER(angle, 0))
                isGreater = true;
        }
        if (isGreater)
            return true;
        //Degenerate case.
        return false;
    }
}


void TestGrahamConvex()
{
    int numPoint = 10000, numCase = 100;
    int cs = 0;
    SrPoint2D* points = new SrPoint2D[numPoint];
    while (numCase--)
    {
        int i;
        for (i = 0; i < numPoint; i++)
        {
            points[i].x = rand() % 1000;
            points[i].y = rand() % 1000;
        }
        SrPoint2D* result = NULL;
        int numResult;
        InitConvex(points, numPoint, numResult, result);
        if (IsConvex(result, numResult))
            printf("Case %d Succeeds!n", ++cs);
        else
            printf("Case %d Fails!n", ++cs);
        delete[]result;
    }
    delete[]points;
}

int main()
{
    TestGrahamConvex();
    return 0;
}