189# ifndef GEOM_ALGORITHMS_H
190# define GEOM_ALGORITHMS_H
240 verts.append(it.get_current_vertex());
252 for (
size_t i = 0; i <
verts.size(); ++i)
296 if (
verts.size() < 3)
300 const size_t n =
verts.size();
301 for (
size_t i = 0; i < n; ++i)
308 const int curr = turn > 0 ? 1 : -1;
311 else if (
sign != curr)
320 for (
size_t i = 0; i <
verts.size() / 2; ++i)
356 const size_t tri,
const size_t third)
360 const size_t tmp = a;
369 template <
typename TriangleArray,
typename GroupFn>
374 edges.
reserve(triangles.size() * 3);
377 const auto & tri = triangles(
tri_idx);
393 while (first < edges.
size())
395 size_t last = first + 1;
397 edges(last).u == edges(first).u
and
398 edges(last).v == edges(first).v)
443 const size_t tmp = u;
449 if (
const auto *ptr = boundary.
search(e); ptr !=
nullptr)
456 template <
typename IndexedTriangle,
typename InCirclePredicate>
470 for (
size_t i = 1; i < n; ++i)
498 const size_t i1 = n + 1;
499 const size_t i2 = n + 2;
509 for (
size_t t = 0; t < work.
size(); ++t)
511 const auto & [a, b, c, alive] = work(t);
515 if (point_in_circumcircle(
pts, a, b, c,
pidx))
547 const size_t tmp = u;
557 for (
size_t t = 0; t < work.
size(); ++t)
564 for (
size_t t = 0; t < work.
size(); ++t)
566 const auto & [a, b, c, alive] = work(t);
570 if (a >= n
or b >= n
or c >= n)
576 out.append(IndexedTriangle{a, b, c});
662 for (
size_t i =
verts.size(); i > 0; --i)
682 if (
Segment curr = it.get_current_segment();
741 Vertex & curr = it.get_current_vertex();
775 <<
"No valid ear found; polygon may be non-simple";
783 <<
"No valid ear found during triangulation";
914 return {a, b,
dist2(a, b)};
928 const size_t l,
const size_t r)
934 for (
size_t i =
l; i <
r; ++i)
935 for (
size_t j = i + 1; j <
r; ++j)
958 const size_t n =
r -
l;
962 const size_t mid =
l + n / 2;
971 for (
size_t i = 0; i <
py.size(); ++i)
995 for (
size_t i = 0; i <
py.size(); ++i)
1002 for (
size_t i = 0; i <
strip.size(); ++i)
1003 for (
size_t j = i + 1; j <
strip.size(); ++j)
1031 <<
"Closest pair requires at least 2 points";
1036 for (
size_t i = 1; i <
px.size(); ++i)
1037 if (
px(i) ==
px(i - 1))
1038 return {
px(i - 1),
px(i), 0};
1184 <<
"MinimumEnclosingCircle: empty point set";
1188 std::random_device rd;
1189 std::mt19937
gen(rd());
1190 for (
size_t i =
pts.size() - 1; i > 0; --i)
1192 std::uniform_int_distribution<size_t> dist(0, i);
1193 if (
const size_t j = dist(
gen); i != j)
1194 std::swap(
pts(i),
pts(j));
1205 for (
const auto & p: points)
1217 for (
size_t j = 0; j < n; ++j)
1230 for (
size_t k = 0;
k < n; ++
k)
1241 for (
size_t i = 1; i <
pts.size(); ++i)
1344 const size_t n =
verts.size();
1364 for (
size_t i = 0; i < n; ++i)
1366 const size_t ni = (i + 1) % n;
1399 const size_t n =
verts.size();
1414 return verts(i).distance_squared_to(
verts((i + 1) % n));
1429 for (
size_t i = 0; i < n; ++i)
1431 const size_t ni = (i + 1) % n;
1594 if (
verts.size() < 3)
1640 if (
out.get_last() == p)
1643 while (
out.size() >= 2)
1653 static_cast<void>(
out.remove_last());
1663 for (
size_t i = 0; i <
pts.size(); ++i)
1667 static_cast<void>(
ret.remove_last());
1677 for (
size_t i = 0; i <
clean.size(); ++i)
1680 if (
ret.size() >= 3)
1712 for (
size_t i = 0; i <
clp.size(); ++i)
1727 for (
size_t j = 0; j <
input.size(); ++j)
1808 return h.dy() > 0
or h.dy() == 0
and h.dx() >= 0;
1814 return a.
dx() * b.
dy() - a.
dy() * b.
dx();
1820 return a.
dx() * b.
dx() + a.
dy() * b.
dy();
1860 if (
out.get_last() == p)
1863 while (
out.size() >= 2)
1873 static_cast<void>(
out.remove_last());
1883 for (
size_t i = 0; i <
pts.size(); ++i)
1887 static_cast<void>(
ret.remove_last());
1897 for (
size_t i = 0; i <
clean.size(); ++i)
1900 if (
ret.size() >= 3)
1924 verts.append(it.get_current_vertex());
1932 for (
size_t i = 0; i <
verts.size(); ++i)
1969 for (
size_t i = 0; i <
hps.size(); ++i)
1980 if (
hp.offset() > last.offset())
1982 static_cast<void>(
unique.remove_last());
1994 for (
size_t i = 0; i <
unique.size(); ++i)
2014 if (
hp.offset() >
dq.get_last().offset())
2027 if (
not dq.is_empty())
2057 while (
not tmp.is_empty())
2076 return (*
this)(
hps);
2158 for (
size_t i = 2; i <
pts.size(); ++i)
2220 for (
size_t i = 0; i <
all.size(); ++i)
2238 const size_t n =
ret.sites.size();
2245 GeomBowyerWatsonUtils::triangulate<IndexedTriangle>(
2248 const size_t ia,
const size_t ib,
2249 const size_t ic,
const size_t ip)
2270 return (*
this)(points);
2364 for (
size_t i = 2; i < s.
size(); ++i)
2365 if (
orientation(s(0).position, s(1).position, s(i).position)
2390 const size_t ia,
const size_t ib,
const size_t ic,
const size_t ip)
2415 if (
det > 0)
return true;
2416 if (
det < 0)
return false;
2420 if (
det < 0)
return true;
2421 if (
det > 0)
return false;
2436 for (
size_t i = 0; i <
input.size(); ++i)
2446 for (
size_t i = 0; i <
all.size(); ++i)
2447 if (
ret.is_empty()
or ret.get_last().position !=
all(i).position)
2464 const size_t n =
ret.sites.size();
2473 for (
size_t i = 0; i < n; ++i)
2475 pts.append(
ret.sites(i).position);
2476 wts.append(
ret.sites(i).weight);
2481 ret.triangles = GeomBowyerWatsonUtils::triangulate<IndexedTriangle>(
2484 const size_t ia,
const size_t ib,
2485 const size_t ic,
const size_t ip)
2538 for (
int i = 0; i < 3; ++i)
2539 if (t.
v[i] ==
id)
return static_cast<size_t>(i);
2546 for (
int i = 0; i < 3; ++i)
2547 if (t.
adj[i] == n)
return static_cast<size_t>(i);
2553 const Tri & t,
const size_t pidx)
2569 const size_t ia,
const size_t ib,
2570 const size_t ic,
const size_t ip)
2577 if (
det > 0)
return true;
2578 if (
det < 0)
return false;
2582 if (
det < 0)
return true;
2583 if (
det > 0)
return false;
2599 const size_t pidx,
const size_t root)
2605 for (
size_t c = 0; c < dag(
cur).children.
size(); ++c)
2623 for (
const unsigned long nb:
tris(
ot).adj)
2625 if (
nb ==
NONE)
continue;
2647 ret.sites.reserve(
all.size());
2648 for (
size_t i = 0; i <
all.size(); ++i)
2649 if (
ret.sites.is_empty()
or ret.sites.get_last() !=
all(i))
2650 ret.sites.append(
all(i));
2653 const size_t n =
ret.sites.size();
2670 for (
size_t i = 1; i < n; ++i)
2689 for (
size_t i = 0; i < n; ++i)
2693 std::random_device rd;
2694 std::mt19937
gen(rd());
2695 for (
size_t i = n - 1; i > 0; --i)
2697 std::uniform_int_distribution<size_t>
dis(0, i);
2698 const size_t j =
dis(
gen);
2699 const size_t tmp = order(i);
2700 order(i) = order(j);
2720 for (
size_t oi = 0;
oi < n; ++
oi)
2722 const size_t pidx = order(
oi);
2763 for (
int e = 0; e < 3; ++e)
2780 const size_t V =
pts.size();
2784 for (
size_t i = 0; i <
V; ++i)
2793 for (
size_t bi = 0;
bi < boundary.
size(); ++
bi)
2795 size_t u = boundary(
bi).u, v = boundary(
bi).v;
2796 const size_t ext = boundary(
bi).ext;
2797 const size_t old_ct = boundary(
bi).old_ct;
2802 const size_t tmp = u;
2829 const size_t u =
tris(
nt).v[0];
2830 const size_t v =
tris(
nt).v[1];
2845 ret.triangles.reserve(2 * n);
2846 for (
size_t t = 0; t <
tris.size(); ++t)
2849 if (
not tr.alive)
continue;
2850 if (
tr.v[0] >= n
or tr.v[1] >= n
or tr.v[2] >= n)
continue;
2864 return (*
this)(points);
2939 size_t lo = 0, hi =
pts.size();
2953 for (
int i = 0; i < 3; ++i)
2954 if (t.
v[i] ==
id)
return static_cast<size_t>(i);
2961 for (
int i = 0; i < 3; ++i)
2962 if (t.
adj[i] == n)
return static_cast<size_t>(i);
2969 const size_t a,
const size_t b)
2971 for (
int i = 0; i < 3; ++i)
2973 const size_t e0 = t.
v[(i + 1) % 3];
2974 const size_t e1 = t.
v[(i + 2) % 3];
2976 return static_cast<size_t>(i);
2986 for (
size_t i = 0; i <
dt_tris.size(); ++i)
2992 {
false,
false,
false},
3001 const size_t first,
const size_t last)
3003 if (last - first != 2)
3006 const auto &
e0 = edges(first);
3007 const auto &
e1 = edges(first + 1);
3009 const size_t t0 =
e0.tri;
3010 const size_t t1 =
e1.tri;
3025 const size_t u,
const size_t v,
3029 for (
size_t t = 0; t <
tris.size(); ++t)
3051 const size_t a,
const size_t b,
3052 const size_t c,
const size_t d)
3091 const size_t p =
ta.v[(
la + 1) % 3];
3092 const size_t q =
ta.v[(
la + 2) % 3];
3101 const bool con_a1 =
ta.constrained[(
la + 1) % 3];
3102 const bool con_a2 =
ta.constrained[(
la + 2) % 3];
3103 const bool con_b1 =
tb.constrained[(
lb + 1) % 3];
3104 const bool con_b2 =
tb.constrained[(
lb + 2) % 3];
3109 const size_t bp =
tb.v[(
lb + 1) % 3];
3110 const size_t bq =
tb.v[(
lb + 2) % 3];
3226 const size_t u,
const size_t v)
3231 for (
size_t t = 0; t <
tris.size(); ++t)
3236 for (
int e = 0; e < 3; ++e)
3238 const size_t w1 =
tris(t).v[(e + 1) % 3];
3239 const size_t w2 =
tris(t).v[(e + 2) % 3];
3247 const size_t nb =
tris(t).adj[e];
3267 const size_t u,
const size_t v)
3269 for (
size_t t = 0; t <
tris.size(); ++t)
3274 if (
const size_t loc = edge_opposite(
tris(t), u, v); loc != NONE)
3275 tris(t).constrained[loc] =
true;
3281 const size_t u,
const size_t v)
3285 if (edge_exists(
tris, u, v,
et,
el))
3287 mark_constrained(
tris, u, v);
3296 if (edge_exists(
tris, u, v,
et,
el))
3298 mark_constrained(
tris, u, v);
3306 if (edge_exists(
tris, u, v,
et,
el))
3307 mark_constrained(
tris, u, v);
3336 if (is_convex_quad(
pts, a, b, c, d))
3353 if (edge_exists(
tris, u, v,
et,
el))
3354 mark_constrained(
tris, u, v);
3370 for (
size_t t = 0; t <
tris.size(); ++t)
3375 for (
int e = 0; e < 3; ++e)
3377 if (
tris(t).constrained[e])
3380 const size_t nb =
tris(t).adj[e];
3381 if (
nb == NONE
or nb <= t)
3388 const size_t lb = adj_of(
tris(
nb), t);
3415 const size_t a =
tris(t).v[e];
3416 const size_t b =
tris(t).v[(e + 1) % 3];
3417 const size_t c =
tris(t).v[(e + 2) % 3];
3419 if (is_convex_quad(
pts, a, b, c, d))
3443 const Segment & s = it.get_curr();
3450 for (
size_t i = 0; i <
con_arr.size(); ++i)
3451 for (
size_t j = i + 1; j <
con_arr.size(); ++j)
3457 return lexicographic_less(a, b);
3462 for (
size_t i = 0; i <
all.size(); ++i)
3478 const Segment & seg = it.get_curr();
3479 const size_t u = find_point_index(sites, seg.
get_src_point());
3480 const size_t v = find_point_index(sites, seg.
get_tgt_point());
3482 if (u == NONE
or v == NONE
or u == v)
3488 for (
size_t i = 0; i < sites.
size(); ++i)
3490 if (i == u
or i == v)
3498 if (
chain.size() > 2)
3504 (sites(a).get_x() -
pu.get_x()) *
3505 (sites(a).get_x() -
pu.get_x()) +
3506 (sites(a).get_y() -
pu.get_y()) *
3507 (sites(a).get_y() -
pu.get_y());
3509 (sites(b).get_x() -
pu.get_x()) *
3510 (sites(b).get_x() -
pu.get_x()) +
3511 (sites(b).get_y() -
pu.get_y()) *
3512 (sites(b).get_y() -
pu.get_y());
3517 for (
size_t i = 0; i + 1 <
chain.size(); ++i)
3523 for (
size_t j = 0; j < result.
size(); ++j)
3524 if (result(j).u == a
and result(j).v == b)
3552 const size_t n =
ret.sites.size();
3570 for (
size_t i = 0; i < n; ++i)
3576 if (triangles.is_empty())
3585 build_adjacency(
tris, triangles);
3600 ret.triangles.reserve(
tris.size());
3601 for (
size_t t = 0; t <
tris.size(); ++t)
3627 const std::initializer_list<Segment>
constraints)
const
3630 for (
const Point & p: points)
3646 for (
size_t i = 0; i < result.
triangles.size(); ++i)
3776 if (
verts.size() < 3)
3792 const Point q(
mid.get_x() - dy,
mid.get_y() + dx);
3801 <<
"Sites and clipped polygons size mismatch";
3806 for (
size_t i = 0; i < sites.
size(); ++i)
3810 cell.
site = sites(i);
3812 ret.append(std::move(cell));
3831 if (
dt.triangles.is_empty())
3836 for (
size_t idx = 0; idx <
dt.triangles.size(); ++idx)
3838 const auto & [i, j,
k] =
dt.triangles(idx);
3845 for (
size_t i = 0; i <
dt.sites.size(); ++i)
3852 const size_t first,
const size_t last)
3854 if (const size_t cnt = last - first; cnt >= 2)
3856 const Point & p1 = centers(edges(first).tri);
3857 const Point & p2 = centers(edges(first + 1).tri);
3858 ret.edges.append(Edge{
3859 edges(first).u, edges(first).v,
3860 p1, p2, false, Point()
3865 const auto & edge = edges(first);
3866 const size_t u = edge.u;
3867 const size_t v = edge.v;
3868 const size_t tri = edge.tri;
3869 const size_t third = edge.third;
3897 for (
size_t s = 0; s <
dt.sites.size(); ++s)
3899 for (
size_t t = 0; t <
dt.triangles.size(); ++t)
3901 const auto & [
ti,
tj,
tk] =
dt.triangles(t);
3907 ret.cells.reserve(
dt.sites.size());
3908 for (
size_t s = 0; s <
dt.sites.size(); ++s)
3912 for (
size_t idx = 0; idx <
incidence(s).size(); ++idx)
3915 const Point site =
dt.sites(s);
3916 if (
verts.size() > 1)
3939 for (
size_t k = 0;
k <
verts.size(); ++
k)
3943 static_cast<void>(
clean.remove_last());
3946 cell.site_index = s;
3948 cell.bounded =
on_hull(s) == 0;
3949 cell.vertices =
clean;
3950 ret.cells.append(std::move(cell));
3975 return (*
this)(delaunay(
il));
3999 const size_t n = sites.
size();
4006 for (
size_t i = 0; i < n; ++i)
4014 for (
size_t di = 0;
di <
dt.sites.size(); ++
di)
4017 for (
size_t oi = 0;
oi < n; ++
oi)
4018 if (
dt.sites(
di) == sites(
oi))
4028 for (
size_t i = 0; i < n; ++i)
4030 for (
size_t t = 0; t <
dt.triangles.size(); ++t)
4046 for (
size_t i = 0; i < n; ++i)
4053 if (sites(j) != sites(i))
4054 hps.append(bisector_halfplane_for_site(sites(i), sites(j)));
4075 return clipped_cells(
vor.sites,
clip);
4101 return clipped_cells(delaunay(
il).sites,
clip);
4114 return indexed_clipped_cells(sites, clipped_cells(sites,
clip));
4127 return clipped_cells_indexed(
vor.sites,
clip);
4140 return clipped_cells_indexed(delaunay(
point_set).sites,
clip);
4154 return clipped_cells_indexed(delaunay(
il).sites,
clip);
4230 return (*
this)(
pts);
4307 if (a->
y != b->
y)
return a->
y > b->
y;
4308 if (a->
x != b->
x)
return a->
x < b->
x;
4310 return a->
id < b->
id;
4325 if (
a !=
o.a)
return a <
o.
a;
4326 if (
b !=
o.b)
return b <
o.b;
4331 static constexpr double kEps = 1e-9;
4346 for (
size_t i = 2; i <
pts.size(); ++i)
4359 std::swap(t.j, t.k);
4373 return (
px +
qx) / 2.0;
4382 const double a = 1.0 /
z0 - 1.0 /
z1;
4383 const double b = -2.0 * (
px /
z0 -
qx /
z1);
4387 if (std::fabs(a) <=
kEps)
4390 double disc = b * b - 4.0 * a * c;
4395 const double sq = std::sqrt(
disc);
4396 const double x1 = (-b - sq) / (2.0 * a);
4397 const double x2 = (-b + sq) / (2.0 * a);
4398 return py <
qy ? std::max(x1, x2) : std::min(x1, x2);
4405 if (head ==
nullptr)
4408 for (
Arc *arc = head; arc !=
nullptr; arc = arc->
next)
4410 const double left = arc->
prev ==
nullptr ?
4411 -std::numeric_limits<double>::infinity() :
4413 const double right = arc->next ==
nullptr ?
4414 std::numeric_limits<double>::infinity() :
4422 while (tail->
next !=
nullptr)
4429 if (arc !=
nullptr and arc->
circle !=
nullptr)
4445 if (arc ==
nullptr or arc->
prev ==
nullptr or arc->
next ==
nullptr)
4449 const size_t ib = arc->
site;
4468 const double det = 2.0 * (
ax * (
by - cy) +
bx * (cy -
ay) + cx * (
ay -
by));
4472 const double a2 =
ax *
ax +
ay *
ay;
4473 const double b2 =
bx *
bx +
by *
by;
4474 const double c2 = cx * cx + cy * cy;
4476 const double ux = (a2 * (
by - cy) + b2 * (cy -
ay) +
c2 * (
ay -
by)) /
det;
4477 const double uy = (a2 * (cx -
bx) + b2 * (
ax - cx) +
c2 * (
bx -
ax)) /
det;
4479 const double r = std::hypot(
ax -
ux,
ay -
uy);
4484 event_pool.append(std::make_unique<Event>());
4488 ev->is_site =
false;
4501 for (
size_t t = 0; t <
dt.triangles.size(); ++t)
4503 const auto & [i, j,
k] =
dt.triangles(t);
4504 const Point & a =
dt.sites(i);
4505 const Point & b =
dt.sites(j);
4511 for (
size_t p = 0; p <
dt.sites.size(); ++p)
4513 if (p == i
or p == j
or p ==
k)
4534 const size_t n = sites.
size();
4544 arc_pool.append(std::make_unique<Arc>());
4551 for (
size_t i = 0; i < n; ++i)
4553 event_pool.append(std::make_unique<Event>());
4563 Arc *head =
nullptr;
4571 if (i == j
or j ==
k or i ==
k)
4576 size_t a = i, b = j, c =
k;
4577 if (a > b) std::swap(a, b);
4578 if (b > c) std::swap(b, c);
4579 if (a > b) std::swap(a, b);
4581 if (
const auto ptr =
seen.insert(
TriKey{a, b, c}); ptr ==
nullptr)
4593 const double ly =
ev->y;
4600 if (head ==
nullptr)
4615 if (right->
next !=
nullptr)
4634 if (left ==
nullptr or right ==
nullptr)
4651 out.triangles = std::move(
tris);
4687 return (*
this)(
pts);
4767 points.
append(it.get_curr());
4769 if (points.
size() <= 1)
4776 unique.append(points(0));
4777 for (
size_t i = 1; i < points.
size(); ++i)
4779 unique.append(points(i));
4803 const size_t n = points.
size();
4810 ret.add_vertex(points(0));
4816 for (
size_t i = 0; i < n; ++i)
4818 while (
lower.size() >= 2 &&
4820 static_cast<void>(
lower.remove_last());
4821 lower.append(points(i));
4826 for (
size_t i = n; i > 0; --i)
4828 const Point & p = points(i - 1);
4829 while (
upper.size() >= 2 &&
4831 static_cast<void>(
upper.remove_last());
4836 static_cast<void>(
lower.remove_last());
4837 static_cast<void>(
upper.remove_last());
4839 for (
size_t i = 0; i <
lower.size(); ++i)
4842 for (
size_t i = 0; i <
upper.size(); ++i)
4845 if (
ret.size() >= 3)
4902 points.
append(it.get_curr());
4904 if (points.
size() <= 1)
4911 unique.append(points(0));
4912 for (
size_t i = 1; i < points.
size(); ++i)
4914 unique.append(points(i));
4938 const size_t n = points.
size();
4945 ret.add_vertex(points(0));
4951 for (
size_t i = 1; i < n; ++i)
4953 if (points(i).get_y() < points(
pivot_idx).get_y())
4959 if (points(i).get_y() == points(
pivot_idx).get_y() &&
4960 points(i).get_x() < points(
pivot_idx).get_x())
4974 for (
size_t i = 1; i < n; ++i)
4986 return pivot.distance_squared_to(a) <
4987 pivot.distance_squared_to(b);
4993 for (
size_t i = 0; i < polar.
size();)
4996 while (j + 1 < polar.
size() &&
5016 for (
size_t i = 1; i <
filtered.size(); ++i)
5019 while (
hull.size() >= 2 &&
5021 static_cast<void>(
hull.remove_last());
5025 for (
size_t i = 0; i <
hull.size(); ++i)
5028 if (
ret.size() >= 3)
5100 if (
const Point & p = it.get_curr(); p.is_right_of(s))
5156 ret.add_vertex(first);
5176 <<
"BruteForceConvexHull: broken chain (degenerate input?)";
5186 if (
ret.size() >= 3)
5262 const Point *current = start;
5266 ret.add_vertex(*current);
5276 if (
next ==
nullptr)
5290 if (
next ==
nullptr)
5295 while (current != start);
5297 if (
ret.size() >= 3)
5362 const Point & p = it.get_curr();
5413 ret.first.append(p);
5418 ret.second.append(p);
5480 return std::make_pair(leftmost, rightmost);
5502 ret.second.append(p);
5522 if (leftmost == rightmost)
5524 ret.add_vertex(leftmost);
5540 ret.add_vertex(it.get_curr());
5542 if (
ret.size() >= 3)
5637 Point a0 =
s1.get_src_point();
5638 Point a1 =
s1.get_tgt_point();
5639 Point b0 =
s2.get_src_point();
5640 Point b1 =
s2.get_tgt_point();
5666 if (
not s1.is_parallel_with(
s2))
5669 out.point =
s1.intersection_with(
s2);
5674 if (
out.overlap.get_src_point() ==
out.overlap.get_tgt_point())
5677 out.point =
out.overlap.get_src_point();
5740 template <
typename Event,
typename CmpEvent>
5800 return std::move(
se.event);
5806 return queue_.min().event;
5828 template <
typename Handler>
5841 template <
typename Handler>
5912 return (
y1 + y2) / 2;
5914 return y1 + (x - x1) * (y2 -
y1) / (x2 - x1);
5973 return static_cast<int>(a.
type) <
static_cast<int>(b.
type);
6050 if (a ==
nullptr)
return b;
6051 if (b ==
nullptr)
return a;
6063 if (
root ==
nullptr)
6068 Node *left =
nullptr;
6069 Node *right =
nullptr;
6072 node->
right = right;
6086 if (
root ==
nullptr)
6093 if (
root->label < label)
6109 if (
root ==
nullptr)
6118 if (
root->label < label)
6142 while (
cur !=
nullptr)
6158 while (
cur !=
nullptr)
6185 for (
size_t i = 0; i <
segs.size(); ++i)
6196 return nodes_(seg) !=
nullptr;
6206 const auto node =
new Node{
6208 static_cast<unsigned>(
rng_()),
6234 while (
cur !=
nullptr)
6235 if (
cur->label < label)
6252 const Node *succ =
nullptr;
6254 while (
cur !=
nullptr)
6284 const size_t i,
const size_t j,
6290 const size_t lo = i < j ? i : j;
6291 const size_t hi = i < j ? j : i;
6292 const size_t key = lo * n + hi;
6299 if (
not hit.intersects())
6315 if (
ip.get_x() < sx)
6334 const size_t n = segments.
size();
6343 for (
size_t i = 0; i < n; ++i)
6346 <<
"Segment " << i <<
" is degenerate (zero length)";
6352 for (
size_t i = 0; i < n; ++i)
6362 auto pair_key = [n](
const size_t i,
const size_t j)
6364 const size_t lo = i < j ? i : j;
6365 const size_t hi = i < j ? j : i;
6378 const size_t lo = i < j ? i : j;
6379 const size_t hi = i < j ? j : i;
6386 while (
eq.has_events())
6393 const size_t idx =
ev.seg_a;
6403 if (
const size_t idx =
ev.seg_a; status.
contains(idx))
6406 const size_t succ = status.
successor(idx);
6414 const size_t a =
ev.seg_a;
6415 const size_t b =
ev.seg_b;
6448 if (block.
size() > 2)
6450 for (
size_t i = 0; i + 1 < block.
size(); ++i)
6451 for (
size_t j = i + 1; j < block.
size(); ++j)
6455 if (block.
size() == 0)
6462 for (
size_t i = 0; i < block.
size(); ++i)
6463 status.
erase(block(i));
6464 for (
size_t i = 0; i < block.
size(); ++i)
6465 status.
insert(block(i), sx);
6467 for (
size_t i = 0; i < block.
size(); ++i)
6469 const size_t s = block(i);
6618 if (
root ==
nullptr)
6625 if (
root->label < label)
6639 if (a ==
nullptr)
return b;
6640 if (b ==
nullptr)
return a;
6652 if (
root ==
nullptr)
6657 Node *left =
nullptr;
6658 Node *right =
nullptr;
6661 node->
right = right;
6676 if (
root ==
nullptr)
6685 if (
root->label < label)
6709 while (
cur !=
nullptr)
6725 while (
cur !=
nullptr)
6752 for (
size_t i = 0; i <
verts.size(); ++i)
6763 return nodes_(edge) !=
nullptr;
6773 const auto node =
new Node{
6775 static_cast<unsigned>(
rng_()),
6798 while (
cur !=
nullptr)
6830 for (
size_t i = 0; i < n; ++i)
6843 const size_t top =
sorted(0);
6848 for (
size_t i = 0; i < n; ++i)
6852 for (
size_t i = top; i !=
bot; i = (i + 1) % n)
6861 for (
size_t i = 2; i < n - 1; ++i)
6863 const size_t curr =
sorted(i);
6868 while (stack.
size() > 1)
6870 const size_t a = stack.
top();
6871 const size_t b = stack.
top(1);
6887 const size_t peek = stack.
top();
6909 const size_t last =
sorted(n - 1);
6910 while (stack.
size() > 1)
6912 const size_t a = stack.
top();
6913 const size_t b = stack.
top(1);
6924 const size_t n = v.
size();
6925 const Point & prev = v((i + n - 1) % n);
6926 const Point & curr = v(i);
6952 const size_t n =
verts.size();
6955 for (
size_t i = 0; i < n; ++i)
6964 for (
size_t i = 0; i < n; ++i)
6970 for (
size_t v = 0; v < n; ++v)
6972 auto &
nbrs = adj(v);
6976 [&
verts, v](
const size_t a,
const size_t b)
6997 using HalfEdge = std::pair<size_t, size_t>;
6999 for (
size_t v = 0; v < n; ++v)
7001 const auto &
nbrs = adj(v);
7002 if (
nbrs.is_empty())
7004 for (
size_t k = 0;
k <
nbrs.size(); ++
k)
7006 const size_t u =
nbrs(
k);
7015 for (
size_t i = 0; i < face.size(); ++i)
7018 const Point & b =
verts(face((i + 1) % face.size()));
7028 const HalfEdge start =
kv;
7029 if (visited.
search(start) !=
nullptr)
7032 HalfEdge
cur = start;
7055 for (
const size_t idx: face)
7066 const size_t n =
verts.size();
7076 for (
size_t i = 0; i < n; ++i)
7084 for (
size_t i = 0; i < n; ++i)
7089 for (
size_t i = 0; i < n; ++i)
7100 for (
size_t i = 0; i < n; ++i)
7106 auto add_diagonal = [&](
const size_t a,
const size_t b)
7110 if ((a + 1) % n == b
or (b + 1) % n == a)
7113 const size_t lo = a < b ? a : b;
7114 const size_t hi = a < b ? b : a;
7124 for (
size_t oi = 0;
oi < n; ++
oi)
7126 const size_t vi = order(
oi);
7127 const size_t prev = (vi + n - 1) % n;
7128 const size_t e_prev = prev;
7129 const size_t e_curr = vi;
7210 if (faces.
size() == 0)
7214 for (
size_t i = 0; i < n; ++i)
7259 if (faces(
fi).
size() < 3)
7264 for (
size_t i = 0; i < faces(
fi).
size(); ++i)
7278 for (
size_t i = 0; i <
mono.size(); ++i)
7280 if (
piece.size() >= 3)
7290 result.
append(it.get_curr());
7307 const size_t n = v.
size();
7314 for (
size_t i = 1; i < n; ++i)
7316 if (v(i).get_y() > v(
top_idx).get_y()
or
7318 v(i).get_x() < v(
top_idx).get_x()))
7321 if (v(i).get_y() < v(
bot_idx).get_y()
or
7323 v(i).get_x() > v(
bot_idx).get_x()))
7331 const size_t next = (i + 1) % n;
7334 if (v(
next).get_y() > v(i).get_y())
7343 const size_t prev = (i + n - 1) % n;
7346 if (v(prev).get_y() > v(i).get_y())
7406 if (
verts.size() < 3)
7418 for (
size_t i = 1; i < v.
size(); ++i)
7419 if (v(i).get_y() < v(
bot).get_y()
or (v(i).get_y() == v(
bot).get_y()
and
7420 v(i).get_x() < v(
bot).get_x()))
7425 for (
size_t i = 0; i < v.
size(); ++i)
7433 const size_t j = (i + 1) % v.
size();
7434 return {v(j).get_x() - v(i).get_x(), v(j).get_y() - v(i).get_y()};
7470 const size_t np =
pv.size();
7471 const size_t nq =
qv.size();
7512 for (
size_t i = 0; i < result.
size(); ++i)
7515 clean.append(result(i));
7518 static_cast<void>(
clean.remove_last());
7521 for (
size_t i = 0; i <
clean.size(); ++i)
7524 if (
ret.size() >= 3)
7591 static constexpr double kEps = 1e-12;
7599 const double bx,
const double by)
noexcept
7622 (1.0 - t) *
ay + t *
by);
7629 for (
size_t i = 0; i < v.
size(); ++i)
7645 double best_p = -std::numeric_limits<double>::infinity();
7646 double best_q = -std::numeric_limits<double>::infinity();
7648 for (
size_t i = 0; i < p.
size(); ++i)
7660 for (
size_t i = 0; i < q.
size(); ++i)
7685 if (
sp.intersects_with(
itq.get_current_segment()))
7692 vp = it.get_current_vertex().to_point();
7699 vq = it.get_current_vertex().to_point();
7707 [[
nodiscard]]
static ClosestSimplexResult
7729 const double abx = b.
vx - a.
vx;
7730 const double aby = b.
vy - a.
vy;
7735 out.reduced_simplex = {a};
7738 out.witness_a = a.
a;
7739 out.witness_b = a.
b;
7745 if (t < 0.0) t = 0.0;
7746 if (t > 1.0) t = 1.0;
7750 out.reduced_simplex = {a};
7753 out.witness_a = a.
a;
7754 out.witness_b = a.
b;
7756 else if (t >= 1.0 -
kEps)
7758 out.reduced_simplex = {b};
7761 out.witness_a = b.
a;
7762 out.witness_b = b.
b;
7766 out.reduced_simplex = {a, b};
7782 const double abx = b.
vx - a.
vx;
7783 const double aby = b.
vy - a.
vy;
7784 const double acx = c.
vx - a.
vx;
7785 const double acy = c.
vy - a.
vy;
7786 const double apx = -a.
vx;
7787 const double apy = -a.
vy;
7791 if (d1 <= 0.0
and d2 <= 0.0)
7793 out.reduced_simplex = {a};
7796 out.witness_a = a.
a;
7797 out.witness_b = a.
b;
7802 const double bpx = -b.
vx;
7803 const double bpy = -b.
vy;
7806 if (d3 >= 0.0
and d4 <= d3)
7808 out.reduced_simplex = {b};
7811 out.witness_a = b.
a;
7812 out.witness_b = b.
b;
7817 const double vc = d1 * d4 - d3 * d2;
7820 const double t = d1 / (d1 - d3);
7821 out.reduced_simplex = {a, b};
7830 const double cpx = -c.
vx;
7831 const double cpy = -c.
vy;
7836 out.reduced_simplex = {c};
7839 out.witness_a = c.
a;
7840 out.witness_b = c.
b;
7845 const double vb = d5 * d2 - d1 *
d6;
7848 const double t = d2 / (d2 -
d6);
7849 out.reduced_simplex = {a, c};
7858 const double va = d3 *
d6 - d5 * d4;
7861 const double t = (d4 - d3) / (d4 - d3 + (d5 -
d6));
7862 const double bcx = c.
vx - b.
vx;
7863 const double bcy = c.
vy - b.
vy;
7864 out.reduced_simplex = {b, c};
7877 out.reduced_simplex = {a, b, c};
7880 out.witness_a = a.
a;
7881 out.witness_b = a.
b;
7882 out.contains_origin =
true;
7889 const double wc =
vc *
inv;
7890 out.reduced_simplex = {a, b, c};
7905 out.contains_origin =
true;
7920 out.distance_squared = 0;
7922 out.closest_on_first =
pv(0);
7923 out.closest_on_second =
pv(0);
7924 out.intersects =
true;
7932 for (
size_t i = 0; i <
pv.size(); ++i)
7935 for (
size_t j = 0; j <
qv.size(); ++j)
7938 const Point & q1 =
qv((j + 1) %
qv.size());
7951 for (
size_t i = 0; i <
qv.size(); ++i)
7954 for (
size_t j = 0; j <
pv.size(); ++j)
7957 const Point & p1 =
pv((j + 1) %
pv.size());
7999 <<
"First polygon must be convex";
8001 <<
"Second polygon must be convex";
8009 r.closest_on_first =
pv(0);
8010 r.closest_on_second =
pv(0);
8011 r.intersects =
true;
8025 std::vector<SupportPoint>
simplex;
8031 double prev_d2 = std::numeric_limits<double>::infinity();
8052 r.closest_on_first =
step.witness_a;
8053 r.closest_on_second =
step.witness_b;
8054 r.intersects =
true;
8055 r.gjk_iterations =
iters + 1;
8125 for (
size_t i = 0; i < points.
size(); ++i)
8126 if (
uniq.is_empty()
or uniq.get_last() != points(i))
8127 uniq.append(points(i));
8165 :
tree(xmin, ymin, xmax, ymax),
8166 bounds_(xmin, ymin, xmax, ymax)
8227 template <
typename Op>
8248 auto recurse = [&](
const auto & self,
8252 const bool split_on_x) ->
void
8276 part.split_on_x = split_on_x;
8290 for (
size_t i = 0; i <
mid; ++i)
8292 for (
size_t i =
mid; i <
sorted.size(); ++i)
8297 Rectangle left(region.get_xmin(), region.get_ymin(),
8298 pivot.get_x(), region.get_ymax());
8300 region.get_xmax(), region.get_ymax());
8301 self(self,
left_pts, left, depth + 1,
false);
8302 self(self,
right_pts, right, depth + 1,
false);
8306 Rectangle bottom(region.get_xmin(), region.get_ymin(),
8307 region.get_xmax(),
pivot.get_y());
8309 region.get_xmax(), region.get_ymax());
8310 self(self,
left_pts, bottom, depth + 1,
true);
8311 self(self,
right_pts, top, depth + 1,
true);
8362 for (
size_t i = 0; i < face.
size(); ++i)
8372 const size_t tmp = u;
8376 return v == u + 1
or (u == 0
and v == n - 1);
8385 const size_t n1 = f1.
size();
8386 const size_t n2 = f2.
size();
8392 if ((
pu1 + 1) % n1 !=
pv1)
8394 const size_t tmp = u;
8401 const size_t prev_u = f1((
pu1 + n1 - 1) % n1);
8402 const size_t next_v = f1((
pv1 + 1) % n1);
8407 const size_t next_u = f2((
pu2 + 1) % n2);
8408 const size_t prev_v = f2((
pv2 + n2 - 1) % n2);
8423 const size_t n1 = f1.
size();
8424 const size_t n2 = f2.
size();
8429 if ((
pu1 + 1) % n1 !=
pv1)
8431 const size_t tmp = u;
8443 for (
size_t k = 0;
k < n1 - 1; ++
k)
8447 for (
size_t k = 0;
k < n2 - 1; ++
k)
8471 const size_t n =
pts.size();
8488 const Triangle & t = it.get_curr();
8490 for (
size_t i = 0; i < n; ++i)
8501 const size_t tmp =
i1;
8510 faces.
append(std::move(face));
8520 const auto f1 = faces(
fi);
8523 const size_t u = f1(
k);
8524 const size_t v = f1((
k + 1) % f1.size());
8531 const auto f2 = faces(
fj);
8542 for (
size_t i = 0; i < faces.
size(); ++i)
8562 for (
size_t k = 0;
k < faces(
fi).
size(); ++
k)
8567 result.
append(std::move(p));
8650 size_t lo = 0, hi = arr.
size();
8652 if (
const size_t mid = lo + (hi - lo) / 2; arr(
mid).get_y() < ymin)
8663 size_t lo = 0, hi = arr.
size();
8665 if (
const size_t mid = lo + (hi - lo) / 2; arr(
mid).get_y() <= ymax)
8675 size_t lo = 0, hi =
n_;
8677 if (
const size_t mid = lo + (hi - lo) / 2;
pts_(
mid).get_x() <
xval)
8687 size_t lo = 0, hi =
n_;
8689 if (
const size_t mid = lo + (hi - lo) / 2;
pts_(
mid).get_x() <=
xval)
8696 void build_node(
const size_t node,
const size_t lo,
const size_t hi)
8704 const size_t mid = lo + (hi - lo) / 2;
8709 const auto & left =
tree_(2 * node).y_sorted;
8710 const auto & right =
tree_(2 * node + 1).y_sorted;
8711 auto & merged =
tree_(node).y_sorted;
8712 merged.reserve(left.size() + right.size());
8714 size_t i = 0, j = 0;
8715 while (i < left.size()
and j < right.size())
8716 if (left(i).get_y() < right(j).get_y()
or
8717 (left(i).get_y() == right(j).get_y() &&
8718 left(i).get_x() <= right(j).get_x()))
8719 merged.append(left(i++));
8721 merged.append(right(j++));
8722 while (i < left.size())
8723 merged.append(left(i++));
8724 while (j < right.size())
8725 merged.append(right(j++));
8729 const size_t qlo,
const size_t qhi,
8738 const auto &
ys =
tree_(node).y_sorted;
8741 for (
size_t i = from; i < to; ++i)
8746 const size_t mid = lo + (hi - lo) / 2;
8757 pts_.append(it.get_curr());
8775 for (
size_t i = 0; i <
tree_sz; ++i)
8816 const size_t hi) ->
size_t
8832 const size_t mid = lo + (hi - lo) / 2;
8879 if (v.
size() < 3)
return false;
8910 nx = -dy *
md / len;
8916 ny = -dx *
md / len;
8969 const size_t n = v.
size();
8974 for (
size_t i = 0; i < n; ++i)
8976 const size_t j = (i + 1) % n;
9005 const size_t n = v.
size();
9011 for (
size_t i = 0; i < n; ++i)
9014 offset_edge(v(i), v((i + 1) % n), distance,
false, a, b);
9021 for (
size_t i = 0; i < n; ++i)
9023 const size_t j = (i + 1) % n;
9035 if (result.
size() >= 3)
9157 for (
size_t i = 1; i <
r.polygons.size(); ++i)
9163 return r.polygons(
best);
9199 return dx * dx + dy * dy;
9206 const size_t n = v.
size();
9207 const bool inward = distance < 0;
9209 if (inward) d = -distance;
9215 for (
size_t i = 0; i < n; ++i)
9218 offset_edge(v(i), v((i + 1) % n), d, inward, a, b);
9234 for (
size_t i = 0; i < n; ++i)
9236 const size_t j = (i + 1) % n;
9270 for (
size_t i = 0; i < n; ++i)
9272 const size_t j = (i + 1) % n;
9320 return (
I.get_x() -
P.get_x()) / dx;
9321 return (
I.get_y() -
P.get_y()) / (
Q.get_y() -
P.get_y());
9328 const size_t n =
raw.size();
9331 for (
size_t i = 0; i < n; ++i)
9334 const Point & a1 =
raw((i + 1) % n);
9337 for (
size_t j = i + 2; j < n; ++j)
9340 if (j == (i + n - 1) % n)
9342 if (i == 0 && j == n - 1)
9346 const Point & b1 =
raw((j + 1) % n);
9366 for (
size_t i = 1; i < idx.
size(); ++i)
9368 const size_t val = idx(i);
9371 while (j > 0 &&
keys(idx(j - 1)) > key)
9373 idx(j) = idx(j - 1);
9385 const size_t n =
raw.size();
9393 for (
size_t k = 0;
k <
m; ++
k)
9402 for (
size_t i = 0; i < n; ++i)
9412 for (
size_t k = 0;
k <
m; ++
k)
9419 else if (
isects(
k).edge_j == i)
9435 for (
size_t p = 0; p < order.
size(); ++p)
9448 iv.is_intersection =
true;
9455 for (
size_t k = 0;
k <
m; ++
k)
9468 for (
size_t i = 0; i <
pts.size(); ++i)
9470 if (
pts.size() >= 3)
9486 const size_t safety =
N * 2 + 10;
9488 for (
size_t s = 0; s <
N; ++s)
9490 if (
not aug(s).is_intersection ||
aug(s).visited)
9498 aug(s).visited =
true;
9500 aug(
aug(s).partner).visited =
true;
9504 const size_t entry = idx;
9510 idx = (idx + 1) %
N;
9512 while (
not aug(idx).is_intersection)
9515 idx = (idx + 1) %
N;
9521 aug(idx).visited =
true;
9524 aug(
aug(idx).partner).visited =
true;
9525 idx =
aug(idx).partner;
9531 while (idx != entry);
9534 if (boundary.
size() >= 3)
9555 if (sa > 0 &&
raw.size() >= 3)
9599 if (dx > 0
and dy >= 0)
return 0;
9601 if (dx < 0
and dy <= 0)
return 2;
9683 if (a.
t != b.
t)
return a.
t < b.
t;
9697 const Point & query)
9700 for (
size_t i = 0; i < n; ++i)
9732 if (
tree_.is_empty())
9734 return tree_.min().edge;
9754 const Point & query)
const
9759 <<
"Query point must be strictly inside the polygon";
9765 const size_t n =
verts.size();
9770 for (
size_t i = 0; i < n; ++i)
9783 for (
size_t ei = 0;
ei < n; ++
ei)
9800 for (
size_t i = 0; i < n; ++i)
9802 for (
size_t i = 0; i < n; ++i)
9805 for (
size_t ei = 0;
ei < n; ++
ei)
9807 const size_t a =
ei, b = (
ei + 1) % n;
9826 for (
size_t ei = 0;
ei < n; ++
ei)
9833 for (
size_t oi = 0;
oi < n; ++
oi)
9835 const size_t vi = order(
oi);
9837 const Point dir = v;
9844 status.
erase(it.get_curr());
9848 status.
insert(it.get_curr(), dir);
9886 for (
size_t i = 0; i <
vis_pts.size(); ++i)
9892 if (result.
size() >= 3)
9897 rv.
append(it.get_current_vertex());
9898 if (
rv.size() >= 2
and rv(0) ==
rv(
rv.size() - 1))
9901 for (
size_t i = 0; i + 1 <
rv.size(); ++i)
9954 for (
size_t i = 0; i <
pts.size(); ++i)
9955 if (
pts(i) == p)
return i;
9967 const Triangle & t = it.get_curr();
9985 bool operator()(
const EdgeKey & a,
const EdgeKey & b)
const
10005 size_t u, v, tri, local;
10011 for (
int e = 0; e < 3; ++e)
10013 size_t u =
tris(
ti).v[(e + 1) % 3];
10014 size_t v =
tris(
ti).v[(e + 2) % 3];
10017 const size_t tmp = u;
10031 return a.tri < b.tri;
10034 for (
size_t i = 0; i + 1 < edges.
size(); ++i)
10035 if (edges(i).u == edges(i + 1).u
and edges(i).v == edges(i + 1).v)
10037 const size_t t1 = edges(i).tri,
l1 = edges(i).local;
10038 const size_t t2 = edges(i + 1).tri,
l2 = edges(i + 1).local;
10067 for (
size_t i = 0; i <
tris.size(); ++i)
10074 const size_t src,
const size_t dst)
10085 for (
size_t i = 0; i <
tris.size(); ++i)
10112 for (
size_t i = 0; i < path.
size() / 2; ++i)
10114 const size_t tmp = path(i);
10115 path(i) = path(path.
size() - 1 - i);
10116 path(path.
size() - 1 - i) =
tmp;
10141 const Point & source,
10142 const Point & target)
const
10147 <<
"Source must be inside the polygon";
10149 <<
"Target must be inside the polygon";
10152 if (source == target)
10160 const Segment seg(source, target);
10161 bool blocked =
false;
10228 for (
size_t i = 0; i + 1 <
sleeve.size(); ++i)
10230 size_t s0 = 0,
s1 = 0;
10233 bool found[3] = {
false,
false,
false};
10234 for (
int a = 0; a < 3; ++a)
10235 for (
unsigned long b:
tris(
sleeve(i + 1)).v)
10242 for (
int a = 0; a < 3; ++a)
10265 Point apex = source;
10268 size_t ai = 0,
li = 0,
ri = 0;
10272 for (
size_t i = 1; i < portals.
size(); ++i)
10274 const Point &
pr = portals(i).right;
10275 const Point &
pl = portals(i).left;
10405 size_t lo = 0, hi =
verts.size();
10420 if (dx > 0
and dy >= 0)
return 0;
10422 if (dx < 0
and dy <= 0)
return 2;
10435 if (q1 !=
q2)
return q1 <
q2;
10464 const size_t n = segments.
size();
10471 ret.faces.append(uf);
10483 for (
size_t i = 0; i < n; ++i)
10485 all_pts.append(segments(i).get_src_point());
10486 all_pts.append(segments(i).get_tgt_point());
10488 for (
size_t i = 0; i <
inters.size(); ++i)
10499 for (
size_t i = 0; i <
all_pts.size(); ++i)
10500 if (
ret.vertices.is_empty()
or ret.vertices.get_last() !=
all_pts(i))
10505 for (
size_t si = 0;
si < n; ++
si)
10507 const Point &
sp = segments(
si).get_src_point();
10508 const Point & tp = segments(
si).get_tgt_point();
10527 for (
size_t i = 0; i <
on_seg.size(); ++i)
10536 const auto &
verts =
ret.vertices;
10549 for (
size_t i = 0; i + 1 <
on_seg.size(); ++i)
10554 const size_t ne =
ret.edges.size();
10555 const size_t nhe = 2 *
ne;
10556 const size_t nv =
ret.vertices.size();
10562 uf.unbounded =
true;
10563 ret.faces.append(uf);
10571 for (
size_t i = 0; i <
ne; ++i)
10573 he.append(HalfEdge{
10574 ret.edges(i).src,
ret.edges(i).tgt,
10577 he.append(HalfEdge{
10578 ret.edges(i).tgt,
ret.edges(i).src,
10586 for (
size_t i = 0; i <
nv; ++i)
10588 for (
size_t h = 0;
h <
nhe; ++
h)
10592 const auto &
verts =
ret.vertices;
10593 for (
size_t v = 0; v <
nv; ++v)
10595 if (
inc(v).size() <= 1)
continue;
10599 -
verts(v).get_x();
10601 -
verts(v).get_y();
10603 -
verts(v).get_x();
10605 -
verts(v).get_y();
10619 for (
size_t v = 0; v <
nv; ++v)
10621 const auto & list =
inc(v);
10622 const size_t sz = list.size();
10623 if (sz == 0)
continue;
10625 for (
size_t i = 0; i < sz; ++i)
10646 const size_t prev = (i == 0) ? sz - 1 : i - 1;
10647 he(
he(list(i)).twin).next = list(prev);
10654 for (
size_t i = 0; i <
nhe; ++i)
10657 for (
size_t h = 0;
h <
nhe; ++
h)
10659 if (visited(
h))
continue;
10662 face.unbounded =
false;
10669 visited(
cur) =
true;
10682 if (signed_area < 0)
10683 face.unbounded =
true;
10685 ret.faces.append(face);
10694 for (
size_t i = 0; i <
ret.faces.size(); ++i)
10695 if (
ret.faces(i).unbounded)
10701 ret.faces(0).unbounded =
true;
10772 auto [sites, triangles] = delaunay(points);
10777 if (triangles.size() == 0)
10783 for (
size_t t = 0; t < triangles.size(); ++t)
10785 const auto & tri = triangles(t);
10786 const Point & a = sites(tri.i);
10787 const Point & b = sites(tri.j);
10788 const Point & c = sites(tri.k);
10802 ret.triangles.append(tri);
10817 return u ==
o.u
and v ==
o.v;
10822 return u <
o.u
or (u ==
o.u
and v <
o.v);
10826 auto make_key = [](
const size_t a,
const size_t b) -> EdgeKey
10828 return a < b ? EdgeKey{a, b} : EdgeKey{b, a};
10833 for (
size_t t = 0; t <
ret.triangles.size(); ++t)
10835 const auto & [i, j,
k] =
ret.triangles(t);
10846 for (
size_t i = 0; i <
all_edges.size();)
11004 const size_t n = sites.
size();
11005 if (n == 0)
return ret;
11012 for (
size_t i = 0; i < n; ++i)
11013 rt_input.append({sites(i).position, sites(i).weight});
11018 if (
rt.triangles.size() == 0)
return ret;
11024 for (
size_t ri = 0;
ri <
rt.sites.size(); ++
ri)
11027 for (
size_t oi = 0;
oi < n; ++
oi)
11028 if (
rt.sites(
ri).position == sites(
oi).position)
11039 for (
size_t t = 0; t <
rt.triangles.size(); ++t)
11041 const auto & tri =
rt.triangles(t);
11058 const size_t first,
const size_t last)
11060 if (last - first == 1)
11063 hull_sites.insert(rt_to_ours(edge_refs(first).u));
11064 hull_sites.insert(rt_to_ours(edge_refs(first).v));
11069 if (last - first >= 2)
11070 for (
size_t a = first; a + 1 < last; ++a)
11071 for (
size_t b = a + 1; b < last; ++b)
11078 pe.unbounded =
false;
11080 ret.edges.append(
pe);
11085 ret.cells.reserve(n);
11086 for (
size_t s = 0; s < n; ++s)
11089 cell.site_index = s;
11091 cell.weight = sites(s).weight;
11092 cell.bounded =
hull_sites.search(s) ==
nullptr;
11095 for (
size_t t = 0; t <
rt.triangles.size(); ++t)
11097 const auto & [i, j,
k] =
rt.triangles(t);
11100 cell.vertices.append(
pcenters(t));
11103 ret.cells.append(cell);
11181 for (
size_t i = 0; i <= n; ++i)
11195 for (
size_t i = 0; i <= n; ++i)
11217 s * a.
get_x() + t * b.get_x(),
11218 s * a.
get_y() + t * b.get_y()
11222 const Point q0 = lerp(p0, p1);
11223 const Point q1 = lerp(p1, p2);
11224 const Point q2 = lerp(p2, p3);
11248 auto update = [&](
const Point & p)
11250 if (p.get_x() <
mnx)
mnx = p.get_x();
11251 if (p.get_x() >
mxx)
mxx = p.get_x();
11252 if (p.get_y() <
mny)
mny = p.get_y();
11253 if (p.get_y() >
mxy)
mxy = p.get_y();
11264 const Point & p2,
const size_t n)
11268 for (
size_t i = 0; i <
pts.size(); ++i)
11276 const Point & p2,
const Point & p3,
const size_t n)
11280 for (
size_t i = 0; i <
pts.size(); ++i)
11425 for (
size_t i = 0, j = v.
size() - 1; i < j; ++i, --j)
11437 for (
size_t i = 0; i <
pts.size(); ++i)
11439 if (
pts.size() >= 3)
11450 return (
I.get_x() -
P.get_x()) / dx;
11451 return (
I.get_y() -
P.get_y()) / (
Q.get_y() -
P.get_y());
11459 const size_t n = poly.
size();
11460 for (
size_t i = 0; i < n; ++i)
11462 const Point & a = poly(i);
11463 const Point & b = poly((i + 1) % n);
11487 for (
size_t i = 1; i < idx.
size(); ++i)
11489 const size_t val = idx(i);
11492 while (j > 0 &&
keys(idx(j - 1)) > key)
11494 idx(j) = idx(j - 1);
11520 const size_t na =
va.size();
11521 const size_t nb =
vb.size();
11527 for (
size_t i = 0; i <
na; ++i)
11533 for (
size_t j = 0; j <
nb; ++j)
11557 for (
size_t k = 0;
k < pairs.
size(); ++
k)
11565 for (
size_t i = 0; i <
na; ++i)
11574 for (
size_t k = 0;
k < pairs.
size(); ++
k)
11575 if (pairs(
k).edge_a == i)
11583 for (
size_t p = 0; p <
edge_pairs.size(); ++p)
11587 for (
size_t p = 0; p < order.
size(); ++p)
11592 in.
pt = pairs(
k).pt;
11593 in.intersect =
true;
11594 in.alpha = pairs(
k).alpha_a;
11600 for (
size_t j = 0; j <
nb; ++j)
11608 for (
size_t k = 0;
k < pairs.
size(); ++
k)
11609 if (pairs(
k).edge_b == j)
11616 for (
size_t p = 0; p <
edge_pairs.size(); ++p)
11620 for (
size_t p = 0; p < order.
size(); ++p)
11625 in.
pt = pairs(
k).pt;
11626 in.intersect =
true;
11627 in.alpha = pairs(
k).alpha_b;
11633 for (
size_t k = 0;
k < pairs.
size(); ++
k)
11645 for (
size_t i = 0; i <
list_a.size(); ++i)
11646 if (
list_a(i).intersect)
11656 for (
size_t j = 0; j <
list_b.size(); ++j)
11657 if (
list_b(j).intersect)
11669 for (
size_t s = 0; s <
list_a.size(); ++s)
11700 idx = (idx + 1) %
cur_list.size();
11704 idx = (idx + 1) %
cur_list.size();
11721 while (! (
on_a && idx == s));
11724 if (boundary.
size() >= 3)
11773 if (!
res.is_empty())
11810 if (!
res.is_empty())
11864 std::ostringstream
os;
11885 return "LINESTRING (" +
11897 auto pt = [](
const Point & p)
11899 return dbl(p.get_x()) +
" " +
dbl(p.get_y());
11901 return "POLYGON ((" + pt(t.
get_p1()) +
", " + pt(t.
get_p2()) +
11910 return "POLYGON ((" +
11911 dbl(
r.get_xmin()) +
" " +
dbl(
r.get_ymin()) +
", " +
11912 dbl(
r.get_xmax()) +
" " +
dbl(
r.get_ymin()) +
", " +
11913 dbl(
r.get_xmax()) +
" " +
dbl(
r.get_ymax()) +
", " +
11914 dbl(
r.get_xmin()) +
" " +
dbl(
r.get_ymax()) +
", " +
11915 dbl(
r.get_xmin()) +
" " +
dbl(
r.get_ymin()) +
"))";
11923 std::ostringstream
os;
11924 os << std::setprecision(15) <<
"POLYGON ((";
11929 const Point & v = it.get_current_vertex();
11962 return R
"({"type":"Point","coordinates":[)" +
11971 return R
"({"type":"LineString","coordinates":[[)" +
11981 auto pt = [](
const Point & p)
11983 return "[" +
dbl(p.get_x()) +
"," +
dbl(p.get_y()) +
"]";
11985 return R
"({"type":"Polygon","coordinates":[[)" +
11993 std::ostringstream
os;
11994 os << std::setprecision(15);
11995 os << R
"({"type":"Polygon","coordinates":[[)";
12000 const Point & v = it.get_current_vertex();
12020 return R
"({"type":"Point","coordinates":[)" +
12131 return {xmin, ymin, xmax, ymax};
12155 const bool split_x)
const
12158 entries_(entry_idx).bbox.get_xmin() +
entries_(entry_idx).bbox.get_xmax() :
12159 entries_(entry_idx).bbox.get_ymin() +
entries_(entry_idx).bbox.get_ymax();
12163 const size_t hi,
const bool split_x)
const
12165 for (
size_t i = lo + 1; i < hi; ++i)
12167 const size_t key = idx(i);
12172 idx(j) = idx(j - 1);
12182 const size_t kth,
const bool split_x)
const
12184 while (hi - lo > 16)
12197 const size_t tmp = idx(
lt);
12203 else if (val >
pivot)
12206 const size_t tmp = idx(
gt);
12216 else if (
kth >=
gt)
12240 for (
size_t i = lo + 1; i < hi; ++i)
12246 const bool split_x = dx >= dy;
12248 const size_t mid = lo + (hi - lo) / 2;
12251 const size_t left =
build(idx, lo,
mid);
12252 const size_t right =
build(idx,
mid, hi);
12320 for (
size_t i = 0; i <
entries_.size(); ++i)
12374 auto dfs = [&](
const auto & self,
const size_t ni,
const size_t depth) ->
size_t
12391 dn.
left = self(self, n.
left, depth + 1);
12427 if (last <= first + 1)
12434 for (
size_t i = first + 1; i < last; ++i)
12472 for (
size_t i = 0; i < n; ++i)
12473 keep.append(
false);
12475 keep(n - 1) =
true;
12480 for (
size_t i = 0; i < n; ++i)
12495 arr.
append(it.get_curr());
12496 return (*
this)(arr, epsilon);
12514 const size_t n =
verts.size();
12521 for (
size_t i = 1; i < n; ++i)
12534 for (
size_t i = 0; i <=
far_idx; ++i)
12538 for (
size_t i =
far_idx; i < n; ++i)
12542 auto s1 = (*this)(
chain1, epsilon);
12543 auto s2 = (*this)(
chain2, epsilon);
12547 for (
size_t i = 0; i <
s1.size(); ++i)
12549 for (
size_t i = 1; i + 1 <
s2.size(); ++i)
12550 merged.append(
s2(i));
12552 if (merged.size() < 3)
12556 for (
size_t i = 0; i < merged.size(); ++i)
12610 for (
size_t i = 0; i < n; ++i)
12612 prev.append(i == 0 ? n : i - 1);
12613 next.append(i == n - 1 ? n : i + 1);
12622 for (
size_t i = 0; i < n; ++i)
12625 for (
size_t i = 1; i + 1 < n; ++i)
12634 auto [idx, area] = heap.
get();
12636 if (
not alive(idx))
12638 if (area !=
areas(idx))
12644 alive(idx) =
false;
12645 const size_t p = prev(idx);
12646 const size_t nx =
next(idx);
12654 if (p < n
and p != 0
and prev(p) < n)
12669 for (
size_t i = 0; i < n; ++i)
12701 arr.
append(it.get_curr());
12719 const size_t n =
verts.size();
12726 for (
size_t i = 0; i < n; ++i)
12728 prev.append(i == 0 ? n - 1 : i - 1);
12729 next.append(i == n - 1 ? 0 : i + 1);
12736 for (
size_t i = 0; i < n; ++i)
12743 size_t remaining = n;
12746 auto [idx, area] = heap.
get();
12748 if (
not alive(idx))
12750 if (area !=
areas(idx))
12755 alive(idx) =
false;
12757 const size_t p = prev(idx);
12758 const size_t nx =
next(idx);
12772 for (
size_t i = 0; i < n; ++i)
12776 if (result.
size() < 3)
12780 for (
size_t i = 0; i < result.
size(); ++i)
12810 const size_t n =
pts.size();
12813 for (
size_t i = 0; i + 1 < n; ++i)
12815 const auto & p =
pts(i);
12816 const auto & q =
pts(i + 1);
12819 p.get_y() *
one_r + q.get_y() *
r));
12821 p.get_y() *
r + q.get_y() *
one_r));
12830 const size_t n =
pts.size();
12833 for (
size_t i = 0; i < n; ++i)
12835 const auto & p =
pts(i);
12836 const auto & q =
pts((i + 1) % n);
12838 p.get_y() *
one_r + q.get_y() *
r));
12840 p.get_y() *
r + q.get_y() *
one_r));
12862 ratio >=
Geom_Number(1, 2)) <<
"ratio must be in (0, 0.5)";
12865 for (
size_t k = 0;
k < iterations; ++
k)
12878 arr.
append(it.get_curr());
12879 return (*
this)(arr, iterations, ratio);
12898 ratio >=
Geom_Number(1, 2)) <<
"ratio must be in (0, 0.5)";
12901 for (
size_t k = 0;
k < iterations; ++
k)
12905 for (
size_t i = 0; i <
cur.size(); ++i)
13064 const auto & [type, index, left, right] =
dag(
ni);
13073 ni = p <
xp ? left : right;
13092 else if (cross > 0)
13114 if (p ==
points(i))
return true;
13158 const auto & [type, index, left, right] = dag(
ni);
13167 ni = p <
xp ? left : right;
13173 segs(index).get_tgt_point(), p);
13176 else if (cross < 0)
13186 const size_t top,
const size_t bottom,
13187 const size_t leftp,
const size_t rightp)
13189 const size_t ti =
traps.size();
13190 const size_t di = dag.
size();
13192 top, bottom, leftp, rightp,
13201 const NodeType type,
const size_t index,
13202 const size_t left,
const size_t right)
13228 const size_t seg_idx,
const size_t start_trap)
13243 if (
rp < rightp
or rp == rightp)
13272 const size_t seg_idx,
const size_t trap_idx)
13281 for (
size_t i = 0; i <
pts.size(); ++i)
13412 const size_t seg_idx,
13420 for (
size_t i = 0; i <
pts.size(); ++i)
13439 for (
size_t ci = 0;
ci <
nc; ++
ci)
13446 const bool is_last = (
ci ==
nc - 1);
13668 ret.segments.reserve(n + 2);
13669 for (
size_t i = 0; i < n; ++i)
13684 for (
size_t i = 0; i < n; ++i)
13686 all_pts.append(
ret.segments(i).get_src_point());
13687 all_pts.append(
ret.segments(i).get_tgt_point());
13692 for (
size_t i = 0; i <
all_pts.size(); ++i)
13694 if (
ret.points.is_empty()
or
13700 ret.num_input_points =
ret.points.size();
13703 if (
ret.points.is_empty())
13716 for (
size_t i = 1; i <
ret.points.size(); ++i)
13718 if (
ret.points(i).get_x() <
mnx)
mnx =
ret.points(i).get_x();
13719 if (
ret.points(i).get_x() >
mxx)
mxx =
ret.points(i).get_x();
13720 if (
ret.points(i).get_y() <
mny)
mny =
ret.points(i).get_y();
13721 if (
ret.points(i).get_y() >
mxy)
mxy =
ret.points(i).get_y();
13734 const size_t bl_idx =
ret.points.size();
13735 ret.points.append(
bl);
13736 ret.points.append(
br);
13737 ret.points.append(
tl);
13738 const size_t tr_idx =
ret.points.size();
13739 ret.points.append(
tr);
13754 if (n == 0)
return ret;
13758 for (
size_t i = 0; i < n; ++i)
13763 std::random_device rd;
13764 std::mt19937
gen(rd());
13765 for (
size_t i = n - 1; i > 0; --i)
13767 std::uniform_int_distribution<size_t>
dis(0, i);
13768 const size_t j =
dis(
gen);
13769 const size_t tmp = order(i);
13770 order(i) = order(j);
13776 for (
size_t oi = 0;
oi < n; ++
oi)
13778 const size_t si = order(
oi);
13788 ret.points,
ret.segments,
ret.trapezoids,
si, start);
13792 ret.trapezoids,
ret.dag,
13796 ret.trapezoids,
ret.dag,
13817 return (*
this)(
segs);
13829 for (
size_t i = 0; i < polygons.
size(); ++i)
13831 const Polygon & poly = polygons(i);
13835 segs.append(it.get_current_segment());
13837 return (*
this)(
segs);
13845#if __cplusplus >= 202002L
13847# include <concepts>
13870 template <
typename T>
13873 {
T(i) } -> std::convertible_to<T>;
13874 { a + b } -> std::convertible_to<T>;
13875 { a - b } -> std::convertible_to<T>;
13876 { a * b } -> std::convertible_to<T>;
13877 { a / b } -> std::convertible_to<T>;
13878 { a == b } -> std::convertible_to<bool>;
13879 { a != b } -> std::convertible_to<bool>;
13880 { a < b } -> std::convertible_to<bool>;
13881 { a <= b } -> std::convertible_to<bool>;
13882 { a > b } -> std::convertible_to<bool>;
13883 { a >= b } -> std::convertible_to<bool>;
13884 { -a } -> std::convertible_to<T>;
13889 "Geom_Number must satisfy GeomNumberType");
13891 "double must satisfy GeomNumberType");
13893 "long long must satisfy GeomNumberType");
Exception handling system with formatted messages for Aleph-w.
#define ah_domain_error_if(C)
Throws std::domain_error if condition holds.
bool operator==(const Time &l, const Time &r)
bool operator<(const Time &l, const Time &r)
Deduplicate sequential Aleph containers in-place.
Axis-aligned bounding box tree for spatial queries.
Rectangle root_bbox() const
Return the root bounding box (union of all entries).
Array< Entry > entries_
Original entries stored in the tree.
size_t build(Array< size_t > &idx, const size_t lo, const size_t hi)
void select_kth_by_center(Array< size_t > &idx, size_t lo, size_t hi, const size_t kth, const bool split_x) const
static bool box_contains_point(const Rectangle &r, const Point &p)
Checks if a rectangle contains a point.
static Rectangle union_bbox(const Rectangle &a, const Rectangle &b)
Computes the union of two axis-aligned bounding boxes.
void query_point_impl(const size_t ni, const Point &p, Array< size_t > &results) const
Geom_Number center_coord(const size_t entry_idx, const bool split_x) const
bool is_empty() const
Whether the tree is empty.
size_t size() const
Number of entries.
void insertion_sort_by_center(Array< size_t > &idx, const size_t lo, const size_t hi, const bool split_x) const
static constexpr size_t NONE
void build(const Array< Entry > &entries)
Build the AABB tree from an array of entries.
DebugSnapshot debug_snapshot() const
Return the full tree structure for visualization/debug.
static bool boxes_overlap(const Rectangle &a, const Rectangle &b)
Checks if two rectangles overlap.
Array< size_t > query(const Rectangle &query) const
Find all entries whose bounding box overlaps the query rectangle.
Array< size_t > query_point(const Point &p) const
Find all entries whose bounding box contains the query point.
Array< Node > nodes_
Internal array storing tree nodes.
void query_impl(const size_t ni, const Rectangle &q, Array< size_t > &results) const
Alpha shape of a point set.
Result operator()(const DynList< Point > &points, const Geom_Number &alpha_squared) const
Compute the α-shape.
Andrew's monotonic chain convex hull algorithm.
Polygon operator()(const DynList< Point > &point_set) const
Compute the convex hull of a point set.
static Array< Point > sorted_unique_points(const DynList< Point > &point_set)
Extracts points from a list, sorts them and removes duplicates.
static Geom_Number turn(const Point &a, const Point &b, const Point &c)
Computes the turn orientation of triple (a, b, c).
Simple dynamic array with automatic resizing and functional operations.
constexpr size_t size() const noexcept
Return the number of elements stored in the stack.
void empty() noexcept
Empties the container.
constexpr bool is_empty() const noexcept
Checks if the container is empty.
T & insert(const T &data)
insert a copy of data at the beginning of the array.
T & append(const T &data)
Append a copy of data
T & get_last() noexcept
return a modifiable reference to the last element.
void reserve(size_t cap)
Reserves cap cells into the array.
Quadratic and cubic Bézier curves with exact rational arithmetic.
static Polygon approximate_quadratic(const Point &p0, const Point &p1, const Point &p2, const size_t n)
Approximate a quadratic Bézier as a polyline (polygon without closing).
static Array< Point > sample_cubic(const Point &p0, const Point &p1, const Point &p2, const Point &p3, const size_t n)
Sample a cubic Bézier into n+1 points.
static Rectangle control_bbox(const Point &p0, const Point &p1, const Point &p2, const Point &p3)
Compute the bounding box of a cubic Bézier's control polygon.
static Array< Point > sample_quadratic(const Point &p0, const Point &p1, const Point &p2, const size_t n)
Sample a quadratic Bézier into n+1 points (t = 0, 1/n, ..., 1).
static SplitResult split_cubic(const Point &p0, const Point &p1, const Point &p2, const Point &p3, const Geom_Number &t)
static Point quadratic(const Point &p0, const Point &p1, const Point &p2, const Geom_Number &t)
Evaluate a quadratic Bézier at parameter t.
static Polygon approximate_cubic(const Point &p0, const Point &p1, const Point &p2, const Point &p3, const size_t n)
Approximate a cubic Bézier as a polyline.
static Point cubic(const Point &p0, const Point &p1, const Point &p2, const Point &p3, const Geom_Number &t)
Evaluate a cubic Bézier at parameter t.
Boolean operations on simple polygons (union, intersection, difference) using the Greiner-Hormann alg...
static void sort_indices_by_alpha(Array< size_t > &idx, const Array< Geom_Number > &keys)
Insertion-sort indices by a key array.
static Array< Point > extract(const Polygon &p)
Extract vertices from a polygon into an Array.
static void reverse_array(Array< Point > &v)
Reverse a vertex array in place.
static Geom_Number edge_param(const Point &P, const Point &Q, const Point &I)
Compute parameter t of point I along directed edge P→Q.
static Array< Polygon > compute_union(const Polygon &a, const Polygon &b)
Array< Polygon > difference(const Polygon &a, const Polygon &b) const
Convenience: difference (a minus b).
Array< Polygon > polygon_union(const Polygon &a, const Polygon &b) const
Convenience: union.
Array< Polygon > intersection(const Polygon &a, const Polygon &b) const
Convenience: intersection.
Array< Polygon > operator()(const Polygon &a, const Polygon &b, Op op) const
Compute a boolean operation on two simple polygons.
static bool point_inside_ccw(const Array< Point > &poly, const Point &p)
Point-in-polygon test on a CCW vertex array (winding number).
static Array< Polygon > greiner_hormann(const Array< Point > &va, const Array< Point > &vb, bool start_entering)
Run Greiner-Hormann and return result polygons.
Op
Type of Boolean operation to perform.
@ INTERSECTION
Find regions common to both polygons.
@ DIFFERENCE
Find regions in the first polygon not covered by the second.
@ UNION
Find regions covered by either polygon.
static Polygon build_poly(const Array< Point > &pts)
Build a polygon from a vertex array (bypasses colinearity merge).
static Array< Polygon > compute_difference(const Polygon &a, const Polygon &b)
static void ensure_ccw(Array< Point > &v)
Ensure vertices are in CCW order.
static Array< Polygon > compute_intersection(const Polygon &a, const Polygon &b)
Brute force convex hull algorithm.
static SegmentSet extreme_edges(const DynList< Point > &point_set)
Polygon operator()(const DynList< Point > &point_set) const
Compute the convex hull of a point set.
static bool are_all_points_on_left(const DynList< Point > &l, const Segment &s)
Checks if all points in a list lie to the left of a segment.
Chaikin corner-cutting subdivision for polygon smoothing.
static Array< Point > smooth_closed_once(const Array< Point > &pts, const Geom_Number &r)
Array< Point > operator()(const DynList< Point > &polyline, size_t iterations, const Geom_Number &ratio=Geom_Number(1, 4)) const
Overload accepting DynList.
static Polygon smooth_polygon(const Polygon &poly, size_t iterations, const Geom_Number &ratio=Geom_Number(1, 4))
Smooth a closed polygon.
Array< Point > operator()(const Array< Point > &polyline, size_t iterations, const Geom_Number &ratio=Geom_Number(1, 4)) const
Smooth an open polyline.
static Array< Point > smooth_open_once(const Array< Point > &pts, const Geom_Number &r)
Closest pair of points via divide and conquer.
static Result make_result(const Point &a, const Point &b)
Construct a Result object for a pair (a,b).
static Result brute_force(const Array< Point > &px, const size_t l, const size_t r)
Brute-force closest pair on a small subarray.
static Geom_Number dist2(const Point &a, const Point &b)
Squared Euclidean distance between two points.
Result operator()(const DynList< Point > &point_set) const
Compute the closest pair of points.
Segment closest_segment(DynList< Point > &point_set) const
Convenience wrapper returning the closest segment.
static Result recurse(const Array< Point > &px, const size_t l, const size_t r, Array< Point > &py)
Recursive divide-and-conquer step.
Constrained Delaunay Triangulation via Sloan's flip-based method.
Result operator()(const std::initializer_list< Point > points, const std::initializer_list< Segment > constraints) const
Convenience overload using initializer lists.
static void flip_edge(Array< Tri > &tris, const size_t tri_a, const size_t tri_b)
Flip the diagonal of the quad formed by tri_a and tri_b.
static bool lexicographic_less(const Point &p1, const Point &p2)
static Array< Point > merge_and_deduplicate(const DynList< Point > &points, const DynList< Segment > &constraints)
Merge input points, constraint endpoints, and constraint-constraint intersection points; then dedupli...
static size_t adj_of(const Tri &t, const size_t n)
Return the local index of the edge shared with neighbor n.
static constexpr size_t NONE
static bool edge_exists(const Array< Tri > &tris, const size_t u, const size_t v, size_t &out_tri, size_t &out_local)
Check if edge (u,v) already exists in the mesh.
static void build_adjacency(Array< Tri > &tris, const Array< IndexedTriangle > &dt_tris)
Build adjacency mesh from DT output.
static Array< IndexedEdge > map_constraints(const Array< Point > &sites, const DynList< Segment > &constraints)
Split constraints at interior collinear points.
static bool is_convex_quad(const Array< Point > &pts, const size_t a, const size_t b, const size_t c, const size_t d)
True if diagonal (a,d) in quad (a,b,d) + (a,d,c) forms a convex quad.
static void lawson_flip(const Array< Point > &pts, Array< Tri > &tris)
Lawson flip pass: restore Delaunay for non-constrained edges.
static void mark_constrained(Array< Tri > &tris, const size_t u, const size_t v)
Mark edge (u,v) as constrained in the mesh.
static void enforce_constraint(Array< Point > &pts, Array< Tri > &tris, const size_t u, const size_t v)
Enforce a single constraint edge (u,v) using Sloan's flip approach.
static DynList< Triangle > as_triangles(const Result &result)
Convert indexed triangulation to geometric triangles.
Result operator()(const DynList< Point > &points, const DynList< Segment > &constraints) const
Compute constrained Delaunay triangulation.
static size_t find_point_index(const Array< Point > &pts, const Point &p)
Find index of point p in sorted unique array, or NONE.
static size_t local_of(const Tri &t, const size_t id)
Return local index (0,1,2) of vertex id in triangle t, or NONE.
static Array< CrossingEdge > find_crossing_edges(const Array< Point > &pts, const Array< Tri > &tris, const size_t u, const size_t v)
static size_t edge_opposite(const Tri &t, const size_t a, const size_t b)
Return the local edge index for edge (a,b) in triangle t.
Decompose a simple polygon into convex parts using Hertel-Mehlhorn.
static Array< Point > extract_vertices(const Polygon &p)
static size_t find_pos(const Array< size_t > &face, size_t v)
Array< Polygon > operator()(const Polygon &poly) const
Decompose polygon poly into convex parts.
static bool can_merge(const Array< Point > &pts, const Array< size_t > &f1, const Array< size_t > &f2, size_t u, size_t v)
Check whether merging faces f1 and f2 across diagonal (u,v) is convex.
static Array< size_t > merge_faces(const Array< size_t > &f1, const Array< size_t > &f2, size_t u, size_t v)
Merge f1 and f2 by removing their shared edge (u,v).
static constexpr size_t NONE
static bool is_polygon_edge(size_t u, size_t v, size_t n)
Distance between two closed convex polygons using GJK.
Result operator()(const Polygon &p, const Polygon &q) const
Compute the distance between two closed convex polygons.
static Result exact_refine_disjoint(const Polygon &p, const Polygon &q, const Array< Point > &pv, const Array< Point > &qv)
static double norm2(const double x, const double y) noexcept
static ClosestSimplexResult closest_from_simplex(const std::vector< SupportPoint > &simplex)
static constexpr size_t kMaxIters
static double dot2(const double ax, const double ay, const double bx, const double by) noexcept
static Point point_from_double(const double x, const double y)
static bool polygons_intersect_or_contain(const Polygon &p, const Polygon &q)
static SupportPoint support(const Array< Point > &p, const Array< Point > &q, const double dx, const double dy)
static Point lerp_point(const Point &a, const Point &b, const double t)
static double to_double(const Geom_Number &v)
static Point centroid_of(const Array< Point > &v)
static constexpr double kEps
Basic exact intersection for closed convex polygons.
Polygon operator()(const Polygon &subject, const Polygon &clip) const
Intersect two closed convex polygons.
static Array< Point > normalize_vertices(const Array< Point > &pts)
static Point line_intersection(const Point &s, const Point &e, const Point &a, const Point &b)
static Geom_Number signed_double_area(const Array< Point > &verts)
static bool inside_half_plane(const Point &p, const Point &a, const Point &b, const bool clip_ccw)
static Array< Point > extract_vertices(const Polygon &poly)
static void push_clean(Array< Point > &out, const Point &p)
static Polygon build_polygon(const Array< Point > &pts)
static bool is_convex(const Array< Point > &verts)
Inward (erosion) and outward (dilation) offset of convex polygons.
static Point line_intersect(const Point &a1, const Point &a2, const Point &b1, const Point &b2)
Line-line intersection of (a1,a2) and (b1,b2).
static bool is_convex(const Array< Point > &v)
static void ensure_ccw(Array< Point > &v)
static Geom_Number signed_double_area(const Array< Point > &v)
static void offset_edge(const Point &a, const Point &b, const Geom_Number &d, const bool inward, Point &oa, Point &ob)
Compute two points on the offset line of edge (a, b) moved by distance d along its outward (or inward...
static Polygon outward(const Polygon &convex_poly, const Geom_Number &distance)
Outward offset (dilation) of a convex polygon.
static Array< Point > extract_verts(const Polygon &poly)
static Polygon inward(const Polygon &convex_poly, const Geom_Number &distance)
Inward offset (erosion) of a convex polygon.
Polygon triangulation using the ear-cutting algorithm.
static Geom_Number signed_double_area(const Polygon &p)
static Array< Point > extract_vertices(const Polygon &p)
static bool diagonal(const Polygon &p, const Vertex &a, const Vertex &b)
Check if segment (a, b) is a valid internal diagonal.
static EarsSet init_ears(const Polygon &p)
Initialize the set of ear vertices.
static void normalize_to_ccw(Polygon &p)
static bool in_cone(const Polygon &p, const Vertex &a, const Vertex &b)
Check if vertex b is inside the cone formed at vertex a.
static bool diagonalize(const Polygon &p, const Segment &s)
Check if a segment is a valid diagonal of the polygon.
DynList< Triangle > operator()(const Polygon &poly) const
Triangulate the polygon.
Exact Delaunay triangulation using the Bowyer-Watson incremental algorithm.
static bool all_collinear(const Array< Point > &pts)
Check if all points in an array are collinear.
Result operator()(const DynList< Point > &point_set) const
Compute Delaunay triangulation of a point set.
static bool point_in_circumcircle(const Array< Point > &pts, const size_t ia, const size_t ib, const size_t ic, const size_t ip)
Predicate to check if a point lies inside the circumcircle of a triangle.
static DynList< Triangle > as_triangles(const Result &result)
Convert indexed triangulation to geometric triangles.
static Array< Point > unique_points(const DynList< Point > &point_set)
Extracts unique points from a list and sorts them lexicographically.
static bool lexicographic_less(const Point &p1, const Point &p2)
Lexicographical comparison for points.
O(n log n) expected-time Delaunay's triangulation.
static bool in_cc(const Array< Point > &pts, const size_t ia, const size_t ib, const size_t ic, const size_t ip)
Standard in-circumcircle test using exact in_circle_determinant.
Result operator()(const std::initializer_list< Point > il) const
static constexpr size_t NONE
static void remap_adj(Array< Tri > &tris, const size_t ot, const size_t nt)
Update neighbor n of the old triangle ot to point to a new triangle nt.
static bool point_in_tri(const Array< Point > &pts, const Tri &t, const size_t pidx)
True if point p is inside or on triangle t (using orientation tests).
static size_t adj_of(const Tri &t, const size_t n)
Return the local index of the edge shared with neighbor n.
static size_t locate(const Array< Point > &pts, const Array< Tri > &tris, const Array< DagNode > &dag, const size_t pidx, const size_t root)
Locate the leaf triangle containing point pidx via a DAG walk.
static size_t local_of(const Tri &t, const size_t id)
Return local index (0,1,2) of vertex id in triangle t, or NONE.
Result operator()(const DynList< Point > &point_set) const
bool has_curr() const noexcept
Return true if the iterator has current item.
void append(Dlink *node) noexcept
Insert node before this.
void insert(Dlink *node) noexcept
Insert node after this.
Douglas-Peucker polyline simplification.
static void dp_recurse(const Array< Point > &pts, const size_t first, const size_t last, const Geom_Number &eps2, Array< bool > &keep)
Polygon simplify_polygon(const Polygon &poly, const Geom_Number &epsilon) const
Simplify a closed polygon.
Array< Point > operator()(const DynList< Point > &polyline, const Geom_Number &epsilon) const
Overload accepting DynList.
Array< Point > operator()(const Array< Point > &polyline, const Geom_Number &epsilon) const
Simplify an open polyline.
Dynamic heap of elements of type T ordered by a comparison functor.
T get()
Alias for getMin().
T & put(const T &item)
Synonym of insert().
Dynamic doubly linked list with O(1) size and bidirectional access.
T & append(const T &item)
Append a copied item at the end of the list.
Iterator on the items of list.
T & get_curr() const
Return the current item.
Dynamic singly linked list with functional programming support.
T & append(const T &item)
T & get_last() const
Return the last item of the list.
Dynamic set backed by balanced binary search trees with automatic memory management.
long position(const Key &key) const
Returns the infix (ordered) position of the key.
Key * insert(const Key &key)
Inserts a key into the dynamic set.
size_t remove(const Key &key)
Removes a key from the dynamic set.
Key * search(const Key &key) const
Find an element in the set.
Very simple queue implemented with a contiguous array.
size_t size() const noexcept
Return the number of elements stored in the stack.
bool is_empty() const noexcept
Return true if stack is empty.
T pop() noexcept
Pop by moving the top of stack.
T & top() noexcept
Return a modifiable reference to stack's top.
T & push(const T &data) noexcept(std::is_nothrow_copy_assignable_v< T >)
Push a copy of data
bool is_empty() const noexcept
Shared Bowyer-Watson core parameterized by in-circle predicate.
static Array< IndexedTriangle > triangulate(Array< Point > pts, const size_t n, InCirclePredicate point_in_circumcircle)
static void toggle_edge(EdgeSet &boundary, size_t u, size_t v)
Shared polygon helpers used across multiple geometry algorithms.
static Geom_Number area(const Array< Point > &verts)
Compute the absolute area of a vertex array.
static Geom_Number signed_double_area(const Polygon &poly)
Compute twice the signed area of a polygon.
static Geom_Number signed_area(const Polygon &poly)
Compute the actual signed area of a polygon.
static bool is_convex(const Array< Point > &verts)
Check if a vertex array forms a convex polygon.
static Geom_Number signed_double_area(const Array< Point > &verts)
Compute twice the signed area (shoelace formula without division).
static void ensure_ccw(Array< Point > &verts)
static Geom_Number signed_area(const Array< Point > &verts)
Compute the actual signed area of a vertex array.
static Geom_Number area(const Polygon &poly)
Compute the absolute area of a polygon.
static Array< Point > extract_vertices(const Polygon &poly)
Extract vertices from a polygon into an array for indexed access.
Serialization utilities for geometry objects.
static std::string to_geojson(const Triangle &t)
Triangle → GeoJSON Polygon.
static std::string to_wkt(const Point3D &p)
Converts a Point3D to WKT format: "POINT Z (x y z)".
static std::string to_geojson(const Polygon &poly)
Polygon → GeoJSON Polygon.
static std::string to_wkt(const Polygon &poly)
Converts a Polygon to WKT format: "POLYGON ((x1 y1, x2 y2, ..., x1 y1))".
static std::string to_wkt(const Rectangle &r)
Converts a Rectangle to WKT format: "POLYGON ((xmin ymin, xmax ymin, xmax ymax, xmin ymax,...
static std::string to_wkt(const Triangle &t)
Converts a Triangle to WKT format: "POLYGON ((x1 y1, x2 y2, x3 y3, x1 y1))".
static std::string dbl(const Geom_Number &n)
Formats a Geom_Number as a string with high precision.
static std::string to_geojson(const Segment &s)
Converts a Segment to a GeoJSON LineString.
static std::string to_wkt(const Segment &s)
Converts a Segment to WKT format: "LINESTRING (x1 y1, x2 y2)".
static std::string to_geojson(const Point &p)
Converts a Point to a GeoJSON geometry object.
static std::string to_wkt(const Point &p)
Converts a Point to WKT format: "POINT (x y)".
static std::string to_geojson(const Point3D &p)
Point3D → GeoJSON Point with Z.
Shared edge-group extraction for triangle meshes.
static void append_edge(Array< EdgeRef > &edges, size_t a, size_t b, const size_t tri, const size_t third)
static void for_each_sorted_edge_group(const TriangleArray &triangles, GroupFn on_group)
Gift wrapping (Jarvis march) convex hull algorithm.
static const Point * get_lowest_point(const DynList< Point > &point_set)
Finds the point with the minimum y-coordinate (and minimum x on ties).
Polygon operator()(const DynList< Point > &point_set) const
Compute the convex hull of a point set.
Graham scan convex hull algorithm.
static Array< Point > sorted_unique_points(const DynList< Point > &point_set)
Extracts points from a list, sorts them and removes duplicates.
static Geom_Number turn(const Point &a, const Point &b, const Point &c)
Computes the turn orientation of triple (a, b, c).
Polygon operator()(const DynList< Point > &point_set) const
Compute the convex hull of a point set.
void next_ne() noexcept
Move the iterator one position forward guaranteeing no exception.
bool has_curr() const noexcept
constexpr bool is_empty() const noexcept
size_t size() const noexcept
Count the number of elements of the list.
Exact bounded intersection of half-planes.
static Point line_intersection(const HalfPlane &a, const HalfPlane &b)
static Geom_Number cross_dir(const HalfPlane &a, const HalfPlane &b)
static Geom_Number dot_dir(const HalfPlane &a, const HalfPlane &b)
static Polygon build_polygon(const Array< Point > &pts)
static bool parallel(const HalfPlane &a, const HalfPlane &b)
static bool upper_half(const HalfPlane &h)
static Geom_Number signed_double_area(const Array< Point > &verts)
static bool same_direction(const HalfPlane &a, const HalfPlane &b)
static void push_clean(Array< Point > &out, const Point &p)
static Array< HalfPlane > from_convex_polygon(const Polygon &poly)
Build half-planes from the edges of a closed convex polygon.
static Array< Point > normalize_vertices(const Array< Point > &pts)
Polygon operator()(const Array< HalfPlane > &halfplanes) const
Intersect half-planes and return bounded feasible polygon.
Spatial point index for O(log n) nearest-neighbor queries.
DebugSnapshot debug_snapshot() const
Build a deterministic geometric partition snapshot for visualization.
bool is_empty() const noexcept
Return true if the tree is empty.
static KDTreePointSearch build(const Array< Point > &points, const Geom_Number &xmin, const Geom_Number &ymin, const Geom_Number &xmax, const Geom_Number &ymax)
Build a balanced KD-tree from a point array.
size_t size() const noexcept
Return the number of points in the tree.
bool contains(const Point &p) const
Check if a point exists in the tree.
static bool lexicographic_less(const Point &a, const Point &b)
KDTreePointSearch(const Geom_Number &xmin, const Geom_Number &ymin, const Geom_Number &xmax, const Geom_Number &ymax)
Construct an empty KD-tree for the given bounding region.
std::optional< Point > nearest(const Point &p) const
Find the nearest neighbor to query point p.
void range(const Geom_Number &xmin, const Geom_Number &ymin, const Geom_Number &xmax, const Geom_Number &ymax, DynList< Point > *out) const
Collect all points inside the given rectangle.
static Array< Point > unique_points(Array< Point > points)
bool insert(const Point &p)
Insert a point. Returns true if inserted, false if duplicate.
void for_each(Op &&op) const
Apply an operation to every point (inorder traversal).
Reusable event-driven sweep line framework.
DynSetTree< SeqEvent, Avl_Tree, CmpSeqEvent > queue_
void enqueue(Event &&e)
Enqueue an event (move version).
bool has_events() const noexcept
True if the event queue is non-empty.
Event dequeue()
Remove and return the next (minimum) event.
void enqueue(const Event &e)
Enqueue an event.
void run(Handler &&handler)
Run the sweep: process every event through handler.
size_t pending() const noexcept
Number of pending events.
void run(Handler &&handler, Array< Event > &out)
Run the sweep, collecting every event into out.
const Event & peek() const
Peek at the next event without removing it.
void clear() noexcept
Discard all pending events.
Smallest circle enclosing a point set (Welzl's algorithm).
Circle operator()(const DynList< Point > &points) const
Compute the minimum enclosing circle of a point set.
static Circle from_two_points(const Point &a, const Point &b)
Smallest circle with a and b on its boundary (diameter).
static Circle from_one_point(const Point &a)
Circle through a single point (degenerate, radius = 0).
static Circle mec_with_two_points(const Array< Point > &pts, size_t n, const Point &p1, const Point &p2)
Minimum enclosing circle with boundary points p1 and p2.
static Circle welzl_iterative(const Array< Point > &pts)
Iterative Welzl: three nested loops.
static Circle from_three_points(const Point &a, const Point &b, const Point &c)
Circumscribed circle through three points.
static Circle mec_with_point(const Array< Point > &pts, size_t n, const Point &p1)
Minimum enclosing circle with boundary point p1.
Circle operator()(std::initializer_list< Point > points) const
Convenience overload accepting an initializer list.
Exact Minkowski sum of two closed convex polygons.
static Array< Point > normalize(Array< Point > v)
Ensure CCW and rotate so that the bottom-most vertex is first.
static Array< Point > extract_vertices(const Polygon &poly)
static Geom_Number signed_double_area(const Array< Point > &v)
Polygon operator()(const Polygon &P, const Polygon &Q) const
Compute the Minkowski sum of two convex polygons.
static Point edge_vec(const Array< Point > &v, const size_t i)
static bool is_convex(const Array< Point > &verts)
static Geom_Number cross(const Point &a, const Point &b)
size_t left_edge_of_point(const Point &p, const Geom_Number &sweep_y) const
static void split_by_label(Node *root, const Geom_Number &label, Node *&left, Node *&right)
size_t successor_for_insert(const size_t edge, const Geom_Number &sweep_y) const
const Array< Point > & verts_
bool contains(const size_t edge) const
static Node * merge(Node *a, Node *b)
static Node * insert_by_label(Node *root, Node *node)
void erase(const size_t edge)
Geom_Number fresh_label_between(const size_t pred, const size_t succ) const
size_t predecessor_for_insert(const size_t edge, const Geom_Number &sweep_y) const
void insert(const size_t edge, const Geom_Number &sweep_y)
static Node * erase_by_label(Node *root, const Geom_Number &label, Node *&removed)
EdgeStatusTree(const Array< Point > &verts)
static void destroy(const Node *n)
O(n log n) triangulation of simple polygons via y-monotone partition + linear-time monotone triangula...
static Geom_Number signed_double_area(const Array< Point > &v)
static void ensure_ccw(Array< Point > &v)
static bool is_below(const Point &a, const Point &b)
Lexicographical "below" comparison.
static bool is_above(const Point &a, const Point &b)
Lexicographical "above" comparison (y primary, x secondary).
static bool edge_status_less(const size_t lhs, const size_t rhs, const Geom_Number &sweep_y, const Array< Point > &verts)
static Array< Array< size_t > > build_faces_from_diagonals(const Array< Point > &verts, const DynSetTree< std::pair< size_t, size_t > > &diagonals)
VertexType
Classification of vertices for y-monotone decomposition.
static bool is_y_monotone(const Array< Point > &v)
Check if a CCW polygon vertex array is y-monotone.
static Array< Array< size_t > > decompose_to_monotone_faces(const Array< Point > &verts)
DynList< Triangle > operator()(Polygon p) const
Triangulate a simple polygon.
static DynList< Triangle > triangulate_monotone(const Array< Point > &verts)
Triangulate a y-monotone polygon given as a vertex array.
static bool edge_goes_down(const Array< Point > &verts, const size_t edge)
Check if the edge starting at edge points downwards.
static VertexType classify_vertex(const Array< Point > &v, const size_t i)
static Array< Point > extract_vertices(const Polygon &p)
static Geom_Number edge_x_at_y(const Array< Point > &verts, const size_t edge, const Geom_Number &y)
static bool regular_interior_right(const Array< Point > &v, const size_t i)
Represents a point in 3D space with exact rational coordinates.
const Geom_Number & get_z() const
Gets the z-coordinate.
const Geom_Number & get_x() const
Gets the x-coordinate.
const Geom_Number & get_y() const
Gets the y-coordinate.
Exact point-in-polygon classification via winding number.
static bool strictly_contains(const Polygon &poly, const Point &p)
Return true only for strict interior points.
static Location locate(const Polygon &poly, const Point &p)
Classify point location with respect to a polygon.
static bool contains(const Polygon &poly, const Point &p)
Return true if the point is inside or on the boundary.
Represents a point with rectangular coordinates in a 2D plane.
const Geom_Number & get_x() const noexcept
Gets the x-coordinate value.
const Geom_Number & get_y() const noexcept
Gets the y-coordinate value.
Point midpoint(const Point &other) const
Calculates the midpoint between this point and another.
bool is_to_left_on_from(const Point &p1, const Point &p2) const
Checks if this point is to the left of or on the line from p1 to p2.
bool is_left_of(const Point &p1, const Point &p2) const
Checks if this point is to the left of the directed line from p1 to p2.
Geom_Number distance_squared_to(const Point &that) const
Calculates the squared Euclidean distance to another point.
Offset (inflate/deflate) an arbitrary simple polygon.
static Result cleanup(const Array< Point > &raw)
Full cleanup pipeline.
Result operator()(const Polygon &poly, const Geom_Number &distance, const JoinType join=JoinType::Miter, const Geom_Number &miter_limit=Geom_Number(2)) const
Offset a simple polygon by the given distance.
static Array< Polygon > extract_contours(Array< AugVertex > &aug)
Walk the augmented polygon and extract valid CCW contours.
JoinType
Strategy for connecting offset edges at vertices.
@ Bevel
Connect edge endpoints with a straight line.
@ Miter
Extend edges until they intersect.
static void sort_by_alpha(Array< size_t > &idx, const Array< Geom_Number > &keys)
Insertion-sort indices by alpha key.
Polygon offset_polygon(const Polygon &poly, const Geom_Number &distance, JoinType join=JoinType::Miter, const Geom_Number &miter_limit=Geom_Number(2)) const
Convenience: return the single largest-area result polygon.
static Polygon build_poly(const Array< Point > &pts)
Build a polygon from a point array.
static void offset_edge(const Point &a, const Point &b, const Geom_Number &d, const bool inward, Point &oa, Point &ob)
Compute two points on the offset line of edge (a→b) shifted by |d| along its outward or inward normal...
static Array< Point > compute_raw_offset(const Array< Point > &v, const Geom_Number &distance, const JoinType join, const Geom_Number &miter_limit)
static Array< Intersection > find_self_intersections(const Array< Point > &raw)
O(n²) scan for all proper self-intersections.
static Geom_Number abs_area(const Polygon &p)
static Array< AugVertex > build_augmented(const Array< Point > &raw, const Array< Intersection > &isects)
Build the augmented vertex list with crossing links.
static Point line_intersect(const Point &a1, const Point &a2, const Point &b1, const Point &b2)
Line–line intersection of (a1,a2) and (b1,b2).
static Geom_Number sq_dist(const Point &a, const Point &b)
Compute the squared distance between two points.
static Geom_Number edge_param(const Point &P, const Point &Q, const Point &I)
Compute parameter t of point I along directed edge P→Q.
Iterator over the edges (segments) of a polygon.
bool has_curr() const
Check if there is a current segment.
A general (irregular) 2D polygon defined by a sequence of vertices.
const Vertex & get_next_vertex(const Vertex &v) const
Get the vertex following v in the polygon.
const Vertex & get_prev_vertex(const Vertex &v) const
Get the vertex preceding v in the polygon.
void remove_vertex(const Vertex &v)
Remove a vertex from the polygon.
PointLocation locate_point(const Point &p) const
Classify a point against this closed polygon.
const Vertex & get_first_vertex() const
Get the first vertex of the polygon.
void add_vertex(const Point &point)
Add a vertex to the polygon.
void close()
Close the polygon.
const bool & is_closed() const
Check if the polygon is closed.
const size_t & size() const
Get the number of vertices.
Power diagram (weighted Voronoi diagram).
static Point power_center(const WeightedSite &a, const WeightedSite &b, const WeightedSite &c)
Compute the power center of three weighted sites.
Result operator()(const Array< WeightedSite > &sites) const
Compute the power diagram.
QuickHull convex hull algorithm.
static std::pair< DynList< Point >, DynList< Point > > partition(const DynList< Point > &point_set, const Point &a, const Point &b)
Split points by the line through directed segment (a,b).
Polygon operator()(const DynList< Point > &point_set) const
Compute the convex hull of a point set.
static std::pair< Point, Point > search_extremes(const DynList< Point > &point_set)
Find leftmost and rightmost points by x coordinate.
static DynList< Point > quick_hull(DynList< Point > &point_set, const Point &a, const Point &b)
Recursive QuickHull step for a directed edge (a,b).
static std::pair< DynList< Point >, DynList< Point > > get_right_points(DynList< Point > &point_set, const Point &a, const Point &b, const Point &c)
Partition points to the right of edges (a,c) and (c,b).
static Point get_farthest_point(const DynList< Point > &point_set, const Segment &s)
Return the point in point_set farthest from segment s.
Static 2D range tree for orthogonal range queries.
void build_node(const size_t node, const size_t lo, const size_t hi)
bool is_empty() const noexcept
DynList< Point > query(const Geom_Number &xmin, const Geom_Number &xmax, const Geom_Number &ymin, const Geom_Number &ymax) const
Query: return all points inside [xmin,xmax] × [ymin,ymax].
Array< Point > pts_
points sorted by x (primary key)
Array< Node > tree_
implicit binary tree (1-indexed)
size_t size() const noexcept
size_t lower_bound_x(const Geom_Number &xval) const
Binary search: first index i in pts_ where pts_(i).get_x() >= xval.
static size_t lower_bound_y(const Array< Point > &arr, const Geom_Number &ymin)
Binary search: first index i in arr where arr(i).get_y() >= ymin.
void query_range(const size_t node, const size_t lo, const size_t hi, const size_t qlo, const size_t qhi, const Geom_Number &ymin, const Geom_Number &ymax, DynList< Point > &out) const
size_t upper_bound_x(const Geom_Number &xval) const
Binary search: first index i in pts_ where pts_(i).get_x() > xval.
void build(const DynList< Point > &points)
Build the range tree from a point set.
static size_t upper_bound_y(const Array< Point > &arr, const Geom_Number &ymax)
Binary search: first index i in arr where arr(i).get_y() > ymax.
DebugSnapshot debug_snapshot() const
Return structural snapshot for visualization/debug.
An axis-aligned rectangle.
const Geom_Number & get_xmin() const
Gets the minimum x-coordinate.
const Geom_Number & get_ymax() const
Gets the maximum y-coordinate.
const Geom_Number & get_ymin() const
Gets the minimum y-coordinate.
const Geom_Number & get_xmax() const
Gets the maximum x-coordinate.
Exact regular triangulation (weighted Delaunay) via Bowyer-Watson.
static bool all_collinear(const Array< WeightedSite > &s)
static bool point_in_power_circle(const Array< Point > &pts, const Array< Geom_Number > &wts, const size_t ia, const size_t ib, const size_t ic, const size_t ip)
Weighted in-circle (power) predicate.
static bool lex_less(const Point &p1, const Point &p2)
static Array< WeightedSite > unique_sites(const Array< WeightedSite > &input)
Result operator()(const Array< WeightedSite > &weighted_sites) const
Compute the regular triangulation of a weighted point set.
Rotating calipers metrics for convex polygons.
static DiameterResult make_diameter(const Point &a, const Point &b)
static Array< Point > extract_vertices(const Polygon &poly)
Extract polygon vertices into an array.
static DiameterResult diameter(const Polygon &poly)
Compute convex polygon diameter (farthest vertex pair).
static WidthResult minimum_width(const Polygon &poly)
Compute the minimum width of a closed convex polygon.
static bool is_convex(const Array< Point > &verts)
Check if a polygon vertex cycle is convex.
Compute the full planar subdivision induced by a set of segments.
static size_t find_vertex(const Array< Point > &verts, const Point &p)
Find the index of point p in a sorted vertex array (binary search).
static bool angle_lt(const Geom_Number &dx1, const Geom_Number &dy1, const Geom_Number &dx2, const Geom_Number &dy2)
Compare two directions from the same origin by angle (exact).
static constexpr size_t NONE
static int angle_quad(const Geom_Number &dx, const Geom_Number &dy)
Angular quadrant of direction (dx, dy).
Result operator()(const Array< Segment > &segments) const
Compute the arrangement of a set of segments.
static bool pt_less(const Point &a, const Point &b)
Lexicographic point comparison.
Dedicated exact intersection for a single pair of segments.
Kind
Types of intersection between two segments.
@ OVERLAP
Intersection over a collinear interval.
@ POINT
Intersection at a single point.
static Segment collinear_overlap(const Segment &s1, const Segment &s2)
static Point point_max_on_axis(const Point &a, const Point &b, const bool vertical_axis)
static bool point_less_on_axis(const Point &a, const Point &b, const bool vertical_axis)
Result operator()(const Segment &s1, const Segment &s2) const
Computes the intersection of two line segments.
static Point point_min_on_axis(const Point &a, const Point &b, const bool vertical_axis)
Represents a line segment between two points.
bool intersects_properly_with(const Segment &s) const
Checks if this segment properly intersects another segment.
Point intersection_with(const Segment &s) const
Computes the intersection point of the infinite lines defined by two segments.
bool intersects_with(const Segment &s) const
Checks if this segment intersects another one (including endpoints and collinear overlap).
Point project(const Point &p) const
Orthogonal projection of a point onto this segment's infinite line, clamped to the segment's endpoint...
Geom_Number distance_to(const Point &p) const
Calculates the Euclidean distance from a point to this segment.
bool contains(const Point &p) const
Checks if a point lies on this segment.
const Point & get_tgt_point() const noexcept
Gets the target point of the segment.
const Point & get_src_point() const noexcept
Gets the source point of the segment.
Compute the shortest Euclidean path between two points inside a simple polygon.
static size_t find_tri(const Array< Point > &pts, const Array< ITri > &tris, const Point &p)
Find a triangle containing point p.
static constexpr size_t NONE
static Array< size_t > find_sleeve(const Array< ITri > &tris, const size_t src, const size_t dst)
BFS on dual graph → sleeve.
static Array< ITri > build_tris(const Array< Point > &pts, const DynList< Triangle > &tl)
Build indexed triangulation with adjacency.
static Geom_Number cross(const Point &a, const Point &b, const Point &c)
Cross product (b-a) x (c-a).
DynList< Point > operator()(const Polygon &polygon, const Point &source, const Point &target) const
Compute the shortest path between two points in a simple polygon.
static bool point_in_triangle(const Array< Point > &pts, const ITri &t, const Point &p)
Point-in-triangle test.
static size_t find_index(const Array< Point > &pts, const Point &p)
Match a Point from a Triangle to an index in pts (exact comparison).
Sweep-line status as a balanced tree with O(log n) updates.
size_t predecessor_for_insert(const size_t seg, const Geom_Number &sx) const
StatusTree(const Array< Segment > &segs)
bool contains(const size_t seg) const
static Node * insert_by_label(Node *root, Node *node)
Geom_Number fresh_label_between(const size_t pred, const size_t succ) const
static void destroy(const Node *n)
size_t predecessor(const size_t seg) const
static Node * merge(Node *a, Node *b)
size_t successor(const size_t seg) const
void insert(const size_t seg, const Geom_Number &sx)
size_t successor_for_insert(const size_t seg, const Geom_Number &sx) const
static void split_by_label(Node *root, const Geom_Number &label, Node *&left, Node *&right)
void erase(const size_t seg)
static Node * erase_by_label(Node *root, const Geom_Number &label, Node *&removed)
const Array< Segment > & segs_
Report all pairwise intersection points among a set of segments.
static bool status_less(const size_t lhs, const size_t rhs, const Geom_Number &sx, const Array< Segment > &segs)
static void check_and_enqueue(const Array< Segment > &segs, const size_t i, const size_t j, const Geom_Number &sx, LineSweepFramework< Event, EventLess > &eq, DynSetTree< size_t, Treap_Rk > &seen_pairs, const size_t n)
Detect an intersection and enqueue it as a future event.
static bool slope_less(const Segment &a, const Segment &b)
static Segment canonicalize(const Segment &s)
Canonicalize a segment so its src is the "left" endpoint.
Array< Intersection > operator()(const Array< Segment > &segments) const
Find all pairwise segment intersection points.
static bool event_less(const Event &a, const Event &b)
Ordering for events in the sweep queue.
static Geom_Number y_at_x(const Segment &s, const Geom_Number &x)
Evaluate the y-coordinate of segment s at x = x.
EventType
Types of events in the Bentley-Ottmann sweep.
O(log n) point location via trapezoidal map with DAG search.
static void split_multiple_trapezoids(Array< Point > &pts, Array< Segment > &segs, Array< Trapezoid > &traps, Array< DagNode > &dag, const size_t seg_idx, const Array< size_t > &crossed)
Handle the case where a segment crosses multiple trapezoids.
LocationType
Classification of a point location query result.
@ TRAPEZOID
The point is strictly inside a trapezoid.
@ ON_POINT
The point coincides with a segment endpoint.
@ ON_SEGMENT
The point lies on a segment boundary.
Result operator()(const Array< Polygon > &polygons) const
Build a trapezoidal map from multiple closed polygons.
static constexpr size_t NONE
static void update_neighbor(Array< Trapezoid > &traps, const size_t neighbor_idx, const size_t old_ti, const size_t new_ti)
Update all neighbor pointers that referenced old_ti to point to new_ti.
static size_t dag_locate(const Array< Point > &pts, const Array< Segment > &segs, const Array< DagNode > &dag, const Point &p, size_t root)
Locate the DAG leaf (trapezoid) containing point p.
static Array< size_t > find_crossed_trapezoids(const Array< Point > &pts, const Array< Segment > &segs, const Array< Trapezoid > &traps, const size_t seg_idx, const size_t start_trap)
Find all trapezoids crossed by segment seg_idx, starting from the trapezoid containing the left endpo...
static bool point_above_segment(const Point &p, const Segment &seg)
Test if point p is above segment seg (left of the directed segment).
static size_t make_trapezoid(Array< Trapezoid > &traps, Array< DagNode > &dag, const size_t top, const size_t bottom, const size_t leftp, const size_t rightp)
Create a new trapezoid and a DAG leaf for it.
Result operator()(const Array< Segment > &input_segments) const
Build a trapezoidal map from a set of non-crossing segments.
static void split_single_trapezoid(Array< Point > &pts, Array< Segment > &segs, Array< Trapezoid > &traps, Array< DagNode > &dag, const size_t seg_idx, const size_t trap_idx)
Handle the case where a segment is fully contained in one trapezoid.
static void replace_leaf(Array< DagNode > &dag, const size_t leaf_idx, const NodeType type, const size_t index, const size_t left, const size_t right)
Replace a DAG leaf with an internal node (X or Y).
Result operator()(const Polygon &polygon) const
Build a trapezoidal map from a closed polygon's edges.
NodeType
Types of nodes in the search DAG.
@ X_NODE
A node that splits the search space by an X-coordinate.
@ LEAF
A leaf node representing a trapezoid.
@ Y_NODE
A node that splits the search space by a segment (Y-direction).
A non-degenerate triangle defined by three points.
const Point & get_p3() const
Gets the third vertex.
const Point & get_p2() const
Gets the second vertex.
const Point & get_p1() const
Gets the first vertex.
A vertex in a polygon's doubly linked vertex list.
const Vertex & next_vertex() const
Get the next vertex in the polygon.
Point to_point() const
Return this vertex as a plain Point value.
BST-based sweep status for active edges, ordered by ray-intersection distance.
void erase(const size_t edge)
size_t min() const
Return the nearest edge (smallest ray_param).
DynSetTree< EdgeKey, Treap, EdgeKeyCmp > tree_
EdgeStatusTree(const Array< Point > &verts, const size_t n, const Point &query)
bool contains(const size_t edge) const
const Array< Point > & verts_
Array< Geom_Number > params_
void insert(const size_t edge, const Point &dir)
Compute the visibility polygon from a point inside a simple polygon.
static bool angle_less(const Point &q, const Point &a, const Point &b)
True if direction (a - q) has a smaller angle than (b - q).
static Point ray_edge_hit(const Point &q, const Point &dir, const Point &e0, const Point &e1)
Intersection point of ray from q through dir with edge (e0, e1).
Polygon operator()(const Polygon &polygon, const Point &query) const
Compute the visibility polygon.
static int angle_quadrant(const Geom_Number &dx, const Geom_Number &dy)
Quadrant of direction (dx, dy): 0 = +x+y, 1 = -x+y, 2 = -x-y, 3 = +x-y.
static Geom_Number ray_param(const Point &q, const Point &dir, const Point &e0, const Point &e1)
Parametric t along ray q + t*(dir-q) for intersection with edge (e0,e1).
Visvalingam-Whyatt polyline simplification.
static Array< Point > simplify_open_array(const Array< Point > &polyline, const Geom_Number &area_threshold)
Array< Point > operator()(const Array< Point > &polyline, const Geom_Number &area_threshold) const
Simplify an open polyline.
static Geom_Number effective_area(const Point &a, const Point &b, const Point &c)
Array< Point > operator()(const DynList< Point > &polyline, const Geom_Number &area_threshold) const
Overload accepting DynList.
static Polygon simplify_polygon(const Polygon &poly, const Geom_Number &area_threshold)
Simplify a closed polygon.
Fortune sweep-line Voronoi construction.
Array< ClippedCell > clipped_cells(const DynList< Point > &pts, const Polygon &clip) const
Compute Voronoi cells clipped to a bounding polygon.
static bool all_collinear(const Array< Point > &pts)
static void enqueue_circle_event(Arc *arc, const double sweepline_y, const Array< Point > &sites, DynBinHeap< Event *, EventCmp > &queue, Array< std::unique_ptr< Event > > &event_pool, size_t &next_event_id)
VoronoiDiagramFromDelaunay voronoi_
Result operator()(const std::initializer_list< Point > il) const
Compute the Voronoi diagram from an initializer list.
DelaunayTriangulationBowyerWatson fallback_
static Arc * locate_arc(Arc *head, const Array< Point > &sites, const double x, const double sweepline_y)
static double as_double(const Geom_Number &v)
static void invalidate_circle_event(Arc *arc)
Array< ClippedCell > clipped_cells(const std::initializer_list< Point > il, const Polygon &clip) const
This is an overloaded member function, provided for convenience. It differs from the above function o...
static DelaunayTriangulationBowyerWatson::IndexedTriangle normalized_triangle(const size_t i, const size_t j, const size_t k, const Array< Point > &sites)
static constexpr double kEps
static double breakpoint_x(const Point &left_site, const Point &right_site, const double sweepline_y)
static bool is_valid_delaunay(const DelaunayTriangulationBowyerWatson::Result &dt)
static DelaunayTriangulationBowyerWatson::Result triangulate_sweep(const Array< Point > &sites)
Result operator()(const DynList< Point > &pts) const
Compute the Voronoi diagram for a list of points.
Voronoi diagram derived as the dual of a Delaunay triangulation.
static Array< ClippedCell > indexed_clipped_cells(const Array< Point > &sites, const Array< Polygon > &polys)
Array< ClippedCell > clipped_cells_indexed(const std::initializer_list< Point > il, const Polygon &clip) const
Convenience overload with initializer-list input.
static Point circumcenter(const Point &a, const Point &b, const Point &c)
Computes the circumcenter of three points.
Result operator()(const DelaunayTriangulationBowyerWatson::Result &dt) const
Build Voronoi from a precomputed Delaunay triangulation.
static HalfPlaneIntersection::HalfPlane bisector_halfplane_for_site(const Point &s, const Point &t)
Creates a half-plane from a directed edge.
static bool is_convex(const Array< Point > &verts)
Checks if a vertex set forms a convex polygon.
static Array< ClippedCell > clipped_cells_indexed(const Array< Point > &sites, const Polygon &clip)
Clip Voronoi cells and return explicit site-indexed records.
Array< Polygon > clipped_cells(const std::initializer_list< Point > il, const Polygon &clip) const
Convenience overload with initializer-list input.
Array< ClippedCell > clipped_cells_indexed(const DynList< Point > &point_set, const Polygon &clip) const
Compute/clip Voronoi cells and return site-indexed records.
static Array< Polygon > clipped_cells(const Array< Point > &sites, const Polygon &clip)
Clip Voronoi cells against a closed convex polygon.
static Array< ClippedCell > clipped_cells_indexed(const Result &vor, const Polygon &clip)
Clip Voronoi cells (from result) into site-indexed records.
static Array< Point > extract_vertices(const Polygon &p)
Extracts vertices from a closed polygon.
Array< Polygon > clipped_cells(const DynList< Point > &point_set, const Polygon &clip) const
Compute Voronoi and clip its cells against a closed convex polygon.
DelaunayTriangulationBowyerWatson delaunay
Internal Delaunay triangulator.
static Array< Polygon > clipped_cells(const Result &vor, const Polygon &clip)
Clip Voronoi cells against a closed convex polygon.
O(n log n) Voronoi diagram construction.
VoronoiDiagramFromDelaunay voronoi_
Internal dual builder.
Array< ClippedCell > clipped_cells(const DynList< Point > &pts, const Polygon &clip) const
Compute Voronoi cells clipped to a bounding polygon.
DelaunayTriangulationRandomizedIncremental delaunay_
Internal Delaunay triangulator.
Result operator()(const std::initializer_list< Point > il) const
Compute the Voronoi diagram from an initializer list.
Result operator()(const DynList< Point > &pts) const
Compute the Voronoi diagram for a list of points.
void for_each(Operation &operation)
Traverse all the container and performs an operation on each element.
2D k-d tree for efficient spatial point operations.
bool insert(const Point &p)
Insert a point into the tree.
constexpr size_t size() const noexcept
Get the number of points in the tree.
static void range(Node *root, const Rectangle &rect, DynList< Point > *q)
Recursively find points within a rectangle.
void for_each(Op &&op) const
Apply an operation to every point in the tree (inorder).
static K2Tree build(Array< Point > points, const Point &pmin, const Point &pmax)
Build a balanced k-d tree from an array of points.
bool contains(const Point &p) const noexcept
Check if the tree contains a point.
std::optional< Point > nearest(const Point &p) const noexcept
Find the nearest point to a query point.
constexpr bool is_empty() const noexcept
Check if the tree is empty.
constexpr size_t size() const noexcept
Returns the number of entries in the table.
__gmp_expr< typename __gmp_resolve_expr< T, V >::value_type, __gmp_binary_expr< __gmp_expr< T, U >, __gmp_expr< V, W >, __gmp_hypot_function > > hypot(const __gmp_expr< T, U > &expr1, const __gmp_expr< V, W > &expr2)
__gmp_expr< T, __gmp_unary_expr< __gmp_expr< T, U >, __gmp_y1_function > > y1(const __gmp_expr< T, U > &expr)
__gmp_expr< mpfr_t, mpfr_t > mpfr_class
__gmp_expr< T, __gmp_binary_expr< __gmp_expr< T, U >, unsigned long int, __gmp_root_function > > root(const __gmp_expr< T, U > &expr, unsigned long int l)
bool vertical
If true, use vertical layout (default).
Singly linked list implementations with head-tail access.
Freq_Node * pred
Predecessor node in level-order traversal.
Main namespace for Aleph-w library functions.
and
Check uniqueness with explicit hash + equality functors.
void in_place_unique(Container &c, Compare cmp={})
Remove duplicates in-place preserving first occurrence order.
Orientation
Classification of three-point orientation.
Geom_Number in_circle_determinant(const Point &a, const Point &b, const Point &c, const Point &p)
Return the exact in-circle determinant for (a,b,c,p).
bool all(Container &container, Operation &operation)
Return true if all elements satisfy a predicate.
bool eq(const C1 &c1, const C2 &c2, Eq e=Eq())
Check equality of two containers using a predicate.
Geom_Number area_of_parallelogram(const Point &a, const Point &b, const Point &c)
Compute the signed area of the parallelogram defined by vectors a->b and a->c.
Itor unique(Itor __first, Itor __last, BinaryPredicate __binary_pred=BinaryPredicate())
Remove consecutive duplicates in place.
bool on_segment(const Segment &s, const Point &p)
Return true if p lies on segment s (exact).
size_t size(Node *root) noexcept
bool segments_intersect(const Segment &s1, const Segment &s2)
Return true if segments s1 and s2 intersect (including endpoints).
Divide_Conquer_DP_Result< Cost > divide_and_conquer_partition_dp(const size_t groups, const size_t n, Transition_Cost_Fn transition_cost, const Cost inf=dp_optimization_detail::default_inf< Cost >())
Optimize partition DP using divide-and-conquer optimization.
std::decay_t< typename HeadC::Item_Type > T
Point segment_intersection_point(const Segment &s1, const Segment &s2)
Compute the exact intersection point of segments s1 and s2.
bool contains(const std::string_view &str, const std::string_view &substr)
Check if substr appears inside str.
DynArray< T > & in_place_sort(DynArray< T > &c, Cmp cmp=Cmp())
Sorts a DynArray in place.
Geom_Number square_root(const Geom_Number &x)
Square root of x (wrapper over mpfr).
double geom_number_to_double(const Geom_Number &n)
Converts a Geom_Number to its double precision representation.
Orientation orientation(const Point &a, const Point &b, const Point &c)
Return the orientation of the triple (a, b, c).
mpq_class Geom_Number
Numeric type used by the geometry module.
SegmentSegmentIntersection::Result segment_segment_intersection(const Segment &s1, const Segment &s2)
Convenience free-function wrapper for SegmentSegmentIntersection.
std::ostream & join(const C &c, const std::string &sep, std::ostream &out)
Join elements of an Aleph-style container into a stream.
void next()
Advance all underlying iterators (bounds-checked).
T sum(const Container &container, const T &init=T{})
Compute sum of all elements.
void quicksort_op(C< T > &a, const Compare &cmp=Compare(), const size_t threshold=Quicksort_Threshold)
Optimized quicksort for containers using operator().
2D polygon representation and geometric operations.
Structure used for debugging and visualizing the tree structure.
size_t user_index
Entry::index.
size_t node_index
Index of the node.
size_t entry_index
Index in entries_ array.
size_t right
Right child index.
size_t left
Left child index.
Rectangle bbox
Bounding box of the node.
size_t depth
Depth in the tree.
bool is_leaf
True if it's a leaf node.
A snapshot of the entire tree for debugging.
size_t root
Index of the root node.
Array< DebugNode > nodes
Array of all nodes in debug format.
An entry in the tree, consisting of a bounding box and a user-defined index.
size_t index
User-defined identifier for the entry.
Rectangle bbox
The axis-aligned bounding box.
Internal tree node structure.
bool is_leaf() const
Checks if this node is a leaf.
Rectangle bbox
Bounding box enclosing all descendants.
size_t entry_idx
Index in entries_ (leaf only).
size_t left
Left child index.
size_t right
Right child index.
Result of an alpha-shape computation.
Array< DelaunayTriangulationBowyerWatson::IndexedTriangle > triangles
Triangles that satisfy the alpha criterion.
Array< Point > sites
Deduplicated and sorted input sites.
Array< Segment > boundary_edges
Segments forming the boundary of the alpha-shape.
Lexicographical comparison for points (x primary, y secondary).
bool operator()(const Point &p1, const Point &p2) const
De Casteljau subdivision: split a cubic Bézier at parameter t into two cubic Béziers (left and right)...
Intersection pair found between edge i of A and edge j of B.
Comparator for segments to ensure strict weak ordering in sets.
bool operator()(const Segment &s1, const Segment &s2) const
Compares two segments lexicographically.
static bool cmp_point(const Point &p1, const Point &p2)
Strict lexicographical comparison for points.
Strict weak order for sorting points by (y,x).
bool operator()(const Point &p1, const Point &p2) const
Strict weak order for sorting points lexicographically by (x,y).
bool operator()(const Point &p1, const Point &p2) const
Geom_Number distance_squared
Collect edges that cross segment (u,v) by scanning all mesh edges.
size_t local
local index in tri (opposite vertex)
Array< IndexedEdge > constrained_edges
Array< IndexedTriangle > triangles
size_t adj[3]
adj[i] = neighbor across edge opposite v[i]
size_t v[3]
CCW vertex indices.
bool constrained[3]
constrained[i] = edge opp v[i] is constrained
std::vector< SupportPoint > reduced_simplex
Holds the results of the distance computation.
Geom_Number distance
The minimum Euclidean distance.
Geom_Number distance_squared
The minimum squared distance.
Point closest_on_second
Witness point on the second polygon.
bool intersects
True if the polygons overlap.
Point closest_on_first
Witness point on the first polygon.
size_t gjk_iterations
Count of iterations performed.
A point in the Minkowski difference and its source components.
Point a
Contributing point from the first polygon.
Point b
Contributing point from the second polygon.
Point v
Resulting point in Minkowski difference (a - b).
double vy
Double-precision cache for faster predicates.
double vx
Double-precision cache for faster predicates.
Represents a triangle by the indices of its three vertices.
size_t j
Index of the second vertex.
size_t i
Index of the first vertex.
size_t k
Index of the third vertex.
The result of a Delaunay triangulation.
Array< IndexedTriangle > triangles
Triangles forming the Delaunay triangulation.
Array< Point > sites
Unique, sorted input points used for triangulation.
size_t v[3]
vertex indices, CCW order
size_t adj[3]
adj[i] = neighbor across edge opposite v[i]
bool operator()(const UndirectedEdge &a, const UndirectedEdge &b) const
Represents a reference to a triangle edge.
size_t v
Index of the second vertex of the edge (u < v).
size_t third
Index of the third vertex in the triangle (not on the edge).
size_t u
Index of the first vertex of the edge (u < v).
size_t tri
Index of the triangle containing this edge.
Lexicographical comparison for points.
bool operator()(const Point &p1, const Point &p2) const
bool outside(const Point &x) const
Geom_Number offset() const
HalfPlane(Point p_, Point q_)
Represents a spatial partition or leaf in the KD-tree.
bool split_on_x
True if split is vertical, false if horizontal.
Rectangle region
Geometric region covered by this node.
Geom_Number split_value
The coordinate value of the split.
Point representative
Point stored in the leaf.
bool is_leaf
True if it's a leaf node.
size_t depth
Depth of the node in the tree.
A complete snapshot of the tree for debugging.
Array< DebugPartition > partitions
List of all internal and leaf partitions.
Array< Point > points
All points stored in the tree.
Rectangle bounds
Global bounding box of the tree.
Comparator for SeqEvent, ensuring a strict weak ordering.
CmpEvent cmp
The user-provided base comparator.
bool operator()(const SeqEvent &a, const SeqEvent &b) const
Internal event wrapper that adds a sequence number for tie-breaking.
size_t seq
Insertion sequence number.
Event event
The user-defined event payload.
Result type: a circle defined by center and squared radius.
Geom_Number radius() const
Exact radius (square root of radius_squared).
bool contains(const Point &p) const
True if point p lies inside or on the circle boundary.
Geom_Number radius_squared
Augmented vertex in the cleaned polygon.
Intersection found between edge i and edge j of the raw polygon.
The result of an offset operation, which may produce multiple polygons.
Array< Polygon > polygons
Set of resulting simple polygons.
size_t size() const noexcept
bool is_empty() const noexcept
Checks if the result contains no polygons.
Iterator over the vertices of a polygon.
Represents a cell in the power diagram.
size_t site_index
Index of the site this cell belongs to.
bool bounded
True if the cell is completely bounded.
Array< Point > vertices
Vertices of the cell in CCW order.
Point site
Coordinates of the site.
Geom_Number weight
Weight of the site.
Represents an edge in the power diagram.
size_t site_v
Index of the second site sharing this edge.
Point tgt
Target point of the edge (if not unbounded).
size_t site_u
Index of the first site sharing this edge.
Point direction
Direction vector for unbounded rays.
bool unbounded
True if the edge is an unbounded ray.
Point src
Source point of the edge.
Complete result of a power diagram computation.
Array< PowerCell > cells
Cells of the power diagram.
Array< Point > vertices
Vertices of the power diagram.
Array< WeightedSite > sites
Sorted and unique input sites.
Array< PowerEdge > edges
Edges of the power diagram.
A site with an associated weight (squared radius).
Geom_Number weight
Weight of the site.
Point position
Coordinates of the site.
Represents a node in the range tree for debugging.
size_t left
Left child tree index.
size_t y_sorted_size
Number of points in the y-sorted secondary structure.
size_t right
Right child tree index.
Geom_Number split_x
X-coordinate used to split this node.
bool is_leaf
True if this is a leaf node.
Geom_Number xmin
Minimum x-coordinate in this subtree.
size_t tree_index
Implicit tree index (1-based).
size_t lo
Lower index into the x-sorted point array.
Geom_Number xmax
Maximum x-coordinate in this subtree.
size_t hi
Upper index into the x-sorted point array.
A complete snapshot of the range tree for debugging.
Array< DebugNode > nodes
List of all nodes in debug format.
Array< Point > x_sorted_points
The underlying x-sorted point array.
Internal node structure containing the y-sorted secondary array.
Array< Point > y_sorted
Points in this node's x-range, sorted by y.
Array< IndexedTriangle > triangles
Array< WeightedSite > sites
A 2-D site with an associated weight (squared radius).
Geom_Number distance_squared
Geom_Number width_squared
Represents an edge in the planar subdivision.
size_t tgt
Index of the target vertex.
size_t src
Index of the source vertex.
size_t seg_idx
Index of the original segment that generated this edge.
Represents a face (region) in the planar subdivision.
bool unbounded
True if this is the infinite outer face.
DynList< size_t > boundary
Ordered indices of vertices forming the face boundary.
Internal half-edge structure used for face traversal.
size_t edge_idx
Index of the corresponding undirected edge.
size_t twin
Index of the twin half-edge going in the opposite direction.
size_t next
Index of the next half-edge in CCW order around the face.
size_t origin
Index of the vertex where this half-edge starts.
size_t face
Index of the face this half-edge borders.
size_t target
Index of the vertex where this half-edge ends.
The complete result of the segment arrangement computation.
Array< Point > vertices
Unique vertices (endpoints and intersection points).
Array< ArrEdge > edges
Sub-segments forming the arrangement graph.
Array< ArrFace > faces
Faces defined by the arrangement.
Detailed result of a segment-segment intersection test.
bool intersects() const noexcept
Returns true if any intersection exists.
Point point
The intersection point (if kind is POINT).
Segment overlap
The overlap segment (if kind is OVERLAP).
Kind kind
Type of intersection.
Functor wrapper for event_less (used by LineSweepFramework).
bool operator()(const Event &a, const Event &b) const
Represents an event in the sweep-line algorithm.
Point pt
Coordinates of the event.
size_t seg_b
Index of the second segment (only for INTERSECTION).
EventType type
Type of the event.
size_t seg_a
Index of the first segment involved.
A single intersection record.
size_t seg_j
Index of second segment (seg_i < seg_j).
Point point
Intersection point.
size_t seg_i
Index of the first segment.
A node in the search Directed Acyclic Graph (DAG).
NodeType type
The type of the node.
size_t index
Index into points, segments, or trapezoids depending on type.
size_t right
Index of the right child (or above/right child).
size_t left
Index of the left child (or below/left child).
Result of a point location query.
size_t segment_index
Index of the segment the point lies on.
size_t trapezoid_index
Index of the containing trapezoid (if any).
size_t point_index
Index of the point the query point coincides with.
LocationType type
Classification of the point location.
The trapezoidal map and DAG search structure.
Array< Trapezoid > trapezoids
All trapezoids created during construction.
Array< Segment > segments
Input segments.
Array< DagNode > dag
The search directed acyclic graph.
LocationResult locate(const Point &p) const
Locate a point in the trapezoidal map.
size_t num_input_points
Number of unique points derived from segments.
Array< Point > points
Unique endpoints of all segments.
size_t num_input_segments
Number of segments originally provided.
size_t dag_root
Index of the DAG root node.
bool contains(const Point &p) const
Check if point is inside any input polygon.
A trapezoid in the decomposition.
size_t leftp
Index into points array (left vertical wall).
size_t dag_leaf
Back-pointer to the corresponding DAG leaf node.
size_t lower_left
Neighbor trapezoid to the lower left.
size_t rightp
Index into points array (right vertical wall).
size_t upper_right
Neighbor trapezoid to the upper right.
bool active
True if this trapezoid is part of the current map.
size_t lower_right
Neighbor trapezoid to the lower right.
size_t upper_left
Neighbor trapezoid to the upper left.
size_t bottom
Index into segments array (lower boundary).
size_t top
Index into segments array (upper boundary).
bool operator()(const EdgeKey &a, const EdgeKey &b) const
Entry in the priority queue for vertex removal.
size_t idx
Index of the vertex in the polyline.
bool operator<(const VWEntry &o) const
Comparison based on effective area.
Geom_Number area
Effective area of the triangle formed with its neighbors.
bool operator>(const VWEntry &o) const
Comparison based on effective area.
Arc in the beach-line state structure.
Arc * next
Next arc in the beach-line.
Event * circle
Potential circle event generated by this arc.
Arc * prev
Previous arc in the beach-line.
size_t site
Site index generating this parabolic arc.
Comparator for prioritizing events in the sweep-line.
bool operator()(const Event *a, const Event *b) const
Forward declaration of beach-line arc.
bool valid
False if the circle event was invalidated.
double y
Y-coordinate of the event.
bool is_site
True if it's a site event, false for circle.
Arc * arc
Corresponding arc in the beach-line (circle events).
size_t id
Unique event ID for deterministic sorting.
double x
X-coordinate of the site or circle center.
Key for identifying a triangle by its vertex indices.
bool operator<(const TriKey &o) const
Represents a Voronoi cell.
size_t site_index
Index of the site this cell belongs to.
Array< Point > vertices
Vertices of the cell in CCW order.
Point site
Coordinates of the site.
bool bounded
True if the cell is completely bounded.
Represents a Voronoi cell clipped by a bounding polygon.
size_t site_index
Index of the site.
Polygon polygon
Clipped convex polygon representing the cell.
Point site
Coordinates of the site.
Represents an edge in the Voronoi diagram.
Point tgt
Target point of the edge (only if not unbounded).
Point src
Source point of the edge.
size_t site_v
Index of the second site sharing this edge.
Point direction
Direction vector for unbounded rays.
size_t site_u
Index of the first site sharing this edge.
bool unbounded
True if the edge is an unbounded ray.
The complete result of a Voronoi diagram computation.
Array< Edge > edges
Voronoi edges.
Array< Point > vertices
Voronoi vertices.
Array< Cell > cells
Voronoi cells.
Array< Point > sites
Input sites.
FooMap m(5, fst_unit_pair_hash, snd_unit_pair_hash)
2D k-d tree implementation for spatial point indexing.
Circular queue implementations backed by arrays.
Stack implementations backed by dynamic or fixed arrays.
Dynamic binary heap with node-based storage.
Dynamic doubly linked list implementation.
Dynamic set implementations based on balanced binary search trees.
Comprehensive sorting algorithms and search utilities for Aleph-w.