From 4c976272c53e5d62afbfccd7c31cd6eaf4a52f54 Mon Sep 17 00:00:00 2001 From: lisyarus Date: Fri, 16 Sep 2022 18:30:39 +0300 Subject: [PATCH] Switch from define-based to tag-based robust predicates switch --- libs/cg/include/psemek/cg/convex/inside.hpp | 12 +++-- .../psemek/cg/convex_hull_2d/andrew.hpp | 14 +++-- .../psemek/cg/convex_hull_2d/graham.hpp | 14 +++-- .../psemek/cg/convex_hull_2d/jarvis.hpp | 12 +++-- .../psemek/cg/convex_hull_2d/quickhull.hpp | 28 ++++++---- .../psemek/cg/convex_hull_3d/quickhull.hpp | 24 +++++---- .../psemek/cg/triangulation/delaunay.hpp | 14 +++-- .../psemek/cg/triangulation/ear_clipping.hpp | 26 +++++++--- .../psemek/cg/triangulation/monotone.hpp | 52 ++++++++++++------- .../psemek/cg/triangulation/triangulation.hpp | 16 ++++-- libs/geom/CMakeLists.txt | 12 ++--- libs/geom/include/psemek/geom/contains.hpp | 16 ++++-- libs/geom/include/psemek/geom/incircle.hpp | 20 +++---- .../geom/include/psemek/geom/intersection.hpp | 48 ++++++++++------- libs/geom/include/psemek/geom/orientation.hpp | 51 ++++++++---------- libs/geom/include/psemek/geom/robust.hpp | 12 +++++ 16 files changed, 232 insertions(+), 139 deletions(-) create mode 100644 libs/geom/include/psemek/geom/robust.hpp diff --git a/libs/cg/include/psemek/cg/convex/inside.hpp b/libs/cg/include/psemek/cg/convex/inside.hpp index f593356c..633ae91f 100644 --- a/libs/cg/include/psemek/cg/convex/inside.hpp +++ b/libs/cg/include/psemek/cg/convex/inside.hpp @@ -7,18 +7,24 @@ namespace psemek::cg { - template - bool inside(geom::point const & p, Body const & body) + template + bool inside(RobustTag robust_tag, geom::point const & p, Body const & body) { auto const & vs = vertices(body); auto const & fs = faces(body); for (auto const & f : fs) { - if (geom::orientation(vs[f[0]], vs[f[1]], vs[f[2]], p) == geom::sign_t::negative) + if (geom::orientation(robust_tag, vs[f[0]], vs[f[1]], vs[f[2]], p) == geom::sign_t::negative) return false; } return true; } + template + bool inside(geom::point const & p, Body const & body) + { + return inside(geom::default_robust_tag, p, body); + } + } diff --git a/libs/cg/include/psemek/cg/convex_hull_2d/andrew.hpp b/libs/cg/include/psemek/cg/convex_hull_2d/andrew.hpp index 2e10c7f7..b1a53dee 100644 --- a/libs/cg/include/psemek/cg/convex_hull_2d/andrew.hpp +++ b/libs/cg/include/psemek/cg/convex_hull_2d/andrew.hpp @@ -9,8 +9,8 @@ namespace psemek::cg { - template - OutIterator andrew_convex_hull(InIterator begin, InIterator end, OutIterator out) + template + OutIterator andrew_convex_hull(RobustTag robust_tag, InIterator begin, InIterator end, OutIterator out) { // need to store iterators to sort them std::vector its(end - begin); @@ -25,7 +25,7 @@ namespace psemek::cg // find lower half of the hull for (auto it = its.begin() + 1; it != its.end(); ++it) { - while (hull.size() >= 2 && orientation(**(hull.end() - 2), **(hull.end() - 1), **it) != geom::sign_t::positive) + while (hull.size() >= 2 && orientation(robust_tag, **(hull.end() - 2), **(hull.end() - 1), **it) != geom::sign_t::positive) hull.pop_back(); hull.push_back(*it); @@ -37,7 +37,7 @@ namespace psemek::cg if (*it == *(hull.end() - 2)) continue; - while (hull.size() >= 2 && orientation(**(hull.end() - 2), **(hull.end() - 1), **it) != geom::sign_t::positive) + while (hull.size() >= 2 && orientation(robust_tag, **(hull.end() - 2), **(hull.end() - 1), **it) != geom::sign_t::positive) hull.pop_back(); hull.push_back(*it); @@ -50,4 +50,10 @@ namespace psemek::cg return std::copy(hull.begin(), hull.end(), out); } + template + OutIterator andrew_convex_hull(InIterator begin, InIterator end, OutIterator out) + { + return andrew_convex_hull(geom::default_robust_tag, begin, end, out); + } + } diff --git a/libs/cg/include/psemek/cg/convex_hull_2d/graham.hpp b/libs/cg/include/psemek/cg/convex_hull_2d/graham.hpp index 9bf4f927..db2d38b0 100644 --- a/libs/cg/include/psemek/cg/convex_hull_2d/graham.hpp +++ b/libs/cg/include/psemek/cg/convex_hull_2d/graham.hpp @@ -9,8 +9,8 @@ namespace psemek::cg { - template - OutIterator graham_convex_hull(InIterator begin, InIterator end, OutIterator out) + template + OutIterator graham_convex_hull(RobustTag robust_tag, InIterator begin, InIterator end, OutIterator out) { // need to store iterators to sort them std::vector its(end - begin); @@ -24,7 +24,7 @@ namespace psemek::cg // sort with respect to angle to leftmost point std::sort(its.begin() + 1, its.end(), [&](auto i1, auto i2){ - auto o = orientation(*its.front(), *i1, *i2); + auto o = orientation(robust_tag, *its.front(), *i1, *i2); // carefully deal with parallel points if (o == geom::sign_t::positive) @@ -41,7 +41,7 @@ namespace psemek::cg for (auto it = its.begin() + 1; it != its.end(); ++it) { - while (hull_end - its.begin() >= 2 && orientation(**(hull_end - 2), **(hull_end - 1), **it) != geom::sign_t::positive) + while (hull_end - its.begin() >= 2 && orientation(robust_tag, **(hull_end - 2), **(hull_end - 1), **it) != geom::sign_t::positive) --hull_end; *hull_end++ = *it; @@ -51,4 +51,10 @@ namespace psemek::cg return std::copy(its.begin(), hull_end, out); } + template + OutIterator graham_convex_hull(InIterator begin, InIterator end, OutIterator out) + { + return graham_convex_hull(geom::default_robust_tag, begin, end, out); + } + } diff --git a/libs/cg/include/psemek/cg/convex_hull_2d/jarvis.hpp b/libs/cg/include/psemek/cg/convex_hull_2d/jarvis.hpp index 36d30ecc..a44888dc 100644 --- a/libs/cg/include/psemek/cg/convex_hull_2d/jarvis.hpp +++ b/libs/cg/include/psemek/cg/convex_hull_2d/jarvis.hpp @@ -5,8 +5,8 @@ namespace psemek::cg { - template - OutIterator jarvis_convex_hull(InIterator begin, InIterator end, OutIterator out) + template + OutIterator jarvis_convex_hull(RobustTag robust_tag, InIterator begin, InIterator end, OutIterator out) { auto first_hull_point = std::min_element(begin, end); auto last_hull_point = first_hull_point; @@ -26,7 +26,7 @@ namespace psemek::cg if (jt == it) continue; if (jt == last_hull_point) continue; - if (orientation(*last_hull_point, *it, *jt) != geom::sign_t::positive) + if (orientation(robust_tag, *last_hull_point, *it, *jt) != geom::sign_t::positive) { is_hull_edge = false; break; @@ -48,4 +48,10 @@ namespace psemek::cg return out; } + template + OutIterator jarvis_convex_hull(InIterator begin, InIterator end, OutIterator out) + { + return jarvis_convex_hull(geom::default_robust_tag, begin, end, out); + } + } diff --git a/libs/cg/include/psemek/cg/convex_hull_2d/quickhull.hpp b/libs/cg/include/psemek/cg/convex_hull_2d/quickhull.hpp index 61bdfa84..0c08f9ea 100644 --- a/libs/cg/include/psemek/cg/convex_hull_2d/quickhull.hpp +++ b/libs/cg/include/psemek/cg/convex_hull_2d/quickhull.hpp @@ -12,8 +12,8 @@ namespace psemek::cg namespace detail { - template - OutIterator quick_hull_recursive_helper(InIterator p1, InIterator p2, ItIterator begin, ItIterator end, OutIterator out) + template + OutIterator quick_hull_recursive_helper(RobustTag robust_tag, InIterator p1, InIterator p2, ItIterator begin, ItIterator end, OutIterator out) { if (begin == end) return *out++ = p1; @@ -21,22 +21,22 @@ namespace psemek::cg // find point in [begin, end) furthest from segment (p1,p2) auto mid = *std::max_element(begin, end, [&](auto it1, auto it2){ // TODO: design a robust predicate - return orientation(*it2, (*it2) + ((*p2) - (*p1)), *it1) == geom::sign_t::positive; + return orientation(robust_tag, *it2, (*it2) + ((*p2) - (*p1)), *it1) == geom::sign_t::positive; }); - auto end1 = std::partition(begin, end, [&](auto it){ return orientation(*p1, *it, *mid) == geom::sign_t::positive; }); - auto end2 = std::partition(end1, end, [&](auto it){ return orientation(*mid, *it, *p2) == geom::sign_t::positive; }); + auto end1 = std::partition(begin, end, [&](auto it){ return orientation(robust_tag, *p1, *it, *mid) == geom::sign_t::positive; }); + auto end2 = std::partition(end1, end, [&](auto it){ return orientation(robust_tag, *mid, *it, *p2) == geom::sign_t::positive; }); - out = quick_hull_recursive_helper(p1, mid, begin, end1, out); - out = quick_hull_recursive_helper(mid, p2, end1, end2, out); + out = quick_hull_recursive_helper(robust_tag, p1, mid, begin, end1, out); + out = quick_hull_recursive_helper(robust_tag, mid, p2, end1, end2, out); return out; } } - template - OutIterator quick_hull(InIterator begin, InIterator end, OutIterator out) + template + OutIterator quick_hull(RobustTag robust_tag, InIterator begin, InIterator end, OutIterator out) { // need to store iterators to move sets of points around std::vector its(end - begin); @@ -52,16 +52,22 @@ namespace psemek::cg { auto p = std::find_if(its.begin() + 1, its.end(), [&](auto it){ return std::all_of(its.begin() + 1, its.end(), [&](auto jt){ - return it == jt || orientation(*its.front(), *it, *jt) == geom::sign_t::negative; + return it == jt || orientation(robust_tag, *its.front(), *it, *jt) == geom::sign_t::negative; }); }); std::iter_swap(its.begin() + 1, p); } // everything is set up, do the recursion - out = detail::quick_hull_recursive_helper(its[0], its[1], its.begin() + 2, its.end(), out); + out = detail::quick_hull_recursive_helper(robust_tag, its[0], its[1], its.begin() + 2, its.end(), out); return *out++ = *(its.begin() + 1); } + template + OutIterator quick_hull(InIterator begin, InIterator end, OutIterator out) + { + return quick_hull(geom::default_robust_tag, begin, end, out); + } + } diff --git a/libs/cg/include/psemek/cg/convex_hull_3d/quickhull.hpp b/libs/cg/include/psemek/cg/convex_hull_3d/quickhull.hpp index 12251022..c85a2742 100644 --- a/libs/cg/include/psemek/cg/convex_hull_3d/quickhull.hpp +++ b/libs/cg/include/psemek/cg/convex_hull_3d/quickhull.hpp @@ -15,8 +15,8 @@ namespace psemek::cg { - template - std::vector> quickhull (std::vector> const & points) + template + std::vector> quickhull(RobustTag robust_tag, std::vector> const & points) { struct face; @@ -72,7 +72,7 @@ namespace psemek::cg return dist012(p1) < dist012(p2); }) - points.begin(); - if (geom::orientation(points[i0], points[i1], points[i2], points[i3]) == geom::sign_t::negative) + if (geom::orientation(robust_tag, points[i0], points[i1], points[i2], points[i3]) == geom::sign_t::negative) std::swap(i2, i3); triangle t0 { i0, i1, i2 }; @@ -89,13 +89,13 @@ namespace psemek::cg if (i == i2) continue; if (i == i3) continue; - if (geom::orientation(points[t0[0]], points[t0[1]], points[t0[2]], points[i]) == geom::sign_t::negative) + if (geom::orientation(robust_tag, points[t0[0]], points[t0[1]], points[t0[2]], points[i]) == geom::sign_t::negative) out0.push_back(i); - else if (geom::orientation(points[t1[0]], points[t1[1]], points[t1[2]], points[i]) == geom::sign_t::negative) + else if (geom::orientation(robust_tag, points[t1[0]], points[t1[1]], points[t1[2]], points[i]) == geom::sign_t::negative) out1.push_back(i); - else if (geom::orientation(points[t2[0]], points[t2[1]], points[t2[2]], points[i]) == geom::sign_t::negative) + else if (geom::orientation(robust_tag, points[t2[0]], points[t2[1]], points[t2[2]], points[i]) == geom::sign_t::negative) out2.push_back(i); - else if (geom::orientation(points[t3[0]], points[t3[1]], points[t3[2]], points[i]) == geom::sign_t::negative) + else if (geom::orientation(robust_tag, points[t3[0]], points[t3[1]], points[t3[2]], points[i]) == geom::sign_t::negative) out3.push_back(i); } @@ -171,7 +171,7 @@ namespace psemek::cg { handle it = h->neigh[i]; triangle const & t = it->tri; - if (it->valid && (geom::orientation(points[t[0]], points[t[1]], points[t[2]], points[max_i]) == geom::sign_t::negative)) + if (it->valid && (geom::orientation(robust_tag, points[t[0]], points[t[1]], points[t[2]], points[max_i]) == geom::sign_t::negative)) { it->valid = false; visible_queue.push(it); @@ -222,7 +222,7 @@ namespace psemek::cg triangle const t { h->tri[(i+1)%3], h->tri[i], max_i }; auto it = std::partition(unassigned.begin(), unassigned_end, [&](Index i){ - return geom::orientation(points[t[0]], points[t[1]], points[t[2]], points[i]) == geom::sign_t::positive; + return geom::orientation(robust_tag, points[t[0]], points[t[1]], points[t[2]], points[i]) == geom::sign_t::positive; }); std::vector out_set(it, unassigned_end); @@ -252,4 +252,10 @@ namespace psemek::cg return result; } + template + auto quickhull(std::vector> const & points) + { + return quickhull(geom::default_robust_tag, points); + } + } diff --git a/libs/cg/include/psemek/cg/triangulation/delaunay.hpp b/libs/cg/include/psemek/cg/triangulation/delaunay.hpp index e21fa522..c664401e 100644 --- a/libs/cg/include/psemek/cg/triangulation/delaunay.hpp +++ b/libs/cg/include/psemek/cg/triangulation/delaunay.hpp @@ -8,8 +8,8 @@ namespace psemek::cg { - template - auto delaunay(InputIterator begin, InputIterator end) + template + auto delaunay(RobustTag robust_tag, InputIterator begin, InputIterator end) { std::vector edge_queue; @@ -58,7 +58,7 @@ namespace psemek::cg auto p3 = e.twin().next().next().origin(); // decide if a flip is needed - if (in_circle(*p0.data(), *p1.data(), *p2.data(), *p3.data()) != geom::sign_t::positive) continue; + if (in_circle(robust_tag, *p0.data(), *p1.data(), *p2.data(), *p3.data()) != geom::sign_t::positive) continue; auto f0 = e.face(); auto f1 = twin.face(); @@ -97,7 +97,13 @@ namespace psemek::cg push(tprev); } }; - return detail::triangulate(begin, end, callback); + return detail::triangulate(robust_tag, begin, end, callback); + } + + template + auto delaunay(InputIterator begin, InputIterator end) + { + return delaunay(geom::default_robust_tag, begin, end); } } diff --git a/libs/cg/include/psemek/cg/triangulation/ear_clipping.hpp b/libs/cg/include/psemek/cg/triangulation/ear_clipping.hpp index afde9669..12eada8c 100644 --- a/libs/cg/include/psemek/cg/triangulation/ear_clipping.hpp +++ b/libs/cg/include/psemek/cg/triangulation/ear_clipping.hpp @@ -9,8 +9,8 @@ namespace psemek::cg { - template - void ear_clipping(Iterator begin, Iterator end, OutputIterator out) + template + OutputIterator ear_clipping(RobustTag robust_tag, Iterator begin, Iterator end, OutputIterator out) { IndexType const count = std::distance(begin, end); @@ -35,7 +35,7 @@ namespace psemek::cg geom::simplex triangle{*(begin + j), *(begin + l), *(begin + h)}; - if (geom::orientation(triangle[0], triangle[1], triangle[2]) == geom::sign_t::negative) + if (geom::orientation(robust_tag, triangle[0], triangle[1], triangle[2]) == geom::sign_t::negative) continue; for (IndexType k = next[h]; k != j; k = next[k]) @@ -61,14 +61,28 @@ namespace psemek::cg break; } } + + return out; } - template - std::vector ear_clipping(Iterator begin, Iterator end) + template + OutputIterator ear_clipping(Iterator begin, Iterator end, OutputIterator out) + { + return ear_clipping(geom::default_robust_tag, begin, end, out); + } + + template + std::vector ear_clipping(RobustTag robust_tag, Iterator begin, Iterator end) { std::vector result; - ear_clipping>, IndexType>(begin, end, std::back_inserter(result)); + ear_clipping(robust_tag, begin, end, std::back_inserter(result)); return result; } + template + auto ear_clipping(Iterator begin, Iterator end) + { + return ear_clipping(geom::default_robust_tag, begin, end); + } + } diff --git a/libs/cg/include/psemek/cg/triangulation/monotone.hpp b/libs/cg/include/psemek/cg/triangulation/monotone.hpp index 10d70c01..aa52fec3 100644 --- a/libs/cg/include/psemek/cg/triangulation/monotone.hpp +++ b/libs/cg/include/psemek/cg/triangulation/monotone.hpp @@ -10,8 +10,8 @@ namespace psemek::cg { - template - void monotone_triangulation(Iterator begin, Iterator end, OutputIterator out) + template + OutputIterator monotone_triangulation(RobustTag robust_tag, Iterator begin, Iterator end, OutputIterator out) { IndexType const count = std::distance(begin, end); IndexType const null = IndexType(-1); @@ -62,23 +62,23 @@ namespace psemek::cg auto const j0 = at(j); auto const j1 = at((j + 1) % count); - auto const oi0 = geom::orientation(i0, i1, j0); - auto const oi1 = geom::orientation(i0, i1, j1); + auto const oi0 = geom::orientation(RobustTag{}, i0, i1, j0); + auto const oi1 = geom::orientation(RobustTag{}, i0, i1, j1); if (oi0 == oi1) return oi0 == geom::sign_t::positive; - return geom::orientation(j0, j1, i0) == geom::sign_t::negative; + return geom::orientation(RobustTag{}, j0, j1, i0) == geom::sign_t::negative; } bool operator()(IndexType e, vertex_id const & v) const { - return geom::orientation(at(e), at((e + 1) % count), at(v.index)) == geom::sign_t::positive; + return geom::orientation(RobustTag{}, at(e), at((e + 1) % count), at(v.index)) == geom::sign_t::positive; } bool operator()(vertex_id const & v, IndexType e) const { - return geom::orientation(at(e), at((e + 1) % count), at(v.index)) == geom::sign_t::negative; + return geom::orientation(RobustTag{}, at(e), at((e + 1) % count), at(v.index)) == geom::sign_t::negative; } }; @@ -110,7 +110,7 @@ namespace psemek::cg bool lp = at(p) < at(e); bool ln = at(n) < at(e); - auto s = geom::orientation(at(p), at(e), at(n)); + auto s = geom::orientation(robust_tag, at(p), at(e), at(n)); if (!lp && !ln && s == geom::sign_t::positive) return event_type::start; @@ -135,10 +135,10 @@ namespace psemek::cg auto f = at(finish(ni)); auto p = at(start[prev[ni]]); bool contains; - if (geom::orientation(f, vi, p) == geom::sign_t::negative) - contains = geom::orientation(f, vi, vj) == geom::sign_t::negative && geom::orientation(vj, vi, p) == geom::sign_t::negative; + if (geom::orientation(robust_tag, f, vi, p) == geom::sign_t::negative) + contains = geom::orientation(robust_tag, f, vi, vj) == geom::sign_t::negative && geom::orientation(robust_tag, vj, vi, p) == geom::sign_t::negative; else - contains = geom::orientation(f, vi, vj) == geom::sign_t::negative || geom::orientation(vj, vi, p) == geom::sign_t::negative; + contains = geom::orientation(robust_tag, f, vi, vj) == geom::sign_t::negative || geom::orientation(robust_tag, vj, vi, p) == geom::sign_t::negative; if (contains) break; ni = twin[prev[ni]]; @@ -151,10 +151,10 @@ namespace psemek::cg auto f = at(finish(nj)); auto p = at(start[prev[nj]]); bool contains; - if (geom::orientation(f, vj, p) == geom::sign_t::negative) - contains = geom::orientation(f, vj, vi) == geom::sign_t::negative && geom::orientation(vi, vj, p) == geom::sign_t::negative; + if (geom::orientation(robust_tag, f, vj, p) == geom::sign_t::negative) + contains = geom::orientation(robust_tag, f, vj, vi) == geom::sign_t::negative && geom::orientation(robust_tag, vi, vj, p) == geom::sign_t::negative; else - contains = geom::orientation(f, vj, vi) == geom::sign_t::negative || geom::orientation(vi, vj, p) == geom::sign_t::negative; + contains = geom::orientation(robust_tag, f, vj, vi) == geom::sign_t::negative || geom::orientation(robust_tag, vi, vj, p) == geom::sign_t::negative; if (contains) break; nj = twin[prev[nj]]; @@ -290,7 +290,7 @@ namespace psemek::cg if (left && stack_left) { - while (stack.size() >= 2 && geom::orientation(at(v), at(stack.back()), at(stack[stack.size() - 2])) == geom::sign_t::negative) + while (stack.size() >= 2 && geom::orientation(robust_tag, at(v), at(stack.back()), at(stack[stack.size() - 2])) == geom::sign_t::negative) { *out++ = v; *out++ = stack[stack.size() - 2]; @@ -301,7 +301,7 @@ namespace psemek::cg } else if (!left && !stack_left) { - while (stack.size() >= 2 && geom::orientation(at(v), at(stack.back()), at(stack[stack.size() - 2])) == geom::sign_t::positive) + while (stack.size() >= 2 && geom::orientation(robust_tag, at(v), at(stack.back()), at(stack[stack.size() - 2])) == geom::sign_t::positive) { *out++ = v; *out++ = stack.back(); @@ -347,14 +347,28 @@ namespace psemek::cg stack.clear(); } + + return out; } - template - std::vector monotone_triangulation(Iterator begin, Iterator end) + template + auto monotone_triangulation(Iterator begin, Iterator end, OutputIterator out) + { + return monotone_triangulation(geom::default_robust_tag, begin, end, out); + } + + template + std::vector monotone_triangulation(RobustTag robust_tag, Iterator begin, Iterator end) { std::vector result; - monotone_triangulation>, IndexType>(begin, end, std::back_inserter(result)); + monotone_triangulation(robust_tag, begin, end, std::back_inserter(result)); return result; } + template + auto monotone_triangulation(Iterator begin, Iterator end) + { + return monotone_triangulation(geom::default_robust_tag, begin, end); + } + } diff --git a/libs/cg/include/psemek/cg/triangulation/triangulation.hpp b/libs/cg/include/psemek/cg/triangulation/triangulation.hpp index e0933a46..4cb4afd7 100644 --- a/libs/cg/include/psemek/cg/triangulation/triangulation.hpp +++ b/libs/cg/include/psemek/cg/triangulation/triangulation.hpp @@ -13,8 +13,8 @@ namespace psemek::cg namespace detail { - template - auto triangulate(Iterator begin, Iterator end, Callback && callback) + template + auto triangulate(RobustTag robust_tag, Iterator begin, Iterator end, Callback && callback) { using point_type = std::decay_t; static_assert(point_type::static_dimension == 2); @@ -93,7 +93,7 @@ namespace psemek::cg bool degenerate = false; // find first hull edge visible from new point - while (orientation(**it, *hp0.data(), *hp1.data()) != geom::sign_t::positive) + while (orientation(robust_tag, **it, *hp0.data(), *hp1.data()) != geom::sign_t::positive) { move_hull_edge(); if (cur_hull_edge == hull_start) @@ -192,7 +192,7 @@ namespace psemek::cg } // until current edge is visible - while (orientation(**it, *hp0.data(), *hp1.data()) == geom::sign_t::positive) + while (orientation(robust_tag, **it, *hp0.data(), *hp1.data()) == geom::sign_t::positive) { // fill new triangle @@ -253,10 +253,16 @@ namespace psemek::cg } + template + auto triangulate(RobustTag robust_tag, InputIterator begin, InputIterator end) + { + return detail::triangulate(robust_tag, begin, end, util::nop); + } + template auto triangulate(InputIterator begin, InputIterator end) { - return detail::triangulate(begin, end, util::nop); + return triangulate(geom::default_robust_tag, begin, end); } } diff --git a/libs/geom/CMakeLists.txt b/libs/geom/CMakeLists.txt index 39105516..af0efa85 100644 --- a/libs/geom/CMakeLists.txt +++ b/libs/geom/CMakeLists.txt @@ -1,19 +1,13 @@ option(PSEMEK_GEOM_ROBUST_PREDICATES "Use robust geometric predicates" OFF) -if(PSEMEK_GEOM_ROBUST_PREDICATES) - find_package(Boost REQUIRED) - find_package(GMP REQUIRED) -endif() +find_package(Boost REQUIRED) +find_package(GMP REQUIRED) file(GLOB_RECURSE PSEMEK_GEOM_HEADERS "include/*.hpp") file(GLOB_RECURSE PSEMEK_GEOM_SOURCES "source/*.cpp") psemek_add_library(psemek-geom ${PSEMEK_GEOM_HEADERS} ${PSEMEK_GEOM_SOURCES}) target_include_directories(psemek-geom PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include") -target_link_libraries(psemek-geom PUBLIC psemek-util) -if(PSEMEK_GEOM_ROBUST_PREDICATES) - target_link_libraries(psemek-geom PUBLIC Boost::boost GMP) - target_compile_definitions(psemek-geom PUBLIC -DPSEMEK_GEOM_ROBUST_PREDICATES=1) -endif() +target_link_libraries(psemek-geom PUBLIC psemek-util Boost::boost GMP) psemek_glob_tests(psemek-geom tests) diff --git a/libs/geom/include/psemek/geom/contains.hpp b/libs/geom/include/psemek/geom/contains.hpp index 579e1fe1..064f48f2 100644 --- a/libs/geom/include/psemek/geom/contains.hpp +++ b/libs/geom/include/psemek/geom/contains.hpp @@ -39,14 +39,20 @@ namespace psemek::geom return true; } - template - bool contains(simplex, 2> const & t, point const & p) + template + bool contains(RobustTag robust_tag, simplex, 2> const & t, point const & p) { return true - && orientation(t[0], t[1], p) != sign_t::negative - && orientation(t[1], t[2], p) != sign_t::negative - && orientation(t[2], t[0], p) != sign_t::negative + && orientation(robust_tag, t[0], t[1], p) != sign_t::negative + && orientation(robust_tag, t[1], t[2], p) != sign_t::negative + && orientation(robust_tag, t[2], t[0], p) != sign_t::negative ; } + template + bool contains(simplex, 2> const & t, point const & p) + { + return contains(default_robust_tag, t, p); + } + } diff --git a/libs/geom/include/psemek/geom/incircle.hpp b/libs/geom/include/psemek/geom/incircle.hpp index ed152fda..604c5dca 100644 --- a/libs/geom/include/psemek/geom/incircle.hpp +++ b/libs/geom/include/psemek/geom/incircle.hpp @@ -3,6 +3,7 @@ #include #include #include +#include #ifdef PSEMEK_GEOM_ROBUST_PREDICATES #include @@ -14,12 +15,7 @@ namespace psemek::geom { template -#ifdef PSEMEK_GEOM_ROBUST_PREDICATES - std::enable_if_t, sign_t> -#else - sign_t -#endif - in_circle(point const & p0, Points const & ... points) + sign_t in_circle(fast_predicate_tag, point const & p0, Points const & ... points) { auto proj = [](vector const & v) { @@ -33,17 +29,16 @@ namespace psemek::geom return u; }; - auto result = orientation(proj(points - p0) ...); + auto result = orientation(fast, proj(points - p0) ...); if constexpr ((N % 2) == 0) return inverse(result); else return result; } -#ifdef PSEMEK_GEOM_ROBUST_PREDICATES template std::enable_if_t, sign_t> - in_circle(point const & p0, point const & p1, point const & p2, point const & p3) + in_circle(robust_predicate_tag, point const & p0, point const & p1, point const & p2, point const & p3) { constexpr T error = std::numeric_limits::epsilon() * T(29) / T(2); @@ -83,6 +78,11 @@ namespace psemek::geom return in_circle(cast(p0), cast(p1), cast(p2), cast(p3)); } } -#endif + + template + sign_t in_circle(Args const & ... args) + { + return in_circle(default_robust_tag, args...); + } } diff --git a/libs/geom/include/psemek/geom/intersection.hpp b/libs/geom/include/psemek/geom/intersection.hpp index 916a8da7..4c28d24b 100644 --- a/libs/geom/include/psemek/geom/intersection.hpp +++ b/libs/geom/include/psemek/geom/intersection.hpp @@ -18,15 +18,21 @@ namespace psemek::geom { + template + bool intersect(RobustTag robust_tag, segment> const & s0, segment> const & s1) + { + auto const o00 = orientation(robust_tag, s0[0], s0[1], s1[0]); + auto const o01 = orientation(robust_tag, s0[0], s0[1], s1[1]); + auto const o10 = orientation(robust_tag, s1[0], s1[1], s0[0]); + auto const o11 = orientation(robust_tag, s1[0], s1[1], s0[1]); + + return ((o00 != o01) || (o00 == sign_t::zero)) && ((o10 != o11) || (o10 == sign_t::zero)); + } + template bool intersect(segment> const & s0, segment> const & s1) { - auto const o00 = orientation(s0[0], s0[1], s1[0]); - auto const o01 = orientation(s0[0], s0[1], s1[1]); - auto const o10 = orientation(s1[0], s1[1], s0[0]); - auto const o11 = orientation(s1[0], s1[1], s0[1]); - - return ((o00 != o01) || (o00 == sign_t::zero)) && ((o10 != o11) || (o10 == sign_t::zero)); + return intersect(default_robust_tag, s0, s1); } // TODO: robust implementation @@ -91,21 +97,27 @@ namespace psemek::geom } } + template + bool intersect(RobustTag robust_tag, triangle> const & t0, triangle> const & t1) + { + if (contains(robust_tag, t0, t1[0]) || contains(robust_tag, t0, t1[1]) || contains(robust_tag, t0, t1[2])) return true; + + if (contains(robust_tag, t1, t0[0]) || contains(robust_tag, t1, t0[1]) || contains(robust_tag, t1, t0[2])) return true; + + if (intersect(robust_tag, simplex{t0[0], t0[1]}, simplex{t1[0], t1[1]})) return true; + if (intersect(robust_tag, simplex{t0[0], t0[1]}, simplex{t1[1], t1[2]})) return true; + if (intersect(robust_tag, simplex{t0[1], t0[2]}, simplex{t1[0], t1[1]})) return true; + if (intersect(robust_tag, simplex{t0[1], t0[2]}, simplex{t1[1], t1[2]})) return true; + if (intersect(robust_tag, simplex{t0[2], t0[0]}, simplex{t1[0], t1[1]})) return true; + if (intersect(robust_tag, simplex{t0[2], t0[0]}, simplex{t1[1], t1[2]})) return true; + + return false; + } + template bool intersect(triangle> const & t0, triangle> const & t1) { - if (contains(t0, t1[0]) || contains(t0, t1[1]) || contains(t0, t1[2])) return true; - - if (contains(t1, t0[0]) || contains(t1, t0[1]) || contains(t1, t0[2])) return true; - - if (intersect(simplex{t0[0], t0[1]}, simplex{t1[0], t1[1]})) return true; - if (intersect(simplex{t0[0], t0[1]}, simplex{t1[1], t1[2]})) return true; - if (intersect(simplex{t0[1], t0[2]}, simplex{t1[0], t1[1]})) return true; - if (intersect(simplex{t0[1], t0[2]}, simplex{t1[1], t1[2]})) return true; - if (intersect(simplex{t0[2], t0[0]}, simplex{t1[0], t1[1]})) return true; - if (intersect(simplex{t0[2], t0[0]}, simplex{t1[1], t1[2]})) return true; - - return false; + return intersect(default_robust_tag, t0, t1); } template diff --git a/libs/geom/include/psemek/geom/orientation.hpp b/libs/geom/include/psemek/geom/orientation.hpp index 63303bdb..9b1ee544 100644 --- a/libs/geom/include/psemek/geom/orientation.hpp +++ b/libs/geom/include/psemek/geom/orientation.hpp @@ -3,10 +3,9 @@ #include #include #include +#include -#ifdef PSEMEK_GEOM_ROBUST_PREDICATES #include -#endif #include #include @@ -14,10 +13,22 @@ namespace psemek::geom { -#ifdef PSEMEK_GEOM_ROBUST_PREDICATES + template + sign_t orientation(fast_predicate_tag, vector const & v1, Vectors const & ... vs) + { + T const d = det(v1, vs...); + + if (d > T{}) + return sign_t::positive; + else if (d < T{}) + return sign_t::negative; + else + return sign_t::zero; + } + template std::enable_if_t, sign_t> - orientation(point const & p0, point const & p1, point const & p2) + orientation(robust_predicate_tag, point const & p0, point const & p1, point const & p2) { constexpr T error = std::numeric_limits::epsilon() * T(5) / T(2); @@ -38,33 +49,9 @@ namespace psemek::geom return orientation(cast(p0), cast(p1), cast(p2)); } } -#endif - - template -#ifdef PSEMEK_GEOM_ROBUST_PREDICATES - std::enable_if_t, sign_t> -#else - sign_t -#endif - orientation(vector const & v1, Vectors const & ... vs) - { - T const d = det(v1, vs...); - - if (d > T{}) - return sign_t::positive; - else if (d < T{}) - return sign_t::negative; - else - return sign_t::zero; - } template -#ifdef PSEMEK_GEOM_ROBUST_PREDICATES - std::enable_if_t, sign_t> -#else - sign_t -#endif - orientation(point const & p0, Points const & ... ps) + sign_t orientation(fast_predicate_tag, point const & p0, Points const & ... ps) { auto o = orientation((ps - p0)...); @@ -78,4 +65,10 @@ namespace psemek::geom } } + template + sign_t orientation(Args const & ... args) + { + return orientation(default_robust_tag, args...); + } + } diff --git a/libs/geom/include/psemek/geom/robust.hpp b/libs/geom/include/psemek/geom/robust.hpp new file mode 100644 index 00000000..64fe8a04 --- /dev/null +++ b/libs/geom/include/psemek/geom/robust.hpp @@ -0,0 +1,12 @@ +#pragma once + +namespace psemek::geom +{ + + constexpr struct fast_predicate_tag{} fast; + + constexpr struct robust_predicate_tag{} robust; + + constexpr auto default_robust_tag = fast; + +}