convex_polygon.hpp¶
规范化凸多边形后以 O(log \(n\)) 判断点的位置，并提供有向直线裁剪与面积重心；支持退化点集和顺逆时针输入。
Complexity: Time: O(n) convex polygon preprocessing, O(log n) point query, O(n) clipping or centroid. Space: O(n), except O(1) for queries/centroid.
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(n) convex polygon preprocessing, O(log n) point query,
/// O(n) clipping or centroid. Space: O(n), except O(1) for queries/centroid.
#include "noya/geometry/geometry_base.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
namespace convex_polygon_internal {
template <class T>
auto turn(const point<T> &a, const point<T> &b, const point<T> &c) {
using W = std::conditional_t<std::numeric_limits<T>::is_integer,
__int128_t, long double>;
W bx = W(b.x) - W(a.x), by = W(b.y) - W(a.y);
W cx = W(c.x) - W(a.x), cy = W(c.y) - W(a.y);
return bx * cy - by * cx;
}
template <class T>
bool on_segment(const point<T> &p, const point<T> &a, const point<T> &b) {
return turn(a, b, p) == 0 && std::min(a.x, b.x) <= p.x &&
p.x <= std::max(a.x, b.x) && std::min(a.y, b.y) <= p.y &&
p.y <= std::max(a.y, b.y);
}
template <class T> point<long double> cast(const point<T> &p) {
return {static_cast<long double>(p.x), static_cast<long double>(p.y)};
}
} // namespace convex_polygon_internal
/// @brief Preprocess a convex boundary for logarithmic point location.
/// The input is in clockwise or counter-clockwise boundary order. Consecutive
/// duplicates, a repeated closing point, and redundant collinear vertices are
/// removed in O(n). Collinear input becomes its two extreme points.
/// For integral coordinates, all determinant products/differences must fit in
/// signed 128 bits; floating coordinates use exact long-double comparisons.
template <class T> class convex_polygon {
public:
convex_polygon() = default;
explicit convex_polygon(const std::vector<point<T>> &pg) { build(pg); }
void build(const std::vector<point<T>> &pg) {
using convex_polygon_internal::turn;
std::vector<point<T>> clean;
clean.reserve(pg.size());
for (const auto &p : pg) {
if (clean.empty() || !(clean.back() == p)) {
clean.push_back(p);
}
}
if (clean.size() > 1 && clean.front() == clean.back()) {
clean.pop_back();
}
int n = int(clean.size());
int winding = 0;
for (int i = 0; i < n; i++) {
auto z = turn(clean[i], clean[(i + 1) % n], clean[(i + 2) % n]);
if (z != 0) {
winding = z > 0 ? 1 : -1;
break;
}
}
if (winding == 0) {
hull.clear();
if (!clean.empty()) {
auto [lo, hi] = std::minmax_element(clean.begin(), clean.end());
hull.push_back(*lo);
if (!(*lo == *hi)) {
hull.push_back(*hi);
}
}
return;
}
if (winding < 0) {
std::reverse(clean.begin(), clean.end());
}
hull.clear();
hull.reserve(clean.size());
for (int i = 0; i < n; i++) {
auto z = turn(clean[(i + n - 1) % n], clean[i], clean[(i + 1) % n]);
assert(z >= 0); // The supplied boundary must be convex.
if (z != 0) {
hull.push_back(clean[i]);
}
}
auto first = std::min_element(hull.begin(), hull.end());
std::rotate(hull.begin(), first, hull.end());
}
const std::vector<point<T>> &vertices() const { return hull; }
int size() const { return int(hull.size()); }
bool empty() const { return hull.empty(); }
/// @brief Return -1 outside, 0 on the boundary, or 1 inside, in O(log n).
/// A degenerate point/segment has boundary points but no interior.
int locate(const point<T> &p) const {
using convex_polygon_internal::turn;
int n = size();
if (n == 0) {
return -1;
}
if (n == 1) {
return hull[0] == p ? 0 : -1;
}
if (n == 2) {
return convex_polygon_internal::on_segment(p, hull[0], hull[1]) ? 0 : -1;
}
auto left = turn(hull[0], hull[1], p);
auto right = turn(hull[0], hull[n - 1], p);
if (left < 0 || right > 0) {
return -1;
}
if (left == 0) {
return convex_polygon_internal::on_segment(p, hull[0], hull[1]) ? 0 : -1;
}
if (right == 0) {
return convex_polygon_internal::on_segment(p, hull[0], hull[n - 1]) ? 0 : -1;
}
int lo = 1, hi = n - 1;
while (hi - lo > 1) {
int mid = lo + (hi - lo) / 2;
if (turn(hull[0], hull[mid], p) >= 0) {
lo = mid;
} else {
hi = mid;
}
}
auto side = turn(hull[lo], hull[hi], p);
return (side > 0) - (side < 0);
}
private:
std::vector<point<T>> hull;
};
/// @brief Clip a convex boundary to the closed left half-plane of a -> b.
/// Accepts either winding and empty/point/segment boundaries. a != b is required.
/// eps is an absolute perpendicular-distance tolerance; the returned boundary
/// keeps the input winding and has no repeated closing point.
template <class T, class U>
std::vector<point<long double>> convex_polygon_cut(
const std::vector<point<T>> &pg, const point<U> &a, const point<U> &b,
long double eps = 1e-12L) {
using convex_polygon_internal::cast;
assert(eps >= 0);
auto origin = cast(a);
auto dir = cast(b) - origin;
long double len = std::hypot(dir.x, dir.y);
assert(len > 0);
dir /= len;
auto side = [&](const point<long double> &p) {
long double z = cross(dir, p - origin);
return std::abs(z) <= eps ? 0.0L : z;
};
std::vector<point<long double>> out;
auto append = [&](const point<long double> &p) {
if (out.empty() || distance(out.back(), p) > eps) {
out.push_back(p);
}
};
for (int i = 0; i < int(pg.size()); i++) {
auto p = cast(pg[i]);
auto q = cast(pg[(i + 1) % pg.size()]);
long double s = side(p), t = side(q);
if (s >= 0) {
append(p);
}
if ((s < 0 && t > 0) || (s > 0 && t < 0)) {
append(p + (q - p) * (s / (s - t)));
}
}
if (out.size() > 1 && distance(out.front(), out.back()) <= eps) {
out.pop_back();
}
return out;
}
/// @brief Return the area centroid of a simple polygon in either winding.
/// Returns nullopt when the absolute doubled signed area is <= eps.
/// Computations use long double and a translated origin to reduce cancellation.
template <class T>
std::optional<point<long double>> polygon_centroid(
const std::vector<point<T>> &pg, long double eps = 0) {
using convex_polygon_internal::cast;
assert(eps >= 0);
if (pg.size() < 3) {
return std::nullopt;
}
auto origin = cast(pg[0]);
point<long double> sum{};
long double area = 0;
for (int i = 1; i + 1 < int(pg.size()); i++) {
auto a = cast(pg[i]) - origin;
auto b = cast(pg[i + 1]) - origin;
long double weight = cross(a, b);
area += weight;
sum += (a + b) * weight;
}
if (std::abs(area) <= eps) {
return std::nullopt;
}
return origin + sum / (3 * area);
}
} // namespace noya
// NOYA_CONVEX_POLYGON_HPP
#ifndef NOYA_CONVEX_POLYGON_HPP
#define NOYA_CONVEX_POLYGON_HPP 1
/// @complexity Time: O(n) convex polygon preprocessing, O(log n) point query,
/// O(n) clipping or centroid. Space: O(n), except O(1) for queries/centroid.
#include "noya/geometry/geometry_base.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
namespace convex_polygon_internal {
template <class T>
auto turn(const point<T> &a, const point<T> &b, const point<T> &c) {
using W = std::conditional_t<std::numeric_limits<T>::is_integer,
__int128_t, long double>;
W bx = W(b.x) - W(a.x), by = W(b.y) - W(a.y);
W cx = W(c.x) - W(a.x), cy = W(c.y) - W(a.y);
return bx * cy - by * cx;
}
template <class T>
bool on_segment(const point<T> &p, const point<T> &a, const point<T> &b) {
return turn(a, b, p) == 0 && std::min(a.x, b.x) <= p.x &&
p.x <= std::max(a.x, b.x) && std::min(a.y, b.y) <= p.y &&
p.y <= std::max(a.y, b.y);
}
template <class T> point<long double> cast(const point<T> &p) {
return {static_cast<long double>(p.x), static_cast<long double>(p.y)};
}
} // namespace convex_polygon_internal
/// @brief Preprocess a convex boundary for logarithmic point location.
/// The input is in clockwise or counter-clockwise boundary order. Consecutive
/// duplicates, a repeated closing point, and redundant collinear vertices are
/// removed in O(n). Collinear input becomes its two extreme points.
/// For integral coordinates, all determinant products/differences must fit in
/// signed 128 bits; floating coordinates use exact long-double comparisons.
template <class T> class convex_polygon {
public:
convex_polygon() = default;
explicit convex_polygon(const std::vector<point<T>> &pg) { build(pg); }
void build(const std::vector<point<T>> &pg) {
using convex_polygon_internal::turn;
std::vector<point<T>> clean;
clean.reserve(pg.size());
for (const auto &p : pg) {
if (clean.empty() || !(clean.back() == p)) {
clean.push_back(p);
}
}
if (clean.size() > 1 && clean.front() == clean.back()) {
clean.pop_back();
}
int n = int(clean.size());
int winding = 0;
for (int i = 0; i < n; i++) {
auto z = turn(clean[i], clean[(i + 1) % n], clean[(i + 2) % n]);
if (z != 0) {
winding = z > 0 ? 1 : -1;
break;
}
}
if (winding == 0) {
hull.clear();
if (!clean.empty()) {
auto [lo, hi] = std::minmax_element(clean.begin(), clean.end());
hull.push_back(*lo);
if (!(*lo == *hi)) {
hull.push_back(*hi);
}
}
return;
}
if (winding < 0) {
std::reverse(clean.begin(), clean.end());
}
hull.clear();
hull.reserve(clean.size());
for (int i = 0; i < n; i++) {
auto z = turn(clean[(i + n - 1) % n], clean[i], clean[(i + 1) % n]);
assert(z >= 0); // The supplied boundary must be convex.
if (z != 0) {
hull.push_back(clean[i]);
}
}
auto first = std::min_element(hull.begin(), hull.end());
std::rotate(hull.begin(), first, hull.end());
}
const std::vector<point<T>> &vertices() const { return hull; }
int size() const { return int(hull.size()); }
bool empty() const { return hull.empty(); }
/// @brief Return -1 outside, 0 on the boundary, or 1 inside, in O(log n).
/// A degenerate point/segment has boundary points but no interior.
int locate(const point<T> &p) const {
using convex_polygon_internal::turn;
int n = size();
if (n == 0) {
return -1;
}
if (n == 1) {
return hull[0] == p ? 0 : -1;
}
if (n == 2) {
return convex_polygon_internal::on_segment(p, hull[0], hull[1]) ? 0 : -1;
}
auto left = turn(hull[0], hull[1], p);
auto right = turn(hull[0], hull[n - 1], p);
if (left < 0 || right > 0) {
return -1;
}
if (left == 0) {
return convex_polygon_internal::on_segment(p, hull[0], hull[1]) ? 0 : -1;
}
if (right == 0) {
return convex_polygon_internal::on_segment(p, hull[0], hull[n - 1]) ? 0 : -1;
}
int lo = 1, hi = n - 1;
while (hi - lo > 1) {
int mid = lo + (hi - lo) / 2;
if (turn(hull[0], hull[mid], p) >= 0) {
lo = mid;
} else {
hi = mid;
}
}
auto side = turn(hull[lo], hull[hi], p);
return (side > 0) - (side < 0);
}
private:
std::vector<point<T>> hull;
};
/// @brief Clip a convex boundary to the closed left half-plane of a -> b.
/// Accepts either winding and empty/point/segment boundaries. a != b is required.
/// eps is an absolute perpendicular-distance tolerance; the returned boundary
/// keeps the input winding and has no repeated closing point.
template <class T, class U>
std::vector<point<long double>> convex_polygon_cut(
const std::vector<point<T>> &pg, const point<U> &a, const point<U> &b,
long double eps = 1e-12L) {
using convex_polygon_internal::cast;
assert(eps >= 0);
auto origin = cast(a);
auto dir = cast(b) - origin;
long double len = std::hypot(dir.x, dir.y);
assert(len > 0);
dir /= len;
auto side = [&](const point<long double> &p) {
long double z = cross(dir, p - origin);
return std::abs(z) <= eps ? 0.0L : z;
};
std::vector<point<long double>> out;
auto append = [&](const point<long double> &p) {
if (out.empty() || distance(out.back(), p) > eps) {
out.push_back(p);
}
};
for (int i = 0; i < int(pg.size()); i++) {
auto p = cast(pg[i]);
auto q = cast(pg[(i + 1) % pg.size()]);
long double s = side(p), t = side(q);
if (s >= 0) {
append(p);
}
if ((s < 0 && t > 0) || (s > 0 && t < 0)) {
append(p + (q - p) * (s / (s - t)));
}
}
if (out.size() > 1 && distance(out.front(), out.back()) <= eps) {
out.pop_back();
}
return out;
}
/// @brief Return the area centroid of a simple polygon in either winding.
/// Returns nullopt when the absolute doubled signed area is <= eps.
/// Computations use long double and a translated origin to reduce cancellation.
template <class T>
std::optional<point<long double>> polygon_centroid(
const std::vector<point<T>> &pg, long double eps = 0) {
using convex_polygon_internal::cast;
assert(eps >= 0);
if (pg.size() < 3) {
return std::nullopt;
}
auto origin = cast(pg[0]);
point<long double> sum{};
long double area = 0;
for (int i = 1; i + 1 < int(pg.size()); i++) {
auto a = cast(pg[i]) - origin;
auto b = cast(pg[i + 1]) - origin;
long double weight = cross(a, b);
area += weight;
sum += (a + b) * weight;
}
if (std::abs(area) <= eps) {
return std::nullopt;
}
return origin + sum / (3 * area);
}
} // namespace noya
#endif // NOYA_CONVEX_POLYGON_HPP
/// @complexity Time: O(n) convex polygon preprocessing, O(log n) point query,
/// O(n) clipping or centroid. Space: O(n), except O(1) for queries/centroid.
/// @complexity Time: O(1) per primitive; O(n log n) for convex hull.
/// Space: O(1) per primitive and O(n) for hull construction.
#include <algorithm>
#include <cmath>
#include <optional>
#include <vector>
namespace noya {
/// @brief Two-dimensional point with vector arithmetic and lexicographic order.
template <class T> struct point {
T x{};
T y{};
point() = default;
point(T x_, T y_) : x(x_), y(y_) {}
point &operator+=(const point &rhs) {
x += rhs.x;
y += rhs.y;
return *this;
}
point &operator-=(const point &rhs) {
x -= rhs.x;
y -= rhs.y;
return *this;
}
point &operator*=(const T &scl) {
x *= scl;
y *= scl;
return *this;
}
point &operator/=(const T &scl) {
x /= scl;
y /= scl;
return *this;
}
friend point operator+(point l, const point &r) { return l += r; }
friend point operator-(point l, const point &r) { return l -= r; }
friend point operator*(point val, const T &scl) { return val *= scl; }
friend point operator*(const T &scl, point val) { return val *= scl; }
friend point operator/(point val, const T &scl) { return val /= scl; }
friend bool operator==(const point &, const point &) = default;
friend bool operator<(const point &l, const point &r) {
return l.x < r.x || (l.x == r.x && l.y < r.y);
}
};
/// @brief Circle represented by a center and a nonnegative radius.
template <class Real> struct circle {
point<Real> o;
Real r1{};
};
/// @brief Return the dot product of two vectors.
template <class T> T dot(const point<T> &a, const point<T> &b) {
return a.x * b.x + a.y * b.y;
}
/// @brief Return the signed cross product of two vectors.
template <class T> T cross(const point<T> &a, const point<T> &b) {
return a.x * b.y - a.y * b.x;
}
/// @brief Return cross(a - o1, b - o1).
template <class T>
T cross(const point<T> &o1, const point<T> &a, const point<T> &b) {
return cross(a - o1, b - o1);
}
/// @brief Return the squared Euclidean norm.
template <class T> T norm2(const point<T> &val) { return dot(val, val); }
/// @brief Compare a value with zero using an optional absolute tolerance.
template <class T> int sign(const T &val, const T &eps = T{}) {
return (val > eps) - (val < -eps);
}
/// @brief Return -1, 0, or 1 for a clockwise, collinear, or counter-clockwise
/// turn.
template <class T>
int orientation(const point<T> &a, const point<T> &b, const point<T> &c,
const T &eps = T{}) {
return sign(cross(a, b, c), eps);
}
/// @brief Test whether p lies on the closed segment [a, b].
template <class T>
bool on_segment(const point<T> &p, const point<T> &a, const point<T> &b,
const T &eps = T{}) {
if (orientation(a, b, p, eps) != 0) {
return false;
}
return std::min(a.x, b.x) - eps <= p.x && p.x <= std::max(a.x, b.x) + eps &&
std::min(a.y, b.y) - eps <= p.y && p.y <= std::max(a.y, b.y) + eps;
}
/// @brief Test whether the closed segments [a, b] and [c, d] intersect.
template <class T>
bool segments_intersect(const point<T> &a, const point<T> &b, const point<T> &c,
const point<T> &d, const T &eps = T{}) {
int s1 = orientation(a, b, c, eps);
int s2 = orientation(a, b, d, eps);
int s3 = orientation(c, d, a, eps);
int s4 = orientation(c, d, b, eps);
if (s1 == 0 && on_segment(c, a, b, eps)) {
return true;
}
if (s2 == 0 && on_segment(d, a, b, eps)) {
return true;
}
if (s3 == 0 && on_segment(a, c, d, eps)) {
return true;
}
if (s4 == 0 && on_segment(b, c, d, eps)) {
return true;
}
return s1 * s2 < 0 && s3 * s4 < 0;
}
/// @brief Intersect the infinite lines through (a, b) and (c, d), returning
/// nullopt when they are parallel or coincident.
template <class T>
std::optional<point<long double>>
line_intersection(const point<T> &a, const point<T> &b, const point<T> &c,
const point<T> &d, long double eps = 0) {
point<long double> lhs{static_cast<long double>(a.x),
static_cast<long double>(a.y)};
point<long double> b1{static_cast<long double>(b.x),
static_cast<long double>(b.y)};
point<long double> z{static_cast<long double>(c.x),
static_cast<long double>(c.y)};
point<long double> d1{static_cast<long double>(d.x),
static_cast<long double>(d.y)};
point<long double> da = b1 - lhs;
point<long double> db = d1 - z;
long double den = cross(da, db);
if (std::abs(den) <= eps) {
return std::nullopt;
}
long double rat = cross(z - lhs, db) / den;
return lhs + da * rat;
}
/// @brief Return the convex hull in counter-clockwise order without repetition.
template <class T>
std::vector<point<T>> convex_hull(std::vector<point<T>> pt, bool kc = false) {
std::sort(pt.begin(), pt.end());
pt.erase(std::unique(pt.begin(), pt.end()), pt.end());
if (pt.size() <= 1) {
return pt;
}
bool col = true;
for (int i = 2; i < int(pt.size()); i++) {
col &= orientation(pt[0], pt[1], pt[i]) == 0;
}
if (kc && col) {
return pt;
}
std::vector<point<T>> lo, hi;
for (const point<T> &p : pt) {
while (lo.size() >= 2) {
int rot = orientation(lo[lo.size() - 2], lo.back(), p);
if (rot > 0 || (kc && rot == 0)) {
break;
}
lo.pop_back();
}
lo.push_back(p);
}
for (auto it = pt.rbegin(); it != pt.rend(); ++it) {
while (hi.size() >= 2) {
int rot = orientation(hi[hi.size() - 2], hi.back(), *it);
if (rot > 0 || (kc && rot == 0)) {
break;
}
hi.pop_back();
}
hi.push_back(*it);
}
lo.pop_back();
hi.pop_back();
lo.insert(lo.end(), hi.begin(), hi.end());
return lo;
}
/// @brief Return twice the signed area of a polygon.
template <class T> T polygon_area2(const std::vector<point<T>> &pg) {
T res{};
for (int i = 0; i < int(pg.size()); i++) {
res += cross(pg[i], pg[(i + 1) % pg.size()]);
}
return res;
}
/// @brief Classify a point relative to a polygon: -1 outside, 0 boundary, 1
/// inside.
template <class T>
int point_in_polygon(const point<T> &p, const std::vector<point<T>> &pg) {
bool in = false;
for (int i = 0; i < int(pg.size()); i++) {
point<T> a = pg[i];
point<T> b = pg[(i + 1) % pg.size()];
if (on_segment(p, a, b)) {
return 0;
}
if (a.y <= p.y && p.y < b.y && orientation(a, b, p) > 0) {
in = !in;
}
if (b.y <= p.y && p.y < a.y && orientation(a, b, p) < 0) {
in = !in;
}
}
return in ? 1 : -1;
}
/// @brief Return the Euclidean distance between two points.
template <class T> long double distance(const point<T> &a, const point<T> &b) {
return std::hypot(
static_cast<long double>(a.x) - static_cast<long double>(b.x),
static_cast<long double>(a.y) - static_cast<long double>(b.y));
}
} // namespace noya
#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
namespace convex_polygon_internal {
template <class T>
auto turn(const point<T> &a, const point<T> &b, const point<T> &c) {
using W = std::conditional_t<std::numeric_limits<T>::is_integer,
__int128_t, long double>;
W bx = W(b.x) - W(a.x), by = W(b.y) - W(a.y);
W cx = W(c.x) - W(a.x), cy = W(c.y) - W(a.y);
return bx * cy - by * cx;
}
template <class T>
bool on_segment(const point<T> &p, const point<T> &a, const point<T> &b) {
return turn(a, b, p) == 0 && std::min(a.x, b.x) <= p.x &&
p.x <= std::max(a.x, b.x) && std::min(a.y, b.y) <= p.y &&
p.y <= std::max(a.y, b.y);
}
template <class T> point<long double> cast(const point<T> &p) {
return {static_cast<long double>(p.x), static_cast<long double>(p.y)};
}
} // namespace convex_polygon_internal
/// @brief Preprocess a convex boundary for logarithmic point location.
/// The input is in clockwise or counter-clockwise boundary order. Consecutive
/// duplicates, a repeated closing point, and redundant collinear vertices are
/// removed in O(n). Collinear input becomes its two extreme points.
/// For integral coordinates, all determinant products/differences must fit in
/// signed 128 bits; floating coordinates use exact long-double comparisons.
template <class T> class convex_polygon {
public:
convex_polygon() = default;
explicit convex_polygon(const std::vector<point<T>> &pg) { build(pg); }
void build(const std::vector<point<T>> &pg) {
using convex_polygon_internal::turn;
std::vector<point<T>> clean;
clean.reserve(pg.size());
for (const auto &p : pg) {
if (clean.empty() || !(clean.back() == p)) {
clean.push_back(p);
}
}
if (clean.size() > 1 && clean.front() == clean.back()) {
clean.pop_back();
}
int n = int(clean.size());
int winding = 0;
for (int i = 0; i < n; i++) {
auto z = turn(clean[i], clean[(i + 1) % n], clean[(i + 2) % n]);
if (z != 0) {
winding = z > 0 ? 1 : -1;
break;
}
}
if (winding == 0) {
hull.clear();
if (!clean.empty()) {
auto [lo, hi] = std::minmax_element(clean.begin(), clean.end());
hull.push_back(*lo);
if (!(*lo == *hi)) {
hull.push_back(*hi);
}
}
return;
}
if (winding < 0) {
std::reverse(clean.begin(), clean.end());
}
hull.clear();
hull.reserve(clean.size());
for (int i = 0; i < n; i++) {
auto z = turn(clean[(i + n - 1) % n], clean[i], clean[(i + 1) % n]);
assert(z >= 0); // The supplied boundary must be convex.
if (z != 0) {
hull.push_back(clean[i]);
}
}
auto first = std::min_element(hull.begin(), hull.end());
std::rotate(hull.begin(), first, hull.end());
}
const std::vector<point<T>> &vertices() const { return hull; }
int size() const { return int(hull.size()); }
bool empty() const { return hull.empty(); }
/// @brief Return -1 outside, 0 on the boundary, or 1 inside, in O(log n).
/// A degenerate point/segment has boundary points but no interior.
int locate(const point<T> &p) const {
using convex_polygon_internal::turn;
int n = size();
if (n == 0) {
return -1;
}
if (n == 1) {
return hull[0] == p ? 0 : -1;
}
if (n == 2) {
return convex_polygon_internal::on_segment(p, hull[0], hull[1]) ? 0 : -1;
}
auto left = turn(hull[0], hull[1], p);
auto right = turn(hull[0], hull[n - 1], p);
if (left < 0 || right > 0) {
return -1;
}
if (left == 0) {
return convex_polygon_internal::on_segment(p, hull[0], hull[1]) ? 0 : -1;
}
if (right == 0) {
return convex_polygon_internal::on_segment(p, hull[0], hull[n - 1]) ? 0 : -1;
}
int lo = 1, hi = n - 1;
while (hi - lo > 1) {
int mid = lo + (hi - lo) / 2;
if (turn(hull[0], hull[mid], p) >= 0) {
lo = mid;
} else {
hi = mid;
}
}
auto side = turn(hull[lo], hull[hi], p);
return (side > 0) - (side < 0);
}
private:
std::vector<point<T>> hull;
};
/// @brief Clip a convex boundary to the closed left half-plane of a -> b.
/// Accepts either winding and empty/point/segment boundaries. a != b is required.
/// eps is an absolute perpendicular-distance tolerance; the returned boundary
/// keeps the input winding and has no repeated closing point.
template <class T, class U>
std::vector<point<long double>> convex_polygon_cut(
const std::vector<point<T>> &pg, const point<U> &a, const point<U> &b,
long double eps = 1e-12L) {
using convex_polygon_internal::cast;
assert(eps >= 0);
auto origin = cast(a);
auto dir = cast(b) - origin;
long double len = std::hypot(dir.x, dir.y);
assert(len > 0);
dir /= len;
auto side = [&](const point<long double> &p) {
long double z = cross(dir, p - origin);
return std::abs(z) <= eps ? 0.0L : z;
};
std::vector<point<long double>> out;
auto append = [&](const point<long double> &p) {
if (out.empty() || distance(out.back(), p) > eps) {
out.push_back(p);
}
};
for (int i = 0; i < int(pg.size()); i++) {
auto p = cast(pg[i]);
auto q = cast(pg[(i + 1) % pg.size()]);
long double s = side(p), t = side(q);
if (s >= 0) {
append(p);
}
if ((s < 0 && t > 0) || (s > 0 && t < 0)) {
append(p + (q - p) * (s / (s - t)));
}
}
if (out.size() > 1 && distance(out.front(), out.back()) <= eps) {
out.pop_back();
}
return out;
}
/// @brief Return the area centroid of a simple polygon in either winding.
/// Returns nullopt when the absolute doubled signed area is <= eps.
/// Computations use long double and a translated origin to reduce cancellation.
template <class T>
std::optional<point<long double>> polygon_centroid(
const std::vector<point<T>> &pg, long double eps = 0) {
using convex_polygon_internal::cast;
assert(eps >= 0);
if (pg.size() < 3) {
return std::nullopt;
}
auto origin = cast(pg[0]);
point<long double> sum{};
long double area = 0;
for (int i = 1; i + 1 < int(pg.size()); i++) {
auto a = cast(pg[i]) - origin;
auto b = cast(pg[i + 1]) - origin;
long double weight = cross(a, b);
area += weight;
sum += (a + b) * weight;
}
if (std::abs(area) <= eps) {
return std::nullopt;
}
return origin + sum / (3 * area);
}
} // namespace noya
// NOYA_CONVEX_POLYGON_HPP