Voronoi Diagram
(geometry/voronoi_diagram.hpp)
- View this file on GitHub
- Last update: 2026-10-05 22:23:07+09:00
- Include:
#include "geometry/voronoi_diagram.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.
Returned lengths or constructed coordinates still use long double.
voronoi_diagram constructs the ordinary Euclidean Voronoi diagram of a set
of distinct exact-coordinate sites. The diagram is the geometric dual of an exact
Delaunay triangulation, so it contains finite segments, unbounded rays, and,
when every site is collinear, full lines.
Combinatorial decisions use exact predicates. Finite vertices and
parametric edge geometry use long double because circumcenters need not be
integral.
Types
enum class VoronoiEdgeKind {
Segment,
Ray,
Line,
};
struct VoronoiEdge {
VoronoiEdgeKind kind;
int first_site;
int second_site;
int first_vertex;
int second_vertex;
Point<long double> point;
Point<long double> direction;
};
struct VoronoiDiagram {
std::vector<Point<long double>> vertices;
std::vector<VoronoiEdge> edges;
std::vector<std::vector<int>> cell_edges;
};
Each edge is the shared one-dimensional boundary of the cells belonging to
first_site and second_site, with first_site < second_site. Its points are
described by
| Kind | Parameter range | Vertex indices | Direction |
|---|---|---|---|
Segment |
$0 \leq t \leq 1$ | Both indices are valid; point is vertices[first_vertex]. |
The displacement from first_vertex to second_vertex. |
Ray |
$0 \leq t$ |
first_vertex is the finite endpoint and second_vertex == -1. |
A unit vector pointing toward infinity. |
Line |
$-\infty < t < \infty$ | Both indices are -1. |
A unit vector along the perpendicular bisector. |
cell_edges[site] contains the indices of all one-dimensional boundary edges
of that site’s cell. They are ordered counterclockwise by the direction from
site to the neighboring site. Two cells that meet only at one Voronoi vertex
do not share an edge and therefore do not appear in each other’s lists.
Function
The exact signature is:
template <ExactCoordinate T>
VoronoiDiagram voronoi_diagram(
const std::vector<Point<T>>& sites
);
| Function | Description | Complexity |
|---|---|---|
voronoi_diagram(sites) |
Returns all finite vertices, parametric edges, and per-site boundary-edge lists. | $O(N\log N)$ time and $O(N)$ memory |
Site indices are their zero-based positions in sites. All sites must be
pairwise distinct. Empty and one-site inputs have no vertices or edges. Two
sites produce one Line. A larger collinear input produces one line between
each consecutive pair of sites.
Cocircular Delaunay triangles are merged into one Voronoi vertex. Consequently, the result does not expose artificial zero-length edges introduced by a choice of Delaunay triangulation.
The degree-four incircle expressions must fit in signed 128-bit arithmetic. This is satisfied, for example, by signed 32-bit input coordinates. As with the other floating-point geometry constructions, callers should use an appropriate tolerance when comparing returned coordinates.
Example
#include "geometry/voronoi_diagram.hpp"
#include <iostream>
#include <vector>
int main() {
using Point = m1une::geometry::Point<long long>;
std::vector<Point> sites;
sites.emplace_back(0, 0);
sites.emplace_back(6, 0);
sites.emplace_back(0, 8);
auto diagram = m1une::geometry::voronoi_diagram(sites);
std::cout << diagram.vertices.size() << '\n'; // 1
std::cout << diagram.edges.size() << '\n'; // 3 rays
for (int edge_index : diagram.cell_edges[0]) {
const auto& edge = diagram.edges[edge_index];
std::cout << edge.first_site << ' ' << edge.second_site << '\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
Verified with
- verify/geometry/centroid.test.cpp
- verify/geometry/geometry_algorithms.test.cpp
- verify/geometry/rational.test.cpp
- verify/geometry/voronoi_diagram.test.cpp
Code
#ifndef M1UNE_GEOMETRY_VORONOI_DIAGRAM_HPP
#define M1UNE_GEOMETRY_VORONOI_DIAGRAM_HPP 1
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <limits>
#include <numeric>
#include <utility>
#include <vector>
#include "euclidean_mst.hpp"
namespace m1une {
namespace geometry {
enum class VoronoiEdgeKind {
Segment,
Ray,
Line,
};
struct VoronoiEdge {
VoronoiEdgeKind kind;
int first_site;
int second_site;
int first_vertex;
int second_vertex;
Point<long double> point;
Point<long double> direction;
};
struct VoronoiDiagram {
std::vector<Point<long double>> vertices;
std::vector<VoronoiEdge> edges;
std::vector<std::vector<int>> cell_edges;
};
namespace voronoi_diagram_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>>& sites,
int origin,
int first,
int second
) {
int first_half = direction_half(sites[origin], sites[first]);
int second_half = direction_half(sites[origin], sites[second]);
if (first_half != second_half) return first_half < second_half;
using W = wide_type<T>;
W first_x = W(sites[first].x) - W(sites[origin].x);
W first_y = W(sites[first].y) - W(sites[origin].y);
W second_x = W(sites[second].x) - W(sites[origin].x);
W second_y = W(sites[second].y) - W(sites[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;
}
template <ExactCoordinate T>
bool cocircular(
const Point<T>& first,
const Point<T>& second,
const Point<T>& third,
const Point<T>& fourth
) {
using W = wide_type<T>;
W ax = W(first.x) - W(fourth.x);
W ay = W(first.y) - W(fourth.y);
W bx = W(second.x) - W(fourth.x);
W by = W(second.y) - W(fourth.y);
W cx = W(third.x) - W(fourth.x);
W cy = W(third.y) - W(fourth.y);
W a_norm = ax * ax + ay * ay;
W b_norm = bx * bx + by * by;
W c_norm = cx * cx + cy * cy;
W determinant =
(bx * cy - by * cx) * a_norm +
(cx * ay - cy * ax) * b_norm +
(ax * by - ay * bx) * c_norm;
return determinant == 0;
}
template <ExactCoordinate T>
Point<long double> circumcenter(
const Point<T>& first,
const Point<T>& second,
const Point<T>& third
) {
long double ax = static_cast<long double>(first.x);
long double ay = static_cast<long double>(first.y);
long double bx = static_cast<long double>(second.x);
long double by = static_cast<long double>(second.y);
long double cx = static_cast<long double>(third.x);
long double cy = static_cast<long double>(third.y);
long double denominator = 2 * (
ax * (by - cy) +
bx * (cy - ay) +
cx * (ay - by)
);
assert(denominator != 0);
long double first_norm = ax * ax + ay * ay;
long double second_norm = bx * bx + by * by;
long double third_norm = cx * cx + cy * cy;
return Point<long double>(
(first_norm * (by - cy) +
second_norm * (cy - ay) +
third_norm * (ay - by)) /
denominator,
(first_norm * (cx - bx) +
second_norm * (ax - cx) +
third_norm * (bx - ax)) /
denominator
);
}
inline Point<long double> unit(Point<long double> direction) {
long double length = norm(direction);
assert(length != 0);
return direction / length;
}
inline int other_site(const VoronoiEdge& edge, int site) {
assert(edge.first_site == site || edge.second_site == site);
return edge.first_site == site ? edge.second_site : edge.first_site;
}
} // namespace voronoi_diagram_detail
// Constructs the ordinary Euclidean Voronoi diagram of distinct exact-coordinate sites.
template <ExactCoordinate T>
VoronoiDiagram voronoi_diagram(const std::vector<Point<T>>& sites) {
namespace detail = voronoi_diagram_detail;
assert(sites.size() <= std::size_t(std::numeric_limits<int>::max()));
const int size = int(sites.size());
std::vector<int> site_order(size);
std::iota(site_order.begin(), site_order.end(), 0);
std::sort(site_order.begin(), site_order.end(), [&](int first, int second) {
return sites[first] < sites[second];
});
for (int index = 1; index < size; ++index) {
assert(sites[site_order[index - 1]] != sites[site_order[index]]);
}
std::vector<std::pair<int, int>> delaunay_edges =
geometry::detail::EuclideanDelaunay<T>(sites).get_edges();
for (auto& [first, second] : delaunay_edges) {
if (first > second) std::swap(first, second);
}
std::sort(delaunay_edges.begin(), delaunay_edges.end());
delaunay_edges.erase(
std::unique(delaunay_edges.begin(), delaunay_edges.end()),
delaunay_edges.end()
);
auto find_edge_index = [&](int first, int second) {
if (first > second) std::swap(first, second);
auto iterator = std::lower_bound(
delaunay_edges.begin(),
delaunay_edges.end(),
std::pair(first, second)
);
if (
iterator == delaunay_edges.end() ||
*iterator != std::pair(first, second)
) {
return -1;
}
return int(iterator - delaunay_edges.begin());
};
std::vector<std::vector<int>> neighbors(size);
for (int index = 0; index < int(delaunay_edges.size()); ++index) {
auto [first, second] = delaunay_edges[index];
neighbors[first].push_back(second);
neighbors[second].push_back(first);
}
for (int site = 0; site < size; ++site) {
std::sort(
neighbors[site].begin(),
neighbors[site].end(),
[&](int first, int second) {
return detail::direction_less(sites, site, first, second);
}
);
}
std::vector<std::array<int, 3>> triangles;
for (int site = 0; site < size; ++site) {
int degree = int(neighbors[site].size());
for (int index = 0; index < degree; ++index) {
int first = neighbors[site][index];
int second = neighbors[site][(index + 1) % degree];
if (orientation(sites[site], sites[first], sites[second]) <= 0) {
continue;
}
if (find_edge_index(first, second) == -1) continue;
std::array<int, 3> triangle{site, first, second};
std::sort(triangle.begin(), triangle.end());
triangles.push_back(triangle);
}
}
std::sort(triangles.begin(), triangles.end());
triangles.erase(
std::unique(triangles.begin(), triangles.end()),
triangles.end()
);
for (auto& triangle : triangles) {
if (orientation(
sites[triangle[0]],
sites[triangle[1]],
sites[triangle[2]]
) < 0) {
std::swap(triangle[1], triangle[2]);
}
}
std::vector<std::array<int, 2>> incident_triangles(
delaunay_edges.size(),
std::array<int, 2>{-1, -1}
);
std::vector<int> incident_count(delaunay_edges.size(), 0);
for (int triangle = 0; triangle < int(triangles.size()); ++triangle) {
for (int side = 0; side < 3; ++side) {
int first = triangles[triangle][side];
int second = triangles[triangle][(side + 1) % 3];
int edge = find_edge_index(first, second);
assert(edge != -1);
assert(incident_count[edge] < 2);
incident_triangles[edge][incident_count[edge]++] = triangle;
}
}
std::vector<int> parent(triangles.size());
std::vector<int> component_size(triangles.size(), 1);
std::iota(parent.begin(), parent.end(), 0);
auto find_root = [&](auto&& self, int vertex) -> int {
if (parent[vertex] == vertex) return vertex;
return parent[vertex] = self(self, parent[vertex]);
};
auto merge = [&](int first, int second) {
first = find_root(find_root, first);
second = find_root(find_root, second);
if (first == second) return;
if (component_size[first] < component_size[second]) {
std::swap(first, second);
}
parent[second] = first;
component_size[first] += component_size[second];
};
for (int edge = 0; edge < int(delaunay_edges.size()); ++edge) {
if (incident_count[edge] != 2) continue;
int first_triangle = incident_triangles[edge][0];
int second_triangle = incident_triangles[edge][1];
const auto& first = triangles[first_triangle];
const auto& second = triangles[second_triangle];
int fourth = second[0];
if (fourth == first[0] || fourth == first[1] || fourth == first[2]) {
fourth = second[1];
}
if (fourth == first[0] || fourth == first[1] || fourth == first[2]) {
fourth = second[2];
}
assert(
fourth != first[0] &&
fourth != first[1] &&
fourth != first[2]
);
if (detail::cocircular(
sites[first[0]],
sites[first[1]],
sites[first[2]],
sites[fourth]
)) {
merge(first_triangle, second_triangle);
}
}
VoronoiDiagram result;
result.cell_edges.resize(size);
std::vector<int> root_vertex(triangles.size(), -1);
std::vector<int> triangle_vertex(triangles.size(), -1);
for (int triangle = 0; triangle < int(triangles.size()); ++triangle) {
int root = find_root(find_root, triangle);
if (root_vertex[root] == -1) {
const auto& sites_on_circle = triangles[triangle];
root_vertex[root] = int(result.vertices.size());
result.vertices.push_back(detail::circumcenter(
sites[sites_on_circle[0]],
sites[sites_on_circle[1]],
sites[sites_on_circle[2]]
));
}
triangle_vertex[triangle] = root_vertex[root];
}
result.edges.reserve(delaunay_edges.size());
for (int edge = 0; edge < int(delaunay_edges.size()); ++edge) {
auto [first_site, second_site] = delaunay_edges[edge];
VoronoiEdge voronoi_edge;
voronoi_edge.first_site = first_site;
voronoi_edge.second_site = second_site;
voronoi_edge.first_vertex = -1;
voronoi_edge.second_vertex = -1;
if (incident_count[edge] == 2) {
int first_vertex =
triangle_vertex[incident_triangles[edge][0]];
int second_vertex =
triangle_vertex[incident_triangles[edge][1]];
if (first_vertex == second_vertex) continue;
if (first_vertex > second_vertex) {
std::swap(first_vertex, second_vertex);
}
voronoi_edge.kind = VoronoiEdgeKind::Segment;
voronoi_edge.first_vertex = first_vertex;
voronoi_edge.second_vertex = second_vertex;
voronoi_edge.point = result.vertices[first_vertex];
voronoi_edge.direction =
result.vertices[second_vertex] - result.vertices[first_vertex];
} else if (incident_count[edge] == 1) {
int triangle = incident_triangles[edge][0];
int third_site = triangles[triangle][0];
if (third_site == first_site || third_site == second_site) {
third_site = triangles[triangle][1];
}
if (third_site == first_site || third_site == second_site) {
third_site = triangles[triangle][2];
}
assert(third_site != first_site && third_site != second_site);
Point<long double> first(sites[first_site]);
Point<long double> second(sites[second_site]);
Point<long double> edge_direction = second - first;
Point<long double> outward;
if (orientation(
sites[first_site],
sites[second_site],
sites[third_site]
) > 0) {
outward = Point<long double>(
edge_direction.y,
-edge_direction.x
);
} else {
outward = Point<long double>(
-edge_direction.y,
edge_direction.x
);
}
voronoi_edge.kind = VoronoiEdgeKind::Ray;
voronoi_edge.first_vertex = triangle_vertex[triangle];
voronoi_edge.point = result.vertices[voronoi_edge.first_vertex];
voronoi_edge.direction = detail::unit(outward);
} else {
assert(incident_count[edge] == 0);
Point<long double> first(sites[first_site]);
Point<long double> second(sites[second_site]);
Point<long double> edge_direction = second - first;
voronoi_edge.kind = VoronoiEdgeKind::Line;
voronoi_edge.point = (first + second) / 2.0L;
voronoi_edge.direction = detail::unit(Point<long double>(
edge_direction.y,
-edge_direction.x
));
}
int voronoi_edge_index = int(result.edges.size());
result.edges.push_back(voronoi_edge);
result.cell_edges[first_site].push_back(voronoi_edge_index);
result.cell_edges[second_site].push_back(voronoi_edge_index);
}
for (int site = 0; site < size; ++site) {
std::sort(
result.cell_edges[site].begin(),
result.cell_edges[site].end(),
[&](int first_edge, int second_edge) {
int first_other =
detail::other_site(result.edges[first_edge], site);
int second_other =
detail::other_site(result.edges[second_edge], site);
return detail::direction_less(
sites,
site,
first_other,
second_other
);
}
);
}
return result;
}
} // namespace geometry
} // namespace m1une
#endif // M1UNE_GEOMETRY_VORONOI_DIAGRAM_HPP#line 1 "geometry/voronoi_diagram.hpp"
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <limits>
#include <numeric>
#include <utility>
#include <vector>
#line 1 "geometry/euclidean_mst.hpp"
#line 10 "geometry/euclidean_mst.hpp"
#include <tuple>
#line 13 "geometry/euclidean_mst.hpp"
#line 1 "ds/dsu/dsu.hpp"
#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 16 "geometry/voronoi_diagram.hpp"
namespace m1une {
namespace geometry {
enum class VoronoiEdgeKind {
Segment,
Ray,
Line,
};
struct VoronoiEdge {
VoronoiEdgeKind kind;
int first_site;
int second_site;
int first_vertex;
int second_vertex;
Point<long double> point;
Point<long double> direction;
};
struct VoronoiDiagram {
std::vector<Point<long double>> vertices;
std::vector<VoronoiEdge> edges;
std::vector<std::vector<int>> cell_edges;
};
namespace voronoi_diagram_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>>& sites,
int origin,
int first,
int second
) {
int first_half = direction_half(sites[origin], sites[first]);
int second_half = direction_half(sites[origin], sites[second]);
if (first_half != second_half) return first_half < second_half;
using W = wide_type<T>;
W first_x = W(sites[first].x) - W(sites[origin].x);
W first_y = W(sites[first].y) - W(sites[origin].y);
W second_x = W(sites[second].x) - W(sites[origin].x);
W second_y = W(sites[second].y) - W(sites[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;
}
template <ExactCoordinate T>
bool cocircular(
const Point<T>& first,
const Point<T>& second,
const Point<T>& third,
const Point<T>& fourth
) {
using W = wide_type<T>;
W ax = W(first.x) - W(fourth.x);
W ay = W(first.y) - W(fourth.y);
W bx = W(second.x) - W(fourth.x);
W by = W(second.y) - W(fourth.y);
W cx = W(third.x) - W(fourth.x);
W cy = W(third.y) - W(fourth.y);
W a_norm = ax * ax + ay * ay;
W b_norm = bx * bx + by * by;
W c_norm = cx * cx + cy * cy;
W determinant =
(bx * cy - by * cx) * a_norm +
(cx * ay - cy * ax) * b_norm +
(ax * by - ay * bx) * c_norm;
return determinant == 0;
}
template <ExactCoordinate T>
Point<long double> circumcenter(
const Point<T>& first,
const Point<T>& second,
const Point<T>& third
) {
long double ax = static_cast<long double>(first.x);
long double ay = static_cast<long double>(first.y);
long double bx = static_cast<long double>(second.x);
long double by = static_cast<long double>(second.y);
long double cx = static_cast<long double>(third.x);
long double cy = static_cast<long double>(third.y);
long double denominator = 2 * (
ax * (by - cy) +
bx * (cy - ay) +
cx * (ay - by)
);
assert(denominator != 0);
long double first_norm = ax * ax + ay * ay;
long double second_norm = bx * bx + by * by;
long double third_norm = cx * cx + cy * cy;
return Point<long double>(
(first_norm * (by - cy) +
second_norm * (cy - ay) +
third_norm * (ay - by)) /
denominator,
(first_norm * (cx - bx) +
second_norm * (ax - cx) +
third_norm * (bx - ax)) /
denominator
);
}
inline Point<long double> unit(Point<long double> direction) {
long double length = norm(direction);
assert(length != 0);
return direction / length;
}
inline int other_site(const VoronoiEdge& edge, int site) {
assert(edge.first_site == site || edge.second_site == site);
return edge.first_site == site ? edge.second_site : edge.first_site;
}
} // namespace voronoi_diagram_detail
// Constructs the ordinary Euclidean Voronoi diagram of distinct exact-coordinate sites.
template <ExactCoordinate T>
VoronoiDiagram voronoi_diagram(const std::vector<Point<T>>& sites) {
namespace detail = voronoi_diagram_detail;
assert(sites.size() <= std::size_t(std::numeric_limits<int>::max()));
const int size = int(sites.size());
std::vector<int> site_order(size);
std::iota(site_order.begin(), site_order.end(), 0);
std::sort(site_order.begin(), site_order.end(), [&](int first, int second) {
return sites[first] < sites[second];
});
for (int index = 1; index < size; ++index) {
assert(sites[site_order[index - 1]] != sites[site_order[index]]);
}
std::vector<std::pair<int, int>> delaunay_edges =
geometry::detail::EuclideanDelaunay<T>(sites).get_edges();
for (auto& [first, second] : delaunay_edges) {
if (first > second) std::swap(first, second);
}
std::sort(delaunay_edges.begin(), delaunay_edges.end());
delaunay_edges.erase(
std::unique(delaunay_edges.begin(), delaunay_edges.end()),
delaunay_edges.end()
);
auto find_edge_index = [&](int first, int second) {
if (first > second) std::swap(first, second);
auto iterator = std::lower_bound(
delaunay_edges.begin(),
delaunay_edges.end(),
std::pair(first, second)
);
if (
iterator == delaunay_edges.end() ||
*iterator != std::pair(first, second)
) {
return -1;
}
return int(iterator - delaunay_edges.begin());
};
std::vector<std::vector<int>> neighbors(size);
for (int index = 0; index < int(delaunay_edges.size()); ++index) {
auto [first, second] = delaunay_edges[index];
neighbors[first].push_back(second);
neighbors[second].push_back(first);
}
for (int site = 0; site < size; ++site) {
std::sort(
neighbors[site].begin(),
neighbors[site].end(),
[&](int first, int second) {
return detail::direction_less(sites, site, first, second);
}
);
}
std::vector<std::array<int, 3>> triangles;
for (int site = 0; site < size; ++site) {
int degree = int(neighbors[site].size());
for (int index = 0; index < degree; ++index) {
int first = neighbors[site][index];
int second = neighbors[site][(index + 1) % degree];
if (orientation(sites[site], sites[first], sites[second]) <= 0) {
continue;
}
if (find_edge_index(first, second) == -1) continue;
std::array<int, 3> triangle{site, first, second};
std::sort(triangle.begin(), triangle.end());
triangles.push_back(triangle);
}
}
std::sort(triangles.begin(), triangles.end());
triangles.erase(
std::unique(triangles.begin(), triangles.end()),
triangles.end()
);
for (auto& triangle : triangles) {
if (orientation(
sites[triangle[0]],
sites[triangle[1]],
sites[triangle[2]]
) < 0) {
std::swap(triangle[1], triangle[2]);
}
}
std::vector<std::array<int, 2>> incident_triangles(
delaunay_edges.size(),
std::array<int, 2>{-1, -1}
);
std::vector<int> incident_count(delaunay_edges.size(), 0);
for (int triangle = 0; triangle < int(triangles.size()); ++triangle) {
for (int side = 0; side < 3; ++side) {
int first = triangles[triangle][side];
int second = triangles[triangle][(side + 1) % 3];
int edge = find_edge_index(first, second);
assert(edge != -1);
assert(incident_count[edge] < 2);
incident_triangles[edge][incident_count[edge]++] = triangle;
}
}
std::vector<int> parent(triangles.size());
std::vector<int> component_size(triangles.size(), 1);
std::iota(parent.begin(), parent.end(), 0);
auto find_root = [&](auto&& self, int vertex) -> int {
if (parent[vertex] == vertex) return vertex;
return parent[vertex] = self(self, parent[vertex]);
};
auto merge = [&](int first, int second) {
first = find_root(find_root, first);
second = find_root(find_root, second);
if (first == second) return;
if (component_size[first] < component_size[second]) {
std::swap(first, second);
}
parent[second] = first;
component_size[first] += component_size[second];
};
for (int edge = 0; edge < int(delaunay_edges.size()); ++edge) {
if (incident_count[edge] != 2) continue;
int first_triangle = incident_triangles[edge][0];
int second_triangle = incident_triangles[edge][1];
const auto& first = triangles[first_triangle];
const auto& second = triangles[second_triangle];
int fourth = second[0];
if (fourth == first[0] || fourth == first[1] || fourth == first[2]) {
fourth = second[1];
}
if (fourth == first[0] || fourth == first[1] || fourth == first[2]) {
fourth = second[2];
}
assert(
fourth != first[0] &&
fourth != first[1] &&
fourth != first[2]
);
if (detail::cocircular(
sites[first[0]],
sites[first[1]],
sites[first[2]],
sites[fourth]
)) {
merge(first_triangle, second_triangle);
}
}
VoronoiDiagram result;
result.cell_edges.resize(size);
std::vector<int> root_vertex(triangles.size(), -1);
std::vector<int> triangle_vertex(triangles.size(), -1);
for (int triangle = 0; triangle < int(triangles.size()); ++triangle) {
int root = find_root(find_root, triangle);
if (root_vertex[root] == -1) {
const auto& sites_on_circle = triangles[triangle];
root_vertex[root] = int(result.vertices.size());
result.vertices.push_back(detail::circumcenter(
sites[sites_on_circle[0]],
sites[sites_on_circle[1]],
sites[sites_on_circle[2]]
));
}
triangle_vertex[triangle] = root_vertex[root];
}
result.edges.reserve(delaunay_edges.size());
for (int edge = 0; edge < int(delaunay_edges.size()); ++edge) {
auto [first_site, second_site] = delaunay_edges[edge];
VoronoiEdge voronoi_edge;
voronoi_edge.first_site = first_site;
voronoi_edge.second_site = second_site;
voronoi_edge.first_vertex = -1;
voronoi_edge.second_vertex = -1;
if (incident_count[edge] == 2) {
int first_vertex =
triangle_vertex[incident_triangles[edge][0]];
int second_vertex =
triangle_vertex[incident_triangles[edge][1]];
if (first_vertex == second_vertex) continue;
if (first_vertex > second_vertex) {
std::swap(first_vertex, second_vertex);
}
voronoi_edge.kind = VoronoiEdgeKind::Segment;
voronoi_edge.first_vertex = first_vertex;
voronoi_edge.second_vertex = second_vertex;
voronoi_edge.point = result.vertices[first_vertex];
voronoi_edge.direction =
result.vertices[second_vertex] - result.vertices[first_vertex];
} else if (incident_count[edge] == 1) {
int triangle = incident_triangles[edge][0];
int third_site = triangles[triangle][0];
if (third_site == first_site || third_site == second_site) {
third_site = triangles[triangle][1];
}
if (third_site == first_site || third_site == second_site) {
third_site = triangles[triangle][2];
}
assert(third_site != first_site && third_site != second_site);
Point<long double> first(sites[first_site]);
Point<long double> second(sites[second_site]);
Point<long double> edge_direction = second - first;
Point<long double> outward;
if (orientation(
sites[first_site],
sites[second_site],
sites[third_site]
) > 0) {
outward = Point<long double>(
edge_direction.y,
-edge_direction.x
);
} else {
outward = Point<long double>(
-edge_direction.y,
edge_direction.x
);
}
voronoi_edge.kind = VoronoiEdgeKind::Ray;
voronoi_edge.first_vertex = triangle_vertex[triangle];
voronoi_edge.point = result.vertices[voronoi_edge.first_vertex];
voronoi_edge.direction = detail::unit(outward);
} else {
assert(incident_count[edge] == 0);
Point<long double> first(sites[first_site]);
Point<long double> second(sites[second_site]);
Point<long double> edge_direction = second - first;
voronoi_edge.kind = VoronoiEdgeKind::Line;
voronoi_edge.point = (first + second) / 2.0L;
voronoi_edge.direction = detail::unit(Point<long double>(
edge_direction.y,
-edge_direction.x
));
}
int voronoi_edge_index = int(result.edges.size());
result.edges.push_back(voronoi_edge);
result.cell_edges[first_site].push_back(voronoi_edge_index);
result.cell_edges[second_site].push_back(voronoi_edge_index);
}
for (int site = 0; site < size; ++site) {
std::sort(
result.cell_edges[site].begin(),
result.cell_edges[site].end(),
[&](int first_edge, int second_edge) {
int first_other =
detail::other_site(result.edges[first_edge], site);
int second_other =
detail::other_site(result.edges[second_edge], site);
return detail::direction_less(
sites,
site,
first_other,
second_other
);
}
);
}
return result;
}
} // namespace geometry
} // namespace m1une