Graham扫描法求凸包

凸包用不严谨的话来讲,给定二维平面上的点集,凸包就是将最外层的点连接起来构成的凸多边型,它能包含点集中所有的点。

如下图(图片来看网络):

图一

算法的主要方法是先排序,然后扫描。

点集排序

  1. 选取基点。选择所有点中y坐标最小的点,如果多个点都是最小值,选x坐标最小的点。
  2. 剩下的点按与基点构成的向量进行排序。计算向量与x轴的夹角大小来排序,可以使用两个向量叉积来判断第二个向量在第一个向量的左边还是右边。只求向量z轴上的值来判断方向即可。

代码如下

/* 
 *        Name:  getdist
 * Description:  计算两点间的距离
 * 
 */
static float 
getdist(POINT p1, POINT p2)
{
    return sqrt((p2.x - p1.x)*(p2.x - p1.x) + (p2.y - p1.y)*(p2.y - p1.y));
}

/* 
 *        Name:  getcross
 * Description:  计算叉积z方向的值,用于判断两个向量的转向
 * 
 */
static float 
getcross(POINT p0, POINT p1, POINT p2)
{
    return (p1.x - p0.x)*(p2.y - p0.y) - (p2.x - p0.x)*(p1.y - p0.y);
}

/* 
 *        Name:  anglecmp
 * Description:  比较两个向量的逆时针转角方向,p1大于p2时返回大于0,等于时返回等于0
 * 
 */
static float
anglecmp(POINT p0, POINT p1, POINT p2)
{
    float cross = getcross(p0, p1, p2);
    if (cross == 0) {
        return getdist(p0, p2) - getdist(p0, p1);
    }
    return cross;
}

/* 
 *        Name:  vectorsort
 * Description:  根据基点进行向量排序, 快速排序递归实现
 * 
 */
static void
vectorsort(POINT *arr, int left, int right)
{
    int i, mid, last;

    if (left >= right) {
        return;
    }

    mid = (left + right)/2;
    swap(arr, left, mid);
    last = left;
    for (i=left + 1; i<=right; i++) {
        if (anglecmp(arr[0], arr[i], arr[left])>0) {
            swap(arr, i, ++last);
        }
    }
    swap(arr, left, last);
    vectorsort(arr, left, last - 1);
    vectorsort(arr, last + 1, right);
}

取基点那一部分在下面的代码中。

点集扫描

基点是必然在凸包上的。除基点外,逆时针进行计算,每相邻两个点(A,B)与上一个点(T)构成的两个向量(TA,TB)使用叉积计算旋转方向。后面那个向量(TB)必须在前面向量(TA)的左边。如果TB在TA右边,说明TA不在凸包上。只求向量z轴上的值来判断方向即可。注意跳过相同的点。

如下图(图片来看网络):

图一

代码如下

/* 
 *        Name:  getconvex
 * Description:  计算凸包
 * 
 */
void
getconvex(POINT *arr, int len, int *n)
{
    int i, base, top;

    /* 小于4个点的不计算凸包 */
    if (len < 4) {
        *n = len;
        return;
    }
    /* 计算基点,交换到0位置 */
    base = 0;
    for (i=0; i<len; i++) {
        if (arr[i].y == arr[base].y && arr[i].x < arr[base].x) {
            base = i;
        } else if (arr[i].y < arr[base].y) {
            base = i;
        }
    }
    swap(arr, base, 0);

    /* 排序 */
    vectorsort(arr, 1, len-1);

    /* 计算凸包 */
    top = 1;
    for (i=2; i<len; i++) {
        while (top>0 && getcross(arr[top-1], arr[top], arr[i])<=0) {
            top--;
        }
        arr[++top] = arr[i];
    }
    *n = top;
}

完整源码

/* grahamscan.c */
/*
 * =====================================================================================
 *
 *       Filename:  grahamscan.c
 *
 *    Description:  graham扫描凸包
 *
 *        Version:  1.0
 *        Created:  2017-03-10 23:23
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  simon
 *
 * =====================================================================================
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct point_s {
    float x;
    float y;
} POINT;

static void
swap(POINT *arr, int a, int b)
{
    POINT temp = arr[a];
    arr[a] = arr[b];
    arr[b] = temp;
}

/* 
 *        Name:  getdist
 * Description:  计算两点间的距离
 * 
 */
static float 
getdist(POINT p1, POINT p2)
{
    return sqrt((p2.x - p1.x)*(p2.x - p1.x) + (p2.y - p1.y)*(p2.y - p1.y));
}

/* 
 *        Name:  getcross
 * Description:  计算叉积z方向的值,用于判断两个向量的转向
 * 
 */
static float 
getcross(POINT p0, POINT p1, POINT p2)
{
    return (p1.x - p0.x)*(p2.y - p0.y) - (p2.x - p0.x)*(p1.y - p0.y);
}

/* 
 *        Name:  getpod
 * Description:  计算两个向量的点积
 * 
 */
static float 
getdot(POINT p0, POINT p1, POINT p2)
{
    return (p1.x - p0.x)*(p2.x - p0.x) + (p1.y - p0.y)*(p2.y - p0.y);
}

/* 
 *        Name:  anglecmp
 * Description:  比较两个向量的逆时针转角方向,p1大于p2时返回大于0,等于时返回等于0
 * 
 */
static float
anglecmp(POINT p0, POINT p1, POINT p2)
{
    float cross = getcross(p0, p1, p2);
    if (cross == 0) {
        return getdist(p0, p2) - getdist(p0, p1);
    }
    return cross;
}

/* 
 *        Name:  vectorsort
 * Description:  根据基点进行向量排序, 快速排序递归实现
 * 
 */
static void
vectorsort(POINT *arr, int left, int right)
{
    int i, mid, last;

    if (left >= right) {
        return;
    }

    mid = (left + right)/2;
    swap(arr, left, mid);
    last = left;
    for (i=left + 1; i<=right; i++) {
        if (anglecmp(arr[0], arr[i], arr[left])>0) {
            swap(arr, i, ++last);
        }
    }
    swap(arr, left, last);
    vectorsort(arr, left, last - 1);
    vectorsort(arr, last + 1, right);
}

/* 
 *        Name:  getconvex
 * Description:  计算凸包
 * 
 */
void
getconvex(POINT *arr, int len, int *n)
{
    int i, base, top;

    /* 小于4个点的不计算凸包 */
    if (len < 4) {
        *n = len;
        return;
    }
    /* 计算基点,交换到0位置 */
    base = 0;
    for (i=0; i<len; i++) {
        if (arr[i].y == arr[base].y && arr[i].x < arr[base].x) {
            base = i;
        } else if (arr[i].y < arr[base].y) {
            base = i;
        }
    }
    swap(arr, base, 0);

    /* 排序 */
    vectorsort(arr, 1, len-1);

    /* 计算凸包 */
    top = 1;
    for (i=2; i<len; i++) {
        while (top>0 && getcross(arr[top-1], arr[top], arr[i])<=0) {
            top--;
        }
        arr[++top] = arr[i];
    }
    *n = top + 1;
}

int
main(int argc, char *argv[])
{
    int n, i, length;
    POINT a[] = {{1.0, 2.0},{2.0, 0.5},{2.5, 1.0},{2.0, 0.0},{4.0, 2.0}};
    length = sizeof(a)/sizeof(POINT);
    getconvex(a, length, &n);
    for (i=0; i<=n; i++) {
        printf("[%f, %f] ", a[i].x, a[i].y);
    }
    printf("\n");
}

grahamscan.c

转载文章请注明出处: http://mingnote.com


Date 2017-03-11(星期六) 18:42 By Ming In 几何算法 Tags 几何算法,算法,