RotateRect.h
#include <vector>/**
* @brief 计算点集最小旋转外接矩形
*/
class RotateRect {
public:enum { CALIPERS_MAXHEIGHT = 0, CALIPERS_MINAREARECT = 1, CALIPERS_MAXDIST = 2 };struct Point {float x, y;};using Points = std::vector<Point>;struct Sizef {float width, height;};struct Rect {float left, top, right, bottom;};explicit RotateRect(const Points& pts);const Point& center() const;const Sizef& size() const;float angle() const;
public:void update();Points toPoints() const;Rect getOutLine() const;private:double crossProduct(const Point& A, const Point& B, const Point& C) const;double calculateDistance(const Point& A, const Point& B) const;Points convexHull() const;/* we will use usual cartesian coordinates */void rotatingCalipers(const Point* pts, int n, int mode, float* out) const;private:const Points& m_inputs;Point m_center;//! returns width and height of the rectangleSizef m_size;//! returns the rotation angle. When the angle is 0, 90, 180, 270 etc., the rectangle becomes an up-right rectangle.float m_angle;
};
RotateRect.cpp
#include "RotateRect.h"
#include <cmath>
#include <assert.h>#ifndef CV_PI
#define CV_PI 3.1415926535897932384626433832795
#endifRotateRect::RotateRect(const RotateRect::Points& pts) : m_inputs(pts), m_angle(0), m_size{ 0, 0 }, m_center{ 0, 0 } {}/** @brief 计算BA与CA之间面积* @param [in] A 点* @param [in] B 点* @param [in] C 点* @returnBA x CA > 0, C在B的逆时针方向,B在C右边BA x CA < 0, C在B的顺时针方向,B在C左边BA x CA = 0, C和B共线,可能正方向,也可能反方向*/
double RotateRect::crossProduct(const RotateRect::Point& A, const RotateRect::Point& B, const RotateRect::Point& C) const {return (B.x - A.x) * (C.y - A.y) - (C.x - A.x) * (B.y - A.y);
}/** @brief 找到x方向最小点为起始点,查找所有点最右边点,到起始点结束* @return 所有凸包点,逆时针方向*/
RotateRect::Points RotateRect::convexHull() const {const auto& points = m_inputs;int n = points.size();if (n <= 3) {return points;}std::vector<Point> hull;int l = 0;for (int i = 1; i < n; i++) {if (points[i].x < points[l].x) {l = i;}}int p = l, q;do {hull.push_back(points[p]);q = (p + 1) % n;for (int i = 0; i < n; i++) {if (crossProduct(points[p], points[i], points[q]) > 0) {q = i;}}p = q;} while (p != l);return hull;
}double RotateRect::calculateDistance(const RotateRect::Point& A, const RotateRect::Point& B) const {return std::sqrt(std::pow(B.x - A.x, 2) + std::pow(B.y - A.y, 2));
}void rotate90CCW(const RotateRect::Point& in, RotateRect::Point& out)
{out.x = -in.y;out.y = in.x;
}void rotate90CW(const RotateRect::Point& in, RotateRect::Point& out)
{out.x = in.y;out.y = -in.x;
}void rotate180(const RotateRect::Point& in, RotateRect::Point& out)
{out.x = -in.x;out.y = -in.y;
}/* return true if first vector is to the right (clockwise) of the second */
bool firstVecIsRight(const RotateRect::Point& vec1, const RotateRect::Point& vec2)
{RotateRect::Point tmp;rotate90CW(vec1, tmp);return tmp.x * vec2.x + tmp.y * vec2.y < 0;
}/** @brief 计算凸包点的最小旋转矩形*//* we will use usual cartesian coordinates */
void RotateRect::rotatingCalipers(const RotateRect::Point* points, int n, int mode, float* out) const
{using Point = RotateRect::Point;float minarea = FLT_MAX;float max_dist = 0;char buffer[32] = {};int i, k;std::vector<float> abuf(n * 3);float* inv_vect_length = abuf.data();Point* vect = (Point*)(inv_vect_length + n);int left = 0, bottom = 0, right = 0, top = 0;int seq[4] = { -1, -1, -1, -1 };Point rot_vect[4];/* rotating calipers sides will always have coordinates(a,b) (-b,a) (-a,-b) (b, -a)*//* this is a first base vector (a,b) initialized by (1,0) */float orientation = 0;float base_a;float base_b = 0;float left_x, right_x, top_y, bottom_y;Point pt0 = points[0];left_x = right_x = pt0.x;top_y = bottom_y = pt0.y;for (i = 0; i < n; i++){double dx, dy;if (pt0.x < left_x)left_x = pt0.x, left = i;if (pt0.x > right_x)right_x = pt0.x, right = i;if (pt0.y > top_y)top_y = pt0.y, top = i;if (pt0.y < bottom_y)bottom_y = pt0.y, bottom = i;Point pt = points[(i + 1) & (i + 1 < n ? -1 : 0)];dx = pt.x - pt0.x;dy = pt.y - pt0.y;vect[i].x = (float)dx;vect[i].y = (float)dy;inv_vect_length[i] = (float)(1. / std::sqrt(dx * dx + dy * dy));pt0 = pt;}// find convex hull orientation{double ax = vect[n - 1].x;double ay = vect[n - 1].y;for (i = 0; i < n; i++){double bx = vect[i].x;double by = vect[i].y;double convexity = ax * by - ay * bx;if (convexity != 0){orientation = (convexity > 0) ? 1.f : (-1.f);break;}ax = bx;ay = by;}assert(orientation != 0);}base_a = orientation;/*****************************************************************************************//* init calipers position */seq[0] = bottom;seq[1] = right;seq[2] = top;seq[3] = left;/*****************************************************************************************//* Main loop - evaluate angles and rotate calipers *//* all of edges will be checked while rotating calipers by 90 degrees */for (k = 0; k < n; k++){/* number of calipers edges, that has minimal angle with edge */int main_element = 0;/* choose minimum angle between calipers side and polygon edge by dot product sign */rot_vect[0] = vect[seq[0]];rotate90CW(vect[seq[1]], rot_vect[1]);rotate180(vect[seq[2]], rot_vect[2]);rotate90CCW(vect[seq[3]], rot_vect[3]);for (i = 1; i < 4; i++){if (firstVecIsRight(rot_vect[i], rot_vect[main_element]))main_element = i;}/*rotate calipers*/{//get next baseint pindex = seq[main_element];float lead_x = vect[pindex].x * inv_vect_length[pindex];float lead_y = vect[pindex].y * inv_vect_length[pindex];switch (main_element){case 0:base_a = lead_x;base_b = lead_y;break;case 1:base_a = lead_y;base_b = -lead_x;break;case 2:base_a = -lead_x;base_b = -lead_y;break;case 3:base_a = -lead_y;base_b = lead_x;break;default:assert("main_element should be 0, 1, 2 or 3" && false);}}/* change base point of main edge */seq[main_element] += 1;seq[main_element] = (seq[main_element] == n) ? 0 : seq[main_element];switch (mode){case RotateRect::CALIPERS_MAXHEIGHT:{/* now main element lies on edge aligned to calipers side *//* find opposite element i.e. transform *//* 0->2, 1->3, 2->0, 3->1 */int opposite_el = main_element ^ 2;float dx = points[seq[opposite_el]].x - points[seq[main_element]].x;float dy = points[seq[opposite_el]].y - points[seq[main_element]].y;float dist;if (main_element & 1)dist = (float)fabs(dx * base_a + dy * base_b);elsedist = (float)fabs(dx * (-base_b) + dy * base_a);if (dist > max_dist)max_dist = dist;}break;case RotateRect::CALIPERS_MINAREARECT:/* find area of rectangle */{float height;float area;/* find vector left-right */float dx = points[seq[1]].x - points[seq[3]].x;float dy = points[seq[1]].y - points[seq[3]].y;/* dotproduct */float width = dx * base_a + dy * base_b;/* find vector left-right */dx = points[seq[2]].x - points[seq[0]].x;dy = points[seq[2]].y - points[seq[0]].y;/* dotproduct */height = -dx * base_b + dy * base_a;area = width * height;if (area <= minarea){float* buf = (float*)buffer;minarea = area;/* leftist point */((int*)buf)[0] = seq[3];buf[1] = base_a;buf[2] = width;buf[3] = base_b;buf[4] = height;/* bottom point */((int*)buf)[5] = seq[0];buf[6] = area;}}break;} /*switch */} /* for */switch (mode){case RotateRect::CALIPERS_MINAREARECT:{float* buf = (float*)buffer;float A1 = buf[1];float B1 = buf[3];float A2 = -buf[3];float B2 = buf[1];float C1 = A1 * points[((int*)buf)[0]].x + points[((int*)buf)[0]].y * B1;float C2 = A2 * points[((int*)buf)[5]].x + points[((int*)buf)[5]].y * B2;float idet = 1.f / (A1 * B2 - A2 * B1);float px = (C1 * B2 - C2 * B1) * idet;float py = (A1 * C2 - A2 * C1) * idet;out[0] = px;out[1] = py;out[2] = A1 * buf[2];out[3] = B1 * buf[2];out[4] = A2 * buf[4];out[5] = B2 * buf[4];}break;case RotateRect::CALIPERS_MAXHEIGHT:{out[0] = max_dist;}break;}
}void RotateRect::update()
{using Point = RotateRect::Point;std::vector<Point> hull;Point out[3];hull = convexHull();int n = hull.size();const Point* hpoints = &hull[0];if (n > 2){rotatingCalipers(hpoints, n, RotateRect::CALIPERS_MINAREARECT, (float*)out);m_center.x = out[0].x + (out[1].x + out[2].x) * 0.5f;m_center.y = out[0].y + (out[1].y + out[2].y) * 0.5f;m_size.width = (float)std::sqrt((double)out[1].x * out[1].x + (double)out[1].y * out[1].y);m_size.height = (float)std::sqrt((double)out[2].x * out[2].x + (double)out[2].y * out[2].y);m_angle = (float)atan2((double)out[1].y, (double)out[1].x);}else if (n == 2){m_center.x = (hpoints[0].x + hpoints[1].x) * 0.5f;m_center.y = (hpoints[0].y + hpoints[1].y) * 0.5f;double dx = hpoints[1].x - hpoints[0].x;double dy = hpoints[1].y - hpoints[0].y;m_size.width = (float)std::sqrt(dx * dx + dy * dy);m_size.height = 0;m_angle = (float)atan2(dy, dx);}else{if (n == 1)m_center = hpoints[0];}m_angle = (float)(m_angle * 180 / CV_PI);
}RotateRect::Points RotateRect::toPoints() const
{RotateRect::Points pt(4);double _angle = m_angle * CV_PI / 180.;float b = (float)cos(_angle) * 0.5f;float a = (float)sin(_angle) * 0.5f;pt[0].x = m_center.x - a * m_size.height - b * m_size.width;pt[0].y = m_center.y + b * m_size.height - a * m_size.width;pt[1].x = m_center.x + a * m_size.height - b * m_size.width;pt[1].y = m_center.y - b * m_size.height - a * m_size.width;pt[2].x = 2 * m_center.x - pt[0].x;pt[2].y = 2 * m_center.y - pt[0].y;pt[3].x = 2 * m_center.x - pt[1].x;pt[3].y = 2 * m_center.y - pt[1].y;return pt;
}const RotateRect::Point& RotateRect::center() const {return m_center;
}const RotateRect::Sizef& RotateRect::size() const {return m_size;
}float RotateRect::angle() const {return m_angle;
}RotateRect::Rect RotateRect::getOutLine() const {if (m_inputs.empty())return { 0, 0, 0, 0 };using Number = std::numeric_limits<float>;RotateRect::Rect rect{ Number::max(), Number::max(), Number::min(), Number::min() };for (const auto& point : m_inputs) {if (point.x < rect.left)rect.left = point.x;if (point.x > rect.right)rect.right = point.x;if (point.y < rect.top)rect.top = point.y;if (point.y > rect.bottom)rect.bottom = point.y;}return rect;
}
main.cpp
RotateRect::Points myPoints = { {233, 86}, {322, 106}, {214, 154}, {307, 176}, {286, 209}, {331, 183}, {346, 319}, {392, 294}, {356, 346}, {419, 311}, {419, 311}, {778, 1031}, {840, 995} };RotateRect myRotateRect(myPoints);myRotateRect.update();