sotanishy's competitive programming library

sotanishy's code snippets for competitive programming

View the Project on GitHub sotanishy/cp-library-cpp

:warning: Delaunay Diagram
(geometry/delaunay_diagram.hpp)

Description

与えられた点の Delaunay 図は，Voronoi 図の双対であり，Voronoi 領域の隣接関係を表す．

Voronoi 図とは，平面に配置された母点に対して，平面の各点がどの母点に最も近いかによって領域分けされた図である．

Delaunay 辺の本数は Euler の多面体定理より $O(n)$ 本である．また，Euclidean MST は Delaunay 図の部分グラフである．

Delaunay 図は平面走査アルゴリズムである Fortune のアルゴリズムで求めることができる．

Operations

Reference

Depends on

Code

#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;
}
Back to top page