幾何の直線・線分

geometry/linear.hpp · GitHub で見る

共通ページ

#include "lib/geometry/linear.hpp"
using namespace geometry;

Segment<std::int64_t> s{IPoint{0, 0}, IPoint{4, 0}};
Segment<std::int64_t> t{IPoint{2, -1}, IPoint{2, 1}};
bool hit = intersects(s, t); // true
auto points = intersection(s, t); // (2,0) の1点
Real d = distance(s, IPoint{5, 1}); // sqrt(2)

Line<T>{a,b} は a・b を通る無限直線、Segment<T>{a,b} は両端を含む線分です。a・b は Point<T> です。引数同士の座標型は揃えて使います。 始点と終点が同じものも点として扱います。整数の交差判定は128 bitの整数で行い、交点・射影・距離は Real または FPoint で返します。

操作	内容・返り値	計算量
Line<T>{a,b}, Segment<T>{a,b}	直線・閉線分の構築。端点は .a,.b	\(O(1)\)
parallel(l,m,eps = epsilon)	方向ベクトルの外積が0か	\(O(1)\)
orthogonal(l,m,eps = epsilon)	方向ベクトルの内積が0か	\(O(1)\)
on_line(l,p,eps = epsilon)	p が直線上か	\(O(1)\)
on_segment(s,p,eps = epsilon)	p が線分上か。端点も含む	\(O(1)\)
intersects(a,b,eps = epsilon)	直線・線分の任意の組で共通点があるか	\(O(1)\)
projection(l,p)	p の直線への射影。退化直線ならその点	\(O(1)\)
reflection(l,p)	直線について p と対称な点。退化直線ならその点を中心に点対称	\(O(1)\)
intersection(l,m,eps = epsilon)	直線同士の唯一の交点を optional<FPoint> で返す。平行・一致直線は nullopt	\(O(1)\)
intersection(s,t,eps = epsilon)	線分同士の共通部分を vector<FPoint> で返す。0点・1点・重なりの両端2点	\(O(1)\)
intersection(l,s,eps = epsilon), intersection(s,l,eps)	直線と線分の共通部分。重なる場合は線分の両端	\(O(1)\)
distance(a,b)	点・直線・線分の任意の組の最短距離。対称な引数順も使用可	\(O(1)\)
distance_squared(a,b)	同じ組の距離二乗。点同士は wide_t<T>、それ以外は Real	\(O(1)\)
perpendicular_bisector(a,b)	垂直二等分線を optional<Line<Real>> で返す。同一点なら nullopt	\(O(1)\)

線分同士の重なりの両端は辞書順です。直線と線分の重なりは線分の a・b の順です。退化した入力でも共通部分が1点なら、その点を返します。 parallel・orthogonal は方向ベクトルに対する判定なので、零ベクトルを含む場合は両方 true です。幾何学的な交差判定には intersects を使います。

整数の場合は EPS を無視します。実数の交差判定では外積・内積に対する絶対 EPS を使い、線分の座標範囲判定には座標単位の EPS を使います。 直線同士・直線と線分・線分同士の距離では、必要なら第3引数に EPS を渡せます。点への距離・射影では EPS で位置を動かしません。 128 bitの中間値と実数精度については点の説明の前提に従います。実数の交点を厳密な有理数として返す関数ではありません。

全ての操作の時間・追加空間は \(O(1)\) です。

依存ライブラリ

• 64 bit キーの安定基数ソート
• 幾何の点・ベクトル・偏角順

検証コード

• test/geometry_circle_line.test.cpp
• test/geometry_circle_polygon.test.cpp
• test/geometry_closest_pair.test.cpp
• test/geometry_convex_cut.test.cpp
• test/geometry_convex_hull.test.cpp
• test/geometry_furthest_pair.test.cpp
• test/geometry_polygon_contains.test.cpp
• test/geometry_projection.test.cpp
• test/geometry_segment_distance.test.cpp

コード

#ifndef GEOMETRY_LINEAR_HPP
#define GEOMETRY_LINEAR_HPP

#include "point.hpp"
#include <algorithm>
#include <array>
#include <vector>

namespace geometry {
template<class T> struct Line { Point<T> a, b; };
template<class T> struct Segment { Point<T> a, b; };

template<class T>
constexpr bool parallel(Line<T> l, Line<T> m, Real eps = epsilon) {
	using W = wide_t<T>;
	return sign((W(l.b.x) - l.a.x) * (W(m.b.y) - m.a.y) - (W(l.b.y) - l.a.y) * (W(m.b.x) - m.a.x), eps) == 0;
}
template<class T>
constexpr bool orthogonal(Line<T> l, Line<T> m, Real eps = epsilon) {
	using W = wide_t<T>;
	return sign((W(l.b.x) - l.a.x) * (W(m.b.x) - m.a.x) + (W(l.b.y) - l.a.y) * (W(m.b.y) - m.a.y), eps) == 0;
}
template<class T>
constexpr bool on_line(Line<T> l, Point<T> p, Real eps = epsilon) {
	return l.a == l.b ? almost_equal(l.a, p, eps) : orientation(l.a, l.b, p, eps) == 0;
}
template<class T>
constexpr bool on_segment(Segment<T> s, Point<T> p, Real eps = epsilon) {
	if (orientation(s.a, s.b, p, eps)) return false;
	if constexpr (integer<T>)
		return std::min(s.a.x, s.b.x) <= p.x && p.x <= std::max(s.a.x, s.b.x)
			&& std::min(s.a.y, s.b.y) <= p.y && p.y <= std::max(s.a.y, s.b.y);
	else return std::min(s.a.x, s.b.x) - eps <= p.x && p.x <= std::max(s.a.x, s.b.x) + eps
		&& std::min(s.a.y, s.b.y) - eps <= p.y && p.y <= std::max(s.a.y, s.b.y) + eps;
}
template<class T>
constexpr bool intersects(Segment<T> s, Segment<T> t, Real eps = epsilon) {
	int a = orientation(s.a, s.b, t.a, eps), b = orientation(s.a, s.b, t.b, eps);
	int c = orientation(t.a, t.b, s.a, eps), d = orientation(t.a, t.b, s.b, eps);
	return (!a && on_segment(s, t.a, eps)) || (!b && on_segment(s, t.b, eps))
		|| (!c && on_segment(t, s.a, eps)) || (!d && on_segment(t, s.b, eps)) || (a * b < 0 && c * d < 0);
}
template<class T>
constexpr bool intersects(Line<T> l, Segment<T> s, Real eps = epsilon) {
	if (l.a == l.b) return on_segment(s, l.a, eps);
	return orientation(l.a, l.b, s.a, eps) * orientation(l.a, l.b, s.b, eps) <= 0;
}
template<class T>
constexpr bool intersects(Segment<T> s, Line<T> l, Real eps = epsilon) { return intersects(l, s, eps); }
template<class T>
constexpr bool intersects(Line<T> l, Line<T> m, Real eps = epsilon) {
	if (l.a == l.b) return on_line(m, l.a, eps);
	if (m.a == m.b) return on_line(l, m.a, eps);
	return !parallel(l, m, eps) || on_line(l, m.a, eps);
}
template<class T>
FPoint projection(Line<T> l, Point<T> p) {
	using W = wide_t<T>;
	W x = W(l.b.x) - l.a.x, y = W(l.b.y) - l.a.y, n = x * x + y * y;
	if (!n) return to_real(l.a);
	Real t = Real(x * (W(p.x) - l.a.x) + y * (W(p.y) - l.a.y)) / Real(n);
	return {Real(l.a.x) + Real(x) * t, Real(l.a.y) + Real(y) * t};
}
template<class T>
FPoint reflection(Line<T> l, Point<T> p) { return projection(l, p) * Real(2) - to_real(p); }
template<class T>
std::optional<FPoint> intersection(Line<T> l, Line<T> m, Real eps = epsilon) {
	if (l.a == l.b) return on_line(m, l.a, eps) ? std::optional(to_real(l.a)) : std::nullopt;
	if (m.a == m.b) return on_line(l, m.a, eps) ? std::optional(to_real(m.a)) : std::nullopt;
	using W = wide_t<T>;
	W x = W(l.b.x) - l.a.x, y = W(l.b.y) - l.a.y, u = W(m.b.x) - m.a.x, v = W(m.b.y) - m.a.y;
	W d = x * v - y * u;
	if (!sign(d, eps)) return std::nullopt;
	Real t = Real((W(m.a.x) - l.a.x) * v - (W(m.a.y) - l.a.y) * u) / Real(d);
	return FPoint{Real(l.a.x) + Real(x) * t, Real(l.a.y) + Real(y) * t};
}
template<class T>
std::vector<FPoint> intersection(Segment<T> s, Segment<T> t, Real eps = epsilon) {
	if (!intersects(s, t, eps)) return {};
	if (auto p = intersection(Line<T>{s.a, s.b}, Line<T>{t.a, t.b}, eps)) return {*p};
	std::array<Point<T>, 4> candidates{s.a, s.b, t.a, t.b};
	std::sort(candidates.begin(), candidates.end());
	std::vector<FPoint> result;
	std::optional<Point<T>> first, last;
	for (auto p : candidates) if (on_segment(s, p, eps) && on_segment(t, p, eps)) {
		if (!first) first = p;
		last = p;
	}
	if (first) {
		result.push_back(to_real(*first));
		if (!almost_equal(*first, *last, eps)) result.push_back(to_real(*last));
	}
	return result;
}
template<class T>
std::vector<FPoint> intersection(Line<T> l, Segment<T> s, Real eps = epsilon) {
	if (!intersects(l, s, eps)) return {};
	if (auto p = intersection(l, Line<T>{s.a, s.b}, eps)) return {*p};
	if (almost_equal(s.a, s.b, eps)) return {to_real(s.a)};
	return {to_real(s.a), to_real(s.b)};
}
template<class T>
std::vector<FPoint> intersection(Segment<T> s, Line<T> l, Real eps = epsilon) { return intersection(l, s, eps); }
template<class T>
Real distance_squared(Line<T> l, Point<T> p) {
	if (l.a == l.b) return Real(distance_squared(l.a, p));
	Real c = Real(cross(l.a, l.b, p));
	return c * c / Real(distance_squared(l.a, l.b));
}
template<class T>
Real distance_squared(Segment<T> s, Point<T> p) {
	using W = wide_t<T>;
	W x = W(s.b.x) - s.a.x, y = W(s.b.y) - s.a.y;
	W n = x * x + y * y, d = (W(p.x) - s.a.x) * x + (W(p.y) - s.a.y) * y;
	if (d <= 0) return Real(distance_squared(s.a, p));
	if (d >= n) return Real(distance_squared(s.b, p));
	Real c = Real(cross(s.a, s.b, p));
	return c * c / Real(n);
}
template<class T> Real distance(Line<T> l, Point<T> p) { return std::sqrt(distance_squared(l, p)); }
template<class T> Real distance(Point<T> p, Line<T> l) { return distance(l, p); }
template<class T> Real distance(Segment<T> s, Point<T> p) { return std::sqrt(distance_squared(s, p)); }
template<class T> Real distance(Point<T> p, Segment<T> s) { return distance(s, p); }
template<class T> Real distance_squared(Point<T> p, Line<T> l) { return distance_squared(l, p); }
template<class T> Real distance_squared(Point<T> p, Segment<T> s) { return distance_squared(s, p); }
template<class T>
Real distance_squared(Line<T> l, Line<T> m, Real eps = epsilon) {
	if (intersects(l, m, eps)) return 0;
	return l.a == l.b ? distance_squared(m, l.a) : distance_squared(l, m.a);
}
template<class T>
Real distance_squared(Line<T> l, Segment<T> s, Real eps = epsilon) {
	if (intersects(l, s, eps)) return 0;
	return l.a == l.b ? distance_squared(s, l.a) : std::min(distance_squared(l, s.a), distance_squared(l, s.b));
}
template<class T>
Real distance_squared(Segment<T> s, Line<T> l, Real eps = epsilon) { return distance_squared(l, s, eps); }
template<class T>
Real distance_squared(Segment<T> s, Segment<T> t, Real eps = epsilon) {
	if (intersects(s, t, eps)) return 0;
	return std::min({distance_squared(s, t.a), distance_squared(s, t.b), distance_squared(t, s.a), distance_squared(t, s.b)});
}
template<class T> Real distance(Line<T> l, Line<T> m, Real eps = epsilon) { return std::sqrt(distance_squared(l, m, eps)); }
template<class T> Real distance(Line<T> l, Segment<T> s, Real eps = epsilon) { return std::sqrt(distance_squared(l, s, eps)); }
template<class T> Real distance(Segment<T> s, Line<T> l, Real eps = epsilon) { return distance(l, s, eps); }
template<class T> Real distance(Segment<T> s, Segment<T> t, Real eps = epsilon) { return std::sqrt(distance_squared(s, t, eps)); }
template<class T>
std::optional<Line<Real>> perpendicular_bisector(Point<T> a, Point<T> b) {
	if (a == b) return std::nullopt;
	using W = wide_t<T>;
	FPoint mid{Real(a.x) / 2 + Real(b.x) / 2, Real(a.y) / 2 + Real(b.y) / 2};
	FPoint d{Real(W(b.y) - a.y), Real(W(a.x) - b.x)};
	return Line<Real>{mid, mid + d};
}
}

#endif

Back to top page