convex_layers.hpp¶
逐层删除凸包边界，返回每个点从 1 开始的层号；输入点须互异，行列式与带权高度比较须能用 128 位有符号整数表示。
Complexity: Time: \(O(n \log^{2} n)\). Space: \(O(n)\).
Source: Library Checker
AC 记录：convex_layers。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O(n log^2 n).
/// Space: O(n).
#include "noya/geometry/geometry_base.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iterator>
#include <utility>
#include <vector>
namespace noya {
namespace convex_layers_detail {
using point_type = point<std::int64_t>;
class left_hull {
struct node {
int bl = 0;
int br = 0;
int low = 0;
int hi = 0;
int l = -1;
int r = -1;
};
std::vector<point_type> pt_;
std::vector<node> nd_;
int rt_ = 0;
bool rev_ = false;
static __int128 turn(const point_type &a, const point_type &b,
const point_type &c) {
__int128 x1 = __int128(b.x) - a.x, y1 = __int128(b.y) - a.y;
__int128 x2 = __int128(c.x) - a.x, y2 = __int128(c.y) - a.y;
return x1 * y2 - y1 * x2;
}
bool leaf(int cur) const { return nd_[cur].l == -1 && nd_[cur].r == -1; }
void pull(int cur) {
int l = nd_[cur].l;
int r = nd_[cur].r;
std::int64_t sy = pt_[nd_[r].low].y;
while (!leaf(l) || !leaf(r)) {
int a = nd_[l].bl;
int b = nd_[l].br;
int c = nd_[r].bl;
int d = nd_[r].br;
if (a != b && turn(pt_[a], pt_[b], pt_[c]) > 0) {
l = nd_[l].l;
} else if (c != d && turn(pt_[b], pt_[c], pt_[d]) > 0) {
r = nd_[r].r;
} else if (a == b) {
r = nd_[r].l;
} else if (c == d) {
l = nd_[l].r;
} else {
__int128 lhs = turn(pt_[a], pt_[b], pt_[c]);
__int128 rhs = turn(pt_[b], pt_[a], pt_[d]);
assert(lhs + rhs >= 0);
// Mirroring both coordinates preserves turns and reverses this comparison.
__int128 cmp = lhs * (__int128(pt_[d].y) - sy) +
rhs * (__int128(pt_[c].y) - sy);
if (lhs + rhs == 0 || (rev_ ? cmp > 0 : cmp < 0)) {
l = nd_[l].r;
} else {
r = nd_[r].l;
}
}
}
nd_[cur].bl = nd_[l].low;
nd_[cur].br = nd_[r].low;
}
void build(int cur, int low, int hi) {
nd_[cur].low = low;
nd_[cur].hi = hi;
if (hi - low == 1) {
nd_[cur].bl = nd_[cur].br = low;
nd_[cur].l = nd_[cur].r = -1;
return;
}
int mid = (low + hi) / 2;
nd_[cur].l = cur + 1;
nd_[cur].r = cur + 2 * (mid - low);
build(nd_[cur].l, low, mid);
build(nd_[cur].r, mid, hi);
pull(cur);
}
// Delete a sorted batch, rebuilding each affected bridge only once.
int erase(int cur, const std::vector<int> &pos, int l, int r) {
if (l == r) return cur;
if (r - l == nd_[cur].hi - nd_[cur].low) return -1;
int cut = nd_[nd_[cur].r].low;
int mid = int(std::lower_bound(pos.begin() + l, pos.begin() + r, cut) - pos.begin());
nd_[cur].l = erase(nd_[cur].l, pos, l, mid);
nd_[cur].r = erase(nd_[cur].r, pos, mid, r);
if (nd_[cur].l == -1) return nd_[cur].r;
if (nd_[cur].r == -1) return nd_[cur].l;
pull(cur);
return cur;
}
void collect(int cur, int low, int hi, std::vector<int> &res) {
if (leaf(cur)) {
res.push_back(nd_[cur].low);
} else if (hi <= nd_[cur].bl) {
collect(nd_[cur].l, low, hi, res);
} else if (low >= nd_[cur].br) {
collect(nd_[cur].r, low, hi, res);
} else {
collect(nd_[cur].l, low, nd_[cur].bl, res);
collect(nd_[cur].r, nd_[cur].br, hi, res);
}
}
public:
explicit left_hull(std::vector<point_type> pt, bool rev = false)
: pt_(std::move(pt)), nd_(pt_.size() * 2), rev_(rev) {
build(0, 0, int(pt_.size()));
}
std::vector<int> hull() {
if (rt_ == -1) {
return {};
}
std::vector<int> res;
collect(rt_, 0, int(pt_.size()) - 1, res);
return res;
}
void erase(const std::vector<int> &pos) {
if (!pos.empty()) rt_ = erase(rt_, pos, 0, int(pos.size()));
}
};
} // namespace convex_layers_detail
/// @brief Return the one-based onion layer of every distinct point. Two
/// decremental hull structures maintain the left and right boundary chains;
/// deleting the current boundary from both structures exposes the next layer.
/// Input points must be distinct. Determinants and weighted height comparisons
/// must fit in signed 128 bits.
inline std::vector<int>
convex_layers(const std::vector<point<std::int64_t>> &in) {
using convex_layers_detail::left_hull;
using point_type = convex_layers_detail::point_type;
int n = int(in.size());
if (n == 0) {
return {};
}
std::vector<int> ord(n);
for (int i = 0; i < n; i++) {
ord[i] = i;
}
std::sort(ord.begin(), ord.end(), [&](int lhs, int rhs) {
return std::pair(in[lhs].y, in[lhs].x) < std::pair(in[rhs].y, in[rhs].x);
});
std::vector<point_type> or1(n);
for (int i = 0; i < n; i++) {
or1[i] = in[ord[i]];
}
left_hull l(or1);
std::vector<point_type> rev(or1.rbegin(), or1.rend());
left_hull r(std::move(rev), true);
std::vector<int> lev(n), ans(n);
int del = 0;
for (int le1 = 1; del < n; le1++) {
auto a = l.hull(), b = r.hull();
std::reverse(b.begin(), b.end());
for (int &idx : b) idx = n - 1 - idx;
std::vector<int> bd;
bd.reserve(a.size() + b.size());
std::set_union(a.begin(), a.end(), b.begin(), b.end(), std::back_inserter(bd));
for (int idx : bd) {
lev[idx] = le1;
del++;
}
l.erase(bd);
std::reverse(bd.begin(), bd.end());
for (int &idx : bd) idx = n - 1 - idx;
r.erase(bd);
}
for (int i = 0; i < n; i++) {
ans[ord[i]] = lev[i];
}
return ans;
}
} // namespace noya
// NOYA_CONVEX_LAYERS_HPP
#ifndef NOYA_CONVEX_LAYERS_HPP
#define NOYA_CONVEX_LAYERS_HPP 1
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O(n log^2 n).
/// Space: O(n).
#include "noya/geometry/geometry_base.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iterator>
#include <utility>
#include <vector>
namespace noya {
namespace convex_layers_detail {
using point_type = point<std::int64_t>;
class left_hull {
struct node {
int bl = 0;
int br = 0;
int low = 0;
int hi = 0;
int l = -1;
int r = -1;
};
std::vector<point_type> pt_;
std::vector<node> nd_;
int rt_ = 0;
bool rev_ = false;
static __int128 turn(const point_type &a, const point_type &b,
const point_type &c) {
__int128 x1 = __int128(b.x) - a.x, y1 = __int128(b.y) - a.y;
__int128 x2 = __int128(c.x) - a.x, y2 = __int128(c.y) - a.y;
return x1 * y2 - y1 * x2;
}
bool leaf(int cur) const { return nd_[cur].l == -1 && nd_[cur].r == -1; }
void pull(int cur) {
int l = nd_[cur].l;
int r = nd_[cur].r;
std::int64_t sy = pt_[nd_[r].low].y;
while (!leaf(l) || !leaf(r)) {
int a = nd_[l].bl;
int b = nd_[l].br;
int c = nd_[r].bl;
int d = nd_[r].br;
if (a != b && turn(pt_[a], pt_[b], pt_[c]) > 0) {
l = nd_[l].l;
} else if (c != d && turn(pt_[b], pt_[c], pt_[d]) > 0) {
r = nd_[r].r;
} else if (a == b) {
r = nd_[r].l;
} else if (c == d) {
l = nd_[l].r;
} else {
__int128 lhs = turn(pt_[a], pt_[b], pt_[c]);
__int128 rhs = turn(pt_[b], pt_[a], pt_[d]);
assert(lhs + rhs >= 0);
// Mirroring both coordinates preserves turns and reverses this comparison.
__int128 cmp = lhs * (__int128(pt_[d].y) - sy) +
rhs * (__int128(pt_[c].y) - sy);
if (lhs + rhs == 0 || (rev_ ? cmp > 0 : cmp < 0)) {
l = nd_[l].r;
} else {
r = nd_[r].l;
}
}
}
nd_[cur].bl = nd_[l].low;
nd_[cur].br = nd_[r].low;
}
void build(int cur, int low, int hi) {
nd_[cur].low = low;
nd_[cur].hi = hi;
if (hi - low == 1) {
nd_[cur].bl = nd_[cur].br = low;
nd_[cur].l = nd_[cur].r = -1;
return;
}
int mid = (low + hi) / 2;
nd_[cur].l = cur + 1;
nd_[cur].r = cur + 2 * (mid - low);
build(nd_[cur].l, low, mid);
build(nd_[cur].r, mid, hi);
pull(cur);
}
// Delete a sorted batch, rebuilding each affected bridge only once.
int erase(int cur, const std::vector<int> &pos, int l, int r) {
if (l == r) return cur;
if (r - l == nd_[cur].hi - nd_[cur].low) return -1;
int cut = nd_[nd_[cur].r].low;
int mid = int(std::lower_bound(pos.begin() + l, pos.begin() + r, cut) - pos.begin());
nd_[cur].l = erase(nd_[cur].l, pos, l, mid);
nd_[cur].r = erase(nd_[cur].r, pos, mid, r);
if (nd_[cur].l == -1) return nd_[cur].r;
if (nd_[cur].r == -1) return nd_[cur].l;
pull(cur);
return cur;
}
void collect(int cur, int low, int hi, std::vector<int> &res) {
if (leaf(cur)) {
res.push_back(nd_[cur].low);
} else if (hi <= nd_[cur].bl) {
collect(nd_[cur].l, low, hi, res);
} else if (low >= nd_[cur].br) {
collect(nd_[cur].r, low, hi, res);
} else {
collect(nd_[cur].l, low, nd_[cur].bl, res);
collect(nd_[cur].r, nd_[cur].br, hi, res);
}
}
public:
explicit left_hull(std::vector<point_type> pt, bool rev = false)
: pt_(std::move(pt)), nd_(pt_.size() * 2), rev_(rev) {
build(0, 0, int(pt_.size()));
}
std::vector<int> hull() {
if (rt_ == -1) {
return {};
}
std::vector<int> res;
collect(rt_, 0, int(pt_.size()) - 1, res);
return res;
}
void erase(const std::vector<int> &pos) {
if (!pos.empty()) rt_ = erase(rt_, pos, 0, int(pos.size()));
}
};
} // namespace convex_layers_detail
/// @brief Return the one-based onion layer of every distinct point. Two
/// decremental hull structures maintain the left and right boundary chains;
/// deleting the current boundary from both structures exposes the next layer.
/// Input points must be distinct. Determinants and weighted height comparisons
/// must fit in signed 128 bits.
inline std::vector<int>
convex_layers(const std::vector<point<std::int64_t>> &in) {
using convex_layers_detail::left_hull;
using point_type = convex_layers_detail::point_type;
int n = int(in.size());
if (n == 0) {
return {};
}
std::vector<int> ord(n);
for (int i = 0; i < n; i++) {
ord[i] = i;
}
std::sort(ord.begin(), ord.end(), [&](int lhs, int rhs) {
return std::pair(in[lhs].y, in[lhs].x) < std::pair(in[rhs].y, in[rhs].x);
});
std::vector<point_type> or1(n);
for (int i = 0; i < n; i++) {
or1[i] = in[ord[i]];
}
left_hull l(or1);
std::vector<point_type> rev(or1.rbegin(), or1.rend());
left_hull r(std::move(rev), true);
std::vector<int> lev(n), ans(n);
int del = 0;
for (int le1 = 1; del < n; le1++) {
auto a = l.hull(), b = r.hull();
std::reverse(b.begin(), b.end());
for (int &idx : b) idx = n - 1 - idx;
std::vector<int> bd;
bd.reserve(a.size() + b.size());
std::set_union(a.begin(), a.end(), b.begin(), b.end(), std::back_inserter(bd));
for (int idx : bd) {
lev[idx] = le1;
del++;
}
l.erase(bd);
std::reverse(bd.begin(), bd.end());
for (int &idx : bd) idx = n - 1 - idx;
r.erase(bd);
}
for (int i = 0; i < n; i++) {
ans[ord[i]] = lev[i];
}
return ans;
}
} // namespace noya
#endif // NOYA_CONVEX_LAYERS_HPP
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O(n log^2 n).
/// Space: O(n).
/// @complexity Time: O(1) per point/segment primitive, O(n) for polygon area or
/// point location, O(n log n) for convex hull. Space: O(1), except O(n) for hull.
#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.
/// Starts at the lexicographically smallest point; kc retains boundary points.
/// Coordinate differences and determinant arithmetic must fit in T.
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;
}
if (kc) {
bool col = true;
for (int i = 2; i < int(pt.size()) && col; i++) {
col = orientation(pt[0], pt[1], pt[i]) == 0;
}
if (col) return pt;
}
std::vector<point<T>> res;
res.reserve(pt.size() + 1);
for (const point<T> &p : pt) {
while (res.size() >= 2) {
int rot = orientation(res[res.size() - 2], res.back(), p);
if (rot > 0 || (kc && rot == 0)) {
break;
}
res.pop_back();
}
res.push_back(p);
}
auto k = res.size(); // Keep the lower chain while building the upper chain.
for (int i = int(pt.size()) - 2; i >= 0; i--) {
while (res.size() > k) {
int rot = orientation(res[res.size() - 2], res.back(), pt[i]);
if (rot > 0 || (kc && rot == 0)) {
break;
}
res.pop_back();
}
res.push_back(pt[i]);
}
res.pop_back();
return res;
}
/// @brief Return twice the signed area of a polygon.
/// Coordinate differences, products and the accumulated area must fit in T.
template <class T> T polygon_area2(const std::vector<point<T>> &pg) {
T res{};
if (pg.size() < 3) return res;
auto prev = pg[1] - pg[0];
for (int i = 2; i < int(pg.size()); i++) {
auto cur = pg[i] - pg[0];
res += cross(prev, cur);
prev = cur;
}
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 <cstdint>
#include <iterator>
#include <utility>
#include <vector>
namespace noya {
namespace convex_layers_detail {
using point_type = point<std::int64_t>;
class left_hull {
struct node {
int bl = 0;
int br = 0;
int low = 0;
int hi = 0;
int l = -1;
int r = -1;
};
std::vector<point_type> pt_;
std::vector<node> nd_;
int rt_ = 0;
bool rev_ = false;
static __int128 turn(const point_type &a, const point_type &b,
const point_type &c) {
__int128 x1 = __int128(b.x) - a.x, y1 = __int128(b.y) - a.y;
__int128 x2 = __int128(c.x) - a.x, y2 = __int128(c.y) - a.y;
return x1 * y2 - y1 * x2;
}
bool leaf(int cur) const { return nd_[cur].l == -1 && nd_[cur].r == -1; }
void pull(int cur) {
int l = nd_[cur].l;
int r = nd_[cur].r;
std::int64_t sy = pt_[nd_[r].low].y;
while (!leaf(l) || !leaf(r)) {
int a = nd_[l].bl;
int b = nd_[l].br;
int c = nd_[r].bl;
int d = nd_[r].br;
if (a != b && turn(pt_[a], pt_[b], pt_[c]) > 0) {
l = nd_[l].l;
} else if (c != d && turn(pt_[b], pt_[c], pt_[d]) > 0) {
r = nd_[r].r;
} else if (a == b) {
r = nd_[r].l;
} else if (c == d) {
l = nd_[l].r;
} else {
__int128 lhs = turn(pt_[a], pt_[b], pt_[c]);
__int128 rhs = turn(pt_[b], pt_[a], pt_[d]);
assert(lhs + rhs >= 0);
// Mirroring both coordinates preserves turns and reverses this comparison.
__int128 cmp = lhs * (__int128(pt_[d].y) - sy) +
rhs * (__int128(pt_[c].y) - sy);
if (lhs + rhs == 0 || (rev_ ? cmp > 0 : cmp < 0)) {
l = nd_[l].r;
} else {
r = nd_[r].l;
}
}
}
nd_[cur].bl = nd_[l].low;
nd_[cur].br = nd_[r].low;
}
void build(int cur, int low, int hi) {
nd_[cur].low = low;
nd_[cur].hi = hi;
if (hi - low == 1) {
nd_[cur].bl = nd_[cur].br = low;
nd_[cur].l = nd_[cur].r = -1;
return;
}
int mid = (low + hi) / 2;
nd_[cur].l = cur + 1;
nd_[cur].r = cur + 2 * (mid - low);
build(nd_[cur].l, low, mid);
build(nd_[cur].r, mid, hi);
pull(cur);
}
// Delete a sorted batch, rebuilding each affected bridge only once.
int erase(int cur, const std::vector<int> &pos, int l, int r) {
if (l == r) return cur;
if (r - l == nd_[cur].hi - nd_[cur].low) return -1;
int cut = nd_[nd_[cur].r].low;
int mid = int(std::lower_bound(pos.begin() + l, pos.begin() + r, cut) - pos.begin());
nd_[cur].l = erase(nd_[cur].l, pos, l, mid);
nd_[cur].r = erase(nd_[cur].r, pos, mid, r);
if (nd_[cur].l == -1) return nd_[cur].r;
if (nd_[cur].r == -1) return nd_[cur].l;
pull(cur);
return cur;
}
void collect(int cur, int low, int hi, std::vector<int> &res) {
if (leaf(cur)) {
res.push_back(nd_[cur].low);
} else if (hi <= nd_[cur].bl) {
collect(nd_[cur].l, low, hi, res);
} else if (low >= nd_[cur].br) {
collect(nd_[cur].r, low, hi, res);
} else {
collect(nd_[cur].l, low, nd_[cur].bl, res);
collect(nd_[cur].r, nd_[cur].br, hi, res);
}
}
public:
explicit left_hull(std::vector<point_type> pt, bool rev = false)
: pt_(std::move(pt)), nd_(pt_.size() * 2), rev_(rev) {
build(0, 0, int(pt_.size()));
}
std::vector<int> hull() {
if (rt_ == -1) {
return {};
}
std::vector<int> res;
collect(rt_, 0, int(pt_.size()) - 1, res);
return res;
}
void erase(const std::vector<int> &pos) {
if (!pos.empty()) rt_ = erase(rt_, pos, 0, int(pos.size()));
}
};
} // namespace convex_layers_detail
/// @brief Return the one-based onion layer of every distinct point. Two
/// decremental hull structures maintain the left and right boundary chains;
/// deleting the current boundary from both structures exposes the next layer.
/// Input points must be distinct. Determinants and weighted height comparisons
/// must fit in signed 128 bits.
inline std::vector<int>
convex_layers(const std::vector<point<std::int64_t>> &in) {
using convex_layers_detail::left_hull;
using point_type = convex_layers_detail::point_type;
int n = int(in.size());
if (n == 0) {
return {};
}
std::vector<int> ord(n);
for (int i = 0; i < n; i++) {
ord[i] = i;
}
std::sort(ord.begin(), ord.end(), [&](int lhs, int rhs) {
return std::pair(in[lhs].y, in[lhs].x) < std::pair(in[rhs].y, in[rhs].x);
});
std::vector<point_type> or1(n);
for (int i = 0; i < n; i++) {
or1[i] = in[ord[i]];
}
left_hull l(or1);
std::vector<point_type> rev(or1.rbegin(), or1.rend());
left_hull r(std::move(rev), true);
std::vector<int> lev(n), ans(n);
int del = 0;
for (int le1 = 1; del < n; le1++) {
auto a = l.hull(), b = r.hull();
std::reverse(b.begin(), b.end());
for (int &idx : b) idx = n - 1 - idx;
std::vector<int> bd;
bd.reserve(a.size() + b.size());
std::set_union(a.begin(), a.end(), b.begin(), b.end(), std::back_inserter(bd));
for (int idx : bd) {
lev[idx] = le1;
del++;
}
l.erase(bd);
std::reverse(bd.begin(), bd.end());
for (int &idx : bd) idx = n - 1 - idx;
r.erase(bd);
}
for (int i = 0; i < n; i++) {
ans[ord[i]] = lev[i];
}
return ans;
}
} // namespace noya
// NOYA_CONVEX_LAYERS_HPP