library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub maspypy/library

:warning: geo/delaunay_triangulation.hpp

Depends on

Code

#include "geo/base.hpp"

// とりあえずコピペしたものを使う

namespace CP_ALGO {
// https://cp-algorithms.com/geometry/delaunay.html
typedef long long ll;

bool ge(const ll& a, const ll& b) { return a >= b; }
bool le(const ll& a, const ll& b) { return a <= b; }
bool eq(const ll& a, const ll& b) { return a == b; }
bool gt(const ll& a, const ll& b) { return a > b; }
bool lt(const ll& a, const ll& b) { return a < b; }
int sgn(const ll& a) { return a >= 0 ? a ? 1 : 0 : -1; }

struct pt {
  ll x, y;
  pt() {}
  pt(ll _x, ll _y) : x(_x), y(_y) {}
  pt operator-(const pt& p) const { return pt(x - p.x, y - p.y); }
  ll cross(const pt& p) const { return x * p.y - y * p.x; }
  ll cross(const pt& a, const pt& b) const { return (a - *this).cross(b - *this); }
  ll dot(const pt& p) const { return x * p.x + y * p.y; }
  ll dot(const pt& a, const pt& b) const { return (a - *this).dot(b - *this); }
  ll sqrLength() const { return this->dot(*this); }
  bool operator==(const pt& p) const { return eq(x, p.x) && eq(y, p.y); }
};

const pt inf_pt = pt(1e18, 1e18);

struct QuadEdge {
  pt origin;
  QuadEdge* rot = nullptr;
  QuadEdge* onext = nullptr;
  bool used = false;
  QuadEdge* rev() const { return rot->rot; }
  QuadEdge* lnext() const { return rot->rev()->onext->rot; }
  QuadEdge* oprev() const { return rot->onext->rot; }
  pt dest() const { return rev()->origin; }
};

QuadEdge* make_edge(pt from, pt to) {
  QuadEdge* e1 = new QuadEdge;
  QuadEdge* e2 = new QuadEdge;
  QuadEdge* e3 = new QuadEdge;
  QuadEdge* e4 = new QuadEdge;
  e1->origin = from;
  e2->origin = to;
  e3->origin = e4->origin = inf_pt;
  e1->rot = e3;
  e2->rot = e4;
  e3->rot = e2;
  e4->rot = e1;
  e1->onext = e1;
  e2->onext = e2;
  e3->onext = e4;
  e4->onext = e3;
  return e1;
}

void splice(QuadEdge* a, QuadEdge* b) {
  swap(a->onext->rot->onext, b->onext->rot->onext);
  swap(a->onext, b->onext);
}

void delete_edge(QuadEdge* e) {
  splice(e, e->oprev());
  splice(e->rev(), e->rev()->oprev());
  delete e->rev()->rot;
  delete e->rev();
  delete e->rot;
  delete e;
}

QuadEdge* connect(QuadEdge* a, QuadEdge* b) {
  QuadEdge* e = make_edge(a->dest(), b->origin);
  splice(e, a->lnext());
  splice(e->rev(), b);
  return e;
}

bool left_of(pt p, QuadEdge* e) { return gt(p.cross(e->origin, e->dest()), 0); }

bool right_of(pt p, QuadEdge* e) { return lt(p.cross(e->origin, e->dest()), 0); }

template <class T>
T det3(T a1, T a2, T a3, T b1, T b2, T b3, T c1, T c2, T c3) {
  return a1 * (b2 * c3 - c2 * b3) - a2 * (b1 * c3 - c1 * b3) + a3 * (b1 * c2 - c1 * b2);
}

bool in_circle(pt a, pt b, pt c, pt d) {
// If there is __int128, calculate directly.
// Otherwise, calculate angles.
#if defined(__LP64__) || defined(_WIN64)
  __int128 det = -det3<__int128>(b.x, b.y, b.sqrLength(), c.x, c.y, c.sqrLength(), d.x, d.y, d.sqrLength());
  det += det3<__int128>(a.x, a.y, a.sqrLength(), c.x, c.y, c.sqrLength(), d.x, d.y, d.sqrLength());
  det -= det3<__int128>(a.x, a.y, a.sqrLength(), b.x, b.y, b.sqrLength(), d.x, d.y, d.sqrLength());
  det += det3<__int128>(a.x, a.y, a.sqrLength(), b.x, b.y, b.sqrLength(), c.x, c.y, c.sqrLength());
  return det > 0;
#else
  auto ang = [](pt l, pt mid, pt r) {
    ll x = mid.dot(l, r);
    ll y = mid.cross(l, r);
    long double res = atan2((long double)x, (long double)y);
    return res;
  };
  long double kek = ang(a, b, c) + ang(c, d, a) - ang(b, c, d) - ang(d, a, b);
  if (kek > 1e-8)
    return true;
  else
    return false;
#endif
}

pair<QuadEdge*, QuadEdge*> build_tr(int l, int r, vector<pt>& p) {
  if (r - l + 1 == 2) {
    QuadEdge* res = make_edge(p[l], p[r]);
    return make_pair(res, res->rev());
  }
  if (r - l + 1 == 3) {
    QuadEdge *a = make_edge(p[l], p[l + 1]), *b = make_edge(p[l + 1], p[r]);
    splice(a->rev(), b);
    int sg = sgn(p[l].cross(p[l + 1], p[r]));
    if (sg == 0) return make_pair(a, b->rev());
    QuadEdge* c = connect(b, a);
    if (sg == 1)
      return make_pair(a, b->rev());
    else
      return make_pair(c->rev(), c);
  }
  int mid = (l + r) / 2;
  QuadEdge *ldo, *ldi, *rdo, *rdi;
  tie(ldo, ldi) = build_tr(l, mid, p);
  tie(rdi, rdo) = build_tr(mid + 1, r, p);
  while (true) {
    if (left_of(rdi->origin, ldi)) {
      ldi = ldi->lnext();
      continue;
    }
    if (right_of(ldi->origin, rdi)) {
      rdi = rdi->rev()->onext;
      continue;
    }
    break;
  }
  QuadEdge* basel = connect(rdi->rev(), ldi);
  auto valid = [&basel](QuadEdge* e) { return right_of(e->dest(), basel); };
  if (ldi->origin == ldo->origin) ldo = basel->rev();
  if (rdi->origin == rdo->origin) rdo = basel;
  while (true) {
    QuadEdge* lcand = basel->rev()->onext;
    if (valid(lcand)) {
      while (in_circle(basel->dest(), basel->origin, lcand->dest(), lcand->onext->dest())) {
        QuadEdge* t = lcand->onext;
        delete_edge(lcand);
        lcand = t;
      }
    }
    QuadEdge* rcand = basel->oprev();
    if (valid(rcand)) {
      while (in_circle(basel->dest(), basel->origin, rcand->dest(), rcand->oprev()->dest())) {
        QuadEdge* t = rcand->oprev();
        delete_edge(rcand);
        rcand = t;
      }
    }
    if (!valid(lcand) && !valid(rcand)) break;
    if (!valid(lcand) || (valid(rcand) && in_circle(lcand->dest(), lcand->origin, rcand->origin, rcand->dest())))
      basel = connect(rcand, basel->rev());
    else
      basel = connect(basel->rev(), lcand->rev());
  }
  return make_pair(ldo, rdo);
}

vector<tuple<pt, pt, pt>> delaunay(vector<pt> p) {
  sort(p.begin(), p.end(), [](const pt& a, const pt& b) { return lt(a.x, b.x) || (eq(a.x, b.x) && lt(a.y, b.y)); });
  auto res = build_tr(0, (int)p.size() - 1, p);
  QuadEdge* e = res.first;
  vector<QuadEdge*> edges = {e};
  while (lt(e->onext->dest().cross(e->dest(), e->origin), 0)) e = e->onext;
  auto add = [&p, &e, &edges]() {
    QuadEdge* curr = e;
    do {
      curr->used = true;
      p.push_back(curr->origin);
      edges.push_back(curr->rev());
      curr = curr->lnext();
    } while (curr != e);
  };
  add();
  p.clear();
  int kek = 0;
  while (kek < (int)edges.size()) {
    if (!(e = edges[kek++])->used) add();
  }
  vector<tuple<pt, pt, pt>> ans;
  for (int i = 0; i < (int)p.size(); i += 3) { ans.push_back(make_tuple(p[i], p[i + 1], p[i + 2])); }
  return ans;
}
} // namespace CP_ALGO

template <typename T>
vc<tuple<int, int, int>> delaunay(vc<Point<T>> point) {
  using P = Point<T>;
  map<P, int> MP;
  vc<CP_ALGO::pt> tmp;
  int n = len(point);
  FOR(i, n) {
    P p = point[i];
    MP[p] = i;
    tmp.eb(CP_ALGO::pt(p.x, p.y));
  }
  auto tri = CP_ALGO::delaunay(tmp);
  vc<tuple<int, int, int>> ANS;
  for (auto [a, b, c]: tri) {
    int i = MP[P(a.x, a.y)];
    int j = MP[P(b.x, b.y)];
    int k = MP[P(c.x, c.y)];
    ANS.eb(i, j, k);
  }
  return ANS;
}
#line 1 "geo/delaunay_triangulation.hpp"

#line 2 "geo/base.hpp"
template <typename T>
struct Point {
  T x, y;

  Point() : x(0), y(0) {}

  template <typename A, typename B>
  Point(A x, B y) : x(x), y(y) {}

  template <typename A, typename B>
  Point(pair<A, B> p) : x(p.fi), y(p.se) {}

  Point operator+=(const Point p) {
    x += p.x, y += p.y;
    return *this;
  }
  Point operator-=(const Point p) {
    x -= p.x, y -= p.y;
    return *this;
  }
  Point operator+(Point p) const { return {x + p.x, y + p.y}; }
  Point operator-(Point p) const { return {x - p.x, y - p.y}; }
  bool operator==(Point p) const { return x == p.x && y == p.y; }
  bool operator!=(Point p) const { return x != p.x || y != p.y; }
  Point operator-() const { return {-x, -y}; }
  Point operator*(T t) const { return {x * t, y * t}; }
  Point operator/(T t) const { return {x / t, y / t}; }

  bool operator<(Point p) const {
    if (x != p.x) return x < p.x;
    return y < p.y;
  }
  T dot(const Point& other) const { return x * other.x + y * other.y; }
  T det(const Point& other) const { return x * other.y - y * other.x; }

  double norm() { return sqrtl(x * x + y * y); }
  double angle() { return atan2(y, x); }

  Point rotate(double theta) {
    static_assert(!is_integral<T>::value);
    double c = cos(theta), s = sin(theta);
    return Point{c * x - s * y, s * x + c * y};
  }
  Point rot90(bool ccw) { return (ccw ? Point{-y, x} : Point{y, -x}); }
};

#ifdef FASTIO
template <typename T>
void rd(Point<T>& p) {
  fastio::rd(p.x), fastio::rd(p.y);
}
template <typename T>
void wt(Point<T>& p) {
  fastio::wt(p.x);
  fastio::wt(' ');
  fastio::wt(p.y);
}
#endif

// A -> B -> C と進むときに、左に曲がるならば +1、右に曲がるならば -1
template <typename T>
int ccw(Point<T> A, Point<T> B, Point<T> C) {
  T x = (B - A).det(C - A);
  if (x > 0) return 1;
  if (x < 0) return -1;
  return 0;
}

template <typename REAL, typename T, typename U>
REAL dist(Point<T> A, Point<U> B) {
  REAL dx = REAL(A.x) - REAL(B.x);
  REAL dy = REAL(A.y) - REAL(B.y);
  return sqrt(dx * dx + dy * dy);
}

// ax+by+c
template <typename T>
struct Line {
  T a, b, c;

  Line(T a, T b, T c) : a(a), b(b), c(c) {}
  Line(Point<T> A, Point<T> B) { a = A.y - B.y, b = B.x - A.x, c = A.x * B.y - A.y * B.x; }
  Line(T x1, T y1, T x2, T y2) : Line(Point<T>(x1, y1), Point<T>(x2, y2)) {}

  template <typename U>
  U eval(Point<U> P) {
    return U(a) * P.x + U(b) * P.y + U(c);
  }

  template <typename U>
  T eval(U x, U y) {
    return a * x + b * y + c;
  }

  // 同じ直線が同じ a,b,c で表現されるようにする
  void normalize() {
    static_assert(is_same_v<T, int> || is_same_v<T, long long>);
    T g = gcd(gcd(abs(a), abs(b)), abs(c));
    a /= g, b /= g, c /= g;
    if (b < 0) { a = -a, b = -b, c = -c; }
    if (b == 0 && a < 0) { a = -a, b = -b, c = -c; }
  }

  bool is_parallel(Line other) { return a * other.b - b * other.a == 0; }
  bool is_orthogonal(Line other) { return a * other.a + b * other.b == 0; }
};

template <typename T>
struct Segment {
  Point<T> A, B;

  Segment(Point<T> A, Point<T> B) : A(A), B(B) {}
  Segment(T x1, T y1, T x2, T y2) : Segment(Point<T>(x1, y1), Point<T>(x2, y2)) {}

  bool contain(Point<T> C) {
    T det = (C - A).det(B - A);
    if (det != 0) return 0;
    return (C - A).dot(B - A) >= 0 && (C - B).dot(A - B) >= 0;
  }

  Line<T> to_Line() { return Line(A, B); }
};

template <typename REAL>
struct Circle {
  Point<REAL> O;
  REAL r;
  Circle() {}
  Circle(Point<REAL> O, REAL r) : O(O), r(r) {}
  Circle(REAL x, REAL y, REAL r) : O(x, y), r(r) {}
  template <typename T>
  bool contain(Point<T> p) {
    REAL dx = p.x - O.x, dy = p.y - O.y;
    return dx * dx + dy * dy <= r * r;
  }
};
#line 3 "geo/delaunay_triangulation.hpp"

// とりあえずコピペしたものを使う

namespace CP_ALGO {
// https://cp-algorithms.com/geometry/delaunay.html
typedef long long ll;

bool ge(const ll& a, const ll& b) { return a >= b; }
bool le(const ll& a, const ll& b) { return a <= b; }
bool eq(const ll& a, const ll& b) { return a == b; }
bool gt(const ll& a, const ll& b) { return a > b; }
bool lt(const ll& a, const ll& b) { return a < b; }
int sgn(const ll& a) { return a >= 0 ? a ? 1 : 0 : -1; }

struct pt {
  ll x, y;
  pt() {}
  pt(ll _x, ll _y) : x(_x), y(_y) {}
  pt operator-(const pt& p) const { return pt(x - p.x, y - p.y); }
  ll cross(const pt& p) const { return x * p.y - y * p.x; }
  ll cross(const pt& a, const pt& b) const { return (a - *this).cross(b - *this); }
  ll dot(const pt& p) const { return x * p.x + y * p.y; }
  ll dot(const pt& a, const pt& b) const { return (a - *this).dot(b - *this); }
  ll sqrLength() const { return this->dot(*this); }
  bool operator==(const pt& p) const { return eq(x, p.x) && eq(y, p.y); }
};

const pt inf_pt = pt(1e18, 1e18);

struct QuadEdge {
  pt origin;
  QuadEdge* rot = nullptr;
  QuadEdge* onext = nullptr;
  bool used = false;
  QuadEdge* rev() const { return rot->rot; }
  QuadEdge* lnext() const { return rot->rev()->onext->rot; }
  QuadEdge* oprev() const { return rot->onext->rot; }
  pt dest() const { return rev()->origin; }
};

QuadEdge* make_edge(pt from, pt to) {
  QuadEdge* e1 = new QuadEdge;
  QuadEdge* e2 = new QuadEdge;
  QuadEdge* e3 = new QuadEdge;
  QuadEdge* e4 = new QuadEdge;
  e1->origin = from;
  e2->origin = to;
  e3->origin = e4->origin = inf_pt;
  e1->rot = e3;
  e2->rot = e4;
  e3->rot = e2;
  e4->rot = e1;
  e1->onext = e1;
  e2->onext = e2;
  e3->onext = e4;
  e4->onext = e3;
  return e1;
}

void splice(QuadEdge* a, QuadEdge* b) {
  swap(a->onext->rot->onext, b->onext->rot->onext);
  swap(a->onext, b->onext);
}

void delete_edge(QuadEdge* e) {
  splice(e, e->oprev());
  splice(e->rev(), e->rev()->oprev());
  delete e->rev()->rot;
  delete e->rev();
  delete e->rot;
  delete e;
}

QuadEdge* connect(QuadEdge* a, QuadEdge* b) {
  QuadEdge* e = make_edge(a->dest(), b->origin);
  splice(e, a->lnext());
  splice(e->rev(), b);
  return e;
}

bool left_of(pt p, QuadEdge* e) { return gt(p.cross(e->origin, e->dest()), 0); }

bool right_of(pt p, QuadEdge* e) { return lt(p.cross(e->origin, e->dest()), 0); }

template <class T>
T det3(T a1, T a2, T a3, T b1, T b2, T b3, T c1, T c2, T c3) {
  return a1 * (b2 * c3 - c2 * b3) - a2 * (b1 * c3 - c1 * b3) + a3 * (b1 * c2 - c1 * b2);
}

bool in_circle(pt a, pt b, pt c, pt d) {
// If there is __int128, calculate directly.
// Otherwise, calculate angles.
#if defined(__LP64__) || defined(_WIN64)
  __int128 det = -det3<__int128>(b.x, b.y, b.sqrLength(), c.x, c.y, c.sqrLength(), d.x, d.y, d.sqrLength());
  det += det3<__int128>(a.x, a.y, a.sqrLength(), c.x, c.y, c.sqrLength(), d.x, d.y, d.sqrLength());
  det -= det3<__int128>(a.x, a.y, a.sqrLength(), b.x, b.y, b.sqrLength(), d.x, d.y, d.sqrLength());
  det += det3<__int128>(a.x, a.y, a.sqrLength(), b.x, b.y, b.sqrLength(), c.x, c.y, c.sqrLength());
  return det > 0;
#else
  auto ang = [](pt l, pt mid, pt r) {
    ll x = mid.dot(l, r);
    ll y = mid.cross(l, r);
    long double res = atan2((long double)x, (long double)y);
    return res;
  };
  long double kek = ang(a, b, c) + ang(c, d, a) - ang(b, c, d) - ang(d, a, b);
  if (kek > 1e-8)
    return true;
  else
    return false;
#endif
}

pair<QuadEdge*, QuadEdge*> build_tr(int l, int r, vector<pt>& p) {
  if (r - l + 1 == 2) {
    QuadEdge* res = make_edge(p[l], p[r]);
    return make_pair(res, res->rev());
  }
  if (r - l + 1 == 3) {
    QuadEdge *a = make_edge(p[l], p[l + 1]), *b = make_edge(p[l + 1], p[r]);
    splice(a->rev(), b);
    int sg = sgn(p[l].cross(p[l + 1], p[r]));
    if (sg == 0) return make_pair(a, b->rev());
    QuadEdge* c = connect(b, a);
    if (sg == 1)
      return make_pair(a, b->rev());
    else
      return make_pair(c->rev(), c);
  }
  int mid = (l + r) / 2;
  QuadEdge *ldo, *ldi, *rdo, *rdi;
  tie(ldo, ldi) = build_tr(l, mid, p);
  tie(rdi, rdo) = build_tr(mid + 1, r, p);
  while (true) {
    if (left_of(rdi->origin, ldi)) {
      ldi = ldi->lnext();
      continue;
    }
    if (right_of(ldi->origin, rdi)) {
      rdi = rdi->rev()->onext;
      continue;
    }
    break;
  }
  QuadEdge* basel = connect(rdi->rev(), ldi);
  auto valid = [&basel](QuadEdge* e) { return right_of(e->dest(), basel); };
  if (ldi->origin == ldo->origin) ldo = basel->rev();
  if (rdi->origin == rdo->origin) rdo = basel;
  while (true) {
    QuadEdge* lcand = basel->rev()->onext;
    if (valid(lcand)) {
      while (in_circle(basel->dest(), basel->origin, lcand->dest(), lcand->onext->dest())) {
        QuadEdge* t = lcand->onext;
        delete_edge(lcand);
        lcand = t;
      }
    }
    QuadEdge* rcand = basel->oprev();
    if (valid(rcand)) {
      while (in_circle(basel->dest(), basel->origin, rcand->dest(), rcand->oprev()->dest())) {
        QuadEdge* t = rcand->oprev();
        delete_edge(rcand);
        rcand = t;
      }
    }
    if (!valid(lcand) && !valid(rcand)) break;
    if (!valid(lcand) || (valid(rcand) && in_circle(lcand->dest(), lcand->origin, rcand->origin, rcand->dest())))
      basel = connect(rcand, basel->rev());
    else
      basel = connect(basel->rev(), lcand->rev());
  }
  return make_pair(ldo, rdo);
}

vector<tuple<pt, pt, pt>> delaunay(vector<pt> p) {
  sort(p.begin(), p.end(), [](const pt& a, const pt& b) { return lt(a.x, b.x) || (eq(a.x, b.x) && lt(a.y, b.y)); });
  auto res = build_tr(0, (int)p.size() - 1, p);
  QuadEdge* e = res.first;
  vector<QuadEdge*> edges = {e};
  while (lt(e->onext->dest().cross(e->dest(), e->origin), 0)) e = e->onext;
  auto add = [&p, &e, &edges]() {
    QuadEdge* curr = e;
    do {
      curr->used = true;
      p.push_back(curr->origin);
      edges.push_back(curr->rev());
      curr = curr->lnext();
    } while (curr != e);
  };
  add();
  p.clear();
  int kek = 0;
  while (kek < (int)edges.size()) {
    if (!(e = edges[kek++])->used) add();
  }
  vector<tuple<pt, pt, pt>> ans;
  for (int i = 0; i < (int)p.size(); i += 3) { ans.push_back(make_tuple(p[i], p[i + 1], p[i + 2])); }
  return ans;
}
} // namespace CP_ALGO

template <typename T>
vc<tuple<int, int, int>> delaunay(vc<Point<T>> point) {
  using P = Point<T>;
  map<P, int> MP;
  vc<CP_ALGO::pt> tmp;
  int n = len(point);
  FOR(i, n) {
    P p = point[i];
    MP[p] = i;
    tmp.eb(CP_ALGO::pt(p.x, p.y));
  }
  auto tri = CP_ALGO::delaunay(tmp);
  vc<tuple<int, int, int>> ANS;
  for (auto [a, b, c]: tri) {
    int i = MP[P(a.x, a.y)];
    int j = MP[P(b.x, b.y)];
    int k = MP[P(c.x, c.y)];
    ANS.eb(i, j, k);
  }
  return ANS;
}
Back to top page