#include "geometry/calipers_box.hpp" #include "geometry/bounding_box.hpp" #include "geometry/convex_hull.hpp" #include "geometry/line2d.hpp" #include "geometry/polygon.hpp" #include "geometry/segment2d.hpp" #include "base/assert.hpp" #include #include #include using namespace std; namespace m2 { namespace { static_assert(numeric_limits::has_infinity, ""); double const kInf = numeric_limits::infinity(); // Checks whether (p1 - p) x (p2 - p) >= 0. bool IsCCW(PointD const & p1, PointD const & p2, PointD const & p, double eps) { return robust::OrientedS(p1, p2, p) > -eps; } PointD Ort(PointD const & p) { return PointD(-p.y, p.x); } // For each facet of the |hull| calls |fn| with the smallest rectangle // containing the hull and with one side collinear to the facet. template void ForEachRect(ConvexHull const & hull, Fn && fn) { ASSERT_GREATER(hull.Size(), 2, ()); size_t j = 0, k = 0, l = 0; for (size_t i = 0; i < hull.Size(); ++i) { auto const ab = hull.SegmentAt(i).Dir(); j = max(j, i + 1); while (DotProduct(ab, hull.SegmentAt(j).Dir()) > CalipersBox::kEps) ++j; k = max(k, j); while (CrossProduct(ab, hull.SegmentAt(k).Dir()) > CalipersBox::kEps) ++k; l = max(l, k); while (DotProduct(ab, hull.SegmentAt(l).Dir()) < -CalipersBox::kEps) ++l; auto const oab = Ort(ab); array const lines = {{Line2D(hull.PointAt(i), ab), Line2D(hull.PointAt(j), oab), Line2D(hull.PointAt(k), ab), Line2D(hull.PointAt(l), oab)}}; vector corners; for (size_t i = 0; i < lines.size(); ++i) { auto const j = (i + 1) % lines.size(); auto result = LineIntersector::Intersect(lines[i], lines[j], CalipersBox::kEps); if (result.m_type == LineIntersector::Result::Type::One) corners.push_back(result.m_point); } if (corners.size() != 4) continue; auto const it = min_element(corners.begin(), corners.end()); rotate(corners.begin(), it, corners.end()); fn(corners); } } } // namespace CalipersBox::CalipersBox(vector const & points) : m_points({}) { ConvexHull hull(points, kEps); if (hull.Size() < 3) { m_points = hull.Points(); return; } double bestArea = kInf; vector bestPoints; ForEachRect(hull, [&](vector const & points) { ASSERT_EQUAL(points.size(), 4, ()); double const area = GetPolygonArea(points.begin(), points.end()); if (area < bestArea) { bestArea = area; bestPoints = points; } }); if (bestPoints.empty()) { BoundingBox bbox(points); auto const min = bbox.Min(); auto const max = bbox.Max(); auto const width = max.x - min.x; auto const height = max.y - min.y; m_points.resize(4); m_points[0] = min; m_points[1] = m_points[0] + PointD(width, 0); m_points[2] = m_points[1] + PointD(0, height); m_points[3] = m_points[0] + PointD(0, height); return; } ASSERT_EQUAL(bestPoints.size(), 4, ()); m_points = bestPoints; } bool CalipersBox::HasPoint(PointD const & p, double eps) const { auto const n = m_points.size(); if (n == 0) return false; if (n == 1) return AlmostEqualAbs(m_points[0], p, eps); if (n == 2) return IsPointOnSegmentEps(p, m_points[0], m_points[1], eps); for (size_t i = 0; i < n; ++i) { auto const & a = m_points[i]; auto const & b = m_points[(i + 1) % n]; if (!IsCCW(b, p, a, eps)) return false; } return true; } } // namespace m2