Linear Objects
(geometry/linear.hpp)
- View this file on GitHub
- Last update: 2026-10-05 22:23:07+09:00
- Include:
#include "geometry/linear.hpp"
Overview
Rational coordinates are supported, including Rational<BigInt>. Predicates
that use exact arithmetic for integral inputs also use exact rational arithmetic
and ignore eps. Constructed coordinates and lengths with long double return
types remain approximations. Complexity bounds count scalar operations; rational
inputs add the underlying arithmetic and gcd costs.
This header provides Line<T>, Segment<T>, and Ray<T> together with
projection, reflection, distances, intersection predicates, segment centroids,
classified intersections, and closest-point witnesses.
Integral predicates use signed 128-bit arithmetic. Constructed coordinates and
distances are returned using long double.
For floating-point coordinates, eps is relative to the directions or segment
lengths involved. Parallelism, orthogonality, and collinearity use scale-aware
determinant or dot-product tests. Segment endpoint checks allow a normalized
parameter error of eps, so uniformly scaling the complete configuration does
not change the result. Ray direction comparisons use the same scale-aware
rules. A degenerate segment has no length scale; its point test therefore uses
eps as an absolute coordinate tolerance.
Types
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;
};
A line requires distinct endpoints. A segment may be degenerate. A ray starts
at origin, contains through, and continues indefinitely in the direction
through - origin; these two points must be distinct.
Intersections among lines, segments, and rays use one common result type:
enum class LinearIntersectionKind {
Empty,
Point,
Segment,
Ray,
Line,
};
struct LinearIntersection {
LinearIntersectionKind kind;
Point<long double> first;
Point<long double> second;
};
The kind describes the complete intersection set. For Point, first and
second are equal. For Segment, they are its two endpoints. For Ray,
first is the origin and second is a distinct point in its forward
direction. For Line, the two points are distinct representatives of the
line. An Empty result has no meaningful point values. Segment endpoints do
not have a canonical order.
linear.hpp provides these overloads:
template <Coordinate T>
LinearIntersection linear_intersection(
const Line<T>& first,
const Line<T>& second,
long double eps = 1e-12L
);
template <Coordinate T>
LinearIntersection linear_intersection(
const Line<T>& line,
const Segment<T>& segment,
long double eps = 1e-12L
);
template <Coordinate T>
LinearIntersection linear_intersection(
const Segment<T>& segment,
const Line<T>& line,
long double eps = 1e-12L
);
template <Coordinate T>
LinearIntersection linear_intersection(
const Segment<T>& first,
const Segment<T>& second,
long double eps = 1e-12L
);
template <Coordinate T>
LinearIntersection linear_intersection(
const Ray<T>& ray,
const Line<T>& line,
long double eps = 1e-12L
);
template <Coordinate T>
LinearIntersection linear_intersection(
const Line<T>& line,
const Ray<T>& ray,
long double eps = 1e-12L
);
template <Coordinate T>
LinearIntersection linear_intersection(
const Ray<T>& ray,
const Segment<T>& segment,
long double eps = 1e-12L
);
template <Coordinate T>
LinearIntersection linear_intersection(
const Segment<T>& segment,
const Ray<T>& ray,
long double eps = 1e-12L
);
template <Coordinate T>
LinearIntersection linear_intersection(
const Ray<T>& first,
const Ray<T>& second,
long double eps = 1e-12L
);
Every overload runs in $O(1)$.
All overloads share one parametric kernel. Internally an object is written as
$P(t)=O+tD$ with one of three parameter domains: every real value for a line,
$[0,1]$ for a segment, or $[0,+\infty)$ for a ray. Nonparallel objects test
the unique crossing parameter against both domains. Collinear objects
intersect their one-dimensional parameter intervals and translate the result
back into Empty, Point, Segment, Ray, or Line. Degenerate segments
are handled as singleton points before a supporting direction is required.
Closest points
struct ClosestPoints {
Point<long double> first;
Point<long double> second;
};
template <Coordinate T>
ClosestPoints closest_points(
const Point<T>& first,
const Point<T>& second
);
template <Coordinate T>
ClosestPoints closest_points(const Line<T>& line, const Point<T>& point);
template <Coordinate T>
ClosestPoints closest_points(const Point<T>& point, const Line<T>& line);
template <Coordinate T>
ClosestPoints closest_points(
const Segment<T>& segment,
const Point<T>& point
);
template <Coordinate T>
ClosestPoints closest_points(
const Point<T>& point,
const Segment<T>& segment
);
template <Coordinate T>
ClosestPoints closest_points(const Ray<T>& ray, const Point<T>& point);
template <Coordinate T>
ClosestPoints closest_points(const Point<T>& point, const Ray<T>& ray);
template <Coordinate T>
ClosestPoints closest_points(
const Line<T>& first,
const Line<T>& second,
long double eps = 1e-12L
);
template <Coordinate T>
ClosestPoints closest_points(
const Line<T>& line,
const Segment<T>& segment,
long double eps = 1e-12L
);
template <Coordinate T>
ClosestPoints closest_points(
const Segment<T>& segment,
const Line<T>& line,
long double eps = 1e-12L
);
template <Coordinate T>
ClosestPoints closest_points(
const Segment<T>& first,
const Segment<T>& second,
long double eps = 1e-12L
);
template <Coordinate T>
ClosestPoints closest_points(
const Line<T>& line,
const Ray<T>& ray,
long double eps = 1e-12L
);
template <Coordinate T>
ClosestPoints closest_points(
const Ray<T>& ray,
const Line<T>& line,
long double eps = 1e-12L
);
template <Coordinate T>
ClosestPoints closest_points(
const Ray<T>& ray,
const Segment<T>& segment,
long double eps = 1e-12L
);
template <Coordinate T>
ClosestPoints closest_points(
const Segment<T>& segment,
const Ray<T>& ray,
long double eps = 1e-12L
);
template <Coordinate T>
ClosestPoints closest_points(
const Ray<T>& first,
const Ray<T>& second,
long double eps = 1e-12L
);
first belongs to the first argument and second belongs to the second. Their
distance is minimal. Intersecting objects return the same point twice.
Reversing the arguments swaps the witnesses; when several minimizing pairs
exist, a deterministic coordinate-based tie-break is used. Every overload
runs in $O(1)$. The linear-object distance overloads return the distance
between these witnesses.
Line and segment operations
| Function | Description | Complexity |
|---|---|---|
centroid(segment) |
Returns the midpoint, which is the centroid under uniform length. | $O(1)$ |
on_line(line, point, eps) |
Tests whether a point lies on an infinite line. | $O(1)$ |
parallel(first, second, eps) |
Tests whether two lines are parallel. | $O(1)$ |
orthogonal(first, second, eps) |
Tests whether two lines are perpendicular. | $O(1)$ |
projection(line, point) |
Returns the orthogonal projection. | $O(1)$ |
reflection(line, point) |
Reflects a point across a line. | $O(1)$ |
intersects(first, second, eps) |
Tests whether two infinite lines share a point, including coincidence. | $O(1)$ |
distance(line, point) |
Point-to-line distance. Both argument orders are supported. | $O(1)$ |
distance(first, second) |
Distance between two infinite lines. | $O(1)$ |
on_segment(segment, point, eps) |
Tests whether a point lies on a closed segment. | $O(1)$ |
projection(segment, point) |
Returns the closest point on the closed segment, including for a degenerate segment. | $O(1)$ |
intersects(first, second, eps) |
Tests whether two closed segments intersect, including overlap. | $O(1)$ |
intersects(line, segment, eps) |
Tests whether an infinite line and a closed segment intersect. Both argument orders are supported. | $O(1)$ |
distance(segment, point) |
Point-to-segment distance. Both argument orders are supported. | $O(1)$ |
distance(first, second) |
Segment-to-segment distance. | $O(1)$ |
distance(line, segment) |
Line-to-segment distance. Both argument orders are supported. | $O(1)$ |
linear_intersection(first, second, eps) |
Classifies the complete intersection as empty, a point, a segment, a ray, or a line. Every line/segment/ray pairing is supported. | $O(1)$ |
For a degenerate segment, linear_intersection treats it as a point. Its
centroid is also that endpoint. Infinite lines have no centroid overload.
Ray operations
| Function | Description | Complexity |
|---|---|---|
on_ray(ray, point, eps) |
Tests whether point belongs to the closed ray. The origin is included. |
$O(1)$ |
projection(ray, point) |
Returns the closest point on the ray, clamped to its origin. | $O(1)$ |
distance(ray, point) |
Returns the point-to-ray distance. Both argument orders are supported. | $O(1)$ |
reflection(line, ray) |
Reflects the complete ray across an infinite line. | $O(1)$ |
reflected_ray(incoming, hit, mirror, eps) |
Creates the outgoing ray after incoming hits mirror at hit. |
$O(1)$ |
intersects(ray, line, eps) |
Tests whether a ray and an infinite line share a point. Both argument orders are supported. | $O(1)$ |
intersects(ray, segment, eps) |
Tests whether a ray and a closed segment share a point. Both argument orders are supported. | $O(1)$ |
intersects(first, second, eps) |
Tests whether two closed rays share a point. | $O(1)$ |
distance(ray, line) |
Minimum ray-to-line distance. Both argument orders are supported. | $O(1)$ |
distance(ray, segment) |
Minimum ray-to-segment distance. Both argument orders are supported. | $O(1)$ |
distance(first, second) |
Minimum distance between two rays. | $O(1)$ |
reflected_ray requires hit to lie on mirror. Its result starts at hit;
the angle of incidence equals the angle of reflection. Collinear overlap counts
as an intersection. A line and ray can intersect in an empty set, a point, or
a ray. A segment and ray can intersect in an empty set, a point, or a segment.
Two rays can intersect in an empty set, a point, a segment, or a ray.
Closest-points example
Segment<long long> segment;
segment.a = Point<long long>(0, 0);
segment.b = Point<long long>(2, 0);
Ray<long long> ray;
ray.origin = Point<long long>(5, 3);
ray.through = Point<long long>(6, 3);
ClosestPoints result = closest_points(segment, ray);
// result.first == (2, 0), result.second == (5, 3)
Example
#include "geometry/linear.hpp"
#include <iostream>
int main() {
using namespace m1une::geometry;
Segment<long long> first;
first.a = Point<long long>(0, 0);
first.b = Point<long long>(2, 2);
Segment<long long> second;
second.a = Point<long long>(0, 2);
second.b = Point<long long>(2, 0);
LinearIntersection result = linear_intersection(first, second);
if (result.kind == LinearIntersectionKind::Point) {
std::cout << result.first.x << " " << result.first.y << "\n";
}
}
Depends on
Required by
- Geometry Bundle
(geometry/all.hpp)
- Circles
(geometry/circle.hpp)
- Circle Coverage Areas
(geometry/circle_coverage_areas.hpp)
- Area of Union of Circles
(geometry/circle_union_area.hpp)
- Convex Decomposition
(geometry/convex_decomposition.hpp)
- Convex Polygons
(geometry/convex_polygon.hpp)
- Half-Plane Intersection
(geometry/half_plane_intersection.hpp)
- Manhattan Segment Intersections
(geometry/manhattan_segment_intersections.hpp)
- Minimum Enclosing Circle
(geometry/minimum_enclosing_circle.hpp)
- Perpendicular Bisector
(geometry/perpendicular_bisector.hpp)
- Polygons
(geometry/polygon.hpp)
- Steiner Convex Decomposition
(geometry/steiner_convex_decomposition.hpp)
Verified with
- verify/geometry/centroid.test.cpp
- verify/geometry/circle_boundary_intersection.test.cpp
- verify/geometry/circle_circle_intersection.test.cpp
- verify/geometry/circle_circle_intersection_area.test.cpp
- verify/geometry/circle_common_tangents.test.cpp
- verify/geometry/circle_coverage_areas.test.cpp
- verify/geometry/circle_filled.test.cpp
- verify/geometry/circle_line_intersection.test.cpp
- verify/geometry/circle_operations.test.cpp
- verify/geometry/circle_polygon_intersection_area.test.cpp
- verify/geometry/circle_ray.test.cpp
- verify/geometry/circle_relation.test.cpp
- verify/geometry/circle_tangent_points.test.cpp
- verify/geometry/circle_union_area.test.cpp
- verify/geometry/circumcircle.test.cpp
- verify/geometry/closest_points.test.cpp
- verify/geometry/convex_decomposition.test.cpp
- verify/geometry/convex_diameter.test.cpp
- verify/geometry/convex_polygon.test.cpp
- verify/geometry/floating_predicates.test.cpp
- verify/geometry/geometry_algorithms.test.cpp
- verify/geometry/half_plane_intersection.test.cpp
- verify/geometry/half_plane_intersection_random.test.cpp
- verify/geometry/incircle.test.cpp
- verify/geometry/is_convex_polygon.test.cpp
- verify/geometry/linear_intersection.test.cpp
- verify/geometry/manhattan_segment_intersections.test.cpp
- verify/geometry/minimum_enclosing_circle.test.cpp
- verify/geometry/perpendicular_bisector.test.cpp
- verify/geometry/point_in_polygon.test.cpp
- verify/geometry/polygon_area.test.cpp
- verify/geometry/polygon_clipping.test.cpp
- verify/geometry/polygon_filled.test.cpp
- verify/geometry/polygon_operations.test.cpp
- verify/geometry/projection.test.cpp
- verify/geometry/rational.test.cpp
- verify/geometry/ray.test.cpp
- verify/geometry/segment_intersection.test.cpp
- verify/geometry/segment_intersection_point.test.cpp
- verify/geometry/steiner_convex_decomposition.test.cpp
- verify/geometry/voronoi_diagram.test.cpp
Code
#ifndef M1UNE_GEOMETRY_LINEAR_HPP
#define M1UNE_GEOMETRY_LINEAR_HPP 1
#include <algorithm>
#include <cassert>
#include <cmath>
#include "point.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
#endif // M1UNE_GEOMETRY_LINEAR_HPP#line 1 "geometry/linear.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#line 1 "geometry/point.hpp"
#line 5 "geometry/point.hpp"
#include <concepts>
#line 7 "geometry/point.hpp"
#include <type_traits>
#line 1 "geometry/detail/floating_predicate.hpp"
namespace m1une {
namespace geometry {
namespace predicate_detail {
template <typename T>
constexpr T absolute(T value) {
return value < T(0) ? -value : value;
}
template <typename T>
constexpr T max_value(T first, T second) {
return first < second ? second : first;
}
template <typename T>
constexpr T vector_scale(T x, T y) {
return max_value(absolute(x), absolute(y));
}
template <bool Exact, typename T>
constexpr int scaled_sign(T value, T scale, long double eps) {
if constexpr (Exact) {
return (value > T(0)) - (value < T(0));
} else {
const T tolerance = T(eps) * scale;
return (value > tolerance) - (value < -tolerance);
}
}
template <bool Exact, typename T>
constexpr T determinant_scale(T ax, T ay, T bx, T by) {
if constexpr (Exact) {
return T(0);
} else {
return vector_scale(ax, ay) * vector_scale(bx, by);
}
}
template <bool Exact, typename T>
constexpr int determinant_sign(
T ax,
T ay,
T bx,
T by,
long double eps
) {
const T determinant = ax * by - ay * bx;
return scaled_sign<Exact>(
determinant,
determinant_scale<Exact>(ax, ay, bx, by),
eps
);
}
template <bool Exact, typename T>
constexpr int orientation_sign(
T direction_x,
T direction_y,
T offset_x,
T offset_y,
long double eps
) {
const T determinant =
direction_x * offset_y - direction_y * offset_x;
T scale = T(0);
if constexpr (!Exact) {
const T direction_scale =
vector_scale(direction_x, direction_y);
scale = direction_scale * max_value(
direction_scale,
vector_scale(offset_x, offset_y)
);
}
return scaled_sign<Exact>(determinant, scale, eps);
}
template <bool Exact, typename T>
constexpr int dot_sign(
T ax,
T ay,
T bx,
T by,
long double eps
) {
const T value = ax * bx + ay * by;
T scale = T(0);
if constexpr (!Exact) {
scale = vector_scale(ax, ay) * vector_scale(bx, by);
}
return scaled_sign<Exact>(value, scale, eps);
}
} // namespace predicate_detail
} // namespace geometry
} // namespace m1une
#line 10 "geometry/point.hpp"
namespace m1une {
namespace geometry {
template <typename T>
concept Coordinate = !std::same_as<std::remove_cv_t<T>, bool> &&
(std::is_arithmetic_v<T> ||
(std::copyable<T> && std::totally_ordered<T> && requires(T a, T b) {
T(0);
T(1);
static_cast<long double>(a);
{ +a } -> std::same_as<T>;
{ -a } -> std::same_as<T>;
{ a + b } -> std::same_as<T>;
{ a - b } -> std::same_as<T>;
{ a * b } -> std::same_as<T>;
{ a / b } -> std::same_as<T>;
{ a += b } -> std::same_as<T&>;
{ a -= b } -> std::same_as<T&>;
}));
// Custom coordinate types keep their own exact arithmetic.
template <typename T>
concept ExactCoordinate = Coordinate<T> && !std::floating_point<T>;
template <Coordinate T>
using wide_type = std::conditional_t<std::integral<T>, __int128_t,
std::conditional_t<std::floating_point<T>, long double, T>>;
template <Coordinate T>
struct Point {
T x;
T y;
constexpr Point() : x(0), y(0) {}
constexpr Point(T x_value, T y_value) : x(x_value), y(y_value) {}
template <Coordinate U>
explicit constexpr Point(const Point<U>& other)
: x(static_cast<T>(other.x)), y(static_cast<T>(other.y)) {}
constexpr Point& operator+=(const Point& other) {
x += other.x;
y += other.y;
return *this;
}
constexpr Point& operator-=(const Point& other) {
x -= other.x;
y -= other.y;
return *this;
}
constexpr Point operator+() const {
return *this;
}
constexpr Point operator-() const {
return Point(-x, -y);
}
friend constexpr Point operator+(Point left, const Point& right) {
return left += right;
}
friend constexpr Point operator-(Point left, const Point& right) {
return left -= right;
}
friend constexpr bool operator==(const Point&, const Point&) = default;
friend constexpr bool operator<(const Point& left, const Point& right) {
if (left.x != right.x) return left.x < right.x;
return left.y < right.y;
}
};
template <Coordinate T>
constexpr Point<long double> centroid(const Point<T>& point) {
return Point<long double>(point);
}
template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(const Point<T>& point, Scalar scalar) {
using Result = std::common_type_t<T, Scalar>;
return Point<Result>(
Result(point.x) * Result(scalar),
Result(point.y) * Result(scalar)
);
}
template <typename Scalar, Coordinate T>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(Scalar scalar, const Point<T>& point) {
return point * scalar;
}
template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator/(const Point<T>& point, Scalar scalar) {
using Result = std::common_type_t<T, Scalar>;
return Point<Result>(
Result(point.x) / Result(scalar),
Result(point.y) / Result(scalar)
);
}
template <Coordinate T>
constexpr wide_type<T> dot(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
return W(a.x) * W(b.x) + W(a.y) * W(b.y);
}
template <Coordinate T>
constexpr wide_type<T> cross(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
return W(a.x) * W(b.y) - W(a.y) * W(b.x);
}
template <Coordinate T>
constexpr wide_type<T> cross(
const Point<T>& origin,
const Point<T>& a,
const Point<T>& b
) {
using W = wide_type<T>;
W ax = W(a.x) - W(origin.x);
W ay = W(a.y) - W(origin.y);
W bx = W(b.x) - W(origin.x);
W by = W(b.y) - W(origin.y);
return ax * by - ay * bx;
}
template <Coordinate T>
constexpr wide_type<T> norm2(const Point<T>& point) {
return dot(point, point);
}
template <Coordinate T>
constexpr wide_type<T> distance2(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
W dx = W(a.x) - W(b.x);
W dy = W(a.y) - W(b.y);
return dx * dx + dy * dy;
}
template <Coordinate T>
long double norm(const Point<T>& point) {
return std::hypot(
static_cast<long double>(point.x),
static_cast<long double>(point.y)
);
}
template <Coordinate T>
long double distance(const Point<T>& a, const Point<T>& b) {
return std::hypot(
static_cast<long double>(a.x) - static_cast<long double>(b.x),
static_cast<long double>(a.y) - static_cast<long double>(b.y)
);
}
template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
(std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> internal_division_point(
const Point<T>& a,
const Point<T>& b,
M m,
N n
) {
long double first_ratio = static_cast<long double>(m);
long double second_ratio = static_cast<long double>(n);
long double denominator = first_ratio + second_ratio;
assert(denominator != 0);
Point<long double> first(a);
Point<long double> direction = Point<long double>(b) - first;
return first + direction * (first_ratio / denominator);
}
template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
(std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> external_division_point(
const Point<T>& a,
const Point<T>& b,
M m,
N n
) {
long double first_ratio = static_cast<long double>(m);
long double second_ratio = static_cast<long double>(n);
long double denominator = first_ratio - second_ratio;
assert(denominator != 0);
Point<long double> first(a);
Point<long double> direction = Point<long double>(b) - first;
return first + direction * (first_ratio / denominator);
}
template <Coordinate T>
constexpr int sign(wide_type<T> value, long double eps = 1e-12L) {
return predicate_detail::scaled_sign<ExactCoordinate<T>>(
value,
wide_type<T>(1),
eps
);
}
template <Coordinate T>
constexpr int orientation(
const Point<T>& a,
const Point<T>& b,
const Point<T>& c,
long double eps = 1e-12L
) {
using W = wide_type<T>;
const W first_x = W(b.x) - W(a.x);
const W first_y = W(b.y) - W(a.y);
const W second_x = W(c.x) - W(a.x);
const W second_y = W(c.y) - W(a.y);
return predicate_detail::orientation_sign<ExactCoordinate<T>>(
first_x,
first_y,
second_x,
second_y,
eps
);
}
template <Coordinate T>
constexpr bool collinear(
const Point<T>& a,
const Point<T>& b,
const Point<T>& c,
long double eps = 1e-12L
) {
return orientation(a, b, c, eps) == 0;
}
template <Coordinate T>
Point<long double> rotate(const Point<T>& point, long double angle) {
long double cosine = std::cos(angle);
long double sine = std::sin(angle);
return Point<long double>(
static_cast<long double>(point.x) * cosine -
static_cast<long double>(point.y) * sine,
static_cast<long double>(point.x) * sine +
static_cast<long double>(point.y) * cosine
);
}
template <Coordinate T>
Point<long double> normalized(const Point<T>& point) {
long double length = norm(point);
assert(length != 0);
return Point<long double>(
static_cast<long double>(point.x) / length,
static_cast<long double>(point.y) / length
);
}
} // namespace geometry
} // namespace m1une
#line 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