geometry/Caliper.cpp

• View this file on GitHub
• Last update: 2023-02-27 10:03:35-03:00

Depends on

• Epsilon.cpp
• bits/stdc++.h
• geometry/Circle2D.cpp
• geometry/GeometryEpsilon.cpp
• geometry/Line2D.cpp
• geometry/Polygon2D.cpp
• geometry/Trigonometry.cpp
• utils/Annotation.cpp

Code

#ifndef _LIB_GEOMETRY_CALIPER
#define _LIB_GEOMETRY_CALIPER
#include "Line2D.cpp"
#include "Polygon2D.cpp"
#include <bits/stdc++.h>

namespace lib {
using namespace std;
namespace geo {
namespace plane {
template <typename T, typename Large = T,
          typename enable_if<!is_integral<T>::value>::type * = nullptr,
          typename enable_if<!is_integral<T>::value>::type * = nullptr>
struct Caliper {
  typedef Point<T, Large> point;
  typedef Line<T, Large> line;
  point p;
  Large ang;
  Caliper(point a, Large alpha) : p(a) {
    ang = remainder(alpha, 2 * trig::PI);
    while (ang < 0)
      ang += 2 * trig::PI;
  }
  Large angle_to(const point &q) const {
    return remainder(arg(q - p) - ang, 2 * trig::PI);
  }
  void rotate(double theta) {
    ang += theta;
    while (ang > 2 * trig::PI)
      ang -= 2 * trig::PI;
    while (ang < 0)
      ang += 2 * trig::PI;
  }
  void move(const point &q) { p = q; }
  point versor() const { return point::polar(1.0, ang); }
  line as_line(Large scale = 1.0) const {
    return line(p, p + versor() * scale);
  }
  friend Large dist(const Caliper &a, const Caliper &b) {
    return dist(a.as_line(), b.p);
  }
};

template <typename T, typename Large = T> struct PolygonCalipers {
  constexpr static Large LIMIT = 4 * acosl(-1);

  typedef Point<T, Large> point;
  typedef Caliper<T, Large> caliper;
  typedef ConvexPolygon<T, Large> polygon;
  typedef pair<int, Large> descriptor;

  polygon poly;
  vector<caliper> calipers;
  vector<int> indices;
  vector<int> walked;
  Large angle_walked;

  PolygonCalipers(const polygon &poly, const vector<descriptor> &descriptors)
      : poly(poly), walked(descriptors.size()), angle_walked(0) {
    indices.reserve(descriptors.size());
    calipers.reserve(descriptors.size());
    for (size_t i = 0; i < descriptors.size(); i++) {
      calipers.emplace_back(poly[descriptors[i].first], descriptors[i].second);
      indices.emplace_back(descriptors[i].first);
    }
  }
  caliper operator[](int i) const { return calipers[i]; }
  int index(int i) const { return indices[i]; }
  bool has_next() const {
    return *min_element(walked.begin(), walked.end()) < poly.size() &&
           angle_walked < LIMIT;
  }
  Large angle_to_next(int i) const {
    int u = indices[i];
    return calipers[i].angle_to(poly[u + 1]);
  }
  void step_(int i) {
    int u = indices[i]++;
    indices[i] %= poly.size();
    calipers[i].move(poly[u + 1]);
    walked[i]++;
  }

  void next() {
    int i = 0;
    Large best = angle_to_next(0);
    for (size_t j = 1; j < calipers.size(); j++) {
      Large cur = angle_to_next(j);
      if (cur < best) {
        best = cur;
        i = j;
      }
    }
    Large alpha = angle_to_next(i);
    for (auto &caliper : calipers)
      caliper.rotate(alpha);
    step_(i);
    angle_walked += alpha;
  }
};
} // namespace plane

} // namespace geo
} // namespace lib

#endif

#line 1 "geometry/Caliper.cpp"


#line 1 "geometry/Line2D.cpp"


#line 1 "geometry/GeometryEpsilon.cpp"


#line 1 "Epsilon.cpp"


#include <bits/stdc++.h>

namespace lib {
using namespace std;

template <typename T = double> struct Epsilon {
  T eps;
  constexpr Epsilon(T eps = 1e-9) : eps(eps) {}

  template <typename G,
            typename enable_if<is_floating_point<G>::value>::type * = nullptr>
  int operator()(G a, G b = 0) const {
    return a + eps < b ? -1 : (b + eps < a ? 1 : 0);
  }

  template <typename G,
            typename enable_if<!is_floating_point<G>::value>::type * = nullptr>
  int operator()(G a, G b = 0) const {
    return a < b ? -1 : (a > b ? 1 : 0);
  }

  template <typename G,
            typename enable_if<is_floating_point<G>::value>::type * = nullptr>
  bool null(G a) const {
    return (*this)(a) == 0;
  }

  template <typename G,
            typename enable_if<!is_floating_point<G>::value>::type * = nullptr>
  bool null(G a) const {
    return a == 0;
  }
};
} // namespace lib


#line 5 "geometry/GeometryEpsilon.cpp"

#define GEOMETRY_EPSILON(T, x)                                                 \
  template <>                                                                  \
  lib::Epsilon<T> *lib::geo::GeometryEpsilon<T>::eps =                         \
      new lib::Epsilon<T>((x));

#define GEOMETRY_COMPARE0(T, x) GeometryEpsilon<T>()((x))
#define GEOMETRY_COMPARE(T, x, y) GeometryEpsilon<T>()((x), (y))

namespace lib {
using namespace std;
namespace geo {
template <typename T> struct GeometryEpsilon {
  static Epsilon<T> *eps;
  template <typename G> int operator()(G a, G b = 0) const {
    return (*eps)(a, b);
  }
};

GEOMETRY_EPSILON(int, 0);
GEOMETRY_EPSILON(long, 0);
GEOMETRY_EPSILON(long long, 0);
} // namespace geo
} // namespace lib


#line 1 "geometry/Trigonometry.cpp"


#line 4 "geometry/Trigonometry.cpp"

namespace lib {
using namespace std;
namespace geo {
namespace trig {
constexpr static long double PI = 3.141592653589793238462643383279502884197169399375105820974944l;
double cos(double x) { return ::cos(x); }
double sin(double x) { return ::sin(x); }
double asin(double x) { return ::asin(x); }
double acos(double x) { return ::acos(x); }
double atan2(double y, double x) { return ::atan2(y, x); }
long double cos(long double x) { return ::cosl(x); }
long double sin(long double x) { return ::sinl(x); }
long double asin(long double x) { return ::asinl(x); }
long double acos(long double x) { return ::acosl(x); }
long double atan2(long double y, long double x) { return ::atan2l(y, x); }
} // namespace trig
} // namespace geo
} // namespace lib


#line 6 "geometry/Line2D.cpp"

namespace lib {
using namespace std;
namespace geo {
namespace plane {
namespace {
template <typename T> bool scalar_between(T a, T o, T b) {
  if (a > b)
    swap(a, b);
  return GEOMETRY_COMPARE(T, a, o) <= 0 && GEOMETRY_COMPARE(T, o, b) <= 0;
}

template <typename T> bool scalar_strictly_between(T a, T o, T b) {
  if (a > b)
    swap(a, b);
  int x = GEOMETRY_COMPARE(T, a, o);
  int y = GEOMETRY_COMPARE(T, o, b);
  return x <= 0 && y <= 0 && (x < 0 || y < 0);
}
} // namespace

template <typename T, typename Large = T> struct Point {
  T x, y;
  Point() : x(0), y(0) {}
  Point(T x, T y) : x(x), y(y) {}
  template <typename G, typename H> explicit operator Point<G, H>() const {
    return Point<G, H>((G)x, (G)y);
  }
  friend Point reversed(const Point &a) { return Point(a.y, a.x); }
  Point &operator+=(const Point &rhs) {
    x += rhs.x, y += rhs.y;
    return *this;
  }
  Point &operator-=(const Point &rhs) {
    x -= rhs.x, y -= rhs.y;
    return *this;
  }
  Point &operator*=(T k) {
    x *= k, y *= k;
    return *this;
  }
  Point &operator/=(T k) {
    x /= k, y /= k;
    return *this;
  }
  Point operator+(const Point &rhs) const {
    Point res = *this;
    return res += rhs;
  }
  Point operator-(const Point &rhs) const {
    Point res = *this;
    return res -= rhs;
  }
  Point operator*(T k) const {
    Point res = *this;
    return res *= k;
  }
  Point operator/(T k) const {
    Point res = *this;
    return res /= k;
  }
  Point operator-() const { return Point(-x, -y); }
  inline friend Point convolve(const Point &a, const Point &b) {
    return Point(a.x * b.x - a.y * b.y, a.x * b.y + b.x * a.y);
  }
  inline friend Large cross(const Point &a, const Point &b) {
    return (Large)a.x * b.y - (Large)a.y * b.x;
  }
  friend Large cross(const Point &a, const Point &b, const Point &c) {
    return cross(b - a, c - a);
  }
  inline friend Large dot(const Point &a, const Point &b) {
    return (Large)a.x * b.x + (Large)a.y * b.y;
  }
  friend int ccw(const Point &u, const Point &v) {
    return GEOMETRY_COMPARE0(Large, cross(u, v));
  }
  friend int ccw(const Point &a, const Point &b, const Point &c) {
    return ccw(b - a, c - a);
  }
  friend int half_ccw(const Point& u, const Point& v) {
    int dot_sgn = GEOMETRY_COMPARE0(Large, dot(u, v));
    int ccw_sgn = ccw(u, v);
    if(dot_sgn == 0) return ccw_sgn ? 1 : 0;
    return dot_sgn * ccw_sgn;
  }
  friend Large norm(const Point &a) { return sqrtl(dot(a, a)); }
  friend Large norm_sq(const Point &a) { return dot(a, a); }
  bool is_null() const { return GEOMETRY_COMPARE0(Large, norm_sq(*this)) == 0; }
  bool is_versor() const {
    return GEOMETRY_COMPARE(Large, norm_sq(*this), (Large)1) == 0;
  }
  static Point polar(Large d, Large theta) {
    return Point(trig::cos(theta) * d, trig::sin(theta) * d);
  }
  friend Point rotate(const Point &a, Large theta) {
    return convolve(a, polar((Large)1, theta));
  }
  friend Point ortho(const Point &a) { return Point(-a.y, a.x); }
  friend Large arg(const Point &a) { return trig::atan2(a.y, a.x); }
  friend Large signed_angle(const Point &v, const Point &w) {
    return remainder(arg(w) - arg(v), 2.0 * trig::PI);
  }
  friend Large angle(const Point &v, const Point &w) {
    return abs(signed_angle(v, w));
  }
  friend Large ccw_angle(const Point &v) {
    Large res = arg(v);
    if (res < 0)
      res += 2.0 * trig::PI;
    return res;
  }
  friend Large ccw_angle(const Point &v, const Point &w) {
    Large res = signed_angle(v, w);
    if (res < 0)
      res += 2.0 * trig::PI;
    return res;
  }
  inline friend Point normalized(const Point &a, Large k) {
    return a.is_null() ? Point() : a / norm(a) * k;
  }
  inline friend Point versor(const Point &a) { return normalized(a, (Large)1); }
  friend bool collinear(const Point &a, const Point &b) {
    return GEOMETRY_COMPARE0(Large, cross(a, b)) == 0;
  }
  friend bool collinear(const Point &a, const Point &b, const Point &c) {
    return collinear(b - a, c - a);
  }
  friend Point project(const Point &a, const Point &v) {
    return v / norm_sq(v) * dot(a, v);
  }
  template <typename G = T,
            typename enable_if<!is_integral<G>::value>::type * = nullptr>
  friend Point reflect(const Point &a, const Point &v) {
    Point n = versor(v);
    return a - n * 2 * dot(n, v);
  }
  friend bool between(const Point &a, const Point &b, const Point &c) {
    return collinear(a, b, c) &&
           GEOMETRY_COMPARE0(Large, dot(a - b, c - b)) <= 0;
  }
  friend bool strictly_between(const Point &a, const Point &b, const Point &c) {
    return collinear(a, b, c) &&
           GEOMETRY_COMPARE0(Large, dot(a - b, c - b)) < 0;
  }
  friend bool collinear_between(const Point a, const Point &o, const Point &b) {
    return scalar_between(a.x, o.x, b.x) && scalar_between(a.y, o.y, b.y);
  }
  friend bool collinear_strictly_between(const Point &a, const Point &o,
                                         const Point &b) {
    return scalar_between(a.x, o.x, b.x) && scalar_between(a.y, o.y, b.y);
  }
  friend Large dist(const Point &a, const Point &b) { return norm(a - b); }
  friend bool operator==(const Point &a, const Point &b) {
    return GEOMETRY_COMPARE(T, a.x, b.x) == 0 &&
           GEOMETRY_COMPARE(T, a.y, b.y) == 0;
  }
  friend bool operator!=(const Point &a, const Point &b) { return !(a == b); }
  friend bool operator<(const Point &a, const Point &b) {
    return tie(a.y, a.x) < tie(b.y, b.x);
  }
  friend bool operator>(const Point &a, const Point &b) {
    return tie(a.y, a.x) > tie(b.y, b.x);
  }
  friend bool operator>=(const Point &a, const Point &b) {
    return tie(a.y, a.x) >= tie(b.y, b.x);
  }
  friend bool operator<=(const Point &a, const Point &b) {
    return tie(a.y, a.x) <= tie(b.y, b.x);
  }
  friend istream &operator>>(istream &in, Point &p) { return in >> p.x >> p.y; }
  friend ostream &operator<<(ostream &out, const Point &p) {
    return out << p.x << " " << p.y;
  }
};

template <typename T, typename Large = T> struct Rectangle {
  typedef Point<T, Large> point;

  T minx, miny, maxx, maxy;
  Rectangle() {
    minx = miny = numeric_limits<T>::max();
    maxx = maxy = numeric_limits<T>::min();
  }

  Rectangle(const initializer_list<point> &points) : Rectangle() {
    for (const auto &p : points) {
      minx = min(minx, p.x);
      maxx = max(maxx, p.x);
      miny = min(miny, p.y);
      maxy = max(maxy, p.y);
    }
  }

  bool contains(const point &p) const {
    return GEOMETRY_COMPARE(T, minx, p.x) <= 0 &&
           GEOMETRY_COMPARE(T, p.x, maxx) <= 0 &&
           GEOMETRY_COMPARE(T, miny, p.y) <= 0 &&
           GEOMETRY_COMPARE(T, p.y, maxy) <= 0;
  }
};

template <typename T, typename Large = T> struct Line {
  typedef Point<T, Large> point;
  typedef Line<T, Large> line;
  point a, b;
  Line(point a, point b) : a(a), b(b) {}
  template <typename G = T,
            typename enable_if<!is_integral<G>::value>::type * = nullptr>
  Line(T A, T B, T C) {
    if (GEOMETRY_COMPARE0(Large, A))
      a = point(-C / A, 0), b = point((-C - B) / A, 1);
    else if (GEOMETRY_COMPARE0(Large, B))
      a = point(0, -C / B), b = point(1, (-C - A) / B);
    else
      assert(false);
  }
  template <typename G, typename H> explicit operator Line<G, H>() const {
    return Line<G, H>(Point<G, H>(a), Point<G, H>(b));
  }
  point direction() const { return b - a; }
  friend point project(const point &p, const line &v) {
    return project(p - v.a, v.b - v.a) + v.a;
  }
  friend bool collinear(const line &u, const line &v) {
    return collinear(u.a, u.b, v.a) && collinear(u.a, u.b, v.b);
  }
  bool contains(const point &p) const { return collinear(a, b, p); }
  friend bool parallel(const line &u, const line &v) {
    return collinear(u.b - u.a, v.b - v.a);
  }
  friend bool opposite(const line &l, const point &p1, const point &p2) {
    int x = GEOMETRY_COMPARE0(Large, cross(p1 - l.a, l.direction()));
    int y = GEOMETRY_COMPARE0(Large, cross(p2 - l.a, l.direction()));
    return x * y <= 0;
  }
  friend pair<point, bool> intersect(const line &l1, const line &l2) {
    Large c1 = cross(l2.a - l1.a, l1.b - l1.a);
    Large c2 = cross(l2.b - l1.a, l1.b - l1.a);
    if (GEOMETRY_COMPARE0(Large, c1 - c2) == 0)
      return {{}, false};
    return {(l2.b * c1 - l2.a * c2) / (c1 - c2), true};
  }
  friend bool has_unique_intersection(const line &l1, const line &l2) {
    return !parallel(l1, l2);
  }
  friend bool has_intersection(const line &l1, const line &l2) {
    return collinear(l1, l2) || has_unique_intersection(l1, l2);
  }
  friend Large dist(const line &l1, const point &p) {
    // TODO: improve this
    return dist(p, project(p, l1));
  }
  friend Large dist(const line &l1, const line &l2) {
    if (has_intersection(l1, l2))
      return 0;
    // TODO: improve this
    return dist(l1.a, project(l1.a, l2));
  }
};

template <typename T, typename Large = T> struct Ray {
  typedef Point<T, Large> point;
  typedef Line<T, Large> line;
  typedef Ray<T, Large> ray;
  point a, b;

  Ray(point a, point direction) : a(a), b(a + direction) {}

  static ray from_points(point a, point b) { return ray(a, b - a); }
  point direction() const { return b - a; }
  point direction_versor() const { return versor(direction()); }

  line as_line() const { return line(a, b); }
  explicit operator line() const { return as_line(); }

  template <typename G, typename H> explicit operator Ray<G, H>() const {
    return Ray<G, H>(Point<G, H>(a), Point<G, H>(b));
  }
  bool contains(const point &p) const {
    return collinear(a, b, p) &&
           GEOMETRY_COMPARE0(Large, dot(p - a, b - a)) >= 0;
  }
  bool strictly_contains(const point &p) const {
    return collinear(a, b, p) &&
           GEOMETRY_COMPARE0(Large, dot(p - a, b - a)) > 0;
  }
  bool collinear_contains(const point &p) const {
    point dir = direction();
    int dx = GEOMETRY_COMPARE0(T, dir.x);
    if (dx == 0)
      return GEOMETRY_COMPARE0(T, dir.y) * GEOMETRY_COMPARE0(T, p.y - a.y) >= 0;
    else
      return dx * GEOMETRY_COMPARE0(T, p.x - a.x) >= 0;
  }
  bool collinear_strictly_contains(const point &p) const {
    point dir = direction();
    int dx = GEOMETRY_COMPARE0(T, dir.x);
    if (dx == 0)
      return GEOMETRY_COMPARE0(T, dir.y) * GEOMETRY_COMPARE0(T, p.y - a.y) > 0;
    else
      return dx * GEOMETRY_COMPARE0(T, p.x - a.x) > 0;
  }
  friend pair<point, bool> intersect(const ray &r, const line &l) {
    auto p = intersect(r.as_line(), l);
    if (!p.second)
      return {{}, false};
    if (!r.collinear_contains(p.first))
      return {{}, false};
    return p;
  }
  friend pair<point, bool> intersect(const ray &a, const ray &b) {
    auto p = intersect(a, b.as_line());
    if (!p.second)
      return {{}, false};
    if (!b.collinear_contains(p.first))
      return {{}, false};
    return p;
  }
  friend bool has_unique_intersection(const ray &r, const line &l) {
    if (!has_unique_intersection(r.as_line(), l))
      return false;
    int x = GEOMETRY_COMPARE0(Large, cross(r.direction(), l.direction()));
    int y = GEOMETRY_COMPARE0(Large, cross(r.a - l.a, l.direction()));
    return x * y <= 0;
  }
  friend bool has_intersection(const ray &r, const line &l) {
    return collinear(r.as_line(), l) || has_unique_intersection(r, l);
  }
  friend bool has_unique_intersection(const ray &r1, const ray &r2) {
    // TODO: not efficient
    return has_unique_intersection(r1, r2.as_line()) &&
           has_unique_intersection(r2, r1.as_line());
  }
  friend bool has_intersection(const ray &r1, const ray &r2) {
    return r1.contains(r2.a) || has_unique_intersection(r1, r2);
  }
  friend Large dist(const ray &r, const point &p) {
    if (GEOMETRY_COMPARE0(Large, dot(r.direction(), p - r.a)) < 0)
      return dist(p, r.a);
    return dist(r.as_line(), p);
  }
  friend Large dist(const ray &r, const line &l) {
    if (has_intersection(r, l))
      return Large(0);
    return dist(l, r.a);
  }
  friend Large dist(const ray &r1, const ray &r2) {
    if (has_intersection(r1, r2))
      return Large(0);
    return min(dist(r1, r2.a), dist(r2, r1.a));
  }
};

template <typename T, typename Large = T> struct Halfplane {
  typedef Point<T, Large> point;
  typedef Line<T, Large> line;
  typedef Ray<T, Large> ray;
  typedef Halfplane<T, Large> halfplane;
  point a, b;

  Halfplane(point a, point direction) : a(a), b(a + direction) {}

  static halfplane from_points(point a, point b) { return halfplane(a, b - a); }
  point direction() const { return b - a; }
  point direction_versor() const { return versor(direction()); }

  line as_line() const { return line(a, b); }
  explicit operator line() const { return as_line(); }

  ray as_ray() const { return ray(a, b); }
  explicit operator ray() const { return as_ray(); }

  template <typename G, typename H> explicit operator Halfplane<G, H>() const {
    return Halfplane<G, H>(Point<G, H>(a), Point<G, H>(b));
  }

  bool contains(const point& p) const {
    return ccw(a, b, p) <= 0;
  }
  bool strictly_contains(const point& p) const {
    return ccw(a, b, p) < 0;
  }
};

template <typename T, typename Large = T> struct Segment {
  typedef Point<T, Large> point;
  typedef Line<T, Large> line;
  typedef Segment<T, Large> segment;
  typedef Ray<T, Large> ray;
  point a, b;

  Segment() {}
  Segment(point a, point b) : a(a), b(b) {}
  line as_line() const { return line(a, b); }
  explicit operator line() const { return as_line(); }
  bool is_degenerate() const { return a == b; }

  template <typename G, typename H> explicit operator Segment<G, H>() const {
    return Segment<G, H>(Point<G, H>(a), Point<G, H>(b));
  }
  bool contains(const point &p) const { return between(a, p, b); }
  bool strictly_contains(const point &p) const {
    return strictly_between(a, p, b);
  }
  bool collinear_contains(const point &p) const {
    return collinear_between(a, p, b);
  }
  bool collinear_strictly_contains(const point &p) const {
    return collinear_strictly_between(a, p, b);
  }
  friend pair<point, bool> intersect(const segment &s, const line &l) {
    auto p = intersect(s.as_line(), l);
    if (!p.second)
      return {{}, false};
    if (!s.collinear_contains(p.first))
      return {{}, false};
    return p;
  }
  friend pair<point, bool> intersect(const segment &s, const ray &r) {
    auto p = intersect(s.as_line(), r.as_line());
    if (!p.second)
      return {{}, false};
    if (!s.collinear_contains(p.first) || !r.collinear_contains(p.first))
      return {{}, false};
    return p;
  }
  friend pair<segment, int> intersect_segment(segment s1, segment s2) {
    if (collinear(s1.as_line(), s2.as_line())) {
      if (s1.a > s1.b)
        swap(s1.a, s1.b);
      if (s2.a > s2.b)
        swap(s2.a, s2.b);
      segment res(max(s1.a, s2.a), min(s1.b, s2.b));
      return {res, int(res.a <= res.b) * 2};
    } else {
      auto p = intersect(s1, s2);
      return {segment(p.first, p.first), p.second};
    }
  }
  friend pair<point, bool> intersect(const segment &s1, const segment &s2) {
    auto p = intersect(s1, s2.as_line());
    if (!p.second)
      return {{}, false};
    if (!s2.collinear_contains(p.first))
      return {{}, false};
    return p;
  }
  friend bool has_unique_intersection(const segment &s, const line &l) {
    if (!has_unique_intersection(s.as_line(), l))
      return false;
    return opposite(l, s.a, s.b);
  }
  friend bool has_intersection(const segment &s, const line &l) {
    return collinear(s.as_line(), l) || has_unique_intersection(s, l);
  }
  friend bool has_unique_intersection(const segment &s, const ray &r) {
    if (!has_unique_intersection(r, s.as_line()))
      return false;
    return opposite(r.as_line(), s.a, s.b);
  }
  friend bool has_intersection(const segment &s, const ray &r) {
    return r.contains(s.a) || r.contains(s.b) || has_unique_intersection(s, r);
  }
  friend bool has_unique_intersection(const segment &s1, const segment &s2) {
    if (!has_unique_intersection(s1.as_line(), s2.as_line()))
      return false;
    return opposite(s2.as_line(), s1.a, s1.b) &&
           opposite(s1.as_line(), s2.a, s2.b);
  }
  friend bool has_intersection(const segment &s1, const segment &s2) {
    return s1.contains(s2.a) || s1.contains(s2.b) ||
           has_unique_intersection(s1, s2);
  }
  friend Large dist(const segment &s, const point &p) {
    if (GEOMETRY_COMPARE0(Large, dot(p - s.a, s.b - s.a)) <= 0)
      return dist(s.a, p);
    if (GEOMETRY_COMPARE0(Large, dot(p - s.b, s.a - s.b)) <= 0)
      return dist(s.b, p);
    return dist(s.as_line(), p);
  }
  friend Large dist(const segment &s, const line &l) {
    if (has_intersection(s, l))
      return Large(0);
    return min(dist(l, s.a), dist(l, s.b));
  }
  friend Large dist(const segment &s, const ray &r) {
    if (has_intersection(s, r))
      return Large(0);
    return min({dist(r, s.a), dist(r, s.b), dist(s, r.a)});
  }
  friend Large dist(const segment &s1, const segment &s2) {
    if (has_intersection(s1, s2))
      return Large(0);
    return min(
        {dist(s1, s2.a), dist(s1, s2.b), dist(s2, s1.a), dist(s2, s1.b)});
  }

  friend bool operator==(const segment &l1, const segment &l2) {
    return tie(l1.a, l1.b) == tie(l2.a, l2.b);
  }
  friend bool operator!=(const segment &l1, const segment &l2) {
    return !(l1 == l2);
  }
  friend bool operator<(const segment &l1, const segment &l2) {
    return tie(l1.a, l1.b) < tie(l2.a, l2.b);
  }
};

template <typename Direction, typename T, typename Large> struct AngleComparator {
  using type = typename Direction::type;
  using point = Point<T, Large>;

  Direction dir;
  AngleComparator() {}
  AngleComparator(Direction dir) : dir(dir) {}
  bool operator()(const type &a, const type &b) const {
    return ccw(dir(a), dir(b)) > 0;
  }
  template <typename Iterator>
  static void sortByAngle(Iterator begin, Iterator end, const Direction& dir = Direction()) {
    AngleComparator cmp(dir);
    begin =
        partition(begin, end, [&dir](const type &p) { return dir(p).is_null(); });
    auto half =
        partition(begin, end, [&dir](const type &p) { return dir(p) > point(); });
    sort(begin, half, cmp);
    sort(half, end, cmp);
  }
  template <typename Iterator>
  static Iterator minByAngle(Iterator begin, Iterator end, const Direction& dir = Direction()) {
    AngleComparator cmp(dir);
    return min_element(begin, end, [&dir, &cmp](const type& a, const type& b) {
      bool part_a = dir(a) > point();
      bool part_b = dir(b) > point();
      if(part_a == part_b)
        return cmp(a, b);
      return part_a > part_b;
    });
  }
};
template <typename Ray> struct RayDirection {
  using point = typename Ray::point;
  using type = Ray;
  point operator()(const type& rhs) const {
    return rhs.direction();
  }
};
template <typename Point> struct PointDirection {
  using type = Point;
  Point pivot;
  PointDirection() : pivot() {}
  PointDirection(Point pivot) : pivot(pivot) {}
  Point operator()(const Point& rhs) const {
    return (rhs - pivot).direction();
  }
};
} // namespace plane

template <typename T, typename Large = T> struct CartesianPlane {
  typedef plane::Point<T, Large> point;
  typedef plane::Line<T, Large> line;
  typedef plane::Rectangle<T, Large> rectangle;
  typedef plane::Segment<T, Large> segment;
  typedef plane::Ray<T, Large> ray;
  typedef plane::Halfplane<T, Large> halfplane;

  template<typename Direction>
  using angle_comparator = plane::AngleComparator<Direction, T, Large>;
};

} // namespace geo
} // namespace lib


#line 1 "geometry/Polygon2D.cpp"


#line 1 "geometry/Circle2D.cpp"


#line 1 "utils/Annotation.cpp"


#line 4 "utils/Annotation.cpp"

namespace lib {
using namespace std;
template <typename T, typename A = void>
struct Note : public T {
private:
    A m_data = A();
    Note(const T& t, const A& a) : T(t), m_data(a) {}
public:
    using T::T;

    static Note make(const T& t, const A& a) {
        return Note(t, a);
    }

    friend A& annotation(Note& note) {
        return note.m_data;
    }
    friend const A& annotation(const Note& note) {
        return note.m_data;
    }

    template<typename C, typename D>
    operator Note<T,A>() const {
        return Note<C, D>(*this, m_data);
    }
};

template <typename T>
struct Note<T, void> : public T {
    using T::T;
    using T::operator=;
    
    Note(const T& a) : T(a) {}
    Note(T &&a): T(std::move(a)) {}
};

template<typename T, typename A>
Note<T, A> make_note(const T& t, const A& a) {
    return Note<T, A>::make(t, a);
}
} // namespace lib


#line 6 "geometry/Circle2D.cpp"

namespace lib {
using namespace std;
namespace geo {
namespace plane {
template <typename T, typename Large = T> struct Barycentric {
  typedef Point<T, Large> point;
  point r1, r2, r3;
  T a, b, c;

  Barycentric(const point &r1, const point &r2, const point &r3, T a = 1,
              T b = 1, T c = 1)
      : r1(r1), r2(r2), r3(r3), a(a), b(b), c(c) {}
  point as_point() const { return (r1 * a + r2 * b + r3 * c) / (a + b + c); }

  static Barycentric centroid(const point &r1, const point &r2,
                              const point &r3) {
    return Barycentric(r1, r2, r3);
  }
  static Barycentric circumcenter(const point &r1, const point &r2,
                                  const point &r3) {
    Large a = norm_sq(r2 - r3), b = norm_sq(r3 - r1), c = norm_sq(r1 - r2);
    return Barycentric(r1, r2, r3, a * (b + c - a), b * (c + a - b),
                       c * (a + b - c));
  }
  static Barycentric incenter(const point &r1, const point &r2,
                              const point &r3) {
    return Barycentric(r1, r2, r3, norm(r2 - r3), norm(r1 - r3), norm(r1 - r2));
  }
  static Barycentric orthocenter(const point &r1, const point &r2,
                                 const point &r3) {
    Large a = norm_sq(r2 - r3), b = norm_sq(r3 - r1), c = norm_sq(r1 - r2);
    return Barycentric(r1, r2, r3, (a + b - c) * (c + a - b),
                       (b + c - a) * (a + b - c), (c + a - b) * (b + c - a));
  }
  static Barycentric excenter(const point &r1, const point &r2,
                              const point &r3) {
    return Barycentric(r1, r2, r3, -norm(r2 - r3), norm(r1 - r3),
                       norm(r1 - r2));
  }
};

template <typename T, typename Large = T> struct Circle {
  typedef Point<T, Large> point;
  typedef Line<T, Large> line;
  typedef Barycentric<Large> bary;
  typedef Segment<T, Large> segment;
  point center;
  T radius;

  Circle(point center, T radius) : center(center), radius(radius) {}
  Circle(const point &p1, const point &p2, const point &p3) {
    center = bary::circumcenter(p1, p2, p3).as_point();
    radius = dist(center, p1);
  }
  Circle(const point &p1, const point &p2) {
    center = (p1 + p2) / 2;
    radius = dist(center, p1);
  }
  bool crosses_x_axis(point p = point()) const {
    auto c = center - p;
    return GEOMETRY_COMPARE0(T, c.y + radius) >= 0 && GEOMETRY_COMPARE0(T, c.y - radius) < 0;
  }
  static Circle incircle(const point &p1, const point &p2, const point &p3) {
    point center = bary::incenter(p1, p2, p3).as_point();
    return Circle(center, dist(line(p1, p2), center));
  }
  friend pair<segment, int> intersect_segment(const Circle &c, const line &l) {
    point H = project(c.center, l);
    Large h = norm(H - c.center);
    if (GEOMETRY_COMPARE(Large, c.radius, h) < 0)
      return {{}, 0};
    Large norma = sqrtl(c.radius + h) * sqrtl(c.radius - h);
    point v = normalized(l.direction(), norma);
    segment res = segment(H - v, H + v);
    return {res, res.is_degenerate() ? 1 : 2};
  }
  friend Large intersection_area(const Circle &a, const Circle &b) {
    Large d = norm(a.center - b.center);
    if (GEOMETRY_COMPARE(Large, a.radius + b.radius, d) <= 0)
      return 0.0;
    if (GEOMETRY_COMPARE(Large, d, abs(a.radius - b.radius)) <= 0) {
      T r = min(a.radius, b.radius);
      return r * r * trig::PI;
    }

    auto compute = [d](Large ra, Large rb) {
      Large sup = rb * rb + d * d - ra * ra;
      Large alpha = trig::acos(sup / (2.0 * rb * d));
      Large s = alpha * rb * rb;
      Large t = rb * rb * trig::sin(alpha) * trig::cos(alpha);
      return s - t;
    };
    return compute(a.radius, b.radius) + compute(b.radius, a.radius);
  }
  static Large intersection_signed_area(T r, const point &a, const point &b) {
    Circle C(point(), r);
    auto ps = intersect_segment(C, line(a, b));
    if (!ps.second)
      return r * r * signed_angle(a, b) / 2;
    auto s = ps.first;
    bool outa = !contains(C, a), outb = !contains(C, b);
    if (outa && outb) {
      segment ab(a, b);
      if (ab.contains(s.a) && ab.contains(s.b))
        return (r * r * (signed_angle(a, b) - signed_angle(s.a, s.b)) +
                cross(s.a, s.b)) /
               2;
      return r * r * signed_angle(a, b) / 2;
    } else if (outa)
      return (r * r * signed_angle(a, s.a) + cross(s.a, b)) / 2;
    else if (outb)
      return (r * r * signed_angle(s.b, b) + cross(a, s.b)) / 2;
    else
      return cross(a, b) / 2;
  }
  friend vector<point> tangents(const Circle &C, const point &p) {
    return _tangents({p, T()}, C, {1});
  }
  friend vector<line> inner_tangents(const Circle& a, const Circle& b) {
    return _tangents(a, b, {-1});
  }
  friend vector<line> outer_tangents(const Circle& a, const Circle& b) {
    return _tangents(a, b, {1});
  }
  friend vector<line> _tangents(const Circle& a, const Circle& b, const initializer_list<int>& r_sgn) {
    vector<line> res;
    for(int r_s : r_sgn) {
      point d = b.center - a.center;
      Large dr = (a.radius - b.radius*r_s), d2 = norm_sq(d), h2 = d2 - dr*dr;
      if(GEOMETRY_COMPARE0(Large, d2) == 0) continue;
      if(GEOMETRY_COMPARE0(Large, h2) < 0) continue;
      for(T sgn : {-1, 1}) {
        point v = (d * dr + ortho(d) * sqrtl(h2) * sgn) / d2;
        res.push_back({a.center + v * a.radius, b.center + v * (b.radius * r_s)});
      }
      if(GEOMETRY_COMPARE0(Large, h2) == 0) res.pop_back();
    }
    return res;
  }
  friend vector<Note<line, int>> angular_tangents(const Circle& a, const vector<Circle>& v, vector<int>& sgn) {
    vector<Note<line, int>> res;
    res.reserve(4 * v.size());
    int i = 0;
    sgn = vector<int>(v.size());
    vector<bool> reversed(4);
    bool null_a = GEOMETRY_COMPARE0(T, a.radius) == 0;

    for(int i = 0; i < v.size(); i++) {
      bool null_i = GEOMETRY_COMPARE0(T, v[i].radius) == 0;
      assert(!null_a || !null_i);
      vector<line> tgts;
      if(null_a || null_i) tgts = _tangents(a, v[i], {1});
      else tgts = _tangents(a, v[i], {+1, -1});
      if(tgts.empty()) continue;

      fill(reversed.begin(), reversed.end(), false);
      int j = 0;
      for(auto& t : tgts) {
        // direct tangents
        if(ccw(t.b - t.a, a.center - t.a) < 0)
          swap(t.a, t.b), reversed[j] = true;
        res.push_back(make_note<line, int>(t, i));
        j++;
      }

      // check signal
      auto it = AngleComparator<RayDirection<line>, T, Large>::minByAngle(tgts.begin(), tgts.end());
      point ta = reversed[it - tgts.begin()] ? it->b : it->a;
      point dir = v[i].center - ta;
      sgn[i] = half_ccw(it->direction(), dir);
    }
    AngleComparator<RayDirection<line>, T, Large>::sortByAngle(res.begin(), res.end());
    return res;
  }
  friend bool contains(const Circle &c, const point &p) {
    return GEOMETRY_COMPARE(Large, dist(p, c.center), c.radius) <= 0;
  }
  friend bool contains(const Circle &c, const segment &s) {
    return GEOMETRY_COMPARE(Large, dist(s.a, c.center), c.radius) <= 0 &&
           GEOMETRY_COMPARE(Large, dist(s.b, c.center), c.radius) <= 0;
  }
  template <typename L>
  friend bool partially_contains(const Circle &c, const L &l) {
    return GEOMETRY_COMPARE(Large, dist(l, c.center), c.radius) <= 0;
  }
  template <typename L>
  friend bool has_unique_intersection(const Circle &c, const L &l) {
    return GEOMETRY_COMPARE(Large, dist(l, c.center), c.radius) == 0;
  }
  template <typename L>
  friend bool has_intersection(const Circle &c, const L &l) {
    return GEOMETRY_COMPARE(Large, dist(l, c.center), c.radius) <= 0;
  }
  friend bool has_intersection(const Circle &c, const segment &s) {
    return GEOMETRY_COMPARE(Large, dist(s, c.center), c.radius) <= 0 &&
           (GEOMETRY_COMPARE(Large, dist(s.a, c.center), c.radius) >= 0 ||
            GEOMETRY_COMPARE(Large, dist(s.b, c.center), c.radius) >= 0);
  }
};
} // namespace plane

template <typename T, typename Large = T>
struct CirclePlane : public CartesianPlane<T, Large> {
  typedef plane::Circle<T, Large> circle;
};

} // namespace geo
} // namespace lib


#line 6 "geometry/Polygon2D.cpp"

namespace lib {
using namespace std;
namespace geo {
namespace plane {

template <typename T, typename Large = T> struct ConvexHullComparator {
  typedef Point<T, Large> point;
  point pivot;
  ConvexHullComparator(point p) : pivot(p) {}
  template <typename G>
  bool operator()(const pair<point, G> &a, const pair<point, G> &b) const {
    int k = ccw(pivot, a.first, b.first);
    if (k == 0)
      return norm_sq(a.first) < norm_sq(b.first);
    return k > 0;
  }
};

template <typename T, typename Large = T> struct Polygon {
  typedef Point<T, Large> point;
  typedef Polygon<T, Large> polygon;
  typedef Circle<T, Large> circle;
  vector<point> p;

  Polygon() {}
  Polygon(const vector<point> &p) : p(p) {}
  template <typename G> Polygon(const vector<pair<point, G>> &g) : p(g.size()) {
    for (size_t i = 0; i < g.size(); i++)
      p[i] = g[i].first;
  }
  template <typename A, typename B> explicit operator Polygon<A, B>() const {
    vector<Point<A, B>> v(p.size());
    for (size_t i = 0; i < p.size(); i++)
      v[i] = Point<A, B>(p[i]);
    return Polygon<A, B>(v);
  }
  inline int index(int i) const {
    if (i >= size())
      i %= size();
    else if (i < 0) {
      i %= size();
      if (i < 0)
        i += size();
    }
    return i;
  }
  inline int size() const { return p.size(); }
  inline point &operator[](int i) { return p[index(i)]; }
  inline point operator[](int i) const { return p[index(i)]; }
  void erase(int i) { p.erase(p.begin() + index(i)); }
  polygon &operator+=(const point &pt) {
    for (auto &q : p)
      q += pt;
    return *this;
  }
  polygon &operator-=(const point &pt) {
    for (auto &q : p)
      q -= pt;
    return *this;
  }
  polygon &operator*=(const Large k) {
    for (auto &q : p)
      q *= k;
    return *this;
  }
  polygon &operator/=(const Large k) {
    for (auto &q : p)
      q /= k;
    return *this;
  }
  polygon operator-() const {
    polygon res = *this;
    for (auto &q : res.p)
      q = -q;
    return res;
  }
  void reserve(int n) { p.reserve(n); }
  bool is_ccw() const {
    int n = size();
    int i = min_element(p.begin(), p.end()) - p.begin();
    return ccw(p[i], p[i + 1], p[i - 1]) >= 0;
  }
  bool is_degenerate() const {
    int n = size();
    if (n < 3)
      return true;
    for (int i = 0; i < n; i++) {
      if (GEOMETRY_COMPARE0(Large, cross(p[i + 2] - p[i], p[i + 1] - p[i])))
        return false;
    }
    return true;
  }
  inline operator vector<point>() const { return p; }

  friend Large double_area(const Polygon &p) {
    int n = p.size();
    Large res = 0;
    for (int i = 0; i < n; i++) {
      res += cross(p[i], p[i + 1]);
    }
    return abs(res);
  }
  friend Large area(const Polygon &p) { return double_area(p) / 2; }
  friend Large perimeter(const Polygon &p) {
    int n = p.size();
    Large res = 0;
    for (int i = 0; i < n; i++)
      res += dist(p[i], p[i + 1]);
    return res;
  }

  int test(const point &p) const {
    const Polygon &poly = *this;
    int n = size();
    int wn = 0;
    for (int i = 0; i < n; i++) {
      if (p == poly[i])
        return 0;
      int j = i + 1;
      if (poly[i].y == p.y && poly[j].y == p.y) {
        if (min(poly[i].x, poly[j].x) <= p.x &&
            p.x <= max(poly[i].x, poly[j].x))
          return 0;
      } else {
        bool below = poly[i].y < p.y;
        if (below != (poly[j].y < p.y)) {
          auto sig = ccw(poly[i], poly[j], p);
          if (sig == 0)
            return 0;
          if (below == (sig > 0))
            wn += below ? 1 : -1;
        }
      }
    }
    return wn == 0 ? 1 : -1;
  }

  template <typename G>
  static vector<pair<point, G>> convex_hull(vector<pair<point, G>> p,
                                            bool keep_border = false) {
    if (p.size() <= 1)
      return p;
    sort(p.begin(), p.end());
    vector<pair<point, G>> res;
    res.reserve(p.size() + 1);
    for (int step = 0; step < 2; step++) {
      auto start = res.size();
      for (auto &q : p) {
        while (res.size() >= start + 2) {
          int sig = ccw(res[res.size() - 2].first, res.back().first, q.first);
          if ((sig == 0 && !keep_border) || sig < 0)
            res.pop_back();
          else
            break;
        }
        res.push_back(q);
      }
      res.pop_back();
      if (step == 0)
        reverse(p.begin(), p.end());
    }
    if (res.size() == 2 && res[0] == res[1])
      res.pop_back();
    return res;
  }

  static polygon convex_hull(const vector<point> &p, bool keep_border = false) {
    vector<pair<point, int>> v(p.size());
    for (size_t i = 0; i < p.size(); i++)
      v[i] = {p[i], i};
    auto res = convex_hull(v, keep_border);
    return polygon(res);
  }

  friend vector<polygon> triangulation(polygon poly) {
    if (poly.size() < 3)
      return {};
    vector<polygon> res;
    int ptr = 0;
    int n;
    while ((n = poly.size()) > 3) {
      for (int &i = ptr;; i++) {
        if (ccw(poly[i - 1], poly[i], poly[i + 1]) > 0) {
          auto trig = polygon({poly[i - 1], poly[i], poly[i + 1]});
          bool good = true;
          for (int j = 0; j < n; j++) {
            good &= trig.test(poly[j]) >= 0;
          }
          if (!good)
            continue;
          poly.erase(i--);
          res.push_back(trig);
          break;
        }
      }
    }
    res.push_back(poly);
    return res;
  }

  friend Large intersection_area(const Polygon &p, const circle &C) {
    Large res = 0;
    int n = p.size();
    for (int i = 0; i < n; i++) {
      res += circle::intersection_signed_area(C.radius, p[i + 1] - C.center,
                                              p[i] - C.center);
    }
    return abs(res);
  }
};

template <typename T, typename Large = T>
struct ConvexPolygon : public Polygon<T, Large> {
  typedef Point<T, Large> point;
  typedef Segment<T, Large> segment;
  typedef Line<T, Large> line;
  typedef Halfplane<T, Large> halfplane;
  typedef Circle<T, Large> circle;
  typedef AngleComparator<PointDirection<point>, T, Large> angle_comparator;
  using Polygon<T, Large>::p;
  int top;
  ConvexPolygon() {}
  ConvexPolygon(const vector<point> &p) : Polygon<T, Large>(p) { normalize(); }
  template <typename G>
  ConvexPolygon(const vector<pair<point, G>> &p) : Polygon<T, Large>(p) {
    normalize();
  }
  void normalize() {
    auto bottom = min_element(p.begin(), p.end());
    rotate(p.begin(), bottom, p.end());
    top = max_element(p.begin(), p.end()) - p.begin();
  }
  ConvexPolygon &operator+=(const point &pt) {
    for (auto &q : p)
      q += pt;
    return *this;
  }
  ConvexPolygon &operator-=(const point &pt) {
    for (auto &q : p)
      q -= pt;
    return *this;
  }
  ConvexPolygon &operator*=(const Large k) {
    for (auto &q : p)
      q *= k;
    return *this;
  }
  ConvexPolygon &operator/=(const Large k) {
    for (auto &q : p)
      q /= k;
    return *this;
  }
  ConvexPolygon operator-() const {
    ConvexPolygon res = *this;
    for (auto &q : res.p)
      q = -q;
    return res;
  }

  int test(const point &q) const {
    if (q < p[0] || q > p[top])
      return 1;
    auto sig = ccw(p[0], p[top], q);
    if (sig == 0) {
      if (q == p[0] || q == p[top])
        return 0;
      return top == 1 || top + 1 == this->size() ? 0 : -1;
    } else if (sig < 0) {
      auto it = lower_bound(p.begin() + 1, p.begin() + top, q);
      return ccw(it[-1], q, it[0]);
    } else {
      auto it = upper_bound(p.rbegin(), p.rend() - top - 1, q);
      auto pit_deref = it == p.rbegin() ? p[0] : it[-1];
      return ccw(*it, q, pit_deref);
    }
  }
  template <typename Function> int extreme(Function direction) const {
    int n = this->size(), left = 0, leftSig;
    const ConvexPolygon &poly = *this;
    auto vertex_cmp = [&poly, direction](int i, int j) {
      return ccw(poly[j] - poly[i], direction(poly[j]));
    };
    auto is_extreme = [n, vertex_cmp](int i, int &iSig) {
      return (iSig = vertex_cmp(i + 1, i)) >= 0 && vertex_cmp(i, i - 1) < 0;
    };
    for (int right = is_extreme(0, leftSig) ? 1 : n; left + 1 < right;) {
      int mid = (left + right) / 2, midSig;
      if (is_extreme(mid, midSig))
        return mid;
      if (leftSig != midSig ? leftSig < midSig
                            : leftSig == vertex_cmp(left, mid))
        right = mid;
      else
        left = mid, leftSig = midSig;
    }
    return poly.index(left);
  }
  void stab_extremes(const line &l, int &left, int &right) const {
    point direction = l.direction();
    right = extreme([&direction](const point &) { return direction; });
    left = extreme([&direction](const point &) { return -direction; });
  }
  friend vector<point> intersect(const ConvexPolygon &poly, const line &l) {
    point direction = l.direction();

    int left, right;
    poly.stab_extremes(l, left, right);
    auto vertex_cmp = [&l, &direction](const point &q) {
      return ccw(q - l.a, direction);
    };
    int rightSig = vertex_cmp(poly[right]), leftSig = vertex_cmp(poly[left]);
    if (rightSig < 0 || leftSig > 0)
      return {};
    auto intersectChain = [&l, &poly, vertex_cmp](int first, int last,
                                                  int firstSig) {
      int n = poly.size();
      while (poly.index(first + 1) != poly.index(last)) {
        int mid = (first + last + (first < last ? 0 : n)) / 2;
        mid = poly.index(mid);
        if (vertex_cmp(poly[mid]) == firstSig)
          first = mid;
        else
          last = mid;
      }
      return intersect(l, line(poly[first], poly[last]));
    };
    return {intersectChain(left, right, leftSig).first,
            intersectChain(right, left, rightSig).first};
  }
  friend bool has_intersection(const ConvexPolygon &p, const line &l) {
    point direction = l.direction();
    int left, right;
    p.stab_extremes(l, left, right);
    auto vertex_cmp = [&l, &direction](const point &q) {
      return ccw(q - l.a, direction);
    };
    int rightSig = vertex_cmp(p[right]), leftSig = vertex_cmp(p[left]);
    if (rightSig < 0 || leftSig > 0)
      return false;
    return true;
  }
  friend Large dist(const ConvexPolygon &p, const line &l) {
    point direction = l.direction();
    int left, right;
    p.stab_extremes(l, left, right);
    auto vertex_cmp = [&l, &direction](const point &q) {
      return ccw(q - l.a, direction);
    };
    int rightSig = vertex_cmp(p[right]), leftSig = vertex_cmp(p[left]);
    if (rightSig < 0 || leftSig > 0) {
      return min(dist(l, p[right]), dist(l, p[left]));
    } else {
      return 0;
    }
  }
  template <typename Function>
  friend void antipodals(const ConvexPolygon &poly, Function f) {
    if (poly.size() <= 1)
      return;
    if (poly.size() == 2)
      return void(f(0, 1));
    auto area = [&poly](int i, int j, int k) {
      return abs(cross(poly[i], poly[j], poly[k]));
    };
    auto func = [f, &poly](int i, int j) {
      return f(poly.index(i), poly.index(j));
    };

    int p = -1;
    int q = 0;
    while (area(p, p + 1, q + 1) > area(p, p + 1, q))
      q++;
    int p0 = 0;
    int q0 = q;
    while (poly.index(q) != p0) {
      p++;
      func(p, q);
      while (area(p, p + 1, q + 1) > area(p, p + 1, q)) {
        q++;
        if (poly.index(p) != poly.index(q0) || poly.index(q) != p0)
          func(p, q);
        else
          return;
      }
      if (area(p, p + 1, q + 1) == area(p, p + 1, q)) {
        if (poly.index(p) != poly.index(q0) || poly.index(q) != p0)
          func(p, q + 1);
        else
          func(p + 1, q);
      }
    }
  }
  friend ConvexPolygon minkowski_sum(const vector<ConvexPolygon> &v) {
    vector<point> vectors;
    point origin;
    for (auto &poly : v) {
      origin += poly[0];
      for (int i = 0; i < poly.size(); i++)
        vectors.push_back(poly[i + 1] - poly[i]);
    }
    angle_comparator::sortByAngle(vectors.begin(), vectors.end());
    auto last = point();
    if (!vectors.empty()) {
      last = vectors.back();
      vectors.pop_back();
    }
    vector<point> res;
    res.push_back(origin);
    for (auto &v : vectors) {
      res.push_back(res.back() + v);
      int n = res.size();
      if (n >= 3 && collinear(res[n - 3], res[n - 2], res[n - 1]))
        res.erase(res.begin() + n - 2);
    }
    int n = res.size();
    if (n >= 3 && collinear(res[n - 2], res[n - 1], res[0]))
      res.pop_back();
    if (res.size() >= 3 && collinear(res.back(), res[0], res[1]))
      res.erase(res.begin());
    return ConvexPolygon(res);
  }
  friend ConvexPolygon minkowski_sum(const ConvexPolygon &a,
                                     const ConvexPolygon &b) {
    vector<ConvexPolygon> v;
    v.push_back(a);
    v.push_back(b);
    return minkowski_sum(v);
  }
  friend ConvexPolygon intersect(const ConvexPolygon &a,
                                 const ConvexPolygon &b) {
    vector<point> candidates;
    auto consider = [&candidates](const ConvexPolygon &a,
                                  const ConvexPolygon &b) {
      for (int i = 0; i < a.size(); i++) {
        if (b.test(a[i]) <= 0)
          candidates.push_back(a[i]);
        segment s(a[i], a[i + 1]);
        vector<point> ps = intersect(b, s.as_line());
        for (auto p : ps) {
          if (s.contains(p))
            candidates.push_back(p);
        }
      }
    };
    consider(a, b);
    consider(b, a);
    auto res = ConvexPolygon(ConvexPolygon::convex_hull(candidates));
    return res;
  }
  friend Large intersection_area_or_dist(const ConvexPolygon &a,
                                         const ConvexPolygon &b) {
    ConvexPolygon inter = intersect(a, b);
    if (inter.size() > 0)
      return max(area(inter), Large(0));
    ConvexPolygon sum = minkowski_sum(a, -b);
    Large res = numeric_limits<Large>::max();
    for (int i = 0; i < sum.size(); i++) {
      res = min(res, dist(segment(sum[i], sum[i + 1]), point()));
    }
    return -res;
  }
  void cut(const halfplane& pl) {
    int n = this->size();
    if(n < 3) return;
    p.push_back(p[0]);

    auto pl_line = pl.as_line();

    vector<point> out;
    bool inside = pl.strictly_contains(p[0]);
    if(inside) out.push_back(p[0]);

    for(int i = 1; i <= n; i++) {
      if(pl.strictly_contains(p[i])) {
        if(!inside) {
          out.push_back(intersect(pl_line, line(p[i-1], p[i])).first);
        }
        out.push_back(p[i]);
        inside = true;
      } else {
        if(inside) {
          out.push_back(intersect(pl_line, line(p[i-1], p[i])).first);
        }
        inside = false;
      }
    }

    if(!out.empty() && out[0] == out.back()) out.pop_back();
    *this = ConvexPolygon(ConvexPolygon::convex_hull(out));
  }
  void cut(const ConvexPolygon &rhs) {
    for(int i = 0; i < rhs.size(); i++) {
      cut(halfplane::from_points(rhs[i], rhs[i+1]));
    }
  }
};
} // namespace plane

template <typename T, typename Large = T>
struct PolygonPlane : public CirclePlane<T, Large> {
  typedef plane::Polygon<T, Large> polygon;
  typedef plane::ConvexPolygon<T, Large> convex_polygon;
};

} // namespace geo
} // namespace lib


#line 6 "geometry/Caliper.cpp"

namespace lib {
using namespace std;
namespace geo {
namespace plane {
template <typename T, typename Large = T,
          typename enable_if<!is_integral<T>::value>::type * = nullptr,
          typename enable_if<!is_integral<T>::value>::type * = nullptr>
struct Caliper {
  typedef Point<T, Large> point;
  typedef Line<T, Large> line;
  point p;
  Large ang;
  Caliper(point a, Large alpha) : p(a) {
    ang = remainder(alpha, 2 * trig::PI);
    while (ang < 0)
      ang += 2 * trig::PI;
  }
  Large angle_to(const point &q) const {
    return remainder(arg(q - p) - ang, 2 * trig::PI);
  }
  void rotate(double theta) {
    ang += theta;
    while (ang > 2 * trig::PI)
      ang -= 2 * trig::PI;
    while (ang < 0)
      ang += 2 * trig::PI;
  }
  void move(const point &q) { p = q; }
  point versor() const { return point::polar(1.0, ang); }
  line as_line(Large scale = 1.0) const {
    return line(p, p + versor() * scale);
  }
  friend Large dist(const Caliper &a, const Caliper &b) {
    return dist(a.as_line(), b.p);
  }
};

template <typename T, typename Large = T> struct PolygonCalipers {
  constexpr static Large LIMIT = 4 * acosl(-1);

  typedef Point<T, Large> point;
  typedef Caliper<T, Large> caliper;
  typedef ConvexPolygon<T, Large> polygon;
  typedef pair<int, Large> descriptor;

  polygon poly;
  vector<caliper> calipers;
  vector<int> indices;
  vector<int> walked;
  Large angle_walked;

  PolygonCalipers(const polygon &poly, const vector<descriptor> &descriptors)
      : poly(poly), walked(descriptors.size()), angle_walked(0) {
    indices.reserve(descriptors.size());
    calipers.reserve(descriptors.size());
    for (size_t i = 0; i < descriptors.size(); i++) {
      calipers.emplace_back(poly[descriptors[i].first], descriptors[i].second);
      indices.emplace_back(descriptors[i].first);
    }
  }
  caliper operator[](int i) const { return calipers[i]; }
  int index(int i) const { return indices[i]; }
  bool has_next() const {
    return *min_element(walked.begin(), walked.end()) < poly.size() &&
           angle_walked < LIMIT;
  }
  Large angle_to_next(int i) const {
    int u = indices[i];
    return calipers[i].angle_to(poly[u + 1]);
  }
  void step_(int i) {
    int u = indices[i]++;
    indices[i] %= poly.size();
    calipers[i].move(poly[u + 1]);
    walked[i]++;
  }

  void next() {
    int i = 0;
    Large best = angle_to_next(0);
    for (size_t j = 1; j < calipers.size(); j++) {
      Large cur = angle_to_next(j);
      if (cur < best) {
        best = cur;
        i = j;
      }
    }
    Large alpha = angle_to_next(i);
    for (auto &caliper : calipers)
      caliper.rotate(alpha);
    step_(i);
    angle_walked += alpha;
  }
};
} // namespace plane

} // namespace geo
} // namespace lib

Back to top page