Minimum Enclosing Circle
(geometry/minimum_enclosing_circle.hpp)
- View this file on GitHub
- Last update: 2026-10-05 22:23:07+09:00
- Include:
#include "geometry/minimum_enclosing_circle.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.
minimum_enclosing_circle(points) finds the unique smallest circle containing
every input point. It uses the randomized incremental algorithm and returns the
circle together with one support set of at most three original points that
defines it.
Duplicate and collinear points are supported.
Interface
struct MinimumEnclosingCircle {
Circle<long double> circle;
std::vector<int> support;
};
template <Coordinate T>
std::optional<MinimumEnclosingCircle> minimum_enclosing_circle(
const std::vector<Point<T>>& points,
long double eps = 1e-12L
);
Result
| Member | Description |
|---|---|
circle |
The minimum enclosing Circle<long double>. |
support |
One to three sorted original indices whose points define the circle. |
The function returns std::nullopt for an empty input. For one or more points,
support is nonempty. It may omit other input points that also lie on the
boundary. The returned circle has filled == true.
For integral coordinates, containment decisions use signed 128-bit arithmetic;
all predicate intermediates must fit in that type. The returned center and
radius are still long double. For floating-point coordinates, eps controls
containment comparisons relative to the squared radius.
The shuffle is deterministic, so repeated calls on the same input produce the same support set.
Operations
| Function | Description | Complexity |
|---|---|---|
minimum_enclosing_circle(points, eps) |
Returns the smallest circle containing all points without modifying the input. | Expected $O(N)$ time and $O(N)$ memory |
Example
#include "geometry/minimum_enclosing_circle.hpp"
#include <cassert>
#include <vector>
int main() {
using Point = m1une::geometry::Point<long long>;
std::vector<Point> points;
points.emplace_back(0, 0);
points.emplace_back(4, 0);
points.emplace_back(2, 1);
auto answer = m1une::geometry::minimum_enclosing_circle(points);
assert(answer.has_value());
assert(answer->circle.center.x == 2.0L);
assert(answer->circle.center.y == 0.0L);
assert(answer->circle.radius == 2.0L);
}
Depends on
- Circles
(geometry/circle.hpp)
- geometry/detail/floating_predicate.hpp
- Linear Objects
(geometry/linear.hpp)
- 2D Point and Predicates
(geometry/point.hpp)
Required by
Verified with
- verify/geometry/centroid.test.cpp
- verify/geometry/geometry_algorithms.test.cpp
- verify/geometry/minimum_enclosing_circle.test.cpp
- verify/geometry/rational.test.cpp
Code
#ifndef M1UNE_GEOMETRY_MINIMUM_ENCLOSING_CIRCLE_HPP
#define M1UNE_GEOMETRY_MINIMUM_ENCLOSING_CIRCLE_HPP 1
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <numeric>
#include <optional>
#include <type_traits>
#include <vector>
#include "circle.hpp"
namespace m1une {
namespace geometry {
struct MinimumEnclosingCircle {
Circle<long double> circle;
std::vector<int> support;
};
namespace internal {
struct MinimumEnclosingCircleSupport {
int size = 0;
std::array<int, 3> index = {-1, -1, -1};
};
template <Coordinate T>
Circle<long double> circle_from_support(
const std::vector<Point<T>>& points,
const MinimumEnclosingCircleSupport& support
) {
assert(1 <= support.size && support.size <= 3);
Point<long double> first(points[support.index[0]]);
if (support.size == 1) return Circle<long double>{first, 0.0L};
Point<long double> second(points[support.index[1]]);
if (support.size == 2) {
Point<long double> center = (first + second) / 2.0L;
return Circle<long double>{center, distance(first, second) / 2.0L};
}
Point<long double> third(points[support.index[2]]);
long double denominator = 2.0L * cross(first, second, third);
assert(denominator != 0.0L);
long double first_norm = norm2(first);
long double second_norm = norm2(second);
long double third_norm = norm2(third);
Point<long double> center(
(
first_norm * (second.y - third.y) +
second_norm * (third.y - first.y) +
third_norm * (first.y - second.y)
) / denominator,
(
first_norm * (third.x - second.x) +
second_norm * (first.x - third.x) +
third_norm * (second.x - first.x)
) / denominator
);
return Circle<long double>{center, distance(center, first)};
}
template <Coordinate T>
void reduce_collinear_support(
const std::vector<Point<T>>& points,
MinimumEnclosingCircleSupport& support
) {
if (support.size != 3) return;
const Point<T>& first = points[support.index[0]];
const Point<T>& second = points[support.index[1]];
const Point<T>& third = points[support.index[2]];
if (cross(first, second, third) != 0) return;
wide_type<T> first_second = distance2(first, second);
wide_type<T> second_third = distance2(second, third);
wide_type<T> first_third = distance2(first, third);
if (first_second >= second_third && first_second >= first_third) {
support.size = 2;
} else if (second_third >= first_second && second_third >= first_third) {
support.size = 2;
support.index[0] = support.index[1];
support.index[1] = support.index[2];
} else {
support.size = 2;
support.index[1] = support.index[2];
}
}
template <Coordinate T>
bool support_contains(
const std::vector<Point<T>>& points,
const MinimumEnclosingCircleSupport& support,
int point_index,
long double eps
) {
if (support.size == 0) return false;
const Point<T>& point = points[point_index];
const Point<T>& first = points[support.index[0]];
if (support.size == 1) return point == first;
if constexpr (ExactCoordinate<T>) {
using W = wide_type<T>;
const Point<T>& second = points[support.index[1]];
if (support.size == 2) {
W ux = W(2) * W(point.x) - W(first.x) - W(second.x);
W uy = W(2) * W(point.y) - W(first.y) - W(second.y);
W dx = W(first.x) - W(second.x);
W dy = W(first.y) - W(second.y);
return ux * ux + uy * uy <= dx * dx + dy * dy;
}
const Point<T>& third = points[support.index[2]];
W ax = W(first.x) - W(point.x);
W ay = W(first.y) - W(point.y);
W bx = W(second.x) - W(point.x);
W by = W(second.y) - W(point.y);
W cx = W(third.x) - W(point.x);
W cy = W(third.y) - W(point.y);
W a_norm = ax * ax + ay * ay;
W b_norm = bx * bx + by * by;
W c_norm = cx * cx + cy * cy;
W determinant =
ax * (by * c_norm - cy * b_norm) -
ay * (bx * c_norm - cx * b_norm) +
a_norm * (bx * cy - by * cx);
W direction = cross(first, second, third);
return direction > 0 ? determinant >= 0 : determinant <= 0;
} else {
Circle<long double> circle = circle_from_support(points, support);
long double squared = distance2(Point<long double>(point), circle.center);
long double radius_squared = circle.radius * circle.radius;
long double tolerance = eps * std::max(1.0L, radius_squared);
return squared <= radius_squared + tolerance;
}
}
inline std::uint64_t minimum_enclosing_circle_random(std::uint64_t& state) {
state ^= state << 7;
state ^= state >> 9;
return state;
}
} // namespace internal
// Returns the unique minimum circle containing every input point.
template <Coordinate T>
std::optional<MinimumEnclosingCircle> minimum_enclosing_circle(
const std::vector<Point<T>>& points,
long double eps = 1e-12L
) {
if (points.empty()) return std::nullopt;
assert(eps >= 0.0L);
const int n = int(points.size());
std::vector<int> order(n);
std::iota(order.begin(), order.end(), 0);
std::uint64_t state = 0x9e3779b97f4a7c15ULL;
for (int index = n - 1; index >= 1; index--) {
int other = int(
internal::minimum_enclosing_circle_random(state) %
std::uint64_t(index + 1)
);
std::swap(order[index], order[other]);
}
internal::MinimumEnclosingCircleSupport support;
for (int first = 0; first < n; first++) {
if (internal::support_contains(points, support, order[first], eps)) {
continue;
}
support.size = 1;
support.index[0] = order[first];
for (int second = 0; second < first; second++) {
if (internal::support_contains(points, support, order[second], eps)) {
continue;
}
support.size = 2;
support.index[0] = order[first];
support.index[1] = order[second];
for (int third = 0; third < second; third++) {
if (internal::support_contains(
points,
support,
order[third],
eps
)) {
continue;
}
support.size = 3;
support.index[0] = order[first];
support.index[1] = order[second];
support.index[2] = order[third];
internal::reduce_collinear_support(points, support);
}
}
}
MinimumEnclosingCircle result;
result.circle = internal::circle_from_support(points, support);
result.support.assign(
support.index.begin(),
support.index.begin() + support.size
);
std::sort(result.support.begin(), result.support.end());
return result;
}
} // namespace geometry
} // namespace m1une
#endif // M1UNE_GEOMETRY_MINIMUM_ENCLOSING_CIRCLE_HPP#line 1 "geometry/minimum_enclosing_circle.hpp"
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <numeric>
#include <optional>
#include <type_traits>
#include <vector>
#line 1 "geometry/circle.hpp"
#line 8 "geometry/circle.hpp"
#include <cstddef>
#include <numbers>
#line 13 "geometry/circle.hpp"
#line 1 "geometry/linear.hpp"
#line 7 "geometry/linear.hpp"
#line 1 "geometry/point.hpp"
#line 5 "geometry/point.hpp"
#include <concepts>
#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 15 "geometry/minimum_enclosing_circle.hpp"
namespace m1une {
namespace geometry {
struct MinimumEnclosingCircle {
Circle<long double> circle;
std::vector<int> support;
};
namespace internal {
struct MinimumEnclosingCircleSupport {
int size = 0;
std::array<int, 3> index = {-1, -1, -1};
};
template <Coordinate T>
Circle<long double> circle_from_support(
const std::vector<Point<T>>& points,
const MinimumEnclosingCircleSupport& support
) {
assert(1 <= support.size && support.size <= 3);
Point<long double> first(points[support.index[0]]);
if (support.size == 1) return Circle<long double>{first, 0.0L};
Point<long double> second(points[support.index[1]]);
if (support.size == 2) {
Point<long double> center = (first + second) / 2.0L;
return Circle<long double>{center, distance(first, second) / 2.0L};
}
Point<long double> third(points[support.index[2]]);
long double denominator = 2.0L * cross(first, second, third);
assert(denominator != 0.0L);
long double first_norm = norm2(first);
long double second_norm = norm2(second);
long double third_norm = norm2(third);
Point<long double> center(
(
first_norm * (second.y - third.y) +
second_norm * (third.y - first.y) +
third_norm * (first.y - second.y)
) / denominator,
(
first_norm * (third.x - second.x) +
second_norm * (first.x - third.x) +
third_norm * (second.x - first.x)
) / denominator
);
return Circle<long double>{center, distance(center, first)};
}
template <Coordinate T>
void reduce_collinear_support(
const std::vector<Point<T>>& points,
MinimumEnclosingCircleSupport& support
) {
if (support.size != 3) return;
const Point<T>& first = points[support.index[0]];
const Point<T>& second = points[support.index[1]];
const Point<T>& third = points[support.index[2]];
if (cross(first, second, third) != 0) return;
wide_type<T> first_second = distance2(first, second);
wide_type<T> second_third = distance2(second, third);
wide_type<T> first_third = distance2(first, third);
if (first_second >= second_third && first_second >= first_third) {
support.size = 2;
} else if (second_third >= first_second && second_third >= first_third) {
support.size = 2;
support.index[0] = support.index[1];
support.index[1] = support.index[2];
} else {
support.size = 2;
support.index[1] = support.index[2];
}
}
template <Coordinate T>
bool support_contains(
const std::vector<Point<T>>& points,
const MinimumEnclosingCircleSupport& support,
int point_index,
long double eps
) {
if (support.size == 0) return false;
const Point<T>& point = points[point_index];
const Point<T>& first = points[support.index[0]];
if (support.size == 1) return point == first;
if constexpr (ExactCoordinate<T>) {
using W = wide_type<T>;
const Point<T>& second = points[support.index[1]];
if (support.size == 2) {
W ux = W(2) * W(point.x) - W(first.x) - W(second.x);
W uy = W(2) * W(point.y) - W(first.y) - W(second.y);
W dx = W(first.x) - W(second.x);
W dy = W(first.y) - W(second.y);
return ux * ux + uy * uy <= dx * dx + dy * dy;
}
const Point<T>& third = points[support.index[2]];
W ax = W(first.x) - W(point.x);
W ay = W(first.y) - W(point.y);
W bx = W(second.x) - W(point.x);
W by = W(second.y) - W(point.y);
W cx = W(third.x) - W(point.x);
W cy = W(third.y) - W(point.y);
W a_norm = ax * ax + ay * ay;
W b_norm = bx * bx + by * by;
W c_norm = cx * cx + cy * cy;
W determinant =
ax * (by * c_norm - cy * b_norm) -
ay * (bx * c_norm - cx * b_norm) +
a_norm * (bx * cy - by * cx);
W direction = cross(first, second, third);
return direction > 0 ? determinant >= 0 : determinant <= 0;
} else {
Circle<long double> circle = circle_from_support(points, support);
long double squared = distance2(Point<long double>(point), circle.center);
long double radius_squared = circle.radius * circle.radius;
long double tolerance = eps * std::max(1.0L, radius_squared);
return squared <= radius_squared + tolerance;
}
}
inline std::uint64_t minimum_enclosing_circle_random(std::uint64_t& state) {
state ^= state << 7;
state ^= state >> 9;
return state;
}
} // namespace internal
// Returns the unique minimum circle containing every input point.
template <Coordinate T>
std::optional<MinimumEnclosingCircle> minimum_enclosing_circle(
const std::vector<Point<T>>& points,
long double eps = 1e-12L
) {
if (points.empty()) return std::nullopt;
assert(eps >= 0.0L);
const int n = int(points.size());
std::vector<int> order(n);
std::iota(order.begin(), order.end(), 0);
std::uint64_t state = 0x9e3779b97f4a7c15ULL;
for (int index = n - 1; index >= 1; index--) {
int other = int(
internal::minimum_enclosing_circle_random(state) %
std::uint64_t(index + 1)
);
std::swap(order[index], order[other]);
}
internal::MinimumEnclosingCircleSupport support;
for (int first = 0; first < n; first++) {
if (internal::support_contains(points, support, order[first], eps)) {
continue;
}
support.size = 1;
support.index[0] = order[first];
for (int second = 0; second < first; second++) {
if (internal::support_contains(points, support, order[second], eps)) {
continue;
}
support.size = 2;
support.index[0] = order[first];
support.index[1] = order[second];
for (int third = 0; third < second; third++) {
if (internal::support_contains(
points,
support,
order[third],
eps
)) {
continue;
}
support.size = 3;
support.index[0] = order[first];
support.index[1] = order[second];
support.index[2] = order[third];
internal::reduce_collinear_support(points, support);
}
}
}
MinimumEnclosingCircle result;
result.circle = internal::circle_from_support(points, support);
result.support.assign(
support.index.begin(),
support.index.begin() + support.size
);
std::sort(result.support.begin(), result.support.end());
return result;
}
} // namespace geometry
} // namespace m1une