凸包计算及连通区最小外接矩形计算

声希Censh

已于 2023-03-21 00:09:05 修改

阅读量739

点赞数 1

分类专栏：图像处理算法文章标签：算法 c++ 开发语言

于 2023-03-21 00:05:53 首次发布

本文链接：https://blog.csdn.net/weixin_43294620/article/details/129679380

版权

图像处理算法专栏收录该内容

3 篇文章 4 订阅

订阅专栏

以下是一个简单的算法，用来计算凸包，并根据凸包的坐标来计算最小外接矩形

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#define _USE_MATH_DEFINES
#include <math.h>

using namespace std;

struct Point{
    double x,y;
};

double cross(Point a,Point b,Point c){
    double x1=b.x-a.x,x2=c.x-a.x;
    double y1=b.y-a.y,y2=c.y-a.y;
    return x1*y2-x2*y1;
}

bool cmp(Point a,Point b){
    if(a.x!=b.x) return a.x<b.x;
    return a.y<b.y;
}

vector<Point> convex_hull(vector<Point> &p){
    int n=p.size(),k=0;
    vector<Point> h(2*n);
    sort(p.begin(),p.end(),cmp);

    for(int i=0;i<n;i++){
        while(k>=2 && cross(h[k-2],h[k-1],p[i])<=0) k--;
        h[k++]=p[i];
    }

    for(int i=n-2,t=k;i>=0;i--){
        while(k>t && cross(h[k-2],h[k-1],p[i])<=0) k--;
        h[k++]=p[i];
    }

    h.resize(k-1);
    return h;
}

double min_bounding_rectangle(vector<Point> &p){
    vector<Point> h=convex_hull(p);
    int n=h.size();
    double ans=1e20;
    for(int i=0;i<n;i++){
        int j=(i+1)%n,k=(j+1)%n;
        while(cross(h[i],h[j],h[k])<cross(h[i],h[j],h[(k+1)%n])){
            k=(k+1)%n;
        }
        double dx=h[i].x-h[j].x,dy=h[i].y-h[j].y;
        double area=fabs(cross(h[i],h[j],h[k]))/sqrt(dx*dx+dy*dy);
        ans=min(ans,area);
    }
    return ans;
}


struct rect //四个点的矩形
{
    point first, second, third, fouth;
    double area, perimeter, aspect;
};
//4计算最小外接矩形，利用凸包的边进行旋转,返回的rect中包含最小外接矩形的:4个顶点、面积、周长、纵横比
static rect minBoundRect(const std::vector<point> &convexHullCoor, int methodID=0)
{//用来计算给定坐标点集合的最小外接矩形，返回最小外接矩形四个顶点坐标以及面积、周长
 //需要返回的外接矩形信息
    rect minBoundRect;
    //当凸包点数小于3时，直接返回0
    if(convexHullCoor.size()<3)
    {
        point zero;
        zero.width=0, zero.height=0;
        minBoundRect.first=zero;
        minBoundRect.second=zero;
        minBoundRect.third=zero;
        minBoundRect.fouth=zero;
        minBoundRect.area=0;
        minBoundRect.perimeter=0;  //严格来说两点的周长不为零
        minBoundRect.aspect=0;  //纵横比
        return minBoundRect;
        }
    //计算方位角
    double *pAngle = new double[convexHullCoor.size()];
    for(int i=0; i<convexHullCoor.size()-1; i++)
    {
        pAngle[i] = std::atan2(convexHullCoor[i+1].height-convexHullCoor[i].height,convexHullCoor[i+1].width-convexHullCoor[i].width);
    }
    //最后一个方位角，利用第一个计算
    pAngle[convexHullCoor.size()-1] = std::atan2(convexHullCoor[0].height-convexHullCoor[convexHullCoor.size()-1].height,convexHullCoor[0].width-convexHullCoor[convexHullCoor.size()-1].width);

    //对pi/2取余，旋转角度
    double *pEdgeAngle = new double[convexHullCoor.size()];
    for(int i=0; i<convexHullCoor.size(); i++)
    {
        pEdgeAngle[i] = fmod(pAngle[i],M_PI/2);
        double tmp = fmod(pAngle[i],M_PI/2);
        tmp<0? pEdgeAngle[i] = M_PI/2+tmp:pEdgeAngle[i] =tmp;
    }
    //对pEdgeAngle排序,升序
    std::sort(pEdgeAngle,pEdgeAngle+convexHullCoor.size());
    //去除重复出现的元素,得到去重之后的数量uniqueNum，即只需要旋转以下角度就可以
    int uniqueNum = std::unique(pEdgeAngle,pEdgeAngle+convexHullCoor.size())-pEdgeAngle; //去重函数返回地址为：去重后最后一个不重复元素地址
    //得到坐标矩阵
    std::vector<std::vector<double>> xy(convexHullCoor.size(),std::vector<double>(2,0));
    for(int i=0; i<convexHullCoor.size(); i++)
    {
        xy[i][0] = convexHullCoor[i].width;
        xy[i][1] = convexHullCoor[i].height;
    }
    double maxValue = DBL_MAX;
    double area = DBL_MAX;
    double perimeter = DBL_MAX;
    double aspect=0;
    std::vector<std::vector<double>> rect(4,std::vector<double>(2,0));
    //将数据旋转相应角度
    std::vector<std::vector<double>> rot(2,std::vector<double>(2,0)); //2×2矩阵，double类型，初始化为0
    for(int i=0; i<uniqueNum; i++)
    {
        //计算rot的值,数据旋转-edgeAngle
        rot[0][0] =  cos(-pEdgeAngle[i]);
        rot[0][1] =  sin(-pEdgeAngle[i]);
        rot[1][0] =  -sin(-pEdgeAngle[i]);
        rot[1][1] =  cos(-pEdgeAngle[i]);
        std::vector<std::vector<double>> xyr = matrixMultiply(xy,rot,convexHullCoor.size(),2,2);
        //计算最值，即当前旋转下外接矩形的顶点坐标
        double minX = DBL_MAX, minY = DBL_MAX;
        double maxX = -DBL_MAX, maxY = -DBL_MAX;  //不能用DBL_MIN,  DBL_MIN表示最小的正数，坐标可能为负
        for(int i=0; i<convexHullCoor.size(); i++)
        {
            if(xyr[i][0]<minX) minX = xyr[i][0];
            if(xyr[i][0]>maxX) maxX = xyr[i][0];
            if(xyr[i][1]<minY) minY = xyr[i][1];
            if(xyr[i][1]>maxY) maxY = xyr[i][1];
        }
        //计算矩形的面积、周长以及纵横比
        double area_t = (maxX-minX)*(maxY-minY);
        double perimeter_t = 2*((maxX-minX)+(maxY-minY));
        double aspect_t;
        (maxX-minX)>(maxY-minY)? aspect_t=(maxX-minX)/(maxY-minY):aspect_t=(maxY-minY)/(maxX-minX);
        //默认取最小面积，即methodID=0，也可以取最小周长
        double currentValue=0;
        if(methodID==0) currentValue = area_t;
        else if(methodID==1) currentValue = perimeter_t;
        if(currentValue<maxValue)
        {
            maxValue = currentValue;
            area = area_t;
            perimeter = perimeter_t;
            aspect = aspect_t;
            rect[0][0]=minX, rect[0][1]=minY, rect[1][0]=maxX, rect[1][1]=minY;
            rect[2][0]=maxX, rect[2][1]=maxY, rect[3][0]=minX, rect[3][1]=maxY;
            //先将rot转置
            rot = transpose(rot);
            rect = matrixMultiply(rect,rot,4,2,2);//反变换，得到未旋转的坐标
        }
    }

    minBoundRect.area = area;
    minBoundRect.perimeter = perimeter;
    minBoundRect.first.width = rect[0][0], minBoundRect.first.height=rect[0][1];
    minBoundRect.second.width = rect[1][0], minBoundRect.second.height=rect[1][1];
    minBoundRect.third.width = rect[2][0], minBoundRect.third.height=rect[2][1];
    minBoundRect.fouth.width = rect[3][0], minBoundRect.fouth.height=rect[3][1];
    minBoundRect.aspect = aspect;
    delete [] pAngle;
    delete [] pEdgeAngle;
    return minBoundRect;
}

//矩阵乘法， double类型数组
std::vector<std::vector<double>> matrixMultiply(std::vector<std::vector<double>>& a, std::vector<std::vector<double>>& b, int m1, int n1, int n2)
{//m1、n1为矩阵a的行数和列数，n2为矩阵b的列数
    std::vector<std::vector<double>> result(m1, std::vector<double>(n2, 0));  // 存储结果的矩阵，初始化为0
    for (int i = 0; i < m1; i++) {
        for (int j = 0; j < n2; j++) {
            for (int k = 0; k < n1; k++) {
                result[i][j] += a[i][k] * b[k][j];  // 矩阵乘法公式
            }
        }
    }

    return result;
}

//矩阵转置， double类型数组
std::vector<std::vector<double>> transpose(std::vector<std::vector<double>>& matrix)
{
    int m = matrix.size(); // 原矩阵行数
    int n = matrix[0].size(); // 原矩阵列数
    std::vector<std::vector<double>> res(n, std::vector<double>(m)); // 初始化转置矩阵
    // 将原矩阵中的每个元素放入对应位置的转置矩阵中
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            res[j][i] = matrix[i][j];
        }
    }
    return res;
}

int main(){
    int n;
    cin>>n;
    vector<Point> p(n);
    for(int i=0;i<n;i++){
        cin>>p[i].x>>p[i].y;
    }
    double ans=minBoundRect(p);
    cout<<ans<<endl;
    return 0;
}

以上代码中，struct Point定义了一个二维坐标系中的点，cross()函数用来计算两个向量的叉积，cmp()函数用来排序，convex_hull()函数用来计算凸包，min_bounding_rectangle()函数用来计算凸包的最小外接矩形，最后在main()函数中输入点的坐标，输出计算结果。

对算法本身进行分析，如对下图所示的边缘区域进行凸包求取的计算。

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-SrKt49F3-1679327903495)(C:\Users\Censh77\AppData\Roaming\Typora\typora-user-images\image-20230315141246136.png)]$

在此之前先了解交叉积的相关知识：

交叉积可以被定义为：

模长：（在这里θ表示两向量之间的夹角（共起点的前提下）（0°≤θ≤180°），它位于这两个矢量所定义的平面上。）

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-KBvkg6Tb-1679327904311)(null)]

c的长度在数值上等于以a，**b，**夹角为θ组成的平行四边形的面积。

方向：a向量与b向量的向量积的方向与这两个向量所在平面垂直，且遵守右手定则。（一个简单的确定满足“右手定则”的结果向量的方向的方法是这样的：若坐标系是满足右手定则的，当右手的四指从a以不超过180度的转角转向b时，竖起的大拇指指向是c的方向。）

已知点 $a(x_1,y_1)$ , $b(x_2,y_2)$ ，则叉积的二维坐标表示为：
$a×b=x_1y_2-x_2y_1$
上述代码中用到了cross函数来求交叉积（向量积、外积），正是用到了上述公式。

double cross(Point a,Point b,Point c){
    double x1=b.x-a.x,x2=c.x-a.x;
    double y1=b.y-a.y,y2=c.y-a.y;
    return x1*y2-x2*y1;
}

cross函数参数中有三个点a、b、c，组成的向量ab和ac如下图所示：

在这里插入图片描述

对于叉积ab×ac

若叉积大于0，则ab逆时针旋转到ac；

若叉积等于0，则a、b、c三点共线；

若叉积小于0，则ab顺时针旋转到ac；

知道了以上知识点，回到凸包函数中：

vector<Point> convex_hull(vector<Point> &p){
    int n=p.size(),k=0;
    vector<Point> h(2*n);
    sort(p.begin(),p.end(),cmp);

    for(int i=0;i<n;i++){
        while(k>=2 && cross(h[k-2],h[k-1],p[i])<=0) k--;
        h[k++]=p[i];
    }

    for(int i=n-2,t=k;i>=0;i--){
        while(k>t && cross(h[k-2],h[k-1],p[i])<=0) k--;
        h[k++]=p[i];
    }

    h.resize(k-1);
    return h;
}

该函数参数为点集，即需要从这些点中找到凸包点，首先将这些点进行排序(sort函数)，排序规则如下，即从左到右，从上往下。

bool cmp(Point a,Point b){
    if(a.x!=b.x) return a.x<b.x;
    return a.y<b.y;
}

然后依次遍历这些点集，根据叉积的性质得到凸包点。部分流程如下图所示。

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-K6lGf40Q-1679327903497)(C:\Users\Censh77\AppData\Roaming\Typora\typora-user-images\image-20230315141746986.png)]$

点1为起始点，依次遍历，上图右侧是两种不同情况下的交叉积，可以看到右上的图表示向量12和向量13的交叉积，然而在图像中y(height)方向于直角坐标系方向不一致，经过坐标变换后可得到变换后的示意图，此时应用右手定则是符合的，所以在图像中交叉积大于零的时候，表示从第一条直线到第二条直线是顺时针旋转的。

convex_hull函数有两个for循环，通过遍历分别得到上下两部分的凸包点集合。

有了凸包点集之后，可以通过凸包点集计算该凸包的最小外接矩形。最小外接矩形的计算函数为minBoundRect。该函数执行以下操作：

1、先计算方位角（方位角是用来进行旋转的）

2、将方位角对 $\pi/2$ 进行取余并去除重复角度

3、将边缘坐标写成矩阵的形式（便于做矩阵乘法）

4、遍历旋转角度，并计算相应旋转矩阵（用来和边缘坐标矩阵相乘），得到旋转后的边缘坐标矩阵

5、计算旋转后的边缘坐标矩阵的最大最小height、width值，此矩形即为当前角度下的外接矩形

6、将上一步计算得到的外接矩形进行反旋转（即乘以相应旋转矩阵的转置矩阵）

7、反旋转得到的矩形就是一个原连通区的外接矩形，根据所使用的方法可以选择取面积最小的矩形或是周长最小的矩形，遍历结束就可以得到最小外接矩形的四个点及相应面积、周长、纵横比等参数

以上计算过程详见代码，其中用到了矩阵乘法matrixMultiply以及矩阵转置transpose，另外需要注意利用取余函数fmod时候的正负值问题。