Update Triangulation_Vonoroi

codebook/8_Geometry/DelaunayTriangulation_dq.cpp

@@ -1,102 +1,82 @@
-struct pddi {
-    pdd p;
-    int id;
-    pddi(double a = 0, double b = 0, int c = -1):p(pdd(a, b)), id(c) {}
-};
+/* Delaunay Triangulation:
+Given a sets of points on 2D plane, find a
+triangulation such that no points will strictly
+inside circumcircle of any triangle. */
 struct Edge {
-    int id;
-    list<Edge>::iterator twin;
-    Edge(int _id = 0):id(_id) {}
+  int id; // oidx[id]
+  list<Edge>::iterator twin;
+  Edge(int _id = 0):id(_id) {}
 };
-int inCircle(pddi a, pddi b, pddi c, pddi p) {
-    if (ori(a.p, b.p, c.p) < 0)
-        swap(b, c);
-    Point a3(a.p), b3(b.p), c3(c.p), p3(p.p);
-    b3 = b3 - a3, c3 = c3 - a3, p3 = p3 - a3;
-    Point f = cross(b3, c3); // normal vector
-    return sign(dot(p3, f)); // check same direction, in: < 0, on: = 0, out: > 0
-}
-int intersection(pddi a, pddi b, pddi c, pddi d) { // seg(a, b) and seg(c, d)
-    return ori(a.p, c.p, b.p) * ori(a.p, b.p, d.p) > 0 && ori(c.p, a.p, d.p) * ori(c.p, d.p, b.p) > 0;
-}
 struct Delaunay { // 0-base
-    list<Edge> head[N]; // graph
-    pddi p[N];
-    int n, rename[N];
-    void init(int _n, pddi _p[]) {
-        n = _n;
-        for (int i = 0; i < n; ++i) head[i].clear();
-        for (int i = 0; i < n; ++i) p[i] = _p[i];
-        sort(p, p + n, [&](pddi a, pddi b){return a.p < b.p;});
-        for (int i = 0; i < n; ++i) rename[p[i].id] = i;
-        divide(0, n - 1);
-    }
-    void addEdge(int u, int v) {
-        head[u].push_front(Edge(v));
-        head[v].push_front(Edge(u));
-        head[u].begin() -> twin = head[v].begin();
-        head[v].begin() -> twin = head[u].begin();
-    }
-    void divide(int l, int r) {
-        if (l == r) return;
-        if (l + 1 == r) return addEdge(l, l + 1);
-        int mid = (l + r) >> 1;
-        divide(l, mid), divide(mid + 1, r);
-        int nowl = l, nowr = r;
-        for (int update = 1; update;) {
-            update = 0;
-            pddi ptL = p[nowl], ptR = p[nowr];
-            for (auto it : head[nowl]) {
-                pddi t = p[it.id];
-                int v = ori(ptR.p, ptL.p, t.p);
-                if (v > 0 || (v == 0 && abs2(ptR.p - t.p) < abs2(ptR.p - ptL.p))) {
-                    nowl = it.id, update = 1;
-                    break;
-                }
-            }
-            if (update)	continue;
-            for (auto it : head[nowr]) {
-                pddi t = p[it.id];
-                int v = ori(ptL.p, ptR.p, t.p);
-                if (v < 0 || (v == 0 && abs2(ptL.p - t.p) < abs2(ptL.p - ptR.p))) {
-                    nowr = it.id, update = 1;
-                    break;
-                }
-            }
+  int n, oidx[N];
+  list<Edge> head[N]; // result udir. graph
+  pll p[N];
+  void init(int _n, pll _p[]) {
+    n = _n, iota(oidx, oidx + n, 0);
+    for (int i = 0; i < n; ++i) head[i].clear();
+    sort(oidx, oidx + n, [&](int a, int b) 
+    { return _p[a] < _p[b]; });
+    for (int i = 0; i < n; ++i) p[i] = _p[oidx[i]];
+    divide(0, n - 1);
+  }
+  void addEdge(int u, int v) {
+    head[u].push_front(Edge(v));
+    head[v].push_front(Edge(u));
+    head[u].begin()->twin = head[v].begin();
+    head[v].begin()->twin = head[u].begin();
+  }
+  void divide(int l, int r) {
+    if (l == r) return;
+    if (l + 1 == r) return addEdge(l, l + 1);
+    int mid = (l + r) >> 1;
+    divide(l, mid), divide(mid + 1, r);
+    int nowl = l, nowr = r;
+    for (int update = 1; update;) {
+      update = 0;
+      pll ptL = p[nowl], ptR = p[nowr];
+      for (auto it : head[nowl]) {
+        pll t = p[it.id];
+        int v = ori(ptR, ptL, t);
+        if (v > 0 || (v == 0 && abs2(ptR - t) < abs2(ptR - ptL))) {
+          nowl = it.id, update = 1;
+          break;
         }
-        addEdge(nowl, nowr); // add tangent
-        while (true) {
-            pddi ptL = p[nowl], ptR = p[nowr];
-            int ch = -1, side = 0;
-            for (auto it : head[nowl])
-                if (ori(ptL.p, ptR.p, p[it.id].p) > 0 && (ch == -1 || inCircle(ptL, ptR, p[ch], p[it.id]) < 0))
-                    ch = it.id, side = -1;
-            for (auto it : head[nowr])
-                if (ori(ptR.p, p[it.id].p, ptL.p) > 0 && (ch == -1 || inCircle(ptL, ptR, p[ch], p[it.id]) < 0))
-                    ch = it.id, side = 1;
-            if (ch == -1) break; // upper common tangent
-            if (side == -1) {
-                for (auto it = head[nowl].begin(); it != head[nowl].end(); )
-                    if (intersection(ptL, p[it -> id], ptR, p[ch]))
-                        head[it -> id].erase(it -> twin), head[nowl].erase(it++);
-                    else ++it;
-                nowl = ch, addEdge(nowl, nowr);
-            } 
-            else {
-                for (auto it = head[nowr].begin(); it != head[nowr].end(); )
-                    if (intersection(ptR, p[it -> id], ptL, p[ch]))
-                        head[it -> id].erase(it -> twin), head[nowr].erase(it++);
-                    else ++it;
-                nowr = ch, addEdge(nowl, nowr);
-            }
+      }
+      if (update)	continue;
+      for (auto it : head[nowr]) {
+        pll t = p[it.id];
+        int v = ori(ptL, ptR, t);
+        if (v < 0 || (v == 0 && abs2(ptL - t) < abs2(ptL - ptR))) {
+          nowr = it.id, update = 1;
+          break;
         }
+      }
     }
-    vector<pii> getEdge() {
-        vector<pii> ret;
-        for (int i = 0; i < n; ++i)
-            for (auto it : head[i])
-                if (it.id >= i)
-                    ret.pb(pii(p[i].id, p[it.id].id));
-        return ret;
+    addEdge(nowl, nowr); // add tangent
+    while (true) {
+      pll ptL = p[nowl], ptR = p[nowr];
+      int ch = -1, side = 0;
+      for (auto it : head[nowl])
+        if (ori(ptL, ptR, p[it.id]) > 0 && (ch == -1 || in_cc({ptL, ptR, p[ch]}, p[it.id])))
+          ch = it.id, side = -1;
+      for (auto it : head[nowr])
+        if (ori(ptR, p[it.id], ptL) > 0 && (ch == -1 || in_cc({ptL, ptR, p[ch]}, p[it.id])))
+          ch = it.id, side = 1;
+      if (ch == -1) break; // upper common tangent
+      if (side == -1) {
+        for (auto it = head[nowl].begin(); it != head[nowl].end(); )
+          if (seg_strict_intersect(ptL, p[it->id], ptR, p[ch]))
+            head[it->id].erase(it->twin), head[nowl].erase(it++);
+          else ++it;
+        nowl = ch, addEdge(nowl, nowr);
+      } 
+      else {
+        for (auto it = head[nowr].begin(); it != head[nowr].end(); )
+          if (seg_strict_intersect(ptR, p[it->id], ptL, p[ch]))
+            head[it->id].erase(it->twin), head[nowr].erase(it++);
+          else ++it;
+        nowr = ch, addEdge(nowl, nowr);
+      }
     }
+  }
 } tool;

codebook/8_Geometry/Triangulation_Vonoroi.cpp

@@ -1,39 +1,12 @@
-vector<Line> ls[N];
-pll arr[N];
-Line make_line(pdd p, Line l) {
-  pdd d = l.Y - l.X; d = perp(d);
-  pdd m = (l.X + l.Y) / 2;
-  l = Line(m, m + d);
-  if (ori(l.X, l.Y, p) < 0) 
-    l = Line(m + d, m);
-  return l;
-}
-double calc_area(int id) {
-  // use to calculate the area of point "strictly in the convex hull"
-  vector<Line> hpi = halfPlaneInter(ls[id]);
-  vector<pdd> ps;
-  for (int i = 0; i < SZ(hpi); ++i)
-    ps.pb(intersect(hpi[i].X, hpi[i].Y, hpi[(i + 1) % SZ(hpi)].X, hpi[(i + 1) % SZ(hpi)].Y));
-  double rt = 0;
-  for (int i = 0; i < SZ(ps); ++i)
-    rt += cross(ps[i], ps[(i + 1) % SZ(ps)]);
-  return fabs(rt) / 2;
-}
-void solve(int n, pii *oarr) {
-  map<pll, int> mp;
+// all coord. is even, you may want to call halfPlaneInter after then
+vector<vector<Line>> vec;
+void build_voronoi_line(int n, pll *arr) {
+  tool.init(n, arr); // Delaunay
+  vec.clear(), vec.resize(n);
   for (int i = 0; i < n; ++i)
-    arr[i] = pll(oarr[i].X, oarr[i].Y), mp[arr[i]] = i;
-  build(n, arr); // Triangulation
-  for (auto *t : triang) {
-    vector<int> p;
-    for (int i = 0; i < 3; ++i)
-      if (mp.find(t->p[i]) != mp.end())
-        p.pb(mp[t->p[i]]);
-    for (int i = 0; i < SZ(p); ++i)
-      for (int j = i + 1; j < SZ(p); ++j) { 
-        Line l(oarr[p[i]], oarr[p[j]]);
-        ls[p[i]].pb(make_line(oarr[p[i]], l));
-        ls[p[j]].pb(make_line(oarr[p[j]], l));
-      }
-  }
+    for (auto e : tool.head[i]) {
+      int u = tool.oidx[i], v = tool.oidx[e.id];
+      pll m = (arr[v] + arr[u]) / 2LL, d = perp(arr[v] - arr[u]);
+      vec[u].pb(Line(m, m + d));
+    }
 }

codebook/8_Geometry/point_in_circle.cpp

@@ -1,12 +1,7 @@
-// return p4 is strictly in circumcircle of tri(p1,p2,p3)
-ll sqr(ll x) { return x * x; }
-bool in_cc(const pll& p1, const pll& p2, const pll& p3, const pll& p4) {
-  ll u11 = p1.X - p4.X; ll u12 = p1.Y - p4.Y;
-  ll u21 = p2.X - p4.X; ll u22 = p2.Y - p4.Y;
-  ll u31 = p3.X - p4.X; ll u32 = p3.Y - p4.Y;
-  ll u13 = sqr(p1.X) - sqr(p4.X) + sqr(p1.Y) - sqr(p4.Y);
-  ll u23 = sqr(p2.X) - sqr(p4.X) + sqr(p2.Y) - sqr(p4.Y);
-  ll u33 = sqr(p3.X) - sqr(p4.X) + sqr(p3.Y) - sqr(p4.Y);
-  __int128 det = (__int128)-u13 * u22 * u31 + (__int128)u12 * u23 * u31 + (__int128)u13 * u21 * u32 - (__int128)u11 * u23 * u32 - (__int128)u12 * u21 * u33 + (__int128)u11 * u22 * u33;
-  return det > eps;
+// return q's relation with circumcircle of tri(p[0],p[1],p[2])
+bool in_cc(const array<pll, 3> &p, pll q) {
+  __int128 det = 0;
+  for (int i = 0; i < 3; ++i) 
+    det += __int128(abs2(p[i]) - abs2(q)) * cross(p[(i + 1) % 3] - q, p[(i + 2) % 3] - q);
+  return det > 0; // in: >0, on: =0, out: <0
 }

codebook/content.tex

@@ -235,9 +235,9 @@ \subsection{3Dpoint*} % test by HDU 3662
 \lstinputlisting{8_Geometry/3Dpoint.cpp}
 \subsection{Convexhull3D*} % test by HDU 3662, Ptz. Summer 2023 olmrgcsi And His Friends' Contest K
 \lstinputlisting{8_Geometry/Convexhull3D.cpp}
-\subsection{DelaunayTriangulation*} % test by Zerojudge b370
-\lstinputlisting{8_Geometry/DelaunayTriangulation.cpp}
-\subsection{Triangulation Vonoroi*} % test by NPSC 2018 territory
+\subsection{DelaunayTriangulation*} % test by Zerojudge b370, TIOJ 1310
+\lstinputlisting{8_Geometry/DelaunayTriangulation_dq.cpp}
+\subsection{Triangulation Vonoroi*} % test by CF 104555 J. Try NPSC 2018 territory?
 \lstinputlisting{8_Geometry/Triangulation_Vonoroi.cpp}
 \subsection{Minkowski Sum*} % test by Zerojudge b398
 \lstinputlisting{8_Geometry/Minkowski_Sum.cpp}