geometry/Circle2D.cpp

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

Depends on

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

Required by

• geometry/Caliper.cpp
• geometry/Polygon2D.cpp

Code

#ifndef _LIB_GEOMETRY_CIRCLE_2D
#define _LIB_GEOMETRY_CIRCLE_2D
#include "../utils/Annotation.cpp"
#include "Line2D.cpp"
#include <bits/stdc++.h>

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

#endif

#line 1 "geometry/Circle2D.cpp"


#line 1 "utils/Annotation.cpp"


#include <bits/stdc++.h>

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 1 "geometry/Line2D.cpp"


#line 1 "geometry/GeometryEpsilon.cpp"


#line 1 "Epsilon.cpp"


#line 4 "Epsilon.cpp"

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 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

Back to top page