计算几何 - 寻找凸包算法

寻找凸包的算法

点集Q的凸包，是一个最小的凸多边形P，满足Q中的每个点都在P的边界上，或者在P的内部。

Graham扫描法: 复杂度O(nlogn)

1. 选取y最小的点，多个y最小的话，选取其中x最小的点，作为p0

2. 剩余的点，按照p0和pi的极角的逆时针排序，编号为p1,p2,...,pm

3. 如果m小于2，表示点数小于3，形不成多边形

4. 设定以空栈S，将p0、p1、p2压入栈中。

5. for i=3 to m

   1. 得到栈顶的2个点pi-1和pi-2，如果t1t0向t0pi转的时候，不是左转，就把顶点t0出栈；

   2. 如果出栈了t0，就继续a，直到栈的顶点不再出栈位置

   3. 将pi入栈

6. return S

图中，从a到f是一步一步选择的过程。

Jarvis步进法：复杂度O(nh)，h是凸包顶点数

先找到最下边结点里最左边的点p0，然后寻找使得p0p1极角最小的点，则p1也是凸包顶点；继续寻找使得p1p2极角最小的点，直到达到最高点pk，上图是p3，此时已经构造好了CH(Q)的右链。为了构造左链，寻找pk+1使得pkpk+1极角最小，但此时x轴啊原x轴的负方向。

除此之外，还有增量法、分治法、剪枝-搜索方法，其中剪枝-搜索方法复杂度为O(nlgh)。

向量知识

确定连续线段是向左转还是向右转

点积等于零是指两向量垂直。

叉积等于零是指两向量平行；叉积大于0，右手定则大拇指朝向自己，向量按逆时针方向首尾相接；叉积小于0，右手定则大拇指朝外，向量按顺时针方向首尾相接。

对于向量p1和p2，叉积是由点(0,0)、p1、p2和p1+p2构成的平行四边形的有向面积。另一种与之等价但更有效的的叉积定义方式是将其看做矩阵行列式：
p1×p2 = x1y2 - x2y1 = - p2×p1
若p1×p2为正，则相对于原点(0,0)来说，p1位于p2顺时针方向；若p1×p2为负，p1位于p2逆时针方向；若为0则方向相同，或相反。

若是相对于点p0(x0,y0)而非原点，则p0p1和p0p2的叉积为(p1-p0)×(p2-p0) = (x1-x0)(y2-y0)-(x2-x0)(y1-y0)。

确定连续线段是向左转还是向右转

对于线段p0p1和p1p2，采用叉积可以避免计算角度，只需简单的计算一下p0p2是位于p0p1的顺时针还是逆时针方向。计算叉积
(p2-p0)×(p1-p0) = (x2-x0)(y1-y0) - (x1-x0)(y2-y0)
若结果为负，p0p2在p0p1的逆时针方向，在p1处左转；结果为正则右转；为0表示三点共线。

Graham扫描法: 复杂度O(nlogn)

http://blog.csdn.net/hqd_acm/article/details/6218420
该博文写的算是通俗易懂，故推荐。
博文代码中有一处小错误，已改正如下。

模板代码

#include <algorithm>
#include <iostream>
#include <vector>
#include <math.h>
using namespace std;
//二维点(或向量)结构体定义
#ifndef _WINDEF_
struct POINT { int x; int y; };
#endif
typedef vector<POINT> PTARRAY;
//判断两个点(或向量)是否相等
bool operator==(const POINT &pt1, const POINT &pt2) {
    return (pt1.x == pt2.x && pt1.y == pt2.y);
}
// 比较两个向量pt1和pt2分别与x轴向量(1, 0)的夹角
bool CompareVector(const POINT &pt1, const POINT &pt2) {
    //求向量的模
    float m1 = sqrt((float)(pt1.x * pt1.x + pt1.y * pt1.y));
    float m2 = sqrt((float)(pt2.x * pt2.x + pt2.y * pt2.y));
    //两个向量分别与(1, 0)求内积
    float v1 = pt1.x / m1, v2 = pt2.x / m2; 
    //通过两个向量对应单位向量横坐标的大小简介比较它们与x轴夹角的大小
    //自己的判断方法： 以下x轴、y轴坐标，指的都是单位向量，从x轴开始逆时针旋转， 
    //y1 < 0 <= y2 ,则向量1的夹角大于向量2的夹角
    // y1,y2 >= 0 , x值小的的夹角大
    // y1,y2 <= 0, x值大的夹角大
    return (v1 > v2 || (v1 == v2 && m1 < m2)); //这一步不理解
}
//计算凸包
void CalcConvexHull(PTARRAY &vecSrc) {
    //点集中至少应有3个点，才能构成多边形
    if (vecSrc.size() < 3) {
        return;
    }
    //查找基点
    POINT ptBase = vecSrc.front(); //将第1个点预设为最小点
    /*在所有点中选取y坐标最小的一点H，当作基点。
      如果存在多个点的y坐标都为最小值，则选取x坐标最小的一点。
      坐标相同的点应排除。
    */
    for (PTARRAY::iterator i = vecSrc.begin() + 1; i != vecSrc.end(); ++i) {
        //如果当前点的y值小于最小点，或y值相等，x值较小
        if (i->y < ptBase.y || (i->y == ptBase.y && i->x < ptBase.x)) {
            //将当前点作为最小点
            ptBase = *i;
        }
    }
    //计算出各点与基点构成的向量
    for (PTARRAY::iterator i = vecSrc.begin(); i != vecSrc.end();) {
        //排除与基点相同的点，避免后面的排序计算中出现除0错误
        if (*i == ptBase) {
            i = vecSrc.erase(i);
        }
        else {
            //方向由基点到目标点
            i->x -= ptBase.x, i->y -= ptBase.y;
            ++i;
        }
    }
    //按各向量与横坐标之间的夹角排序
    sort(vecSrc.begin(), vecSrc.end(), &CompareVector);
    //删除相同的向量
    vecSrc.erase(unique(vecSrc.begin(), vecSrc.end()), vecSrc.end());//  unique是啥用法？
    //计算得到首尾依次相联的向量
    for (PTARRAY::reverse_iterator ri = vecSrc.rbegin();
        ri != vecSrc.rend() - 1; ++ri) {
        PTARRAY::reverse_iterator riNext = ri + 1;
        //向量三角形计算公式
        ri->x -= riNext->x, ri->y -= riNext->y;
    }
    //依次删除不在凸包上的向量
    for (PTARRAY::iterator i = vecSrc.begin() + 1; i != vecSrc.end(); ++i) {
        //回溯删除旋转方向相反的向量，使用外积判断旋转方向
        for (PTARRAY::iterator iLast = i - 1; iLast != vecSrc.begin();) {
            int v1 = i->x * iLast->y, v2 = i->y * iLast->x;
            //如果叉积小于0，则没有逆向旋转
            //如果叉积等于0，还需判断方向是否相逆
            if (v1 < v2 || (v1 == v2 && i->x * iLast->x > 0 &&
                i->y * iLast->y > 0)) {
                break;
            }
            //删除前一个向量后，需更新当前向量，与前面的向量首尾相连
            //向量三角形计算公式
            i->x += iLast->x, i->y += iLast->y;
            iLast = (i = vecSrc.erase(iLast)) - 1;
        }
    }
    //将所有首尾相连的向量依次累加，换算成坐标
    vecSrc.front().x += ptBase.x, vecSrc.front().y += ptBase.y;
    for (PTARRAY::iterator i = vecSrc.begin() + 1; i != vecSrc.end(); ++i) {
        i->x += (i - 1)->x, i->y += (i - 1)->y;
    }
    //添加基点，全部的凸包计算完成
    vecSrc.push_back(ptBase);
}

int main(void) {
    int nPtCnt = 10; //生成的随机点数
    PTARRAY vecSrc; //初始点集 和  最后的凸包都存在vecSrc中
    for (int i = 0; i < nPtCnt; ++i) {
        POINT ptIn = { rand() % 16, rand() % 16 };
        vecSrc.push_back(ptIn);
        cout << ptIn.x << ", " << ptIn.y << endl;
    }
    CalcConvexHull(vecSrc);
    cout << "\nConvex Hull:\n";
    for (PTARRAY::iterator i = vecSrc.begin(); i != vecSrc.end(); ++i) {
        cout << i->x << ", " << i->y << endl;
    }
    return 0;
}

代码语法(STL)解读

vector.push_back();

拓展练习

简单变形题1 - 判断凸多边形、凹多边形

http://www.acmerblog.com/hdu-2108-shape-of-hdu-3252.html

拓展阅读

计算几何基础——矢量和叉积 && 叉积、线段相交判断、凸包