lisyarus/psemek
1
0
Fork 1
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Switch from define-based to tag-based robust predicates switch

Browse source
This commit is contained in:
Nikita Lisitsa 2022-09-16 18:30:39 +03:00
parent 46458abd71
commit 4c976272c5
16 changed files with 232 additions and 139 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

12
libs/cg/include/psemek/cg/convex/inside.hpp
View file

@ -7,18 +7,24 @@
namespace psemek::cg
{
template <typename T, std::size_t N, typename Body>
bool inside(geom::point<T, N> const & p, Body const & body)
template <typename RobustTag, typename T, std::size_t N, typename Body>
bool inside(RobustTag robust_tag, geom::point<T, N> 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 <typename T, std::size_t N, typename Body>
bool inside(geom::point<T, N> const & p, Body const & body)
{
return inside(geom::default_robust_tag, p, body);
}
}

14
libs/cg/include/psemek/cg/convex_hull_2d/andrew.hpp
View file

@ -9,8 +9,8 @@
namespace psemek::cg
{
template <typename InIterator, typename OutIterator>
OutIterator andrew_convex_hull(InIterator begin, InIterator end, OutIterator out)
template <typename RobustTag, typename InIterator, typename OutIterator>
OutIterator andrew_convex_hull(RobustTag robust_tag, InIterator begin, InIterator end, OutIterator out)
{
// need to store iterators to sort them
std::vector<InIterator> 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 <typename InIterator, typename OutIterator>
OutIterator andrew_convex_hull(InIterator begin, InIterator end, OutIterator out)
{
return andrew_convex_hull(geom::default_robust_tag, begin, end, out);
}
}

14
libs/cg/include/psemek/cg/convex_hull_2d/graham.hpp
View file

@ -9,8 +9,8 @@
namespace psemek::cg
{
template <typename InIterator, typename OutIterator>
OutIterator graham_convex_hull(InIterator begin, InIterator end, OutIterator out)
template <typename RobustTag, typename InIterator, typename OutIterator>
OutIterator graham_convex_hull(RobustTag robust_tag, InIterator begin, InIterator end, OutIterator out)
{
// need to store iterators to sort them
std::vector<InIterator> 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 <typename InIterator, typename OutIterator>
OutIterator graham_convex_hull(InIterator begin, InIterator end, OutIterator out)
{
return graham_convex_hull(geom::default_robust_tag, begin, end, out);
}
}

12
libs/cg/include/psemek/cg/convex_hull_2d/jarvis.hpp
View file

@ -5,8 +5,8 @@
namespace psemek::cg
{
template <typename InIterator, typename OutIterator>
OutIterator jarvis_convex_hull(InIterator begin, InIterator end, OutIterator out)
template <typename RobustTag, typename InIterator, typename OutIterator>
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 <typename InIterator, typename OutIterator>
OutIterator jarvis_convex_hull(InIterator begin, InIterator end, OutIterator out)
{
return jarvis_convex_hull(geom::default_robust_tag, begin, end, out);
}
}

28
libs/cg/include/psemek/cg/convex_hull_2d/quickhull.hpp
View file

@ -12,8 +12,8 @@ namespace psemek::cg
namespace detail
{
template <typename InIterator, typename ItIterator, typename OutIterator>
OutIterator quick_hull_recursive_helper(InIterator p1, InIterator p2, ItIterator begin, ItIterator end, OutIterator out)
template <typename RobustTag, typename InIterator, typename ItIterator, typename OutIterator>
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 <typename InIterator, typename OutIterator>
OutIterator quick_hull(InIterator begin, InIterator end, OutIterator out)
template <typename RobustTag, typename InIterator, typename OutIterator>
OutIterator quick_hull(RobustTag robust_tag, InIterator begin, InIterator end, OutIterator out)
{
// need to store iterators to move sets of points around
std::vector<InIterator> 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 <typename InIterator, typename OutIterator>
OutIterator quick_hull(InIterator begin, InIterator end, OutIterator out)
{
return quick_hull(geom::default_robust_tag, begin, end, out);
}
}

24
libs/cg/include/psemek/cg/convex_hull_3d/quickhull.hpp
View file

@ -15,8 +15,8 @@
namespace psemek::cg
{
template <typename T, typename Index = std::size_t>
std::vector<geom::triangle<Index>> quickhull (std::vector<geom::point<T, 3ul>> const & points)
template <typename RobustTag, typename T, typename Index = std::size_t>
std::vector<geom::triangle<Index>> quickhull(RobustTag robust_tag, std::vector<geom::point<T, 3ul>> 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<Index> out_set(it, unassigned_end);
@ -252,4 +252,10 @@ namespace psemek::cg
return result;
}
template <typename T, typename Index = std::size_t>
auto quickhull(std::vector<geom::point<T, 3ul>> const & points)
{
return quickhull(geom::default_robust_tag, points);
}
}

14
libs/cg/include/psemek/cg/triangulation/delaunay.hpp
View file

@ -8,8 +8,8 @@
namespace psemek::cg
{
template <typename Index = std::size_t, typename InputIterator>
auto delaunay(InputIterator begin, InputIterator end)
template <typename Index = std::size_t, typename RobustTag, typename InputIterator>
auto delaunay(RobustTag robust_tag, InputIterator begin, InputIterator end)
{
std::vector<Index> 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<Index>(robust_tag, begin, end, callback);
}
template <typename Index = std::size_t, typename InputIterator>
auto delaunay(InputIterator begin, InputIterator end)
{
return delaunay<Index>(geom::default_robust_tag, begin, end);
}
}

26
libs/cg/include/psemek/cg/triangulation/ear_clipping.hpp
View file

@ -9,8 +9,8 @@
namespace psemek::cg
{
template <typename Iterator, typename OutputIterator, typename IndexType = std::size_t>
void ear_clipping(Iterator begin, Iterator end, OutputIterator out)
template <typename IndexType = std::size_t, typename RobustTag, typename Iterator, typename OutputIterator>
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 <typename IndexType, typename Iterator>
std::vector<IndexType> ear_clipping(Iterator begin, Iterator end)
template <typename IndexType = std::size_t, typename Iterator, typename OutputIterator>
OutputIterator ear_clipping(Iterator begin, Iterator end, OutputIterator out)
{
return ear_clipping<IndexType>(geom::default_robust_tag, begin, end, out);
}
template <typename IndexType, typename RobustTag, typename Iterator>
std::vector<IndexType> ear_clipping(RobustTag robust_tag, Iterator begin, Iterator end)
{
std::vector<IndexType> result;
ear_clipping<Iterator, std::back_insert_iterator<std::vector<IndexType>>, IndexType>(begin, end, std::back_inserter(result));
ear_clipping<IndexType>(robust_tag, begin, end, std::back_inserter(result));
return result;
}
template <typename IndexType, typename Iterator>
auto ear_clipping(Iterator begin, Iterator end)
{
return ear_clipping<IndexType>(geom::default_robust_tag, begin, end);
}
}

52
libs/cg/include/psemek/cg/triangulation/monotone.hpp
View file

@ -10,8 +10,8 @@
namespace psemek::cg
{
template <typename Iterator, typename OutputIterator, typename IndexType = std::size_t>
void monotone_triangulation(Iterator begin, Iterator end, OutputIterator out)
template <typename IndexType = std::size_t, typename RobustTag, typename Iterator, typename OutputIterator>
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 <typename IndexType, typename Iterator>
std::vector<IndexType> monotone_triangulation(Iterator begin, Iterator end)
template <typename IndexType = std::size_t, typename Iterator, typename OutputIterator>
auto monotone_triangulation(Iterator begin, Iterator end, OutputIterator out)
{
return monotone_triangulation<IndexType>(geom::default_robust_tag, begin, end, out);
}
template <typename IndexType, typename RobustTag, typename Iterator>
std::vector<IndexType> monotone_triangulation(RobustTag robust_tag, Iterator begin, Iterator end)
{
std::vector<IndexType> result;
monotone_triangulation<Iterator, std::back_insert_iterator<std::vector<IndexType>>, IndexType>(begin, end, std::back_inserter(result));
monotone_triangulation<IndexType>(robust_tag, begin, end, std::back_inserter(result));
return result;
}
template <typename IndexType, typename Iterator>
auto monotone_triangulation(Iterator begin, Iterator end)
{
return monotone_triangulation<IndexType>(geom::default_robust_tag, begin, end);
}
}

16
libs/cg/include/psemek/cg/triangulation/triangulation.hpp
View file

@ -13,8 +13,8 @@ namespace psemek::cg
namespace detail
{
template <typename Index = std::size_t, typename Iterator, typename Callback>
auto triangulate(Iterator begin, Iterator end, Callback && callback)
template <typename Index, typename RobustTag, typename Iterator, typename Callback>
auto triangulate(RobustTag robust_tag, Iterator begin, Iterator end, Callback && callback)
{
using point_type = std::decay_t<decltype(*begin)>;
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 <typename Index = std::size_t, typename RobustTag, typename InputIterator>
auto triangulate(RobustTag robust_tag, InputIterator begin, InputIterator end)
{
return detail::triangulate<Index>(robust_tag, begin, end, util::nop);
}
template <typename Index = std::size_t, typename InputIterator>
auto triangulate(InputIterator begin, InputIterator end)
{
return detail::triangulate(begin, end, util::nop);
return triangulate<Index>(geom::default_robust_tag, begin, end);
}
}

12
libs/geom/CMakeLists.txt
View file

@ -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)

16
libs/geom/include/psemek/geom/contains.hpp
View file

@ -39,14 +39,20 @@ namespace psemek::geom
return true;
}
template <typename T>
bool contains(simplex<point<T, 2>, 2> const & t, point<T, 2> const & p)
template <typename RobustTag, typename T>
bool contains(RobustTag robust_tag, simplex<point<T, 2>, 2> const & t, point<T, 2> 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 <typename T>
bool contains(simplex<point<T, 2>, 2> const & t, point<T, 2> const & p)
{
return contains(default_robust_tag, t, p);
}
}

20
libs/geom/include/psemek/geom/incircle.hpp
View file

@ -3,6 +3,7 @@
#include <psemek/geom/sign.hpp>
#include <psemek/geom/point.hpp>
#include <psemek/geom/orientation.hpp>
#include <psemek/geom/robust.hpp>
#ifdef PSEMEK_GEOM_ROBUST_PREDICATES
#include <boost/multiprecision/gmp.hpp>
@ -14,12 +15,7 @@ namespace psemek::geom
{
template <typename T, std::size_t N, typename ... Points>
#ifdef PSEMEK_GEOM_ROBUST_PREDICATES
std::enable_if_t<!std::is_floating_point_v<T>, sign_t>
#else
sign_t
#endif
in_circle(point<T, N> const & p0, Points const & ... points)
sign_t in_circle(fast_predicate_tag, point<T, N> const & p0, Points const & ... points)
{
auto proj = [](vector<T, N> 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 <typename T>
std::enable_if_t<std::is_floating_point_v<T>, sign_t>
in_circle(point<T, 2> const & p0, point<T, 2> const & p1, point<T, 2> const & p2, point<T, 2> const & p3)
in_circle(robust_predicate_tag, point<T, 2> const & p0, point<T, 2> const & p1, point<T, 2> const & p2, point<T, 2> const & p3)
{
constexpr T error = std::numeric_limits<T>::epsilon() * T(29) / T(2);
@ -83,6 +78,11 @@ namespace psemek::geom
return in_circle(cast<exact_type>(p0), cast<exact_type>(p1), cast<exact_type>(p2), cast<exact_type>(p3));
}
}
#endif
template <typename ... Args>
sign_t in_circle(Args const & ... args)
{
return in_circle(default_robust_tag, args...);
}
}

48
libs/geom/include/psemek/geom/intersection.hpp
View file

@ -18,15 +18,21 @@
namespace psemek::geom
{
template <typename RobustTag, typename T>
bool intersect(RobustTag robust_tag, segment<point<T, 2>> const & s0, segment<point<T, 2>> 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 <typename T>
bool intersect(segment<point<T, 2>> const & s0, segment<point<T, 2>> 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 <typename RobustTag, typename T>
bool intersect(RobustTag robust_tag, triangle<point<T, 2>> const & t0, triangle<point<T, 2>> 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 <typename T>
bool intersect(triangle<point<T, 2>> const & t0, triangle<point<T, 2>> 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 <typename T, std::size_t N>

51
libs/geom/include/psemek/geom/orientation.hpp
View file

@ -3,10 +3,9 @@
#include <psemek/geom/vector.hpp>
#include <psemek/geom/point.hpp>
#include <psemek/geom/sign.hpp>
#include <psemek/geom/robust.hpp>
#ifdef PSEMEK_GEOM_ROBUST_PREDICATES
#include <boost/multiprecision/gmp.hpp>
#endif
#include <limits>
#include <type_traits>
@ -14,10 +13,22 @@
namespace psemek::geom
{
#ifdef PSEMEK_GEOM_ROBUST_PREDICATES
template <typename T, std::size_t N, typename ... Vectors>
sign_t orientation(fast_predicate_tag, vector<T, N> 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 <typename T>
std::enable_if_t<std::is_floating_point_v<T>, sign_t>
orientation(point<T, 2> const & p0, point<T, 2> const & p1, point<T, 2> const & p2)
orientation(robust_predicate_tag, point<T, 2> const & p0, point<T, 2> const & p1, point<T, 2> const & p2)
{
constexpr T error = std::numeric_limits<T>::epsilon() * T(5) / T(2);
@ -38,33 +49,9 @@ namespace psemek::geom
return orientation(cast<exact_type>(p0), cast<exact_type>(p1), cast<exact_type>(p2));
}
}
#endif
template <typename T, std::size_t N, typename ... Vectors>
#ifdef PSEMEK_GEOM_ROBUST_PREDICATES
std::enable_if_t<!std::is_floating_point_v<T>, sign_t>
#else
sign_t
#endif
orientation(vector<T, N> 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 <typename T, std::size_t N, typename ... Points>
#ifdef PSEMEK_GEOM_ROBUST_PREDICATES
std::enable_if_t<!std::is_floating_point_v<T>, sign_t>
#else
sign_t
#endif
orientation(point<T, N> const & p0, Points const & ... ps)
sign_t orientation(fast_predicate_tag, point<T, N> const & p0, Points const & ... ps)
{
auto o = orientation((ps - p0)...);
@ -78,4 +65,10 @@ namespace psemek::geom
}
}
template <typename ... Args>
sign_t orientation(Args const & ... args)
{
return orientation(default_robust_tag, args...);
}
}

12
libs/geom/include/psemek/geom/robust.hpp Normal file
View file

@ -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;
}