From b8082f3600d8a18ed5a601ac519eeb6d7c79a5a6 Mon Sep 17 00:00:00 2001 From: baluteshih Date: Sat, 21 Oct 2023 01:36:22 +0800 Subject: [PATCH] Update Triangulation_Vonoroi --- .../8_Geometry/DelaunayTriangulation_dq.cpp | 168 ++++++++---------- codebook/8_Geometry/Triangulation_Vonoroi.cpp | 47 ++--- codebook/8_Geometry/point_in_circle.cpp | 17 +- codebook/content.tex | 6 +- 4 files changed, 93 insertions(+), 145 deletions(-) diff --git a/codebook/8_Geometry/DelaunayTriangulation_dq.cpp b/codebook/8_Geometry/DelaunayTriangulation_dq.cpp index fc2b3d86..5cac71ff 100644 --- a/codebook/8_Geometry/DelaunayTriangulation_dq.cpp +++ b/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::iterator twin; - Edge(int _id = 0):id(_id) {} + int id; // oidx[id] + list::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 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 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 getEdge() { - vector 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; diff --git a/codebook/8_Geometry/Triangulation_Vonoroi.cpp b/codebook/8_Geometry/Triangulation_Vonoroi.cpp index 0b0afebd..1a93193a 100644 --- a/codebook/8_Geometry/Triangulation_Vonoroi.cpp +++ b/codebook/8_Geometry/Triangulation_Vonoroi.cpp @@ -1,39 +1,12 @@ -vector 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 hpi = halfPlaneInter(ls[id]); - vector 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 mp; +// all coord. is even, you may want to call halfPlaneInter after then +vector> 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 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)); + } } diff --git a/codebook/8_Geometry/point_in_circle.cpp b/codebook/8_Geometry/point_in_circle.cpp index 15f954b3..c3edccc3 100644 --- a/codebook/8_Geometry/point_in_circle.cpp +++ b/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 &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 } diff --git a/codebook/content.tex b/codebook/content.tex index 5e71a103..757f5dfd 100644 --- a/codebook/content.tex +++ b/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}