Add monotone chain simple polygon triangulation

This commit is contained in:
Nikita Lisitsa 2022-09-15 14:42:13 +03:00
parent 2e44978ff3
commit b2a714f88a

View file

@ -0,0 +1,360 @@
#pragma once
#include <psemek/geom/point.hpp>
#include <psemek/geom/orientation.hpp>
#include <psemek/geom/contains.hpp>
#include <vector>
#include <set>
namespace psemek::cg
{
template <typename Iterator, typename OutputIterator, typename IndexType = std::size_t>
void monotone_triangulation(Iterator begin, Iterator end, OutputIterator out)
{
IndexType const count = std::distance(begin, end);
IndexType const null = IndexType(-1);
auto at = [begin](auto i){ return *(begin + i); };
// Initialize dcel
std::vector<IndexType> start(count, 0);
std::vector<IndexType> next(count, 0);
std::vector<IndexType> prev(count, 0);
std::vector<IndexType> twin(count, 0);
for (IndexType i = 0; i < count; ++i)
{
start[i] = i;
next[i] = (i + 1) % count;
prev[i] = (i + count - 1) % count;
twin[i] = null;
}
auto finish = [&](auto i)
{
return start[next[i]];
};
struct vertex_id
{
IndexType index;
};
struct sweepline_compare
{
struct is_transparent{};
Iterator begin;
IndexType count;
auto const & at(IndexType i) const
{
return *(begin + i);
}
bool operator()(IndexType i, IndexType j) const
{
auto const i0 = at(i);
auto const i1 = at((i + 1) % count);
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);
if (oi0 == oi1)
return oi0 == geom::sign_t::positive;
return geom::orientation(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;
}
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;
}
};
// Edges intersecting sweepline
std::set<IndexType, sweepline_compare> sweepline(sweepline_compare{begin, count});
std::vector<IndexType> events(count, 0);
std::vector<IndexType> helper(count, 0);
for (IndexType i = 0; i < count; ++i)
events[i] = i;
std::sort(events.begin(), events.end(), [&](auto i, auto j){ return at(i) < at(j); });
enum class event_type
{
start,
end,
split,
merge,
regular,
};
auto classify = [&](auto e) -> event_type
{
auto p = (e + count - 1) % count;
auto n = (e + 1) % count;
bool lp = at(p) < at(e);
bool ln = at(n) < at(e);
auto s = geom::orientation(at(p), at(e), at(n));
if (!lp && !ln && s == geom::sign_t::positive)
return event_type::start;
else if (lp && ln && s == geom::sign_t::positive)
return event_type::end;
else if (!lp && !ln && s == geom::sign_t::negative)
return event_type::split;
else if (lp && ln && s == geom::sign_t::negative)
return event_type::merge;
else
return event_type::regular;
};
auto connect = [&](auto i, auto j)
{
auto const vi = at(i);
auto const vj = at(j);
auto ni = i;
while (true)
{
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;
else
contains = geom::orientation(f, vi, vj) == geom::sign_t::negative || geom::orientation(vj, vi, p) == geom::sign_t::negative;
if (contains)
break;
ni = twin[prev[ni]];
}
auto pi = prev[ni];
auto nj = j;
while (true)
{
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;
else
contains = geom::orientation(f, vj, vi) == geom::sign_t::negative || geom::orientation(vi, vj, p) == geom::sign_t::negative;
if (contains)
break;
nj = twin[prev[nj]];
}
auto pj = prev[nj];
IndexType eij = start.size();
IndexType eji = eij + 1;
start.emplace_back();
start.emplace_back();
next.emplace_back();
next.emplace_back();
prev.emplace_back();
prev.emplace_back();
twin.emplace_back();
twin.emplace_back();
start[eij] = i;
next[eij] = nj;
prev[eij] = pi;
twin[eij] = eji;
start[eji] = j;
next[eji] = ni;
prev[eji] = pj;
twin[eji] = eij;
next[pi] = eij;
prev[nj] = eij;
next[pj] = eji;
prev[ni] = eji;
};
// Partition into monotone parts
for (auto e : events)
{
auto p = (e + count - 1) % count;
switch (classify(e))
{
case event_type::start:
sweepline.insert(e);
helper[e] = e;
break;
case event_type::end:
if (auto h = helper[p]; classify(h) == event_type::merge)
connect(e, h);
sweepline.erase(p);
break;
case event_type::split:
{
auto l = *std::prev(sweepline.upper_bound(vertex_id{e}));
connect(e, helper[l]);
helper[l] = e;
}
sweepline.insert(e);
helper[e] = e;
break;
case event_type::merge:
if (auto h = helper[p]; classify(h) == event_type::merge)
connect(e, h);
sweepline.erase(p);
{
auto l = *std::prev(sweepline.upper_bound(vertex_id{e}));
if (auto h = helper[l]; classify(h) == event_type::merge)
connect(e, h);
helper[l] = e;
}
break;
case event_type::regular:
if (at(e) < at((e + 1) % count))
{
if (auto h = helper[p]; classify(h) == event_type::merge)
connect(e, h);
sweepline.erase(p);
sweepline.insert(e);
helper[e] = e;
}
else
{
auto l = *std::prev(sweepline.upper_bound(vertex_id{e}));
if (auto h = helper[l]; classify(h) == event_type::merge)
connect(e, h);
helper[l] = e;
}
break;
}
}
// Triangulate monotone parts
std::vector<bool> used(start.size(), false);
std::vector<IndexType> stack;
for (IndexType e = 0; e < start.size(); ++e)
{
if (used[e]) continue;
auto le = e;
used[e] = true;
// Mark the whole part as used and find leftmost point
for (auto ee = next[e]; ee != e; ee = next[ee])
{
used[ee] = true;
if (at(start[ee]) < at(start[le]))
le = ee;
}
bool stack_left = true;
auto re = prev[le];
stack.push_back(start[le]);
le = next[le];
while (true)
{
bool last = (le == re);
IndexType v;
bool left;
if (at(start[le]) < at(start[re]))
{
v = start[le];
le = next[le];
left = true;
}
else
{
v = start[re];
re = prev[re];
left = false;
}
if (left && stack_left)
{
while (stack.size() >= 2 && geom::orientation(at(v), at(stack.back()), at(stack[stack.size() - 2])) == geom::sign_t::negative)
{
*out++ = v;
*out++ = stack[stack.size() - 2];
*out++ = stack.back();
stack.pop_back();
}
stack.push_back(v);
}
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)
{
*out++ = v;
*out++ = stack.back();
*out++ = stack[stack.size() - 2];
stack.pop_back();
}
stack.push_back(v);
}
else if (left && !stack_left)
{
auto last = stack.back();
while (stack.size() >= 2)
{
*out++ = v;
*out++ = stack.back();
*out++ = stack[stack.size() - 2];
stack.pop_back();
}
stack.pop_back();
stack.push_back(last);
stack.push_back(v);
}
else if (!left && stack_left)
{
auto last = stack.back();
while (stack.size() >= 2)
{
*out++ = v;
*out++ = stack[stack.size() - 2];
*out++ = stack.back();
stack.pop_back();
}
stack.pop_back();
stack.push_back(last);
stack.push_back(v);
}
stack_left = left;
if (last)
break;
}
stack.clear();
}
}
template <typename IndexType, typename Iterator>
std::vector<IndexType> monotone_triangulation(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));
return result;
}
}