Convex Decomposition
(geometry/convex_decomposition.hpp)
- View this file on GitHub
- Last update: 2026-10-06 02:48:54+09:00
- Include:
#include "geometry/convex_decomposition.hpp"
Overview
Rational coordinates are supported, including Rational<BigInt>. Predicates
that use exact arithmetic for integral inputs also use exact rational arithmetic
and ignore eps. Constructed coordinates and lengths with long double return
types remain approximations. Complexity bounds count scalar operations; rational
inputs add the underlying arithmetic and gcd costs.
This header partitions a simple polygon without holes into convex polygons. It provides two choices:
-
convex_decompositiontriangulates the polygon and applies the Hertel–Mehlhorn diagonal-removal algorithm. It is the routine to use under ordinary contest constraints. -
minimum_convex_decompositionuses Keil–Snoeyink dynamic programming over visible diagonals and narrowest-pair stacks. It returns the smallest possible number of pieces. Its running time is sensitive to the number of reflex vertices rather than being quartic in the total vertex count.
Both algorithms return an exact geometric partition: the pieces have disjoint
interiors and their union is exactly the input polygon. The approximation
guarantee of convex_decomposition concerns only how many pieces it returns
compared with the minimum possible number.
Both algorithms use the no-Steiner-point model: every output vertex is an
input vertex, apart from removing redundant boundary vertices. This restriction
is important when interpreting the piece-count guarantee and the word
“minimum.” For a floating-point decomposition that may introduce Steiner
vertices and guarantees fewer than twice the unrestricted optimum, see
steiner_convex_decomposition.
Functions
template <Coordinate T>
std::optional<std::vector<std::vector<Point<T>>>>
convex_decomposition(
std::vector<Point<T>> polygon,
long double eps = 1e-12L
);
template <Coordinate T>
std::optional<std::vector<std::vector<Point<T>>>>
minimum_convex_decomposition(
std::vector<Point<T>> polygon,
long double eps = 1e-12L
);
| Function | Result | Time | Memory |
|---|---|---|---|
convex_decomposition(polygon, eps) |
An exact partition containing at most four times the minimum number of convex pieces. | $O(N^2)$ | $O(N)$ |
minimum_convex_decomposition(polygon, eps) |
A decomposition with the minimum possible number of pieces. | $O(N + \min\lbrace NR^2, R^4 \rbrace)$ | $O(N + \min\lbrace NR^2, R^4 \rbrace)$ |
Here $N$ is the number of vertices after cleanup and $R$ is the number of reflex vertices (vertices whose interior angle is greater than $\pi$).
The bound for the exact routine is easy to misread. The implementation chooses between these two Keil–Snoeyink paths:
- The direct narrowest-pair DP stores only subproblems with at least one reflex endpoint and takes $O(NR^2)$ time and memory.
- When $R^2 < N$, biased-decomposition reduction retains the $O(R^2)$ vertices that can be endpoints of reflex extensions and diagonal extensions, and runs the same DP on that reduced polygon. Including reduction and reconstruction, this path takes $O(N + R^4)$ time and memory.
The minimum of those paths is therefore $O(N + \min\lbrace NR^2, R^4 \rbrace)$. The reduction is internal: output pieces still contain the original polygon boundary, including vertices omitted from the DP instance. $R$ is not the number of vertices in the reduced polygon.
Input and output rules
The input may be clockwise or counterclockwise. A repeated closing point,
consecutive duplicates, and redundant collinear boundary vertices are removed.
The polygon must be simple and have no holes. The exact routine treats
simplicity as a precondition so that it does not add an $O(N^2)$ validation step
to its reflex-sensitive bound. convex_decomposition validates simplicity.
The return value is nullopt when fewer than three effective vertices remain,
the area is zero, validation fails where performed, or construction fails.
Every returned polygon is counterclockwise, does not repeat its first point at the end, and is convex in the non-strict sense. The pieces have disjoint interiors and their union is the original closed polygon.
“Minimum” means the minimum number of non-strictly convex pieces, not that each piece is strictly convex. Collinear vertices can therefore occur on an output boundary. The approximation factor compares against the same no-Steiner-point optimum.
For standard integral coordinate types up to 64 bits, the exact routine performs no floating-point conversion. It scans the cleaned polygon once and selects the smallest integer type that is provably wide enough for its projective visibility predicates and rational ray-intersection comparisons:
| Largest coordinate magnitude | Predicate type |
|---|---|
| Fewer than $2^{30}$ | built-in __int128_t
|
| Fewer than $2^{62}$ | Int256 |
| Otherwise | Int512 |
The largest intermediate has magnitude below $2^{4k+6}$ when coordinate
magnitudes use at most $k$ bits. Thus full-range 32-bit coordinates are not
dispatched to __int128_t, while the common bound
$|x|,|y|\le 10^9$ is.
Selection uses absolute coordinates rather than edge lengths because
homogeneous line coefficients are formed before translation-dependent terms
cancel; translating a polygon far from the origin can therefore select a wider
type. The scan costs $O(N)$ and does not change the overall bound.
Input turns still use the geometry module’s ordinary widened integer type, so
results are exact as long as those ordinary predicates do not overflow. eps
has no effect on integral predicates.
For floating-point coordinates, visibility and ray shooting use long double,
and eps controls predicate tolerance.
Example
#include "geometry/convex_decomposition.hpp"
#include <iostream>
#include <vector>
int main() {
using Point = m1une::geometry::Point<long long>;
std::vector<Point> polygon;
polygon.emplace_back(0, 0);
polygon.emplace_back(5, 0);
polygon.emplace_back(5, 2);
polygon.emplace_back(2, 2);
polygon.emplace_back(2, 5);
polygon.emplace_back(0, 5);
auto parts = m1une::geometry::convex_decomposition(polygon);
if (!parts.has_value()) return 0;
std::cout << parts->size() << "\n"; // 2
}
Depends on
- Circles
(geometry/circle.hpp)
- geometry/detail/floating_predicate.hpp
- Linear Objects
(geometry/linear.hpp)
- 2D Point and Predicates
(geometry/point.hpp)
- Polygons
(geometry/polygon.hpp)
- utilities/detail/fixed_int.hpp
- Int256
(utilities/int256.hpp)
- Int512
(utilities/int512.hpp)
Required by
Verified with
- verify/geometry/centroid.test.cpp
- verify/geometry/convex_decomposition.test.cpp
- verify/geometry/geometry_algorithms.test.cpp
- verify/geometry/rational.test.cpp
Code
#ifndef M1UNE_GEOMETRY_CONVEX_DECOMPOSITION_HPP
#define M1UNE_GEOMETRY_CONVEX_DECOMPOSITION_HPP 1
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <deque>
#include <limits>
#include <map>
#include <memory>
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>
#include "../utilities/int256.hpp"
#include "../utilities/int512.hpp"
#include "polygon.hpp"
namespace m1une {
namespace geometry {
namespace convex_decomposition_detail {
using Index = std::size_t;
using IndexPolygon = std::vector<Index>;
enum class ExactPredicateWidth { Int128, Int256, Int512 };
template <std::integral T>
std::size_t coordinate_magnitude_bits(
const std::vector<Point<T>>& polygon
) {
using Unsigned = std::make_unsigned_t<T>;
std::size_t result = 0;
auto update = [&](T coordinate) {
Unsigned magnitude = static_cast<Unsigned>(coordinate);
if constexpr (std::signed_integral<T>) {
if (coordinate < 0) magnitude = Unsigned(0) - magnitude;
}
std::size_t bits = 0;
while (magnitude != 0) {
++bits;
magnitude >>= 1;
}
result = std::max(result, bits);
};
for (const Point<T>& point : polygon) {
update(point.x);
update(point.y);
}
return result;
}
template <std::integral T>
ExactPredicateWidth select_exact_predicate_width(
const std::vector<Point<T>>& polygon
) {
static_assert(sizeof(T) <= sizeof(std::uint64_t));
const std::size_t coordinate_bits =
coordinate_magnitude_bits(polygon);
// Cross-multiplied ray parameters and projective visibility
// determinants have magnitude below 2^(4 * coordinate_bits + 6).
const std::size_t required_bits = 4 * coordinate_bits + 7;
if (required_bits <= 128) return ExactPredicateWidth::Int128;
if (required_bits <= 256) return ExactPredicateWidth::Int256;
assert(required_bits <= 512);
return ExactPredicateWidth::Int512;
}
template <typename Number, Coordinate T>
Number predicate_number(T value) {
return static_cast<Number>(value);
}
template <Coordinate T, typename Number>
int predicate_sign(const Number& value, long double eps) {
if constexpr (ExactCoordinate<T>) {
return (value > 0) - (value < 0);
} else {
return (value > eps) - (value < -eps);
}
}
template <Coordinate T>
std::optional<std::vector<Point<T>>> prepare_polygon(
std::vector<Point<T>> polygon,
long double eps
) {
polygon =
polygon_detail::clean_polygon_vertices(std::move(polygon), eps);
if (polygon.size() < 3) return std::nullopt;
const int area_sign = sign<T>(polygon_area2(polygon), eps);
if (area_sign == 0 || !is_simple_polygon(polygon, eps)) {
return std::nullopt;
}
if (area_sign < 0) std::reverse(polygon.begin(), polygon.end());
return polygon;
}
template <Coordinate T>
bool is_weakly_convex(
const std::vector<Point<T>>& polygon,
long double eps
) {
if (polygon.size() < 3) return false;
for (Index index = 0; index < polygon.size(); ++index) {
if (
orientation(
polygon[index],
polygon[(index + 1) % polygon.size()],
polygon[(index + 2) % polygon.size()],
eps
) < 0
) {
return false;
}
}
return true;
}
template <Coordinate T>
std::optional<std::vector<IndexPolygon>> triangulate_indices(
const std::vector<Point<T>>& polygon,
long double eps
) {
const Index size = polygon.size();
std::vector<Index> previous(size), next(size);
std::vector<bool> active(size, true);
for (Index index = 0; index < size; ++index) {
previous[index] = (index + size - 1) % size;
next[index] = (index + 1) % size;
}
auto is_ear = [&](Index index) {
if (!active[index]) return false;
const Index first = previous[index];
const Index third = next[index];
if (
orientation(
polygon[first], polygon[index], polygon[third], eps
) <= 0
) {
return false;
}
for (Index other = 0; other < size; ++other) {
if (
!active[other] || other == first || other == index ||
other == third
) {
continue;
}
if (
polygon_detail::in_ccw_triangle(
polygon[other],
polygon[first],
polygon[index],
polygon[third],
eps
)
) {
return false;
}
}
return true;
};
std::deque<Index> ears;
for (Index index = 0; index < size; ++index) {
if (is_ear(index)) ears.push_back(index);
}
std::vector<IndexPolygon> triangles;
triangles.reserve(size - 2);
Index remaining = size;
while (remaining > 3) {
while (!ears.empty() && !is_ear(ears.front())) {
ears.pop_front();
}
if (ears.empty()) return std::nullopt;
const Index ear = ears.front();
ears.pop_front();
const Index first = previous[ear];
const Index third = next[ear];
triangles.push_back(IndexPolygon{first, ear, third});
active[ear] = false;
next[first] = third;
previous[third] = first;
--remaining;
if (is_ear(first)) ears.push_back(first);
if (is_ear(third)) ears.push_back(third);
}
Index first = 0;
while (first < size && !active[first]) ++first;
if (first == size) return std::nullopt;
const Index second = next[first];
const Index third = next[second];
if (
third == first || next[third] != first ||
orientation(
polygon[first], polygon[second], polygon[third], eps
) <= 0
) {
return std::nullopt;
}
triangles.push_back(IndexPolygon{first, second, third});
return triangles;
}
inline Index find_root(std::vector<Index>& parent, Index index) {
Index root = index;
while (parent[root] != root) root = parent[root];
while (parent[index] != index) {
const Index next = parent[index];
parent[index] = root;
index = next;
}
return root;
}
inline std::optional<IndexPolygon> merge_across_edge(
const IndexPolygon& first,
const IndexPolygon& second,
Index edge_first,
Index edge_second
) {
for (Index first_position = 0;
first_position < first.size();
++first_position) {
const Index first_next = (first_position + 1) % first.size();
const Index from = first[first_position];
const Index to = first[first_next];
if (
!(
(from == edge_first && to == edge_second) ||
(from == edge_second && to == edge_first)
)
) {
continue;
}
for (Index second_position = 0;
second_position < second.size();
++second_position) {
const Index second_next =
(second_position + 1) % second.size();
if (
second[second_position] != to ||
second[second_next] != from
) {
continue;
}
IndexPolygon merged;
merged.reserve(first.size() + second.size() - 2);
for (Index position = first_next;
position != first_position;
position = (position + 1) % first.size()) {
merged.push_back(first[position]);
}
for (Index position = second_next;
position != second_position;
position = (position + 1) % second.size()) {
merged.push_back(second[position]);
}
return merged;
}
}
return std::nullopt;
}
template <Coordinate T>
bool is_weakly_convex(
const IndexPolygon& polygon,
const std::vector<Point<T>>& points,
long double eps
) {
for (Index index = 0; index < polygon.size(); ++index) {
if (
orientation(
points[polygon[index]],
points[polygon[(index + 1) % polygon.size()]],
points[polygon[(index + 2) % polygon.size()]],
eps
) < 0
) {
return false;
}
}
return true;
}
template <Coordinate T>
std::vector<Point<T>> materialize(
const IndexPolygon& indices,
const std::vector<Point<T>>& points,
long double eps
) {
std::vector<Point<T>> polygon;
polygon.reserve(indices.size());
for (const Index index : indices) polygon.push_back(points[index]);
return polygon_detail::clean_polygon_vertices(std::move(polygon), eps);
}
struct Diagonal {
Index first;
Index second;
};
template <Coordinate T>
std::optional<std::vector<Point<T>>> prepare_minimum_polygon(
std::vector<Point<T>> polygon,
long double eps
) {
if (polygon.size() >= 2 && polygon.front() == polygon.back()) {
polygon.pop_back();
}
std::vector<Point<T>> distinct;
distinct.reserve(polygon.size());
for (const Point<T>& point : polygon) {
if (distinct.empty() || distinct.back() != point) {
distinct.push_back(point);
}
}
if (distinct.size() >= 2 && distinct.front() == distinct.back()) {
distinct.pop_back();
}
if (distinct.size() < 3) return std::nullopt;
const int original_size = static_cast<int>(distinct.size());
std::vector<int> previous(original_size), next(original_size);
std::vector<bool> removed(original_size, false);
std::deque<int> candidates;
for (int index = 0; index < original_size; ++index) {
previous[index] = (index + original_size - 1) % original_size;
next[index] = (index + 1) % original_size;
candidates.push_back(index);
}
int remaining = original_size;
while (!candidates.empty() && remaining >= 3) {
const int index = candidates.front();
candidates.pop_front();
if (removed[index]) continue;
const int before = previous[index];
const int after = next[index];
if (
orientation(
distinct[before], distinct[index], distinct[after], eps
) != 0 ||
sign<T>(
dot(
distinct[index] - distinct[before],
distinct[after] - distinct[index]
),
eps
) < 0
) {
continue;
}
removed[index] = true;
next[before] = after;
previous[after] = before;
--remaining;
candidates.push_back(before);
candidates.push_back(after);
}
if (remaining < 3) return std::nullopt;
std::vector<Point<T>> cleaned;
cleaned.reserve(static_cast<Index>(remaining));
int first = 0;
while (removed[first]) ++first;
int index = first;
do {
cleaned.push_back(distinct[index]);
index = next[index];
} while (index != first);
const int area_sign = sign<T>(polygon_area2(cleaned), eps);
if (area_sign == 0) return std::nullopt;
if (area_sign < 0) std::reverse(cleaned.begin(), cleaned.end());
return cleaned;
}
template <Coordinate T, typename Number>
class BiasedPolygonReduction {
private:
struct Vector {
Number x;
Number y;
};
struct Fraction {
Number numerator;
Number denominator;
};
struct Chain {
int first_edge;
int edge_count;
};
public:
struct Result {
std::vector<Point<T>> polygon;
std::vector<Index> original_index;
};
BiasedPolygonReduction(
const std::vector<Point<T>>& polygon,
long double eps
)
: polygon_(polygon),
eps_(eps),
size_(static_cast<int>(polygon.size())),
marked_edge_(polygon.size(), false) {
for (int index = 0; index < size_; ++index) {
if (
orientation(
polygon_[(index + size_ - 1) % size_],
polygon_[index],
polygon_[(index + 1) % size_],
eps_
) < 0
) {
reflex_vertices_.push_back(index);
}
}
build_chains();
}
Result run() {
for (const int reflex : reflex_vertices_) {
marked_edge_[(reflex + size_ - 1) % size_] = true;
marked_edge_[reflex] = true;
}
std::vector<std::vector<int>> extension_endpoints(size_);
for (const int reflex : reflex_vertices_) {
const Vector direction = vector_between(
reflex, (reflex + 1) % size_
);
for (const int edge : ray_shoot(reflex, direction)) {
marked_edge_[edge] = true;
extension_endpoints[reflex].push_back(edge);
extension_endpoints[reflex].push_back((edge + 1) % size_);
}
}
for (Index first = 0; first < reflex_vertices_.size(); ++first) {
for (Index second = first + 1;
second < reflex_vertices_.size();
++second) {
const int first_vertex = reflex_vertices_[first];
const int second_vertex = reflex_vertices_[second];
mark_ray(
first_vertex,
vector_between(first_vertex, second_vertex)
);
mark_ray(
second_vertex,
vector_between(second_vertex, first_vertex)
);
}
}
for (const int reflex : reflex_vertices_) {
auto& endpoints = extension_endpoints[reflex];
std::sort(endpoints.begin(), endpoints.end());
endpoints.erase(
std::unique(endpoints.begin(), endpoints.end()),
endpoints.end()
);
for (const int endpoint : endpoints) {
if (endpoint == reflex) continue;
mark_ray(
reflex, vector_between(reflex, endpoint)
);
}
}
Result result;
for (int index = 0; index < size_; ++index) {
if (
marked_edge_[index] ||
marked_edge_[(index + size_ - 1) % size_]
) {
result.polygon.push_back(polygon_[index]);
result.original_index.push_back(static_cast<Index>(index));
}
}
return result;
}
private:
const std::vector<Point<T>>& polygon_;
long double eps_;
int size_;
std::vector<int> reflex_vertices_;
std::vector<Chain> chains_;
std::vector<bool> marked_edge_;
Vector vector_between(int first, int second) const {
return {
predicate_number<Number>(polygon_[first].x) -
predicate_number<Number>(polygon_[second].x),
predicate_number<Number>(polygon_[first].y) -
predicate_number<Number>(polygon_[second].y)
};
}
Number vector_cross(const Vector& first, const Vector& second) const {
return first.x * second.y - first.y * second.x;
}
Number vector_dot(const Vector& first, const Vector& second) const {
return first.x * second.x + first.y * second.y;
}
int vector_sign(const Number& value) const {
return predicate_sign<T>(value, eps_);
}
int quadrant(const Vector& direction) const {
const int x_sign = vector_sign(direction.x);
const int y_sign = vector_sign(direction.y);
if (y_sign >= 0) return x_sign >= 0 ? 0 : 1;
return x_sign < 0 ? 2 : 3;
}
void build_chains() {
// Within one quadrant, monotonically turning edge directions span at
// most pi/2. This gives the bitonic sidedness needed by ray shooting
// without evaluating an angle.
int first_edge = 0;
Vector previous_direction = vector_between(1, 0);
int current_quadrant = quadrant(previous_direction);
for (int edge = 1; edge < size_; ++edge) {
const Vector direction = vector_between(
(edge + 1) % size_, edge
);
const int direction_quadrant = quadrant(direction);
if (
vector_sign(vector_cross(previous_direction, direction)) <
0 ||
direction_quadrant != current_quadrant
) {
chains_.push_back(Chain{
first_edge, edge - first_edge
});
first_edge = edge;
current_quadrant = direction_quadrant;
}
previous_direction = direction;
}
chains_.push_back(Chain{first_edge, size_ - first_edge});
}
int side(
int origin,
const Vector& direction,
int vertex
) const {
return vector_sign(
vector_cross(
direction,
vector_between(vertex % size_, origin)
)
);
}
int edge_side(const Vector& direction, int edge) const {
return vector_sign(
vector_cross(
direction,
vector_between((edge + 1) % size_, edge)
)
);
}
void crossing_on_monotone_part(
const Chain& chain,
int origin,
const Vector& direction,
int first_position,
int last_position,
std::vector<int>& candidates
) const {
if (first_position >= last_position) return;
const int first_sign = side(
origin,
direction,
chain.first_edge + first_position
);
const int last_sign = side(
origin,
direction,
chain.first_edge + last_position
);
if (first_sign == 0) {
candidates.push_back(
chain.first_edge + first_position
);
}
if (last_sign == 0) {
candidates.push_back(
chain.first_edge + last_position - 1
);
}
if (first_sign == 0 || last_sign == 0 || first_sign == last_sign) {
return;
}
int low = first_position;
int high = last_position;
while (high - low > 1) {
const int middle = (low + high) / 2;
const int middle_sign = side(
origin,
direction,
chain.first_edge + middle
);
if (middle_sign == 0 || middle_sign != first_sign) {
high = middle;
} else {
low = middle;
}
}
candidates.push_back(chain.first_edge + high - 1);
}
void chain_candidates(
const Chain& chain,
int origin,
const Vector& direction,
std::vector<int>& candidates
) const {
int split = 0;
const int first_derivative =
edge_side(direction, chain.first_edge);
const int last_derivative = edge_side(
direction,
chain.first_edge + chain.edge_count - 1
);
if (
first_derivative != 0 && last_derivative != 0 &&
first_derivative != last_derivative
) {
int low = 0;
int high = chain.edge_count - 1;
while (high - low > 1) {
const int middle = (low + high) / 2;
const int middle_sign = edge_side(
direction, chain.first_edge + middle
);
if (
middle_sign == 0 ||
middle_sign != first_derivative
) {
high = middle;
} else {
low = middle;
}
}
split = high;
}
if (split == 0) {
crossing_on_monotone_part(
chain,
origin,
direction,
0,
chain.edge_count,
candidates
);
} else {
crossing_on_monotone_part(
chain, origin, direction, 0, split, candidates
);
crossing_on_monotone_part(
chain,
origin,
direction,
split,
chain.edge_count,
candidates
);
}
}
std::vector<int> ray_shoot(
int origin,
const Vector& direction
) const {
std::vector<int> candidates;
for (const Chain& chain : chains_) {
chain_candidates(chain, origin, direction, candidates);
}
std::sort(candidates.begin(), candidates.end());
candidates.erase(
std::unique(candidates.begin(), candidates.end()),
candidates.end()
);
if constexpr (ExactCoordinate<T>) {
return exact_ray_shoot(origin, direction, candidates);
} else {
return floating_ray_shoot(origin, direction, candidates);
}
}
std::vector<int> exact_ray_shoot(
int origin,
const Vector& direction,
const std::vector<int>& candidates
) const {
// Intersections are rational parameters along the ray. Keep their
// numerators and denominators and compare them by cross products.
std::optional<Fraction> best;
std::vector<int> result;
const Number direction_norm2 = vector_dot(direction, direction);
for (int edge : candidates) {
edge %= size_;
const Vector first = vector_between(edge, origin);
const Vector second = vector_between(
(edge + 1) % size_, origin
);
const Vector edge_direction = vector_between(
(edge + 1) % size_, edge
);
Number denominator = vector_cross(direction, edge_direction);
Fraction parameter{0, 1};
bool valid = false;
if (denominator == 0) {
if (vector_cross(direction, first) != 0) continue;
const Number first_parameter =
vector_dot(first, direction);
const Number second_parameter =
vector_dot(second, direction);
if (first_parameter > 0) {
parameter = Fraction{
first_parameter, direction_norm2
};
valid = true;
}
if (
second_parameter > 0 &&
(!valid || second_parameter < parameter.numerator)
) {
parameter = Fraction{
second_parameter, direction_norm2
};
valid = true;
}
} else {
Number numerator = vector_cross(first, edge_direction);
Number segment_numerator = vector_cross(first, direction);
if (denominator < 0) {
denominator = -denominator;
numerator = -numerator;
segment_numerator = -segment_numerator;
}
if (
numerator <= 0 || segment_numerator < 0 ||
segment_numerator > denominator
) {
continue;
}
parameter = Fraction{numerator, denominator};
valid = true;
}
if (!valid) continue;
if (!best.has_value()) {
best = parameter;
result.assign(1, edge);
continue;
}
const Number left =
parameter.numerator * best->denominator;
const Number right =
best->numerator * parameter.denominator;
if (left < right) {
best = parameter;
result.assign(1, edge);
} else if (left == right) {
result.push_back(edge);
}
}
return result;
}
std::vector<int> floating_ray_shoot(
int origin,
const Vector& direction,
const std::vector<int>& candidates
) const {
long double best = std::numeric_limits<long double>::infinity();
std::vector<int> result;
const long double direction_norm2 = vector_dot(
direction, direction
);
for (int edge : candidates) {
edge %= size_;
const Vector first = vector_between(edge, origin);
const Vector second = vector_between(
(edge + 1) % size_, origin
);
const Vector edge_direction = vector_between(
(edge + 1) % size_, edge
);
const long double denominator = vector_cross(
direction, edge_direction
);
long double parameter = -1;
if (std::fabs(denominator) <= eps_) {
if (std::fabs(vector_cross(direction, first)) > eps_) {
continue;
}
const long double first_parameter =
vector_dot(first, direction) / direction_norm2;
const long double second_parameter =
vector_dot(second, direction) / direction_norm2;
if (first_parameter > eps_) parameter = first_parameter;
if (
second_parameter > eps_ &&
(parameter < 0 || second_parameter < parameter)
) {
parameter = second_parameter;
}
} else {
parameter =
vector_cross(first, edge_direction) / denominator;
const long double segment_parameter =
vector_cross(first, direction) / denominator;
if (
parameter <= eps_ || segment_parameter < -eps_ ||
segment_parameter > 1 + eps_
) {
continue;
}
}
if (parameter < 0) continue;
if (parameter + eps_ < best) {
best = parameter;
result.assign(1, edge);
} else if (std::fabs(parameter - best) <= eps_) {
result.push_back(edge);
}
}
return result;
}
void mark_ray(int origin, const Vector& direction) {
for (const int edge : ray_shoot(origin, direction)) {
marked_edge_[edge] = true;
}
}
};
template <Coordinate T, typename Number>
class KeilSnoeyinkDecomposition {
private:
static constexpr int infinity = 100000000;
static constexpr int bad = 1000000000;
struct SolutionNode;
using NodePointer = std::shared_ptr<const SolutionNode>;
struct SolutionNode {
NodePointer first;
NodePointer second;
Diagonal diagonal{0, 0};
bool is_diagonal = false;
};
struct NarrowPair {
int first;
int second;
NodePointer solution;
};
// Two independently popped views of the same narrowest-pair stack.
// Pairs are appended from back to front in the terminology of the
// Keil--Snoeyink paper.
struct PairDeque {
std::vector<NarrowPair> pairs;
int front = -1;
int back = 0;
bool empty_front() const { return front < 0; }
bool more_front() const { return front > 0; }
bool empty_back() const {
return back >= static_cast<int>(pairs.size());
}
bool more_back() const {
return back + 1 < static_cast<int>(pairs.size());
}
const NarrowPair& front_pair() const { return pairs[front]; }
const NarrowPair& under_front() const {
return pairs[front - 1];
}
const NarrowPair& back_pair() const { return pairs[back]; }
const NarrowPair& under_back() const {
return pairs[back + 1];
}
void pop_front() { --front; }
void pop_back() { ++back; }
void restore() {
front = static_cast<int>(pairs.size()) - 1;
back = 0;
}
void clear() {
pairs.clear();
restore();
}
void push(int first, int second, NodePointer solution = nullptr) {
pairs.push_back(NarrowPair{first, second, std::move(solution)});
restore();
}
void push_narrow(
int first,
int second,
NodePointer solution
) {
if (!empty_front() && first <= front_pair().first) return;
while (!empty_front() && front_pair().second >= second) {
pairs.pop_back();
--front;
}
push(first, second, std::move(solution));
}
};
struct State {
int weight = bad;
PairDeque pairs;
};
using ProjectiveNumber = Number;
struct Homogeneous {
ProjectiveNumber w;
ProjectiveNumber x;
ProjectiveNumber y;
Homogeneous negated() const { return {-w, -x, -y}; }
ProjectiveNumber side(const Point<T>& point) const {
return w + x * predicate_number<ProjectiveNumber>(point.x) +
y * predicate_number<ProjectiveNumber>(point.y);
}
static Homogeneous line(
const Point<T>& first,
const Point<T>& second
) {
const ProjectiveNumber ax =
predicate_number<ProjectiveNumber>(first.x);
const ProjectiveNumber ay =
predicate_number<ProjectiveNumber>(first.y);
const ProjectiveNumber bx =
predicate_number<ProjectiveNumber>(second.x);
const ProjectiveNumber by =
predicate_number<ProjectiveNumber>(second.y);
return {ax * by - ay * bx, ay - by, bx - ax};
}
static ProjectiveNumber determinant(
const Homogeneous& first,
const Homogeneous& second,
const Homogeneous& third
) {
return
first.w *
(second.x * third.y - second.y * third.x) -
second.w *
(first.x * third.y - first.y * third.x) +
third.w *
(first.x * second.y - first.y * second.x);
}
};
// ElGindy--Avis visibility-polygon scan, specialized to a polygon
// vertex. Only original polygon vertices are returned; artificial lid
// intersections are discarded.
class VisibilityPolygon {
private:
enum VertexType { right_lid, left_lid, right_wall, left_wall };
public:
VisibilityPolygon(
const std::vector<Point<T>>& polygon,
long double eps
)
: polygon_(polygon), eps_(eps), size_(polygon.size()) {}
std::vector<Index> build(Index origin_index) {
origin_ = &polygon_[origin_index];
vertices_.clear();
types_.clear();
left_lid_index_ = not_saved;
right_lid_index_ = not_saved;
int vertex = static_cast<int>(origin_index);
push(vertex++, right_wall);
do {
push(vertex++, right_wall);
if (vertex >= static_cast<int>(size_ + origin_index)) {
break;
}
Homogeneous edge = line(vertex - 1, vertex);
if (left(edge, *origin_)) continue;
if (!left(edge, point(vertex - 2))) {
vertex = exit_right_bay(
vertex,
top(),
Homogeneous{1, 0, 0}
);
push(vertex++, right_lid);
continue;
}
save_lid();
for (;;) {
if (
orientation(
*origin_, top(), point(vertex), eps_
) > 0
) {
if (
orientation(
point(vertex),
point(vertex + 1),
*origin_,
eps_
) < 0
) {
++vertex;
} else if (left(edge, point(vertex + 1))) {
vertex = exit_left_bay(
vertex,
point(vertex),
line(left_lid_index_, left_lid_index_ - 1)
) + 1;
} else {
restore_lid();
push(vertex++, left_wall);
break;
}
edge = line(vertex - 1, vertex);
} else if (!left(edge, top())) {
if (right_lid_index_ == not_saved) break;
vertex = exit_right_bay(
vertex, top(), edge.negated()
);
push(vertex++, right_lid);
break;
} else {
save_lid();
}
}
} while (vertex < static_cast<int>(size_ + origin_index));
std::vector<Index> result;
while (!vertices_.empty()) {
while (
!vertices_.empty() &&
(types_.back() == right_lid ||
types_.back() == left_lid)
) {
vertices_.pop_back();
types_.pop_back();
}
if (vertices_.empty()) break;
result.push_back(
static_cast<Index>(vertices_.back()) % size_
);
vertices_.pop_back();
types_.pop_back();
}
return result;
}
private:
static constexpr int not_saved = -1;
const std::vector<Point<T>>& polygon_;
long double eps_;
Index size_;
const Point<T>* origin_ = nullptr;
std::vector<int> vertices_;
std::vector<VertexType> types_;
int left_lid_index_ = not_saved;
int right_lid_index_ = not_saved;
const Point<T>& point(int index) const {
int reduced = index % static_cast<int>(size_);
if (reduced < 0) reduced += static_cast<int>(size_);
return polygon_[static_cast<Index>(reduced)];
}
const Point<T>& top() const { return point(vertices_.back()); }
Homogeneous line(int first, int second) const {
return Homogeneous::line(point(first), point(second));
}
bool left(
const Homogeneous& line_value,
const Point<T>& point_value
) const {
return
predicate_sign<T>(line_value.side(point_value), eps_) > 0;
}
bool right(
const Homogeneous& line_value,
const Point<T>& point_value
) const {
return
predicate_sign<T>(line_value.side(point_value), eps_) < 0;
}
bool clockwise(
const Homogeneous& first,
const Homogeneous& second,
const Homogeneous& third
) const {
return predicate_sign<T>(
Homogeneous::determinant(first, second, third), eps_
) < 0;
}
void push(int index, VertexType type) {
vertices_.push_back(index);
types_.push_back(type);
}
void save_lid() {
if (types_.back() == left_wall) {
vertices_.pop_back();
types_.pop_back();
}
left_lid_index_ = vertices_.back();
vertices_.pop_back();
types_.pop_back();
if (types_.back() == right_lid) {
right_lid_index_ = vertices_.back();
vertices_.pop_back();
types_.pop_back();
} else {
right_lid_index_ = not_saved;
}
}
void restore_lid() {
if (right_lid_index_ != not_saved) {
push(right_lid_index_, right_lid);
}
push(left_lid_index_, left_lid);
}
int exit_right_bay(
int vertex,
const Point<T>& bottom,
const Homogeneous& lid
) const {
int winding = 0;
const Homogeneous mouth = Homogeneous::line(*origin_, bottom);
bool current_left = false;
while (++vertex < 3 * static_cast<int>(size_)) {
const bool last_left = current_left;
current_left = left(mouth, point(vertex));
if (
current_left != last_left &&
(orientation(
point(vertex - 1),
point(vertex),
*origin_,
eps_
) > 0) == current_left
) {
if (!current_left) {
--winding;
} else if (winding++ == 0) {
const Homogeneous edge = line(vertex - 1, vertex);
if (
left(edge, bottom) &&
!clockwise(mouth, edge, lid)
) {
return vertex - 1;
}
}
}
}
return vertex;
}
int exit_left_bay(
int vertex,
const Point<T>& bottom,
const Homogeneous& lid
) const {
int winding = 0;
const Homogeneous mouth = Homogeneous::line(*origin_, bottom);
bool current_right = false;
while (++vertex < 3 * static_cast<int>(size_)) {
const bool last_right = current_right;
current_right = right(mouth, point(vertex));
if (
current_right != last_right &&
(orientation(
point(vertex - 1),
point(vertex),
*origin_,
eps_
) < 0) == current_right
) {
if (!current_right) {
++winding;
} else if (winding-- == 0) {
const Homogeneous edge = line(vertex - 1, vertex);
if (
right(edge, bottom) &&
!clockwise(mouth, edge, lid)
) {
return vertex - 1;
}
}
}
}
return vertex;
}
};
public:
KeilSnoeyinkDecomposition(
const std::vector<Point<T>>& polygon,
long double eps
)
: polygon_(polygon),
eps_(eps),
size_(polygon.size()),
reflex_(size_, false),
remapped_(size_) {}
std::optional<std::vector<Diagonal>> run() {
initialize_reflex_vertices();
allocate_states();
initialize_visibility();
initialize_base_subproblems();
for (int span = 3; span < static_cast<int>(size_); ++span) {
for (const int first : reflex_vertices_) {
const int last = first + span;
if (last >= static_cast<int>(size_)) break;
if (!visible(first, last)) continue;
initialize_pairs(first, last);
if (reflex_[last]) {
for (int middle = first + 1; middle < last; ++middle) {
type_a(first, middle, last);
}
} else {
for (const int middle : reflex_vertices_) {
if (middle <= first) continue;
if (middle >= last - 1) break;
type_a(first, middle, last);
}
type_a(first, last - 1, last);
}
}
for (const int last : reflex_vertices_) {
if (last < span) continue;
const int first = last - span;
if (reflex_[first] || !visible(first, last)) continue;
initialize_pairs(first, last);
type_b(first, first + 1, last);
for (const int middle : reflex_vertices_) {
if (middle <= first + 1) continue;
if (middle >= last) break;
type_b(first, middle, last);
}
}
}
if (weight(0, static_cast<int>(size_) - 1) >= infinity) {
return std::nullopt;
}
std::vector<Diagonal> diagonals;
const PairDeque& root = state(0, static_cast<int>(size_) - 1).pairs;
if (root.pairs.empty()) return std::nullopt;
flatten(root.pairs.front().solution, diagonals);
for (Diagonal& diagonal : diagonals) {
if (diagonal.second < diagonal.first) {
std::swap(diagonal.first, diagonal.second);
}
}
std::sort(
diagonals.begin(),
diagonals.end(),
[](const Diagonal& first, const Diagonal& second) {
return std::pair(first.first, first.second) <
std::pair(second.first, second.second);
}
);
diagonals.erase(
std::unique(
diagonals.begin(),
diagonals.end(),
[](const Diagonal& first, const Diagonal& second) {
return
first.first == second.first &&
first.second == second.second;
}
),
diagonals.end()
);
if (
static_cast<int>(diagonals.size()) !=
weight(0, static_cast<int>(size_) - 1)
) {
return std::nullopt;
}
return diagonals;
}
private:
const std::vector<Point<T>>& polygon_;
long double eps_;
Index size_;
std::vector<bool> reflex_;
std::vector<int> reflex_vertices_;
std::vector<int> remapped_;
std::vector<std::vector<State>> states_;
void initialize_reflex_vertices() {
for (int index = 0; index < static_cast<int>(size_); ++index) {
reflex_[index] = index == 0 || orientation(
polygon_[(index + static_cast<int>(size_) - 1) % size_],
polygon_[index],
polygon_[(index + 1) % size_],
eps_
) < 0;
if (reflex_[index]) reflex_vertices_.push_back(index);
}
}
void allocate_states() {
int next = 0;
for (const int index : reflex_vertices_) remapped_[index] = next++;
for (int index = 0; index < static_cast<int>(size_); ++index) {
if (!reflex_[index]) remapped_[index] = next++;
}
const int reflex_count = static_cast<int>(reflex_vertices_.size());
states_.resize(size_);
for (int index = 0; index < static_cast<int>(size_); ++index) {
states_[index].resize(reflex_[index] ? size_ : reflex_count);
}
}
State& state(int first, int last) {
assert(first <= last);
assert(reflex_[first] || reflex_[last]);
return states_[first][remapped_[last]];
}
const State& state(int first, int last) const {
assert(first <= last);
assert(reflex_[first] || reflex_[last]);
return states_[first][remapped_[last]];
}
int weight(int first, int last) const {
return state(first, last).weight;
}
bool visible(int first, int last) const {
return weight(first, last) < bad;
}
void initialize_visibility() {
VisibilityPolygon visibility_polygon(polygon_, eps_);
for (const int reflex : reflex_vertices_) {
for (const Index visible_vertex :
visibility_polygon.build(static_cast<Index>(reflex))) {
int first = reflex;
int last = static_cast<int>(visible_vertex);
if (last < first) std::swap(first, last);
if (first == last) continue;
state(first, last).weight = infinity;
}
}
// The visibility scan treats the closing edge like every other
// boundary edge, but make the DP convention explicit.
state(0, static_cast<int>(size_) - 1).weight = infinity;
}
void initialize_base_subproblems() {
const int size = static_cast<int>(size_);
for (const int index : reflex_vertices_) {
if (index + 1 < size) state(index, index + 1).weight = 0;
if (index > 0) state(index - 1, index).weight = 0;
if (index + 2 < size && visible(index, index + 2)) {
State& base = state(index, index + 2);
base.weight = 0;
base.pairs.clear();
base.pairs.push(index + 1, index + 1, nullptr);
}
if (index >= 2 && visible(index - 2, index)) {
State& base = state(index - 2, index);
base.weight = 0;
base.pairs.clear();
base.pairs.push(index - 1, index - 1, nullptr);
}
}
}
void initialize_pairs(int first, int last) {
state(first, last).pairs.clear();
}
void update(
int first,
int last,
int new_weight,
int pair_first,
int pair_second,
NodePointer solution
) {
State& subproblem = state(first, last);
if (new_weight > subproblem.weight) return;
if (new_weight < subproblem.weight) {
subproblem.weight = new_weight;
subproblem.pairs.clear();
}
subproblem.pairs.push_narrow(
pair_first, pair_second, std::move(solution)
);
}
NodePointer concatenate(NodePointer first, NodePointer second) const {
if (!first) return second;
if (!second) return first;
SolutionNode node;
node.first = std::move(first);
node.second = std::move(second);
return std::make_shared<const SolutionNode>(std::move(node));
}
NodePointer append_diagonal(
NodePointer solution,
int first,
int last
) const {
SolutionNode leaf;
leaf.diagonal = Diagonal{
static_cast<Index>(first), static_cast<Index>(last)
};
leaf.is_diagonal = true;
return concatenate(
std::move(solution),
std::make_shared<const SolutionNode>(std::move(leaf))
);
}
NodePointer any_solution(int first, int last) const {
const PairDeque& pairs = state(first, last).pairs;
assert(!pairs.pairs.empty());
return pairs.pairs.front().solution;
}
void type_a(int first, int middle, int last) {
if (!visible(first, middle)) return;
int top = middle;
int new_weight = weight(first, middle);
NodePointer solution;
if (last - middle > 1) {
if (!visible(middle, last)) return;
new_weight += weight(middle, last) + 1;
solution = any_solution(middle, last);
solution = append_diagonal(
std::move(solution), middle, last
);
}
bool use_first_diagonal = false;
if (middle - first > 1) {
PairDeque& pairs = state(first, middle).pairs;
if (pairs.empty_back()) return;
if (
orientation(
polygon_[last],
polygon_[middle],
polygon_[pairs.back_pair().second],
eps_
) <= 0
) {
while (
pairs.more_back() &&
orientation(
polygon_[last],
polygon_[middle],
polygon_[pairs.under_back().second],
eps_
) <= 0
) {
pairs.pop_back();
}
if (
!pairs.empty_back() &&
orientation(
polygon_[last],
polygon_[first],
polygon_[pairs.back_pair().first],
eps_
) >= 0
) {
top = pairs.back_pair().first;
} else {
++new_weight;
use_first_diagonal = true;
}
} else {
++new_weight;
use_first_diagonal = true;
}
solution = concatenate(
pairs.back_pair().solution, std::move(solution)
);
if (use_first_diagonal) {
solution = append_diagonal(
std::move(solution), first, middle
);
}
}
update(
first,
last,
new_weight,
top,
middle,
std::move(solution)
);
}
void type_b(int first, int middle, int last) {
if (!visible(middle, last)) return;
int top = middle;
int new_weight = weight(middle, last);
NodePointer solution;
if (middle - first > 1) {
if (!visible(first, middle)) return;
new_weight += weight(first, middle) + 1;
solution = any_solution(first, middle);
solution = append_diagonal(
std::move(solution), first, middle
);
}
bool use_second_diagonal = false;
if (last - middle > 1) {
PairDeque& pairs = state(middle, last).pairs;
if (pairs.empty_front()) return;
if (
orientation(
polygon_[first],
polygon_[middle],
polygon_[pairs.front_pair().first],
eps_
) >= 0
) {
while (
pairs.more_front() &&
orientation(
polygon_[first],
polygon_[middle],
polygon_[pairs.under_front().first],
eps_
) >= 0
) {
pairs.pop_front();
}
if (
!pairs.empty_front() &&
orientation(
polygon_[first],
polygon_[last],
polygon_[pairs.front_pair().second],
eps_
) <= 0
) {
top = pairs.front_pair().second;
} else {
++new_weight;
use_second_diagonal = true;
}
} else {
++new_weight;
use_second_diagonal = true;
}
solution = concatenate(
std::move(solution), pairs.front_pair().solution
);
if (use_second_diagonal) {
solution = append_diagonal(
std::move(solution), middle, last
);
}
}
update(
first,
last,
new_weight,
middle,
top,
std::move(solution)
);
}
void flatten(
const NodePointer& node,
std::vector<Diagonal>& diagonals
) const {
if (!node) return;
if (node->is_diagonal) {
diagonals.push_back(node->diagonal);
return;
}
flatten(node->first, diagonals);
flatten(node->second, diagonals);
}
};
template <Coordinate T>
std::optional<std::vector<IndexPolygon>> faces_from_diagonals(
const std::vector<Point<T>>& polygon,
const std::vector<Diagonal>& diagonals,
long double eps
) {
const Index size = polygon.size();
std::vector<std::vector<Index>> adjacent(size);
auto add_edge = [&](Index first, Index second) {
adjacent[first].push_back(second);
adjacent[second].push_back(first);
};
for (Index index = 0; index < size; ++index) {
add_edge(index, (index + 1) % size);
}
for (const Diagonal& diagonal : diagonals) {
if (
diagonal.first >= size || diagonal.second >= size ||
diagonal.first == diagonal.second
) {
return std::nullopt;
}
add_edge(diagonal.first, diagonal.second);
}
for (Index origin = 0; origin < size; ++origin) {
auto upper_half = [&](Index target) {
const auto direction = polygon[target] - polygon[origin];
const int y_sign = sign<T>(direction.y, eps);
return
y_sign > 0 ||
(y_sign == 0 && sign<T>(direction.x, eps) >= 0);
};
std::sort(
adjacent[origin].begin(),
adjacent[origin].end(),
[&](Index first, Index second) {
const bool first_upper = upper_half(first);
const bool second_upper = upper_half(second);
if (first_upper != second_upper) return first_upper;
const auto first_direction =
polygon[first] - polygon[origin];
const auto second_direction =
polygon[second] - polygon[origin];
const int turn = sign<T>(
cross(first_direction, second_direction), eps
);
if (turn != 0) return turn > 0;
return
norm2(first_direction) < norm2(second_direction);
}
);
adjacent[origin].erase(
std::unique(
adjacent[origin].begin(), adjacent[origin].end()
),
adjacent[origin].end()
);
}
std::vector<std::vector<bool>> visited(size);
Index directed_edge_count = 0;
for (Index index = 0; index < size; ++index) {
visited[index].assign(adjacent[index].size(), false);
directed_edge_count += adjacent[index].size();
}
std::vector<IndexPolygon> result;
for (Index start = 0; start < size; ++start) {
for (Index start_position = 0;
start_position < adjacent[start].size();
++start_position) {
if (visited[start][start_position]) continue;
IndexPolygon face;
Index current = start;
Index position = start_position;
for (Index guard = 0;; ++guard) {
if (guard > directed_edge_count) return std::nullopt;
if (visited[current][position]) {
if (current != start || position != start_position) {
return std::nullopt;
}
break;
}
visited[current][position] = true;
face.push_back(current);
const Index next = adjacent[current][position];
const auto reverse = std::find(
adjacent[next].begin(), adjacent[next].end(), current
);
if (reverse == adjacent[next].end()) return std::nullopt;
const Index reverse_position =
static_cast<Index>(reverse - adjacent[next].begin());
current = next;
position =
(reverse_position + adjacent[current].size() - 1) %
adjacent[current].size();
}
if (face.size() < 3) continue;
std::vector<Point<T>> points;
points.reserve(face.size());
for (const Index index : face) points.push_back(polygon[index]);
if (sign<T>(polygon_area2(points), eps) > 0) {
if (!is_weakly_convex(face, polygon, eps)) {
return std::nullopt;
}
result.push_back(std::move(face));
}
}
}
if (result.size() != diagonals.size() + 1) return std::nullopt;
return result;
}
} // namespace convex_decomposition_detail
template <Coordinate T>
std::optional<std::vector<std::vector<Point<T>>>>
convex_decomposition(
std::vector<Point<T>> polygon,
long double eps = 1e-12L
) {
auto prepared = convex_decomposition_detail::prepare_polygon(
std::move(polygon), eps
);
if (!prepared.has_value()) return std::nullopt;
polygon = std::move(*prepared);
if (convex_decomposition_detail::is_weakly_convex(polygon, eps)) {
return std::vector<std::vector<Point<T>>>{std::move(polygon)};
}
auto triangulation =
convex_decomposition_detail::triangulate_indices(polygon, eps);
if (!triangulation.has_value()) return std::nullopt;
using convex_decomposition_detail::Index;
using convex_decomposition_detail::IndexPolygon;
const Index triangle_count = triangulation->size();
const Index absent = std::numeric_limits<Index>::max();
struct Owners {
Index first = std::numeric_limits<Index>::max();
Index second = std::numeric_limits<Index>::max();
};
std::map<std::pair<Index, Index>, Owners> edge_owners;
for (Index triangle = 0; triangle < triangle_count; ++triangle) {
for (Index edge = 0; edge < 3; ++edge) {
Index first = (*triangulation)[triangle][edge];
Index second = (*triangulation)[triangle][(edge + 1) % 3];
if (second < first) std::swap(first, second);
Owners& owners = edge_owners[{first, second}];
if (owners.first == absent) {
owners.first = triangle;
} else {
owners.second = triangle;
}
}
}
std::vector<Index> parent(triangle_count);
std::vector<IndexPolygon> pieces = std::move(*triangulation);
for (Index index = 0; index < triangle_count; ++index) {
parent[index] = index;
}
for (const auto& [edge, owners] : edge_owners) {
if (owners.second == absent) continue;
Index first_root =
convex_decomposition_detail::find_root(parent, owners.first);
Index second_root =
convex_decomposition_detail::find_root(parent, owners.second);
if (first_root == second_root) continue;
auto merged = convex_decomposition_detail::merge_across_edge(
pieces[first_root],
pieces[second_root],
edge.first,
edge.second
);
if (
!merged.has_value() ||
!convex_decomposition_detail::is_weakly_convex(
*merged, polygon, eps
)
) {
continue;
}
pieces[first_root] = std::move(*merged);
pieces[second_root].clear();
parent[second_root] = first_root;
}
std::vector<std::vector<Point<T>>> result;
for (Index index = 0; index < triangle_count; ++index) {
if (convex_decomposition_detail::find_root(parent, index) != index) {
continue;
}
result.push_back(convex_decomposition_detail::materialize(
pieces[index], polygon, eps
));
}
return result;
}
template <Coordinate T>
std::optional<std::vector<std::vector<Point<T>>>>
minimum_convex_decomposition(
std::vector<Point<T>> polygon,
long double eps = 1e-12L
) {
auto prepared = convex_decomposition_detail::prepare_minimum_polygon(
std::move(polygon), eps
);
if (!prepared.has_value()) return std::nullopt;
polygon = std::move(*prepared);
if (convex_decomposition_detail::is_weakly_convex(polygon, eps)) {
return std::vector<std::vector<Point<T>>>{std::move(polygon)};
}
std::size_t first_reflex = polygon.size();
for (std::size_t index = 0; index < polygon.size(); ++index) {
if (
orientation(
polygon[(index + polygon.size() - 1) % polygon.size()],
polygon[index],
polygon[(index + 1) % polygon.size()],
eps
) < 0
) {
first_reflex = index;
break;
}
}
assert(first_reflex < polygon.size());
std::rotate(
polygon.begin(), polygon.begin() + first_reflex, polygon.end()
);
std::size_t reflex_count = 0;
for (std::size_t index = 0; index < polygon.size(); ++index) {
if (
orientation(
polygon[(index + polygon.size() - 1) % polygon.size()],
polygon[index],
polygon[(index + 1) % polygon.size()],
eps
) < 0
) {
++reflex_count;
}
}
auto solve = [&]<typename Number>()
-> std::optional<std::vector<std::vector<Point<T>>>> {
std::vector<Point<T>> dynamic_polygon = polygon;
std::vector<convex_decomposition_detail::Index> original_index(
polygon.size()
);
for (std::size_t index = 0; index < polygon.size(); ++index) {
original_index[index] = index;
}
if (
reflex_count > 0 &&
polygon.size() / reflex_count > reflex_count
) {
convex_decomposition_detail::BiasedPolygonReduction<T, Number>
reduction(polygon, eps);
auto reduced = reduction.run();
dynamic_polygon = std::move(reduced.polygon);
original_index = std::move(reduced.original_index);
}
convex_decomposition_detail::KeilSnoeyinkDecomposition<T, Number>
solver(dynamic_polygon, eps);
auto diagonals = solver.run();
if (!diagonals.has_value()) return std::nullopt;
for (auto& diagonal : *diagonals) {
diagonal.first = original_index[diagonal.first];
diagonal.second = original_index[diagonal.second];
}
auto index_decomposition =
convex_decomposition_detail::faces_from_diagonals(
polygon, *diagonals, eps
);
if (!index_decomposition.has_value()) return std::nullopt;
std::vector<std::vector<Point<T>>> result;
result.reserve(index_decomposition->size());
for (const auto& indices : *index_decomposition) {
result.push_back(convex_decomposition_detail::materialize(
indices, polygon, eps
));
}
return result;
};
if constexpr (std::integral<T>) {
using convex_decomposition_detail::ExactPredicateWidth;
switch (
convex_decomposition_detail::select_exact_predicate_width(
polygon
)
) {
case ExactPredicateWidth::Int128:
return solve.template operator()<__int128_t>();
case ExactPredicateWidth::Int256:
return solve.template operator()<
m1une::utilities::Int256
>();
case ExactPredicateWidth::Int512:
return solve.template operator()<
m1une::utilities::Int512
>();
}
assert(false);
return std::nullopt;
} else {
return solve.template operator()<wide_type<T>>();
}
}
} // namespace geometry
} // namespace m1une
#endif // M1UNE_GEOMETRY_CONVEX_DECOMPOSITION_HPP#line 1 "geometry/convex_decomposition.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <deque>
#include <limits>
#include <map>
#include <memory>
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>
#line 1 "utilities/int256.hpp"
#include <string>
#include <string_view>
#line 1 "utilities/detail/fixed_int.hpp"
#line 5 "utilities/detail/fixed_int.hpp"
#include <array>
#include <concepts>
#line 9 "utilities/detail/fixed_int.hpp"
#include <istream>
#include <ostream>
#include <stdexcept>
#line 16 "utilities/detail/fixed_int.hpp"
namespace m1une {
namespace utilities {
namespace detail {
// A signed two's-complement integer whose arithmetic wraps modulo 2^Bits.
// Public aliases select contest-friendly fixed widths in int*.hpp.
template <std::size_t Bits>
class FixedInt {
static_assert(Bits >= 64);
static_assert(Bits % 64 == 0);
private:
static constexpr std::size_t limb_count = Bits / 64;
using LimbArray = std::array<std::uint64_t, limb_count>;
public:
static constexpr std::size_t bit_width = Bits;
constexpr FixedInt() = default;
template <std::integral Integer>
constexpr FixedInt(Integer value) {
static_assert(sizeof(Integer) <= sizeof(std::uint64_t));
if constexpr (std::signed_integral<Integer>) {
const std::uint64_t extension =
value < 0 ? ~std::uint64_t(0) : std::uint64_t(0);
limbs_.fill(extension);
limbs_[0] = static_cast<std::uint64_t>(
static_cast<std::int64_t>(value)
);
} else {
limbs_[0] = static_cast<std::uint64_t>(value);
}
}
explicit FixedInt(std::string_view text) { read(text); }
FixedInt& operator=(std::string_view text) {
read(text);
return *this;
}
void read(std::string_view text) {
if (text.empty()) {
throw std::invalid_argument("empty fixed-width integer");
}
const bool negative = text.front() == '-';
std::size_t position =
(text.front() == '-' || text.front() == '+') ? 1 : 0;
if (position == text.size()) {
throw std::invalid_argument("invalid fixed-width integer");
}
FixedInt result;
for (; position < text.size(); ++position) {
const char digit = text[position];
if (digit < '0' || digit > '9') {
throw std::invalid_argument("invalid fixed-width integer");
}
result.multiply_unsigned_small(10);
result += FixedInt(static_cast<unsigned>(digit - '0'));
}
*this = negative ? -result : result;
}
constexpr bool is_zero() const {
for (const std::uint64_t limb : limbs_) {
if (limb != 0) return false;
}
return true;
}
constexpr bool is_negative() const {
return (limbs_.back() >> 63) != 0;
}
constexpr int sign() const {
if (is_zero()) return 0;
return is_negative() ? -1 : 1;
}
constexpr FixedInt operator+() const { return *this; }
constexpr FixedInt operator-() const {
FixedInt result;
result.limbs_ = limbs_;
negate_unsigned(result.limbs_);
return result;
}
constexpr FixedInt& operator+=(const FixedInt& other) {
__uint128_t carry = 0;
for (std::size_t index = 0; index < limb_count; ++index) {
const __uint128_t current =
__uint128_t(limbs_[index]) + other.limbs_[index] + carry;
limbs_[index] = static_cast<std::uint64_t>(current);
carry = current >> 64;
}
return *this;
}
constexpr FixedInt& operator-=(const FixedInt& other) {
return *this += -other;
}
constexpr FixedInt& operator*=(const FixedInt& other) {
LimbArray product{};
for (std::size_t first = 0; first < limb_count; ++first) {
__uint128_t carry = 0;
for (
std::size_t second = 0;
first + second < limb_count;
++second
) {
const std::size_t position = first + second;
const __uint128_t current =
__uint128_t(limbs_[first]) * other.limbs_[second] +
product[position] + carry;
product[position] = static_cast<std::uint64_t>(current);
carry = current >> 64;
}
}
limbs_ = product;
return *this;
}
constexpr FixedInt& multiply_small(std::uint64_t value) {
multiply_unsigned_small(value);
return *this;
}
constexpr FixedInt& operator/=(const FixedInt& other) {
return *this = divmod(*this, other).first;
}
constexpr FixedInt& operator%=(const FixedInt& other) {
return *this = divmod(*this, other).second;
}
std::string to_string() const {
if (is_zero()) return "0";
const bool negative = is_negative();
LimbArray magnitude = unsigned_magnitude();
constexpr std::size_t chunk_capacity =
(Bits * 30103 / 100000 + 9) / 9;
std::array<std::uint32_t, chunk_capacity> chunks{};
std::size_t chunk_count = 0;
while (!magnitude_is_zero(magnitude)) {
chunks[chunk_count++] = static_cast<std::uint32_t>(
divide_unsigned_by_small(magnitude, 1000000000)
);
}
std::string result;
result.reserve(chunk_count * 9 + negative);
if (negative) result.push_back('-');
result += std::to_string(chunks[chunk_count - 1]);
char digits[9];
while (--chunk_count != 0) {
std::uint32_t chunk = chunks[chunk_count - 1];
for (int index = 8; index >= 0; --index) {
digits[index] = static_cast<char>('0' + chunk % 10);
chunk /= 10;
}
result.append(digits, 9);
}
return result;
}
friend constexpr std::pair<FixedInt, FixedInt> divmod(
const FixedInt& dividend,
const FixedInt& divisor
) {
if (divisor.is_zero()) {
throw std::domain_error("fixed-width integer division by zero");
}
const bool quotient_negative =
dividend.is_negative() != divisor.is_negative();
const bool remainder_negative = dividend.is_negative();
auto [quotient_limbs, remainder_limbs] = divide_unsigned(
dividend.unsigned_magnitude(), divisor.unsigned_magnitude()
);
FixedInt quotient;
FixedInt remainder;
quotient.limbs_ = quotient_limbs;
remainder.limbs_ = remainder_limbs;
if (quotient_negative) quotient = -quotient;
if (remainder_negative) remainder = -remainder;
return std::make_pair(quotient, remainder);
}
friend constexpr std::pair<FixedInt, std::int64_t> divmod_small(
const FixedInt& dividend,
std::uint32_t divisor
) {
if (divisor == 0) {
throw std::domain_error("fixed-width integer division by zero");
}
LimbArray quotient_limbs = dividend.unsigned_magnitude();
const std::uint64_t unsigned_remainder =
divide_unsigned_by_small(quotient_limbs, divisor);
FixedInt quotient;
quotient.limbs_ = quotient_limbs;
if (dividend.is_negative()) quotient = -quotient;
const std::int64_t remainder = dividend.is_negative()
? -std::int64_t(unsigned_remainder)
: std::int64_t(unsigned_remainder);
return std::make_pair(quotient, remainder);
}
friend constexpr std::int64_t mod_small(
const FixedInt& dividend,
std::uint32_t divisor
) {
if (divisor == 0) {
throw std::domain_error("fixed-width integer division by zero");
}
const bool negative = dividend.is_negative();
const std::uint64_t remainder = negative
? remainder_unsigned_by_small(
dividend.unsigned_magnitude(),
divisor
)
: remainder_unsigned_by_small(
dividend.limbs_, divisor
);
return negative ? -std::int64_t(remainder)
: std::int64_t(remainder);
}
friend constexpr FixedInt operator+(
FixedInt first,
const FixedInt& second
) {
return first += second;
}
friend constexpr FixedInt operator-(
FixedInt first,
const FixedInt& second
) {
return first -= second;
}
friend constexpr FixedInt operator*(
FixedInt first,
const FixedInt& second
) {
return first *= second;
}
friend constexpr FixedInt operator/(
FixedInt first,
const FixedInt& second
) {
return first /= second;
}
friend constexpr FixedInt operator%(
FixedInt first,
const FixedInt& second
) {
return first %= second;
}
friend constexpr bool operator==(
const FixedInt& first,
const FixedInt& second
) = default;
friend constexpr bool operator<(
const FixedInt& first,
const FixedInt& second
) {
const bool first_negative = first.is_negative();
const bool second_negative = second.is_negative();
if (first_negative != second_negative) return first_negative;
return compare_unsigned(first.limbs_, second.limbs_) < 0;
}
friend constexpr bool operator!=(
const FixedInt& first,
const FixedInt& second
) {
return !(first == second);
}
friend constexpr bool operator>(
const FixedInt& first,
const FixedInt& second
) {
return second < first;
}
friend constexpr bool operator<=(
const FixedInt& first,
const FixedInt& second
) {
return !(second < first);
}
friend constexpr bool operator>=(
const FixedInt& first,
const FixedInt& second
) {
return !(first < second);
}
friend std::ostream& operator<<(
std::ostream& output,
const FixedInt& value
) {
return output << value.to_string();
}
friend std::istream& operator>>(
std::istream& input,
FixedInt& value
) {
std::string text;
if (input >> text) value.read(text);
return input;
}
private:
LimbArray limbs_{};
constexpr LimbArray unsigned_magnitude() const {
LimbArray result = limbs_;
if (is_negative()) negate_unsigned(result);
return result;
}
constexpr void multiply_unsigned_small(std::uint64_t value) {
__uint128_t carry = 0;
for (std::size_t index = 0; index < limb_count; ++index) {
const __uint128_t current =
__uint128_t(limbs_[index]) * value + carry;
limbs_[index] = static_cast<std::uint64_t>(current);
carry = current >> 64;
}
}
static constexpr void negate_unsigned(LimbArray& value) {
for (std::uint64_t& limb : value) limb = ~limb;
for (std::size_t index = 0; index < limb_count; ++index) {
if (++value[index] != 0) break;
}
}
static constexpr int compare_unsigned(
const LimbArray& first,
const LimbArray& second
) {
for (std::size_t offset = 0; offset < limb_count; ++offset) {
const std::size_t index = limb_count - 1 - offset;
if (first[index] != second[index]) {
return first[index] < second[index] ? -1 : 1;
}
}
return 0;
}
static constexpr void subtract_unsigned(
LimbArray& first,
const LimbArray& second
) {
std::uint64_t borrow = 0;
for (std::size_t index = 0; index < limb_count; ++index) {
const std::uint64_t previous = first[index];
first[index] -= second[index] + borrow;
const bool addition_overflow =
borrow != 0 && second[index] == ~std::uint64_t(0);
borrow = addition_overflow ||
previous < second[index] + borrow;
}
}
static constexpr void shift_left_one(LimbArray& value) {
std::uint64_t carry = 0;
for (std::size_t index = 0; index < limb_count; ++index) {
const std::uint64_t next_carry = value[index] >> 63;
value[index] = (value[index] << 1) | carry;
carry = next_carry;
}
}
static constexpr std::pair<LimbArray, LimbArray> divide_unsigned(
const LimbArray& dividend,
const LimbArray& divisor
) {
LimbArray quotient{};
LimbArray remainder{};
for (std::size_t offset = 0; offset < Bits; ++offset) {
const std::size_t bit = Bits - 1 - offset;
shift_left_one(remainder);
remainder[0] |=
(dividend[bit / 64] >> (bit % 64)) & std::uint64_t(1);
if (compare_unsigned(remainder, divisor) >= 0) {
subtract_unsigned(remainder, divisor);
quotient[bit / 64] |= std::uint64_t(1) << (bit % 64);
}
}
return std::make_pair(quotient, remainder);
}
static bool magnitude_is_zero(const LimbArray& value) {
for (const std::uint64_t limb : value) {
if (limb != 0) return false;
}
return true;
}
static constexpr std::uint64_t divide_unsigned_by_small(
LimbArray& value,
std::uint32_t divisor
) {
std::uint64_t remainder = 0;
for (std::size_t offset = 0; offset < limb_count; ++offset) {
const std::size_t index = limb_count - 1 - offset;
const std::uint64_t high =
(remainder << 32) | (value[index] >> 32);
const std::uint64_t quotient_high = high / divisor;
remainder = high % divisor;
const std::uint64_t low =
(remainder << 32) | std::uint32_t(value[index]);
const std::uint64_t quotient_low = low / divisor;
remainder = low % divisor;
value[index] = (quotient_high << 32) | quotient_low;
}
return remainder;
}
static constexpr std::uint64_t remainder_unsigned_by_small(
const LimbArray& value,
std::uint32_t divisor
) {
std::uint64_t remainder = 0;
for (std::size_t offset = 0; offset < limb_count; ++offset) {
const std::size_t index = limb_count - 1 - offset;
remainder = ((remainder << 32) | (value[index] >> 32)) % divisor;
remainder =
((remainder << 32) | std::uint32_t(value[index])) % divisor;
}
return remainder;
}
};
} // namespace detail
} // namespace utilities
} // namespace m1une
#line 8 "utilities/int256.hpp"
namespace m1une {
namespace utilities {
using Int256 = detail::FixedInt<256>;
using i256 = Int256;
inline Int256 parse_int256(std::string_view text) {
return Int256(text);
}
inline std::string to_string(const Int256& value) {
return value.to_string();
}
} // namespace utilities
} // namespace m1une
#line 1 "utilities/int512.hpp"
#line 6 "utilities/int512.hpp"
#line 8 "utilities/int512.hpp"
namespace m1une {
namespace utilities {
using Int512 = detail::FixedInt<512>;
using i512 = Int512;
inline Int512 parse_int512(std::string_view text) {
return Int512(text);
}
inline std::string to_string(const Int512& value) {
return value.to_string();
}
} // namespace utilities
} // namespace m1une
#line 1 "geometry/polygon.hpp"
#line 10 "geometry/polygon.hpp"
#include <numbers>
#line 14 "geometry/polygon.hpp"
#line 1 "geometry/circle.hpp"
#line 13 "geometry/circle.hpp"
#line 1 "geometry/linear.hpp"
#line 7 "geometry/linear.hpp"
#line 1 "geometry/point.hpp"
#line 8 "geometry/point.hpp"
#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 9 "geometry/linear.hpp"
namespace m1une {
namespace geometry {
template <Coordinate T>
struct Line {
Point<T> a;
Point<T> b;
};
template <Coordinate T>
struct Segment {
Point<T> a;
Point<T> b;
};
template <Coordinate T>
struct Ray {
Point<T> origin;
Point<T> through;
};
enum class LinearIntersectionKind {
Empty,
Point,
Segment,
Ray,
Line,
};
struct LinearIntersection {
LinearIntersectionKind kind;
Point<long double> first;
Point<long double> second;
};
struct ClosestPoints {
Point<long double> first;
Point<long double> second;
};
namespace linear_intersection_detail {
inline LinearIntersection make_empty() {
const Point<long double> zero;
return LinearIntersection{
LinearIntersectionKind::Empty,
zero,
zero,
};
}
template <Coordinate T>
LinearIntersection make_point(const Point<T>& point) {
const Point<long double> converted(point);
return LinearIntersection{
LinearIntersectionKind::Point,
converted,
converted,
};
}
template <Coordinate T>
LinearIntersection make_object(
LinearIntersectionKind kind,
const Point<T>& first,
const Point<T>& second
) {
return LinearIntersection{
kind,
Point<long double>(first),
Point<long double>(second),
};
}
} // namespace linear_intersection_detail
template <Coordinate T>
constexpr Point<long double> centroid(const Segment<T>& segment) {
return Point<long double>(
(
static_cast<long double>(segment.a.x) +
static_cast<long double>(segment.b.x)
) / 2,
(
static_cast<long double>(segment.a.y) +
static_cast<long double>(segment.b.y)
) / 2
);
}
template <Coordinate T>
bool on_line(
const Line<T>& line,
const Point<T>& point,
long double eps = 1e-12L
) {
assert(line.a != line.b);
return orientation(line.a, line.b, point, eps) == 0;
}
template <Coordinate T>
bool parallel(const Line<T>& first, const Line<T>& second, long double eps = 1e-12L) {
using W = wide_type<T>;
W first_x = W(first.b.x) - W(first.a.x);
W first_y = W(first.b.y) - W(first.a.y);
W second_x = W(second.b.x) - W(second.a.x);
W second_y = W(second.b.y) - W(second.a.y);
return predicate_detail::determinant_sign<ExactCoordinate<T>>(
first_x,
first_y,
second_x,
second_y,
eps
) == 0;
}
template <Coordinate T>
bool orthogonal(const Line<T>& first, const Line<T>& second, long double eps = 1e-12L) {
using W = wide_type<T>;
W first_x = W(first.b.x) - W(first.a.x);
W first_y = W(first.b.y) - W(first.a.y);
W second_x = W(second.b.x) - W(second.a.x);
W second_y = W(second.b.y) - W(second.a.y);
return predicate_detail::dot_sign<ExactCoordinate<T>>(
first_x,
first_y,
second_x,
second_y,
eps
) == 0;
}
template <Coordinate T>
Point<long double> projection(const Line<T>& line, const Point<T>& point) {
assert(line.a != line.b);
Point<long double> a(line.a);
Point<long double> direction(
static_cast<long double>(line.b.x) - static_cast<long double>(line.a.x),
static_cast<long double>(line.b.y) - static_cast<long double>(line.a.y)
);
Point<long double> offset(
static_cast<long double>(point.x) - a.x,
static_cast<long double>(point.y) - a.y
);
long double ratio = dot(offset, direction) / dot(direction, direction);
return a + direction * ratio;
}
template <Coordinate T>
Point<long double> reflection(const Line<T>& line, const Point<T>& point) {
Point<long double> projected = projection(line, point);
return projected * 2.0L - Point<long double>(point);
}
template <Coordinate T>
bool intersects(
const Line<T>& first,
const Line<T>& second,
long double eps = 1e-12L
) {
return !parallel(first, second, eps) || on_line(first, second.a, eps);
}
template <Coordinate T>
bool on_segment(
const Segment<T>& segment,
const Point<T>& point,
long double eps = 1e-12L
) {
if (orientation(segment.a, segment.b, point, eps) != 0) return false;
using W = wide_type<T>;
const W direction_x = W(segment.b.x) - W(segment.a.x);
const W direction_y = W(segment.b.y) - W(segment.a.y);
if (direction_x == W(0) && direction_y == W(0)) {
if constexpr (ExactCoordinate<T>) {
return point == segment.a;
} else {
return
predicate_detail::absolute(W(point.x) - W(segment.a.x)) <= eps &&
predicate_detail::absolute(W(point.y) - W(segment.a.y)) <= eps;
}
}
const W offset_x = W(point.x) - W(segment.a.x);
const W offset_y = W(point.y) - W(segment.a.y);
const W projection =
offset_x * direction_x + offset_y * direction_y;
const W length_squared =
direction_x * direction_x + direction_y * direction_y;
return
predicate_detail::scaled_sign<ExactCoordinate<T>>(
projection,
length_squared,
eps
) >= 0 &&
predicate_detail::scaled_sign<ExactCoordinate<T>>(
projection - length_squared,
length_squared,
eps
) <= 0;
}
template <Coordinate T>
Point<long double> projection(
const Segment<T>& segment,
const Point<T>& point
) {
const Point<long double> first(segment.a);
const Point<long double> direction =
Point<long double>(segment.b) - first;
const long double length_squared = dot(direction, direction);
if (length_squared == 0) return first;
const long double ratio = std::clamp(
dot(Point<long double>(point) - first, direction) / length_squared,
0.0L,
1.0L
);
return first + direction * ratio;
}
template <Coordinate T>
bool intersects(
const Segment<T>& first,
const Segment<T>& second,
long double eps = 1e-12L
) {
int abc = orientation(first.a, first.b, second.a, eps);
int abd = orientation(first.a, first.b, second.b, eps);
int cda = orientation(second.a, second.b, first.a, eps);
int cdb = orientation(second.a, second.b, first.b, eps);
if (abc == 0 && on_segment(first, second.a, eps)) return true;
if (abd == 0 && on_segment(first, second.b, eps)) return true;
if (cda == 0 && on_segment(second, first.a, eps)) return true;
if (cdb == 0 && on_segment(second, first.b, eps)) return true;
return abc * abd < 0 && cda * cdb < 0;
}
template <Coordinate T>
bool intersects(
const Line<T>& line,
const Segment<T>& segment,
long double eps = 1e-12L
) {
int first_side = orientation(line.a, line.b, segment.a, eps);
int second_side = orientation(line.a, line.b, segment.b, eps);
return first_side == 0 || second_side == 0 || first_side != second_side;
}
template <Coordinate T>
bool intersects(
const Segment<T>& segment,
const Line<T>& line,
long double eps = 1e-12L
) {
return intersects(line, segment, eps);
}
namespace linear_parameter_detail {
template <Coordinate T>
struct Parameters {
wide_type<T> denominator;
wide_type<T> denominator_scale;
wide_type<T> first_numerator;
wide_type<T> second_numerator;
};
template <Coordinate T>
Parameters<T> parameters(
const Point<T>& first_origin,
const Point<T>& first_through,
const Point<T>& second_origin,
const Point<T>& second_through
) {
using W = wide_type<T>;
W first_x = W(first_through.x) - W(first_origin.x);
W first_y = W(first_through.y) - W(first_origin.y);
W second_x = W(second_through.x) - W(second_origin.x);
W second_y = W(second_through.y) - W(second_origin.y);
W offset_x = W(second_origin.x) - W(first_origin.x);
W offset_y = W(second_origin.y) - W(first_origin.y);
return Parameters<T>{
first_x * second_y - first_y * second_x,
predicate_detail::determinant_scale<ExactCoordinate<T>>(
first_x,
first_y,
second_x,
second_y
),
offset_x * second_y - offset_y * second_x,
offset_x * first_y - offset_y * first_x
};
}
template <Coordinate T>
int denominator_sign(const Parameters<T>& values, long double eps) {
return predicate_detail::scaled_sign<ExactCoordinate<T>>(
values.denominator,
values.denominator_scale,
eps
);
}
template <Coordinate T>
bool ratio_nonnegative(
wide_type<T> numerator,
wide_type<T> denominator,
long double eps
) {
const int numerator_sign =
predicate_detail::scaled_sign<ExactCoordinate<T>>(
numerator,
predicate_detail::absolute(denominator),
eps
);
const int denominator_direction =
(denominator > 0) - (denominator < 0);
return
numerator_sign == 0 ||
numerator_sign == denominator_direction;
}
template <Coordinate T>
bool ratio_in_unit_interval(
wide_type<T> numerator,
wide_type<T> denominator,
long double eps
) {
const auto scale = predicate_detail::absolute(denominator);
const int start_sign =
predicate_detail::scaled_sign<ExactCoordinate<T>>(
numerator,
scale,
eps
);
const int finish_sign =
predicate_detail::scaled_sign<ExactCoordinate<T>>(
numerator - denominator,
scale,
eps
);
if (denominator > 0) {
return start_sign >= 0 && finish_sign <= 0;
}
return start_sign <= 0 && finish_sign >= 0;
}
} // namespace linear_parameter_detail
template <Coordinate T>
bool on_ray(
const Ray<T>& ray,
const Point<T>& point,
long double eps = 1e-12L
) {
assert(ray.origin != ray.through);
if (orientation(ray.origin, ray.through, point, eps) != 0) return false;
using W = wide_type<T>;
W direction_x = W(ray.through.x) - W(ray.origin.x);
W direction_y = W(ray.through.y) - W(ray.origin.y);
W offset_x = W(point.x) - W(ray.origin.x);
W offset_y = W(point.y) - W(ray.origin.y);
const W projection =
direction_x * offset_x + direction_y * offset_y;
const W length_squared =
direction_x * direction_x + direction_y * direction_y;
return predicate_detail::scaled_sign<ExactCoordinate<T>>(
projection,
length_squared,
eps
) >= 0;
}
template <Coordinate T>
Point<long double> projection(const Ray<T>& ray, const Point<T>& point) {
assert(ray.origin != ray.through);
Point<long double> origin(ray.origin);
Point<long double> direction =
Point<long double>(ray.through) - origin;
Point<long double> offset = Point<long double>(point) - origin;
long double ratio = dot(offset, direction) / dot(direction, direction);
if (ratio < 0) ratio = 0;
return origin + direction * ratio;
}
template <Coordinate T>
Ray<long double> reflection(const Line<T>& line, const Ray<T>& ray) {
assert(ray.origin != ray.through);
return Ray<long double>{
reflection(line, ray.origin),
reflection(line, ray.through)
};
}
template <Coordinate T>
Ray<long double> reflected_ray(
const Ray<T>& incoming,
const Point<T>& hit,
const Line<T>& mirror,
long double eps = 1e-12L
) {
assert(incoming.origin != incoming.through);
assert(on_line(mirror, hit, eps));
Point<T> translated = hit + (incoming.through - incoming.origin);
return Ray<long double>{
Point<long double>(hit),
reflection(mirror, translated)
};
}
template <Coordinate T>
bool intersects(
const Ray<T>& ray,
const Line<T>& line,
long double eps = 1e-12L
) {
assert(ray.origin != ray.through);
assert(line.a != line.b);
linear_parameter_detail::Parameters<T> values =
linear_parameter_detail::parameters(
ray.origin,
ray.through,
line.a,
line.b
);
if (linear_parameter_detail::denominator_sign(values, eps) == 0) {
return on_line(line, ray.origin, eps);
}
return linear_parameter_detail::ratio_nonnegative<T>(
values.first_numerator,
values.denominator,
eps
);
}
template <Coordinate T>
bool intersects(
const Line<T>& line,
const Ray<T>& ray,
long double eps = 1e-12L
) {
return intersects(ray, line, eps);
}
template <Coordinate T>
bool intersects(
const Ray<T>& ray,
const Segment<T>& segment,
long double eps = 1e-12L
) {
assert(ray.origin != ray.through);
if (segment.a == segment.b) return on_ray(ray, segment.a, eps);
linear_parameter_detail::Parameters<T> values =
linear_parameter_detail::parameters(
ray.origin,
ray.through,
segment.a,
segment.b
);
if (linear_parameter_detail::denominator_sign(values, eps) == 0) {
if (orientation(ray.origin, ray.through, segment.a, eps) != 0) {
return false;
}
return on_ray(ray, segment.a, eps) ||
on_ray(ray, segment.b, eps) ||
on_segment(segment, ray.origin, eps);
}
return linear_parameter_detail::ratio_nonnegative<T>(
values.first_numerator,
values.denominator,
eps
) &&
linear_parameter_detail::ratio_in_unit_interval<T>(
values.second_numerator,
values.denominator,
eps
);
}
template <Coordinate T>
bool intersects(
const Segment<T>& segment,
const Ray<T>& ray,
long double eps = 1e-12L
) {
return intersects(ray, segment, eps);
}
template <Coordinate T>
bool intersects(
const Ray<T>& first,
const Ray<T>& second,
long double eps = 1e-12L
) {
assert(first.origin != first.through);
assert(second.origin != second.through);
linear_parameter_detail::Parameters<T> values =
linear_parameter_detail::parameters(
first.origin,
first.through,
second.origin,
second.through
);
if (linear_parameter_detail::denominator_sign(values, eps) == 0) {
if (orientation(first.origin, first.through, second.origin, eps) != 0) {
return false;
}
return on_ray(first, second.origin, eps) ||
on_ray(second, first.origin, eps);
}
return linear_parameter_detail::ratio_nonnegative<T>(
values.first_numerator,
values.denominator,
eps
) &&
linear_parameter_detail::ratio_nonnegative<T>(
values.second_numerator,
values.denominator,
eps
);
}
namespace linear_intersection_detail {
enum class Domain {
Line,
Segment,
Ray,
};
template <Coordinate T>
struct ParametricObject {
Point<T> origin;
Point<T> through;
Domain domain;
};
template <Coordinate T>
ParametricObject<T> parametric_object(const Line<T>& line) {
assert(line.a != line.b);
return ParametricObject<T>{line.a, line.b, Domain::Line};
}
template <Coordinate T>
ParametricObject<T> parametric_object(const Segment<T>& segment) {
return ParametricObject<T>{segment.a, segment.b, Domain::Segment};
}
template <Coordinate T>
ParametricObject<T> parametric_object(const Ray<T>& ray) {
assert(ray.origin != ray.through);
return ParametricObject<T>{ray.origin, ray.through, Domain::Ray};
}
template <Coordinate T>
bool contains(
const ParametricObject<T>& object,
const Point<T>& point,
long double eps
) {
if (object.domain == Domain::Line) {
return on_line(Line<T>{object.origin, object.through}, point, eps);
}
if (object.domain == Domain::Segment) {
return on_segment(
Segment<T>{object.origin, object.through},
point,
eps
);
}
return on_ray(Ray<T>{object.origin, object.through}, point, eps);
}
template <Coordinate T>
bool accepts_parameter(
Domain domain,
wide_type<T> numerator,
wide_type<T> denominator,
long double eps
) {
if (domain == Domain::Line) return true;
if (domain == Domain::Ray) {
return linear_parameter_detail::ratio_nonnegative<T>(
numerator,
denominator,
eps
);
}
return linear_parameter_detail::ratio_in_unit_interval<T>(
numerator,
denominator,
eps
);
}
template <Coordinate T>
Point<long double> point_at_ratio(
const ParametricObject<T>& object,
wide_type<T> numerator,
wide_type<T> denominator
) {
const long double ratio =
static_cast<long double>(numerator) /
static_cast<long double>(denominator);
const Point<long double> origin(object.origin);
const Point<long double> direction =
Point<long double>(object.through) - origin;
return origin + direction * ratio;
}
template <Coordinate T>
struct AxisProjection {
bool use_x;
bool negate;
wide_type<T> operator()(const Point<T>& point) const {
const wide_type<T> value = use_x
? wide_type<T>(point.x)
: wide_type<T>(point.y);
return negate ? -value : value;
}
};
template <Coordinate T>
AxisProjection<T> axis_projection(const ParametricObject<T>& object) {
using W = wide_type<T>;
const W direction_x = W(object.through.x) - W(object.origin.x);
const W direction_y = W(object.through.y) - W(object.origin.y);
const bool use_x =
predicate_detail::absolute(direction_x) >=
predicate_detail::absolute(direction_y);
const W component = use_x ? direction_x : direction_y;
assert(component != W(0));
return AxisProjection<T>{use_x, component < W(0)};
}
template <Coordinate T>
struct ParameterInterval {
bool has_lower;
bool has_upper;
wide_type<T> lower;
wide_type<T> upper;
};
template <Coordinate T>
ParameterInterval<T> parameter_interval(
const ParametricObject<T>& object,
const AxisProjection<T>& projection
) {
using W = wide_type<T>;
const W origin = projection(object.origin);
const W through = projection(object.through);
if (object.domain == Domain::Line) {
return ParameterInterval<T>{false, false, W(0), W(0)};
}
if (object.domain == Domain::Segment) {
return ParameterInterval<T>{
true,
true,
std::min(origin, through),
std::max(origin, through),
};
}
if (origin < through) {
return ParameterInterval<T>{true, false, origin, W(0)};
}
return ParameterInterval<T>{false, true, W(0), origin};
}
template <Coordinate T>
ParameterInterval<T> intersect_intervals(
ParameterInterval<T> first,
const ParameterInterval<T>& second
) {
if (
second.has_lower &&
(!first.has_lower || first.lower < second.lower)
) {
first.has_lower = true;
first.lower = second.lower;
}
if (
second.has_upper &&
(!first.has_upper || second.upper < first.upper)
) {
first.has_upper = true;
first.upper = second.upper;
}
return first;
}
template <Coordinate T>
Point<long double> point_at_projection(
const ParametricObject<T>& object,
const AxisProjection<T>& projection,
long double target
) {
const long double origin =
static_cast<long double>(projection(object.origin));
const long double through =
static_cast<long double>(projection(object.through));
const long double ratio = (target - origin) / (through - origin);
const Point<long double> point(object.origin);
const Point<long double> direction =
Point<long double>(object.through) - point;
return point + direction * ratio;
}
template <Coordinate T>
LinearIntersection collinear_intersection(
const ParametricObject<T>& first,
const ParametricObject<T>& second,
long double eps
) {
using W = wide_type<T>;
const AxisProjection<T> projection = axis_projection(first);
const ParameterInterval<T> first_interval =
parameter_interval(first, projection);
const ParameterInterval<T> second_interval =
parameter_interval(second, projection);
const ParameterInterval<T> common =
intersect_intervals(first_interval, second_interval);
W scale = predicate_detail::absolute(
projection(first.through) - projection(first.origin)
);
scale = std::max(
scale,
predicate_detail::absolute(
projection(second.through) - projection(second.origin)
)
);
if (common.has_lower && common.has_upper) {
const int order = predicate_detail::scaled_sign<ExactCoordinate<T>>(
common.lower - common.upper,
scale,
eps
);
if (order > 0) return make_empty();
if (order == 0) {
const long double coordinate =
(
static_cast<long double>(common.lower) +
static_cast<long double>(common.upper)
) / 2.0L;
return make_point(
point_at_projection(first, projection, coordinate)
);
}
return make_object(
LinearIntersectionKind::Segment,
point_at_projection(
first,
projection,
static_cast<long double>(common.lower)
),
point_at_projection(
first,
projection,
static_cast<long double>(common.upper)
)
);
}
const Point<long double> direction =
Point<long double>(first.through) -
Point<long double>(first.origin);
if (common.has_lower) {
const Point<long double> origin = point_at_projection(
first,
projection,
static_cast<long double>(common.lower)
);
return make_object(
LinearIntersectionKind::Ray,
origin,
origin + direction
);
}
if (common.has_upper) {
const Point<long double> origin = point_at_projection(
first,
projection,
static_cast<long double>(common.upper)
);
return make_object(
LinearIntersectionKind::Ray,
origin,
origin - direction
);
}
return make_object(
LinearIntersectionKind::Line,
first.origin,
first.through
);
}
template <Coordinate T>
LinearIntersection intersect(
const ParametricObject<T>& first,
const ParametricObject<T>& second,
long double eps
) {
const bool first_degenerate = first.origin == first.through;
const bool second_degenerate = second.origin == second.through;
if (first_degenerate) {
assert(first.domain == Domain::Segment);
if (contains(second, first.origin, eps)) {
return make_point(first.origin);
}
return make_empty();
}
if (second_degenerate) {
assert(second.domain == Domain::Segment);
if (contains(first, second.origin, eps)) {
return make_point(second.origin);
}
return make_empty();
}
const linear_parameter_detail::Parameters<T> values =
linear_parameter_detail::parameters(
first.origin,
first.through,
second.origin,
second.through
);
if (linear_parameter_detail::denominator_sign(values, eps) != 0) {
if (
!accepts_parameter<T>(
first.domain,
values.first_numerator,
values.denominator,
eps
) ||
!accepts_parameter<T>(
second.domain,
values.second_numerator,
values.denominator,
eps
)
) {
return make_empty();
}
return make_point(
point_at_ratio(
first,
values.first_numerator,
values.denominator
)
);
}
if (
orientation(
first.origin,
first.through,
second.origin,
eps
) != 0
) {
return make_empty();
}
return collinear_intersection(first, second, eps);
}
} // namespace linear_intersection_detail
template <Coordinate T>
LinearIntersection linear_intersection(
const Line<T>& first,
const Line<T>& second,
long double eps = 1e-12L
) {
return linear_intersection_detail::intersect(
linear_intersection_detail::parametric_object(first),
linear_intersection_detail::parametric_object(second),
eps
);
}
template <Coordinate T>
LinearIntersection linear_intersection(
const Line<T>& line,
const Segment<T>& segment,
long double eps = 1e-12L
) {
return linear_intersection_detail::intersect(
linear_intersection_detail::parametric_object(line),
linear_intersection_detail::parametric_object(segment),
eps
);
}
template <Coordinate T>
LinearIntersection linear_intersection(
const Segment<T>& segment,
const Line<T>& line,
long double eps = 1e-12L
) {
return linear_intersection_detail::intersect(
linear_intersection_detail::parametric_object(segment),
linear_intersection_detail::parametric_object(line),
eps
);
}
template <Coordinate T>
LinearIntersection linear_intersection(
const Segment<T>& first,
const Segment<T>& second,
long double eps = 1e-12L
) {
return linear_intersection_detail::intersect(
linear_intersection_detail::parametric_object(first),
linear_intersection_detail::parametric_object(second),
eps
);
}
template <Coordinate T>
LinearIntersection linear_intersection(
const Ray<T>& ray,
const Line<T>& line,
long double eps = 1e-12L
) {
return linear_intersection_detail::intersect(
linear_intersection_detail::parametric_object(ray),
linear_intersection_detail::parametric_object(line),
eps
);
}
template <Coordinate T>
LinearIntersection linear_intersection(
const Line<T>& line,
const Ray<T>& ray,
long double eps = 1e-12L
) {
return linear_intersection_detail::intersect(
linear_intersection_detail::parametric_object(line),
linear_intersection_detail::parametric_object(ray),
eps
);
}
template <Coordinate T>
LinearIntersection linear_intersection(
const Ray<T>& ray,
const Segment<T>& segment,
long double eps = 1e-12L
) {
return linear_intersection_detail::intersect(
linear_intersection_detail::parametric_object(ray),
linear_intersection_detail::parametric_object(segment),
eps
);
}
template <Coordinate T>
LinearIntersection linear_intersection(
const Segment<T>& segment,
const Ray<T>& ray,
long double eps = 1e-12L
) {
return linear_intersection_detail::intersect(
linear_intersection_detail::parametric_object(segment),
linear_intersection_detail::parametric_object(ray),
eps
);
}
template <Coordinate T>
LinearIntersection linear_intersection(
const Ray<T>& first,
const Ray<T>& second,
long double eps = 1e-12L
) {
return linear_intersection_detail::intersect(
linear_intersection_detail::parametric_object(first),
linear_intersection_detail::parametric_object(second),
eps
);
}
namespace closest_points_detail {
inline ClosestPoints reversed(const ClosestPoints& result) {
return ClosestPoints{result.second, result.first};
}
inline bool point_less(
const Point<long double>& first,
const Point<long double>& second
) {
if (first.x != second.x) return first.x < second.x;
return first.y < second.y;
}
inline ClosestPoints common_point(const LinearIntersection& intersection) {
assert(intersection.kind != LinearIntersectionKind::Empty);
Point<long double> point = intersection.first;
if (intersection.kind == LinearIntersectionKind::Segment) {
if (point_less(intersection.second, point)) {
point = intersection.second;
}
} else if (intersection.kind == LinearIntersectionKind::Line) {
const Line<long double> line{
intersection.first,
intersection.second
};
point = projection(line, Point<long double>(0, 0));
}
return ClosestPoints{point, point};
}
inline long double separation2(const ClosestPoints& result) {
return distance2(result.first, result.second);
}
inline bool canonical_less(
const ClosestPoints& first,
const ClosestPoints& second
) {
Point<long double> first_start = first.first;
Point<long double> first_finish = first.second;
if (point_less(first_finish, first_start)) {
std::swap(first_start, first_finish);
}
Point<long double> second_start = second.first;
Point<long double> second_finish = second.second;
if (point_less(second_finish, second_start)) {
std::swap(second_start, second_finish);
}
if (point_less(first_start, second_start)) return true;
if (point_less(second_start, first_start)) return false;
return point_less(first_finish, second_finish);
}
inline void consider(ClosestPoints& best, const ClosestPoints& candidate) {
const long double best_distance = separation2(best);
const long double candidate_distance = separation2(candidate);
if (
candidate_distance < best_distance ||
(
candidate_distance == best_distance &&
canonical_less(candidate, best)
)
) {
best = candidate;
}
}
} // namespace closest_points_detail
template <Coordinate T>
ClosestPoints closest_points(
const Point<T>& first,
const Point<T>& second
) {
return ClosestPoints{
Point<long double>(first),
Point<long double>(second),
};
}
template <Coordinate T>
ClosestPoints closest_points(
const Line<T>& line,
const Point<T>& point
) {
return ClosestPoints{
projection(line, point),
Point<long double>(point),
};
}
template <Coordinate T>
ClosestPoints closest_points(
const Point<T>& point,
const Line<T>& line
) {
return closest_points_detail::reversed(closest_points(line, point));
}
template <Coordinate T>
ClosestPoints closest_points(
const Segment<T>& segment,
const Point<T>& point
) {
return ClosestPoints{
projection(segment, point),
Point<long double>(point),
};
}
template <Coordinate T>
ClosestPoints closest_points(
const Point<T>& point,
const Segment<T>& segment
) {
return closest_points_detail::reversed(closest_points(segment, point));
}
template <Coordinate T>
ClosestPoints closest_points(
const Ray<T>& ray,
const Point<T>& point
) {
return ClosestPoints{
projection(ray, point),
Point<long double>(point),
};
}
template <Coordinate T>
ClosestPoints closest_points(
const Point<T>& point,
const Ray<T>& ray
) {
return closest_points_detail::reversed(closest_points(ray, point));
}
template <Coordinate T>
ClosestPoints closest_points(
const Line<T>& first,
const Line<T>& second,
long double eps = 1e-12L
) {
const LinearIntersection intersection =
linear_intersection(first, second, eps);
if (intersection.kind != LinearIntersectionKind::Empty) {
return closest_points_detail::common_point(intersection);
}
ClosestPoints result = closest_points(first, second.a);
closest_points_detail::consider(
result,
closest_points(first.a, second)
);
return result;
}
template <Coordinate T>
ClosestPoints closest_points(
const Line<T>& line,
const Segment<T>& segment,
long double eps = 1e-12L
) {
const LinearIntersection intersection =
linear_intersection(line, segment, eps);
if (intersection.kind != LinearIntersectionKind::Empty) {
return closest_points_detail::common_point(intersection);
}
ClosestPoints result = closest_points(line, segment.a);
closest_points_detail::consider(
result,
closest_points(line, segment.b)
);
return result;
}
template <Coordinate T>
ClosestPoints closest_points(
const Segment<T>& segment,
const Line<T>& line,
long double eps = 1e-12L
) {
return closest_points_detail::reversed(
closest_points(line, segment, eps)
);
}
template <Coordinate T>
ClosestPoints closest_points(
const Segment<T>& first,
const Segment<T>& second,
long double eps = 1e-12L
) {
const LinearIntersection intersection =
linear_intersection(first, second, eps);
if (intersection.kind != LinearIntersectionKind::Empty) {
return closest_points_detail::common_point(intersection);
}
ClosestPoints result = closest_points(first, second.a);
closest_points_detail::consider(
result,
closest_points(first, second.b)
);
closest_points_detail::consider(
result,
closest_points(first.a, second)
);
closest_points_detail::consider(
result,
closest_points(first.b, second)
);
return result;
}
template <Coordinate T>
ClosestPoints closest_points(
const Line<T>& line,
const Ray<T>& ray,
long double eps = 1e-12L
) {
const LinearIntersection intersection =
linear_intersection(line, ray, eps);
if (intersection.kind != LinearIntersectionKind::Empty) {
return closest_points_detail::common_point(intersection);
}
return closest_points(line, ray.origin);
}
template <Coordinate T>
ClosestPoints closest_points(
const Ray<T>& ray,
const Line<T>& line,
long double eps = 1e-12L
) {
return closest_points_detail::reversed(closest_points(line, ray, eps));
}
template <Coordinate T>
ClosestPoints closest_points(
const Ray<T>& ray,
const Segment<T>& segment,
long double eps = 1e-12L
) {
const LinearIntersection intersection =
linear_intersection(ray, segment, eps);
if (intersection.kind != LinearIntersectionKind::Empty) {
return closest_points_detail::common_point(intersection);
}
ClosestPoints result = closest_points(ray, segment.a);
closest_points_detail::consider(
result,
closest_points(ray, segment.b)
);
closest_points_detail::consider(
result,
closest_points(ray.origin, segment)
);
return result;
}
template <Coordinate T>
ClosestPoints closest_points(
const Segment<T>& segment,
const Ray<T>& ray,
long double eps = 1e-12L
) {
return closest_points_detail::reversed(
closest_points(ray, segment, eps)
);
}
template <Coordinate T>
ClosestPoints closest_points(
const Ray<T>& first,
const Ray<T>& second,
long double eps = 1e-12L
) {
const LinearIntersection intersection =
linear_intersection(first, second, eps);
if (intersection.kind != LinearIntersectionKind::Empty) {
return closest_points_detail::common_point(intersection);
}
ClosestPoints result = closest_points(first, second.origin);
closest_points_detail::consider(
result,
closest_points(first.origin, second)
);
return result;
}
template <Coordinate T>
long double distance(const Line<T>& line, const Point<T>& point) {
const ClosestPoints result = closest_points(line, point);
return geometry::distance(result.first, result.second);
}
template <Coordinate T>
long double distance(const Point<T>& point, const Line<T>& line) {
return distance(line, point);
}
template <Coordinate T>
long double distance(const Segment<T>& segment, const Point<T>& point) {
const ClosestPoints result = closest_points(segment, point);
return geometry::distance(result.first, result.second);
}
template <Coordinate T>
long double distance(const Point<T>& point, const Segment<T>& segment) {
return distance(segment, point);
}
template <Coordinate T>
long double distance(const Ray<T>& ray, const Point<T>& point) {
const ClosestPoints result = closest_points(ray, point);
return geometry::distance(result.first, result.second);
}
template <Coordinate T>
long double distance(const Point<T>& point, const Ray<T>& ray) {
return distance(ray, point);
}
template <Coordinate T>
long double distance(const Line<T>& first, const Line<T>& second) {
const ClosestPoints result = closest_points(first, second);
return geometry::distance(result.first, result.second);
}
template <Coordinate T>
long double distance(const Line<T>& line, const Segment<T>& segment) {
const ClosestPoints result = closest_points(line, segment);
return geometry::distance(result.first, result.second);
}
template <Coordinate T>
long double distance(const Segment<T>& segment, const Line<T>& line) {
return distance(line, segment);
}
template <Coordinate T>
long double distance(const Segment<T>& first, const Segment<T>& second) {
const ClosestPoints result = closest_points(first, second);
return geometry::distance(result.first, result.second);
}
template <Coordinate T>
long double distance(const Line<T>& line, const Ray<T>& ray) {
const ClosestPoints result = closest_points(line, ray);
return geometry::distance(result.first, result.second);
}
template <Coordinate T>
long double distance(const Ray<T>& ray, const Line<T>& line) {
return distance(line, ray);
}
template <Coordinate T>
long double distance(const Ray<T>& ray, const Segment<T>& segment) {
const ClosestPoints result = closest_points(ray, segment);
return geometry::distance(result.first, result.second);
}
template <Coordinate T>
long double distance(const Segment<T>& segment, const Ray<T>& ray) {
return distance(ray, segment);
}
template <Coordinate T>
long double distance(const Ray<T>& first, const Ray<T>& second) {
const ClosestPoints result = closest_points(first, second);
return geometry::distance(result.first, result.second);
}
} // namespace geometry
} // namespace m1une
#line 15 "geometry/circle.hpp"
namespace m1une {
namespace geometry {
template <Coordinate T>
struct Circle {
Point<T> center;
T radius;
bool filled = true;
};
enum class PointInCircle {
Outside = 0,
Boundary = 1,
Inside = 2,
};
enum class CircleRelation {
Separate,
ExternallyTangent,
Intersecting,
InternallyTangent,
Contained,
Coincident,
};
enum class AngularCoverageKind {
Empty,
Point,
Arc,
Full,
};
struct AngularCoverage {
AngularCoverageKind kind = AngularCoverageKind::Empty;
long double begin = 0.0L;
long double end = 0.0L;
};
struct CircleContact {
Point<long double> point;
long double first_argument = 0.0L;
long double second_argument = 0.0L;
};
enum class CircleContactKind {
Empty,
Point,
TwoPoints,
Coincident,
};
struct CircleCircleIntersection {
CircleRelation relation = CircleRelation::Separate;
CircleContactKind contact_kind = CircleContactKind::Empty;
std::array<CircleContact, 2> contacts;
AngularCoverage first_inside_second;
AngularCoverage second_inside_first;
constexpr int contact_count() const noexcept {
if (contact_kind == CircleContactKind::Point) return 1;
if (contact_kind == CircleContactKind::TwoPoints) return 2;
return 0;
}
};
struct CircleLinearContact {
Point<long double> point;
long double circle_argument = 0.0L;
long double linear_parameter = 0.0L;
};
struct CircleLinearIntersection {
int contact_count = 0;
std::array<CircleLinearContact, 2> contacts;
};
namespace circle_detail {
inline int compare(long double first, long double second, long double eps) {
if (first < second - eps) return -1;
if (first > second + eps) return 1;
return 0;
}
inline bool close(
const Point<long double>& first,
const Point<long double>& second,
long double eps
) {
return geometry::distance(first, second) <= eps;
}
inline void push_unique(
std::vector<Point<long double>>& points,
const Point<long double>& point,
long double eps
) {
for (const Point<long double>& existing : points) {
if (close(existing, point, eps)) return;
}
points.push_back(point);
}
inline bool same_line(
const Line<long double>& first,
const Line<long double>& second,
long double eps
) {
Point<long double> first_direction = first.b - first.a;
Point<long double> second_direction = second.b - second.a;
if (std::fabs(cross(first_direction, second_direction)) > eps) {
return false;
}
return std::fabs(cross(first_direction, second.a - first.a)) <= eps;
}
inline Line<long double> tangent_line(
const Point<long double>& contact,
Point<long double> normal,
long double eps
) {
Point<long double> direction(-normal.y, normal.x);
if (
direction.x < -eps ||
(std::fabs(direction.x) <= eps && direction.y < 0)
) {
direction = -direction;
}
return Line<long double>{contact, contact + direction};
}
inline long double circular_segment_angle_term(
long double angle,
long double sine,
long double cosine
) {
if (angle >= 0.01L) return angle - sine * cosine;
const long double squared = angle * angle;
return angle * squared * (
2.0L / 3.0L +
squared * (
-2.0L / 15.0L +
squared * (4.0L / 315.0L - squared * 2.0L / 2835.0L)
)
);
}
inline long double segment_disk_signed_area(
const Point<long double>& first,
const Point<long double>& second,
long double radius,
long double eps
) {
const Point<long double> direction = second - first;
const long double quadratic = dot(direction, direction);
if (quadratic == 0.0L || radius == 0.0L) return 0.0L;
std::vector<long double> cuts = {0.0L, 1.0L};
const long double linear = 2.0L * dot(first, direction);
const long double constant = dot(first, first) - radius * radius;
const long double discriminant =
linear * linear - 4.0L * quadratic * constant;
const long double tolerance = eps * std::max({
1.0L,
std::fabs(linear * linear),
std::fabs(4.0L * quadratic * constant)
});
if (discriminant >= -tolerance) {
const long double root = std::sqrt(std::max(0.0L, discriminant));
const long double first_ratio =
(-linear - root) / (2.0L * quadratic);
const long double second_ratio =
(-linear + root) / (2.0L * quadratic);
if (eps < first_ratio && first_ratio < 1.0L - eps) {
cuts.push_back(first_ratio);
}
if (eps < second_ratio && second_ratio < 1.0L - eps) {
cuts.push_back(second_ratio);
}
}
std::sort(cuts.begin(), cuts.end());
cuts.erase(
std::unique(
cuts.begin(),
cuts.end(),
[eps](long double left, long double right) {
return std::fabs(left - right) <= eps;
}
),
cuts.end()
);
long double result = 0.0L;
for (std::size_t index = 1; index < cuts.size(); ++index) {
const long double left = cuts[index - 1];
const long double right = cuts[index];
const Point<long double> a = first + direction * left;
const Point<long double> b = first + direction * right;
const Point<long double> middle =
first + direction * ((left + right) / 2.0L);
if (norm(middle) <= radius + eps) {
result += cross(a, b) / 2.0L;
} else {
result +=
radius * radius * std::atan2(cross(a, b), dot(a, b)) /
2.0L;
}
}
return result;
}
} // namespace circle_detail
template <Coordinate T>
constexpr Point<long double> centroid(const Circle<T>& circle) {
assert(circle.radius >= 0);
return Point<long double>(circle.center);
}
template <Coordinate T>
constexpr long double circle_circumference(const Circle<T>& circle) {
assert(circle.radius >= 0);
return
2.0L * std::numbers::pi_v<long double> *
static_cast<long double>(circle.radius);
}
template <Coordinate T>
constexpr long double circle_area(const Circle<T>& circle) {
assert(circle.radius >= 0);
const long double radius = static_cast<long double>(circle.radius);
return std::numbers::pi_v<long double> * radius * radius;
}
inline long double normalize_circle_argument(long double argument) {
const long double full = 2.0L * std::numbers::pi_v<long double>;
argument = std::fmod(argument, full);
if (argument < 0.0L) argument += full;
if (argument == full) argument = 0.0L;
return argument;
}
template <Coordinate T>
Point<long double> circle_point_at(
const Circle<T>& circle,
long double argument
) {
assert(circle.radius >= 0);
const long double radius = static_cast<long double>(circle.radius);
return Point<long double>(circle.center) + Point<long double>(
radius * std::cos(argument),
radius * std::sin(argument)
);
}
inline long double angular_measure(const AngularCoverage& coverage) {
if (
coverage.kind == AngularCoverageKind::Empty ||
coverage.kind == AngularCoverageKind::Point
) {
return 0.0L;
}
if (coverage.kind == AngularCoverageKind::Full) {
return 2.0L * std::numbers::pi_v<long double>;
}
assert(coverage.kind == AngularCoverageKind::Arc);
assert(coverage.begin <= coverage.end);
return coverage.end - coverage.begin;
}
template <Coordinate T>
long double circle_arc_length(
const Circle<T>& circle,
const AngularCoverage& coverage
) {
assert(circle.radius >= 0);
return static_cast<long double>(circle.radius) * angular_measure(coverage);
}
template <Coordinate C, Coordinate P>
PointInCircle point_in_circle(
const Circle<C>& circle,
const Point<P>& point,
long double eps = 1e-12L
) {
assert(circle.radius >= 0);
assert(eps >= 0.0L);
if constexpr (ExactCoordinate<C> && ExactCoordinate<P>) {
using W = std::common_type_t<wide_type<C>, wide_type<P>>;
const W dx = W(point.x) - W(circle.center.x);
const W dy = W(point.y) - W(circle.center.y);
const W radius = W(circle.radius);
const W squared_distance = dx * dx + dy * dy;
const W squared_radius = radius * radius;
if (squared_distance < squared_radius) return PointInCircle::Inside;
if (squared_distance > squared_radius) return PointInCircle::Outside;
return PointInCircle::Boundary;
} else {
const int relation = circle_detail::compare(
geometry::distance(
Point<long double>(circle.center),
Point<long double>(point)
),
static_cast<long double>(circle.radius),
eps
);
if (relation < 0) return PointInCircle::Inside;
if (relation > 0) return PointInCircle::Outside;
return PointInCircle::Boundary;
}
}
template <Coordinate C, Coordinate P>
bool contains(
const Circle<C>& circle,
const Point<P>& point,
long double eps = 1e-12L
) {
const PointInCircle relation = point_in_circle(circle, point, eps);
return circle.filled
? relation != PointInCircle::Outside
: relation == PointInCircle::Boundary;
}
template <Coordinate C, Coordinate P>
bool on_circle(
const Circle<C>& circle,
const Point<P>& point,
long double eps = 1e-12L
) {
assert(circle.radius >= 0);
assert(eps >= 0.0L);
if constexpr (ExactCoordinate<C> && ExactCoordinate<P>) {
using W = std::common_type_t<wide_type<C>, wide_type<P>>;
const W dx = W(point.x) - W(circle.center.x);
const W dy = W(point.y) - W(circle.center.y);
const W radius = W(circle.radius);
return dx * dx + dy * dy == radius * radius;
} else {
return circle_detail::compare(
geometry::distance(
Point<long double>(circle.center),
Point<long double>(point)
),
static_cast<long double>(circle.radius),
eps
) == 0;
}
}
template <Coordinate C, Coordinate P>
long double circle_argument(
const Circle<C>& circle,
const Point<P>& point
) {
assert(circle.radius >= 0);
return normalize_circle_argument(std::atan2(
static_cast<long double>(point.y) -
static_cast<long double>(circle.center.y),
static_cast<long double>(point.x) -
static_cast<long double>(circle.center.x)
));
}
template <Coordinate C, Coordinate P>
bool intersects(
const Circle<C>& circle,
const Point<P>& point,
long double eps = 1e-12L
) {
return contains(circle, point, eps);
}
template <Coordinate P, Coordinate C>
bool intersects(
const Point<P>& point,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return intersects(circle, point, eps);
}
template <Coordinate A, Coordinate B>
Circle<long double> circle_from_diameter(
const Point<A>& first,
const Point<B>& second
) {
Point<long double> a(first);
Point<long double> b(second);
Point<long double> center = (a + b) / 2.0L;
return Circle<long double>{center, geometry::distance(a, b) / 2.0L};
}
template <Coordinate T>
std::optional<Circle<long double>> incircle(
const Point<T>& first,
const Point<T>& second,
const Point<T>& third,
long double eps = 1e-12L
) {
assert(eps >= 0.0L);
if (orientation(first, second, third, eps) == 0) return std::nullopt;
long double opposite_first = geometry::distance(second, third);
long double opposite_second = geometry::distance(third, first);
long double opposite_third = geometry::distance(first, second);
long double perimeter =
opposite_first + opposite_second + opposite_third;
Point<long double> center =
(Point<long double>(first) * opposite_first +
Point<long double>(second) * opposite_second +
Point<long double>(third) * opposite_third) /
perimeter;
long double doubled_area = std::fabs(
static_cast<long double>(cross(first, second, third))
);
return Circle<long double>{center, doubled_area / perimeter};
}
template <Coordinate T>
std::optional<Circle<long double>> circumcircle(
const Point<T>& first,
const Point<T>& second,
const Point<T>& third,
long double eps = 1e-12L
) {
assert(eps >= 0.0L);
if (orientation(first, second, third, eps) == 0) return std::nullopt;
Point<long double> origin(first);
Point<long double> u = Point<long double>(second) - origin;
Point<long double> v = Point<long double>(third) - origin;
long double denominator = 2.0L * cross(u, v);
long double u_norm = norm2(u);
long double v_norm = norm2(v);
Point<long double> offset(
(u_norm * v.y - v_norm * u.y) / denominator,
(u.x * v_norm - v.x * u_norm) / denominator
);
Point<long double> center = origin + offset;
return Circle<long double>{center, norm(offset)};
}
template <Coordinate A, Coordinate B>
CircleRelation circle_relation(
const Circle<A>& first,
const Circle<B>& second,
long double eps = 1e-12L
) {
assert(first.radius >= 0);
assert(second.radius >= 0);
assert(eps >= 0.0L);
if constexpr (ExactCoordinate<A> && ExactCoordinate<B>) {
using W = std::common_type_t<wide_type<A>, wide_type<B>>;
W dx = W(second.center.x) - W(first.center.x);
W dy = W(second.center.y) - W(first.center.y);
W squared_distance = dx * dx + dy * dy;
W first_radius = W(first.radius);
W second_radius = W(second.radius);
W sum = first_radius + second_radius;
W difference = first_radius - second_radius;
if (difference < 0) difference = -difference;
if (squared_distance == 0 && difference == 0) {
return CircleRelation::Coincident;
}
if (squared_distance > sum * sum) return CircleRelation::Separate;
if (squared_distance == sum * sum) {
return CircleRelation::ExternallyTangent;
}
if (squared_distance < difference * difference) {
return CircleRelation::Contained;
}
if (squared_distance == difference * difference) {
return CircleRelation::InternallyTangent;
}
return CircleRelation::Intersecting;
} else {
long double center_distance = geometry::distance(
Point<long double>(first.center),
Point<long double>(second.center)
);
long double first_radius = static_cast<long double>(first.radius);
long double second_radius = static_cast<long double>(second.radius);
long double sum = first_radius + second_radius;
long double difference = std::fabs(first_radius - second_radius);
if (
center_distance <= eps &&
difference <= eps
) {
return CircleRelation::Coincident;
}
int outer = circle_detail::compare(center_distance, sum, eps);
if (outer > 0) return CircleRelation::Separate;
if (outer == 0) return CircleRelation::ExternallyTangent;
int inner = circle_detail::compare(center_distance, difference, eps);
if (inner < 0) return CircleRelation::Contained;
if (inner == 0) return CircleRelation::InternallyTangent;
return CircleRelation::Intersecting;
}
}
template <Coordinate C, Coordinate L>
CircleLinearIntersection circle_boundary_intersection(
const Circle<C>& circle,
const Line<L>& line,
long double eps = 1e-12L
) {
assert(circle.radius >= 0);
assert(line.a != line.b);
assert(eps >= 0.0L);
const Point<long double> center(circle.center);
const Point<long double> origin(line.a);
const Point<long double> direction =
Point<long double>(line.b) - origin;
const long double squared_length = dot(direction, direction);
const long double length = std::sqrt(squared_length);
const long double foot_parameter =
dot(center - origin, direction) / squared_length;
const Point<long double> foot =
origin + direction * foot_parameter;
const long double distance_to_line = geometry::distance(center, foot);
const long double radius = static_cast<long double>(circle.radius);
const int relation =
circle_detail::compare(distance_to_line, radius, eps);
CircleLinearIntersection result;
if (relation > 0) return result;
if (relation == 0) {
result.contact_count = 1;
result.contacts[0] = CircleLinearContact{
foot,
circle_argument(circle, foot),
foot_parameter
};
return result;
}
const long double offset = std::sqrt(std::max(
0.0L,
radius * radius - distance_to_line * distance_to_line
));
const long double parameter_offset = offset / length;
result.contact_count = 2;
for (int index = 0; index < 2; ++index) {
const long double parameter = foot_parameter +
(index == 0 ? -parameter_offset : parameter_offset);
const Point<long double> point = origin + direction * parameter;
result.contacts[index] = CircleLinearContact{
point,
circle_argument(circle, point),
parameter
};
}
return result;
}
template <Coordinate L, Coordinate C>
CircleLinearIntersection circle_boundary_intersection(
const Line<L>& line,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return circle_boundary_intersection(circle, line, eps);
}
template <Coordinate C, Coordinate R>
CircleLinearIntersection circle_boundary_intersection(
const Circle<C>& circle,
const Ray<R>& ray,
long double eps = 1e-12L
) {
assert(ray.origin != ray.through);
const Line<R> line{ray.origin, ray.through};
const CircleLinearIntersection line_result =
circle_boundary_intersection(circle, line, eps);
CircleLinearIntersection result;
for (int index = 0; index < line_result.contact_count; ++index) {
CircleLinearContact contact = line_result.contacts[index];
if (contact.linear_parameter < -eps) continue;
if (std::fabs(contact.linear_parameter) <= eps) {
contact.linear_parameter = 0.0L;
contact.point = Point<long double>(ray.origin);
contact.circle_argument = circle_argument(circle, contact.point);
}
result.contacts[result.contact_count++] = contact;
}
return result;
}
template <Coordinate R, Coordinate C>
CircleLinearIntersection circle_boundary_intersection(
const Ray<R>& ray,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return circle_boundary_intersection(circle, ray, eps);
}
template <Coordinate C, Coordinate S>
CircleLinearIntersection circle_boundary_intersection(
const Circle<C>& circle,
const Segment<S>& segment,
long double eps = 1e-12L
) {
assert(circle.radius >= 0);
assert(eps >= 0.0L);
CircleLinearIntersection result;
if (segment.a == segment.b) {
if (on_circle(circle, segment.a, eps)) {
const Point<long double> point(segment.a);
result.contact_count = 1;
result.contacts[0] = CircleLinearContact{
point,
circle_argument(circle, point),
0.0L
};
}
return result;
}
const Line<S> line{segment.a, segment.b};
const CircleLinearIntersection line_result =
circle_boundary_intersection(circle, line, eps);
for (int index = 0; index < line_result.contact_count; ++index) {
CircleLinearContact contact = line_result.contacts[index];
if (
contact.linear_parameter < -eps ||
contact.linear_parameter > 1.0L + eps
) {
continue;
}
if (std::fabs(contact.linear_parameter) <= eps) {
contact.linear_parameter = 0.0L;
contact.point = Point<long double>(segment.a);
} else if (std::fabs(contact.linear_parameter - 1.0L) <= eps) {
contact.linear_parameter = 1.0L;
contact.point = Point<long double>(segment.b);
}
contact.circle_argument = circle_argument(circle, contact.point);
result.contacts[result.contact_count++] = contact;
}
return result;
}
template <Coordinate S, Coordinate C>
CircleLinearIntersection circle_boundary_intersection(
const Segment<S>& segment,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return circle_boundary_intersection(circle, segment, eps);
}
template <Coordinate A, Coordinate B>
CircleCircleIntersection circle_boundary_intersection(
const Circle<A>& first,
const Circle<B>& second,
long double eps = 1e-12L
) {
assert(first.radius >= 0);
assert(second.radius >= 0);
assert(eps >= 0.0L);
const long double full = 2.0L * std::numbers::pi_v<long double>;
const long double first_radius = static_cast<long double>(first.radius);
const long double second_radius = static_cast<long double>(second.radius);
CircleCircleIntersection result;
result.relation = circle_relation(first, second, eps);
auto point_coverage = [](long double argument) {
return AngularCoverage{
AngularCoverageKind::Point,
argument,
argument
};
};
auto full_coverage = [full]() {
return AngularCoverage{AngularCoverageKind::Full, 0.0L, full};
};
if (result.relation == CircleRelation::Coincident) {
if (first_radius == 0.0L) {
const Point<long double> point(first.center);
result.contact_kind = CircleContactKind::Point;
result.contacts[0] = CircleContact{point, 0.0L, 0.0L};
result.first_inside_second = point_coverage(0.0L);
result.second_inside_first = point_coverage(0.0L);
} else {
result.contact_kind = CircleContactKind::Coincident;
result.first_inside_second = full_coverage();
result.second_inside_first = full_coverage();
}
return result;
}
if (result.relation == CircleRelation::Separate) return result;
if (result.relation == CircleRelation::Contained) {
if (first_radius < second_radius) {
result.first_inside_second = first_radius == 0.0L
? point_coverage(0.0L)
: full_coverage();
} else {
result.second_inside_first = second_radius == 0.0L
? point_coverage(0.0L)
: full_coverage();
}
return result;
}
const Point<long double> first_center(first.center);
const Point<long double> second_center(second.center);
const Point<long double> center_direction = second_center - first_center;
const long double center_distance = norm(center_direction);
const Point<long double> unit = center_direction / center_distance;
const long double along =
(first_radius * first_radius - second_radius * second_radius +
center_distance * center_distance) /
(2.0L * center_distance);
const Point<long double> base = first_center + unit * along;
if (
result.relation == CircleRelation::ExternallyTangent ||
result.relation == CircleRelation::InternallyTangent
) {
const long double first_argument =
circle_argument(first, base);
const long double second_argument =
circle_argument(second, base);
result.contact_kind = CircleContactKind::Point;
result.contacts[0] = CircleContact{
base,
first_argument,
second_argument
};
result.first_inside_second = point_coverage(first_argument);
result.second_inside_first = point_coverage(second_argument);
if (result.relation == CircleRelation::InternallyTangent) {
if (first_radius < second_radius && first_radius > 0.0L) {
result.first_inside_second = full_coverage();
} else if (
second_radius < first_radius && second_radius > 0.0L
) {
result.second_inside_first = full_coverage();
}
}
return result;
}
assert(result.relation == CircleRelation::Intersecting);
const long double height = std::sqrt(std::max(
0.0L,
first_radius * first_radius - along * along
));
const Point<long double> perpendicular(-unit.y, unit.x);
const Point<long double> first_point = base - perpendicular * height;
const Point<long double> second_point = base + perpendicular * height;
result.contact_kind = CircleContactKind::TwoPoints;
result.contacts[0] = CircleContact{
first_point,
circle_argument(first, first_point),
circle_argument(second, first_point)
};
result.contacts[1] = CircleContact{
second_point,
circle_argument(first, second_point),
circle_argument(second, second_point)
};
const long double first_begin = result.contacts[0].first_argument;
long double first_end = result.contacts[1].first_argument;
if (first_end <= first_begin) first_end += full;
result.first_inside_second = AngularCoverage{
AngularCoverageKind::Arc,
first_begin,
first_end
};
const long double second_begin = result.contacts[1].second_argument;
long double second_end = result.contacts[0].second_argument;
if (second_end <= second_begin) second_end += full;
result.second_inside_first = AngularCoverage{
AngularCoverageKind::Arc,
second_begin,
second_end
};
return result;
}
template <Coordinate C, Coordinate L>
bool intersects(
const Circle<C>& circle,
const Line<L>& line,
long double eps = 1e-12L
) {
if (circle.filled) {
const Line<long double> converted{
Point<long double>(line.a),
Point<long double>(line.b)
};
return contains(
circle,
projection(converted, Point<long double>(circle.center)),
eps
);
}
return circle_boundary_intersection(circle, line, eps).contact_count > 0;
}
template <Coordinate C, Coordinate L>
bool intersects(
const Line<L>& line,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return intersects(circle, line, eps);
}
template <Coordinate C, Coordinate R>
bool intersects(
const Circle<C>& circle,
const Ray<R>& ray,
long double eps = 1e-12L
) {
if (circle.filled) {
const Ray<long double> converted{
Point<long double>(ray.origin),
Point<long double>(ray.through)
};
return contains(
circle,
projection(converted, Point<long double>(circle.center)),
eps
);
}
return circle_boundary_intersection(circle, ray, eps).contact_count > 0;
}
template <Coordinate C, Coordinate R>
bool intersects(
const Ray<R>& ray,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return intersects(circle, ray, eps);
}
template <Coordinate C, Coordinate S>
bool intersects(
const Circle<C>& circle,
const Segment<S>& segment,
long double eps = 1e-12L
) {
if (circle.filled) {
const Segment<long double> converted{
Point<long double>(segment.a),
Point<long double>(segment.b)
};
return contains(
circle,
projection(converted, Point<long double>(circle.center)),
eps
);
}
return
circle_boundary_intersection(circle, segment, eps).contact_count > 0;
}
template <Coordinate C, Coordinate S>
bool intersects(
const Segment<S>& segment,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return intersects(circle, segment, eps);
}
template <Coordinate A, Coordinate B>
bool intersects(
const Circle<A>& first,
const Circle<B>& second,
long double eps = 1e-12L
) {
assert(first.radius >= 0);
assert(second.radius >= 0);
assert(eps >= 0.0L);
if (first.filled && second.filled) {
if constexpr (ExactCoordinate<A> && ExactCoordinate<B>) {
using W = std::common_type_t<wide_type<A>, wide_type<B>>;
const W dx = W(second.center.x) - W(first.center.x);
const W dy = W(second.center.y) - W(first.center.y);
const W radius = W(first.radius) + W(second.radius);
return dx * dx + dy * dy <= radius * radius;
} else {
const long double center_distance = geometry::distance(
Point<long double>(first.center),
Point<long double>(second.center)
);
const long double radius_sum =
static_cast<long double>(first.radius) +
static_cast<long double>(second.radius);
return circle_detail::compare(
center_distance,
radius_sum,
eps
) <= 0;
}
}
if (first.filled != second.filled) {
const long double center_distance = geometry::distance(
Point<long double>(first.center),
Point<long double>(second.center)
);
const long double boundary_radius = first.filled
? static_cast<long double>(second.radius)
: static_cast<long double>(first.radius);
const long double filled_radius = first.filled
? static_cast<long double>(first.radius)
: static_cast<long double>(second.radius);
return circle_detail::compare(
std::fabs(center_distance - boundary_radius),
filled_radius,
eps
) <= 0;
}
CircleRelation relation = circle_relation(first, second, eps);
return
relation == CircleRelation::ExternallyTangent ||
relation == CircleRelation::Intersecting ||
relation == CircleRelation::InternallyTangent ||
relation == CircleRelation::Coincident;
}
template <Coordinate R, Coordinate H, Coordinate C>
Ray<long double> reflected_ray(
const Ray<R>& incoming,
const Point<H>& hit,
const Circle<C>& circle,
long double eps = 1e-12L
) {
assert(incoming.origin != incoming.through);
assert(eps >= 0.0L);
assert(static_cast<long double>(circle.radius) > eps);
assert(
std::fabs(
geometry::distance(
Point<long double>(hit),
Point<long double>(circle.center)
) -
static_cast<long double>(circle.radius)
) <= eps
);
Point<long double> hit_point(hit);
Point<long double> normal = normalized(
hit_point - Point<long double>(circle.center)
);
Point<long double> incoming_direction =
Point<long double>(incoming.through) -
Point<long double>(incoming.origin);
Point<long double> outgoing_direction =
incoming_direction - normal * (2.0L * dot(incoming_direction, normal));
return Ray<long double>{hit_point, hit_point + outgoing_direction};
}
template <Coordinate C, Coordinate P>
std::vector<Point<long double>> tangent_points(
const Circle<C>& circle,
const Point<P>& point,
long double eps = 1e-12L
) {
assert(circle.radius >= 0);
assert(eps >= 0.0L);
Point<long double> center(circle.center);
Point<long double> external(point);
Point<long double> direction = external - center;
long double squared_distance = dot(direction, direction);
long double radius = static_cast<long double>(circle.radius);
if (radius == 0.0L) return {center};
long double center_distance = std::sqrt(squared_distance);
int relation = circle_detail::compare(center_distance, radius, eps);
if (relation < 0) return {};
if (relation == 0) {
return {center + direction * (radius / center_distance)};
}
Point<long double> base =
center + direction * (radius * radius / squared_distance);
long double scale =
radius * std::sqrt(std::max(
0.0L,
squared_distance - radius * radius
)) /
squared_distance;
Point<long double> perpendicular(-direction.y, direction.x);
Point<long double> first = base - perpendicular * scale;
Point<long double> second = base + perpendicular * scale;
if (second < first) std::swap(first, second);
return {first, second};
}
template <Coordinate A, Coordinate B>
std::vector<Line<long double>> common_tangents(
const Circle<A>& first,
const Circle<B>& second,
long double eps = 1e-12L
) {
assert(first.radius >= 0);
assert(second.radius >= 0);
assert(eps >= 0.0L);
Point<long double> first_center(first.center);
Point<long double> second_center(second.center);
Point<long double> direction = second_center - first_center;
long double squared_distance = dot(direction, direction);
long double center_distance = std::sqrt(squared_distance);
if (center_distance <= eps) return {};
long double first_radius = static_cast<long double>(first.radius);
long double second_radius = static_cast<long double>(second.radius);
std::vector<Line<long double>> result;
for (int second_side : {1, -1}) {
long double difference =
first_radius - second_side * second_radius;
int relation = circle_detail::compare(
std::fabs(difference),
center_distance,
eps
);
if (relation > 0) continue;
long double perpendicular_length = relation == 0 ? 0.0L : std::sqrt(
std::max(0.0L, squared_distance - difference * difference)
);
int choices = perpendicular_length <= eps ? 1 : 2;
for (int choice = 0; choice < choices; ++choice) {
long double side = choice == 0 ? -1.0L : 1.0L;
Point<long double> normal =
direction * (difference / squared_distance) +
Point<long double>(-direction.y, direction.x) *
(side * perpendicular_length / squared_distance);
normal = normalized(normal);
Point<long double> contact =
first_center + normal * first_radius;
Line<long double> tangent =
circle_detail::tangent_line(contact, normal, eps);
bool duplicate = false;
for (const Line<long double>& existing : result) {
if (circle_detail::same_line(existing, tangent, eps)) {
duplicate = true;
break;
}
}
if (!duplicate) result.push_back(tangent);
}
}
std::sort(
result.begin(),
result.end(),
[](const Line<long double>& left, const Line<long double>& right) {
if (left.a != right.a) return left.a < right.a;
return left.b < right.b;
}
);
return result;
}
template <Coordinate A, Coordinate B>
std::vector<Point<long double>> common_tangent_points(
const Circle<A>& first,
const Circle<B>& second,
long double eps = 1e-12L
) {
std::vector<Point<long double>> result;
for (const Line<long double>& line : common_tangents(first, second, eps)) {
circle_detail::push_unique(result, line.a, eps);
}
std::sort(result.begin(), result.end());
return result;
}
// These area functions use the enclosed disks, independent of `filled`.
template <Coordinate A, Coordinate B>
long double circle_circle_intersection_area(
const Circle<A>& first,
const Circle<B>& second,
long double eps = 1e-12L
) {
assert(first.radius >= 0);
assert(second.radius >= 0);
assert(eps >= 0.0L);
const long double first_radius = static_cast<long double>(first.radius);
const long double second_radius = static_cast<long double>(second.radius);
const CircleRelation relation = circle_relation(first, second, eps);
if (
relation == CircleRelation::Separate ||
relation == CircleRelation::ExternallyTangent
) {
return 0.0L;
}
if (
relation == CircleRelation::Contained ||
relation == CircleRelation::InternallyTangent ||
relation == CircleRelation::Coincident
) {
const long double radius = std::min(first_radius, second_radius);
return std::numbers::pi_v<long double> * radius * radius;
}
const long double center_distance = geometry::distance(
Point<long double>(first.center),
Point<long double>(second.center)
);
const long double first_cosine = std::clamp(
(
(center_distance - second_radius) *
(center_distance + second_radius) +
first_radius * first_radius
) / (2.0L * center_distance * first_radius),
-1.0L,
1.0L
);
const long double second_cosine = std::clamp(
(
(center_distance - first_radius) *
(center_distance + first_radius) +
second_radius * second_radius
) / (2.0L * center_distance * second_radius),
-1.0L,
1.0L
);
const long double radicand =
(-center_distance + first_radius + second_radius) *
(center_distance + first_radius - second_radius) *
(center_distance - first_radius + second_radius) *
(center_distance + first_radius + second_radius);
const long double height =
std::sqrt(std::max(0.0L, radicand)) / (2.0L * center_distance);
const long double first_sine =
std::clamp(height / first_radius, 0.0L, 1.0L);
const long double second_sine =
std::clamp(height / second_radius, 0.0L, 1.0L);
const long double first_angle = std::atan2(first_sine, first_cosine);
const long double second_angle = std::atan2(second_sine, second_cosine);
return
first_radius * first_radius *
circle_detail::circular_segment_angle_term(
first_angle,
first_sine,
first_cosine
) +
second_radius * second_radius *
circle_detail::circular_segment_angle_term(
second_angle,
second_sine,
second_cosine
);
}
template <Coordinate C, Coordinate P>
long double circle_polygon_intersection_area(
const Circle<C>& circle,
const std::vector<Point<P>>& polygon,
long double eps = 1e-12L
) {
assert(circle.radius >= 0);
assert(eps >= 0.0L);
if (polygon.empty() || circle.radius == 0) return 0.0L;
const Point<long double> center(circle.center);
const long double radius = static_cast<long double>(circle.radius);
long double result = 0.0L;
for (std::size_t index = 0; index < polygon.size(); ++index) {
const Point<long double> first =
Point<long double>(polygon[index]) - center;
const Point<long double> second =
Point<long double>(polygon[(index + 1) % polygon.size()]) - center;
result += circle_detail::segment_disk_signed_area(
first,
second,
radius,
eps
);
}
return std::fabs(result);
}
namespace circle_detail {
template <Coordinate T>
Point<long double> point_toward(
const Circle<T>& circle,
const Point<long double>& target
) {
assert(circle.radius >= 0);
const Point<long double> center(circle.center);
const long double radius = static_cast<long double>(circle.radius);
const Point<long double> direction = target - center;
const long double length = norm(direction);
if (length == 0.0L) {
return center + Point<long double>(-radius, 0.0L);
}
return center + direction * (radius / length);
}
inline void consider(
ClosestPoints& best,
const ClosestPoints& candidate
) {
closest_points_detail::consider(best, candidate);
}
} // namespace circle_detail
template <Coordinate C, Coordinate P>
ClosestPoints closest_points(
const Circle<C>& circle,
const Point<P>& point,
long double eps = 1e-12L
) {
const Point<long double> converted(point);
if (circle.filled && contains(circle, converted, eps)) {
return ClosestPoints{converted, converted};
}
return ClosestPoints{
circle_detail::point_toward(circle, converted),
converted
};
}
template <Coordinate P, Coordinate C>
ClosestPoints closest_points(
const Point<P>& point,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return closest_points_detail::reversed(
closest_points(circle, point, eps)
);
}
template <Coordinate C, Coordinate L>
ClosestPoints closest_points(
const Circle<C>& circle,
const Line<L>& line,
long double eps = 1e-12L
) {
const Line<long double> converted{
Point<long double>(line.a),
Point<long double>(line.b)
};
const Point<long double> point = projection(
converted,
Point<long double>(circle.center)
);
if (circle.filled) return closest_points(circle, point, eps);
const CircleLinearIntersection common =
circle_boundary_intersection(circle, line, eps);
if (common.contact_count > 0) {
return ClosestPoints{
common.contacts[0].point,
common.contacts[0].point
};
}
return ClosestPoints{
circle_detail::point_toward(circle, point),
point
};
}
template <Coordinate L, Coordinate C>
ClosestPoints closest_points(
const Line<L>& line,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return closest_points_detail::reversed(closest_points(circle, line, eps));
}
template <Coordinate C, Coordinate R>
ClosestPoints closest_points(
const Circle<C>& circle,
const Ray<R>& ray,
long double eps = 1e-12L
) {
const Ray<long double> converted{
Point<long double>(ray.origin),
Point<long double>(ray.through)
};
const Point<long double> point = projection(
converted,
Point<long double>(circle.center)
);
if (circle.filled) return closest_points(circle, point, eps);
const CircleLinearIntersection common =
circle_boundary_intersection(circle, ray, eps);
if (common.contact_count > 0) {
return ClosestPoints{
common.contacts[0].point,
common.contacts[0].point
};
}
return ClosestPoints{
circle_detail::point_toward(circle, point),
point
};
}
template <Coordinate R, Coordinate C>
ClosestPoints closest_points(
const Ray<R>& ray,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return closest_points_detail::reversed(closest_points(circle, ray, eps));
}
template <Coordinate C, Coordinate S>
ClosestPoints closest_points(
const Circle<C>& circle,
const Segment<S>& segment,
long double eps = 1e-12L
) {
const Segment<long double> converted{
Point<long double>(segment.a),
Point<long double>(segment.b)
};
const Point<long double> center(circle.center);
const Point<long double> projected = projection(converted, center);
if (circle.filled) return closest_points(circle, projected, eps);
const CircleLinearIntersection common =
circle_boundary_intersection(circle, segment, eps);
if (common.contact_count > 0) {
return ClosestPoints{
common.contacts[0].point,
common.contacts[0].point
};
}
ClosestPoints result{
circle_detail::point_toward(circle, projected),
projected
};
for (const Point<long double>& point : {converted.a, converted.b}) {
circle_detail::consider(
result,
ClosestPoints{circle_detail::point_toward(circle, point), point}
);
}
return result;
}
template <Coordinate S, Coordinate C>
ClosestPoints closest_points(
const Segment<S>& segment,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return closest_points_detail::reversed(
closest_points(circle, segment, eps)
);
}
template <Coordinate A, Coordinate B>
ClosestPoints closest_points(
const Circle<A>& first,
const Circle<B>& second,
long double eps = 1e-12L
) {
if (first.filled && !second.filled) {
return closest_points_detail::reversed(
closest_points(second, first, eps)
);
}
if (!first.filled && second.filled) {
const ClosestPoints center_result =
closest_points(first, second.center, eps);
if (contains(second, center_result.first, eps)) {
return ClosestPoints{center_result.first, center_result.first};
}
const ClosestPoints filled_result =
closest_points(second, center_result.first, eps);
return ClosestPoints{center_result.first, filled_result.first};
}
if (first.filled && second.filled) {
assert(first.radius >= 0);
assert(second.radius >= 0);
const Point<long double> first_center(first.center);
const Point<long double> second_center(second.center);
Point<long double> direction = second_center - first_center;
const long double center_distance = norm(direction);
if (center_distance == 0.0L) {
return ClosestPoints{first_center, first_center};
}
direction = direction / center_distance;
const long double first_radius =
static_cast<long double>(first.radius);
const long double second_radius =
static_cast<long double>(second.radius);
if (intersects(first, second, eps)) {
const long double left = std::max(
-first_radius,
center_distance - second_radius
);
const long double right = std::min(
first_radius,
center_distance + second_radius
);
const Point<long double> common =
first_center + direction * ((left + right) / 2.0L);
return ClosestPoints{common, common};
}
return ClosestPoints{
first_center + direction * first_radius,
second_center - direction * second_radius
};
}
const CircleCircleIntersection common =
circle_boundary_intersection(first, second, eps);
if (common.contact_count() > 0) {
return ClosestPoints{
common.contacts[0].point,
common.contacts[0].point
};
}
const Point<long double> first_center(first.center);
const Point<long double> second_center(second.center);
if (circle_relation(first, second, eps) == CircleRelation::Coincident) {
const Point<long double> point =
circle_detail::point_toward(first, first_center);
return ClosestPoints{point, point};
}
Point<long double> direction = second_center - first_center;
const long double center_distance = norm(direction);
if (center_distance == 0.0L) {
const Point<long double> first_point =
circle_detail::point_toward(first, first_center);
const Point<long double> second_point =
circle_detail::point_toward(second, second_center);
return ClosestPoints{first_point, second_point};
}
direction = direction / center_distance;
const long double first_radius = static_cast<long double>(first.radius);
const long double second_radius = static_cast<long double>(second.radius);
ClosestPoints result{
first_center + direction * first_radius,
second_center + direction * second_radius
};
for (const long double first_sign : {-1.0L, 1.0L}) {
for (const long double second_sign : {-1.0L, 1.0L}) {
circle_detail::consider(
result,
ClosestPoints{
first_center + direction * (first_sign * first_radius),
second_center + direction * (second_sign * second_radius)
}
);
}
}
return result;
}
template <Coordinate C, Coordinate P>
long double distance(const Circle<C>& circle, const Point<P>& point) {
const ClosestPoints result = closest_points(circle, point);
return geometry::distance(result.first, result.second);
}
template <Coordinate P, Coordinate C>
long double distance(const Point<P>& point, const Circle<C>& circle) {
return distance(circle, point);
}
template <Coordinate C, Coordinate L>
long double distance(const Circle<C>& circle, const Line<L>& line) {
const ClosestPoints result = closest_points(circle, line);
return geometry::distance(result.first, result.second);
}
template <Coordinate L, Coordinate C>
long double distance(const Line<L>& line, const Circle<C>& circle) {
return distance(circle, line);
}
template <Coordinate C, Coordinate R>
long double distance(const Circle<C>& circle, const Ray<R>& ray) {
const ClosestPoints result = closest_points(circle, ray);
return geometry::distance(result.first, result.second);
}
template <Coordinate R, Coordinate C>
long double distance(const Ray<R>& ray, const Circle<C>& circle) {
return distance(circle, ray);
}
template <Coordinate C, Coordinate S>
long double distance(const Circle<C>& circle, const Segment<S>& segment) {
const ClosestPoints result = closest_points(circle, segment);
return geometry::distance(result.first, result.second);
}
template <Coordinate S, Coordinate C>
long double distance(const Segment<S>& segment, const Circle<C>& circle) {
return distance(circle, segment);
}
template <Coordinate A, Coordinate B>
long double distance(const Circle<A>& first, const Circle<B>& second) {
const ClosestPoints result = closest_points(first, second);
return geometry::distance(result.first, result.second);
}
} // namespace geometry
} // namespace m1une
#line 16 "geometry/polygon.hpp"
namespace m1une {
namespace geometry {
enum class PointInPolygon {
Outside = 0,
Boundary = 1,
Inside = 2,
};
template <Coordinate T>
struct Polygon {
std::vector<Point<T>> vertices;
bool filled = true;
};
template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
Polygon<std::common_type_t<T, Scalar>> operator*(
const Polygon<T>& polygon,
Scalar scalar
) {
using Result = std::common_type_t<T, Scalar>;
Polygon<Result> scaled;
scaled.vertices.reserve(polygon.vertices.size());
for (const Point<T>& point : polygon.vertices) {
scaled.vertices.push_back(point * scalar);
}
scaled.filled = polygon.filled;
return scaled;
}
template <typename Scalar, Coordinate T>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
Polygon<std::common_type_t<T, Scalar>> operator*(
Scalar scalar,
const Polygon<T>& polygon
) {
return polygon * scalar;
}
struct ParameterInterval {
long double begin = 0.0L;
long double end = 0.0L;
};
template <Coordinate T>
constexpr Point<long double> centroid(
const std::array<Point<T>, 3>& triangle
) {
return Point<long double>(
(
static_cast<long double>(triangle[0].x) +
static_cast<long double>(triangle[1].x) +
static_cast<long double>(triangle[2].x)
) / 3,
(
static_cast<long double>(triangle[0].y) +
static_cast<long double>(triangle[1].y) +
static_cast<long double>(triangle[2].y)
) / 3
);
}
namespace polygon_detail {
template <Coordinate T>
std::vector<Point<T>> clean_polygon_vertices(
std::vector<Point<T>> polygon,
long double eps
) {
if (
polygon.size() >= 2 &&
polygon.front() == polygon.back()
) {
polygon.pop_back();
}
std::vector<Point<T>> deduplicated;
for (const Point<T>& point : polygon) {
if (deduplicated.empty() || deduplicated.back() != point) {
deduplicated.push_back(point);
}
}
if (
deduplicated.size() >= 2 &&
deduplicated.front() == deduplicated.back()
) {
deduplicated.pop_back();
}
bool changed = true;
while (changed && deduplicated.size() >= 3) {
changed = false;
std::vector<Point<T>> cleaned;
std::size_t size = deduplicated.size();
for (std::size_t index = 0; index < size; ++index) {
const Point<T>& previous =
deduplicated[(index + size - 1) % size];
const Point<T>& current = deduplicated[index];
const Point<T>& next =
deduplicated[(index + 1) % size];
if (
orientation(previous, current, next, eps) == 0 &&
sign<T>(dot(current - previous, next - current), eps) >= 0
) {
changed = true;
} else {
cleaned.push_back(current);
}
}
deduplicated = std::move(cleaned);
}
return deduplicated;
}
template <Coordinate T>
bool in_ccw_triangle(
const Point<T>& point,
const Point<T>& first,
const Point<T>& second,
const Point<T>& third,
long double eps
) {
return
orientation(first, second, point, eps) >= 0 &&
orientation(second, third, point, eps) >= 0 &&
orientation(third, first, point, eps) >= 0;
}
} // namespace polygon_detail
template <Coordinate T>
wide_type<T> polygon_area2(const std::vector<Point<T>>& polygon) {
wide_type<T> result = 0;
std::size_t n = polygon.size();
for (std::size_t i = 0; i < n; i++) {
result += cross(polygon[i], polygon[(i + 1) % n]);
}
return result;
}
template <Coordinate T>
long double polygon_area(const std::vector<Point<T>>& polygon) {
return std::fabs(static_cast<long double>(polygon_area2(polygon))) / 2;
}
template <Coordinate T>
std::optional<Point<long double>> polygon_centroid(
const std::vector<Point<T>>& polygon,
long double eps = 1e-12L
) {
if (polygon.size() < 3) return std::nullopt;
wide_type<T> signed_area2 = polygon_area2(polygon);
if (sign<T>(signed_area2, eps) == 0) return std::nullopt;
long double x_numerator = 0;
long double y_numerator = 0;
std::size_t size = polygon.size();
for (std::size_t index = 0; index < size; ++index) {
const Point<T>& current = polygon[index];
const Point<T>& next = polygon[(index + 1) % size];
long double weight = static_cast<long double>(cross(current, next));
x_numerator +=
(static_cast<long double>(current.x) +
static_cast<long double>(next.x)) *
weight;
y_numerator +=
(static_cast<long double>(current.y) +
static_cast<long double>(next.y)) *
weight;
}
long double denominator =
3.0L * static_cast<long double>(signed_area2);
return Point<long double>(
x_numerator / denominator,
y_numerator / denominator
);
}
template <Coordinate T>
std::optional<Point<long double>> centroid(
const std::vector<Point<T>>& polygon,
long double eps = 1e-12L
) {
return polygon_centroid(polygon, eps);
}
template <Coordinate T>
std::optional<Point<long double>> polygon_center_of_gravity(
const std::vector<Point<T>>& polygon,
long double eps = 1e-12L
) {
return polygon_centroid(polygon, eps);
}
template <Coordinate T>
bool is_simple_polygon(
const std::vector<Point<T>>& polygon,
long double eps = 1e-12L
) {
if (polygon.size() < 3) return false;
std::size_t size = polygon.size();
for (std::size_t index = 0; index < size; ++index) {
const Point<T>& previous = polygon[(index + size - 1) % size];
const Point<T>& current = polygon[index];
const Point<T>& next = polygon[(index + 1) % size];
if (current == next) return false;
if (
orientation(previous, current, next, eps) == 0 &&
sign<T>(dot(current - previous, next - current), eps) < 0
) {
return false;
}
}
for (std::size_t first_index = 0; first_index < size; ++first_index) {
Segment<T> first{
polygon[first_index],
polygon[(first_index + 1) % size]
};
for (
std::size_t second_index = first_index + 1;
second_index < size;
++second_index
) {
bool adjacent =
second_index == first_index + 1 ||
(first_index == 0 && second_index + 1 == size);
if (adjacent) continue;
Segment<T> second{
polygon[second_index],
polygon[(second_index + 1) % size]
};
if (intersects(first, second, eps)) return false;
}
}
return true;
}
template <Coordinate T>
std::optional<std::vector<std::array<Point<T>, 3>>> triangulate_polygon(
std::vector<Point<T>> polygon,
long double eps = 1e-12L
) {
polygon =
polygon_detail::clean_polygon_vertices(std::move(polygon), eps);
if (polygon.size() < 3) return std::nullopt;
wide_type<T> signed_area2 = polygon_area2(polygon);
if (sign<T>(signed_area2, eps) == 0) return std::nullopt;
if (!is_simple_polygon(polygon, eps)) return std::nullopt;
if (sign<T>(signed_area2, eps) < 0) {
std::reverse(polygon.begin(), polygon.end());
}
std::vector<std::size_t> remaining(polygon.size());
for (std::size_t index = 0; index < polygon.size(); ++index) {
remaining[index] = index;
}
std::vector<std::array<Point<T>, 3>> result;
result.reserve(polygon.size() - 2);
while (remaining.size() > 3) {
bool found_ear = false;
std::size_t size = remaining.size();
for (std::size_t position = 0; position < size; ++position) {
std::size_t previous_index =
remaining[(position + size - 1) % size];
std::size_t current_index = remaining[position];
std::size_t next_index =
remaining[(position + 1) % size];
const Point<T>& previous = polygon[previous_index];
const Point<T>& current = polygon[current_index];
const Point<T>& next = polygon[next_index];
if (orientation(previous, current, next, eps) <= 0) continue;
bool contains_vertex = false;
for (std::size_t other_index : remaining) {
if (
other_index == previous_index ||
other_index == current_index ||
other_index == next_index
) {
continue;
}
if (
polygon_detail::in_ccw_triangle(
polygon[other_index],
previous,
current,
next,
eps
)
) {
contains_vertex = true;
break;
}
}
if (contains_vertex) continue;
std::array<Point<T>, 3> triangle;
triangle[0] = previous;
triangle[1] = current;
triangle[2] = next;
result.push_back(std::move(triangle));
remaining.erase(
remaining.begin() +
static_cast<std::ptrdiff_t>(position)
);
found_ear = true;
break;
}
if (!found_ear) return std::nullopt;
}
std::array<Point<T>, 3> triangle;
triangle[0] = polygon[remaining[0]];
triangle[1] = polygon[remaining[1]];
triangle[2] = polygon[remaining[2]];
if (orientation(triangle[0], triangle[1], triangle[2], eps) <= 0) {
return std::nullopt;
}
result.push_back(std::move(triangle));
return result;
}
template <Coordinate T>
PointInPolygon point_in_polygon(
const std::vector<Point<T>>& polygon,
const Point<T>& point,
long double eps = 1e-12L
) {
bool inside = false;
std::size_t n = polygon.size();
for (std::size_t i = 0; i < n; i++) {
const Point<T>& a = polygon[i];
const Point<T>& b = polygon[(i + 1) % n];
if (on_segment(Segment<T>{a, b}, point, eps)) {
return PointInPolygon::Boundary;
}
if (a.y <= point.y) {
if (point.y < b.y && orientation(a, b, point, eps) > 0) {
inside = !inside;
}
} else if (b.y <= point.y && orientation(a, b, point, eps) < 0) {
inside = !inside;
}
}
return inside ? PointInPolygon::Inside : PointInPolygon::Outside;
}
template <Coordinate T, Coordinate P>
PointInPolygon point_in_polygon(
const Polygon<T>& polygon,
const Point<P>& point,
long double eps = 1e-12L
) {
assert(polygon.vertices.size() >= 3);
if constexpr (std::is_same_v<T, P>) {
return point_in_polygon(polygon.vertices, point, eps);
} else {
std::vector<Point<long double>> vertices;
vertices.reserve(polygon.vertices.size());
for (const Point<T>& vertex : polygon.vertices) {
vertices.emplace_back(vertex);
}
return point_in_polygon(vertices, Point<long double>(point), eps);
}
}
template <Coordinate T, Coordinate P>
bool contains(
const Polygon<T>& polygon,
const Point<P>& point,
long double eps = 1e-12L
) {
const PointInPolygon relation = point_in_polygon(polygon, point, eps);
return polygon.filled
? relation != PointInPolygon::Outside
: relation == PointInPolygon::Boundary;
}
namespace polygon_clip_detail {
struct Event {
long double parameter;
bool toggle;
};
inline bool close_parameter(
long double first,
long double second,
long double eps
) {
return std::fabs(first - second) <= eps * std::max({
1.0L,
std::fabs(first),
std::fabs(second)
});
}
inline long double parameter_on_line(
const Point<long double>& origin,
const Point<long double>& direction,
const Point<long double>& point
) {
return dot(point - origin, direction) / dot(direction, direction);
}
inline std::vector<Event> grouped_events(
std::vector<Event> events,
long double eps
) {
std::sort(
events.begin(),
events.end(),
[](const Event& first, const Event& second) {
return first.parameter < second.parameter;
}
);
std::vector<Event> result;
for (const Event& event : events) {
if (
result.empty() ||
!close_parameter(result.back().parameter, event.parameter, eps)
) {
result.push_back(event);
} else {
result.back().toggle = result.back().toggle != event.toggle;
}
}
return result;
}
inline std::vector<ParameterInterval> merge_intervals(
std::vector<ParameterInterval> intervals,
long double eps
) {
for (ParameterInterval& interval : intervals) {
if (interval.end < interval.begin) {
std::swap(interval.begin, interval.end);
}
}
std::sort(
intervals.begin(),
intervals.end(),
[](const ParameterInterval& first, const ParameterInterval& second) {
if (first.begin != second.begin) return first.begin < second.begin;
return first.end < second.end;
}
);
std::vector<ParameterInterval> result;
for (const ParameterInterval& interval : intervals) {
if (
result.empty() ||
interval.begin > result.back().end + eps * std::max({
1.0L,
std::fabs(interval.begin),
std::fabs(result.back().end)
})
) {
result.push_back(interval);
} else if (result.back().end < interval.end) {
result.back().end = interval.end;
}
}
return result;
}
template <Coordinate T>
std::vector<ParameterInterval> clip_line(
const Point<long double>& origin,
const Point<long double>& direction,
const Polygon<T>& polygon,
long double eps
) {
assert(polygon.vertices.size() >= 3);
assert(direction != Point<long double>());
std::vector<Event> events;
std::vector<ParameterInterval> intervals;
events.reserve(polygon.vertices.size() * 2);
intervals.reserve(polygon.vertices.size());
const Line<long double> line{origin, origin + direction};
for (std::size_t index = 0; index < polygon.vertices.size(); ++index) {
const Point<long double> first(polygon.vertices[index]);
const Point<long double> second(
polygon.vertices[(index + 1) % polygon.vertices.size()]
);
assert(first != second);
const Segment<long double> edge{first, second};
const LinearIntersection intersection =
linear_intersection(line, edge, eps);
if (intersection.kind == LinearIntersectionKind::Empty) continue;
if (intersection.kind == LinearIntersectionKind::Segment) {
intervals.push_back(ParameterInterval{
parameter_on_line(origin, direction, intersection.first),
parameter_on_line(origin, direction, intersection.second)
});
continue;
}
assert(intersection.kind == LinearIntersectionKind::Point);
const Point<long double> point = intersection.first;
const Point<long double> edge_direction = second - first;
const long double edge_parameter =
dot(point - first, edge_direction) /
dot(edge_direction, edge_direction);
bool toggle = false;
if (edge_parameter <= eps) {
toggle = orientation(line.a, line.b, second, eps) > 0;
} else if (edge_parameter >= 1.0L - eps) {
toggle = orientation(line.a, line.b, first, eps) > 0;
} else {
toggle = true;
}
events.push_back(Event{
parameter_on_line(origin, direction, point),
toggle
});
}
const std::vector<Event> grouped = grouped_events(std::move(events), eps);
if (polygon.filled) {
bool inside = false;
for (std::size_t index = 0; index < grouped.size(); ++index) {
if (index > 0 && inside) {
intervals.push_back(ParameterInterval{
grouped[index - 1].parameter,
grouped[index].parameter
});
}
intervals.push_back(ParameterInterval{
grouped[index].parameter,
grouped[index].parameter
});
inside = inside != grouped[index].toggle;
}
} else {
for (const Event& event : grouped) {
intervals.push_back(ParameterInterval{
event.parameter,
event.parameter
});
}
}
return merge_intervals(std::move(intervals), eps);
}
inline std::vector<ParameterInterval> restrict_domain(
const std::vector<ParameterInterval>& intervals,
long double lower,
long double upper,
long double eps
) {
std::vector<ParameterInterval> result;
result.reserve(intervals.size());
for (const ParameterInterval& interval : intervals) {
long double begin = std::max(interval.begin, lower);
long double end = std::min(interval.end, upper);
if (
end < begin &&
!close_parameter(begin, end, eps)
) {
continue;
}
if (end < begin) {
const long double middle = (begin + end) / 2.0L;
begin = middle;
end = middle;
}
result.push_back(ParameterInterval{begin, end});
}
return merge_intervals(std::move(result), eps);
}
inline std::vector<Event> grouped_circle_events(
std::vector<Event> events,
long double eps
) {
std::vector<Event> result = grouped_events(std::move(events), eps);
const long double full = 2.0L * std::numbers::pi_v<long double>;
if (
result.size() >= 2 &&
close_parameter(result.front().parameter + full,
result.back().parameter, eps)
) {
result.front().parameter = 0.0L;
result.front().toggle =
result.front().toggle != result.back().toggle;
result.pop_back();
}
return result;
}
template <Coordinate C, Coordinate T>
std::vector<AngularCoverage> clip_circle(
const Circle<C>& circle,
const Polygon<T>& polygon,
long double eps
) {
assert(circle.radius >= 0);
assert(polygon.vertices.size() >= 3);
if (circle.radius == 0) {
return contains(polygon, circle.center, eps)
? std::vector<AngularCoverage>{AngularCoverage{
AngularCoverageKind::Point,
0.0L,
0.0L
}}
: std::vector<AngularCoverage>();
}
std::vector<Event> events;
events.reserve(polygon.vertices.size() * 2);
for (std::size_t index = 0; index < polygon.vertices.size(); ++index) {
const Point<long double> first(polygon.vertices[index]);
const Point<long double> second(
polygon.vertices[(index + 1) % polygon.vertices.size()]
);
assert(first != second);
const Point<long double> direction = second - first;
const Segment<long double> edge{first, second};
const CircleLinearIntersection intersection =
circle_boundary_intersection(circle, edge, eps);
for (
int contact_index = 0;
contact_index < intersection.contact_count;
++contact_index
) {
const CircleLinearContact& contact =
intersection.contacts[contact_index];
const Point<long double> radial =
contact.point - Point<long double>(circle.center);
bool toggle = false;
if (contact.linear_parameter <= eps) {
toggle = predicate_detail::dot_sign<false>(
radial.x,
radial.y,
direction.x,
direction.y,
eps
) < 0;
} else if (contact.linear_parameter >= 1.0L - eps) {
toggle = predicate_detail::dot_sign<false>(
radial.x,
radial.y,
-direction.x,
-direction.y,
eps
) < 0;
} else {
toggle = predicate_detail::dot_sign<false>(
radial.x,
radial.y,
direction.x,
direction.y,
eps
) != 0;
}
events.push_back(Event{contact.circle_argument, toggle});
}
}
const std::vector<Event> grouped =
grouped_circle_events(std::move(events), eps);
const long double full = 2.0L * std::numbers::pi_v<long double>;
if (grouped.empty()) {
return contains(polygon, circle_point_at(circle, 0.0L), eps)
? std::vector<AngularCoverage>{AngularCoverage{
AngularCoverageKind::Full,
0.0L,
full
}}
: std::vector<AngularCoverage>();
}
std::vector<bool> inside_gap(grouped.size(), false);
if (polygon.filled) {
const long double wrap_middle = normalize_circle_argument(
(grouped.back().parameter + grouped.front().parameter + full) /
2.0L
);
bool inside = point_in_polygon(
polygon,
circle_point_at(circle, wrap_middle),
eps
) != PointInPolygon::Outside;
for (std::size_t index = 0; index < grouped.size(); ++index) {
inside = inside != grouped[index].toggle;
inside_gap[index] = inside;
}
}
if (
polygon.filled &&
std::all_of(
inside_gap.begin(),
inside_gap.end(),
[](bool inside) { return inside; }
)
) {
return {AngularCoverage{AngularCoverageKind::Full, 0.0L, full}};
}
std::vector<AngularCoverage> result;
for (std::size_t index = 0; index < grouped.size(); ++index) {
const std::size_t previous =
(index + grouped.size() - 1) % grouped.size();
if (!inside_gap[previous] && !inside_gap[index]) {
result.push_back(AngularCoverage{
AngularCoverageKind::Point,
grouped[index].parameter,
grouped[index].parameter
});
}
if (!inside_gap[previous] && inside_gap[index]) {
std::size_t finish = index;
while (inside_gap[finish]) {
finish = (finish + 1) % grouped.size();
}
long double end = grouped[finish].parameter;
if (end <= grouped[index].parameter) end += full;
result.push_back(AngularCoverage{
AngularCoverageKind::Arc,
grouped[index].parameter,
end
});
}
}
std::sort(
result.begin(),
result.end(),
[](const AngularCoverage& first, const AngularCoverage& second) {
return first.begin < second.begin;
}
);
return result;
}
} // namespace polygon_clip_detail
template <Coordinate L, Coordinate T>
std::vector<ParameterInterval> clip(
const Line<L>& line,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
assert(line.a != line.b);
assert(eps >= 0.0L);
const Point<long double> origin(line.a);
const Point<long double> direction =
Point<long double>(line.b) - origin;
return polygon_clip_detail::clip_line(
origin,
direction,
polygon,
eps
);
}
template <Coordinate R, Coordinate T>
std::vector<ParameterInterval> clip(
const Ray<R>& ray,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
assert(ray.origin != ray.through);
assert(eps >= 0.0L);
const Point<long double> origin(ray.origin);
const Point<long double> direction =
Point<long double>(ray.through) - origin;
return polygon_clip_detail::restrict_domain(
polygon_clip_detail::clip_line(origin, direction, polygon, eps),
0.0L,
std::numeric_limits<long double>::infinity(),
eps
);
}
template <Coordinate S, Coordinate T>
std::vector<ParameterInterval> clip(
const Segment<S>& segment,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
assert(eps >= 0.0L);
if (segment.a == segment.b) {
return contains(polygon, segment.a, eps)
? std::vector<ParameterInterval>{ParameterInterval{0.0L, 0.0L}}
: std::vector<ParameterInterval>();
}
const Point<long double> origin(segment.a);
const Point<long double> direction =
Point<long double>(segment.b) - origin;
return polygon_clip_detail::restrict_domain(
polygon_clip_detail::clip_line(origin, direction, polygon, eps),
0.0L,
1.0L,
eps
);
}
template <Coordinate C, Coordinate T>
std::vector<AngularCoverage> clip(
const Circle<C>& circle,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
assert(eps >= 0.0L);
return polygon_clip_detail::clip_circle(circle, polygon, eps);
}
template <Coordinate T>
bool intersects(
const Ray<T>& ray,
const std::vector<Point<T>>& polygon,
long double eps = 1e-12L
) {
assert(polygon.size() >= 3);
if (point_in_polygon(polygon, ray.origin, eps) != PointInPolygon::Outside) {
return true;
}
Polygon<T> region{polygon};
return !clip(ray, region, eps).empty();
}
template <Coordinate T>
bool intersects(
const std::vector<Point<T>>& polygon,
const Ray<T>& ray,
long double eps = 1e-12L
) {
return intersects(ray, polygon, eps);
}
template <Coordinate T>
long double distance(
const Ray<T>& ray,
const std::vector<Point<T>>& polygon
) {
assert(polygon.size() >= 3);
if (intersects(ray, polygon)) return 0;
long double result = std::numeric_limits<long double>::infinity();
std::size_t size = polygon.size();
for (std::size_t index = 0; index < size; ++index) {
result = std::min(
result,
distance(
ray,
Segment<T>{
polygon[index],
polygon[(index + 1) % size]
}
)
);
}
return result;
}
template <Coordinate T>
long double distance(
const std::vector<Point<T>>& polygon,
const Ray<T>& ray
) {
return distance(ray, polygon);
}
template <Coordinate T>
bool intersects(
const std::vector<Point<T>>& first,
const std::vector<Point<T>>& second,
long double eps = 1e-12L
) {
assert(first.size() >= 3);
assert(second.size() >= 3);
std::size_t first_size = first.size();
std::size_t second_size = second.size();
for (
std::size_t first_index = 0;
first_index < first_size;
++first_index
) {
Segment<T> first_edge{
first[first_index],
first[(first_index + 1) % first_size]
};
for (
std::size_t second_index = 0;
second_index < second_size;
++second_index
) {
Segment<T> second_edge{
second[second_index],
second[(second_index + 1) % second_size]
};
if (intersects(first_edge, second_edge, eps)) return true;
}
}
return
point_in_polygon(first, second.front(), eps) !=
PointInPolygon::Outside ||
point_in_polygon(second, first.front(), eps) !=
PointInPolygon::Outside;
}
template <Coordinate T>
long double distance(
const std::vector<Point<T>>& first,
const std::vector<Point<T>>& second
) {
assert(first.size() >= 3);
assert(second.size() >= 3);
if (intersects(first, second)) return 0;
long double result = std::numeric_limits<long double>::infinity();
std::size_t first_size = first.size();
std::size_t second_size = second.size();
for (
std::size_t first_index = 0;
first_index < first_size;
++first_index
) {
Segment<T> first_edge{
first[first_index],
first[(first_index + 1) % first_size]
};
for (
std::size_t second_index = 0;
second_index < second_size;
++second_index
) {
Segment<T> second_edge{
second[second_index],
second[(second_index + 1) % second_size]
};
result = std::min(result, distance(first_edge, second_edge));
}
}
return result;
}
template <Coordinate T>
wide_type<T> polygon_area2(const Polygon<T>& polygon) {
return polygon_area2(polygon.vertices);
}
template <Coordinate T>
long double polygon_area(const Polygon<T>& polygon) {
return polygon_area(polygon.vertices);
}
template <Coordinate T>
std::optional<Point<long double>> polygon_centroid(
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
return polygon_centroid(polygon.vertices, eps);
}
template <Coordinate T>
std::optional<Point<long double>> centroid(
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
return polygon_centroid(polygon.vertices, eps);
}
namespace polygon_detail {
template <Coordinate T>
Segment<long double> edge(const Polygon<T>& polygon, std::size_t index) {
return Segment<long double>{
Point<long double>(polygon.vertices[index]),
Point<long double>(
polygon.vertices[(index + 1) % polygon.vertices.size()]
)
};
}
template <Coordinate T>
ClosestPoints closest_boundary_point(
const Polygon<T>& polygon,
const Point<long double>& point
) {
assert(polygon.vertices.size() >= 3);
ClosestPoints result = closest_points(edge(polygon, 0), point);
for (std::size_t index = 1; index < polygon.vertices.size(); ++index) {
closest_points_detail::consider(
result,
closest_points(edge(polygon, index), point)
);
}
return result;
}
template <Coordinate T, class Object>
ClosestPoints closest_boundary_object(
const Polygon<T>& polygon,
const Object& object
) {
assert(polygon.vertices.size() >= 3);
ClosestPoints result = closest_points(edge(polygon, 0), object);
for (std::size_t index = 1; index < polygon.vertices.size(); ++index) {
closest_points_detail::consider(
result,
closest_points(edge(polygon, index), object)
);
}
return result;
}
template <Coordinate A, Coordinate B>
ClosestPoints closest_boundaries(
const Polygon<A>& first,
const Polygon<B>& second
) {
assert(first.vertices.size() >= 3);
assert(second.vertices.size() >= 3);
ClosestPoints result = closest_points(edge(first, 0), edge(second, 0));
for (
std::size_t first_index = 0;
first_index < first.vertices.size();
++first_index
) {
for (
std::size_t second_index = 0;
second_index < second.vertices.size();
++second_index
) {
closest_points_detail::consider(
result,
closest_points(
edge(first, first_index),
edge(second, second_index)
)
);
}
}
return result;
}
} // namespace polygon_detail
template <Coordinate T, Coordinate P>
ClosestPoints closest_points(
const Polygon<T>& polygon,
const Point<P>& point,
long double eps = 1e-12L
) {
assert(polygon.vertices.size() >= 3);
const Point<long double> converted(point);
if (polygon.filled && contains(polygon, point, eps)) {
return ClosestPoints{converted, converted};
}
return polygon_detail::closest_boundary_point(polygon, converted);
}
template <Coordinate P, Coordinate T>
ClosestPoints closest_points(
const Point<P>& point,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
return closest_points_detail::reversed(
closest_points(polygon, point, eps)
);
}
template <Coordinate T, Coordinate S>
ClosestPoints closest_points(
const Polygon<T>& polygon,
const Segment<S>& segment,
long double eps = 1e-12L
) {
assert(polygon.vertices.size() >= 3);
const Segment<long double> converted{
Point<long double>(segment.a),
Point<long double>(segment.b)
};
if (polygon.filled) {
if (contains(polygon, segment.a, eps)) {
const Point<long double> point(segment.a);
return ClosestPoints{point, point};
}
if (contains(polygon, segment.b, eps)) {
const Point<long double> point(segment.b);
return ClosestPoints{point, point};
}
}
return polygon_detail::closest_boundary_object(polygon, converted);
}
template <Coordinate S, Coordinate T>
ClosestPoints closest_points(
const Segment<S>& segment,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
return closest_points_detail::reversed(
closest_points(polygon, segment, eps)
);
}
template <Coordinate T, Coordinate R>
ClosestPoints closest_points(
const Polygon<T>& polygon,
const Ray<R>& ray,
long double eps = 1e-12L
) {
assert(polygon.vertices.size() >= 3);
const Ray<long double> converted{
Point<long double>(ray.origin),
Point<long double>(ray.through)
};
if (polygon.filled && contains(polygon, ray.origin, eps)) {
const Point<long double> point(ray.origin);
return ClosestPoints{point, point};
}
return polygon_detail::closest_boundary_object(polygon, converted);
}
template <Coordinate R, Coordinate T>
ClosestPoints closest_points(
const Ray<R>& ray,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
return closest_points_detail::reversed(
closest_points(polygon, ray, eps)
);
}
template <Coordinate A, Coordinate B>
ClosestPoints closest_points(
const Polygon<A>& first,
const Polygon<B>& second,
long double eps = 1e-12L
) {
assert(first.vertices.size() >= 3);
assert(second.vertices.size() >= 3);
ClosestPoints result = polygon_detail::closest_boundaries(first, second);
if (geometry::distance(result.first, result.second) <= eps) return result;
if (first.filled) {
for (const Point<B>& vertex : second.vertices) {
if (contains(first, vertex, eps)) {
const Point<long double> point(vertex);
return ClosestPoints{point, point};
}
}
}
if (second.filled) {
for (const Point<A>& vertex : first.vertices) {
if (contains(second, vertex, eps)) {
const Point<long double> point(vertex);
return ClosestPoints{point, point};
}
}
}
return result;
}
template <Coordinate C, Coordinate T>
ClosestPoints closest_points(
const Circle<C>& circle,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
assert(polygon.vertices.size() >= 3);
ClosestPoints result = closest_points(
circle,
polygon_detail::edge(polygon, 0),
eps
);
for (std::size_t index = 1; index < polygon.vertices.size(); ++index) {
closest_points_detail::consider(
result,
closest_points(circle, polygon_detail::edge(polygon, index), eps)
);
}
if (geometry::distance(result.first, result.second) <= eps) return result;
if (polygon.filled) {
Point<long double> member(circle.center);
if (!circle.filled) {
member = circle_detail::point_toward(circle, member);
}
if (contains(polygon, member, eps)) {
return ClosestPoints{member, member};
}
}
return result;
}
template <Coordinate T, Coordinate C>
ClosestPoints closest_points(
const Polygon<T>& polygon,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return closest_points_detail::reversed(
closest_points(circle, polygon, eps)
);
}
template <Coordinate T, Coordinate P>
bool intersects(
const Polygon<T>& polygon,
const Point<P>& point,
long double eps = 1e-12L
) {
return contains(polygon, point, eps);
}
template <Coordinate P, Coordinate T>
bool intersects(
const Point<P>& point,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
return intersects(polygon, point, eps);
}
template <Coordinate T, Coordinate S>
bool intersects(
const Polygon<T>& polygon,
const Segment<S>& segment,
long double eps = 1e-12L
) {
const ClosestPoints result = closest_points(polygon, segment, eps);
return geometry::distance(result.first, result.second) <= eps;
}
template <Coordinate S, Coordinate T>
bool intersects(
const Segment<S>& segment,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
return intersects(polygon, segment, eps);
}
template <Coordinate T, Coordinate R>
bool intersects(
const Polygon<T>& polygon,
const Ray<R>& ray,
long double eps = 1e-12L
) {
const ClosestPoints result = closest_points(polygon, ray, eps);
return geometry::distance(result.first, result.second) <= eps;
}
template <Coordinate R, Coordinate T>
bool intersects(
const Ray<R>& ray,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
return intersects(polygon, ray, eps);
}
template <Coordinate A, Coordinate B>
bool intersects(
const Polygon<A>& first,
const Polygon<B>& second,
long double eps = 1e-12L
) {
const ClosestPoints result = closest_points(first, second, eps);
return geometry::distance(result.first, result.second) <= eps;
}
template <Coordinate C, Coordinate T>
bool intersects(
const Circle<C>& circle,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
const ClosestPoints result = closest_points(circle, polygon, eps);
return geometry::distance(result.first, result.second) <= eps;
}
template <Coordinate T, Coordinate C>
bool intersects(
const Polygon<T>& polygon,
const Circle<C>& circle,
long double eps = 1e-12L
) {
return intersects(circle, polygon, eps);
}
template <Coordinate A, Coordinate B>
long double distance(
const Polygon<A>& first,
const Polygon<B>& second
) {
const ClosestPoints result = closest_points(first, second);
return geometry::distance(result.first, result.second);
}
template <Coordinate C, Coordinate T>
long double distance(
const Circle<C>& circle,
const Polygon<T>& polygon
) {
const ClosestPoints result = closest_points(circle, polygon);
return geometry::distance(result.first, result.second);
}
template <Coordinate T, Coordinate C>
long double distance(
const Polygon<T>& polygon,
const Circle<C>& circle
) {
return distance(circle, polygon);
}
template <Coordinate C, Coordinate T>
long double circle_polygon_intersection_area(
const Circle<C>& circle,
const Polygon<T>& polygon,
long double eps = 1e-12L
) {
return circle_polygon_intersection_area(circle, polygon.vertices, eps);
}
template <Coordinate T, Coordinate P>
long double distance(
const Polygon<T>& polygon,
const Point<P>& point
) {
const ClosestPoints result = closest_points(polygon, point);
return geometry::distance(result.first, result.second);
}
template <Coordinate P, Coordinate T>
long double distance(
const Point<P>& point,
const Polygon<T>& polygon
) {
return distance(polygon, point);
}
template <Coordinate T, Coordinate S>
long double distance(
const Polygon<T>& polygon,
const Segment<S>& segment
) {
const ClosestPoints result = closest_points(polygon, segment);
return geometry::distance(result.first, result.second);
}
template <Coordinate S, Coordinate T>
long double distance(
const Segment<S>& segment,
const Polygon<T>& polygon
) {
return distance(polygon, segment);
}
template <Coordinate T, Coordinate R>
long double distance(
const Polygon<T>& polygon,
const Ray<R>& ray
) {
const ClosestPoints result = closest_points(polygon, ray);
return geometry::distance(result.first, result.second);
}
template <Coordinate R, Coordinate T>
long double distance(
const Ray<R>& ray,
const Polygon<T>& polygon
) {
return distance(polygon, ray);
}
} // namespace geometry
} // namespace m1une
#line 21 "geometry/convex_decomposition.hpp"
namespace m1une {
namespace geometry {
namespace convex_decomposition_detail {
using Index = std::size_t;
using IndexPolygon = std::vector<Index>;
enum class ExactPredicateWidth { Int128, Int256, Int512 };
template <std::integral T>
std::size_t coordinate_magnitude_bits(
const std::vector<Point<T>>& polygon
) {
using Unsigned = std::make_unsigned_t<T>;
std::size_t result = 0;
auto update = [&](T coordinate) {
Unsigned magnitude = static_cast<Unsigned>(coordinate);
if constexpr (std::signed_integral<T>) {
if (coordinate < 0) magnitude = Unsigned(0) - magnitude;
}
std::size_t bits = 0;
while (magnitude != 0) {
++bits;
magnitude >>= 1;
}
result = std::max(result, bits);
};
for (const Point<T>& point : polygon) {
update(point.x);
update(point.y);
}
return result;
}
template <std::integral T>
ExactPredicateWidth select_exact_predicate_width(
const std::vector<Point<T>>& polygon
) {
static_assert(sizeof(T) <= sizeof(std::uint64_t));
const std::size_t coordinate_bits =
coordinate_magnitude_bits(polygon);
// Cross-multiplied ray parameters and projective visibility
// determinants have magnitude below 2^(4 * coordinate_bits + 6).
const std::size_t required_bits = 4 * coordinate_bits + 7;
if (required_bits <= 128) return ExactPredicateWidth::Int128;
if (required_bits <= 256) return ExactPredicateWidth::Int256;
assert(required_bits <= 512);
return ExactPredicateWidth::Int512;
}
template <typename Number, Coordinate T>
Number predicate_number(T value) {
return static_cast<Number>(value);
}
template <Coordinate T, typename Number>
int predicate_sign(const Number& value, long double eps) {
if constexpr (ExactCoordinate<T>) {
return (value > 0) - (value < 0);
} else {
return (value > eps) - (value < -eps);
}
}
template <Coordinate T>
std::optional<std::vector<Point<T>>> prepare_polygon(
std::vector<Point<T>> polygon,
long double eps
) {
polygon =
polygon_detail::clean_polygon_vertices(std::move(polygon), eps);
if (polygon.size() < 3) return std::nullopt;
const int area_sign = sign<T>(polygon_area2(polygon), eps);
if (area_sign == 0 || !is_simple_polygon(polygon, eps)) {
return std::nullopt;
}
if (area_sign < 0) std::reverse(polygon.begin(), polygon.end());
return polygon;
}
template <Coordinate T>
bool is_weakly_convex(
const std::vector<Point<T>>& polygon,
long double eps
) {
if (polygon.size() < 3) return false;
for (Index index = 0; index < polygon.size(); ++index) {
if (
orientation(
polygon[index],
polygon[(index + 1) % polygon.size()],
polygon[(index + 2) % polygon.size()],
eps
) < 0
) {
return false;
}
}
return true;
}
template <Coordinate T>
std::optional<std::vector<IndexPolygon>> triangulate_indices(
const std::vector<Point<T>>& polygon,
long double eps
) {
const Index size = polygon.size();
std::vector<Index> previous(size), next(size);
std::vector<bool> active(size, true);
for (Index index = 0; index < size; ++index) {
previous[index] = (index + size - 1) % size;
next[index] = (index + 1) % size;
}
auto is_ear = [&](Index index) {
if (!active[index]) return false;
const Index first = previous[index];
const Index third = next[index];
if (
orientation(
polygon[first], polygon[index], polygon[third], eps
) <= 0
) {
return false;
}
for (Index other = 0; other < size; ++other) {
if (
!active[other] || other == first || other == index ||
other == third
) {
continue;
}
if (
polygon_detail::in_ccw_triangle(
polygon[other],
polygon[first],
polygon[index],
polygon[third],
eps
)
) {
return false;
}
}
return true;
};
std::deque<Index> ears;
for (Index index = 0; index < size; ++index) {
if (is_ear(index)) ears.push_back(index);
}
std::vector<IndexPolygon> triangles;
triangles.reserve(size - 2);
Index remaining = size;
while (remaining > 3) {
while (!ears.empty() && !is_ear(ears.front())) {
ears.pop_front();
}
if (ears.empty()) return std::nullopt;
const Index ear = ears.front();
ears.pop_front();
const Index first = previous[ear];
const Index third = next[ear];
triangles.push_back(IndexPolygon{first, ear, third});
active[ear] = false;
next[first] = third;
previous[third] = first;
--remaining;
if (is_ear(first)) ears.push_back(first);
if (is_ear(third)) ears.push_back(third);
}
Index first = 0;
while (first < size && !active[first]) ++first;
if (first == size) return std::nullopt;
const Index second = next[first];
const Index third = next[second];
if (
third == first || next[third] != first ||
orientation(
polygon[first], polygon[second], polygon[third], eps
) <= 0
) {
return std::nullopt;
}
triangles.push_back(IndexPolygon{first, second, third});
return triangles;
}
inline Index find_root(std::vector<Index>& parent, Index index) {
Index root = index;
while (parent[root] != root) root = parent[root];
while (parent[index] != index) {
const Index next = parent[index];
parent[index] = root;
index = next;
}
return root;
}
inline std::optional<IndexPolygon> merge_across_edge(
const IndexPolygon& first,
const IndexPolygon& second,
Index edge_first,
Index edge_second
) {
for (Index first_position = 0;
first_position < first.size();
++first_position) {
const Index first_next = (first_position + 1) % first.size();
const Index from = first[first_position];
const Index to = first[first_next];
if (
!(
(from == edge_first && to == edge_second) ||
(from == edge_second && to == edge_first)
)
) {
continue;
}
for (Index second_position = 0;
second_position < second.size();
++second_position) {
const Index second_next =
(second_position + 1) % second.size();
if (
second[second_position] != to ||
second[second_next] != from
) {
continue;
}
IndexPolygon merged;
merged.reserve(first.size() + second.size() - 2);
for (Index position = first_next;
position != first_position;
position = (position + 1) % first.size()) {
merged.push_back(first[position]);
}
for (Index position = second_next;
position != second_position;
position = (position + 1) % second.size()) {
merged.push_back(second[position]);
}
return merged;
}
}
return std::nullopt;
}
template <Coordinate T>
bool is_weakly_convex(
const IndexPolygon& polygon,
const std::vector<Point<T>>& points,
long double eps
) {
for (Index index = 0; index < polygon.size(); ++index) {
if (
orientation(
points[polygon[index]],
points[polygon[(index + 1) % polygon.size()]],
points[polygon[(index + 2) % polygon.size()]],
eps
) < 0
) {
return false;
}
}
return true;
}
template <Coordinate T>
std::vector<Point<T>> materialize(
const IndexPolygon& indices,
const std::vector<Point<T>>& points,
long double eps
) {
std::vector<Point<T>> polygon;
polygon.reserve(indices.size());
for (const Index index : indices) polygon.push_back(points[index]);
return polygon_detail::clean_polygon_vertices(std::move(polygon), eps);
}
struct Diagonal {
Index first;
Index second;
};
template <Coordinate T>
std::optional<std::vector<Point<T>>> prepare_minimum_polygon(
std::vector<Point<T>> polygon,
long double eps
) {
if (polygon.size() >= 2 && polygon.front() == polygon.back()) {
polygon.pop_back();
}
std::vector<Point<T>> distinct;
distinct.reserve(polygon.size());
for (const Point<T>& point : polygon) {
if (distinct.empty() || distinct.back() != point) {
distinct.push_back(point);
}
}
if (distinct.size() >= 2 && distinct.front() == distinct.back()) {
distinct.pop_back();
}
if (distinct.size() < 3) return std::nullopt;
const int original_size = static_cast<int>(distinct.size());
std::vector<int> previous(original_size), next(original_size);
std::vector<bool> removed(original_size, false);
std::deque<int> candidates;
for (int index = 0; index < original_size; ++index) {
previous[index] = (index + original_size - 1) % original_size;
next[index] = (index + 1) % original_size;
candidates.push_back(index);
}
int remaining = original_size;
while (!candidates.empty() && remaining >= 3) {
const int index = candidates.front();
candidates.pop_front();
if (removed[index]) continue;
const int before = previous[index];
const int after = next[index];
if (
orientation(
distinct[before], distinct[index], distinct[after], eps
) != 0 ||
sign<T>(
dot(
distinct[index] - distinct[before],
distinct[after] - distinct[index]
),
eps
) < 0
) {
continue;
}
removed[index] = true;
next[before] = after;
previous[after] = before;
--remaining;
candidates.push_back(before);
candidates.push_back(after);
}
if (remaining < 3) return std::nullopt;
std::vector<Point<T>> cleaned;
cleaned.reserve(static_cast<Index>(remaining));
int first = 0;
while (removed[first]) ++first;
int index = first;
do {
cleaned.push_back(distinct[index]);
index = next[index];
} while (index != first);
const int area_sign = sign<T>(polygon_area2(cleaned), eps);
if (area_sign == 0) return std::nullopt;
if (area_sign < 0) std::reverse(cleaned.begin(), cleaned.end());
return cleaned;
}
template <Coordinate T, typename Number>
class BiasedPolygonReduction {
private:
struct Vector {
Number x;
Number y;
};
struct Fraction {
Number numerator;
Number denominator;
};
struct Chain {
int first_edge;
int edge_count;
};
public:
struct Result {
std::vector<Point<T>> polygon;
std::vector<Index> original_index;
};
BiasedPolygonReduction(
const std::vector<Point<T>>& polygon,
long double eps
)
: polygon_(polygon),
eps_(eps),
size_(static_cast<int>(polygon.size())),
marked_edge_(polygon.size(), false) {
for (int index = 0; index < size_; ++index) {
if (
orientation(
polygon_[(index + size_ - 1) % size_],
polygon_[index],
polygon_[(index + 1) % size_],
eps_
) < 0
) {
reflex_vertices_.push_back(index);
}
}
build_chains();
}
Result run() {
for (const int reflex : reflex_vertices_) {
marked_edge_[(reflex + size_ - 1) % size_] = true;
marked_edge_[reflex] = true;
}
std::vector<std::vector<int>> extension_endpoints(size_);
for (const int reflex : reflex_vertices_) {
const Vector direction = vector_between(
reflex, (reflex + 1) % size_
);
for (const int edge : ray_shoot(reflex, direction)) {
marked_edge_[edge] = true;
extension_endpoints[reflex].push_back(edge);
extension_endpoints[reflex].push_back((edge + 1) % size_);
}
}
for (Index first = 0; first < reflex_vertices_.size(); ++first) {
for (Index second = first + 1;
second < reflex_vertices_.size();
++second) {
const int first_vertex = reflex_vertices_[first];
const int second_vertex = reflex_vertices_[second];
mark_ray(
first_vertex,
vector_between(first_vertex, second_vertex)
);
mark_ray(
second_vertex,
vector_between(second_vertex, first_vertex)
);
}
}
for (const int reflex : reflex_vertices_) {
auto& endpoints = extension_endpoints[reflex];
std::sort(endpoints.begin(), endpoints.end());
endpoints.erase(
std::unique(endpoints.begin(), endpoints.end()),
endpoints.end()
);
for (const int endpoint : endpoints) {
if (endpoint == reflex) continue;
mark_ray(
reflex, vector_between(reflex, endpoint)
);
}
}
Result result;
for (int index = 0; index < size_; ++index) {
if (
marked_edge_[index] ||
marked_edge_[(index + size_ - 1) % size_]
) {
result.polygon.push_back(polygon_[index]);
result.original_index.push_back(static_cast<Index>(index));
}
}
return result;
}
private:
const std::vector<Point<T>>& polygon_;
long double eps_;
int size_;
std::vector<int> reflex_vertices_;
std::vector<Chain> chains_;
std::vector<bool> marked_edge_;
Vector vector_between(int first, int second) const {
return {
predicate_number<Number>(polygon_[first].x) -
predicate_number<Number>(polygon_[second].x),
predicate_number<Number>(polygon_[first].y) -
predicate_number<Number>(polygon_[second].y)
};
}
Number vector_cross(const Vector& first, const Vector& second) const {
return first.x * second.y - first.y * second.x;
}
Number vector_dot(const Vector& first, const Vector& second) const {
return first.x * second.x + first.y * second.y;
}
int vector_sign(const Number& value) const {
return predicate_sign<T>(value, eps_);
}
int quadrant(const Vector& direction) const {
const int x_sign = vector_sign(direction.x);
const int y_sign = vector_sign(direction.y);
if (y_sign >= 0) return x_sign >= 0 ? 0 : 1;
return x_sign < 0 ? 2 : 3;
}
void build_chains() {
// Within one quadrant, monotonically turning edge directions span at
// most pi/2. This gives the bitonic sidedness needed by ray shooting
// without evaluating an angle.
int first_edge = 0;
Vector previous_direction = vector_between(1, 0);
int current_quadrant = quadrant(previous_direction);
for (int edge = 1; edge < size_; ++edge) {
const Vector direction = vector_between(
(edge + 1) % size_, edge
);
const int direction_quadrant = quadrant(direction);
if (
vector_sign(vector_cross(previous_direction, direction)) <
0 ||
direction_quadrant != current_quadrant
) {
chains_.push_back(Chain{
first_edge, edge - first_edge
});
first_edge = edge;
current_quadrant = direction_quadrant;
}
previous_direction = direction;
}
chains_.push_back(Chain{first_edge, size_ - first_edge});
}
int side(
int origin,
const Vector& direction,
int vertex
) const {
return vector_sign(
vector_cross(
direction,
vector_between(vertex % size_, origin)
)
);
}
int edge_side(const Vector& direction, int edge) const {
return vector_sign(
vector_cross(
direction,
vector_between((edge + 1) % size_, edge)
)
);
}
void crossing_on_monotone_part(
const Chain& chain,
int origin,
const Vector& direction,
int first_position,
int last_position,
std::vector<int>& candidates
) const {
if (first_position >= last_position) return;
const int first_sign = side(
origin,
direction,
chain.first_edge + first_position
);
const int last_sign = side(
origin,
direction,
chain.first_edge + last_position
);
if (first_sign == 0) {
candidates.push_back(
chain.first_edge + first_position
);
}
if (last_sign == 0) {
candidates.push_back(
chain.first_edge + last_position - 1
);
}
if (first_sign == 0 || last_sign == 0 || first_sign == last_sign) {
return;
}
int low = first_position;
int high = last_position;
while (high - low > 1) {
const int middle = (low + high) / 2;
const int middle_sign = side(
origin,
direction,
chain.first_edge + middle
);
if (middle_sign == 0 || middle_sign != first_sign) {
high = middle;
} else {
low = middle;
}
}
candidates.push_back(chain.first_edge + high - 1);
}
void chain_candidates(
const Chain& chain,
int origin,
const Vector& direction,
std::vector<int>& candidates
) const {
int split = 0;
const int first_derivative =
edge_side(direction, chain.first_edge);
const int last_derivative = edge_side(
direction,
chain.first_edge + chain.edge_count - 1
);
if (
first_derivative != 0 && last_derivative != 0 &&
first_derivative != last_derivative
) {
int low = 0;
int high = chain.edge_count - 1;
while (high - low > 1) {
const int middle = (low + high) / 2;
const int middle_sign = edge_side(
direction, chain.first_edge + middle
);
if (
middle_sign == 0 ||
middle_sign != first_derivative
) {
high = middle;
} else {
low = middle;
}
}
split = high;
}
if (split == 0) {
crossing_on_monotone_part(
chain,
origin,
direction,
0,
chain.edge_count,
candidates
);
} else {
crossing_on_monotone_part(
chain, origin, direction, 0, split, candidates
);
crossing_on_monotone_part(
chain,
origin,
direction,
split,
chain.edge_count,
candidates
);
}
}
std::vector<int> ray_shoot(
int origin,
const Vector& direction
) const {
std::vector<int> candidates;
for (const Chain& chain : chains_) {
chain_candidates(chain, origin, direction, candidates);
}
std::sort(candidates.begin(), candidates.end());
candidates.erase(
std::unique(candidates.begin(), candidates.end()),
candidates.end()
);
if constexpr (ExactCoordinate<T>) {
return exact_ray_shoot(origin, direction, candidates);
} else {
return floating_ray_shoot(origin, direction, candidates);
}
}
std::vector<int> exact_ray_shoot(
int origin,
const Vector& direction,
const std::vector<int>& candidates
) const {
// Intersections are rational parameters along the ray. Keep their
// numerators and denominators and compare them by cross products.
std::optional<Fraction> best;
std::vector<int> result;
const Number direction_norm2 = vector_dot(direction, direction);
for (int edge : candidates) {
edge %= size_;
const Vector first = vector_between(edge, origin);
const Vector second = vector_between(
(edge + 1) % size_, origin
);
const Vector edge_direction = vector_between(
(edge + 1) % size_, edge
);
Number denominator = vector_cross(direction, edge_direction);
Fraction parameter{0, 1};
bool valid = false;
if (denominator == 0) {
if (vector_cross(direction, first) != 0) continue;
const Number first_parameter =
vector_dot(first, direction);
const Number second_parameter =
vector_dot(second, direction);
if (first_parameter > 0) {
parameter = Fraction{
first_parameter, direction_norm2
};
valid = true;
}
if (
second_parameter > 0 &&
(!valid || second_parameter < parameter.numerator)
) {
parameter = Fraction{
second_parameter, direction_norm2
};
valid = true;
}
} else {
Number numerator = vector_cross(first, edge_direction);
Number segment_numerator = vector_cross(first, direction);
if (denominator < 0) {
denominator = -denominator;
numerator = -numerator;
segment_numerator = -segment_numerator;
}
if (
numerator <= 0 || segment_numerator < 0 ||
segment_numerator > denominator
) {
continue;
}
parameter = Fraction{numerator, denominator};
valid = true;
}
if (!valid) continue;
if (!best.has_value()) {
best = parameter;
result.assign(1, edge);
continue;
}
const Number left =
parameter.numerator * best->denominator;
const Number right =
best->numerator * parameter.denominator;
if (left < right) {
best = parameter;
result.assign(1, edge);
} else if (left == right) {
result.push_back(edge);
}
}
return result;
}
std::vector<int> floating_ray_shoot(
int origin,
const Vector& direction,
const std::vector<int>& candidates
) const {
long double best = std::numeric_limits<long double>::infinity();
std::vector<int> result;
const long double direction_norm2 = vector_dot(
direction, direction
);
for (int edge : candidates) {
edge %= size_;
const Vector first = vector_between(edge, origin);
const Vector second = vector_between(
(edge + 1) % size_, origin
);
const Vector edge_direction = vector_between(
(edge + 1) % size_, edge
);
const long double denominator = vector_cross(
direction, edge_direction
);
long double parameter = -1;
if (std::fabs(denominator) <= eps_) {
if (std::fabs(vector_cross(direction, first)) > eps_) {
continue;
}
const long double first_parameter =
vector_dot(first, direction) / direction_norm2;
const long double second_parameter =
vector_dot(second, direction) / direction_norm2;
if (first_parameter > eps_) parameter = first_parameter;
if (
second_parameter > eps_ &&
(parameter < 0 || second_parameter < parameter)
) {
parameter = second_parameter;
}
} else {
parameter =
vector_cross(first, edge_direction) / denominator;
const long double segment_parameter =
vector_cross(first, direction) / denominator;
if (
parameter <= eps_ || segment_parameter < -eps_ ||
segment_parameter > 1 + eps_
) {
continue;
}
}
if (parameter < 0) continue;
if (parameter + eps_ < best) {
best = parameter;
result.assign(1, edge);
} else if (std::fabs(parameter - best) <= eps_) {
result.push_back(edge);
}
}
return result;
}
void mark_ray(int origin, const Vector& direction) {
for (const int edge : ray_shoot(origin, direction)) {
marked_edge_[edge] = true;
}
}
};
template <Coordinate T, typename Number>
class KeilSnoeyinkDecomposition {
private:
static constexpr int infinity = 100000000;
static constexpr int bad = 1000000000;
struct SolutionNode;
using NodePointer = std::shared_ptr<const SolutionNode>;
struct SolutionNode {
NodePointer first;
NodePointer second;
Diagonal diagonal{0, 0};
bool is_diagonal = false;
};
struct NarrowPair {
int first;
int second;
NodePointer solution;
};
// Two independently popped views of the same narrowest-pair stack.
// Pairs are appended from back to front in the terminology of the
// Keil--Snoeyink paper.
struct PairDeque {
std::vector<NarrowPair> pairs;
int front = -1;
int back = 0;
bool empty_front() const { return front < 0; }
bool more_front() const { return front > 0; }
bool empty_back() const {
return back >= static_cast<int>(pairs.size());
}
bool more_back() const {
return back + 1 < static_cast<int>(pairs.size());
}
const NarrowPair& front_pair() const { return pairs[front]; }
const NarrowPair& under_front() const {
return pairs[front - 1];
}
const NarrowPair& back_pair() const { return pairs[back]; }
const NarrowPair& under_back() const {
return pairs[back + 1];
}
void pop_front() { --front; }
void pop_back() { ++back; }
void restore() {
front = static_cast<int>(pairs.size()) - 1;
back = 0;
}
void clear() {
pairs.clear();
restore();
}
void push(int first, int second, NodePointer solution = nullptr) {
pairs.push_back(NarrowPair{first, second, std::move(solution)});
restore();
}
void push_narrow(
int first,
int second,
NodePointer solution
) {
if (!empty_front() && first <= front_pair().first) return;
while (!empty_front() && front_pair().second >= second) {
pairs.pop_back();
--front;
}
push(first, second, std::move(solution));
}
};
struct State {
int weight = bad;
PairDeque pairs;
};
using ProjectiveNumber = Number;
struct Homogeneous {
ProjectiveNumber w;
ProjectiveNumber x;
ProjectiveNumber y;
Homogeneous negated() const { return {-w, -x, -y}; }
ProjectiveNumber side(const Point<T>& point) const {
return w + x * predicate_number<ProjectiveNumber>(point.x) +
y * predicate_number<ProjectiveNumber>(point.y);
}
static Homogeneous line(
const Point<T>& first,
const Point<T>& second
) {
const ProjectiveNumber ax =
predicate_number<ProjectiveNumber>(first.x);
const ProjectiveNumber ay =
predicate_number<ProjectiveNumber>(first.y);
const ProjectiveNumber bx =
predicate_number<ProjectiveNumber>(second.x);
const ProjectiveNumber by =
predicate_number<ProjectiveNumber>(second.y);
return {ax * by - ay * bx, ay - by, bx - ax};
}
static ProjectiveNumber determinant(
const Homogeneous& first,
const Homogeneous& second,
const Homogeneous& third
) {
return
first.w *
(second.x * third.y - second.y * third.x) -
second.w *
(first.x * third.y - first.y * third.x) +
third.w *
(first.x * second.y - first.y * second.x);
}
};
// ElGindy--Avis visibility-polygon scan, specialized to a polygon
// vertex. Only original polygon vertices are returned; artificial lid
// intersections are discarded.
class VisibilityPolygon {
private:
enum VertexType { right_lid, left_lid, right_wall, left_wall };
public:
VisibilityPolygon(
const std::vector<Point<T>>& polygon,
long double eps
)
: polygon_(polygon), eps_(eps), size_(polygon.size()) {}
std::vector<Index> build(Index origin_index) {
origin_ = &polygon_[origin_index];
vertices_.clear();
types_.clear();
left_lid_index_ = not_saved;
right_lid_index_ = not_saved;
int vertex = static_cast<int>(origin_index);
push(vertex++, right_wall);
do {
push(vertex++, right_wall);
if (vertex >= static_cast<int>(size_ + origin_index)) {
break;
}
Homogeneous edge = line(vertex - 1, vertex);
if (left(edge, *origin_)) continue;
if (!left(edge, point(vertex - 2))) {
vertex = exit_right_bay(
vertex,
top(),
Homogeneous{1, 0, 0}
);
push(vertex++, right_lid);
continue;
}
save_lid();
for (;;) {
if (
orientation(
*origin_, top(), point(vertex), eps_
) > 0
) {
if (
orientation(
point(vertex),
point(vertex + 1),
*origin_,
eps_
) < 0
) {
++vertex;
} else if (left(edge, point(vertex + 1))) {
vertex = exit_left_bay(
vertex,
point(vertex),
line(left_lid_index_, left_lid_index_ - 1)
) + 1;
} else {
restore_lid();
push(vertex++, left_wall);
break;
}
edge = line(vertex - 1, vertex);
} else if (!left(edge, top())) {
if (right_lid_index_ == not_saved) break;
vertex = exit_right_bay(
vertex, top(), edge.negated()
);
push(vertex++, right_lid);
break;
} else {
save_lid();
}
}
} while (vertex < static_cast<int>(size_ + origin_index));
std::vector<Index> result;
while (!vertices_.empty()) {
while (
!vertices_.empty() &&
(types_.back() == right_lid ||
types_.back() == left_lid)
) {
vertices_.pop_back();
types_.pop_back();
}
if (vertices_.empty()) break;
result.push_back(
static_cast<Index>(vertices_.back()) % size_
);
vertices_.pop_back();
types_.pop_back();
}
return result;
}
private:
static constexpr int not_saved = -1;
const std::vector<Point<T>>& polygon_;
long double eps_;
Index size_;
const Point<T>* origin_ = nullptr;
std::vector<int> vertices_;
std::vector<VertexType> types_;
int left_lid_index_ = not_saved;
int right_lid_index_ = not_saved;
const Point<T>& point(int index) const {
int reduced = index % static_cast<int>(size_);
if (reduced < 0) reduced += static_cast<int>(size_);
return polygon_[static_cast<Index>(reduced)];
}
const Point<T>& top() const { return point(vertices_.back()); }
Homogeneous line(int first, int second) const {
return Homogeneous::line(point(first), point(second));
}
bool left(
const Homogeneous& line_value,
const Point<T>& point_value
) const {
return
predicate_sign<T>(line_value.side(point_value), eps_) > 0;
}
bool right(
const Homogeneous& line_value,
const Point<T>& point_value
) const {
return
predicate_sign<T>(line_value.side(point_value), eps_) < 0;
}
bool clockwise(
const Homogeneous& first,
const Homogeneous& second,
const Homogeneous& third
) const {
return predicate_sign<T>(
Homogeneous::determinant(first, second, third), eps_
) < 0;
}
void push(int index, VertexType type) {
vertices_.push_back(index);
types_.push_back(type);
}
void save_lid() {
if (types_.back() == left_wall) {
vertices_.pop_back();
types_.pop_back();
}
left_lid_index_ = vertices_.back();
vertices_.pop_back();
types_.pop_back();
if (types_.back() == right_lid) {
right_lid_index_ = vertices_.back();
vertices_.pop_back();
types_.pop_back();
} else {
right_lid_index_ = not_saved;
}
}
void restore_lid() {
if (right_lid_index_ != not_saved) {
push(right_lid_index_, right_lid);
}
push(left_lid_index_, left_lid);
}
int exit_right_bay(
int vertex,
const Point<T>& bottom,
const Homogeneous& lid
) const {
int winding = 0;
const Homogeneous mouth = Homogeneous::line(*origin_, bottom);
bool current_left = false;
while (++vertex < 3 * static_cast<int>(size_)) {
const bool last_left = current_left;
current_left = left(mouth, point(vertex));
if (
current_left != last_left &&
(orientation(
point(vertex - 1),
point(vertex),
*origin_,
eps_
) > 0) == current_left
) {
if (!current_left) {
--winding;
} else if (winding++ == 0) {
const Homogeneous edge = line(vertex - 1, vertex);
if (
left(edge, bottom) &&
!clockwise(mouth, edge, lid)
) {
return vertex - 1;
}
}
}
}
return vertex;
}
int exit_left_bay(
int vertex,
const Point<T>& bottom,
const Homogeneous& lid
) const {
int winding = 0;
const Homogeneous mouth = Homogeneous::line(*origin_, bottom);
bool current_right = false;
while (++vertex < 3 * static_cast<int>(size_)) {
const bool last_right = current_right;
current_right = right(mouth, point(vertex));
if (
current_right != last_right &&
(orientation(
point(vertex - 1),
point(vertex),
*origin_,
eps_
) < 0) == current_right
) {
if (!current_right) {
++winding;
} else if (winding-- == 0) {
const Homogeneous edge = line(vertex - 1, vertex);
if (
right(edge, bottom) &&
!clockwise(mouth, edge, lid)
) {
return vertex - 1;
}
}
}
}
return vertex;
}
};
public:
KeilSnoeyinkDecomposition(
const std::vector<Point<T>>& polygon,
long double eps
)
: polygon_(polygon),
eps_(eps),
size_(polygon.size()),
reflex_(size_, false),
remapped_(size_) {}
std::optional<std::vector<Diagonal>> run() {
initialize_reflex_vertices();
allocate_states();
initialize_visibility();
initialize_base_subproblems();
for (int span = 3; span < static_cast<int>(size_); ++span) {
for (const int first : reflex_vertices_) {
const int last = first + span;
if (last >= static_cast<int>(size_)) break;
if (!visible(first, last)) continue;
initialize_pairs(first, last);
if (reflex_[last]) {
for (int middle = first + 1; middle < last; ++middle) {
type_a(first, middle, last);
}
} else {
for (const int middle : reflex_vertices_) {
if (middle <= first) continue;
if (middle >= last - 1) break;
type_a(first, middle, last);
}
type_a(first, last - 1, last);
}
}
for (const int last : reflex_vertices_) {
if (last < span) continue;
const int first = last - span;
if (reflex_[first] || !visible(first, last)) continue;
initialize_pairs(first, last);
type_b(first, first + 1, last);
for (const int middle : reflex_vertices_) {
if (middle <= first + 1) continue;
if (middle >= last) break;
type_b(first, middle, last);
}
}
}
if (weight(0, static_cast<int>(size_) - 1) >= infinity) {
return std::nullopt;
}
std::vector<Diagonal> diagonals;
const PairDeque& root = state(0, static_cast<int>(size_) - 1).pairs;
if (root.pairs.empty()) return std::nullopt;
flatten(root.pairs.front().solution, diagonals);
for (Diagonal& diagonal : diagonals) {
if (diagonal.second < diagonal.first) {
std::swap(diagonal.first, diagonal.second);
}
}
std::sort(
diagonals.begin(),
diagonals.end(),
[](const Diagonal& first, const Diagonal& second) {
return std::pair(first.first, first.second) <
std::pair(second.first, second.second);
}
);
diagonals.erase(
std::unique(
diagonals.begin(),
diagonals.end(),
[](const Diagonal& first, const Diagonal& second) {
return
first.first == second.first &&
first.second == second.second;
}
),
diagonals.end()
);
if (
static_cast<int>(diagonals.size()) !=
weight(0, static_cast<int>(size_) - 1)
) {
return std::nullopt;
}
return diagonals;
}
private:
const std::vector<Point<T>>& polygon_;
long double eps_;
Index size_;
std::vector<bool> reflex_;
std::vector<int> reflex_vertices_;
std::vector<int> remapped_;
std::vector<std::vector<State>> states_;
void initialize_reflex_vertices() {
for (int index = 0; index < static_cast<int>(size_); ++index) {
reflex_[index] = index == 0 || orientation(
polygon_[(index + static_cast<int>(size_) - 1) % size_],
polygon_[index],
polygon_[(index + 1) % size_],
eps_
) < 0;
if (reflex_[index]) reflex_vertices_.push_back(index);
}
}
void allocate_states() {
int next = 0;
for (const int index : reflex_vertices_) remapped_[index] = next++;
for (int index = 0; index < static_cast<int>(size_); ++index) {
if (!reflex_[index]) remapped_[index] = next++;
}
const int reflex_count = static_cast<int>(reflex_vertices_.size());
states_.resize(size_);
for (int index = 0; index < static_cast<int>(size_); ++index) {
states_[index].resize(reflex_[index] ? size_ : reflex_count);
}
}
State& state(int first, int last) {
assert(first <= last);
assert(reflex_[first] || reflex_[last]);
return states_[first][remapped_[last]];
}
const State& state(int first, int last) const {
assert(first <= last);
assert(reflex_[first] || reflex_[last]);
return states_[first][remapped_[last]];
}
int weight(int first, int last) const {
return state(first, last).weight;
}
bool visible(int first, int last) const {
return weight(first, last) < bad;
}
void initialize_visibility() {
VisibilityPolygon visibility_polygon(polygon_, eps_);
for (const int reflex : reflex_vertices_) {
for (const Index visible_vertex :
visibility_polygon.build(static_cast<Index>(reflex))) {
int first = reflex;
int last = static_cast<int>(visible_vertex);
if (last < first) std::swap(first, last);
if (first == last) continue;
state(first, last).weight = infinity;
}
}
// The visibility scan treats the closing edge like every other
// boundary edge, but make the DP convention explicit.
state(0, static_cast<int>(size_) - 1).weight = infinity;
}
void initialize_base_subproblems() {
const int size = static_cast<int>(size_);
for (const int index : reflex_vertices_) {
if (index + 1 < size) state(index, index + 1).weight = 0;
if (index > 0) state(index - 1, index).weight = 0;
if (index + 2 < size && visible(index, index + 2)) {
State& base = state(index, index + 2);
base.weight = 0;
base.pairs.clear();
base.pairs.push(index + 1, index + 1, nullptr);
}
if (index >= 2 && visible(index - 2, index)) {
State& base = state(index - 2, index);
base.weight = 0;
base.pairs.clear();
base.pairs.push(index - 1, index - 1, nullptr);
}
}
}
void initialize_pairs(int first, int last) {
state(first, last).pairs.clear();
}
void update(
int first,
int last,
int new_weight,
int pair_first,
int pair_second,
NodePointer solution
) {
State& subproblem = state(first, last);
if (new_weight > subproblem.weight) return;
if (new_weight < subproblem.weight) {
subproblem.weight = new_weight;
subproblem.pairs.clear();
}
subproblem.pairs.push_narrow(
pair_first, pair_second, std::move(solution)
);
}
NodePointer concatenate(NodePointer first, NodePointer second) const {
if (!first) return second;
if (!second) return first;
SolutionNode node;
node.first = std::move(first);
node.second = std::move(second);
return std::make_shared<const SolutionNode>(std::move(node));
}
NodePointer append_diagonal(
NodePointer solution,
int first,
int last
) const {
SolutionNode leaf;
leaf.diagonal = Diagonal{
static_cast<Index>(first), static_cast<Index>(last)
};
leaf.is_diagonal = true;
return concatenate(
std::move(solution),
std::make_shared<const SolutionNode>(std::move(leaf))
);
}
NodePointer any_solution(int first, int last) const {
const PairDeque& pairs = state(first, last).pairs;
assert(!pairs.pairs.empty());
return pairs.pairs.front().solution;
}
void type_a(int first, int middle, int last) {
if (!visible(first, middle)) return;
int top = middle;
int new_weight = weight(first, middle);
NodePointer solution;
if (last - middle > 1) {
if (!visible(middle, last)) return;
new_weight += weight(middle, last) + 1;
solution = any_solution(middle, last);
solution = append_diagonal(
std::move(solution), middle, last
);
}
bool use_first_diagonal = false;
if (middle - first > 1) {
PairDeque& pairs = state(first, middle).pairs;
if (pairs.empty_back()) return;
if (
orientation(
polygon_[last],
polygon_[middle],
polygon_[pairs.back_pair().second],
eps_
) <= 0
) {
while (
pairs.more_back() &&
orientation(
polygon_[last],
polygon_[middle],
polygon_[pairs.under_back().second],
eps_
) <= 0
) {
pairs.pop_back();
}
if (
!pairs.empty_back() &&
orientation(
polygon_[last],
polygon_[first],
polygon_[pairs.back_pair().first],
eps_
) >= 0
) {
top = pairs.back_pair().first;
} else {
++new_weight;
use_first_diagonal = true;
}
} else {
++new_weight;
use_first_diagonal = true;
}
solution = concatenate(
pairs.back_pair().solution, std::move(solution)
);
if (use_first_diagonal) {
solution = append_diagonal(
std::move(solution), first, middle
);
}
}
update(
first,
last,
new_weight,
top,
middle,
std::move(solution)
);
}
void type_b(int first, int middle, int last) {
if (!visible(middle, last)) return;
int top = middle;
int new_weight = weight(middle, last);
NodePointer solution;
if (middle - first > 1) {
if (!visible(first, middle)) return;
new_weight += weight(first, middle) + 1;
solution = any_solution(first, middle);
solution = append_diagonal(
std::move(solution), first, middle
);
}
bool use_second_diagonal = false;
if (last - middle > 1) {
PairDeque& pairs = state(middle, last).pairs;
if (pairs.empty_front()) return;
if (
orientation(
polygon_[first],
polygon_[middle],
polygon_[pairs.front_pair().first],
eps_
) >= 0
) {
while (
pairs.more_front() &&
orientation(
polygon_[first],
polygon_[middle],
polygon_[pairs.under_front().first],
eps_
) >= 0
) {
pairs.pop_front();
}
if (
!pairs.empty_front() &&
orientation(
polygon_[first],
polygon_[last],
polygon_[pairs.front_pair().second],
eps_
) <= 0
) {
top = pairs.front_pair().second;
} else {
++new_weight;
use_second_diagonal = true;
}
} else {
++new_weight;
use_second_diagonal = true;
}
solution = concatenate(
std::move(solution), pairs.front_pair().solution
);
if (use_second_diagonal) {
solution = append_diagonal(
std::move(solution), middle, last
);
}
}
update(
first,
last,
new_weight,
middle,
top,
std::move(solution)
);
}
void flatten(
const NodePointer& node,
std::vector<Diagonal>& diagonals
) const {
if (!node) return;
if (node->is_diagonal) {
diagonals.push_back(node->diagonal);
return;
}
flatten(node->first, diagonals);
flatten(node->second, diagonals);
}
};
template <Coordinate T>
std::optional<std::vector<IndexPolygon>> faces_from_diagonals(
const std::vector<Point<T>>& polygon,
const std::vector<Diagonal>& diagonals,
long double eps
) {
const Index size = polygon.size();
std::vector<std::vector<Index>> adjacent(size);
auto add_edge = [&](Index first, Index second) {
adjacent[first].push_back(second);
adjacent[second].push_back(first);
};
for (Index index = 0; index < size; ++index) {
add_edge(index, (index + 1) % size);
}
for (const Diagonal& diagonal : diagonals) {
if (
diagonal.first >= size || diagonal.second >= size ||
diagonal.first == diagonal.second
) {
return std::nullopt;
}
add_edge(diagonal.first, diagonal.second);
}
for (Index origin = 0; origin < size; ++origin) {
auto upper_half = [&](Index target) {
const auto direction = polygon[target] - polygon[origin];
const int y_sign = sign<T>(direction.y, eps);
return
y_sign > 0 ||
(y_sign == 0 && sign<T>(direction.x, eps) >= 0);
};
std::sort(
adjacent[origin].begin(),
adjacent[origin].end(),
[&](Index first, Index second) {
const bool first_upper = upper_half(first);
const bool second_upper = upper_half(second);
if (first_upper != second_upper) return first_upper;
const auto first_direction =
polygon[first] - polygon[origin];
const auto second_direction =
polygon[second] - polygon[origin];
const int turn = sign<T>(
cross(first_direction, second_direction), eps
);
if (turn != 0) return turn > 0;
return
norm2(first_direction) < norm2(second_direction);
}
);
adjacent[origin].erase(
std::unique(
adjacent[origin].begin(), adjacent[origin].end()
),
adjacent[origin].end()
);
}
std::vector<std::vector<bool>> visited(size);
Index directed_edge_count = 0;
for (Index index = 0; index < size; ++index) {
visited[index].assign(adjacent[index].size(), false);
directed_edge_count += adjacent[index].size();
}
std::vector<IndexPolygon> result;
for (Index start = 0; start < size; ++start) {
for (Index start_position = 0;
start_position < adjacent[start].size();
++start_position) {
if (visited[start][start_position]) continue;
IndexPolygon face;
Index current = start;
Index position = start_position;
for (Index guard = 0;; ++guard) {
if (guard > directed_edge_count) return std::nullopt;
if (visited[current][position]) {
if (current != start || position != start_position) {
return std::nullopt;
}
break;
}
visited[current][position] = true;
face.push_back(current);
const Index next = adjacent[current][position];
const auto reverse = std::find(
adjacent[next].begin(), adjacent[next].end(), current
);
if (reverse == adjacent[next].end()) return std::nullopt;
const Index reverse_position =
static_cast<Index>(reverse - adjacent[next].begin());
current = next;
position =
(reverse_position + adjacent[current].size() - 1) %
adjacent[current].size();
}
if (face.size() < 3) continue;
std::vector<Point<T>> points;
points.reserve(face.size());
for (const Index index : face) points.push_back(polygon[index]);
if (sign<T>(polygon_area2(points), eps) > 0) {
if (!is_weakly_convex(face, polygon, eps)) {
return std::nullopt;
}
result.push_back(std::move(face));
}
}
}
if (result.size() != diagonals.size() + 1) return std::nullopt;
return result;
}
} // namespace convex_decomposition_detail
template <Coordinate T>
std::optional<std::vector<std::vector<Point<T>>>>
convex_decomposition(
std::vector<Point<T>> polygon,
long double eps = 1e-12L
) {
auto prepared = convex_decomposition_detail::prepare_polygon(
std::move(polygon), eps
);
if (!prepared.has_value()) return std::nullopt;
polygon = std::move(*prepared);
if (convex_decomposition_detail::is_weakly_convex(polygon, eps)) {
return std::vector<std::vector<Point<T>>>{std::move(polygon)};
}
auto triangulation =
convex_decomposition_detail::triangulate_indices(polygon, eps);
if (!triangulation.has_value()) return std::nullopt;
using convex_decomposition_detail::Index;
using convex_decomposition_detail::IndexPolygon;
const Index triangle_count = triangulation->size();
const Index absent = std::numeric_limits<Index>::max();
struct Owners {
Index first = std::numeric_limits<Index>::max();
Index second = std::numeric_limits<Index>::max();
};
std::map<std::pair<Index, Index>, Owners> edge_owners;
for (Index triangle = 0; triangle < triangle_count; ++triangle) {
for (Index edge = 0; edge < 3; ++edge) {
Index first = (*triangulation)[triangle][edge];
Index second = (*triangulation)[triangle][(edge + 1) % 3];
if (second < first) std::swap(first, second);
Owners& owners = edge_owners[{first, second}];
if (owners.first == absent) {
owners.first = triangle;
} else {
owners.second = triangle;
}
}
}
std::vector<Index> parent(triangle_count);
std::vector<IndexPolygon> pieces = std::move(*triangulation);
for (Index index = 0; index < triangle_count; ++index) {
parent[index] = index;
}
for (const auto& [edge, owners] : edge_owners) {
if (owners.second == absent) continue;
Index first_root =
convex_decomposition_detail::find_root(parent, owners.first);
Index second_root =
convex_decomposition_detail::find_root(parent, owners.second);
if (first_root == second_root) continue;
auto merged = convex_decomposition_detail::merge_across_edge(
pieces[first_root],
pieces[second_root],
edge.first,
edge.second
);
if (
!merged.has_value() ||
!convex_decomposition_detail::is_weakly_convex(
*merged, polygon, eps
)
) {
continue;
}
pieces[first_root] = std::move(*merged);
pieces[second_root].clear();
parent[second_root] = first_root;
}
std::vector<std::vector<Point<T>>> result;
for (Index index = 0; index < triangle_count; ++index) {
if (convex_decomposition_detail::find_root(parent, index) != index) {
continue;
}
result.push_back(convex_decomposition_detail::materialize(
pieces[index], polygon, eps
));
}
return result;
}
template <Coordinate T>
std::optional<std::vector<std::vector<Point<T>>>>
minimum_convex_decomposition(
std::vector<Point<T>> polygon,
long double eps = 1e-12L
) {
auto prepared = convex_decomposition_detail::prepare_minimum_polygon(
std::move(polygon), eps
);
if (!prepared.has_value()) return std::nullopt;
polygon = std::move(*prepared);
if (convex_decomposition_detail::is_weakly_convex(polygon, eps)) {
return std::vector<std::vector<Point<T>>>{std::move(polygon)};
}
std::size_t first_reflex = polygon.size();
for (std::size_t index = 0; index < polygon.size(); ++index) {
if (
orientation(
polygon[(index + polygon.size() - 1) % polygon.size()],
polygon[index],
polygon[(index + 1) % polygon.size()],
eps
) < 0
) {
first_reflex = index;
break;
}
}
assert(first_reflex < polygon.size());
std::rotate(
polygon.begin(), polygon.begin() + first_reflex, polygon.end()
);
std::size_t reflex_count = 0;
for (std::size_t index = 0; index < polygon.size(); ++index) {
if (
orientation(
polygon[(index + polygon.size() - 1) % polygon.size()],
polygon[index],
polygon[(index + 1) % polygon.size()],
eps
) < 0
) {
++reflex_count;
}
}
auto solve = [&]<typename Number>()
-> std::optional<std::vector<std::vector<Point<T>>>> {
std::vector<Point<T>> dynamic_polygon = polygon;
std::vector<convex_decomposition_detail::Index> original_index(
polygon.size()
);
for (std::size_t index = 0; index < polygon.size(); ++index) {
original_index[index] = index;
}
if (
reflex_count > 0 &&
polygon.size() / reflex_count > reflex_count
) {
convex_decomposition_detail::BiasedPolygonReduction<T, Number>
reduction(polygon, eps);
auto reduced = reduction.run();
dynamic_polygon = std::move(reduced.polygon);
original_index = std::move(reduced.original_index);
}
convex_decomposition_detail::KeilSnoeyinkDecomposition<T, Number>
solver(dynamic_polygon, eps);
auto diagonals = solver.run();
if (!diagonals.has_value()) return std::nullopt;
for (auto& diagonal : *diagonals) {
diagonal.first = original_index[diagonal.first];
diagonal.second = original_index[diagonal.second];
}
auto index_decomposition =
convex_decomposition_detail::faces_from_diagonals(
polygon, *diagonals, eps
);
if (!index_decomposition.has_value()) return std::nullopt;
std::vector<std::vector<Point<T>>> result;
result.reserve(index_decomposition->size());
for (const auto& indices : *index_decomposition) {
result.push_back(convex_decomposition_detail::materialize(
indices, polygon, eps
));
}
return result;
};
if constexpr (std::integral<T>) {
using convex_decomposition_detail::ExactPredicateWidth;
switch (
convex_decomposition_detail::select_exact_predicate_width(
polygon
)
) {
case ExactPredicateWidth::Int128:
return solve.template operator()<__int128_t>();
case ExactPredicateWidth::Int256:
return solve.template operator()<
m1une::utilities::Int256
>();
case ExactPredicateWidth::Int512:
return solve.template operator()<
m1une::utilities::Int512
>();
}
assert(false);
return std::nullopt;
} else {
return solve.template operator()<wide_type<T>>();
}
}
} // namespace geometry
} // namespace m1une