Add Delaunay triangulation and graph templates

content/geometry/algorithms/convex-hull-3d.cpp

@@ -1,112 +1,101 @@
 /*
- *Description:* 3D Convex Hull using randomized incremental algorithm
+ *Description:* 3D Convex Hull using randomized incremental algorithm. Exact integer arithmetic. Sets `degenerate=true` if less than 4 non-coplanar points.
  *Time:* $O(n log n)$ expected
  *Status:* Tested
  */
 template <typename T> struct convex_hull_3d {
   using pt = point3d<T>;
-  const ld EPS = 1e-9;
   struct edge;
   struct face {
     int a, b, c;
-    pt pa, pb, pc, q;
+    pt q;
     edge *e1, *e2, *e3;
     vector<int> points;
     int dead = 1e9;
     face(int a, int b, int c, pt pa, pt pb, pt pc)
-        : a(a), b(b), c(c), pa(pa), pb(pb), pc(pc) {
-      q = ((pb - pa)^(pc - pa));
+        : a(a), b(b), c(c) {
+      q = (pb - pa) ^ (pc - pa);
       e1 = e2 = e3 = nullptr;
     }
   };
-
-  struct edge {
-    edge *rev; face *f;
-  };
+  struct edge { edge *rev; face *f; };
 
   vector<pt> p;
   vector<face *> f;
+  bool degenerate = false;
 
   static void glue(face *F1, face *F2, edge *&e1, edge *&e2) {
-    e1 = new edge, e2 = new edge;
-    e1->rev = e2, e2->rev = e1;
-    e1->f = F2, e2->f = F1;
+    e1 = new edge; e2 = new edge;
+    e1->rev = e2; e2->rev = e1;
+    e1->f = F2; e2->f = F1;
   }
 
   void prepare() {
     int n = (int)p.size();
-    mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
+    mt19937 rng(42);
     shuffle(p.begin(), p.end(), rng);
     vector<int> ve = {0};
-    for (int i = 1; i < n; i++) {
+    for (int i = 1; i < n && (int)ve.size() < 4; i++) {
       if (ve.size() == 1) {
-        if ((p[ve[0]] - p[i]).abs() > EPS) ve.push_back(i);
+        if ((p[ve[0]] - p[i]).nonzero()) ve.push_back(i);
       } else if (ve.size() == 2) {
-        if ((p[ve[1]] - p[ve[0]]).cross(p[i] - p[ve[0]]).abs() > EPS)
+        if (((p[ve[1]] - p[ve[0]]) ^ (p[i] - p[ve[0]])).nonzero())
+          ve.push_back(i);
+      } else {
+        if ((p[i] - p[ve[0]]).dot((p[ve[1]] - p[ve[0]]) ^ (p[ve[2]] - p[ve[0]])) != 0)
           ve.push_back(i);
-      } else if (abs((p[i] - p[ve[0]])
-                .dot((p[ve[1]] - p[ve[0]])
-                .cross(p[ve[2]] - p[ve[0]]))) > EPS) {
-        ve.push_back(i);
-        break;
       }
     }
-    assert((int)ve.size() == 4);
+    if ((int)ve.size() < 4) { degenerate = true; return; }
     vector<pt> ve2;
     for (int i : ve) ve2.push_back(p[i]);
     reverse(ve.begin(), ve.end());
     for (int i : ve) p.erase(p.begin() + i);
     p.insert(p.begin(), ve2.begin(), ve2.end());
   }
 
-  convex_hull_3d(vector<pt> points) : p(move(points)) {
-    int n = p.size();
+  convex_hull_3d() {}
+  convex_hull_3d(vector<pt> points) : p(move(points)) { build_hull(); }
+
+  void build_hull() {
+    int n = (int)p.size();
     prepare();
+    if (degenerate) return;
     vector<face *> new_face(n, nullptr);
     vector<vector<face *>> conflict(n);
-
     auto add_face = [&](int a, int b, int c) {
       face *F = new face(a, b, c, p[a], p[b], p[c]);
-      f.push_back(F);
-      return F;
+      f.push_back(F); return F;
     };
-
     face *F1 = add_face(0, 1, 2);
     face *F2 = add_face(0, 2, 1);
     glue(F1, F2, F1->e1, F2->e3);
     glue(F1, F2, F1->e2, F2->e2);
     glue(F1, F2, F1->e3, F2->e1);
-
     for (int i = 3; i < n; i++) {
       for (face *F : {F1, F2}) {
         T Q = (p[i] - p[F->a]).dot(F->q);
-        if (Q > EPS) conflict[i].push_back(F);
-        if (Q >= -EPS) F->points.push_back(i);
+        if (Q > 0) conflict[i].push_back(F);
+        if (Q >= 0) F->points.push_back(i);
       }
     }
-
     for (int i = 3; i < n; i++) {
-      for (face *F : conflict[i])
-        F->dead = min(F->dead, i);
+      for (face *F : conflict[i]) F->dead = min(F->dead, i);
       int v = -1;
       for (face *F : conflict[i]) {
-        if (F->dead != i)
-          continue;
+        if (F->dead != i) continue;
         int parr[3] = {F->a, F->b, F->c};
         edge *earr[3] = {F->e1, F->e2, F->e3};
         for (int j = 0; j < 3; j++) {
-          int j2 = (j + 1)%3;
+          int j2 = (j + 1) % 3;
           if (earr[j]->f->dead > i) {
             face *Fn = new_face[parr[j]] = add_face(parr[j], parr[j2], i);
             set_union(F->points.begin(), F->points.end(),
                       earr[j]->f->points.begin(), earr[j]->f->points.end(),
                       back_inserter(Fn->points));
             Fn->points.erase(
                 stable_partition(Fn->points.begin(), Fn->points.end(),
-                                 [&](int k) {
-                                   return k > i and
-                                          (p[k] - p[Fn->a]).dot(Fn->q) > EPS;
-                                 }),
+                                 [&](int k) { return k > i && (p[k] - p[Fn->a]).dot(Fn->q) > 0; }),
                 Fn->points.end());
             for (int k : Fn->points) conflict[k].push_back(Fn);
             earr[j]->rev->f = Fn;
@@ -122,26 +111,24 @@ template <typename T> struct convex_hull_3d {
         v = u;
       }
     }
-    f.erase(remove_if(f.begin(), f.end(), [&](face *F) { return F->dead < n; }), f.end());
+    f.erase(remove_if(f.begin(), f.end(),
+                       [&](face *F) { return F->dead < n; }), f.end());
   }
 
-  int num_faces() const { return f.size(); }
-  T area() {
-    T res = 0;
-    for (face *F : f)
-      res += F->q.abs() / 2.0;
+  int num_faces() const { return (int)f.size(); }
+  double area() const {
+    double res = 0;
+    for (face *F : f) res += F->q.abs() / 2.0;
     return res;
   }
-  T volume() {
-    T vol = 0;
-    for (face *F : f)
-      vol += F->pa.dot(F->q);
-    return abs(vol) / 6.0;
+  double volume() const {
+    double vol = 0;
+    for (face *F : f) vol += p[F->a].dot(F->q);
+    return std::abs(vol) / 6.0;
   }
-  vector<array<int, 3>> faces() {
+  vector<array<int, 3>> faces() const {
     vector<array<int, 3>> res;
-    for (face *F : f)
-      res.push_back({F->a, F->b, F->c});
+    for (face *F : f) res.push_back({F->a, F->b, F->c});
     return res;
   }
 };

content/geometry/algorithms/delaunay-graph.cpp

@@ -0,0 +1,41 @@
+/*
+ *Description:* Planar graph of a Delaunay Triangulation with greedy walk and locate. `walk` finds nearest site. `locate` finds all equidistant nearest sites via BFS: 1 site $=$ Voronoi region, 2 $=$ edge, 3+ $=$ vertex (cocircular).
+ *Time:* Build $O(n log n)$ | Walk $O(sqrt(n))$ expected | Locate $O(sqrt(n))$ expected
+ *Status:* Tested
+ *Dependencies:* point2d, delaunay
+ */
+template <typename T> struct delaunay_graph : delaunay<T> {
+  vector<vector<int>> adj;
+  delaunay_graph() {}
+  delaunay_graph(vector<T> &px, vector<T> &py) { build(px, py); }
+  void build(vector<T> &px, vector<T> &py) {
+    delaunay<T>::build(px, py);
+    adj.assign(this->n, {});
+    for (auto &t : this->tri)
+      for (int i = 0; i < 3; i++) for (int j = i+1; j < 3; j++) {
+        adj[t[i]].push_back(t[j]); adj[t[j]].push_back(t[i]);
+      }
+    for (auto &v : adj) { sort(v.begin(), v.end()); v.erase(unique(v.begin(), v.end()), v.end()); }
+  }
+  int walk(vector<T> &px, vector<T> &py, T qx, T qy) const {
+    auto d2 = [&](int i) { return (qx-px[i])*(qx-px[i]) + (qy-py[i])*(qy-py[i]); };
+    int cur = 0; T best = d2(0);
+    for (bool ok = true; ok; ) { ok = false;
+      for (int nb : adj[cur]) { T d = d2(nb); if (d < best) { best=d; cur=nb; ok=true; break; } }
+    }
+    return cur;
+  }
+  vector<int> locate(vector<T> &px, vector<T> &py, T qx, T qy) const {
+    auto d2 = [&](int i) { return (qx-px[i])*(qx-px[i]) + (qy-py[i])*(qy-py[i]); };
+    int v = walk(px, py, qx, qy);
+    T best = d2(v);
+    vector<int> res = {v};
+    vector<bool> vis(this->n, false);
+    vis[v] = true;
+    for (int qi = 0; qi < (int)res.size(); qi++)
+      for (int nb : adj[res[qi]])
+        if (!vis[nb]) { vis[nb] = true; if (d2(nb) == best) res.push_back(nb); }
+    sort(res.begin(), res.end());
+    return res;
+  }
+};

content/geometry/algorithms/delaunay.cpp

@@ -0,0 +1,105 @@
+/*
+ *Description:* Delaunay Triangulation using divide and conquer with quad-edges. Each circumcircle contains none of the input points. If all points are collinear, `tri` is empty. Uses `__int128` for circumcircle test (coords up to $10^9$).
+ *Time:* $O(n log n)$
+ *Status:* Tested
+ *Dependencies:* point2d
+ */
+template <typename T> struct delaunay {
+  using pt = point2d<T>;
+  using lll = __int128_t;
+  struct Quad {
+    Quad *rot, *o;
+    pt p; int id; bool mark;
+    pt& F() { return rot->rot->p; }
+    int& Fi() { return rot->rot->id; }
+    Quad*& r() { return rot->rot; }
+    Quad* prev() { return rot->o->rot; }
+    Quad* next() { return r()->prev(); }
+  };
+  Quad* pool;
+  bool circ(pt p, pt a, pt b, pt c) {
+    lll p2=p.norm(), A=a.norm()-p2, B=b.norm()-p2, C=c.norm()-p2;
+    return p.cross(a,b)*C + p.cross(b,c)*A + p.cross(c,a)*B > 0;
+  }
+  Quad* makeEdge(pt orig, int oi, pt dest, int di) {
+    auto nq = [&]() { return new Quad{nullptr,nullptr,{}, -1, false}; };
+    Quad* r;
+    if (pool) { r = pool; pool = r->o; }
+    else { r = nq(); r->rot = nq(); r->rot->rot = nq(); r->rot->rot->rot = nq(); }
+    r->r()->r() = r; r->rot->rot->rot->rot = r;
+    pt arb{(T)1e18, (T)1e18};
+    Quad* c = r;
+    for (int i = 0; i < 4; i++) {
+      c = c->rot; c->p = arb; c->id = -1;
+      c->o = (i & 1) ? c : c->r();
+    }
+    r->p = orig; r->id = oi; r->F() = dest; r->Fi() = di;
+    return r;
+  }
+  void splice(Quad* a, Quad* b) {
+    swap(a->o->rot->o, b->o->rot->o); swap(a->o, b->o);
+  }
+  Quad* connect(Quad* a, Quad* b) {
+    Quad* q = makeEdge(a->F(), a->Fi(), b->p, b->id);
+    splice(q, a->next()); splice(q->r(), b); return q;
+  }
+  pair<Quad*,Quad*> rec(vector<pt>& s, vector<int>& ids, int lo, int hi) {
+    int sz = hi - lo;
+    if (sz <= 3) {
+      Quad* a = makeEdge(s[lo],ids[lo], s[lo+1],ids[lo+1]);
+      Quad* b = makeEdge(s[lo+1],ids[lo+1], s[lo+sz-1],ids[lo+sz-1]);
+      if (sz == 2) return {a, a->r()};
+      splice(a->r(), b);
+      auto side = s[lo].cross(s[lo+1], s[lo+2]);
+      Quad* c = side ? connect(b, a) : nullptr;
+      return {side<0 ? c->r() : a, side<0 ? c : b->r()};
+    }
+    int mid = lo + sz/2;
+    auto [ra, A] = rec(s, ids, lo, mid);
+    auto [B, rb] = rec(s, ids, mid, hi);
+    while ((B->p.cross(A->F(), A->p)<0 && (A=A->next())) ||
+           (A->p.cross(B->F(), B->p)>0 && (B=B->r()->o)));
+    Quad* base = connect(B->r(), A);
+    if (A->p.x==ra->p.x && A->p.y==ra->p.y) ra = base->r();
+    if (B->p.x==rb->p.x && B->p.y==rb->p.y) rb = base;
+    auto valid = [&](Quad* e) { return e->F().cross(base->F(), base->p) > 0; };
+    for (;;) {
+      Quad* LC = base->r()->o;
+      if (valid(LC)) while (circ(LC->o->F(), base->F(), base->p, LC->F()))
+        { Quad* t=LC->o; splice(LC,LC->prev()); splice(LC->r(),LC->r()->prev()); LC->o=pool; pool=LC; LC=t; }
+      Quad* RC = base->prev();
+      if (valid(RC)) while (circ(RC->prev()->F(), base->F(), base->p, RC->F()))
+        { Quad* t=RC->prev(); splice(RC,RC->prev()); splice(RC->r(),RC->r()->prev()); RC->o=pool; pool=RC; RC=t; }
+      if (!valid(LC) && !valid(RC)) break;
+      if (!valid(LC) || (valid(RC) && circ(RC->F(), base->F(), base->p, LC->F())))
+        base = connect(RC, base->r());
+      else base = connect(base->r(), LC->r());
+    }
+    return {ra, rb};
+  }
+  int n;
+  vector<array<int, 3>> tri;
+  delaunay() : pool(nullptr), n(0) {}
+  delaunay(vector<T> &px, vector<T> &py) : pool(nullptr), n(0) { build(px, py); }
+  void build(vector<T> &px, vector<T> &py) {
+    n = (int)px.size(); tri.clear(); pool = nullptr;
+    if (n < 2) return;
+    vector<int> ord(n); iota(ord.begin(), ord.end(), 0);
+    sort(ord.begin(), ord.end(), [&](int a, int b) {
+      return px[a]<px[b] || (px[a]==px[b] && py[a]<py[b]);
+    });
+    vector<pt> sp(n); vector<int> si(n);
+    for (int i = 0; i < n; i++) { sp[i] = {px[ord[i]], py[ord[i]]}; si[i] = ord[i]; }
+    Quad* e = rec(sp, si, 0, n).first;
+    vector<Quad*> q = {e};
+    int qi = 0;
+    while (e->o->F().cross(e->F(), e->p) < 0) e = e->o;
+    auto add = [&](Quad* e) {
+      Quad* c = e; vector<int> face;
+      do { c->mark=1; face.push_back(c->id); q.push_back(c->r()); c=c->next(); } while (c!=e);
+      if ((int)face.size()==3) tri.push_back({face[0],face[1],face[2]});
+    };
+    add(e);
+    while (qi < (int)q.size()) if (!(e=q[qi++])->mark) add(e);
+  }
+};

content/geometry/primitives/point-3d.cpp

@@ -51,6 +51,7 @@ template <typename T> struct point3d {
     return dot(a ^ b);
   }
 
+  bool nonzero() const { return x != 0 || y != 0 || z != 0; }
   T norm() const { return dot(*this); }
   T sq_dist(const point3d<T> &b) const { return (*this - b).norm(); }
   double abs() const { return sqrt((double)norm()); }