description: Return the minimum circle and one boundary support set. Randomly ordering the points makes the expected number of constraint rebuilds linear: an outside point must belong to the new boundary, reducing the remaining problem successively to circles through one, two, then three fixed points.¶
minimum_enclosing_circle.hpp¶
Return the minimum circle and one boundary support set. Randomly ordering the points makes the expected number of constraint rebuilds linear: an outside point must belong to the new boundary, reducing the remaining problem successively to circles through one, two, then three fixed points.
Verified by minimum_enclosing_circle.
求覆盖所有点的最小圆及其边界支撑点；适合最小覆盖半径问题。
Implementation¶
#ifndef NOYA_MINIMUM_ENCLOSING_CIRCLE_HPP
#define NOYA_MINIMUM_ENCLOSING_CIRCLE_HPP 1
/// @complexity Time: Expected O(n).
/// Space: O(n) shuffled points.
#include "noya/geometry_base.hpp"
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <random>
#include <type_traits>
#include <vector>
namespace noya {
template <class T> struct minimum_enclosing_circle_result {
circle<long double> value;
std::array<point<T>, 3> support{};
int support_size = 0;
};
namespace minimum_enclosing_circle_internal {
inline circle<long double> from_diameter(const point<long double> &first,
const point<long double> &second) {
point<long double> center = (first + second) / 2.0L;
return {center, distance(center, first)};
}
inline circle<long double> from_three(const point<long double> &first,
const point<long double> &second,
const point<long double> &third) {
long double denominator = 2.0L * cross(first, second, third);
if (std::abs(denominator) <= 1e-24L) {
circle<long double> result = from_diameter(first, second);
circle<long double> candidate = from_diameter(first, third);
if (candidate.radius > result.radius) {
result = candidate;
}
candidate = from_diameter(second, third);
if (candidate.radius > result.radius) {
result = candidate;
}
return result;
}
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 {center, distance(center, first)};
}
inline bool contains(const circle<long double> &value,
const point<long double> &candidate) {
long double squared_radius = value.radius * value.radius;
long double tolerance = 1e-12L * std::max(1.0L, squared_radius);
return norm2(candidate - value.center) <= squared_radius + tolerance;
}
} // namespace minimum_enclosing_circle_internal
/// @brief Return the minimum circle and one boundary support set. Randomly
/// ordering the points makes the expected number of constraint rebuilds
/// linear: an outside point must belong to the new boundary, reducing the
/// remaining problem successively to circles through one, two, then three
/// fixed points.
template <class T>
minimum_enclosing_circle_result<T>
minimum_enclosing_circle_with_support(const std::vector<point<T>> &input,
std::uint64_t seed = 712367821ULL) {
assert(!input.empty());
struct indexed_point {
point<long double> value;
point<T> original;
};
std::vector<indexed_point> points;
points.reserve(input.size());
for (const point<T> &value : input) {
points.push_back({{static_cast<long double>(value.x),
static_cast<long double>(value.y)},
value});
}
std::mt19937_64 random(seed);
std::shuffle(points.begin(), points.end(), random);
using namespace minimum_enclosing_circle_internal;
minimum_enclosing_circle_result<T> answer;
answer.value = {points[0].value, 0};
answer.support[0] = points[0].original;
answer.support_size = 1;
auto answer_contains = [&](const indexed_point &candidate) {
if constexpr (!std::is_integral_v<T>) {
return contains(answer.value, candidate.value);
} else {
using wide = __int128;
const auto &a = answer.support[0];
const auto &p = candidate.original;
if (answer.support_size == 1) {
return p == a;
}
const auto &b = answer.support[1];
if (answer.support_size == 2) {
wide dx = wide(2) * p.x - a.x - b.x;
wide dy = wide(2) * p.y - a.y - b.y;
wide diameter2 = wide(a.x - b.x) * (a.x - b.x) +
wide(a.y - b.y) * (a.y - b.y);
return dx * dx + dy * dy <= diameter2;
}
const auto &c = answer.support[2];
wide bx = wide(b.x) - a.x;
wide by = wide(b.y) - a.y;
wide cx = wide(c.x) - a.x;
wide cy = wide(c.y) - a.y;
wide px = wide(p.x) - a.x;
wide py = wide(p.y) - a.y;
wide determinant =
(px * px + py * py) * (bx * cy - by * cx) -
(bx * bx + by * by) * (px * cy - py * cx) +
(cx * cx + cy * cy) * (px * by - py * bx);
return determinant * (bx * cy - by * cx) <= 0;
}
};
for (int i = 1; i < int(points.size()); i++) {
if (answer_contains(points[i])) {
continue;
}
answer.value = {points[i].value, 0};
answer.support[0] = points[i].original;
answer.support_size = 1;
for (int j = 0; j < i; j++) {
if (answer_contains(points[j])) {
continue;
}
answer.value = from_diameter(points[i].value, points[j].value);
answer.support[0] = points[i].original;
answer.support[1] = points[j].original;
answer.support_size = 2;
for (int k = 0; k < j; k++) {
if (!answer_contains(points[k])) {
long double twice_area =
cross(points[i].value, points[j].value, points[k].value);
if (std::abs(twice_area) <= 1e-24L) {
int first = i;
int second = j;
long double best = norm2(points[i].value - points[j].value);
for (auto [left, right] :
{std::pair{i, k}, std::pair{j, k}}) {
long double candidate =
norm2(points[left].value - points[right].value);
if (candidate > best) {
best = candidate;
first = left;
second = right;
}
}
answer.value =
from_diameter(points[first].value, points[second].value);
answer.support[0] = points[first].original;
answer.support[1] = points[second].original;
answer.support_size = 2;
} else {
answer.value =
from_three(points[i].value, points[j].value, points[k].value);
answer.support[0] = points[i].original;
answer.support[1] = points[j].original;
answer.support[2] = points[k].original;
answer.support_size = 3;
}
}
}
}
}
return answer;
}
/// @brief Return the minimum circle containing all points in expected O(n).
template <class T>
circle<long double>
minimum_enclosing_circle(const std::vector<point<T>> &input,
std::uint64_t seed = 712367821ULL) {
return minimum_enclosing_circle_with_support(input, seed).value;
}
/// @brief Mark exactly which integral points lie on the minimum circle. Once
/// the randomized construction supplies two or three support points, equality
/// with their circle is tested after clearing all denominators: a doubled
/// midpoint equation for a diameter, or the integer circumcircle determinant
/// for three non-collinear points.
template <class T>
std::vector<bool>
minimum_enclosing_circle_boundary(const std::vector<point<T>> &input,
std::uint64_t seed = 712367821ULL) {
static_assert(std::is_integral_v<T>);
auto answer = minimum_enclosing_circle_with_support(input, seed);
std::vector<bool> result(input.size());
using wide = __int128;
const auto &a = answer.support[0];
if (answer.support_size == 1) {
for (int i = 0; i < int(input.size()); i++) {
result[i] = input[i] == a;
}
} else if (answer.support_size == 2) {
const auto &b = answer.support[1];
wide radius4 = wide(a.x - b.x) * (a.x - b.x) +
wide(a.y - b.y) * (a.y - b.y);
for (int i = 0; i < int(input.size()); i++) {
wide dx = wide(2) * input[i].x - a.x - b.x;
wide dy = wide(2) * input[i].y - a.y - b.y;
result[i] = dx * dx + dy * dy == radius4;
}
} else {
const auto &b = answer.support[1];
const auto &c = answer.support[2];
wide bx = wide(b.x) - a.x;
wide by = wide(b.y) - a.y;
wide cx = wide(c.x) - a.x;
wide cy = wide(c.y) - a.y;
wide bc = bx * cy - by * cx;
wide bn = bx * bx + by * by;
wide cn = cx * cx + cy * cy;
for (int i = 0; i < int(input.size()); i++) {
wide px = wide(input[i].x) - a.x;
wide py = wide(input[i].y) - a.y;
wide pn = px * px + py * py;
result[i] = pn * bc - bn * (px * cy - py * cx) +
cn * (px * by - py * bx) ==
0;
}
}
return result;
}
} // namespace noya
#endif // NOYA_MINIMUM_ENCLOSING_CIRCLE_HPP
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <optional>
#include <random>
#include <type_traits>
#include <vector>
/// @complexity Time: Expected O(n).
/// Space: O(n) shuffled points.
/// @complexity Time: O(1) per primitive; O(n log n) for convex hull.
/// Space: O(1) per primitive and O(n) for hull construction.
namespace noya {
/// @brief Two-dimensional point with vector arithmetic and lexicographic order.
template <class T> struct point {
T x{};
T y{};
point() = default;
point(T x_, T y_) : x(x_), y(y_) {}
point &operator+=(const point &other) {
x += other.x;
y += other.y;
return *this;
}
point &operator-=(const point &other) {
x -= other.x;
y -= other.y;
return *this;
}
point &operator*=(const T &scale) {
x *= scale;
y *= scale;
return *this;
}
point &operator/=(const T &scale) {
x /= scale;
y /= scale;
return *this;
}
friend point operator+(point left, const point &right) {
return left += right;
}
friend point operator-(point left, const point &right) {
return left -= right;
}
friend point operator*(point value, const T &scale) { return value *= scale; }
friend point operator*(const T &scale, point value) { return value *= scale; }
friend point operator/(point value, const T &scale) { return value /= scale; }
friend bool operator==(const point &, const point &) = default;
friend bool operator<(const point &left, const point &right) {
return left.x < right.x || (left.x == right.x && left.y < right.y);
}
};
/// @brief Circle represented by a center and a nonnegative radius.
template <class Real> struct circle {
point<Real> center;
Real radius{};
};
/// @brief Return the dot product of two vectors.
template <class T> T dot(const point<T> &a, const point<T> &b) {
return a.x * b.x + a.y * b.y;
}
/// @brief Return the signed cross product of two vectors.
template <class T> T cross(const point<T> &a, const point<T> &b) {
return a.x * b.y - a.y * b.x;
}
/// @brief Return cross(a - origin, b - origin).
template <class T>
T cross(const point<T> &origin, const point<T> &a, const point<T> &b) {
return cross(a - origin, b - origin);
}
/// @brief Return the squared Euclidean norm.
template <class T> T norm2(const point<T> &value) { return dot(value, value); }
/// @brief Compare a value with zero using an optional absolute tolerance.
template <class T> int sign(const T &value, const T &epsilon = T{}) {
return (value > epsilon) - (value < -epsilon);
}
/// @brief Return -1, 0, or 1 for a clockwise, collinear, or counter-clockwise
/// turn.
template <class T>
int orientation(const point<T> &a, const point<T> &b, const point<T> &c,
const T &epsilon = T{}) {
return sign(cross(a, b, c), epsilon);
}
/// @brief Test whether p lies on the closed segment [a, b].
template <class T>
bool on_segment(const point<T> &p, const point<T> &a, const point<T> &b,
const T &epsilon = T{}) {
if (orientation(a, b, p, epsilon) != 0) {
return false;
}
return std::min(a.x, b.x) - epsilon <= p.x &&
p.x <= std::max(a.x, b.x) + epsilon &&
std::min(a.y, b.y) - epsilon <= p.y &&
p.y <= std::max(a.y, b.y) + epsilon;
}
/// @brief Test whether the closed segments [a, b] and [c, d] intersect.
template <class T>
bool segments_intersect(const point<T> &a, const point<T> &b, const point<T> &c,
const point<T> &d, const T &epsilon = T{}) {
int ab_c = orientation(a, b, c, epsilon);
int ab_d = orientation(a, b, d, epsilon);
int cd_a = orientation(c, d, a, epsilon);
int cd_b = orientation(c, d, b, epsilon);
if (ab_c == 0 && on_segment(c, a, b, epsilon)) {
return true;
}
if (ab_d == 0 && on_segment(d, a, b, epsilon)) {
return true;
}
if (cd_a == 0 && on_segment(a, c, d, epsilon)) {
return true;
}
if (cd_b == 0 && on_segment(b, c, d, epsilon)) {
return true;
}
return ab_c * ab_d < 0 && cd_a * cd_b < 0;
}
/// @brief Intersect the infinite lines through (a, b) and (c, d), returning
/// nullopt when they are parallel or coincident.
template <class T>
std::optional<point<long double>>
line_intersection(const point<T> &a, const point<T> &b, const point<T> &c,
const point<T> &d, long double epsilon = 0) {
point<long double> first{static_cast<long double>(a.x),
static_cast<long double>(a.y)};
point<long double> second{static_cast<long double>(b.x),
static_cast<long double>(b.y)};
point<long double> third{static_cast<long double>(c.x),
static_cast<long double>(c.y)};
point<long double> fourth{static_cast<long double>(d.x),
static_cast<long double>(d.y)};
point<long double> direction_a = second - first;
point<long double> direction_b = fourth - third;
long double denominator = cross(direction_a, direction_b);
if (std::abs(denominator) <= epsilon) {
return std::nullopt;
}
long double ratio = cross(third - first, direction_b) / denominator;
return first + direction_a * ratio;
}
/// @brief Return the convex hull in counter-clockwise order without repetition.
template <class T>
std::vector<point<T>> convex_hull(std::vector<point<T>> points,
bool keep_collinear = false) {
std::sort(points.begin(), points.end());
points.erase(std::unique(points.begin(), points.end()), points.end());
if (points.size() <= 1) {
return points;
}
bool all_collinear = true;
for (int i = 2; i < int(points.size()); i++) {
all_collinear &= orientation(points[0], points[1], points[i]) == 0;
}
if (keep_collinear && all_collinear) {
return points;
}
std::vector<point<T>> lower, upper;
for (const point<T> &p : points) {
while (lower.size() >= 2) {
int turn = orientation(lower[lower.size() - 2], lower.back(), p);
if (turn > 0 || (keep_collinear && turn == 0)) {
break;
}
lower.pop_back();
}
lower.push_back(p);
}
for (auto it = points.rbegin(); it != points.rend(); ++it) {
while (upper.size() >= 2) {
int turn = orientation(upper[upper.size() - 2], upper.back(), *it);
if (turn > 0 || (keep_collinear && turn == 0)) {
break;
}
upper.pop_back();
}
upper.push_back(*it);
}
lower.pop_back();
upper.pop_back();
lower.insert(lower.end(), upper.begin(), upper.end());
return lower;
}
/// @brief Return twice the signed area of a polygon.
template <class T> T polygon_area2(const std::vector<point<T>> &polygon) {
T result{};
for (int i = 0; i < int(polygon.size()); i++) {
result += cross(polygon[i], polygon[(i + 1) % polygon.size()]);
}
return result;
}
/// @brief Classify a point relative to a polygon: -1 outside, 0 boundary, 1
/// inside.
template <class T>
int point_in_polygon(const point<T> &p, const std::vector<point<T>> &polygon) {
bool inside = false;
for (int i = 0; i < int(polygon.size()); i++) {
point<T> a = polygon[i];
point<T> b = polygon[(i + 1) % polygon.size()];
if (on_segment(p, a, b)) {
return 0;
}
if (a.y <= p.y && p.y < b.y && orientation(a, b, p) > 0) {
inside = !inside;
}
if (b.y <= p.y && p.y < a.y && orientation(a, b, p) < 0) {
inside = !inside;
}
}
return inside ? 1 : -1;
}
/// @brief Return the Euclidean distance between two points.
template <class 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));
}
} // namespace noya
namespace noya {
template <class T> struct minimum_enclosing_circle_result {
circle<long double> value;
std::array<point<T>, 3> support{};
int support_size = 0;
};
namespace minimum_enclosing_circle_internal {
inline circle<long double> from_diameter(const point<long double> &first,
const point<long double> &second) {
point<long double> center = (first + second) / 2.0L;
return {center, distance(center, first)};
}
inline circle<long double> from_three(const point<long double> &first,
const point<long double> &second,
const point<long double> &third) {
long double denominator = 2.0L * cross(first, second, third);
if (std::abs(denominator) <= 1e-24L) {
circle<long double> result = from_diameter(first, second);
circle<long double> candidate = from_diameter(first, third);
if (candidate.radius > result.radius) {
result = candidate;
}
candidate = from_diameter(second, third);
if (candidate.radius > result.radius) {
result = candidate;
}
return result;
}
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 {center, distance(center, first)};
}
inline bool contains(const circle<long double> &value,
const point<long double> &candidate) {
long double squared_radius = value.radius * value.radius;
long double tolerance = 1e-12L * std::max(1.0L, squared_radius);
return norm2(candidate - value.center) <= squared_radius + tolerance;
}
} // namespace minimum_enclosing_circle_internal
/// @brief Return the minimum circle and one boundary support set. Randomly
/// ordering the points makes the expected number of constraint rebuilds
/// linear: an outside point must belong to the new boundary, reducing the
/// remaining problem successively to circles through one, two, then three
/// fixed points.
template <class T>
minimum_enclosing_circle_result<T>
minimum_enclosing_circle_with_support(const std::vector<point<T>> &input,
std::uint64_t seed = 712367821ULL) {
assert(!input.empty());
struct indexed_point {
point<long double> value;
point<T> original;
};
std::vector<indexed_point> points;
points.reserve(input.size());
for (const point<T> &value : input) {
points.push_back({{static_cast<long double>(value.x),
static_cast<long double>(value.y)},
value});
}
std::mt19937_64 random(seed);
std::shuffle(points.begin(), points.end(), random);
using namespace minimum_enclosing_circle_internal;
minimum_enclosing_circle_result<T> answer;
answer.value = {points[0].value, 0};
answer.support[0] = points[0].original;
answer.support_size = 1;
auto answer_contains = [&](const indexed_point &candidate) {
if constexpr (!std::is_integral_v<T>) {
return contains(answer.value, candidate.value);
} else {
using wide = __int128;
const auto &a = answer.support[0];
const auto &p = candidate.original;
if (answer.support_size == 1) {
return p == a;
}
const auto &b = answer.support[1];
if (answer.support_size == 2) {
wide dx = wide(2) * p.x - a.x - b.x;
wide dy = wide(2) * p.y - a.y - b.y;
wide diameter2 = wide(a.x - b.x) * (a.x - b.x) +
wide(a.y - b.y) * (a.y - b.y);
return dx * dx + dy * dy <= diameter2;
}
const auto &c = answer.support[2];
wide bx = wide(b.x) - a.x;
wide by = wide(b.y) - a.y;
wide cx = wide(c.x) - a.x;
wide cy = wide(c.y) - a.y;
wide px = wide(p.x) - a.x;
wide py = wide(p.y) - a.y;
wide determinant =
(px * px + py * py) * (bx * cy - by * cx) -
(bx * bx + by * by) * (px * cy - py * cx) +
(cx * cx + cy * cy) * (px * by - py * bx);
return determinant * (bx * cy - by * cx) <= 0;
}
};
for (int i = 1; i < int(points.size()); i++) {
if (answer_contains(points[i])) {
continue;
}
answer.value = {points[i].value, 0};
answer.support[0] = points[i].original;
answer.support_size = 1;
for (int j = 0; j < i; j++) {
if (answer_contains(points[j])) {
continue;
}
answer.value = from_diameter(points[i].value, points[j].value);
answer.support[0] = points[i].original;
answer.support[1] = points[j].original;
answer.support_size = 2;
for (int k = 0; k < j; k++) {
if (!answer_contains(points[k])) {
long double twice_area =
cross(points[i].value, points[j].value, points[k].value);
if (std::abs(twice_area) <= 1e-24L) {
int first = i;
int second = j;
long double best = norm2(points[i].value - points[j].value);
for (auto [left, right] :
{std::pair{i, k}, std::pair{j, k}}) {
long double candidate =
norm2(points[left].value - points[right].value);
if (candidate > best) {
best = candidate;
first = left;
second = right;
}
}
answer.value =
from_diameter(points[first].value, points[second].value);
answer.support[0] = points[first].original;
answer.support[1] = points[second].original;
answer.support_size = 2;
} else {
answer.value =
from_three(points[i].value, points[j].value, points[k].value);
answer.support[0] = points[i].original;
answer.support[1] = points[j].original;
answer.support[2] = points[k].original;
answer.support_size = 3;
}
}
}
}
}
return answer;
}
/// @brief Return the minimum circle containing all points in expected O(n).
template <class T>
circle<long double>
minimum_enclosing_circle(const std::vector<point<T>> &input,
std::uint64_t seed = 712367821ULL) {
return minimum_enclosing_circle_with_support(input, seed).value;
}
/// @brief Mark exactly which integral points lie on the minimum circle. Once
/// the randomized construction supplies two or three support points, equality
/// with their circle is tested after clearing all denominators: a doubled
/// midpoint equation for a diameter, or the integer circumcircle determinant
/// for three non-collinear points.
template <class T>
std::vector<bool>
minimum_enclosing_circle_boundary(const std::vector<point<T>> &input,
std::uint64_t seed = 712367821ULL) {
static_assert(std::is_integral_v<T>);
auto answer = minimum_enclosing_circle_with_support(input, seed);
std::vector<bool> result(input.size());
using wide = __int128;
const auto &a = answer.support[0];
if (answer.support_size == 1) {
for (int i = 0; i < int(input.size()); i++) {
result[i] = input[i] == a;
}
} else if (answer.support_size == 2) {
const auto &b = answer.support[1];
wide radius4 = wide(a.x - b.x) * (a.x - b.x) +
wide(a.y - b.y) * (a.y - b.y);
for (int i = 0; i < int(input.size()); i++) {
wide dx = wide(2) * input[i].x - a.x - b.x;
wide dy = wide(2) * input[i].y - a.y - b.y;
result[i] = dx * dx + dy * dy == radius4;
}
} else {
const auto &b = answer.support[1];
const auto &c = answer.support[2];
wide bx = wide(b.x) - a.x;
wide by = wide(b.y) - a.y;
wide cx = wide(c.x) - a.x;
wide cy = wide(c.y) - a.y;
wide bc = bx * cy - by * cx;
wide bn = bx * bx + by * by;
wide cn = cx * cx + cy * cy;
for (int i = 0; i < int(input.size()); i++) {
wide px = wide(input[i].x) - a.x;
wide py = wide(input[i].y) - a.y;
wide pn = px * px + py * py;
result[i] = pn * bc - bn * (px * cy - py * cx) +
cn * (px * by - py * bx) ==
0;
}
}
return result;
}
} // namespace noya