C++ 计算凸包点的最小旋转矩形

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::vectorPoint; 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 rectangle Sizef 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 #endif RotateRect::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 点 * return BA 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::vectorPoint 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::vectorfloat 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 base int 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); else dist (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::vectorPoint 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_limitsfloat; 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.cppRotateRect::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();创作不易小小的支持一下吧