C++旋转卡壳法求最小面积外接矩形

旋转卡壳基本概念介绍：(86条消息) 旋转卡壳详解_大学要有梦想的博客-CSDN博客

OpenCV里面有现成的计算最小面积外接矩形的方法，但是由于我装了好久也没装上opencv，最后还是决定自己实现。

求多边形最小面积外接矩形的基本思路是：

1. 求多边形凸包，避免对内部的一些点求外接矩形；

2. 以凸包的一条边为底边，求凸包在该底边上的最上点up，最左点left，最右点right；

3. 计算外接矩形的四个顶点及面积；

4. 重复2-3，找出面积最小的外接矩形。

up：利用向量叉积求面积的原理计算最上点，向量叉积等于两向量围成的四边形的面积；

s = a x b = 底边长 * 高；

left / right：利用向量点积求投影的原理计算最左点和最右点，最左/右点为投影向量最大的点。

θ    所以有，

#pragma once
#include<cmath>
#include<vector>
#include<tuple>

using namespace std;


class Point
{
public:
	double x;    
	double y;    
	double z;    // z坐标（默认为0，如果需要三维点则给z赋值）

	Point(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}

	Point(const Point& p)
	{
		x = p.x;
		y = p.y;
		z = p.z;
	}

	Point operator = (const Point& p)
	{
		this->x = p.x;
		this->y = p.y;
		this->z = p.z;
		return *this;
	}
	Point operator - (const Point& p)
	{
		return Point(x - p.x, y - p.y, z - p.z);
	}

	Point operator + (const Point& p)
	{
		return Point(x + p.x, y + p.y, z + p.z);
	}

	Point operator * (const double& a)
	{
		return Point(x * a, y * a, z * a);
	}

	Point operator / (const double& a)
	{
		return Point(x / a, y / a, z / a);
	}

	//向量点乘
	double operator * (const Point& p1)
	{
		return x * p1.x + y * p1.y + z * p1.z;
	}

	bool operator == (const Point& p)  
	{
		bool result = fabs(x - p.x) < EPS && fabs(y - p.y) < EPS && fabs(z - p.z) < EPS;
		return result;
	}

private:


};
typedef Point Vector;

typedef vector<Point> Polygon;


//计算凸包
vector<Point> convexhull(const vector<Point>& points)
{
	vector<Point> result;

	if (points.size() < 3)
		return result;

	vector<Point> tmp_points = points;
	// 首先将所有点按照字典序排序
	sort(tmp_points.begin(), tmp_points.end(), cmp);

	// 上凸包
	vector<Point> upper_hull;
	// 存入第一个和第二个点
	upper_hull.push_back(tmp_points[0]);
	upper_hull.push_back(tmp_points[1]);
	for (int i = 2; i < tmp_points.size(); ++i)
	{
		upper_hull.push_back(tmp_points[i]);
		while (upper_hull.size() > 2 && cross(upper_hull[upper_hull.size() - 2] - upper_hull[upper_hull.size() - 3], upper_hull[upper_hull.size() - 1] - upper_hull[upper_hull.size() - 3]).z >= 0)
		{
			upper_hull.erase(upper_hull.end() - 2);
		}
	}
	// 下凸包
	vector<Point> lower_hull;
	// 存入倒数第一第二个点
	lower_hull.push_back(tmp_points[tmp_points.size() - 1]);
	lower_hull.push_back(tmp_points[tmp_points.size() - 2]);
	for (int i = tmp_points.size() - 3; i >= 0; --i)
	{
		lower_hull.push_back(tmp_points[i]);
		while (lower_hull.size() > 2 && cross(lower_hull[lower_hull.size() - 2] - lower_hull[lower_hull.size() - 3], lower_hull[lower_hull.size() - 1] - lower_hull[lower_hull.size() - 3]).z >= 0)
		{
			lower_hull.erase(lower_hull.end() - 2);
		}
	}
	// 删除重复点
	lower_hull.erase(lower_hull.begin());
	lower_hull.erase(lower_hull.end() - 1);

	// 合并上下凸包
	upper_hull.insert(upper_hull.end(), lower_hull.begin(), lower_hull.end());

	result = upper_hull;

	return result;
}



//计算最小外接矩形
PolygonMinAreaBoundRectangle(Polygon& points)
{
	vector<pair<double, Polygon>> boundRec;

	Polygon boundary = convexhull(points);
	boundary.push_back(boundary[0]);

	for (int i = 0; i < boundary.size() - 1; i++)
	{
		Line baseLine(boundary[i], boundary[i + 1]);
		Vector p1 = boundary[i + 1] - boundary[i];


		//计算最上顶点
		Point up = boundary[i];
		double uparea = DBL_MIN;	
		for (int j = 0; j < boundary.size()-1; j++)
		{
			Vector p2 = boundary[j] - boundary[i];
			double areaj = abs(p1.x * p2.y - p2.x * p1.y);
			if (uparea < areaj)
			{
				up = boundary[j];
				uparea = areaj;
			}
		}

		//计算最右边顶点
		Point right = boundary[i + 1];
		double rightshadow = DBL_MIN;
		for (int j = 0; j < boundary.size() - 1; j++)
		{
			Vector p3 = boundary[j] - boundary[i];
			double dot = (p1.x * p3.x + p1.y * p3.y) / norm(p1);
			if (rightshadow < dot)
			{
				right = boundary[j];
				rightshadow = dot;
			}
		}


		//计算最左边顶点
		Point left = boundary[i];
		double leftshadow = DBL_MIN;
		Vector p1Inv = boundary[i] - boundary[i + 1];
		for (int j = 0; j < boundary.size() - 1; j++)
		{
			Vector p4 = boundary[j] - boundary[i + 1];
			double dot = (p1Inv.x * p4.x + p1Inv.y * p4.y) / norm(p1Inv);
			if (leftshadow < dot)
			{
				left = boundary[j];
				leftshadow = dot;
			}
		}


		double dis = norm(p1);
		Vector Rvec = right - boundary[i];
		double R = (Rvec.x * p1.x + Rvec.y * p1.y) / dis;
		Vector Lvec = left - boundary[i +1];
		double L = (Lvec.x * p1Inv.x + Lvec.y * p1Inv.y) / dis;
		double length = R + L - dis;

		Vector upvec = up - boundary[i];
		double height = (upvec.x * p1.y - upvec.y * p1.x) / dis;

		double rec_area = length * height;
		Polygon rec;
		Point right_down = boundary[i] + p1 / dis * R;
		Point left_down = boundary[i + 1] + p1Inv / dis * L;

		double upshadowlen = (upvec.x * p1.x + upvec.y * p1.y) / dis;
		Vector upshadow = (p1 / dis) * upshadowlen;
		Vector heightVec = upvec - upshadow;
		Point right_up = right_down + (heightVec / norm(heightVec)) * height;
		Point left_up = left_down + (heightVec / norm(heightVec)) * height;

		rec.push_back(left_up);
		rec.push_back(right_up);
		rec.push_back(right_down);
		rec.push_back(left_down);
		rec.push_back(left_up);

		pair<double, Polygon> m_rec(rec_area, rec);
		boundRec.push_back(m_rec);

	}

	sort(boundRec.begin(), boundRec.end(), [](const pair<double, Polygon> a, const pair<double, Polygon> b) {
		return a.first < b.first;
		});

	return boundRec[0].second;
}