旋转卡壳算法求最小外接矩形代码

旋转卡壳原理：旋转卡壳详解_大学要有梦想的博客-CSDN博客_旋转卡壳

思路：

1.选择卡壳算法用于求凸多边形的最小外接矩形

1.多边形最小的外接矩形一定是以多边形的的一条边为底的一部分

2.通过这条边算出以这条边为底的最上点、最右点和最左点。然后算出以这条边为底的矩形面积

重点：

计算最上点的原理是：

利用到向量叉积求面积的原理，因为求凸包前，所有点是按角度顺序排列好的，向量叉积代表两向量围成的矩形面积，底边不变，高越大，则面积越大，而因为点是按角度排序，所以面积会先变大，后变小，在突然变小的那个地方，说明那个点是最上点。

计算最右点的原理是：

利用向量求点积投影的原理，由于是凸包且点是按角度顺序排列好的，用第一条边对底边投影的大小减第二条边对底边投影的大小，一旦变成正的，说明点是最右点。

计算最左点的原理是：

利用向量求点积投影的原理，由于是凸包且点是按角度顺序排列好的，以最右点那个点作为左起始点计算，用第一条边对底边投影的大小减第二条边对底边投影的大小，一旦变成负的，说明点是最左点。

3.最终通过计算以多边形其他边为底的矩形面积，如果比较全部矩形面积，最小的那个就是外接矩形。

代码：

 /**MABR:Minimum Area Bounding Rectangle
     * @brief :获取由凸包传来的点集，计算最小外接矩形，采用旋转卡壳算法
     * @param[in] ：凸包传来的点集,逆时针顺序的点集
     * @param[out] ：
     * @return ：返回最小外接矩形的四个点集。顺时针方向
     */
    CUBECOM_API vector<iCoord2d> GetMABRByRotatingCaliper(vector<iCoord2d>& points, double eps = 1.0e-6);



 //平面向量叉积
    double Vec2dCross(const iVec2d& fvec, const iVec2d& svec)
    {
      return fvec[0] * svec[1] - fvec[1] * svec[0];
    }
    //平面向量点积
    double Vec2dDot(const iVec2d& fvec, const iVec2d& svec)
    {
      return fvec[0] * svec[0] + fvec[1] * svec[1];
    }
    //点顺着给定的向量偏移scale距离
    void GetPointByVecAndScale(iCoord2d& pt, iVec2d vec, const double scale)
    {   
      //归一化
      vec.NormalizeSelf();
      vec *= scale;
      pt += vec;
    }

CUBECOM_API vector<iCoord2d> CubeHelper::GetMABRByRotatingCaliper(vector<iCoord2d> &points, double eps/* = 1.0e-6*/)
 {
    //归一化
    auto cmpFunc = [=](double x) 
    {
        if (fabs(x) < eps)
            return 0;
        return x > 0 ? 1 : -1;
    };
   
    //初始化用1e18，表示一个大值
    double min_s = 1e18;
   
    size_t ptSize = points.size();
    //首尾连接
    points.push_back(points[0]);
    //左上、左下、右上、右下
    iCoord2d luPt, ldPt, ruPt, rdPt;
    int upPt = 1, rightPt = 1, leftPt = 1;
    for (size_t i = 0; i < ptSize; ++i)
    {
        // 最上面的点,利用叉积求面积，凸包中的点是顺序的，所以以其中一条边为底，
        // 面积会出现由小变大再变小，会出现最高点。
        while (cmpFunc(fabs(CubeHelper::Vec2dCross(points[upPt] - points[i], points[i + 1] - points[i])) -
                       fabs(CubeHelper::Vec2dCross(points[upPt + 1] - points[i], points[i + 1] - points[i]))) <= 0)
        {
            upPt = (upPt + 1) % ptSize;
        }

        // 最右边的点，利用点积的投影，向右边投影越大，凸包中的点是顺序的，
        // 每次相减，由负变正瞬间，说明是相对于这条边最右边的点
        while (cmpFunc(CubeHelper::Vec2dDot(points[rightPt] - points[i], points[i + 1] - points[i]) -
                       CubeHelper::Vec2dDot(points[rightPt + 1] - points[i], points[i + 1] - points[i])) <= 0)
        {
            rightPt = (rightPt + 1) % ptSize;
        }

        if (!i)
            leftPt = rightPt;

        // 最左边的点利用点积的投影，向左边投影越大，凸包中的点是顺序的，
        // 每次相减，由正变负瞬间，说明是相对于这条边最左边的点
        while (cmpFunc(CubeHelper::Vec2dDot(points[leftPt] - points[i], points[i + 1] - points[i]) -
                       CubeHelper::Vec2dDot(points[leftPt + 1] - points[i], points[i + 1] - points[i])) >= 0)
        {
            leftPt = (leftPt + 1) % ptSize;
        }
        //获取底边的长度
        double dis  = points[i].DistPtPt(points[i + 1]);
        if (dis==0.0) continue;
        
        //获取矩形长度，通过点积求投影/底边
        double R      = CubeHelper::Vec2dDot(points[rightPt] - points[i], points[i + 1] - points[i]) / dis;
        double L      = CubeHelper::Vec2dDot(points[leftPt] - points[i], points[i + 1] - points[i]) / dis;
        double length = R - L;
        //获取矩形高度，通过叉积求三角形面积/底边
        double height   = fabs(CubeHelper::Vec2dCross(points[upPt] - points[i], points[i + 1] - points[i])) / dis;
        double rectArea  = length * height;
        if (rectArea < min_s)
        {
            min_s  = rectArea;
            //右下点      
            rdPt = points[i] + (points[i + 1] - points[i]) * (R / dis);
            //右上点
            ruPt = rdPt + (points[rightPt] - rdPt) * (height / points[rightPt].DistPtPt(rdPt));
            //左上点
            luPt = ruPt + (points[i] - rdPt) * (length / R);
            //左下点
            ldPt = luPt + (rdPt - ruPt);
        }
    }
    //存储矩形四个点
    vector<iCoord2d> rectPts {ldPt, luPt, ruPt, rdPt};
    return rectPts;
 }