circle_geometry.hpp¶
Circle intersections, line intersections, tangent points, common tangents, and circumcircles.
求圆与圆/直线交点、切点、公切线和三点外接圆；适合圆相关构造题。
Implementation¶
#ifndef NOYA_CIRCLE_GEOMETRY_HPP
#define NOYA_CIRCLE_GEOMETRY_HPP 1
/// @complexity Time: O(1) per primitive.
/// Space: O(1), excluding returned points.
#include "noya/geometry_base.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <optional>
#include <utility>
#include <vector>
namespace noya {
namespace circle_geometry_internal {
template <class T> point<long double> cast_point(const point<T> &value) {
return {static_cast<long double>(value.x),
static_cast<long double>(value.y)};
}
inline point<long double> perpendicular(const point<long double> &value) {
return {-value.y, value.x};
}
inline bool close_points(const point<long double> &first,
const point<long double> &second,
long double epsilon) {
return norm2(first - second) <= epsilon * epsilon;
}
} // namespace circle_geometry_internal
/// @brief Intersect two circles: nullopt means coincident circles, while a
/// vector of size 0, 1, or 2 represents no, tangent, or two intersections.
template <class FirstReal, class SecondReal>
std::optional<std::vector<point<long double>>>
circle_intersections(const circle<FirstReal> &first,
const circle<SecondReal> &second,
long double epsilon = 1e-12L) {
assert(first.radius >= FirstReal{} && second.radius >= SecondReal{});
using circle_geometry_internal::cast_point;
using circle_geometry_internal::perpendicular;
point<long double> first_center = cast_point(first.center);
point<long double> second_center = cast_point(second.center);
point<long double> delta = second_center - first_center;
long double distance_squared = norm2(delta);
long double first_radius = static_cast<long double>(first.radius);
long double second_radius = static_cast<long double>(second.radius);
long double scale =
std::max({1.0L, first_radius, second_radius,
std::sqrt(distance_squared)});
long double tolerance = epsilon * scale;
if (distance_squared <= tolerance * tolerance) {
if (std::abs(first_radius - second_radius) <= tolerance) {
return std::nullopt;
}
return std::vector<point<long double>>{};
}
long double center_distance = std::sqrt(distance_squared);
if (center_distance > first_radius + second_radius + tolerance ||
center_distance < std::abs(first_radius - second_radius) - tolerance) {
return std::vector<point<long double>>{};
}
long double along =
(distance_squared + first_radius * first_radius -
second_radius * second_radius) /
(2 * center_distance);
long double height_squared = first_radius * first_radius - along * along;
if (height_squared < 0) {
height_squared = 0;
}
point<long double> direction = delta / center_distance;
point<long double> middle = first_center + direction * along;
if (height_squared <= tolerance * tolerance) {
return std::vector<point<long double>>{middle};
}
point<long double> offset =
perpendicular(direction) * std::sqrt(height_squared);
return std::vector<point<long double>>{middle + offset, middle - offset};
}
/// @brief Intersect an infinite line through first and second with a circle,
/// returning zero, one, or two points.
template <class PointReal, class CircleReal>
std::vector<point<long double>>
line_circle_intersections(const point<PointReal> &first,
const point<PointReal> &second,
const circle<CircleReal> &value,
long double epsilon = 1e-12L) {
assert(value.radius >= CircleReal{});
using circle_geometry_internal::cast_point;
point<long double> a = cast_point(first);
point<long double> b = cast_point(second);
point<long double> center = cast_point(value.center);
point<long double> direction = b - a;
long double length_squared = norm2(direction);
assert(length_squared > 0);
long double projection = dot(center - a, direction) / length_squared;
point<long double> closest = a + direction * projection;
long double radius = static_cast<long double>(value.radius);
long double remaining = radius * radius - norm2(closest - center);
long double tolerance = epsilon * std::max(1.0L, radius);
if (remaining < -tolerance * tolerance) {
return {};
}
if (remaining <= tolerance * tolerance) {
return {closest};
}
point<long double> offset =
direction * std::sqrt(remaining / length_squared);
return {closest - offset, closest + offset};
}
/// @brief Return the tangent points from an external/on-circle point; an
/// interior point has no tangent point.
template <class PointReal, class CircleReal>
std::vector<point<long double>>
point_circle_tangents(const point<PointReal> &external,
const circle<CircleReal> &value,
long double epsilon = 1e-12L) {
assert(value.radius > CircleReal{});
using circle_geometry_internal::cast_point;
using circle_geometry_internal::perpendicular;
point<long double> center = cast_point(value.center);
point<long double> relative = cast_point(external) - center;
long double distance_squared = norm2(relative);
long double radius = static_cast<long double>(value.radius);
long double radius_squared = radius * radius;
long double tolerance = epsilon * std::max(1.0L, radius_squared);
if (distance_squared < radius_squared - tolerance) {
return {};
}
if (std::abs(distance_squared - radius_squared) <= tolerance) {
return {cast_point(external)};
}
long double along = radius_squared / distance_squared;
long double across =
radius * std::sqrt(distance_squared - radius_squared) / distance_squared;
point<long double> middle = center + relative * along;
point<long double> offset = perpendicular(relative) * across;
return {middle + offset, middle - offset};
}
/// @brief Return all common tangents as pairs of contact points; nullopt means
/// equal coincident circles and therefore infinitely many tangents.
template <class FirstReal, class SecondReal>
std::optional<std::vector<
std::pair<point<long double>, point<long double>>>>
common_circle_tangents(const circle<FirstReal> &first,
const circle<SecondReal> &second,
long double epsilon = 1e-12L) {
assert(first.radius >= FirstReal{} && second.radius >= SecondReal{});
using circle_geometry_internal::cast_point;
using circle_geometry_internal::close_points;
using circle_geometry_internal::perpendicular;
point<long double> first_center = cast_point(first.center);
point<long double> second_center = cast_point(second.center);
point<long double> delta = second_center - first_center;
long double distance_squared = norm2(delta);
long double first_radius = static_cast<long double>(first.radius);
long double second_radius = static_cast<long double>(second.radius);
long double scale =
std::max({1.0L, first_radius, second_radius,
std::sqrt(distance_squared)});
long double tolerance = epsilon * scale;
if (distance_squared <= tolerance * tolerance) {
if (std::abs(first_radius - second_radius) <= tolerance) {
return std::nullopt;
}
return std::vector<
std::pair<point<long double>, point<long double>>>{};
}
std::vector<std::pair<point<long double>, point<long double>>> result;
for (int radius_sign : {-1, 1}) {
long double signed_second_radius = radius_sign * second_radius;
long double radius_difference = first_radius - signed_second_radius;
long double height_squared =
distance_squared - radius_difference * radius_difference;
if (height_squared < -tolerance * tolerance) {
continue;
}
height_squared = std::max(0.0L, height_squared);
long double height = std::sqrt(height_squared);
for (int side : {-1, 1}) {
point<long double> normal =
(delta * radius_difference +
perpendicular(delta) * (height * side)) /
distance_squared;
std::pair<point<long double>, point<long double>> tangent = {
first_center + normal * first_radius,
second_center + normal * signed_second_radius};
bool duplicate = false;
for (const auto &existing : result) {
duplicate |= close_points(existing.first, tangent.first, tolerance) &&
close_points(existing.second, tangent.second, tolerance);
}
if (!duplicate) {
result.push_back(tangent);
}
if (height == 0) {
break;
}
}
}
return result;
}
/// @brief Return the circumcircle through three non-collinear points, or
/// nullopt when the points are collinear.
template <class T>
std::optional<circle<long double>>
circumcircle(const point<T> &first, const point<T> &second,
const point<T> &third, long double epsilon = 1e-18L) {
using circle_geometry_internal::cast_point;
point<long double> a = cast_point(first);
point<long double> b = cast_point(second);
point<long double> c = cast_point(third);
long double denominator = 2 * cross(a, b, c);
if (std::abs(denominator) <= epsilon) {
return std::nullopt;
}
long double a_norm = norm2(a);
long double b_norm = norm2(b);
long double c_norm = norm2(c);
point<long double> center{
(a_norm * (b.y - c.y) + b_norm * (c.y - a.y) +
c_norm * (a.y - b.y)) /
denominator,
(a_norm * (c.x - b.x) + b_norm * (a.x - c.x) +
c_norm * (b.x - a.x)) /
denominator};
return circle<long double>{center, distance(center, a)};
}
} // namespace noya
#endif // NOYA_CIRCLE_GEOMETRY_HPP
#include <algorithm>
#include <cassert>
#include <cmath>
#include <optional>
#include <utility>
#include <vector>
/// @complexity Time: O(1) per primitive.
/// Space: O(1), excluding returned 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 {
namespace circle_geometry_internal {
template <class T> point<long double> cast_point(const point<T> &value) {
return {static_cast<long double>(value.x),
static_cast<long double>(value.y)};
}
inline point<long double> perpendicular(const point<long double> &value) {
return {-value.y, value.x};
}
inline bool close_points(const point<long double> &first,
const point<long double> &second,
long double epsilon) {
return norm2(first - second) <= epsilon * epsilon;
}
} // namespace circle_geometry_internal
/// @brief Intersect two circles: nullopt means coincident circles, while a
/// vector of size 0, 1, or 2 represents no, tangent, or two intersections.
template <class FirstReal, class SecondReal>
std::optional<std::vector<point<long double>>>
circle_intersections(const circle<FirstReal> &first,
const circle<SecondReal> &second,
long double epsilon = 1e-12L) {
assert(first.radius >= FirstReal{} && second.radius >= SecondReal{});
using circle_geometry_internal::cast_point;
using circle_geometry_internal::perpendicular;
point<long double> first_center = cast_point(first.center);
point<long double> second_center = cast_point(second.center);
point<long double> delta = second_center - first_center;
long double distance_squared = norm2(delta);
long double first_radius = static_cast<long double>(first.radius);
long double second_radius = static_cast<long double>(second.radius);
long double scale =
std::max({1.0L, first_radius, second_radius,
std::sqrt(distance_squared)});
long double tolerance = epsilon * scale;
if (distance_squared <= tolerance * tolerance) {
if (std::abs(first_radius - second_radius) <= tolerance) {
return std::nullopt;
}
return std::vector<point<long double>>{};
}
long double center_distance = std::sqrt(distance_squared);
if (center_distance > first_radius + second_radius + tolerance ||
center_distance < std::abs(first_radius - second_radius) - tolerance) {
return std::vector<point<long double>>{};
}
long double along =
(distance_squared + first_radius * first_radius -
second_radius * second_radius) /
(2 * center_distance);
long double height_squared = first_radius * first_radius - along * along;
if (height_squared < 0) {
height_squared = 0;
}
point<long double> direction = delta / center_distance;
point<long double> middle = first_center + direction * along;
if (height_squared <= tolerance * tolerance) {
return std::vector<point<long double>>{middle};
}
point<long double> offset =
perpendicular(direction) * std::sqrt(height_squared);
return std::vector<point<long double>>{middle + offset, middle - offset};
}
/// @brief Intersect an infinite line through first and second with a circle,
/// returning zero, one, or two points.
template <class PointReal, class CircleReal>
std::vector<point<long double>>
line_circle_intersections(const point<PointReal> &first,
const point<PointReal> &second,
const circle<CircleReal> &value,
long double epsilon = 1e-12L) {
assert(value.radius >= CircleReal{});
using circle_geometry_internal::cast_point;
point<long double> a = cast_point(first);
point<long double> b = cast_point(second);
point<long double> center = cast_point(value.center);
point<long double> direction = b - a;
long double length_squared = norm2(direction);
assert(length_squared > 0);
long double projection = dot(center - a, direction) / length_squared;
point<long double> closest = a + direction * projection;
long double radius = static_cast<long double>(value.radius);
long double remaining = radius * radius - norm2(closest - center);
long double tolerance = epsilon * std::max(1.0L, radius);
if (remaining < -tolerance * tolerance) {
return {};
}
if (remaining <= tolerance * tolerance) {
return {closest};
}
point<long double> offset =
direction * std::sqrt(remaining / length_squared);
return {closest - offset, closest + offset};
}
/// @brief Return the tangent points from an external/on-circle point; an
/// interior point has no tangent point.
template <class PointReal, class CircleReal>
std::vector<point<long double>>
point_circle_tangents(const point<PointReal> &external,
const circle<CircleReal> &value,
long double epsilon = 1e-12L) {
assert(value.radius > CircleReal{});
using circle_geometry_internal::cast_point;
using circle_geometry_internal::perpendicular;
point<long double> center = cast_point(value.center);
point<long double> relative = cast_point(external) - center;
long double distance_squared = norm2(relative);
long double radius = static_cast<long double>(value.radius);
long double radius_squared = radius * radius;
long double tolerance = epsilon * std::max(1.0L, radius_squared);
if (distance_squared < radius_squared - tolerance) {
return {};
}
if (std::abs(distance_squared - radius_squared) <= tolerance) {
return {cast_point(external)};
}
long double along = radius_squared / distance_squared;
long double across =
radius * std::sqrt(distance_squared - radius_squared) / distance_squared;
point<long double> middle = center + relative * along;
point<long double> offset = perpendicular(relative) * across;
return {middle + offset, middle - offset};
}
/// @brief Return all common tangents as pairs of contact points; nullopt means
/// equal coincident circles and therefore infinitely many tangents.
template <class FirstReal, class SecondReal>
std::optional<std::vector<
std::pair<point<long double>, point<long double>>>>
common_circle_tangents(const circle<FirstReal> &first,
const circle<SecondReal> &second,
long double epsilon = 1e-12L) {
assert(first.radius >= FirstReal{} && second.radius >= SecondReal{});
using circle_geometry_internal::cast_point;
using circle_geometry_internal::close_points;
using circle_geometry_internal::perpendicular;
point<long double> first_center = cast_point(first.center);
point<long double> second_center = cast_point(second.center);
point<long double> delta = second_center - first_center;
long double distance_squared = norm2(delta);
long double first_radius = static_cast<long double>(first.radius);
long double second_radius = static_cast<long double>(second.radius);
long double scale =
std::max({1.0L, first_radius, second_radius,
std::sqrt(distance_squared)});
long double tolerance = epsilon * scale;
if (distance_squared <= tolerance * tolerance) {
if (std::abs(first_radius - second_radius) <= tolerance) {
return std::nullopt;
}
return std::vector<
std::pair<point<long double>, point<long double>>>{};
}
std::vector<std::pair<point<long double>, point<long double>>> result;
for (int radius_sign : {-1, 1}) {
long double signed_second_radius = radius_sign * second_radius;
long double radius_difference = first_radius - signed_second_radius;
long double height_squared =
distance_squared - radius_difference * radius_difference;
if (height_squared < -tolerance * tolerance) {
continue;
}
height_squared = std::max(0.0L, height_squared);
long double height = std::sqrt(height_squared);
for (int side : {-1, 1}) {
point<long double> normal =
(delta * radius_difference +
perpendicular(delta) * (height * side)) /
distance_squared;
std::pair<point<long double>, point<long double>> tangent = {
first_center + normal * first_radius,
second_center + normal * signed_second_radius};
bool duplicate = false;
for (const auto &existing : result) {
duplicate |= close_points(existing.first, tangent.first, tolerance) &&
close_points(existing.second, tangent.second, tolerance);
}
if (!duplicate) {
result.push_back(tangent);
}
if (height == 0) {
break;
}
}
}
return result;
}
/// @brief Return the circumcircle through three non-collinear points, or
/// nullopt when the points are collinear.
template <class T>
std::optional<circle<long double>>
circumcircle(const point<T> &first, const point<T> &second,
const point<T> &third, long double epsilon = 1e-18L) {
using circle_geometry_internal::cast_point;
point<long double> a = cast_point(first);
point<long double> b = cast_point(second);
point<long double> c = cast_point(third);
long double denominator = 2 * cross(a, b, c);
if (std::abs(denominator) <= epsilon) {
return std::nullopt;
}
long double a_norm = norm2(a);
long double b_norm = norm2(b);
long double c_norm = norm2(c);
point<long double> center{
(a_norm * (b.y - c.y) + b_norm * (c.y - a.y) +
c_norm * (a.y - b.y)) /
denominator,
(a_norm * (c.x - b.x) + b_norm * (a.x - c.x) +
c_norm * (b.x - a.x)) /
denominator};
return circle<long double>{center, distance(center, a)};
}
} // namespace noya