Convex Layers
(geometry/convex_layers.hpp)
- View this file on GitHub
- Last update: 2026-10-05 22:23:07+09:00
- Include:
#include "geometry/convex_layers.hpp"
Overview
convex_layers computes the onion decomposition of a point set. Layer 1
contains every point on the boundary of the original convex hull. After removing
that boundary, layer 2 contains the next convex-hull boundary, and so on until
no points remain.
The implementation maintains both sides of the convex hull under deletions with segment-tree-shaped bridge structures. This is substantially faster than rebuilding a convex hull for every layer when the input has many nested layers.
Function
template <Coordinate T>
std::vector<int> convex_layers(const std::vector<Point<T>>& points);
| Function | Description | Complexity |
|---|---|---|
convex_layers(points) |
Returns the 1-based removal layer of every point in original input order. Does not mutate points. |
$O(N\log^2 N)$ time and $O(N)$ memory. |
All collinear points on a hull edge belong to the same layer. An empty input returns an empty vector. Duplicate coordinates are allowed and receive the same layer.
Integral calculations use wide_type<T>, which is signed 128-bit arithmetic.
Coordinate differences, cross products, and the bridge-comparison intermediate
products must fit that type. Floating-point coordinates use the geometry
module’s default orientation tolerance.
Example
#include "geometry/convex_layers.hpp"
#include <iostream>
#include <vector>
int main() {
using Point = m1une::geometry::Point<long long>;
std::vector<Point> points;
points.emplace_back(0, 0);
points.emplace_back(4, 0);
points.emplace_back(4, 4);
points.emplace_back(0, 4);
points.emplace_back(2, 2);
std::vector<int> layer = m1une::geometry::convex_layers(points);
for (int value : layer) std::cout << value << "\n";
// 1, 1, 1, 1, 2
}
Depends on
Required by
Verified with
- verify/geometry/centroid.test.cpp
- verify/geometry/convex_layers.test.cpp
- verify/geometry/geometry_algorithms.test.cpp
- verify/geometry/rational.test.cpp
Code
#ifndef M1UNE_GEOMETRY_CONVEX_LAYERS_HPP
#define M1UNE_GEOMETRY_CONVEX_LAYERS_HPP 1
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <utility>
#include <vector>
#include "point.hpp"
namespace m1une {
namespace geometry {
namespace convex_layers_detail {
template <Coordinate T>
struct LayerPoint {
wide_type<T> x;
wide_type<T> y;
};
template <Coordinate T>
wide_type<T> layer_cross(
const LayerPoint<T>& first,
const LayerPoint<T>& second,
const LayerPoint<T>& third
) {
return
(second.x - first.x) * (third.y - first.y) -
(second.y - first.y) * (third.x - first.x);
}
template <Coordinate T>
class DecrementalHull {
private:
struct Node {
int left_bound;
int right_bound;
int bridge_left;
int bridge_right;
int left_child;
int right_child;
};
std::vector<LayerPoint<T>> points;
std::vector<Node> nodes;
int root;
bool is_leaf(int node) const {
return nodes[node].left_child == -1 && nodes[node].right_child == -1;
}
void pull(int node) {
int left = nodes[node].left_child;
int right = nodes[node].right_child;
assert(left != -1 && right != -1);
using Wide = wide_type<T>;
const Wide split_y = points[nodes[right].left_bound].y;
while (!is_leaf(left) || !is_leaf(right)) {
const int a = nodes[left].bridge_left;
const int b = nodes[left].bridge_right;
const int c = nodes[right].bridge_left;
const int d = nodes[right].bridge_right;
if (
a != b &&
sign<T>(layer_cross<T>(points[a], points[b], points[c])) > 0
) {
left = nodes[left].left_child;
} else if (
c != d &&
sign<T>(layer_cross<T>(points[b], points[c], points[d])) > 0
) {
right = nodes[right].right_child;
} else if (a == b) {
right = nodes[right].left_child;
} else if (c == d) {
left = nodes[left].right_child;
} else {
const Wide first =
layer_cross<T>(points[a], points[b], points[c]);
const Wide second =
layer_cross<T>(points[b], points[a], points[d]);
const Wide sum = first + second;
assert(sign<T>(sum) >= 0);
const Wide comparison =
first * points[d].y + second * points[c].y - split_y * sum;
if (sign<T>(sum) == 0 || sign<T>(comparison) < 0) {
left = nodes[left].right_child;
} else {
right = nodes[right].left_child;
}
}
}
nodes[node].bridge_left = nodes[left].left_bound;
nodes[node].bridge_right = nodes[right].left_bound;
}
void build(int node, int left, int right) {
nodes[node].left_bound = left;
nodes[node].right_bound = right;
if (right - left == 1) {
nodes[node].bridge_left = left;
nodes[node].bridge_right = left;
nodes[node].left_child = -1;
nodes[node].right_child = -1;
return;
}
const int middle = (left + right) / 2;
nodes[node].left_child = node + 1;
nodes[node].right_child = node + 2 * (middle - left);
build(nodes[node].left_child, left, middle);
build(nodes[node].right_child, middle, right);
pull(node);
}
int erase(int node, int position) {
if (
position < nodes[node].left_bound ||
nodes[node].right_bound <= position
) {
return node;
}
if (nodes[node].right_bound - nodes[node].left_bound == 1) return -1;
nodes[node].left_child = erase(nodes[node].left_child, position);
nodes[node].right_child = erase(nodes[node].right_child, position);
if (nodes[node].left_child == -1) return nodes[node].right_child;
if (nodes[node].right_child == -1) return nodes[node].left_child;
pull(node);
return node;
}
void collect(
int node,
int left,
int right,
std::vector<int>& result
) const {
if (is_leaf(node)) {
result.push_back(nodes[node].left_bound);
} else if (right <= nodes[node].bridge_left) {
collect(nodes[node].left_child, left, right, result);
} else if (nodes[node].bridge_right <= left) {
collect(nodes[node].right_child, left, right, result);
} else {
assert(
left <= nodes[node].bridge_left &&
nodes[node].bridge_right <= right
);
collect(
nodes[node].left_child,
left,
nodes[node].bridge_left,
result
);
collect(
nodes[node].right_child,
nodes[node].bridge_right,
right,
result
);
}
}
public:
explicit DecrementalHull(std::vector<LayerPoint<T>> ordered_points)
: points(std::move(ordered_points)),
nodes(2 * points.size()),
root(points.empty() ? -1 : 0) {
if (!points.empty()) build(0, 0, int(points.size()));
}
std::vector<int> hull() const {
std::vector<int> result;
if (root != -1) collect(root, 0, int(points.size()) - 1, result);
return result;
}
void erase(int position) {
assert(root != -1);
assert(0 <= position && position < int(points.size()));
root = erase(root, position);
}
};
} // namespace convex_layers_detail
template <Coordinate T>
std::vector<int> convex_layers(const std::vector<Point<T>>& points) {
const int n = int(points.size());
if (n == 0) return {};
struct IndexedPoint {
Point<T> point;
int original_index;
};
std::vector<IndexedPoint> indexed;
indexed.reserve(n);
for (int index = 0; index < n; index++) {
indexed.push_back(IndexedPoint{points[index], index});
}
std::sort(
indexed.begin(),
indexed.end(),
[](const IndexedPoint& first, const IndexedPoint& second) {
if (first.point.y != second.point.y) {
return first.point.y < second.point.y;
}
if (first.point.x != second.point.x) {
return first.point.x < second.point.x;
}
return first.original_index < second.original_index;
}
);
std::vector<Point<T>> ordered;
std::vector<int> position(n);
ordered.reserve(n);
for (const IndexedPoint& item : indexed) {
if (ordered.empty() || !(ordered.back() == item.point)) {
ordered.push_back(item.point);
}
position[item.original_index] = int(ordered.size()) - 1;
}
using LayerPoint = convex_layers_detail::LayerPoint<T>;
using Wide = wide_type<T>;
std::vector<LayerPoint> left_points;
left_points.reserve(ordered.size());
for (const Point<T>& point : ordered) {
left_points.push_back(LayerPoint{Wide(point.x), Wide(point.y)});
}
convex_layers_detail::DecrementalHull<T> left_hull(
std::move(left_points)
);
std::vector<LayerPoint> reversed;
reversed.reserve(ordered.size());
for (auto iterator = ordered.rbegin(); iterator != ordered.rend(); ++iterator) {
reversed.push_back(LayerPoint{-Wide(iterator->x), -Wide(iterator->y)});
}
convex_layers_detail::DecrementalHull<T> right_hull(std::move(reversed));
const int distinct_count = int(ordered.size());
std::vector<int> layer_by_position(distinct_count, 0);
std::vector<int> selected_in_layer(distinct_count, 0);
int remaining = distinct_count;
for (int layer = 1; remaining > 0; layer++) {
std::vector<int> boundary;
auto add_boundary = [&](int ordered_position) {
if (selected_in_layer[ordered_position] == layer) return;
selected_in_layer[ordered_position] = layer;
boundary.push_back(ordered_position);
};
for (int ordered_position : left_hull.hull()) {
add_boundary(ordered_position);
}
for (int reversed_position : right_hull.hull()) {
add_boundary(distinct_count - 1 - reversed_position);
}
assert(!boundary.empty());
for (int ordered_position : boundary) {
layer_by_position[ordered_position] = layer;
left_hull.erase(ordered_position);
right_hull.erase(distinct_count - 1 - ordered_position);
remaining--;
}
}
std::vector<int> result(n);
for (int index = 0; index < n; index++) {
result[index] = layer_by_position[position[index]];
}
return result;
}
} // namespace geometry
} // namespace m1une
#endif // M1UNE_GEOMETRY_CONVEX_LAYERS_HPP#line 1 "geometry/convex_layers.hpp"
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <utility>
#include <vector>
#line 1 "geometry/point.hpp"
#include <cmath>
#include <concepts>
#line 7 "geometry/point.hpp"
#include <type_traits>
#line 1 "geometry/detail/floating_predicate.hpp"
namespace m1une {
namespace geometry {
namespace predicate_detail {
template <typename T>
constexpr T absolute(T value) {
return value < T(0) ? -value : value;
}
template <typename T>
constexpr T max_value(T first, T second) {
return first < second ? second : first;
}
template <typename T>
constexpr T vector_scale(T x, T y) {
return max_value(absolute(x), absolute(y));
}
template <bool Exact, typename T>
constexpr int scaled_sign(T value, T scale, long double eps) {
if constexpr (Exact) {
return (value > T(0)) - (value < T(0));
} else {
const T tolerance = T(eps) * scale;
return (value > tolerance) - (value < -tolerance);
}
}
template <bool Exact, typename T>
constexpr T determinant_scale(T ax, T ay, T bx, T by) {
if constexpr (Exact) {
return T(0);
} else {
return vector_scale(ax, ay) * vector_scale(bx, by);
}
}
template <bool Exact, typename T>
constexpr int determinant_sign(
T ax,
T ay,
T bx,
T by,
long double eps
) {
const T determinant = ax * by - ay * bx;
return scaled_sign<Exact>(
determinant,
determinant_scale<Exact>(ax, ay, bx, by),
eps
);
}
template <bool Exact, typename T>
constexpr int orientation_sign(
T direction_x,
T direction_y,
T offset_x,
T offset_y,
long double eps
) {
const T determinant =
direction_x * offset_y - direction_y * offset_x;
T scale = T(0);
if constexpr (!Exact) {
const T direction_scale =
vector_scale(direction_x, direction_y);
scale = direction_scale * max_value(
direction_scale,
vector_scale(offset_x, offset_y)
);
}
return scaled_sign<Exact>(determinant, scale, eps);
}
template <bool Exact, typename T>
constexpr int dot_sign(
T ax,
T ay,
T bx,
T by,
long double eps
) {
const T value = ax * bx + ay * by;
T scale = T(0);
if constexpr (!Exact) {
scale = vector_scale(ax, ay) * vector_scale(bx, by);
}
return scaled_sign<Exact>(value, scale, eps);
}
} // namespace predicate_detail
} // namespace geometry
} // namespace m1une
#line 10 "geometry/point.hpp"
namespace m1une {
namespace geometry {
template <typename T>
concept Coordinate = !std::same_as<std::remove_cv_t<T>, bool> &&
(std::is_arithmetic_v<T> ||
(std::copyable<T> && std::totally_ordered<T> && requires(T a, T b) {
T(0);
T(1);
static_cast<long double>(a);
{ +a } -> std::same_as<T>;
{ -a } -> std::same_as<T>;
{ a + b } -> std::same_as<T>;
{ a - b } -> std::same_as<T>;
{ a * b } -> std::same_as<T>;
{ a / b } -> std::same_as<T>;
{ a += b } -> std::same_as<T&>;
{ a -= b } -> std::same_as<T&>;
}));
// Custom coordinate types keep their own exact arithmetic.
template <typename T>
concept ExactCoordinate = Coordinate<T> && !std::floating_point<T>;
template <Coordinate T>
using wide_type = std::conditional_t<std::integral<T>, __int128_t,
std::conditional_t<std::floating_point<T>, long double, T>>;
template <Coordinate T>
struct Point {
T x;
T y;
constexpr Point() : x(0), y(0) {}
constexpr Point(T x_value, T y_value) : x(x_value), y(y_value) {}
template <Coordinate U>
explicit constexpr Point(const Point<U>& other)
: x(static_cast<T>(other.x)), y(static_cast<T>(other.y)) {}
constexpr Point& operator+=(const Point& other) {
x += other.x;
y += other.y;
return *this;
}
constexpr Point& operator-=(const Point& other) {
x -= other.x;
y -= other.y;
return *this;
}
constexpr Point operator+() const {
return *this;
}
constexpr Point operator-() const {
return Point(-x, -y);
}
friend constexpr Point operator+(Point left, const Point& right) {
return left += right;
}
friend constexpr Point operator-(Point left, const Point& right) {
return left -= right;
}
friend constexpr bool operator==(const Point&, const Point&) = default;
friend constexpr bool operator<(const Point& left, const Point& right) {
if (left.x != right.x) return left.x < right.x;
return left.y < right.y;
}
};
template <Coordinate T>
constexpr Point<long double> centroid(const Point<T>& point) {
return Point<long double>(point);
}
template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(const Point<T>& point, Scalar scalar) {
using Result = std::common_type_t<T, Scalar>;
return Point<Result>(
Result(point.x) * Result(scalar),
Result(point.y) * Result(scalar)
);
}
template <typename Scalar, Coordinate T>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(Scalar scalar, const Point<T>& point) {
return point * scalar;
}
template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator/(const Point<T>& point, Scalar scalar) {
using Result = std::common_type_t<T, Scalar>;
return Point<Result>(
Result(point.x) / Result(scalar),
Result(point.y) / Result(scalar)
);
}
template <Coordinate T>
constexpr wide_type<T> dot(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
return W(a.x) * W(b.x) + W(a.y) * W(b.y);
}
template <Coordinate T>
constexpr wide_type<T> cross(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
return W(a.x) * W(b.y) - W(a.y) * W(b.x);
}
template <Coordinate T>
constexpr wide_type<T> cross(
const Point<T>& origin,
const Point<T>& a,
const Point<T>& b
) {
using W = wide_type<T>;
W ax = W(a.x) - W(origin.x);
W ay = W(a.y) - W(origin.y);
W bx = W(b.x) - W(origin.x);
W by = W(b.y) - W(origin.y);
return ax * by - ay * bx;
}
template <Coordinate T>
constexpr wide_type<T> norm2(const Point<T>& point) {
return dot(point, point);
}
template <Coordinate T>
constexpr wide_type<T> distance2(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
W dx = W(a.x) - W(b.x);
W dy = W(a.y) - W(b.y);
return dx * dx + dy * dy;
}
template <Coordinate T>
long double norm(const Point<T>& point) {
return std::hypot(
static_cast<long double>(point.x),
static_cast<long double>(point.y)
);
}
template <Coordinate 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)
);
}
template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
(std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> internal_division_point(
const Point<T>& a,
const Point<T>& b,
M m,
N n
) {
long double first_ratio = static_cast<long double>(m);
long double second_ratio = static_cast<long double>(n);
long double denominator = first_ratio + second_ratio;
assert(denominator != 0);
Point<long double> first(a);
Point<long double> direction = Point<long double>(b) - first;
return first + direction * (first_ratio / denominator);
}
template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
(std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> external_division_point(
const Point<T>& a,
const Point<T>& b,
M m,
N n
) {
long double first_ratio = static_cast<long double>(m);
long double second_ratio = static_cast<long double>(n);
long double denominator = first_ratio - second_ratio;
assert(denominator != 0);
Point<long double> first(a);
Point<long double> direction = Point<long double>(b) - first;
return first + direction * (first_ratio / denominator);
}
template <Coordinate T>
constexpr int sign(wide_type<T> value, long double eps = 1e-12L) {
return predicate_detail::scaled_sign<ExactCoordinate<T>>(
value,
wide_type<T>(1),
eps
);
}
template <Coordinate T>
constexpr int orientation(
const Point<T>& a,
const Point<T>& b,
const Point<T>& c,
long double eps = 1e-12L
) {
using W = wide_type<T>;
const W first_x = W(b.x) - W(a.x);
const W first_y = W(b.y) - W(a.y);
const W second_x = W(c.x) - W(a.x);
const W second_y = W(c.y) - W(a.y);
return predicate_detail::orientation_sign<ExactCoordinate<T>>(
first_x,
first_y,
second_x,
second_y,
eps
);
}
template <Coordinate T>
constexpr bool collinear(
const Point<T>& a,
const Point<T>& b,
const Point<T>& c,
long double eps = 1e-12L
) {
return orientation(a, b, c, eps) == 0;
}
template <Coordinate T>
Point<long double> rotate(const Point<T>& point, long double angle) {
long double cosine = std::cos(angle);
long double sine = std::sin(angle);
return Point<long double>(
static_cast<long double>(point.x) * cosine -
static_cast<long double>(point.y) * sine,
static_cast<long double>(point.x) * sine +
static_cast<long double>(point.y) * cosine
);
}
template <Coordinate T>
Point<long double> normalized(const Point<T>& point) {
long double length = norm(point);
assert(length != 0);
return Point<long double>(
static_cast<long double>(point.x) / length,
static_cast<long double>(point.y) / length
);
}
} // namespace geometry
} // namespace m1une
#line 11 "geometry/convex_layers.hpp"
namespace m1une {
namespace geometry {
namespace convex_layers_detail {
template <Coordinate T>
struct LayerPoint {
wide_type<T> x;
wide_type<T> y;
};
template <Coordinate T>
wide_type<T> layer_cross(
const LayerPoint<T>& first,
const LayerPoint<T>& second,
const LayerPoint<T>& third
) {
return
(second.x - first.x) * (third.y - first.y) -
(second.y - first.y) * (third.x - first.x);
}
template <Coordinate T>
class DecrementalHull {
private:
struct Node {
int left_bound;
int right_bound;
int bridge_left;
int bridge_right;
int left_child;
int right_child;
};
std::vector<LayerPoint<T>> points;
std::vector<Node> nodes;
int root;
bool is_leaf(int node) const {
return nodes[node].left_child == -1 && nodes[node].right_child == -1;
}
void pull(int node) {
int left = nodes[node].left_child;
int right = nodes[node].right_child;
assert(left != -1 && right != -1);
using Wide = wide_type<T>;
const Wide split_y = points[nodes[right].left_bound].y;
while (!is_leaf(left) || !is_leaf(right)) {
const int a = nodes[left].bridge_left;
const int b = nodes[left].bridge_right;
const int c = nodes[right].bridge_left;
const int d = nodes[right].bridge_right;
if (
a != b &&
sign<T>(layer_cross<T>(points[a], points[b], points[c])) > 0
) {
left = nodes[left].left_child;
} else if (
c != d &&
sign<T>(layer_cross<T>(points[b], points[c], points[d])) > 0
) {
right = nodes[right].right_child;
} else if (a == b) {
right = nodes[right].left_child;
} else if (c == d) {
left = nodes[left].right_child;
} else {
const Wide first =
layer_cross<T>(points[a], points[b], points[c]);
const Wide second =
layer_cross<T>(points[b], points[a], points[d]);
const Wide sum = first + second;
assert(sign<T>(sum) >= 0);
const Wide comparison =
first * points[d].y + second * points[c].y - split_y * sum;
if (sign<T>(sum) == 0 || sign<T>(comparison) < 0) {
left = nodes[left].right_child;
} else {
right = nodes[right].left_child;
}
}
}
nodes[node].bridge_left = nodes[left].left_bound;
nodes[node].bridge_right = nodes[right].left_bound;
}
void build(int node, int left, int right) {
nodes[node].left_bound = left;
nodes[node].right_bound = right;
if (right - left == 1) {
nodes[node].bridge_left = left;
nodes[node].bridge_right = left;
nodes[node].left_child = -1;
nodes[node].right_child = -1;
return;
}
const int middle = (left + right) / 2;
nodes[node].left_child = node + 1;
nodes[node].right_child = node + 2 * (middle - left);
build(nodes[node].left_child, left, middle);
build(nodes[node].right_child, middle, right);
pull(node);
}
int erase(int node, int position) {
if (
position < nodes[node].left_bound ||
nodes[node].right_bound <= position
) {
return node;
}
if (nodes[node].right_bound - nodes[node].left_bound == 1) return -1;
nodes[node].left_child = erase(nodes[node].left_child, position);
nodes[node].right_child = erase(nodes[node].right_child, position);
if (nodes[node].left_child == -1) return nodes[node].right_child;
if (nodes[node].right_child == -1) return nodes[node].left_child;
pull(node);
return node;
}
void collect(
int node,
int left,
int right,
std::vector<int>& result
) const {
if (is_leaf(node)) {
result.push_back(nodes[node].left_bound);
} else if (right <= nodes[node].bridge_left) {
collect(nodes[node].left_child, left, right, result);
} else if (nodes[node].bridge_right <= left) {
collect(nodes[node].right_child, left, right, result);
} else {
assert(
left <= nodes[node].bridge_left &&
nodes[node].bridge_right <= right
);
collect(
nodes[node].left_child,
left,
nodes[node].bridge_left,
result
);
collect(
nodes[node].right_child,
nodes[node].bridge_right,
right,
result
);
}
}
public:
explicit DecrementalHull(std::vector<LayerPoint<T>> ordered_points)
: points(std::move(ordered_points)),
nodes(2 * points.size()),
root(points.empty() ? -1 : 0) {
if (!points.empty()) build(0, 0, int(points.size()));
}
std::vector<int> hull() const {
std::vector<int> result;
if (root != -1) collect(root, 0, int(points.size()) - 1, result);
return result;
}
void erase(int position) {
assert(root != -1);
assert(0 <= position && position < int(points.size()));
root = erase(root, position);
}
};
} // namespace convex_layers_detail
template <Coordinate T>
std::vector<int> convex_layers(const std::vector<Point<T>>& points) {
const int n = int(points.size());
if (n == 0) return {};
struct IndexedPoint {
Point<T> point;
int original_index;
};
std::vector<IndexedPoint> indexed;
indexed.reserve(n);
for (int index = 0; index < n; index++) {
indexed.push_back(IndexedPoint{points[index], index});
}
std::sort(
indexed.begin(),
indexed.end(),
[](const IndexedPoint& first, const IndexedPoint& second) {
if (first.point.y != second.point.y) {
return first.point.y < second.point.y;
}
if (first.point.x != second.point.x) {
return first.point.x < second.point.x;
}
return first.original_index < second.original_index;
}
);
std::vector<Point<T>> ordered;
std::vector<int> position(n);
ordered.reserve(n);
for (const IndexedPoint& item : indexed) {
if (ordered.empty() || !(ordered.back() == item.point)) {
ordered.push_back(item.point);
}
position[item.original_index] = int(ordered.size()) - 1;
}
using LayerPoint = convex_layers_detail::LayerPoint<T>;
using Wide = wide_type<T>;
std::vector<LayerPoint> left_points;
left_points.reserve(ordered.size());
for (const Point<T>& point : ordered) {
left_points.push_back(LayerPoint{Wide(point.x), Wide(point.y)});
}
convex_layers_detail::DecrementalHull<T> left_hull(
std::move(left_points)
);
std::vector<LayerPoint> reversed;
reversed.reserve(ordered.size());
for (auto iterator = ordered.rbegin(); iterator != ordered.rend(); ++iterator) {
reversed.push_back(LayerPoint{-Wide(iterator->x), -Wide(iterator->y)});
}
convex_layers_detail::DecrementalHull<T> right_hull(std::move(reversed));
const int distinct_count = int(ordered.size());
std::vector<int> layer_by_position(distinct_count, 0);
std::vector<int> selected_in_layer(distinct_count, 0);
int remaining = distinct_count;
for (int layer = 1; remaining > 0; layer++) {
std::vector<int> boundary;
auto add_boundary = [&](int ordered_position) {
if (selected_in_layer[ordered_position] == layer) return;
selected_in_layer[ordered_position] = layer;
boundary.push_back(ordered_position);
};
for (int ordered_position : left_hull.hull()) {
add_boundary(ordered_position);
}
for (int reversed_position : right_hull.hull()) {
add_boundary(distinct_count - 1 - reversed_position);
}
assert(!boundary.empty());
for (int ordered_position : boundary) {
layer_by_position[ordered_position] = layer;
left_hull.erase(ordered_position);
right_hull.erase(distinct_count - 1 - ordered_position);
remaining--;
}
}
std::vector<int> result(n);
for (int index = 0; index < n; index++) {
result[index] = layer_by_position[position[index]];
}
return result;
}
} // namespace geometry
} // namespace m1une