sotanishy's code snippets for competitive programming
View the Project on GitHub sotanishy/cp-library-cpp
#include "geometry/delaunay_diagram.hpp"
与えられた点の Delaunay 図は,Voronoi 図の双対であり,Voronoi 領域の隣接関係を表す.
Voronoi 図とは,平面に配置された母点に対して,平面の各点がどの母点に最も近いかによって領域分けされた図である.
Delaunay 辺の本数は Euler の多面体定理より $O(n)$ 本である.また,Euclidean MST は Delaunay 図の部分グラフである.
Delaunay 図は平面走査アルゴリズムである Fortune のアルゴリズムで求めることができる.
vector<pair<int, int>> delaunay_diagram(vector<Vec> pts)
pts
#pragma once #include <queue> #include <set> #include <utility> #include <vector> #include "geometry.hpp" #include "triangle.hpp" std::vector<std::pair<int, int>> delaunay_diagram(std::vector<Vec> pts) { const int n = pts.size(); constexpr T INF = 1e20; // stores delaunay edges std::vector<std::pair<int, int>> edges; // the sweep line moves from left to right // x coordinate of the sweepline static T sweepx; // check if all points are collinear bool ok = false; for (int i = 2; i < n; ++i) { if (ccw(pts[0], pts[1], pts[i]) != 0) { ok = true; break; } } // handle degenerate cases if (!ok) { // connect all points by a path std::vector<std::pair<Vec, int>> pi(n); for (int i = 0; i < n; ++i) pi[i] = {pts[i], i}; std::ranges::sort(pi, {}, [&](auto& p) { return std::make_pair(p.first.real(), p.first.imag()); }); for (int k = 0; k < n - 1; ++k) { edges.emplace_back(pi[k].second, pi[k + 1].second); } return edges; } // represents a parabola given by the focus p // and the focus q of the next parabola struct Parabola { Vec p; mutable Vec q; int i; mutable int id = 0; Parabola(const Vec& p, const Vec& q, int i) : p(p), q(q), i(i) {} T gety(T x) const { if (q.imag() == INF) return INF; x += eps; Vec m = (p + q) / T(2); Vec dir = rot(p - m, PI / 2); T D = (x - p.real()) * (x - q.real()); return m.imag() + ((m.real() - x) * dir.real() + std::sqrt(D) * std::abs(dir)) / dir.imag(); } bool operator<(T y) const { return gety(sweepx) < y; } bool operator<(const Parabola& o) const { return gety(sweepx) < o.gety(sweepx); } }; // maintains a list of parabola std::multiset<Parabola, std::less<>> beach; using iterator = decltype(beach)::iterator; // represents an event // id >= 0: point event // id < 0: vertex event struct Event { T x; int id; iterator it; Event(T x, int id, iterator it) : x(x), id(id), it(it) {} bool operator>(const Event& e) const { return x > e.x; } }; // maintains a list of events std::priority_queue<Event, std::vector<Event>, std::greater<>> pq; // true if the i-th vertex event is still valid std::vector<bool> valid(1); int k = 1; // next id // create vertex events auto update = [&](const iterator& it) { if (it->i <= -1 || it == beach.begin()) return; // sentinels valid[it->id] = false; auto a = std::prev(it); // handle a vertex event that occurs // when 3 parabolas at (a->p, it->p, it->q) coincide if (eq(ccw(a->p, it->p, it->q), 0)) return; // never coincide when collinear // create a new event it->id = k++; valid.push_back(true); Vec c = circumcenter(a->p, it->p, it->q); // coordinate of the vertex event T x = c.real() + std::abs(c - it->p); if (leq(sweepx, x) && leq(it->gety(x), a->gety(x))) { pq.emplace(x, -it->id, it); } }; auto add_edge = [&](int i, int j) { if (i != -1 && j != -1) edges.emplace_back(i, j); }; // rotate all points by a radian // so that x coordinates are distinct for (auto& p : pts) p = rot(p, 1); // sentinel beach.emplace(Vec(-INF, INF), Vec(-INF * 2, INF), -1); // add all point events for (int i = 0; i < n; ++i) { pq.emplace(pts[i].real(), i, beach.end()); } while (!pq.empty()) { auto e = pq.top(); pq.pop(); sweepx = e.x; if (e.id >= 0) { // point event int i = e.id; auto p = pts[i]; // (c) -> (c, p, c) auto c = beach.lower_bound(p.imag()); auto b = beach.insert(c, Parabola(p, c->p, i)); auto a = beach.insert(b, Parabola(c->p, p, c->i)); add_edge(i, c->i); update(a); update(b); update(c); } else if (valid[-e.id]) { // vertex event // (a, e, b) -> (a, b) auto a = std::prev(e.it); auto b = std::next(e.it); beach.erase(e.it); a->q = b->p; add_edge(a->i, b->i); update(a); update(b); } } return edges; }
#line 2 "geometry/delaunay_diagram.hpp" #include <queue> #include <set> #include <utility> #include <vector> #line 2 "geometry/geometry.hpp" #include <algorithm> #include <cassert> #include <cmath> #include <complex> #include <iostream> #include <numbers> #include <numeric> #line 10 "geometry/geometry.hpp" // note that if T is of an integer type, std::abs does not work using T = double; using Vec = std::complex<T>; std::istream& operator>>(std::istream& is, Vec& p) { T x, y; is >> x >> y; p = {x, y}; return is; } T dot(const Vec& a, const Vec& b) { return (std::conj(a) * b).real(); } T cross(const Vec& a, const Vec& b) { return (std::conj(a) * b).imag(); } constexpr T PI = std::numbers::pi_v<T>; constexpr T eps = 1e-10; inline bool eq(T a, T b) { return std::abs(a - b) <= eps; } inline bool eq(Vec a, Vec b) { return std::abs(a - b) <= eps; } inline bool lt(T a, T b) { return a < b - eps; } inline bool leq(T a, T b) { return a <= b + eps; } struct Line { Vec p1, p2; Line() = default; Line(const Vec& p1, const Vec& p2) : p1(p1), p2(p2) {} Vec dir() const { return p2 - p1; } }; struct Segment : Line { using Line::Line; }; struct Circle { Vec c; T r; Circle() = default; Circle(const Vec& c, T r) : c(c), r(r) {} }; using Polygon = std::vector<Vec>; Vec rot(const Vec& a, T ang) { return a * Vec(std::cos(ang), std::sin(ang)); } Vec perp(const Vec& a) { return Vec(-a.imag(), a.real()); } Vec projection(const Line& l, const Vec& p) { return l.p1 + dot(p - l.p1, l.dir()) * l.dir() / std::norm(l.dir()); } Vec reflection(const Line& l, const Vec& p) { return T(2) * projection(l, p) - p; } // 0: collinear // 1: counter-clockwise // -1: clockwise int ccw(const Vec& a, const Vec& b, const Vec& c) { if (eq(cross(b - a, c - a), 0)) return 0; if (lt(cross(b - a, c - a), 0)) return -1; return 1; } void sort_by_arg(std::vector<Vec>& pts) { std::ranges::sort(pts, [&](auto& p, auto& q) { if ((p.imag() < 0) != (q.imag() < 0)) return (p.imag() < 0); if (cross(p, q) == 0) { if (p == Vec(0, 0)) return !(q.imag() < 0 || (q.imag() == 0 && q.real() > 0)); if (q == Vec(0, 0)) return (p.imag() < 0 || (p.imag() == 0 && p.real() > 0)); return (p.real() > q.real()); } return (cross(p, q) > 0); }); } #line 3 "geometry/intersect.hpp" bool intersect(const Segment& s, const Vec& p) { Vec u = s.p1 - p, v = s.p2 - p; return eq(cross(u, v), 0) && leq(dot(u, v), 0); } // 0: outside // 1: on the border // 2: inside int intersect(const Polygon& poly, const Vec& p) { const int n = poly.size(); bool in = 0; for (int i = 0; i < n; ++i) { auto a = poly[i] - p, b = poly[(i + 1) % n] - p; if (eq(cross(a, b), 0) && (lt(dot(a, b), 0) || eq(dot(a, b), 0))) return 1; if (a.imag() > b.imag()) std::swap(a, b); if (leq(a.imag(), 0) && lt(0, b.imag()) && lt(cross(a, b), 0)) in ^= 1; } return in ? 2 : 0; } int intersect(const Segment& s, const Segment& t) { auto a = s.p1, b = s.p2; auto c = t.p1, d = t.p2; if (ccw(a, b, c) != ccw(a, b, d) && ccw(c, d, a) != ccw(c, d, b)) return 2; if (intersect(s, c) || intersect(s, d) || intersect(t, a) || intersect(t, b)) return 1; return 0; } // true if they have positive area in common or touch on the border bool intersect(const Polygon& poly1, const Polygon& poly2) { const int n = poly1.size(); const int m = poly2.size(); for (int i = 0; i < n; ++i) { for (int j = 0; j < m; ++j) { if (intersect(Segment(poly1[i], poly1[(i + 1) % n]), Segment(poly2[j], poly2[(j + 1) % m]))) { return true; } } } return intersect(poly1, poly2[0]) || intersect(poly2, poly1[0]); } // 0: inside // 1: inscribe // 2: intersect // 3: circumscribe // 4: outside int intersect(const Circle& c1, const Circle& c2) { T d = std::abs(c1.c - c2.c); if (lt(d, std::abs(c2.r - c1.r))) return 0; if (eq(d, std::abs(c2.r - c1.r))) return 1; if (eq(c1.r + c2.r, d)) return 3; if (lt(c1.r + c2.r, d)) return 4; return 2; } #line 4 "geometry/dist.hpp" T dist(const Line& l, const Vec& p) { return std::abs(cross(p - l.p1, l.dir())) / std::abs(l.dir()); } T dist(const Segment& s, const Vec& p) { if (lt(dot(p - s.p1, s.dir()), 0)) return std::abs(p - s.p1); if (lt(dot(p - s.p2, -s.dir()), 0)) return std::abs(p - s.p2); return std::abs(cross(p - s.p1, s.dir())) / std::abs(s.dir()); } T dist(const Segment& s, const Segment& t) { if (intersect(s, t)) return T(0); return std::min( {dist(s, t.p1), dist(s, t.p2), dist(t, s.p1), dist(t, s.p2)}); } #line 4 "geometry/intersection.hpp" Vec intersection(const Line& l, const Line& m) { assert(!eq(cross(l.dir(), m.dir()), 0)); // not parallel Vec r = m.p1 - l.p1; return l.p1 + cross(m.dir(), r) / cross(m.dir(), l.dir()) * l.dir(); } std::vector<Vec> intersection(const Circle& c, const Line& l) { T d = dist(l, c.c); if (lt(c.r, d)) return {}; // no intersection Vec e1 = l.dir() / std::abs(l.dir()); Vec e2 = perp(e1); if (ccw(c.c, l.p1, l.p2) == 1) e2 *= -1; if (eq(c.r, d)) return {c.c + d * e2}; // tangent T t = std::sqrt(c.r * c.r - d * d); return {c.c + d * e2 + t * e1, c.c + d * e2 - t * e1}; } std::vector<Vec> intersection(const Circle& c1, const Circle& c2) { T d = std::abs(c1.c - c2.c); if (lt(c1.r + c2.r, d)) return {}; // outside Vec e1 = (c2.c - c1.c) / std::abs(c2.c - c1.c); Vec e2 = perp(e1); if (lt(d, std::abs(c2.r - c1.r))) return {}; // contain if (eq(d, std::abs(c2.r - c1.r))) return {c1.c + c1.r * e1}; // tangent T x = (c1.r * c1.r - c2.r * c2.r + d * d) / (2 * d); T y = std::sqrt(c1.r * c1.r - x * x); return {c1.c + x * e1 + y * e2, c1.c + x * e1 - y * e2}; } T area_intersection(const Circle& c1, const Circle& c2) { T d = std::abs(c2.c - c1.c); if (leq(c1.r + c2.r, d)) return 0; // outside if (leq(d, std::abs(c2.r - c1.r))) { // inside T r = std::min(c1.r, c2.r); return PI * r * r; } T ans = 0; T a; a = std::acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2 * c1.r * d)); ans += c1.r * c1.r * (a - std::sin(a) * std::cos(a)); a = std::acos((c2.r * c2.r + d * d - c1.r * c1.r) / (2 * c2.r * d)); ans += c2.r * c2.r * (a - std::sin(a) * std::cos(a)); return ans; } #line 4 "geometry/bisector.hpp" Line bisector(const Segment& s) { auto m = (s.p1 + s.p2) / T(2); return Line(m, m + Vec(-s.dir().imag(), s.dir().real())); } std::pair<Line, Line> bisector(const Line& l, const Line& m) { // parallel if (eq(cross(l.dir(), m.dir()), 0)) { auto n = Line(l.p1, l.p1 + perp(l.dir())); auto p = intersection(n, m); auto m = (l.p1 + p) / T(2); return {Line(m, m + l.dir()), Line()}; } auto p = intersection(l, m); T ang = (std::arg(l.dir()) + std::arg(m.dir())) / T(2); auto b1 = Line(p, p + std::polar(T(1), ang)); auto b2 = Line(p, p + std::polar(T(1), ang + PI / T(2))); return {b1, b2}; } #line 5 "geometry/triangle.hpp" Vec centroid(const Vec& A, const Vec& B, const Vec& C) { assert(ccw(A, B, C) != 0); return (A + B + C) / T(3); } Vec incenter(const Vec& A, const Vec& B, const Vec& C) { assert(ccw(A, B, C) != 0); T a = std::abs(B - C); T b = std::abs(C - A); T c = std::abs(A - B); return (a * A + b * B + c * C) / (a + b + c); } Vec circumcenter(const Vec& A, const Vec& B, const Vec& C) { assert(ccw(A, B, C) != 0); return intersection(bisector(Segment(A, B)), bisector(Segment(A, C))); } // large error but beautiful // Vec circumcenter(const Vec& A, const Vec& B, const Vec& C) { // assert(ccw(A, B, C) != 0); // Vec p = C - B, q = A - C, r = B - A; // T a = std::norm(p) * dot(q, r); // T b = std::norm(q) * dot(r, p); // T c = std::norm(r) * dot(p, q); // return (a*A + b*B + c*C) / (a + b + c); // } #line 9 "geometry/delaunay_diagram.hpp" std::vector<std::pair<int, int>> delaunay_diagram(std::vector<Vec> pts) { const int n = pts.size(); constexpr T INF = 1e20; // stores delaunay edges std::vector<std::pair<int, int>> edges; // the sweep line moves from left to right // x coordinate of the sweepline static T sweepx; // check if all points are collinear bool ok = false; for (int i = 2; i < n; ++i) { if (ccw(pts[0], pts[1], pts[i]) != 0) { ok = true; break; } } // handle degenerate cases if (!ok) { // connect all points by a path std::vector<std::pair<Vec, int>> pi(n); for (int i = 0; i < n; ++i) pi[i] = {pts[i], i}; std::ranges::sort(pi, {}, [&](auto& p) { return std::make_pair(p.first.real(), p.first.imag()); }); for (int k = 0; k < n - 1; ++k) { edges.emplace_back(pi[k].second, pi[k + 1].second); } return edges; } // represents a parabola given by the focus p // and the focus q of the next parabola struct Parabola { Vec p; mutable Vec q; int i; mutable int id = 0; Parabola(const Vec& p, const Vec& q, int i) : p(p), q(q), i(i) {} T gety(T x) const { if (q.imag() == INF) return INF; x += eps; Vec m = (p + q) / T(2); Vec dir = rot(p - m, PI / 2); T D = (x - p.real()) * (x - q.real()); return m.imag() + ((m.real() - x) * dir.real() + std::sqrt(D) * std::abs(dir)) / dir.imag(); } bool operator<(T y) const { return gety(sweepx) < y; } bool operator<(const Parabola& o) const { return gety(sweepx) < o.gety(sweepx); } }; // maintains a list of parabola std::multiset<Parabola, std::less<>> beach; using iterator = decltype(beach)::iterator; // represents an event // id >= 0: point event // id < 0: vertex event struct Event { T x; int id; iterator it; Event(T x, int id, iterator it) : x(x), id(id), it(it) {} bool operator>(const Event& e) const { return x > e.x; } }; // maintains a list of events std::priority_queue<Event, std::vector<Event>, std::greater<>> pq; // true if the i-th vertex event is still valid std::vector<bool> valid(1); int k = 1; // next id // create vertex events auto update = [&](const iterator& it) { if (it->i <= -1 || it == beach.begin()) return; // sentinels valid[it->id] = false; auto a = std::prev(it); // handle a vertex event that occurs // when 3 parabolas at (a->p, it->p, it->q) coincide if (eq(ccw(a->p, it->p, it->q), 0)) return; // never coincide when collinear // create a new event it->id = k++; valid.push_back(true); Vec c = circumcenter(a->p, it->p, it->q); // coordinate of the vertex event T x = c.real() + std::abs(c - it->p); if (leq(sweepx, x) && leq(it->gety(x), a->gety(x))) { pq.emplace(x, -it->id, it); } }; auto add_edge = [&](int i, int j) { if (i != -1 && j != -1) edges.emplace_back(i, j); }; // rotate all points by a radian // so that x coordinates are distinct for (auto& p : pts) p = rot(p, 1); // sentinel beach.emplace(Vec(-INF, INF), Vec(-INF * 2, INF), -1); // add all point events for (int i = 0; i < n; ++i) { pq.emplace(pts[i].real(), i, beach.end()); } while (!pq.empty()) { auto e = pq.top(); pq.pop(); sweepx = e.x; if (e.id >= 0) { // point event int i = e.id; auto p = pts[i]; // (c) -> (c, p, c) auto c = beach.lower_bound(p.imag()); auto b = beach.insert(c, Parabola(p, c->p, i)); auto a = beach.insert(b, Parabola(c->p, p, c->i)); add_edge(i, c->i); update(a); update(b); update(c); } else if (valid[-e.id]) { // vertex event // (a, e, b) -> (a, b) auto a = std::prev(e.it); auto b = std::next(e.it); beach.erase(e.it); a->q = b->p; add_edge(a->i, b->i); update(a); update(b); } } return edges; }