凸包多边形最小外切矩形算法

邱秋 • 2017年10月13日 • 阅读:4275 • algorithm


其实我对算法不是很在行, 但是项目中有用到某种算法 来实现某种功能, 也得硬着头皮来实现. 这是很早之前的一个项目了, 要计算一个凸包多边形最小外切矩形 . 遇到这种情况肯定是束手无策.. 在翻了一些资料之后. 终于完成了.

先说说项目要干嘛:

有这么一个 Desktop app, 其连接到外接摄像头之后, 通过摄像头来捕捉图形, 通过计算来实现某种功能

比如我们做 app 的时候, 可以利用摄像头来识别银行卡, 再通过 OCR 来 识别卡上的数字, 达到快速填写银行卡号的目的.

再比如银行需要客户身份证复印件的时候, 只需要把 你的身份证 放到摄像头下点下鼠标, 就能完成电子档的复印件, 而抛弃扫描仪.

再比如去停车场停车的时候, 那个 led 显示器上能快速显示你的车牌号.  这些都是通过高清摄像头来捕捉头像, 然后进行一系列的计算 , 再通过 OCR 识别 来做到的


如下图, 从图片中取出绿色圆部分.
首先我们拿到图片之后, 经过 图片去灰度---> 二值化---> 查找边缘--->得到最大的 blob, 通过坐标计算, 最后就能 裁剪出圆的部分.

那么对象是较为复杂的图形呢, 比如 三角形, 五角星, 不规则多边形 改如何去处理呢.

任何一张图片他最终的形状是矩形, 那么我们就可以通过 计算不规则多边形的最小外切矩形, 然后通过角度摆正 90° , 就能拿到想要的图形.

凸多边形的最小包围矩形至少有一条边与多边形的一条边共线。

暴力算法

遍历每一条边构造包围矩形比较面积大小。就是构造包围矩形,取到投影点到边以及垂直边上取距离最远两点距离即可

/* min value */
#define FLT_MIN 1.175494351e-38F
/* max value */
#define FLT_MAX 3.402823466e+38F

struct OBB {
    Point u[2]; //x, y轴
    Point c; //中心点
    float e[2]; //半长,半宽
};

float MinAreaRec(Point *pts, int ptsNum, OBB &obb) {

    float minArea = FLT_MAX;

    for(int i = 0, j = ptsNum - 1; i < ptsNum; j = i, i++) {//遍历边
        Point u0 = pts[i] - pts[j]; //构造边
        u0 = u0/Length(u0);
        Point u1 = Point(0-u0.y, u0.x); //与u0垂直
        float min0 = 0.0f, max0 = 0.0f, min1 = 0.0f, max1 = 0.0f;

        for(int k = 0; k < ptsNum; k++) {//遍历点
            Point d = pts[k] - pts[j];
            //投影在u0
            float dot = Dot(d, u0);
            if(dot < min0) min0 = dot;
            if(dot > max0) max0 = dot;
            //投影在u1
            dot = Dot(d, u1);
            if(dot < min1) min1 = dot;
            if(dot > max1) max1 = dot;
        }

        float area = (max0 - min0) * (max1 - min1);
        if( area < minArea ) {
            minArea = area;
            obb.c = pts[j] + ( u0 * (max0 + min0) + u1 * (max1 + min1) )*0.5f;
            obb.u[0] = u0;
            obb.u[1] = u1;
            obb.e[0] = (max0 - min0)*0.5f;
            obb.e[1] = (max1 - min1)*0.5f;
        }
    }

    return minArea;

}

Point * GetOBBPoints(OBB obb) {//获取OBB四个顶点坐标
    Point *pts = new Point [4];
    pts[0] = obb.c + ( obb.u[0] * obb.e[0] + obb.u[1] * obb.e[1] );
    pts[1] = obb.c + ( obb.u[0] * obb.e[0] - obb.u[1] * obb.e[1] );
    pts[2] = obb.c - ( obb.u[0] * obb.e[0] + obb.u[1] * obb.e[1] );
    pts[3] = obb.c + ( obb.u[1] * obb.e[1] - obb.u[0] * obb.e[0] );
    return pts;

}

旋转卡尺(旋转卡壳)算法

使用旋转卡尺算法可将计算凸多边形的最小包围矩形的时间消耗减少很多..

取坐标上两极值点构成平行线,旋转两线,当线与多边形一条边重合时,计算构成矩形面积。

继续旋转,直至旋转角度超过 90 度。取最小面积。

该算法仅对凸体有效(暴力法对凸体凹体均有效),因此需要先计算凸体,该算法的时间复杂度受限于凸体的计算过程

float Cos(Point v, Point p1) {
    float dot = Dot(v, p1);
    float cos = dot/(Length(v)*Length(p1));
    return cos;
}

void Setu0u1(Point e, Point &u0, Point &u1) {
    //以e方向为x轴方向,设定xy轴
    u0 = e / Length(e);
    u1 = Point( 0 - u0.y, u0.x);
}

int GetMinAngleIndex(int imin1, int imax0, int imax1, int imin0,
                     Point* e, Point u0, Point u1) {
    //返回旋转角度最小(cos值最大)的点的下标
    int imin_angle_index = 0;

    float cos = 0, maxCos = FLT_MIN;
    cos = Cos(e[imin1], u0);
    if(cos > maxCos){maxCos = cos; imin_angle_index = imin1; }

    cos = Cos(e[imax0], u1);
    if(cos > maxCos){maxCos = cos; imin_angle_index = imax0; }
    cos = Cos(e[imax1], Point(0-u0.x, 0-u0.y));
    if(cos > maxCos){maxCos = cos; imin_angle_index = imax1; }
    cos = Cos(e[imin0], Point(0-u1.x, 0-u1.y));
    if(cos > maxCos){maxCos = cos; imin_angle_index = imin0; }
    return imin_angle_index;

}

void SetMinMax(Point*pts, int i, int iu, Point u0, Point u1,
               float &max0, float &min0, float &max1,
               int & new_imax0, int &new_imin0, int &new_imax1)

{
    //找到x轴投影最大最小,y轴投影最大的长度(y轴最小则是重合边上点,长度为0)
    //以及极值点在pts中的下标
    Point d =  pts[i] - pts[iu];
    float dist0 = Dot( d, u0);
    if(dist0 > max0){ max0 = dist0; new_imax0 = i; }
    if(dist0 < min0){ min0 = dist0; new_imin0 = i; }

    float dist1 = Dot( d, u1);
    if(dist1 > max1){ max1 = dist1; new_imax1 = i; }
}

float MinAreaRec2(Point *pts, int ptsNum, OBB &obb) {//旋转卡壳算法
    //必须是凸包
    float minArea = FLT_MAX;
    //初始化边e
    Point *e = new Point[ ptsNum ];
    for(int i = 0; i < ptsNum; i++) {
        e[i] = pts[(i+1)%ptsNum] - pts[i];
    }
    int iu = 0; //以e[0]为重合边
    //初始化u0 u1
    Point u0, u1;
    Setu0u1(e[iu], u0, u1);

    int imax0 = 0, imax1 = 0, imin0 = 0, imin1 = 0;
    float min0 = FLT_MAX, max0 = FLT_MIN,
          max1 = FLT_MIN, min1 = 0; //min1其实可以不需要设定的,始终为0
                                    //只是为了理解方便加上
                                    //要去掉则需要把下方用到的min1都改为0

    //求三个极值坐标
    for( int i = 0; i < ptsNum; i++) {
        SetMinMax(pts, i, iu, u0, u1,
               max0, min0, max1,
               imax0, imin0, imax1) ;
    }

    for(int i = 0; i < ptsNum ; i++) {
        int iminangle = 0;
        iminangle = GetMinAngleIndex((iu+1)%ptsNum, imax0, imax1, imin0, e, u0, u1);
        if(iminangle == 0)break; //旋转回了初始点 没必要继续

        if(iminangle == imax0){imax0 = (iu + 1)%ptsNum; iu = iminangle; }
        else if(iminangle == imax1){imax1 = (iu + 1)%ptsNum; iu = iminangle; }
        else if(iminangle == imin0){imin0 = (iu + 1)%ptsNum; iu = iminangle; }
        else if(iminangle == (iu+1)%ptsNum){iu = (iu+1)%ptsNum; }

        Setu0u1(e[iu], u0, u1); //重设u0u1

        //维护三个极值点
        int new_imax0 = imax0, new_imax1 = imax1, new_imin0 = imin0;
        min0 =FLT_MAX, max0 = FLT_MIN, max1 = FLT_MIN;

        //确定原先imax0在新坐标系中是什么极值
        SetMinMax(pts, imax0, iu, u0, u1,
               max0, min0, max1,
               new_imax0, new_imin0, new_imax1) ;
        //确定原先imax1在新坐标系中是什么极值
        SetMinMax(pts, imax1, iu, u0, u1,
               max0, min0, max1,
               new_imax0, new_imin0, new_imax1) ;
        //确定原先imin0在新坐标系中是什么极值
        SetMinMax(pts, imin0, iu, u0, u1,
               max0, min0, max1,
               new_imax0, new_imin0, new_imax1) ;

        imax0 = new_imax0;
        imax1 = new_imax1;
        imin0 = new_imin0;
        //维护完毕

        //求面积 设置obb
        float area = (max0 - min0)*(max1 - min1);
        if(area < minArea) {
            minArea = area;

            obb.e[0] = (max0 - min0)*0.5f;
            obb.e[1] = (max1 - min1)*0.5f;
            obb.u[0] = u0;
            obb.u[1] = u1;
            obb.c = pts[iu] + ( u0 * (max0 + min0) + u1 * (max1 + min1) )*0.5f;
        }
    }
    return minArea;
}

拿到了四个点, 其他的就好说了, 从主体上裁切 出矩形, 然后平行到 90° 就 OK 了.

[完]

我,秦始皇,打钱!