Delaunay Triangulation
(geometry/delaunay_triangulation.hpp)

• View this file on GitHub
• Last update: 2026-10-05 22:23:07+09:00
• Include: #include "geometry/delaunay_triangulation.hpp"

Overview

T may be a built-in integer or an exact coordinate class such as math::Rational<long long> or math::Rational<utilities::BigInt>. Rational intermediate arithmetic uses T without floating-point conversion. All intermediate fractions must be representable. Listed complexities count scalar operations; rational arithmetic adds its gcd and integer arithmetic costs.

delaunay_triangulation constructs one Delaunay triangulation of distinct two-dimensional exact-coordinate points. It returns both the undirected edges and the counterclockwise triangular faces, using the original zero-based point indices.

The implementation uses divide and conquer. All topological decisions are made with exact orientation and incircle predicates.

Type

struct DelaunayTriangulation {
    std::vector<std::pair<int, int>> edges;
    std::vector<std::array<int, 3>> triangles;
};

Every edge is stored as (first, second) with first < second, and edges is lexicographically sorted without duplicates.

Every triangle is counterclockwise. Its smallest vertex index is stored first, and triangles is lexicographically sorted without duplicates. Every side of a triangle appears in edges.

Function

The exact signature is:

template <ExactCoordinate T>
DelaunayTriangulation delaunay_triangulation(
    const std::vector<Point<T>>& points
);

Function	Description	Complexity
delaunay_triangulation(points)	Returns one Delaunay triangulation of the indexed points.	$O(N\log N)$ time and $O(N)$ memory

The input points must be pairwise distinct. Empty and one-point inputs have no edges or triangles. Two points produce one edge. A larger collinear input has no triangles and connects consecutive points in lexicographic order.

When four or more points are cocircular, the Delaunay triangulation is not unique. The function returns one valid choice of diagonals; callers should not depend on which valid choice is selected.

The degree-four incircle expressions must fit in signed 128-bit arithmetic. This is satisfied by the coordinate constraints of Library Checker’s Euclidean MST problem.

Example

#include "geometry/delaunay_triangulation.hpp"

#include <iostream>
#include <vector>

int main() {
    using Point = m1une::geometry::Point<long long>;
    std::vector<Point> points;
    points.emplace_back(0, 0);
    points.emplace_back(4, 0);
    points.emplace_back(0, 3);
    points.emplace_back(1, 1);

    auto triangulation =
        m1une::geometry::delaunay_triangulation(points);
    for (const auto& triangle : triangulation.triangles) {
        std::cout << triangle[0] << ' '
                  << triangle[1] << ' '
                  << triangle[2] << '\n';
    }
}

Depends on

• DSU (Disjoint Set Union) (ds/dsu/dsu.hpp)
• geometry/detail/floating_predicate.hpp
• Euclidean Minimum Spanning Tree (geometry/euclidean_mst.hpp)
• 2D Point and Predicates (geometry/point.hpp)

Required by

• Geometry Bundle (geometry/all.hpp)

Verified with

• verify/geometry/centroid.test.cpp
• verify/geometry/delaunay_triangulation.test.cpp
• verify/geometry/geometry_algorithms.test.cpp
• verify/geometry/rational.test.cpp

Code

#ifndef M1UNE_GEOMETRY_DELAUNAY_TRIANGULATION_HPP
#define M1UNE_GEOMETRY_DELAUNAY_TRIANGULATION_HPP 1

#include <algorithm>
#include <array>
#include <cassert>
#include <concepts>
#include <utility>
#include <vector>

#include "euclidean_mst.hpp"

namespace m1une {
namespace geometry {

struct DelaunayTriangulation {
    std::vector<std::pair<int, int>> edges;
    std::vector<std::array<int, 3>> triangles;
};

namespace delaunay_triangulation_detail {

template <ExactCoordinate T>
int direction_half(
    const Point<T>& origin,
    const Point<T>& destination
) {
    using W = wide_type<T>;
    W x = W(destination.x) - W(origin.x);
    W y = W(destination.y) - W(origin.y);
    return y > 0 || (y == 0 && x >= 0) ? 0 : 1;
}

template <ExactCoordinate T>
bool direction_less(
    const std::vector<Point<T>>& points,
    int origin,
    int first,
    int second
) {
    int first_half = direction_half(points[origin], points[first]);
    int second_half = direction_half(points[origin], points[second]);
    if (first_half != second_half) return first_half < second_half;

    using W = wide_type<T>;
    W first_x = W(points[first].x) - W(points[origin].x);
    W first_y = W(points[first].y) - W(points[origin].y);
    W second_x = W(points[second].x) - W(points[origin].x);
    W second_y = W(points[second].y) - W(points[origin].y);
    W product = first_x * second_y - first_y * second_x;
    if (product != 0) return product > 0;

    W first_norm = first_x * first_x + first_y * first_y;
    W second_norm = second_x * second_x + second_y * second_y;
    if (first_norm != second_norm) return first_norm < second_norm;
    return first < second;
}

inline void rotate_minimum_first(std::array<int, 3>& triangle) {
    int minimum = int(std::min_element(triangle.begin(), triangle.end()) -
                      triangle.begin());
    std::rotate(
        triangle.begin(),
        triangle.begin() + minimum,
        triangle.end()
    );
}

}  // namespace delaunay_triangulation_detail

// Constructs one Delaunay triangulation of distinct exact-coordinate points.
template <ExactCoordinate T>
DelaunayTriangulation delaunay_triangulation(
    const std::vector<Point<T>>& points
) {
    namespace detail = delaunay_triangulation_detail;

    DelaunayTriangulation result;
    geometry::detail::EuclideanDelaunay<T> builder(points);
    assert(!builder.has_duplicates());
    result.edges = builder.get_edges();
    for (auto& [first, second] : result.edges) {
        if (first > second) std::swap(first, second);
    }
    std::sort(result.edges.begin(), result.edges.end());
    result.edges.erase(
        std::unique(result.edges.begin(), result.edges.end()),
        result.edges.end()
    );

    std::vector<std::vector<int>> neighbors(points.size());
    for (auto [first, second] : result.edges) {
        neighbors[first].push_back(second);
        neighbors[second].push_back(first);
    }
    for (int point = 0; point < int(points.size()); ++point) {
        std::sort(
            neighbors[point].begin(),
            neighbors[point].end(),
            [&](int first, int second) {
                return detail::direction_less(points, point, first, second);
            }
        );
    }

    auto has_edge = [&](int first, int second) {
        if (first > second) std::swap(first, second);
        return std::binary_search(
            result.edges.begin(),
            result.edges.end(),
            std::pair(first, second)
        );
    };
    result.triangles.reserve(result.edges.size());
    for (int point = 0; point < int(points.size()); ++point) {
        int degree = int(neighbors[point].size());
        for (int index = 0; index < degree; ++index) {
            int first = neighbors[point][index];
            int second = neighbors[point][(index + 1) % degree];
            if (orientation(points[point], points[first], points[second]) <= 0) {
                continue;
            }
            if (!has_edge(first, second)) continue;
            std::array<int, 3> triangle{point, first, second};
            detail::rotate_minimum_first(triangle);
            result.triangles.push_back(triangle);
        }
    }
    std::sort(result.triangles.begin(), result.triangles.end());
    result.triangles.erase(
        std::unique(result.triangles.begin(), result.triangles.end()),
        result.triangles.end()
    );
    return result;
}

}  // namespace geometry
}  // namespace m1une

#endif  // M1UNE_GEOMETRY_DELAUNAY_TRIANGULATION_HPP

#line 1 "geometry/delaunay_triangulation.hpp"



#include <algorithm>
#include <array>
#include <cassert>
#include <concepts>
#include <utility>
#include <vector>

#line 1 "geometry/euclidean_mst.hpp"



#line 6 "geometry/euclidean_mst.hpp"
#include <cmath>
#line 8 "geometry/euclidean_mst.hpp"
#include <cstddef>
#include <limits>
#include <tuple>
#line 13 "geometry/euclidean_mst.hpp"

#line 1 "ds/dsu/dsu.hpp"



#line 5 "ds/dsu/dsu.hpp"
#include <numeric>
#line 8 "ds/dsu/dsu.hpp"

namespace m1une {
namespace ds {

struct Dsu {
   private:
    int _n;
    // parent_or_size[i] is the parent of i if it's >= 0.
    // If it's < 0, then i is a root and -parent_or_size[i] is the size of the group.
    std::vector<int> parent_or_size;

    // Returns {new leader, absorbed leader}. The absorbed leader is -1 when
    // both vertices already belong to the same component.
    std::pair<int, int> merge_leaders(int a, int b) {
        int x = leader(a), y = leader(b);
        if (x == y) return {x, -1};
        if (-parent_or_size[x] < -parent_or_size[y]) std::swap(x, y);
        parent_or_size[x] += parent_or_size[y];
        parent_or_size[y] = x;
        return {x, y};
    }

   public:
    Dsu() : _n(0) {}
    explicit Dsu(int n) : _n(n), parent_or_size(n, -1) {}

    // Merges the group containing 'a' with the group containing 'b'.
    // Returns the leader of the merged group.
    int merge(int a, int b) {
        return merge_leaders(a, b).first;
    }

    // Invokes callback(new_leader, absorbed_leader) after an actual merge.
    // Returns the leader of the merged group.
    template <class Callback>
    int merge(int a, int b, Callback&& callback) {
        std::pair<int, int> merged = merge_leaders(a, b);
        if (merged.second != -1) callback(merged.first, merged.second);
        return merged.first;
    }

    // Returns true if 'a' and 'b' belong to the same group.
    bool same(int a, int b) {
        return leader(a) == leader(b);
    }

    // Returns the leader (representative) of the group containing 'a'.
    int leader(int a) {
        if (parent_or_size[a] < 0) return a;
        // Path compression
        return parent_or_size[a] = leader(parent_or_size[a]);
    }

    // Returns the size of the group containing 'a'.
    int size(int a) {
        return -parent_or_size[leader(a)];
    }

    // Returns a list of all groups, where each group is a vector of its elements.
    std::vector<std::vector<int>> groups() {
        std::vector<int> leader_buf(_n), group_size(_n);
        for (int i = 0; i < _n; i++) {
            leader_buf[i] = leader(i);
            group_size[leader_buf[i]]++;
        }
        std::vector<std::vector<int>> result(_n);
        for (int i = 0; i < _n; i++) {
            result[i].reserve(group_size[i]);
        }
        for (int i = 0; i < _n; i++) {
            result[leader_buf[i]].push_back(i);
        }
        result.erase(std::remove_if(result.begin(), result.end(), [&](const std::vector<int>& v) { return v.empty(); }),
                     result.end());
        return result;
    }
};

}  // namespace ds
}  // namespace m1une


#line 1 "geometry/point.hpp"



#line 7 "geometry/point.hpp"
#include <type_traits>

#line 1 "geometry/detail/floating_predicate.hpp"



namespace m1une {
namespace geometry {
namespace predicate_detail {

template <typename T>
constexpr T absolute(T value) {
    return value < T(0) ? -value : value;
}

template <typename T>
constexpr T max_value(T first, T second) {
    return first < second ? second : first;
}

template <typename T>
constexpr T vector_scale(T x, T y) {
    return max_value(absolute(x), absolute(y));
}

template <bool Exact, typename T>
constexpr int scaled_sign(T value, T scale, long double eps) {
    if constexpr (Exact) {
        return (value > T(0)) - (value < T(0));
    } else {
        const T tolerance = T(eps) * scale;
        return (value > tolerance) - (value < -tolerance);
    }
}

template <bool Exact, typename T>
constexpr T determinant_scale(T ax, T ay, T bx, T by) {
    if constexpr (Exact) {
        return T(0);
    } else {
        return vector_scale(ax, ay) * vector_scale(bx, by);
    }
}

template <bool Exact, typename T>
constexpr int determinant_sign(
    T ax,
    T ay,
    T bx,
    T by,
    long double eps
) {
    const T determinant = ax * by - ay * bx;
    return scaled_sign<Exact>(
        determinant,
        determinant_scale<Exact>(ax, ay, bx, by),
        eps
    );
}

template <bool Exact, typename T>
constexpr int orientation_sign(
    T direction_x,
    T direction_y,
    T offset_x,
    T offset_y,
    long double eps
) {
    const T determinant =
        direction_x * offset_y - direction_y * offset_x;
    T scale = T(0);
    if constexpr (!Exact) {
        const T direction_scale =
            vector_scale(direction_x, direction_y);
        scale = direction_scale * max_value(
            direction_scale,
            vector_scale(offset_x, offset_y)
        );
    }
    return scaled_sign<Exact>(determinant, scale, eps);
}

template <bool Exact, typename T>
constexpr int dot_sign(
    T ax,
    T ay,
    T bx,
    T by,
    long double eps
) {
    const T value = ax * bx + ay * by;
    T scale = T(0);
    if constexpr (!Exact) {
        scale = vector_scale(ax, ay) * vector_scale(bx, by);
    }
    return scaled_sign<Exact>(value, scale, eps);
}

}  // namespace predicate_detail
}  // namespace geometry
}  // namespace m1une


#line 10 "geometry/point.hpp"

namespace m1une {
namespace geometry {

template <typename T>
concept Coordinate = !std::same_as<std::remove_cv_t<T>, bool> &&
    (std::is_arithmetic_v<T> ||
     (std::copyable<T> && std::totally_ordered<T> && requires(T a, T b) {
         T(0);
         T(1);
         static_cast<long double>(a);
         { +a } -> std::same_as<T>;
         { -a } -> std::same_as<T>;
         { a + b } -> std::same_as<T>;
         { a - b } -> std::same_as<T>;
         { a * b } -> std::same_as<T>;
         { a / b } -> std::same_as<T>;
         { a += b } -> std::same_as<T&>;
         { a -= b } -> std::same_as<T&>;
     }));

// Custom coordinate types keep their own exact arithmetic.
template <typename T>
concept ExactCoordinate = Coordinate<T> && !std::floating_point<T>;

template <Coordinate T>
using wide_type = std::conditional_t<std::integral<T>, __int128_t,
    std::conditional_t<std::floating_point<T>, long double, T>>;

template <Coordinate T>
struct Point {
    T x;
    T y;

    constexpr Point() : x(0), y(0) {}
    constexpr Point(T x_value, T y_value) : x(x_value), y(y_value) {}

    template <Coordinate U>
    explicit constexpr Point(const Point<U>& other)
        : x(static_cast<T>(other.x)), y(static_cast<T>(other.y)) {}

    constexpr Point& operator+=(const Point& other) {
        x += other.x;
        y += other.y;
        return *this;
    }

    constexpr Point& operator-=(const Point& other) {
        x -= other.x;
        y -= other.y;
        return *this;
    }

    constexpr Point operator+() const {
        return *this;
    }

    constexpr Point operator-() const {
        return Point(-x, -y);
    }

    friend constexpr Point operator+(Point left, const Point& right) {
        return left += right;
    }

    friend constexpr Point operator-(Point left, const Point& right) {
        return left -= right;
    }

    friend constexpr bool operator==(const Point&, const Point&) = default;

    friend constexpr bool operator<(const Point& left, const Point& right) {
        if (left.x != right.x) return left.x < right.x;
        return left.y < right.y;
    }
};

template <Coordinate T>
constexpr Point<long double> centroid(const Point<T>& point) {
    return Point<long double>(point);
}

template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(const Point<T>& point, Scalar scalar) {
    using Result = std::common_type_t<T, Scalar>;
    return Point<Result>(
        Result(point.x) * Result(scalar),
        Result(point.y) * Result(scalar)
    );
}

template <typename Scalar, Coordinate T>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(Scalar scalar, const Point<T>& point) {
    return point * scalar;
}

template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator/(const Point<T>& point, Scalar scalar) {
    using Result = std::common_type_t<T, Scalar>;
    return Point<Result>(
        Result(point.x) / Result(scalar),
        Result(point.y) / Result(scalar)
    );
}

template <Coordinate T>
constexpr wide_type<T> dot(const Point<T>& a, const Point<T>& b) {
    using W = wide_type<T>;
    return W(a.x) * W(b.x) + W(a.y) * W(b.y);
}

template <Coordinate T>
constexpr wide_type<T> cross(const Point<T>& a, const Point<T>& b) {
    using W = wide_type<T>;
    return W(a.x) * W(b.y) - W(a.y) * W(b.x);
}

template <Coordinate T>
constexpr wide_type<T> cross(
    const Point<T>& origin,
    const Point<T>& a,
    const Point<T>& b
) {
    using W = wide_type<T>;
    W ax = W(a.x) - W(origin.x);
    W ay = W(a.y) - W(origin.y);
    W bx = W(b.x) - W(origin.x);
    W by = W(b.y) - W(origin.y);
    return ax * by - ay * bx;
}

template <Coordinate T>
constexpr wide_type<T> norm2(const Point<T>& point) {
    return dot(point, point);
}

template <Coordinate T>
constexpr wide_type<T> distance2(const Point<T>& a, const Point<T>& b) {
    using W = wide_type<T>;
    W dx = W(a.x) - W(b.x);
    W dy = W(a.y) - W(b.y);
    return dx * dx + dy * dy;
}

template <Coordinate T>
long double norm(const Point<T>& point) {
    return std::hypot(
        static_cast<long double>(point.x),
        static_cast<long double>(point.y)
    );
}

template <Coordinate T>
long double distance(const Point<T>& a, const Point<T>& b) {
    return std::hypot(
        static_cast<long double>(a.x) - static_cast<long double>(b.x),
        static_cast<long double>(a.y) - static_cast<long double>(b.y)
    );
}

template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
         (std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> internal_division_point(
    const Point<T>& a,
    const Point<T>& b,
    M m,
    N n
) {
    long double first_ratio = static_cast<long double>(m);
    long double second_ratio = static_cast<long double>(n);
    long double denominator = first_ratio + second_ratio;
    assert(denominator != 0);
    Point<long double> first(a);
    Point<long double> direction = Point<long double>(b) - first;
    return first + direction * (first_ratio / denominator);
}

template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
         (std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> external_division_point(
    const Point<T>& a,
    const Point<T>& b,
    M m,
    N n
) {
    long double first_ratio = static_cast<long double>(m);
    long double second_ratio = static_cast<long double>(n);
    long double denominator = first_ratio - second_ratio;
    assert(denominator != 0);
    Point<long double> first(a);
    Point<long double> direction = Point<long double>(b) - first;
    return first + direction * (first_ratio / denominator);
}

template <Coordinate T>
constexpr int sign(wide_type<T> value, long double eps = 1e-12L) {
    return predicate_detail::scaled_sign<ExactCoordinate<T>>(
        value,
        wide_type<T>(1),
        eps
    );
}

template <Coordinate T>
constexpr int orientation(
    const Point<T>& a,
    const Point<T>& b,
    const Point<T>& c,
    long double eps = 1e-12L
) {
    using W = wide_type<T>;
    const W first_x = W(b.x) - W(a.x);
    const W first_y = W(b.y) - W(a.y);
    const W second_x = W(c.x) - W(a.x);
    const W second_y = W(c.y) - W(a.y);
    return predicate_detail::orientation_sign<ExactCoordinate<T>>(
        first_x,
        first_y,
        second_x,
        second_y,
        eps
    );
}

template <Coordinate T>
constexpr bool collinear(
    const Point<T>& a,
    const Point<T>& b,
    const Point<T>& c,
    long double eps = 1e-12L
) {
    return orientation(a, b, c, eps) == 0;
}

template <Coordinate T>
Point<long double> rotate(const Point<T>& point, long double angle) {
    long double cosine = std::cos(angle);
    long double sine = std::sin(angle);
    return Point<long double>(
        static_cast<long double>(point.x) * cosine -
            static_cast<long double>(point.y) * sine,
        static_cast<long double>(point.x) * sine +
            static_cast<long double>(point.y) * cosine
    );
}

template <Coordinate T>
Point<long double> normalized(const Point<T>& point) {
    long double length = norm(point);
    assert(length != 0);
    return Point<long double>(
        static_cast<long double>(point.x) / length,
        static_cast<long double>(point.y) / length
    );
}

}  // namespace geometry
}  // namespace m1une


#line 16 "geometry/euclidean_mst.hpp"

namespace m1une {
namespace geometry {

template <class T>
struct EuclideanMstEdge {
    int from;
    int to;
    T squared_distance;
};

template <class T>
struct EuclideanMst {
    long double cost;
    std::vector<EuclideanMstEdge<T>> edges;
};

namespace detail {

template <ExactCoordinate T>
class EuclideanDelaunay {
   private:
    using W = wide_type<T>;

    struct InternalPoint {
        W x;
        W y;

        friend bool operator==(const InternalPoint&, const InternalPoint&) = default;
    };

    struct Edge {
        int to;
        int ccw;
        int cw;
        int reverse;
        bool enabled = false;
    };

    std::vector<int> open_addresses;
    std::vector<InternalPoint> points;
    std::vector<Edge> edges;
    std::vector<int> duplicate_representative;

    static InternalPoint subtract(const InternalPoint& a, const InternalPoint& b) {
        return InternalPoint{a.x - b.x, a.y - b.y};
    }

    static W cross_product(const InternalPoint& a, const InternalPoint& b) {
        return a.x * b.y - a.y * b.x;
    }

    static W squared_norm(const InternalPoint& point) {
        return point.x * point.x + point.y * point.y;
    }

    static bool inside_circumcircle(
        InternalPoint a,
        InternalPoint b,
        InternalPoint c,
        const InternalPoint& d
    ) {
        a = subtract(a, d);
        b = subtract(b, d);
        c = subtract(c, d);
        W determinant = cross_product(b, c) * squared_norm(a)
                      + cross_product(c, a) * squared_norm(b)
                      + cross_product(a, b) * squared_norm(c);
        return determinant > 0;
    }

    int get_open_address() {
        if (open_addresses.empty()) {
            edges.push_back(Edge());
            return int(edges.size()) - 1;
        }
        int result = open_addresses.back();
        open_addresses.pop_back();
        return result;
    }

    std::pair<int, int> add_edge(int from, int to) {
        int forward = get_open_address();
        int backward = get_open_address();
        edges[forward].to = to;
        edges[forward].ccw = forward;
        edges[forward].cw = forward;
        edges[forward].reverse = backward;
        edges[forward].enabled = true;
        edges[backward].to = from;
        edges[backward].ccw = backward;
        edges[backward].cw = backward;
        edges[backward].reverse = forward;
        edges[backward].enabled = true;
        return {forward, backward};
    }

    void erase_directed_edge(int edge) {
        int ccw = edges[edge].ccw;
        int cw = edges[edge].cw;
        edges[ccw].cw = cw;
        edges[cw].ccw = ccw;
        edges[edge].enabled = false;
    }

    void erase_edge(int edge) {
        int reverse = edges[edge].reverse;
        erase_directed_edge(edge);
        erase_directed_edge(reverse);
        open_addresses.push_back(edge);
        open_addresses.push_back(reverse);
    }

    void insert_ccw_after(int edge, int position) {
        int next = edges[position].ccw;
        edges[edge].ccw = next;
        edges[next].cw = edge;
        edges[edge].cw = position;
        edges[position].ccw = edge;
    }

    void insert_cw_after(int edge, int position) {
        int next = edges[position].cw;
        edges[edge].cw = next;
        edges[next].ccw = edge;
        edges[edge].ccw = position;
        edges[position].cw = edge;
    }

    int orientation(int a, int b, int c) const {
        InternalPoint ab = subtract(points[b], points[a]);
        InternalPoint ac = subtract(points[c], points[a]);
        W value = cross_product(ab, ac);
        return (value > 0) - (value < 0);
    }

    std::pair<int, int> go_next(int edge) const {
        int vertex = edges[edge].to;
        int next_edge = edges[edges[edge].reverse].ccw;
        return {vertex, next_edge};
    }

    std::pair<int, int> go_previous(int edge) const {
        int vertex = edges[edges[edge].cw].to;
        int next_edge = edges[edges[edge].cw].reverse;
        return {vertex, next_edge};
    }

    std::tuple<int, int, int, int> lower_tangent(
        int left_vertex,
        int left_edge,
        int right_vertex,
        int right_edge
    ) const {
        while (true) {
            auto [next_left_vertex, next_left_edge] = go_previous(left_edge);
            if (orientation(right_vertex, left_vertex, next_left_vertex) > 0) {
                left_vertex = next_left_vertex;
                left_edge = next_left_edge;
                continue;
            }
            auto [next_right_vertex, next_right_edge] = go_next(right_edge);
            if (orientation(left_vertex, right_vertex, next_right_vertex) < 0) {
                right_vertex = next_right_vertex;
                right_edge = next_right_edge;
                continue;
            }
            break;
        }
        return {left_vertex, left_edge, right_vertex, right_edge};
    }

    std::pair<int, int> extreme_vertex(int vertex, int edge, bool minimum) const {
        std::pair<int, int> result = {vertex, edge};
        int current_vertex = vertex;
        int current_edge = edge;
        do {
            std::tie(current_vertex, current_edge) = go_next(current_edge);
            std::pair<int, int> candidate = {current_vertex, current_edge};
            if ((minimum && candidate < result) || (!minimum && result < candidate)) {
                result = candidate;
            }
        } while (current_edge != edge);
        return result;
    }

    bool inside_circumcircle(int a, int b, int c, int d) const {
        return inside_circumcircle(points[a], points[b], points[c], points[d]);
    }

    std::pair<int, int> merge_triangulations(
        int left_vertex,
        int left_edge,
        int right_vertex,
        int right_edge
    ) {
        std::tie(left_vertex, left_edge) = extreme_vertex(left_vertex, left_edge, false);
        std::tie(right_vertex, right_edge) = extreme_vertex(right_vertex, right_edge, true);

        auto [lower_left, lower_left_edge, lower_right, lower_right_edge]
            = lower_tangent(left_vertex, left_edge, right_vertex, right_edge);
        auto [upper_right, upper_right_edge, upper_left, upper_left_edge]
            = lower_tangent(right_vertex, right_edge, left_vertex, left_edge);
        lower_right_edge = edges[lower_right_edge].cw;
        upper_right_edge = edges[upper_right_edge].cw;

        auto [base, reverse_base] = add_edge(lower_left, lower_right);
        insert_cw_after(base, lower_left_edge);
        insert_ccw_after(reverse_base, lower_right_edge);
        if (lower_left == upper_left) upper_left_edge = base;
        if (lower_right == upper_right) upper_right_edge = reverse_base;

        int left = lower_left;
        int left_candidate = lower_left_edge;
        int right = lower_right;
        int right_candidate = lower_right_edge;
        while (left != upper_left || right != upper_right) {
            int next_left = edges[left_candidate].to;
            int next_right = edges[right_candidate].to;
            int next_left_candidate = edges[left_candidate].ccw;
            int next_right_candidate = edges[right_candidate].cw;

            if (left_candidate != upper_left_edge && next_left_candidate != base) {
                int second_left = edges[next_left_candidate].to;
                if (inside_circumcircle(left, right, next_left, second_left)) {
                    erase_edge(left_candidate);
                    left_candidate = next_left_candidate;
                    continue;
                }
            }

            if (right_candidate != upper_right_edge && next_right_candidate != reverse_base) {
                int second_right = edges[next_right_candidate].to;
                if (inside_circumcircle(next_right, left, right, second_right)) {
                    erase_edge(right_candidate);
                    right_candidate = next_right_candidate;
                    continue;
                }
            }

            bool choose_left = right_candidate == upper_right_edge;
            if (left_candidate != upper_left_edge && right_candidate != upper_right_edge) {
                if (orientation(left, right, next_right) < 0) {
                    choose_left = true;
                } else if (orientation(next_left, left, right) < 0) {
                    choose_left = false;
                } else {
                    choose_left = inside_circumcircle(left, right, next_right, next_left);
                }
            }

            if (choose_left) {
                next_left_candidate = edges[edges[left_candidate].reverse].ccw;
                auto [new_base, new_reverse_base] = add_edge(next_left, right);
                insert_cw_after(new_base, next_left_candidate);
                insert_ccw_after(new_reverse_base, right_candidate);
                left_candidate = next_left_candidate;
                left = next_left;
            } else {
                next_right_candidate = edges[edges[right_candidate].reverse].cw;
                auto [new_reverse_base, new_base] = add_edge(next_right, left);
                insert_ccw_after(new_reverse_base, next_right_candidate);
                insert_cw_after(new_base, left_candidate);
                right_candidate = next_right_candidate;
                right = next_right;
            }
        }
        return {lower_left, base};
    }

    std::pair<int, int> solve_range(int left, int right) {
        if (right - left == 2) {
            auto [forward, backward] = add_edge(left, left + 1);
            (void)backward;
            return {left, forward};
        }
        if (right - left == 3) {
            int middle = left + 1;
            int last = left + 2;
            auto [first_middle, middle_first] = add_edge(left, middle);
            auto [middle_last, last_middle] = add_edge(middle, last);
            int direction = orientation(left, middle, last);
            if (direction == 0) {
                insert_ccw_after(middle_first, middle_last);
                return {left, first_middle};
            }

            auto [first_last, last_first] = add_edge(left, last);
            if (direction > 0) {
                insert_cw_after(first_middle, first_last);
                insert_cw_after(middle_last, middle_first);
                insert_cw_after(last_first, last_middle);
                return {left, first_middle};
            }
            insert_ccw_after(first_middle, first_last);
            insert_ccw_after(middle_last, middle_first);
            insert_ccw_after(last_first, last_middle);
            return {middle, middle_first};
        }

        int middle = (left + right) / 2;
        auto [left_vertex, left_edge] = solve_range(left, middle);
        auto [right_vertex, right_edge] = solve_range(middle, right);
        return merge_triangulations(left_vertex, left_edge, right_vertex, right_edge);
    }

    void solve() {
        int size = int(points.size());
        if (size <= 1) return;

        std::vector<int> order(size);
        for (int i = 0; i < size; i++) order[i] = i;
        std::stable_sort(order.begin(), order.end(), [&](int left, int right) {
            if (points[left].x != points[right].x) {
                return points[left].x < points[right].x;
            }
            return points[left].y < points[right].y;
        });

        std::vector<InternalPoint> original_points = points;
        duplicate_representative.assign(size, 0);
        int unique_size = 0;
        for (int i = 0; i < size; i++) {
            int vertex = order[i];
            if (i == 0 || !(original_points[order[unique_size - 1]] == original_points[vertex])) {
                order[unique_size] = vertex;
                points[unique_size] = original_points[vertex];
                unique_size++;
                duplicate_representative[vertex] = vertex;
            } else {
                duplicate_representative[vertex] = order[unique_size - 1];
            }
        }

        if (unique_size >= 2) solve_range(0, unique_size);
        points.swap(original_points);
        for (auto& edge : edges) edge.to = order[edge.to];
    }

   public:
    explicit EuclideanDelaunay(const std::vector<Point<T>>& input_points) {
        assert(input_points.size() <= std::size_t(std::numeric_limits<int>::max()));
        points.reserve(input_points.size());
        edges.reserve(std::size_t(6) * input_points.size());
        for (const auto& point : input_points) {
            points.push_back(InternalPoint{W(point.x), W(point.y)});
        }
        solve();
    }

    bool has_duplicates() const {
        for (
            int vertex = 0;
            vertex < int(duplicate_representative.size());
            ++vertex
        ) {
            if (duplicate_representative[vertex] != vertex) return true;
        }
        return false;
    }

    std::vector<std::pair<int, int>> get_edges() const {
        std::vector<std::pair<int, int>> result;
        result.reserve(edges.size() / 2 + duplicate_representative.size());
        for (int edge = 0; edge < int(edges.size()); edge++) {
            if (!edges[edge].enabled) continue;
            int reverse = edges[edge].reverse;
            if (edge < reverse) continue;
            result.emplace_back(edges[edge].to, edges[reverse].to);
        }
        for (int vertex = 0; vertex < int(duplicate_representative.size()); vertex++) {
            if (duplicate_representative[vertex] != vertex) {
                result.emplace_back(vertex, duplicate_representative[vertex]);
            }
        }
        return result;
    }
};

}  // namespace detail

// Returns O(n) Delaunay edges containing a Euclidean minimum spanning tree.
template <ExactCoordinate T>
std::vector<EuclideanMstEdge<wide_type<T>>> euclidean_mst_edges(
    const std::vector<Point<T>>& points
) {
    using W = wide_type<T>;
    auto delaunay_edges = detail::EuclideanDelaunay<T>(points).get_edges();
    std::vector<EuclideanMstEdge<W>> result;
    result.reserve(delaunay_edges.size());
    for (auto [from, to] : delaunay_edges) {
        result.push_back(EuclideanMstEdge<W>{from, to, distance2(points[from], points[to])});
    }
    return result;
}

// Returns a Euclidean minimum spanning tree.
template <ExactCoordinate T>
EuclideanMst<wide_type<T>> euclidean_mst(const std::vector<Point<T>>& points) {
    using W = wide_type<T>;
    auto candidates = euclidean_mst_edges(points);
    std::sort(candidates.begin(), candidates.end(), [](const auto& left, const auto& right) {
        if (left.squared_distance != right.squared_distance) {
            return left.squared_distance < right.squared_distance;
        }
        if (left.from != right.from) return left.from < right.from;
        return left.to < right.to;
    });

    m1une::ds::Dsu dsu(int(points.size()));
    EuclideanMst<W> result;
    result.cost = 0;
    result.edges.reserve(points.empty() ? 0 : points.size() - 1);
    for (const auto& edge : candidates) {
        if (dsu.same(edge.from, edge.to)) continue;
        dsu.merge(edge.from, edge.to);
        result.cost += std::sqrt(static_cast<long double>(edge.squared_distance));
        result.edges.push_back(edge);
        if (result.edges.size() + 1 == points.size()) break;
    }
    assert(points.empty() || result.edges.size() + 1 == points.size());
    return result;
}

}  // namespace geometry
}  // namespace m1une


#line 12 "geometry/delaunay_triangulation.hpp"

namespace m1une {
namespace geometry {

struct DelaunayTriangulation {
    std::vector<std::pair<int, int>> edges;
    std::vector<std::array<int, 3>> triangles;
};

namespace delaunay_triangulation_detail {

template <ExactCoordinate T>
int direction_half(
    const Point<T>& origin,
    const Point<T>& destination
) {
    using W = wide_type<T>;
    W x = W(destination.x) - W(origin.x);
    W y = W(destination.y) - W(origin.y);
    return y > 0 || (y == 0 && x >= 0) ? 0 : 1;
}

template <ExactCoordinate T>
bool direction_less(
    const std::vector<Point<T>>& points,
    int origin,
    int first,
    int second
) {
    int first_half = direction_half(points[origin], points[first]);
    int second_half = direction_half(points[origin], points[second]);
    if (first_half != second_half) return first_half < second_half;

    using W = wide_type<T>;
    W first_x = W(points[first].x) - W(points[origin].x);
    W first_y = W(points[first].y) - W(points[origin].y);
    W second_x = W(points[second].x) - W(points[origin].x);
    W second_y = W(points[second].y) - W(points[origin].y);
    W product = first_x * second_y - first_y * second_x;
    if (product != 0) return product > 0;

    W first_norm = first_x * first_x + first_y * first_y;
    W second_norm = second_x * second_x + second_y * second_y;
    if (first_norm != second_norm) return first_norm < second_norm;
    return first < second;
}

inline void rotate_minimum_first(std::array<int, 3>& triangle) {
    int minimum = int(std::min_element(triangle.begin(), triangle.end()) -
                      triangle.begin());
    std::rotate(
        triangle.begin(),
        triangle.begin() + minimum,
        triangle.end()
    );
}

}  // namespace delaunay_triangulation_detail

// Constructs one Delaunay triangulation of distinct exact-coordinate points.
template <ExactCoordinate T>
DelaunayTriangulation delaunay_triangulation(
    const std::vector<Point<T>>& points
) {
    namespace detail = delaunay_triangulation_detail;

    DelaunayTriangulation result;
    geometry::detail::EuclideanDelaunay<T> builder(points);
    assert(!builder.has_duplicates());
    result.edges = builder.get_edges();
    for (auto& [first, second] : result.edges) {
        if (first > second) std::swap(first, second);
    }
    std::sort(result.edges.begin(), result.edges.end());
    result.edges.erase(
        std::unique(result.edges.begin(), result.edges.end()),
        result.edges.end()
    );

    std::vector<std::vector<int>> neighbors(points.size());
    for (auto [first, second] : result.edges) {
        neighbors[first].push_back(second);
        neighbors[second].push_back(first);
    }
    for (int point = 0; point < int(points.size()); ++point) {
        std::sort(
            neighbors[point].begin(),
            neighbors[point].end(),
            [&](int first, int second) {
                return detail::direction_less(points, point, first, second);
            }
        );
    }

    auto has_edge = [&](int first, int second) {
        if (first > second) std::swap(first, second);
        return std::binary_search(
            result.edges.begin(),
            result.edges.end(),
            std::pair(first, second)
        );
    };
    result.triangles.reserve(result.edges.size());
    for (int point = 0; point < int(points.size()); ++point) {
        int degree = int(neighbors[point].size());
        for (int index = 0; index < degree; ++index) {
            int first = neighbors[point][index];
            int second = neighbors[point][(index + 1) % degree];
            if (orientation(points[point], points[first], points[second]) <= 0) {
                continue;
            }
            if (!has_edge(first, second)) continue;
            std::array<int, 3> triangle{point, first, second};
            detail::rotate_minimum_first(triangle);
            result.triangles.push_back(triangle);
        }
    }
    std::sort(result.triangles.begin(), result.triangles.end());
    result.triangles.erase(
        std::unique(result.triangles.begin(), result.triangles.end()),
        result.triangles.end()
    );
    return result;
}

}  // namespace geometry
}  // namespace m1une

Back to top page