cp-library

競技プログラミング用 C++ ライブラリ

View the Project on GitHub qandori/cp-library

:heavy_check_mark: 多角形の面積・内外判定・切断

geometry/polygon.hpp · GitHub で見る

共通ページ 

#include "lib/geometry/polygon.hpp"

using namespace geometry;
std::vector<IPoint> polygon{ {0, 0}, {4, 0}, {4, 3}, {0, 3} };
auto twice = signed_area2(polygon); // 24
auto where = contains(polygon, IPoint(2, 1)); // Location::inside
auto center = centroid(polygon); // (2, 1.5)
auto left = convex_cut(polygon, Line<std::int64_t>{ {2, -1}, {2, 4} });
auto left_area = area(left); // 6

多角形は境界に沿った頂点列 vector<Point<T>> で渡します。末尾に始点を重複させる必要はありません。 通常は反時計回りを使いますが、面積・周長・重心・一般の内外判定は時計回りでも使えます。 T は符号付き整数または浮動小数点数とし、整数なら IPoint、実数なら FPoint が使えます。

操作 内容・返り値 計算量
signed_area2(polygon) 符号付き面積の2倍を wide_t<T> で返す。反時計回りが正 \(O(n)\)
area(polygon) 面積の絶対値を Real で返す \(O(n)\)
perimeter(polygon) 最後から最初への辺を含む周長を Real で返す \(O(n)\)
is_convex(polygon,strict = false,eps = epsilon) 単純多角形が凸か判定。strict = true では同一直線上の連続する3頂点も認めない \(O(n)\)
contains(polygon,p,eps = epsilon) 一般の単純多角形に対する Location::outside・boundary・inside \(O(n)\)
centroid(polygon,eps = epsilon) 面積についての重心を optional<FPoint> で返す。面積0は nullopt \(O(n)\)
convex_cut(polygon,line,eps = epsilon) 凸多角形を有向直線 line.a -> line.b の左側で切り、残る頂点列を vector<FPoint> で返す。直線上も残す \(O(n)\)

is_convex は入力が単純であることを前提にします。自己交差の検査はしません。0・1・2頂点は、strict = false なら凸、true なら凸でないものとします。 contains は頂点・辺上を boundary と区別します。空の多角形は常に outside、単点・線分には内部がなく、そこに乗る点だけが boundary です。

空・単点の面積と周長は0です。2点の周長は線分長の2倍です。 centroid は頂点の単純平均ではなく、囲まれた領域の重心です。 convex_cut には凸な境界順の頂点列と、異なる2点からなる直線を渡します。切断後が空・単点・線分になることもあります。

整数の符号判定と signed_area2 は浮動小数点に変換せず計算します。計算途中の積・和は符号付き128 bitに収めてください。 area・perimeter・centroid・切断後の座標は実数です。浮動小数点の eps は絶対誤差のしきい値なので、座標の大きさに合わせて変更してください。整数では eps を使いません。centroid の面積と退化判定も、整数入力なら128 bitで計算してから実数に変換します。

頂点数を \(n\) とすると、各操作は \(O(n)\) 時間です。convex_cut は返り値を含め \(O(n)\) 空間、それ以外は返り値を除き \(O(1)\) 空間です。

依存ライブラリ

検証コード

コード

#ifndef GEOMETRY_POLYGON_HPP
#define GEOMETRY_POLYGON_HPP

#include "linear.hpp"
#include <cassert>
#include <cmath>
#include <optional>
#include <vector>

namespace geometry {
template <class T>
wide_t<T> signed_area2(const std::vector<Point<T>> &polygon) {
	wide_t<T> result = 0;
	for (std::size_t i = 2; i < polygon.size(); ++i) result += cross(polygon[0], polygon[i - 1], polygon[i]);
	return result;
}

template <class T>
Real area(const std::vector<Point<T>> &polygon) { return std::abs(Real(signed_area2(polygon))) / 2; }

template <class T>
Real perimeter(const std::vector<Point<T>> &polygon) {
	Real result = 0;
	for (std::size_t i = 0; i < polygon.size(); ++i) result += distance(polygon[i], polygon[(i + 1) % polygon.size()]);
	return result;
}

template <class T>
bool is_convex(const std::vector<Point<T>> &polygon, bool strict = false, Real eps = epsilon) {
	if (polygon.size() < 3) return !strict;
	int direction = 0;
	for (std::size_t i = 0; i < polygon.size(); ++i) {
		int turn = orientation(polygon[i], polygon[(i + 1) % polygon.size()], polygon[(i + 2) % polygon.size()], eps);
		if (!turn) { if (strict) return false; }
		else if (direction && turn != direction) return false;
		else direction = turn;
	}
	return true;
}

template <class T>
Location contains(const std::vector<Point<T>> &polygon, Point<T> p, Real eps = epsilon) {
	int winding = 0;
	for (std::size_t i = 0; i < polygon.size(); ++i) {
		const auto &a = polygon[i], &b = polygon[(i + 1) % polygon.size()];
		if (on_segment(Segment<T>{a, b}, p, eps)) return Location::boundary;
		int ay = sign(wide_t<T>(a.y) - p.y, eps), by = sign(wide_t<T>(b.y) - p.y, eps);
		int turn = orientation(a, b, p, eps);
		if (ay <= 0 && by > 0 && turn > 0) ++winding;
		if (ay > 0 && by <= 0 && turn < 0) --winding;
	}
	return winding ? Location::inside : Location::outside;
}

template <class T>
std::optional<FPoint> centroid(const std::vector<Point<T>> &polygon, Real eps = epsilon) {
	if (polygon.size() < 3) return std::nullopt;
	using W = wide_t<T>;
	FPoint origin = to_real(polygon[0]), sum{};
	W total = 0;
	for (std::size_t i = 2; i < polygon.size(); ++i) {
		FPoint a{Real(W(polygon[i - 1].x) - polygon[0].x), Real(W(polygon[i - 1].y) - polygon[0].y)};
		FPoint b{Real(W(polygon[i].x) - polygon[0].x), Real(W(polygon[i].y) - polygon[0].y)};
		W weight = cross(polygon[0], polygon[i - 1], polygon[i]);
		sum += (a + b) * Real(weight);
		total += weight;
	}
	if (!sign(total, eps)) return std::nullopt;
	return origin + sum / (3 * Real(total));
}

template <class T>
std::vector<FPoint> convex_cut(const std::vector<Point<T>> &polygon, Line<T> line, Real eps = epsilon) {
	assert(line.a != line.b);
	std::vector<FPoint> result;
	result.reserve(polygon.size() + 1);
	auto append = [&](FPoint p) { if (result.empty() || result.back() != p) result.push_back(p); };
	for (std::size_t i = 0; i < polygon.size(); ++i) {
		const auto &a = polygon[i], &b = polygon[(i + 1) % polygon.size()];
		auto ca = cross(line.a, line.b, a), cb = cross(line.a, line.b, b);
		int sa = sign(ca, eps), sb = sign(cb, eps);
		if (sa >= 0) append(to_real(a));
		if (sa * sb < 0) append(to_real(a) + (to_real(b) - to_real(a)) * (Real(ca) / (Real(ca) - Real(cb))));
	}
	if (result.size() > 1 && result.front() == result.back()) result.pop_back();
	return result;
}
}

#endif
Back to top page