凸包问题(Graham扫描法)

 

 凸包（百度百科）：

     点集Q的凸包(convex hull)是指一个最小凸多边形，满足Q中的点或者在多边形边上或者在其内。右图中由红色线段表示的多边形就是点集Q={p0,p1,...p12}的凸包。

 

顶点个数n

（1）排序：

    在点集Q中找最左下方的点p0，就是x坐标和y坐标都最小的点，其余的点计算它们的极坐标幅角，以幅角的非降序顺序来排序，如果有幅角相同的，最接近p0的优先。（这是《ACM程序设计培训教程》上的排序方法，由于看不懂什么是幅角，于是上网找了另一种排序方法，效果是一样的，具体见附）

（2）扫描：

    堆栈初始化为chs（n） = {pn-1, p0},栈顶指针ps初始化为1，按排序好的点，从p1到pn-1扫描，如果pi （i=1...n-1）和chs栈顶的两个元素（chs[sp-1],chs[sp]）的三角形符号是正的，说明chs[sp-1],chs[sp],pi是逆时针方向，形成的边是凸边，所以可以把pi加入栈，sp++，继续扫描后面的pi+1；如果三角形符号是负的，那chs[sp]不可能是凸包的极点，出栈，sp--，而pi仍旧是pi。下一轮仍旧是chs[sp-1],chs[sp], pi.

 

附：矢量叉积

   设有点p0(x0,y0),p1(x1,y1),p2(x2,y2).（p0p1）,（p0p2）是共p0的两条向量，叉积d = （p0p1）x(p0p2) = (x1-x0)*(y2-y0) - (x2-x0)*(y1-y0)

   叉积的一个非常重要性质是可以通过它的符号判断两矢量相互之间的顺逆时针关系：
若 d > 0 , 则(p0p1)在(p0p2)的顺时针方向。
若 d < 0 , 则(p0p1)在(p0p2)的逆时针方向。(图示方向)
若 d = 0 , 则(p0p1)与(p0p2)共线，但可能同向也可能反向。

                                       

代码：

 

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <conio.h>
#define N 20
#define inf 1e-6

typedef struct
{
    int x;
    int y;
}point;
point points[N]; //点集

point chs[N];     //栈

int sp; //栈顶指针


//计算两点之间距离

double dis(point a, point b)
{
    return sqrt((a.x - b.x) * (a.x - b.x) * 1.0 + (a.y - b.y) * (a.y - b.y));
}

//通过矢量叉积求极角关系（p0p1）(p0p2)

//k > 0 ,p0p1在p0p2顺时针方向上

double multi(point p0, point p1, point p2)    
{
    return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
}

int cmp(const void *p, const void *q)
{
    point a = *(point *)p;
    point b = *(point *)q;
    double k = multi(points[0], a, b);
    
    if(k < -inf)
        return 1;
    else if(fabs(k) < inf && (dis(a, points[0]) - dis(b, points[0])) > inf) //两点在同一直线上的话，用最近的

        return 1;
    else return -1;
}

void convex_hull(int n)
{
    int i, k, d;
    int miny = points[0].y, index = 0;

    for(i=1; i<n; i++) //找最左下顶点

    {
        if(points[i].y < miny) //找到y坐标最小的点

        {
            miny = points[i].y;
            index = i;
        }
        else if(points[i].y == miny && points[i].x < points[index].x) //相同的话找到x最小的

        {
            index = i;
        }
    }
    //把最左下顶点放到第一个

    point temp;
    temp = points[index];
    points[index] = points[0];
    points[0] = temp;

    qsort(points+1, n-1, sizeof(points[0]), cmp); //p[1:n-1]按相对p[0]的极角从小到大排序

    chs[0] = points[n-1];
    chs[1] = points[0];
    sp = 1;
    k = 1;
    while(k <= n-1)
    {
        double d = multi(chs[sp], chs[sp-1], points[k]);
        if(d <= 0)
        {
            sp++;
            chs[sp] = points[k];
            k++;
        }
        else sp--;
    }
}

int main()
{
    int i, j, n;
    float sum, area;
    freopen("in.txt", "r", stdin);
    scanf("%d", &n);

    for(i=0; i<n; i++)
        scanf("%d%d", &points[i].x, &points[i].y);
    convex_hull(n);
    
    for(i=0; i<=sp; i++)
        printf("(%d, %d) \n ", chs[i].x, chs[i].y);
    sum = 0;
    for(i=1; i<=sp; i++) //求周长

        sum += dis(chs[i-1], chs[i]);
    sum += dis(chs[0], chs[sp]);
    printf("周长：%f\n", sum);

    area = 0;
    for(i=1; i<sp; i++) //求面积

    {
        double k = multi(chs[0], chs[i], chs[i+1]); //三角形面积是向量叉积的1/2

        area += k;
    }
    area = fabs(area) / 2;
    printf("面积：%f\n", area);

    getch();
    return 0;
}

参考：http://hi.baidu.com/fandywang_jlu/blog/item/aa4abd18e49762b74bedbc21.html