凸包・凸多角形・最近点対・最遠点対

geometry/convex.hpp · GitHub で見る

共通ページ

#include "lib/geometry/convex.hpp"

using namespace geometry;
std::vector<IPoint> points{ {0, 0}, {4, 0}, {4, 3}, {0, 3}, {1, 1} };
auto hull = convex_hull(points);
auto where = convex_contains(hull, IPoint(2, 1)); // Location::inside
auto near = closest_pair(points);
if (near) {
    auto p = points[near->first];
    auto q = points[near->second];
    auto squared = near->squared_distance;
}
auto far = furthest_pair(points);
auto sum = minkowski_sum(hull, std::vector<IPoint>{ {0, 0}, {1, 0} });

点群には vector<Point<T>> を渡します。整数なら IPoint、実数なら FPoint を使えます。 convex_hull が返す境界には内点・重複点を含めず、始点を末尾に再び追加しません。 左端、同じ横座標なら最も下の頂点から始まり、反時計回りに並びます。

操作	内容・返り値	計算量
convex_hull_indices(points,keep_collinear = false,eps = epsilon)	凸包上の点の元の添字を vector<int> で返す。同じ座標が複数あれば、そのうちの1点を使う	\(O(n\log n)\)
convex_hull(points,keep_collinear = false,eps = epsilon)	凸包の頂点列を返す	\(O(n\log n)\)
convex_contains(convex,p,eps = epsilon)	厳密凸・反時計回りの頂点列に対する Location::outside・boundary・inside	\(O(\log h)\)
convex_diameter(convex,eps = epsilon)	厳密凸・反時計回りの頂点列の最遠頂点対を optional<PointPair<T>> で返す。添字は渡した頂点列内	\(O(h)\)
diameter_squared(convex,eps = epsilon)	同じ凸多角形の直径の2乗を wide_t<T> で返す。空・単点は0	\(O(h)\)
furthest_pair(points,eps = epsilon)	一般の点群の最遠点対。添字は元の点群内	\(O(n\log n)\)
closest_pair(points)	一般の点群の最近点対。添字は元の点群内。重複点があれば距離0の異なる添字を返す	\(O(n\log n)\)
minkowski_sum(a,b)	2つの凸多角形の Minkowski 和の頂点列を返す	\(O(n+m)\)
PointPair<T>::first / second	点対の添字 int。同距離の候補が複数あれば任意の1組	\(O(1)\)
PointPair<T>::squared_distance	点対の距離の2乗。整数なら符号付き128 bit、実数なら Real	\(O(1)\)

keep_collinear = false では、凸包の辺の途中の点を除きます。true なら辺上の点も残します。 全点が同一直線上の場合は、既定では両端のみ、true では全ての異なる座標を辞書順に返します。全点が同一なら1点、入力が空なら空配列です。

convex_contains と convex_diameter には、3頂点以上なら、同一直線上の連続する3頂点を含まない反時計回りの凸多角形を渡します。 前処理として convex_hull(points) を使えばこの形になります。convex_hull(points,true) の結果をそのまま渡せるとは限りません。 凸であることの検査は省略し、先頭の頂点は任意です。0・1・2頂点もそのまま扱えます。

closest_pair・furthest_pair は2点未満なら nullopt を返します。 convex_diameter は空なら nullopt、単点なら first = second = 0、距離0を返します。 2点以上の一般点群では、全て同じ座標でも、異なる添字の点対を返します。

Minkowski 和は、a の各点と b の各点を足して得られる点群の凸包です。 入力は凸な境界順であれば時計回り・辺上の余分な点・末尾の始点の重複を認めます。空集合との和は空、単点との和は平行移動です。 返り値は辺上の余分な点を除いた反時計回りの凸包で、始点の選び方は convex_hull と同じです。 各入力が凸かどうかは検査しません。非凸な点群は、先にそれぞれ convex_hull に通します。出力座標の和も T に収めてください。

整数では符号判定・距離の2乗を厳密計算し、途中の積・和を符号付き128 bitに収めます。 ソートと重複点の除去には、浮動小数点でも誤差しきい値を使いません。Minkowski 和の向き・辺方向の比較と、辺上の点の除去も0との比較で行い、小さい凸多角形もそのまま扱います。 その他の実数の判定には eps を使いますが、最近点対には使わず距離の大小をそのまま比較します。 座標は有限値とし、実数演算の途中の値も有限に収めてください。

点群の大きさを \(n\)、凸多角形の頂点数を \(h\) とします。 凸包・一般の最近点対・最遠点対は \(O(n\log n)\) 時間・\(O(n)\) 追加空間、凸多角形の内外判定は \(O(\log h)\) 時間・\(O(1)\) 追加空間、直径は \(O(h)\) 時間・\(O(1)\) 追加空間です。 Minkowski 和は入力の大きさを \(n,m\) として \(O(n+m)\) 時間・空間です。空・単点の判定は \(O(1)\) 時間です。

依存ライブラリ

• Closest Pair Core
• 整数点の凸包（補助部品）
• 64 bit キーの安定基数ソート
• SplitMix64
• 幾何の直線・線分
• 幾何の点・ベクトル・偏角順
• 多角形の面積・内外判定・切断

検証コード

• test/geometry_closest_pair.test.cpp
• test/geometry_convex_hull.test.cpp
• test/geometry_furthest_pair.test.cpp

コード

#ifndef GEOMETRY_CONVEX_HPP
#define GEOMETRY_CONVEX_HPP

#include "polygon.hpp"
#include "../components/closest_pair.hpp"
#include "../components/packed_convex_hull.hpp"
#include <algorithm>
#include <cassert>
#include <limits>
#include <numeric>
#include <optional>
#include <vector>

namespace geometry {
template <class T>
struct PointPair {
	int first, second;
	wide_t<T> squared_distance;
};

template <class T>
std::vector<int> convex_hull_indices(const std::vector<Point<T>> &points, bool keep_collinear = false, Real eps = epsilon) {
	assert(points.size() <= std::size_t(std::numeric_limits<int>::max()));
	if constexpr (integer<T> && sizeof(T) <= 8) {
		if (points.size() >= 4) {
			T min_x = points[0].x, max_x = min_x, min_y = points[0].y, max_y = min_y;
			for (auto p : points) {
				min_x = std::min(min_x, p.x), max_x = std::max(max_x, p.x);
				min_y = std::min(min_y, p.y), max_y = std::max(max_y, p.y);
			}
			if (wide_t<T>(max_x) - min_x <= INT32_MAX && wide_t<T>(max_y) - min_y <= INT32_MAX) {
				std::vector<packed_convex_hull_detail::Entry> order(points.size());
				for (int i = 0; i < int(points.size()); ++i)
					order[i] = {(std::uint64_t(wide_t<T>(points[i].x) - min_x) << 32) | std::uint32_t(wide_t<T>(points[i].y) - min_y), i};
				auto hull = packed_convex_hull_detail::build(std::move(order), keep_collinear);
				std::vector<int> result;
				result.reserve(hull.size());
				for (auto value : hull) result.push_back(value.index);
				return result;
			}
		}
	}
	std::vector<int> order(points.size());
	std::iota(order.begin(), order.end(), 0);
	std::sort(order.begin(), order.end(), [&](int a, int b) { return points[a] < points[b]; });
	order.erase(std::unique(order.begin(), order.end(), [&](int a, int b) { return points[a] == points[b]; }), order.end());
	if (order.size() <= 2) return order;
	if (keep_collinear && std::all_of(order.begin() + 2, order.end(), [&](int i) {
		return !orientation(points[order[0]], points[order[1]], points[i], eps);
	})) return order;
	std::vector<int> hull;
	hull.reserve(order.size() * 2);
	auto append = [&](int i, std::size_t base) {
		while (hull.size() >= base) {
			int turn = orientation(points[hull[hull.size() - 2]], points[hull.back()], points[i], eps);
			if (turn > 0 || (keep_collinear && !turn)) break;
			hull.pop_back();
		}
		hull.push_back(i);
	};
	for (int i : order) append(i, 2);
	std::size_t base = hull.size() + 1;
	for (std::size_t i = order.size() - 1; i-- > 0;) append(order[i], base);
	hull.pop_back();
	return hull;
}

template <class T>
std::vector<Point<T>> convex_hull(std::vector<Point<T>> &&points, bool keep_collinear = false, Real eps = epsilon) {
	assert(points.size() <= std::size_t(std::numeric_limits<int>::max()));
	if constexpr (integer<T> && sizeof(T) <= 8) {
		if (points.size() >= 4) {
			T min_x = points[0].x, max_x = min_x, min_y = points[0].y, max_y = min_y;
			for (auto p : points) {
				min_x = std::min(min_x, p.x), max_x = std::max(max_x, p.x);
				min_y = std::min(min_y, p.y), max_y = std::max(max_y, p.y);
			}
			if (wide_t<T>(max_x) - min_x <= INT32_MAX && wide_t<T>(max_y) - min_y <= INT32_MAX) {
				std::vector<std::uint64_t> order(points.size());
				for (int i = 0; i < int(points.size()); ++i)
					order[i] = (std::uint64_t(wide_t<T>(points[i].x) - min_x) << 32) | std::uint32_t(wide_t<T>(points[i].y) - min_y);
				auto hull = packed_convex_hull_detail::build(std::move(order), keep_collinear);
				points.clear();
				for (auto value : hull) points.emplace_back(T(wide_t<T>(min_x) + (value >> 32)), T(wide_t<T>(min_y) + std::uint32_t(value)));
				return std::move(points);
			}
		}
	}
	std::sort(points.begin(), points.end());
	points.erase(std::unique(points.begin(), points.end()), points.end());
	if (points.size() <= 2 || (keep_collinear && std::all_of(points.begin() + 2, points.end(), [&](auto p) {
		return !orientation(points[0], points[1], p, eps);
	}))) return std::move(points);
	std::vector<Point<T>> hull;
	hull.reserve(points.size() * 2);
	auto append = [&](Point<T> p, std::size_t base) {
		while (hull.size() >= base) {
			int turn = orientation(hull[hull.size() - 2], hull.back(), p, eps);
			if (turn > 0 || (keep_collinear && !turn)) break;
			hull.pop_back();
		}
		hull.push_back(p);
	};
	for (auto p : points) append(p, 2);
	std::size_t base = hull.size() + 1;
	for (std::size_t i = points.size() - 1; i-- > 0;) append(points[i], base);
	hull.pop_back();
	return hull;
}

template <class T>
std::vector<Point<T>> convex_hull(const std::vector<Point<T>> &points, bool keep_collinear = false, Real eps = epsilon) {
	return convex_hull(std::vector<Point<T>>(points), keep_collinear, eps);
}

template <class T>
Location convex_contains(const std::vector<Point<T>> &polygon, Point<T> p, Real eps = epsilon) {
	if (polygon.empty()) return Location::outside;
	if (polygon.size() <= 2) return on_segment(Segment<T>{polygon.front(), polygon.back()}, p, eps) ? Location::boundary : Location::outside;
	int first = orientation(polygon[0], polygon[1], p, eps), last = orientation(polygon[0], polygon.back(), p, eps);
	if (first < 0 || last > 0) return Location::outside;
	if (!first) return on_segment(Segment<T>{polygon[0], polygon[1]}, p, eps) ? Location::boundary : Location::outside;
	if (!last) return on_segment(Segment<T>{polygon[0], polygon.back()}, p, eps) ? Location::boundary : Location::outside;
	std::size_t l = 1, r = polygon.size() - 1;
	while (r - l > 1) {
		std::size_t m = (l + r) / 2;
		if (orientation(polygon[0], polygon[m], p, eps) >= 0) l = m;
		else r = m;
	}
	int side = orientation(polygon[l], polygon[r], p, eps);
	if (side > 0) return Location::inside;
	return !side && on_segment(Segment<T>{polygon[l], polygon[r]}, p, eps) ? Location::boundary : Location::outside;
}

template <class T>
std::optional<PointPair<T>> convex_diameter(const std::vector<Point<T>> &polygon, Real eps = epsilon) {
	assert(polygon.size() <= std::size_t(std::numeric_limits<int>::max()));
	int n = int(polygon.size());
	if (!n) return std::nullopt;
	if (n == 1) return PointPair<T>{0, 0, 0};
	PointPair<T> result{0, 1, distance_squared(polygon[0], polygon[1])};
	auto update = [&](int a, int b) {
		auto d = distance_squared(polygon[a], polygon[b]);
		if (d > result.squared_distance) result = {a, b, d};
	};
	for (int i = 0, j = 1; i < n; ++i) {
		int next = (i + 1) % n;
		while (sign(cross(polygon[i], polygon[next], polygon[(j + 1) % n]) - cross(polygon[i], polygon[next], polygon[j]), eps) > 0) j = (j + 1) % n;
		update(i, j);
		update(next, j);
	}
	return result;
}

template <class T>
wide_t<T> diameter_squared(const std::vector<Point<T>> &polygon, Real eps = epsilon) {
	auto result = convex_diameter(polygon, eps);
	return result ? result->squared_distance : wide_t<T>(0);
}

template <class T>
std::optional<PointPair<T>> furthest_pair(const std::vector<Point<T>> &points, Real eps = epsilon) {
	assert(points.size() <= std::size_t(std::numeric_limits<int>::max()));
	if (points.size() < 2) return std::nullopt;
	if constexpr (integer<T> && sizeof(T) <= 8) {
		if (points.size() >= 64) {
			T min_x = points[0].x, max_x = min_x, min_y = points[0].y, max_y = min_y;
			for (auto p : points) {
				min_x = std::min(min_x, p.x), max_x = std::max(max_x, p.x);
				min_y = std::min(min_y, p.y), max_y = std::max(max_y, p.y);
			}
			if (wide_t<T>(max_x) - min_x <= INT32_MAX && wide_t<T>(max_y) - min_y <= INT32_MAX) {
				auto pack = [&](Point<T> p) {
					return (std::uint64_t(wide_t<T>(p.x) - min_x) << 32) | std::uint32_t(wide_t<T>(p.y) - min_y);
				};
				std::vector<std::uint64_t> order(points.size());
				for (std::size_t i = 0; i < points.size(); ++i) order[i] = pack(points[i]);
				auto hull = packed_convex_hull_detail::build(std::move(order), false);
				if (hull.size() == 1) return PointPair<T>{0, 1, 0};
				auto [a, b] = packed_convex_hull_detail::diameter(hull);
				int first = -1, second = -1;
				for (int i = 0; i < int(points.size()) && (first < 0 || second < 0); ++i) {
					auto value = pack(points[i]);
					if (value == a) first = i;
					if (value == b) second = i;
				}
				return PointPair<T>{first, second, distance_squared(points[first], points[second])};
			}
		}
	}
	auto polygon = convex_hull(points, false, eps);
	if (polygon.size() == 1) return PointPair<T>{0, 1, 0};
	auto result = *convex_diameter(polygon, eps);
	auto a = polygon[result.first], b = polygon[result.second];
	result.first = result.second = -1;
	for (int i = 0; i < int(points.size()); ++i) {
		if (points[i] == a) result.first = i;
		if (points[i] == b) result.second = i;
		if (result.first >= 0 && result.second >= 0) break;
	}
	return result;
}

template <class T>
std::vector<Point<T>> minkowski_sum(std::vector<Point<T>> a, std::vector<Point<T>> b) {
	auto normalize = [&](std::vector<Point<T>> &p) {
		p.erase(std::unique(p.begin(), p.end()), p.end());
		if (p.size() > 1 && p.front() == p.back()) p.pop_back();
		if (p.size() <= 1) return;
		int direction = sign(signed_area2(p), 0);
		if (!direction) {
			auto [lo, hi] = std::minmax_element(p.begin(), p.end());
			Point<T> first = *lo, last = *hi;
			p = {first};
			if (first != last) p.push_back(last);
		} else {
			if (direction < 0) std::reverse(p.begin(), p.end());
			std::vector<Point<T>> strict;
			strict.reserve(p.size());
			for (std::size_t i = 0; i < p.size(); ++i)
				if (orientation(p[(i + p.size() - 1) % p.size()], p[i], p[(i + 1) % p.size()], 0)) strict.push_back(p[i]);
			p = std::move(strict);
		}
		if (!p.empty()) std::rotate(p.begin(), std::min_element(p.begin(), p.end(), [](auto x, auto y) { return x.y != y.y ? x.y < y.y : x.x < y.x; }), p.end());
	};
	normalize(a), normalize(b);
	if (a.empty() || b.empty()) return {};
	std::vector<Point<T>> result;
	result.reserve(a.size() + b.size());
	for (std::size_t i = 0, j = 0; i < a.size() || j < b.size();) {
		result.push_back(a[i % a.size()] + b[j % b.size()]);
		if (i == a.size() || a.size() == 1) { ++j; i = a.size(); continue; }
		if (j == b.size() || b.size() == 1) { ++i; j = b.size(); continue; }
		const auto &p = a[i], &q = a[(i + 1) % a.size()], &r = b[j], &s = b[(j + 1) % b.size()];
		wide_t<T> ax = wide_t<T>(q.x) - p.x, ay = wide_t<T>(q.y) - p.y, bx = wide_t<T>(s.x) - r.x, by = wide_t<T>(s.y) - r.y;
		bool ha = ay < 0 || (ay == 0 && ax < 0), hb = by < 0 || (by == 0 && bx < 0);
		int turn = ha != hb ? (ha ? -1 : 1) : sign(ax * by - ay * bx, 0);
		if (turn >= 0) ++i;
		if (turn <= 0) ++j;
	}
	if (result.size() > 1 && result.front() == result.back()) result.pop_back();
	std::rotate(result.begin(), std::min_element(result.begin(), result.end()), result.end());
	return result;
}

template <class T>
std::optional<PointPair<T>> closest_pair(const std::vector<Point<T>> &points) {
	assert(points.size() <= std::size_t(std::numeric_limits<int>::max()));
	int n = int(points.size());
	if (n < 2) return std::nullopt;
	if (n <= 32) {
		PointPair<T> result{0, 1, distance_squared(points[0], points[1])};
		for (int i = 1; i < n; ++i) for (int j = 0; j < i; ++j) {
			auto d = distance_squared(points[i], points[j]);
			if (d < result.squared_distance) result = {i, j, d};
		}
		return result;
	}
	if constexpr (integer<T> && sizeof(T) <= 8) {
		if (n >= 64) {
			T min_x = points[0].x, max_x = min_x, min_y = points[0].y, max_y = min_y;
			for (auto p : points) {
				min_x = std::min(min_x, p.x), max_x = std::max(max_x, p.x);
				min_y = std::min(min_y, p.y), max_y = std::max(max_y, p.y);
			}
			if (wide_t<T>(max_x) - min_x <= INT32_MAX && wide_t<T>(max_y) - min_y <= INT32_MAX) {
				auto pack = [&](Point<T> p) {
					return (std::uint64_t(wide_t<T>(p.y) - min_y) << 32) | std::uint32_t(wide_t<T>(p.x) - min_x);
				};
				std::vector<std::uint64_t> packed(n);
				for (int i = 0; i < n; ++i) packed[i] = pack(points[i]);
				if (auto pair = closest_pair_detail::packed(std::move(packed), std::uint32_t(wide_t<T>(max_x) - min_x), std::uint64_t(n) * 32)) {
					int a = -1, b = -1;
					for (int i = 0; i < n && (a < 0 || b < 0); ++i) {
						auto value = pack(points[i]);
						if (value == pair->first && a < 0) a = i;
						else if (value == pair->second && b < 0) b = i;
					}
					return PointPair<T>{a, b, distance_squared(points[a], points[b])};
				}
			}
		}
	}
	struct Entry { Point<T> point; int index; };
	std::vector<Entry> sorted, buffer(n);
	sorted.reserve(n);
	for (int i = 0; i < n; ++i) sorted.push_back({points[i], i});
	std::sort(sorted.begin(), sorted.end(), [](const auto &a, const auto &b) { return a.point < b.point; });
	for (int i = 1; i < n; ++i)
		if (sorted[i - 1].point == sorted[i].point) return PointPair<T>{sorted[i - 1].index, sorted[i].index, 0};
	PointPair<T> result{sorted[0].index, sorted[1].index, distance_squared(sorted[0].point, sorted[1].point)};
	auto update = [&](const Entry &a, const Entry &b) {
		auto d = distance_squared(a.point, b.point);
		if (d < result.squared_distance) result = {a.index, b.index, d};
	};
	auto by_y = [](const auto &a, const auto &b) { return a.point.y != b.point.y ? a.point.y < b.point.y : a.point.x < b.point.x; };
	auto solve = [&](auto &&self, int l, int r) -> void {
		if (r - l <= 3) {
			for (int i = l; i < r; ++i) for (int j = i + 1; j < r; ++j) update(sorted[i], sorted[j]);
			std::sort(sorted.begin() + l, sorted.begin() + r, by_y);
			return;
		}
		int m = (l + r) / 2;
		T middle_x = sorted[m].point.x;
		self(self, l, m), self(self, m, r);
		std::merge(sorted.begin() + l, sorted.begin() + m, sorted.begin() + m, sorted.begin() + r, buffer.begin() + l, by_y);
		std::copy(buffer.begin() + l, buffer.begin() + r, sorted.begin() + l);
		int count = 0;
		for (int i = l; i < r; ++i) {
			wide_t<T> dx = wide_t<T>(sorted[i].point.x) - middle_x;
			if (dx * dx >= result.squared_distance) continue;
			for (int j = count; j-- > 0;) {
				wide_t<T> dy = wide_t<T>(sorted[i].point.y) - buffer[l + j].point.y;
				if (dy * dy >= result.squared_distance) break;
				update(sorted[i], buffer[l + j]);
			}
			buffer[l + count++] = sorted[i];
		}
	};
	solve(solve, 0, n);
	return result;
}
}

#endif

Back to top page