convex_layers.hpp¶
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.
Verified by convex_layers.
逐层剥离凸包，返回每个不同点所属的洋葱层编号；用于凸包分层和点集深度问题。
Implementation¶
#ifndef NOYA_CONVEX_LAYERS_HPP
#define NOYA_CONVEX_LAYERS_HPP 1
/// @complexity Time: O(n log^2 n).
/// Space: O(n).
#include "noya/geometry_base.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <set>
#include <utility>
#include <vector>
namespace noya {
namespace convex_layers_detail {
using point_type = point<std::int64_t>;
class left_hull {
struct node {
int bridge_left = 0;
int bridge_right = 0;
int low = 0;
int high = 0;
int left = -1;
int right = -1;
};
std::vector<point_type> points_;
std::vector<node> nodes_;
int root_ = 0;
bool leaf(int current) const {
return nodes_[current].left == -1 && nodes_[current].right == -1;
}
void pull(int current) {
int left = nodes_[current].left;
int right = nodes_[current].right;
std::int64_t split_y = points_[nodes_[right].low].y;
while (!leaf(left) || !leaf(right)) {
int a = nodes_[left].bridge_left;
int b = nodes_[left].bridge_right;
int c = nodes_[right].bridge_left;
int d = nodes_[right].bridge_right;
if (a != b && cross(points_[a], points_[b], points_[c]) > 0) {
left = nodes_[left].left;
} else if (c != d && cross(points_[b], points_[c], points_[d]) > 0) {
right = nodes_[right].right;
} else if (a == b) {
right = nodes_[right].left;
} else if (c == d) {
left = nodes_[left].right;
} else {
std::int64_t first = cross(points_[a], points_[b], points_[c]);
std::int64_t second = cross(points_[b], points_[a], points_[d]);
assert(first + second >= 0);
if (first + second == 0 ||
first * points_[d].y + second * points_[c].y <
split_y * (first + second)) {
left = nodes_[left].right;
} else {
right = nodes_[right].left;
}
}
}
nodes_[current].bridge_left = nodes_[left].low;
nodes_[current].bridge_right = nodes_[right].low;
}
void build(int current, int low, int high) {
nodes_[current].low = low;
nodes_[current].high = high;
if (high - low == 1) {
nodes_[current].bridge_left = nodes_[current].bridge_right = low;
nodes_[current].left = nodes_[current].right = -1;
return;
}
int middle = (low + high) / 2;
nodes_[current].left = current + 1;
nodes_[current].right = current + 2 * (middle - low);
build(nodes_[current].left, low, middle);
build(nodes_[current].right, middle, high);
pull(current);
}
int erase(int current, int low, int high) {
if (current == -1 || high <= nodes_[current].low ||
nodes_[current].high <= low) {
return current;
}
if (low <= nodes_[current].low && nodes_[current].high <= high) {
return -1;
}
nodes_[current].left = erase(nodes_[current].left, low, high);
nodes_[current].right = erase(nodes_[current].right, low, high);
if (nodes_[current].left == -1) {
return nodes_[current].right;
}
if (nodes_[current].right == -1) {
return nodes_[current].left;
}
pull(current);
return current;
}
void collect(int current, int low, int high, std::vector<int> &result) {
if (leaf(current)) {
result.push_back(nodes_[current].low);
} else if (high <= nodes_[current].bridge_left) {
collect(nodes_[current].left, low, high, result);
} else if (low >= nodes_[current].bridge_right) {
collect(nodes_[current].right, low, high, result);
} else {
collect(nodes_[current].left, low, nodes_[current].bridge_left, result);
collect(nodes_[current].right, nodes_[current].bridge_right, high, result);
}
}
public:
explicit left_hull(std::vector<point_type> points)
: points_(std::move(points)), nodes_(points_.size() * 2) {
build(0, 0, int(points_.size()));
}
std::vector<int> hull() {
if (root_ == -1) {
return {};
}
std::vector<int> result;
collect(root_, 0, int(points_.size()) - 1, result);
return result;
}
void erase(int position) { root_ = erase(root_, position, position + 1); }
};
} // 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.
inline std::vector<int>
convex_layers(const std::vector<point<std::int64_t>> &input) {
using convex_layers_detail::left_hull;
using point_type = convex_layers_detail::point_type;
int n = int(input.size());
if (n == 0) {
return {};
}
std::vector<int> order(n);
for (int i = 0; i < n; i++) {
order[i] = i;
}
std::sort(order.begin(), order.end(), [&](int first, int second) {
return std::pair(input[first].y, input[first].x) <
std::pair(input[second].y, input[second].x);
});
std::vector<point_type> sorted(n);
for (int i = 0; i < n; i++) {
sorted[i] = input[order[i]];
}
left_hull left(sorted);
std::vector<point_type> reflected(sorted.rbegin(), sorted.rend());
for (auto &value : reflected) {
value.x = -value.x;
value.y = -value.y;
}
left_hull right(std::move(reflected));
std::vector<int> sorted_layer(n), answer(n);
int removed = 0;
for (int layer = 1; removed < n; layer++) {
std::set<int> boundary;
for (int index : left.hull()) {
boundary.insert(index);
}
for (int index : right.hull()) {
boundary.insert(n - 1 - index);
}
for (int index : boundary) {
sorted_layer[index] = layer;
removed++;
left.erase(index);
right.erase(n - 1 - index);
}
}
for (int i = 0; i < n; i++) {
answer[order[i]] = sorted_layer[i];
}
return answer;
}
} // namespace noya
#endif // NOYA_CONVEX_LAYERS_HPP
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <optional>
#include <set>
#include <utility>
#include <vector>
/// @complexity Time: O(n log^2 n).
/// Space: O(n).
/// @complexity Time: O(1) per primitive; O(n log n) for convex hull.
/// Space: O(1) per primitive and O(n) for hull construction.
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 &other) {
x += other.x;
y += other.y;
return *this;
}
point &operator-=(const point &other) {
x -= other.x;
y -= other.y;
return *this;
}
point &operator*=(const T &scale) {
x *= scale;
y *= scale;
return *this;
}
point &operator/=(const T &scale) {
x /= scale;
y /= scale;
return *this;
}
friend point operator+(point left, const point &right) {
return left += right;
}
friend point operator-(point left, const point &right) {
return left -= right;
}
friend point operator*(point value, const T &scale) { return value *= scale; }
friend point operator*(const T &scale, point value) { return value *= scale; }
friend point operator/(point value, const T &scale) { return value /= scale; }
friend bool operator==(const point &, const point &) = default;
friend bool operator<(const point &left, const point &right) {
return left.x < right.x || (left.x == right.x && left.y < right.y);
}
};
/// @brief Circle represented by a center and a nonnegative radius.
template <class Real> struct circle {
point<Real> center;
Real radius{};
};
/// @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 - origin, b - origin).
template <class T>
T cross(const point<T> &origin, const point<T> &a, const point<T> &b) {
return cross(a - origin, b - origin);
}
/// @brief Return the squared Euclidean norm.
template <class T> T norm2(const point<T> &value) { return dot(value, value); }
/// @brief Compare a value with zero using an optional absolute tolerance.
template <class T> int sign(const T &value, const T &epsilon = T{}) {
return (value > epsilon) - (value < -epsilon);
}
/// @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 &epsilon = T{}) {
return sign(cross(a, b, c), epsilon);
}
/// @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 &epsilon = T{}) {
if (orientation(a, b, p, epsilon) != 0) {
return false;
}
return std::min(a.x, b.x) - epsilon <= p.x &&
p.x <= std::max(a.x, b.x) + epsilon &&
std::min(a.y, b.y) - epsilon <= p.y &&
p.y <= std::max(a.y, b.y) + epsilon;
}
/// @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 &epsilon = T{}) {
int ab_c = orientation(a, b, c, epsilon);
int ab_d = orientation(a, b, d, epsilon);
int cd_a = orientation(c, d, a, epsilon);
int cd_b = orientation(c, d, b, epsilon);
if (ab_c == 0 && on_segment(c, a, b, epsilon)) {
return true;
}
if (ab_d == 0 && on_segment(d, a, b, epsilon)) {
return true;
}
if (cd_a == 0 && on_segment(a, c, d, epsilon)) {
return true;
}
if (cd_b == 0 && on_segment(b, c, d, epsilon)) {
return true;
}
return ab_c * ab_d < 0 && cd_a * cd_b < 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 epsilon = 0) {
point<long double> first{static_cast<long double>(a.x),
static_cast<long double>(a.y)};
point<long double> second{static_cast<long double>(b.x),
static_cast<long double>(b.y)};
point<long double> third{static_cast<long double>(c.x),
static_cast<long double>(c.y)};
point<long double> fourth{static_cast<long double>(d.x),
static_cast<long double>(d.y)};
point<long double> direction_a = second - first;
point<long double> direction_b = fourth - third;
long double denominator = cross(direction_a, direction_b);
if (std::abs(denominator) <= epsilon) {
return std::nullopt;
}
long double ratio = cross(third - first, direction_b) / denominator;
return first + direction_a * ratio;
}
/// @brief Return the convex hull in counter-clockwise order without repetition.
template <class T>
std::vector<point<T>> convex_hull(std::vector<point<T>> points,
bool keep_collinear = false) {
std::sort(points.begin(), points.end());
points.erase(std::unique(points.begin(), points.end()), points.end());
if (points.size() <= 1) {
return points;
}
bool all_collinear = true;
for (int i = 2; i < int(points.size()); i++) {
all_collinear &= orientation(points[0], points[1], points[i]) == 0;
}
if (keep_collinear && all_collinear) {
return points;
}
std::vector<point<T>> lower, upper;
for (const point<T> &p : points) {
while (lower.size() >= 2) {
int turn = orientation(lower[lower.size() - 2], lower.back(), p);
if (turn > 0 || (keep_collinear && turn == 0)) {
break;
}
lower.pop_back();
}
lower.push_back(p);
}
for (auto it = points.rbegin(); it != points.rend(); ++it) {
while (upper.size() >= 2) {
int turn = orientation(upper[upper.size() - 2], upper.back(), *it);
if (turn > 0 || (keep_collinear && turn == 0)) {
break;
}
upper.pop_back();
}
upper.push_back(*it);
}
lower.pop_back();
upper.pop_back();
lower.insert(lower.end(), upper.begin(), upper.end());
return lower;
}
/// @brief Return twice the signed area of a polygon.
template <class T> T polygon_area2(const std::vector<point<T>> &polygon) {
T result{};
for (int i = 0; i < int(polygon.size()); i++) {
result += cross(polygon[i], polygon[(i + 1) % polygon.size()]);
}
return result;
}
/// @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>> &polygon) {
bool inside = false;
for (int i = 0; i < int(polygon.size()); i++) {
point<T> a = polygon[i];
point<T> b = polygon[(i + 1) % polygon.size()];
if (on_segment(p, a, b)) {
return 0;
}
if (a.y <= p.y && p.y < b.y && orientation(a, b, p) > 0) {
inside = !inside;
}
if (b.y <= p.y && p.y < a.y && orientation(a, b, p) < 0) {
inside = !inside;
}
}
return inside ? 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
namespace noya {
namespace convex_layers_detail {
using point_type = point<std::int64_t>;
class left_hull {
struct node {
int bridge_left = 0;
int bridge_right = 0;
int low = 0;
int high = 0;
int left = -1;
int right = -1;
};
std::vector<point_type> points_;
std::vector<node> nodes_;
int root_ = 0;
bool leaf(int current) const {
return nodes_[current].left == -1 && nodes_[current].right == -1;
}
void pull(int current) {
int left = nodes_[current].left;
int right = nodes_[current].right;
std::int64_t split_y = points_[nodes_[right].low].y;
while (!leaf(left) || !leaf(right)) {
int a = nodes_[left].bridge_left;
int b = nodes_[left].bridge_right;
int c = nodes_[right].bridge_left;
int d = nodes_[right].bridge_right;
if (a != b && cross(points_[a], points_[b], points_[c]) > 0) {
left = nodes_[left].left;
} else if (c != d && cross(points_[b], points_[c], points_[d]) > 0) {
right = nodes_[right].right;
} else if (a == b) {
right = nodes_[right].left;
} else if (c == d) {
left = nodes_[left].right;
} else {
std::int64_t first = cross(points_[a], points_[b], points_[c]);
std::int64_t second = cross(points_[b], points_[a], points_[d]);
assert(first + second >= 0);
if (first + second == 0 ||
first * points_[d].y + second * points_[c].y <
split_y * (first + second)) {
left = nodes_[left].right;
} else {
right = nodes_[right].left;
}
}
}
nodes_[current].bridge_left = nodes_[left].low;
nodes_[current].bridge_right = nodes_[right].low;
}
void build(int current, int low, int high) {
nodes_[current].low = low;
nodes_[current].high = high;
if (high - low == 1) {
nodes_[current].bridge_left = nodes_[current].bridge_right = low;
nodes_[current].left = nodes_[current].right = -1;
return;
}
int middle = (low + high) / 2;
nodes_[current].left = current + 1;
nodes_[current].right = current + 2 * (middle - low);
build(nodes_[current].left, low, middle);
build(nodes_[current].right, middle, high);
pull(current);
}
int erase(int current, int low, int high) {
if (current == -1 || high <= nodes_[current].low ||
nodes_[current].high <= low) {
return current;
}
if (low <= nodes_[current].low && nodes_[current].high <= high) {
return -1;
}
nodes_[current].left = erase(nodes_[current].left, low, high);
nodes_[current].right = erase(nodes_[current].right, low, high);
if (nodes_[current].left == -1) {
return nodes_[current].right;
}
if (nodes_[current].right == -1) {
return nodes_[current].left;
}
pull(current);
return current;
}
void collect(int current, int low, int high, std::vector<int> &result) {
if (leaf(current)) {
result.push_back(nodes_[current].low);
} else if (high <= nodes_[current].bridge_left) {
collect(nodes_[current].left, low, high, result);
} else if (low >= nodes_[current].bridge_right) {
collect(nodes_[current].right, low, high, result);
} else {
collect(nodes_[current].left, low, nodes_[current].bridge_left, result);
collect(nodes_[current].right, nodes_[current].bridge_right, high, result);
}
}
public:
explicit left_hull(std::vector<point_type> points)
: points_(std::move(points)), nodes_(points_.size() * 2) {
build(0, 0, int(points_.size()));
}
std::vector<int> hull() {
if (root_ == -1) {
return {};
}
std::vector<int> result;
collect(root_, 0, int(points_.size()) - 1, result);
return result;
}
void erase(int position) { root_ = erase(root_, position, position + 1); }
};
} // 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.
inline std::vector<int>
convex_layers(const std::vector<point<std::int64_t>> &input) {
using convex_layers_detail::left_hull;
using point_type = convex_layers_detail::point_type;
int n = int(input.size());
if (n == 0) {
return {};
}
std::vector<int> order(n);
for (int i = 0; i < n; i++) {
order[i] = i;
}
std::sort(order.begin(), order.end(), [&](int first, int second) {
return std::pair(input[first].y, input[first].x) <
std::pair(input[second].y, input[second].x);
});
std::vector<point_type> sorted(n);
for (int i = 0; i < n; i++) {
sorted[i] = input[order[i]];
}
left_hull left(sorted);
std::vector<point_type> reflected(sorted.rbegin(), sorted.rend());
for (auto &value : reflected) {
value.x = -value.x;
value.y = -value.y;
}
left_hull right(std::move(reflected));
std::vector<int> sorted_layer(n), answer(n);
int removed = 0;
for (int layer = 1; removed < n; layer++) {
std::set<int> boundary;
for (int index : left.hull()) {
boundary.insert(index);
}
for (int index : right.hull()) {
boundary.insert(n - 1 - index);
}
for (int index : boundary) {
sorted_layer[index] = layer;
removed++;
left.erase(index);
right.erase(n - 1 - index);
}
}
for (int i = 0; i < n; i++) {
answer[order[i]] = sorted_layer[i];
}
return answer;
}
} // namespace noya