kd_tree.hpp¶
维护低维点集的空间划分，用于矩形范围搜索或最近邻类剪枝；适合维数很小的几何查询。
Complexity: Time: Expected \(O(n \log n)\) build; query time is output/pruning dependent, worst \(O(n)\). Space: \(O(n)\).
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: Expected O(n log n) build; query time is output/pruning dependent, worst O(n).
/// Space: O(n).
#include "noya/geometry/geometry_base.hpp"
#include <algorithm>
#include <cassert>
#include <numeric>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
/// @brief Static two-dimensional KD-tree for rectangle reporting and nearest
/// neighbor queries.
template <class T, class Distance = std::conditional_t<std::is_integral_v<T>,
__int128_t, long double>>
struct kd_tree {
using point_type = point<T>;
using distance_type = Distance;
struct node {
int pid = -1;
int ls = -1;
int rs = -1;
T xl{};
T xr{};
T yl{};
T yr{};
int min_pid = -1;
};
std::vector<point_type> pt;
std::vector<node> tr;
int rt = -1;
kd_tree() = default;
explicit kd_tree(const std::vector<point_type> &a) { build(a); }
/// @brief Rebuild the balanced tree in O(n log n) expected time.
void build(const std::vector<point_type> &a) {
pt = a;
tr.clear();
rt = -1;
std::vector<int> ord(pt.size());
std::iota(ord.begin(), ord.end(), 0);
if (!ord.empty()) {
tr.reserve(ord.size());
rt = build_range(ord, 0, int(ord.size()), 0);
}
}
/// @brief Return the number of stored points.
int size() const { return int(pt.size()); }
/// @brief Return point indices in [l, r) x [dn, top), stopping
/// after lim results when lim is nonnegative.
std::vector<int> range_query(const T &l, const T &r, const T &dn,
const T &top, int lim = -1) const {
assert(!(r < l) && !(top < dn));
int li0 = lim < 0 ? size() : lim;
std::vector<int> res;
if (li0 == 0 || l == r || dn == top)
return res;
res.reserve(std::min(size(), li0));
range_query_at(rt, int(tr.size()), l, r, dn, top, li0, res);
return res;
}
/// @brief Return (point index, squared distance) for a nearest point; the
/// smallest index wins ties, and -1 denotes an empty tree. Squared distances
/// must be representable in distance_type.
std::pair<int, distance_type> nearest(const point_type &q) const {
int bid = -1;
distance_type bd{};
nearest_at(rt, q, bid, bd);
return {bid, bd};
}
private:
static distance_type square_difference(const T &arr, const T &b) {
distance_type dif =
static_cast<distance_type>(arr) - static_cast<distance_type>(b);
return dif * dif;
}
distance_type point_distance(int idx, const point_type &q) const {
return square_difference(pt[idx].x, q.x) +
square_difference(pt[idx].y, q.y);
}
distance_type box_distance(const node &cur, const point_type &q) const {
distance_type res{};
if (q.x < cur.xl) {
res += square_difference(q.x, cur.xl);
} else if (cur.xr < q.x) {
res += square_difference(q.x, cur.xr);
}
if (q.y < cur.yl) {
res += square_difference(q.y, cur.yl);
} else if (cur.yr < q.y) {
res += square_difference(q.y, cur.yr);
}
return res;
}
int build_range(std::vector<int> &ord, int l, int r, int dep) {
int mid = (l + r) / 2;
bool dim = (dep & 1) == 0;
auto cmp = [&](int arr, int b) {
if (dim) {
if (pt[arr].x != pt[b].x) {
return pt[arr].x < pt[b].x;
}
if (pt[arr].y != pt[b].y) {
return pt[arr].y < pt[b].y;
}
} else {
if (pt[arr].y != pt[b].y) {
return pt[arr].y < pt[b].y;
}
if (pt[arr].x != pt[b].x) {
return pt[arr].x < pt[b].x;
}
}
return arr < b;
};
std::nth_element(ord.begin() + l, ord.begin() + mid, ord.begin() + r, cmp);
int id = int(tr.size());
int pid = ord[mid];
tr.push_back(
{pid, -1, -1, pt[pid].x, pt[pid].x, pt[pid].y, pt[pid].y, pid});
if (l < mid) {
tr[id].ls = build_range(ord, l, mid, dep + 1);
extend_box(tr[id], tr[tr[id].ls]);
}
if (mid + 1 < r) {
tr[id].rs = build_range(ord, mid + 1, r, dep + 1);
extend_box(tr[id], tr[tr[id].rs]);
}
return id;
}
static void extend_box(node &tar, const node &src) {
tar.xl = std::min(tar.xl, src.xl);
tar.xr = std::max(tar.xr, src.xr);
tar.yl = std::min(tar.yl, src.yl);
tar.yr = std::max(tar.yr, src.yr);
tar.min_pid = std::min(tar.min_pid, src.min_pid);
}
static bool box_disjoint(const node &cur, const T &l, const T &r, const T &dn,
const T &top) {
return cur.xr < l || !(cur.xl < r) || cur.yr < dn || !(cur.yl < top);
}
static bool box_contained(const node &cur, const T &l, const T &r,
const T &dn, const T &top) {
return !(cur.xl < l) && cur.xr < r && !(cur.yl < dn) && cur.yr < top;
}
void range_query_at(int id, int end, const T &l, const T &r, const T &dn,
const T &top, int li0, std::vector<int> &res) const {
if (id == -1 || int(res.size()) == li0 ||
box_disjoint(tr[id], l, r, dn, top)) {
return;
}
if (box_contained(tr[id], l, r, dn, top)) {
// Preorder construction makes every subtree a contiguous node range.
int stop = id + std::min(end - id, li0 - int(res.size()));
for (int cur = id; cur < stop; ++cur)
res.push_back(tr[cur].pid);
return;
}
int idx = tr[id].pid;
const point_type &val = pt[idx];
if (!(val.x < l) && val.x < r && !(val.y < dn) && val.y < top) {
res.push_back(idx);
}
int left_end = tr[id].rs == -1 ? end : tr[id].rs;
range_query_at(tr[id].ls, left_end, l, r, dn, top, li0, res);
range_query_at(tr[id].rs, end, l, r, dn, top, li0, res);
}
void nearest_at(int id, const point_type &q, int &bid,
distance_type &bd) const {
if (id == -1) {
return;
}
int idx = tr[id].pid;
distance_type dis = point_distance(idx, q);
if (bid == -1 || dis < bd || (dis == bd && idx < bid)) {
bid = idx;
bd = dis;
}
int arr = tr[id].ls;
int b = tr[id].rs;
distance_type d1 = arr == -1 ? distance_type{} : box_distance(tr[arr], q);
distance_type d2 = b == -1 ? distance_type{} : box_distance(tr[b], q);
if (arr == -1 ||
(b != -1 &&
(d2 < d1 || (d2 == d1 && tr[b].min_pid < tr[arr].min_pid)))) {
std::swap(arr, b);
std::swap(d1, d2);
}
if (arr != -1 && (d1 < bd || (d1 == bd && tr[arr].min_pid < bid))) {
nearest_at(arr, q, bid, bd);
}
if (b != -1 && (d2 < bd || (d2 == bd && tr[b].min_pid < bid))) {
nearest_at(b, q, bid, bd);
}
}
};
} // namespace noya
// NOYA_KD_TREE_HPP
#ifndef NOYA_KD_TREE_HPP
#define NOYA_KD_TREE_HPP 1
/// @complexity Time: Expected O(n log n) build; query time is output/pruning dependent, worst O(n).
/// Space: O(n).
#include "noya/geometry/geometry_base.hpp"
#include <algorithm>
#include <cassert>
#include <numeric>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
/// @brief Static two-dimensional KD-tree for rectangle reporting and nearest
/// neighbor queries.
template <class T, class Distance = std::conditional_t<std::is_integral_v<T>,
__int128_t, long double>>
struct kd_tree {
using point_type = point<T>;
using distance_type = Distance;
struct node {
int pid = -1;
int ls = -1;
int rs = -1;
T xl{};
T xr{};
T yl{};
T yr{};
int min_pid = -1;
};
std::vector<point_type> pt;
std::vector<node> tr;
int rt = -1;
kd_tree() = default;
explicit kd_tree(const std::vector<point_type> &a) { build(a); }
/// @brief Rebuild the balanced tree in O(n log n) expected time.
void build(const std::vector<point_type> &a) {
pt = a;
tr.clear();
rt = -1;
std::vector<int> ord(pt.size());
std::iota(ord.begin(), ord.end(), 0);
if (!ord.empty()) {
tr.reserve(ord.size());
rt = build_range(ord, 0, int(ord.size()), 0);
}
}
/// @brief Return the number of stored points.
int size() const { return int(pt.size()); }
/// @brief Return point indices in [l, r) x [dn, top), stopping
/// after lim results when lim is nonnegative.
std::vector<int> range_query(const T &l, const T &r, const T &dn,
const T &top, int lim = -1) const {
assert(!(r < l) && !(top < dn));
int li0 = lim < 0 ? size() : lim;
std::vector<int> res;
if (li0 == 0 || l == r || dn == top)
return res;
res.reserve(std::min(size(), li0));
range_query_at(rt, int(tr.size()), l, r, dn, top, li0, res);
return res;
}
/// @brief Return (point index, squared distance) for a nearest point; the
/// smallest index wins ties, and -1 denotes an empty tree. Squared distances
/// must be representable in distance_type.
std::pair<int, distance_type> nearest(const point_type &q) const {
int bid = -1;
distance_type bd{};
nearest_at(rt, q, bid, bd);
return {bid, bd};
}
private:
static distance_type square_difference(const T &arr, const T &b) {
distance_type dif =
static_cast<distance_type>(arr) - static_cast<distance_type>(b);
return dif * dif;
}
distance_type point_distance(int idx, const point_type &q) const {
return square_difference(pt[idx].x, q.x) +
square_difference(pt[idx].y, q.y);
}
distance_type box_distance(const node &cur, const point_type &q) const {
distance_type res{};
if (q.x < cur.xl) {
res += square_difference(q.x, cur.xl);
} else if (cur.xr < q.x) {
res += square_difference(q.x, cur.xr);
}
if (q.y < cur.yl) {
res += square_difference(q.y, cur.yl);
} else if (cur.yr < q.y) {
res += square_difference(q.y, cur.yr);
}
return res;
}
int build_range(std::vector<int> &ord, int l, int r, int dep) {
int mid = (l + r) / 2;
bool dim = (dep & 1) == 0;
auto cmp = [&](int arr, int b) {
if (dim) {
if (pt[arr].x != pt[b].x) {
return pt[arr].x < pt[b].x;
}
if (pt[arr].y != pt[b].y) {
return pt[arr].y < pt[b].y;
}
} else {
if (pt[arr].y != pt[b].y) {
return pt[arr].y < pt[b].y;
}
if (pt[arr].x != pt[b].x) {
return pt[arr].x < pt[b].x;
}
}
return arr < b;
};
std::nth_element(ord.begin() + l, ord.begin() + mid, ord.begin() + r, cmp);
int id = int(tr.size());
int pid = ord[mid];
tr.push_back(
{pid, -1, -1, pt[pid].x, pt[pid].x, pt[pid].y, pt[pid].y, pid});
if (l < mid) {
tr[id].ls = build_range(ord, l, mid, dep + 1);
extend_box(tr[id], tr[tr[id].ls]);
}
if (mid + 1 < r) {
tr[id].rs = build_range(ord, mid + 1, r, dep + 1);
extend_box(tr[id], tr[tr[id].rs]);
}
return id;
}
static void extend_box(node &tar, const node &src) {
tar.xl = std::min(tar.xl, src.xl);
tar.xr = std::max(tar.xr, src.xr);
tar.yl = std::min(tar.yl, src.yl);
tar.yr = std::max(tar.yr, src.yr);
tar.min_pid = std::min(tar.min_pid, src.min_pid);
}
static bool box_disjoint(const node &cur, const T &l, const T &r, const T &dn,
const T &top) {
return cur.xr < l || !(cur.xl < r) || cur.yr < dn || !(cur.yl < top);
}
static bool box_contained(const node &cur, const T &l, const T &r,
const T &dn, const T &top) {
return !(cur.xl < l) && cur.xr < r && !(cur.yl < dn) && cur.yr < top;
}
void range_query_at(int id, int end, const T &l, const T &r, const T &dn,
const T &top, int li0, std::vector<int> &res) const {
if (id == -1 || int(res.size()) == li0 ||
box_disjoint(tr[id], l, r, dn, top)) {
return;
}
if (box_contained(tr[id], l, r, dn, top)) {
// Preorder construction makes every subtree a contiguous node range.
int stop = id + std::min(end - id, li0 - int(res.size()));
for (int cur = id; cur < stop; ++cur)
res.push_back(tr[cur].pid);
return;
}
int idx = tr[id].pid;
const point_type &val = pt[idx];
if (!(val.x < l) && val.x < r && !(val.y < dn) && val.y < top) {
res.push_back(idx);
}
int left_end = tr[id].rs == -1 ? end : tr[id].rs;
range_query_at(tr[id].ls, left_end, l, r, dn, top, li0, res);
range_query_at(tr[id].rs, end, l, r, dn, top, li0, res);
}
void nearest_at(int id, const point_type &q, int &bid,
distance_type &bd) const {
if (id == -1) {
return;
}
int idx = tr[id].pid;
distance_type dis = point_distance(idx, q);
if (bid == -1 || dis < bd || (dis == bd && idx < bid)) {
bid = idx;
bd = dis;
}
int arr = tr[id].ls;
int b = tr[id].rs;
distance_type d1 = arr == -1 ? distance_type{} : box_distance(tr[arr], q);
distance_type d2 = b == -1 ? distance_type{} : box_distance(tr[b], q);
if (arr == -1 ||
(b != -1 &&
(d2 < d1 || (d2 == d1 && tr[b].min_pid < tr[arr].min_pid)))) {
std::swap(arr, b);
std::swap(d1, d2);
}
if (arr != -1 && (d1 < bd || (d1 == bd && tr[arr].min_pid < bid))) {
nearest_at(arr, q, bid, bd);
}
if (b != -1 && (d2 < bd || (d2 == bd && tr[b].min_pid < bid))) {
nearest_at(b, q, bid, bd);
}
}
};
} // namespace noya
#endif // NOYA_KD_TREE_HPP
/// @complexity Time: Expected O(n log n) build; query time is output/pruning dependent, worst O(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 <numeric>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
/// @brief Static two-dimensional KD-tree for rectangle reporting and nearest
/// neighbor queries.
template <class T, class Distance = std::conditional_t<std::is_integral_v<T>,
__int128_t, long double>>
struct kd_tree {
using point_type = point<T>;
using distance_type = Distance;
struct node {
int pid = -1;
int ls = -1;
int rs = -1;
T xl{};
T xr{};
T yl{};
T yr{};
int min_pid = -1;
};
std::vector<point_type> pt;
std::vector<node> tr;
int rt = -1;
kd_tree() = default;
explicit kd_tree(const std::vector<point_type> &a) { build(a); }
/// @brief Rebuild the balanced tree in O(n log n) expected time.
void build(const std::vector<point_type> &a) {
pt = a;
tr.clear();
rt = -1;
std::vector<int> ord(pt.size());
std::iota(ord.begin(), ord.end(), 0);
if (!ord.empty()) {
tr.reserve(ord.size());
rt = build_range(ord, 0, int(ord.size()), 0);
}
}
/// @brief Return the number of stored points.
int size() const { return int(pt.size()); }
/// @brief Return point indices in [l, r) x [dn, top), stopping
/// after lim results when lim is nonnegative.
std::vector<int> range_query(const T &l, const T &r, const T &dn,
const T &top, int lim = -1) const {
assert(!(r < l) && !(top < dn));
int li0 = lim < 0 ? size() : lim;
std::vector<int> res;
if (li0 == 0 || l == r || dn == top)
return res;
res.reserve(std::min(size(), li0));
range_query_at(rt, int(tr.size()), l, r, dn, top, li0, res);
return res;
}
/// @brief Return (point index, squared distance) for a nearest point; the
/// smallest index wins ties, and -1 denotes an empty tree. Squared distances
/// must be representable in distance_type.
std::pair<int, distance_type> nearest(const point_type &q) const {
int bid = -1;
distance_type bd{};
nearest_at(rt, q, bid, bd);
return {bid, bd};
}
private:
static distance_type square_difference(const T &arr, const T &b) {
distance_type dif =
static_cast<distance_type>(arr) - static_cast<distance_type>(b);
return dif * dif;
}
distance_type point_distance(int idx, const point_type &q) const {
return square_difference(pt[idx].x, q.x) +
square_difference(pt[idx].y, q.y);
}
distance_type box_distance(const node &cur, const point_type &q) const {
distance_type res{};
if (q.x < cur.xl) {
res += square_difference(q.x, cur.xl);
} else if (cur.xr < q.x) {
res += square_difference(q.x, cur.xr);
}
if (q.y < cur.yl) {
res += square_difference(q.y, cur.yl);
} else if (cur.yr < q.y) {
res += square_difference(q.y, cur.yr);
}
return res;
}
int build_range(std::vector<int> &ord, int l, int r, int dep) {
int mid = (l + r) / 2;
bool dim = (dep & 1) == 0;
auto cmp = [&](int arr, int b) {
if (dim) {
if (pt[arr].x != pt[b].x) {
return pt[arr].x < pt[b].x;
}
if (pt[arr].y != pt[b].y) {
return pt[arr].y < pt[b].y;
}
} else {
if (pt[arr].y != pt[b].y) {
return pt[arr].y < pt[b].y;
}
if (pt[arr].x != pt[b].x) {
return pt[arr].x < pt[b].x;
}
}
return arr < b;
};
std::nth_element(ord.begin() + l, ord.begin() + mid, ord.begin() + r, cmp);
int id = int(tr.size());
int pid = ord[mid];
tr.push_back(
{pid, -1, -1, pt[pid].x, pt[pid].x, pt[pid].y, pt[pid].y, pid});
if (l < mid) {
tr[id].ls = build_range(ord, l, mid, dep + 1);
extend_box(tr[id], tr[tr[id].ls]);
}
if (mid + 1 < r) {
tr[id].rs = build_range(ord, mid + 1, r, dep + 1);
extend_box(tr[id], tr[tr[id].rs]);
}
return id;
}
static void extend_box(node &tar, const node &src) {
tar.xl = std::min(tar.xl, src.xl);
tar.xr = std::max(tar.xr, src.xr);
tar.yl = std::min(tar.yl, src.yl);
tar.yr = std::max(tar.yr, src.yr);
tar.min_pid = std::min(tar.min_pid, src.min_pid);
}
static bool box_disjoint(const node &cur, const T &l, const T &r, const T &dn,
const T &top) {
return cur.xr < l || !(cur.xl < r) || cur.yr < dn || !(cur.yl < top);
}
static bool box_contained(const node &cur, const T &l, const T &r,
const T &dn, const T &top) {
return !(cur.xl < l) && cur.xr < r && !(cur.yl < dn) && cur.yr < top;
}
void range_query_at(int id, int end, const T &l, const T &r, const T &dn,
const T &top, int li0, std::vector<int> &res) const {
if (id == -1 || int(res.size()) == li0 ||
box_disjoint(tr[id], l, r, dn, top)) {
return;
}
if (box_contained(tr[id], l, r, dn, top)) {
// Preorder construction makes every subtree a contiguous node range.
int stop = id + std::min(end - id, li0 - int(res.size()));
for (int cur = id; cur < stop; ++cur)
res.push_back(tr[cur].pid);
return;
}
int idx = tr[id].pid;
const point_type &val = pt[idx];
if (!(val.x < l) && val.x < r && !(val.y < dn) && val.y < top) {
res.push_back(idx);
}
int left_end = tr[id].rs == -1 ? end : tr[id].rs;
range_query_at(tr[id].ls, left_end, l, r, dn, top, li0, res);
range_query_at(tr[id].rs, end, l, r, dn, top, li0, res);
}
void nearest_at(int id, const point_type &q, int &bid,
distance_type &bd) const {
if (id == -1) {
return;
}
int idx = tr[id].pid;
distance_type dis = point_distance(idx, q);
if (bid == -1 || dis < bd || (dis == bd && idx < bid)) {
bid = idx;
bd = dis;
}
int arr = tr[id].ls;
int b = tr[id].rs;
distance_type d1 = arr == -1 ? distance_type{} : box_distance(tr[arr], q);
distance_type d2 = b == -1 ? distance_type{} : box_distance(tr[b], q);
if (arr == -1 ||
(b != -1 &&
(d2 < d1 || (d2 == d1 && tr[b].min_pid < tr[arr].min_pid)))) {
std::swap(arr, b);
std::swap(d1, d2);
}
if (arr != -1 && (d1 < bd || (d1 == bd && tr[arr].min_pid < bid))) {
nearest_at(arr, q, bid, bd);
}
if (b != -1 && (d2 < bd || (d2 == bd && tr[b].min_pid < bid))) {
nearest_at(b, q, bid, bd);
}
}
};
} // namespace noya
// NOYA_KD_TREE_HPP