清风明月的笔记

旋转卡壳求点集的最小面积外接矩形

旋转卡壳算法的历史
旋转卡壳算法的实现
完整源码

旋转卡壳算法的历史

先转摘一段网上关于旋转卡壳算法的历史

1978年， M.I. Shamos's Ph.D. 的论文"Computational Geometry"标志着计算机科学的这一领域的诞生。当时他发表成果的是一个寻找凸多边形直径的一个非常简单的算法，即根据多边形的一对点距离的最大值来确定。后来直径演化为由一对对踵点对来确定。 Shamos提出了一个简单的 O(n) 时间的算法来确定一个凸 n 角形的对踵点对。因为他们最多只有 3n/2 对，直径可以在 O(n) 时间内算出。如同Toussaint后来提出的， Shamos的算法就像绕着多边形旋转一对卡壳。因此就有了术语“旋转卡壳”。 1983年， Toussaint发表了一篇论文，其中用同样的技术来解决许多问题。从此，基于此模型的新算法就确立了，解决了许多问题。

他们包括：

计算距离
- 凸多边形直径
- 凸多边形宽
- 凸多边形间最大距离
- 凸多边形间最小距离
外接矩形
- 最小面积外接矩形
- 最小周长外接矩形
三角剖分
- 洋葱三角剖分
- 螺旋三角剖分
四边形剖分
- 凸多边形属性
- 合并凸包
- 找共切线
- 凸多边形交
- 临界切线
- 凸多边形矢量和
最薄截面
- 最薄横截带

旋转卡壳算法的实现

旋转卡壳是几何计算中非常重要的算法。这是仅用于实现计算点集的最小面积外接矩形。

求点集的外接矩形需要使用到凸包。有一个很重要的结论，最小覆盖矩形必有一条边和凸包的一条边重合。

算法的步骤如下：

计算点集的凸包
从一条边开始旋转扫描，使用旋转卡壳逆时针求其余三条边的点。求值顺序是下->右->上->左。
- 右顶点可以使用点积求出。很明显|a||b|cosθ的值最大时为右顶点。即向量a在b上的投影最长的时候。
- 上顶点可以使用叉积求出。叉积可以表示两个向量的平行四边形面积。当这两个向量构成的平行四边形面积最大时，即是上顶点。
- 左顶点跟右顶点一样，也可以使用点积求出。但是左顶点的值应该是投影最小的时候。
计算外接矩形的面积
返回2直接遍历所有边

矩形面积计算方法

假设有如下图所求ABCDEF的凸包。

设当前扫描的边为AB，可以得到下点：A，右点:C，上点：E，左点: F

1、可以通过叉积计算出矩形的高EG

设S为AE和AB的叉积，所以S即是这两边构成的平行四边形面积

S = AB*EG

EG = S/AB

2、通过如下图所示过程可以计算出矩形的长CH

$图一$

3、矩形面积 = EG*CH

代码如下

void 
rotatingcalipers(POINT *arr, int len, POINT *rectangle)
{
    int top, down, right=1, up=0, left=0, downlast, rightlast, uplast, leftlast;
    float area = FLT_MAX, dist, X, Y, k;
    POINT temp;

    getconvex(arr, len, &top);
    arr[++top] = arr[0];

    for (down=0; down<top; down++) {
        // find right
        while (getdot(arr[down], arr[down+1], arr[right]) <= getdot(arr[down], arr[down+1], arr[right+1])) {
            right = (right + 1)%top;
        }

        // find up
        if (down == 0) {
            up = right;
        }
        while (getcross(arr[down], arr[down+1], arr[up]) <= getcross(arr[down], arr[down+1], arr[up+1])) {
            up = (up + 1)%top;
        }

        // find down
        if (down == 0) {
            left = up;
        }
        while (getdot(arr[down], arr[down+1], arr[left]) >= getdot(arr[down], arr[down+1], arr[left+1])) {
            left = (left + 1)%top;
        }

        dist = getdist(arr[down], arr[down + 1]);
        X = getcross(arr[down], arr[down + 1], arr[up])/dist;
        temp.x = arr[right].x + arr[down].x - arr[left].x;
        temp.y = arr[right].y + arr[down].y - arr[left].y;
        Y = getdot(arr[down], arr[down + 1], temp);

        if (area > X*Y) {
            area = X*Y;
            downlast = down;
            rightlast = right;
            uplast = up;
            leftlast = left;
        }
    }

    // 计算外接矩形
    if (arr[downlast + 1].y == arr[downlast].y) {
        rectangle[0].x = arr[leftlast].x;
        rectangle[0].y = arr[downlast].y;

        rectangle[1].x = arr[rightlast].x;
        rectangle[1].y = arr[downlast].y;

        rectangle[2].x = arr[rightlast].x;
        rectangle[2].y = arr[uplast].y;

        rectangle[3].x = arr[leftlast].x;
        rectangle[3].y = arr[uplast].y;

    } else if (arr[downlast + 1].x == arr[downlast].x) {
        rectangle[0].x = arr[downlast].x;
        rectangle[0].y = arr[leftlast].y;

        rectangle[1].x = arr[downlast].x;
        rectangle[1].y = arr[rightlast].y;

        rectangle[2].x = arr[uplast].x;
        rectangle[2].y = arr[rightlast].y;

        rectangle[3].x = arr[uplast].x;
        rectangle[3].y = arr[leftlast].y;

    } else {
        k = (arr[downlast + 1].y - arr[downlast].y)/(arr[downlast + 1].x - arr[downlast].x);

        rectangle[0].x = (k*arr[leftlast].y + arr[leftlast].x - k*arr[downlast].y + k*k*arr[downlast].x)/(k*k + 1.0);
        rectangle[0].y = k*rectangle[0].x + arr[downlast].y - k*arr[downlast].x;

        rectangle[1].x = (k*arr[rightlast].y + arr[rightlast].x - k*arr[downlast].y + k*k*arr[downlast].x)/(k*k + 1.0);
        rectangle[1].y = k*rectangle[1].x + arr[downlast].y - k*arr[downlast].x;

        rectangle[2].x = (k*arr[rightlast].y + arr[rightlast].x - k*arr[uplast].y + k*k*arr[uplast].x)/(k*k + 1.0);
        rectangle[2].y = k*rectangle[2].x + arr[uplast].y - k*arr[uplast].x;

        rectangle[3].x = (k*arr[leftlast].y + arr[leftlast].x - k*arr[uplast].y + k*k*arr[uplast].x)/(k*k + 1.0);
        rectangle[3].y = k*rectangle[3].x + arr[uplast].y - k*arr[uplast].x;
    }
}

完整源码

/* rotatingcalipers.c */
/*
 * =====================================================================================
 *
 *       Filename:  rotatingcalipers.c
 *
 *    Description:  旋转卡壳计算最小面积矩形
 *
 *        Version:  1.0
 *        Created:  2017-03-12 21:52
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  simon
 *
 * =====================================================================================
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.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;
}

void 
rotatingcalipers(POINT *arr, int len, POINT *rectangle)
{
    int top, down, right=1, up=0, left=0, downlast, rightlast, uplast, leftlast;
    float area = FLT_MAX, dist, X, Y, k;
    POINT temp;

    getconvex(arr, len, &top);
    arr[++top] = arr[0];

    for (down=0; down<top; down++) {
        // find right
        while (getdot(arr[down], arr[down+1], arr[right]) <= getdot(arr[down], arr[down+1], arr[right+1])) {
            right = (right + 1)%top;
        }

        // find up
        if (down == 0) {
            up = right;
        }
        while (getcross(arr[down], arr[down+1], arr[up]) <= getcross(arr[down], arr[down+1], arr[up+1])) {
            up = (up + 1)%top;
        }

        // find down
        if (down == 0) {
            left = up;
        }
        while (getdot(arr[down], arr[down+1], arr[left]) >= getdot(arr[down], arr[down+1], arr[left+1])) {
            left = (left + 1)%top;
        }

        dist = getdist(arr[down], arr[down + 1]);
        X = getcross(arr[down], arr[down + 1], arr[up])/dist;
        temp.x = arr[right].x + arr[down].x - arr[left].x;
        temp.y = arr[right].y + arr[down].y - arr[left].y;
        Y = getdot(arr[down], arr[down + 1], temp);

        if (area > X*Y) {
            area = X*Y;
            downlast = down;
            rightlast = right;
            uplast = up;
            leftlast = left;
        }
    }

    // 计算外接矩形
    if (arr[downlast + 1].y == arr[downlast].y) {
        rectangle[0].x = arr[leftlast].x;
        rectangle[0].y = arr[downlast].y;

        rectangle[1].x = arr[rightlast].x;
        rectangle[1].y = arr[downlast].y;

        rectangle[2].x = arr[rightlast].x;
        rectangle[2].y = arr[uplast].y;

        rectangle[3].x = arr[leftlast].x;
        rectangle[3].y = arr[uplast].y;

    } else if (arr[downlast + 1].x == arr[downlast].x) {
        rectangle[0].x = arr[downlast].x;
        rectangle[0].y = arr[leftlast].y;

        rectangle[1].x = arr[downlast].x;
        rectangle[1].y = arr[rightlast].y;

        rectangle[2].x = arr[uplast].x;
        rectangle[2].y = arr[rightlast].y;

        rectangle[3].x = arr[uplast].x;
        rectangle[3].y = arr[leftlast].y;

    } else {
        k = (arr[downlast + 1].y - arr[downlast].y)/(arr[downlast + 1].x - arr[downlast].x);

        rectangle[0].x = (k*arr[leftlast].y + arr[leftlast].x - k*arr[downlast].y + k*k*arr[downlast].x)/(k*k + 1.0);
        rectangle[0].y = k*rectangle[0].x + arr[downlast].y - k*arr[downlast].x;

        rectangle[1].x = (k*arr[rightlast].y + arr[rightlast].x - k*arr[downlast].y + k*k*arr[downlast].x)/(k*k + 1.0);
        rectangle[1].y = k*rectangle[1].x + arr[downlast].y - k*arr[downlast].x;

        rectangle[2].x = (k*arr[rightlast].y + arr[rightlast].x - k*arr[uplast].y + k*k*arr[uplast].x)/(k*k + 1.0);
        rectangle[2].y = k*rectangle[2].x + arr[uplast].y - k*arr[uplast].x;

        rectangle[3].x = (k*arr[leftlast].y + arr[leftlast].x - k*arr[uplast].y + k*k*arr[uplast].x)/(k*k + 1.0);
        rectangle[3].y = k*rectangle[3].x + arr[uplast].y - k*arr[uplast].x;
    }
}

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}};
    POINT b[4];
    length = sizeof(a)/sizeof(POINT);
    rotatingcalipers(a, length, b);
    for (i=0; i<4; i++) {
        printf("[%f, %f] ", b[i].x, b[i].y);
    }
    printf("\n");
}

rotatingcalipers.c

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

旋转卡壳求点集的最小面积外接矩形

旋转卡壳算法的历史

旋转卡壳算法的实现

完整源码

Date 2017-03-12(星期日) 21:56 By Ming In 几何算法 Tags 几何算法，算法,