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

Re-implement monotone chain triangulation to output a DCEL

Browse source
This commit is contained in:
Nikita Lisitsa 2022-09-19 23:35:33 +03:00
parent 92ad22f656
commit 7a57c1c1d8
2 changed files with 202 additions and 137 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

56
examples/triangulation.cpp
View file

@ -64,14 +64,20 @@ struct triangulation_app
{
prof::profiler prof("triangulate");
auto dcel = cg::ear_clipping<std::uint16_t>(geom::fast, points_.begin(), points_.end());
// auto dcel = cg::ear_clipping<std::uint16_t>(geom::fast, points_.begin(), points_.end());
auto dcel = cg::monotone_triangulation<std::uint16_t>(geom::fast, points_.begin(), points_.end());
edges_ = cg::edge_mesh(dcel);
triangles_ = cg::triangle_mesh(dcel);
}
prof::dump();
log::info() << (triangles_.size() / 3) << " triangles";
log::info() << edges_.size() << " edges";
log::info() << triangles_.size() << " triangles";
log::info() << "Euler chi: " << (points_.size() - edges_.size() + triangles_.size());
closest_points_.resize(points_.size());
std::iota(closest_points_.begin(), closest_points_.end(), std::uint16_t{0});
}
void on_left_button_down() override
@ -125,8 +131,7 @@ struct triangulation_app
};
random::generator rng;
for (auto const & e : edges_)
edge(e[0], e[1], gfx::black);
for (auto const & t : triangles_)
{
gfx::color_rgba c;
@ -138,14 +143,52 @@ struct triangulation_app
painter_.triangle(points_[t[0]], points_[t[1]], points_[t[2]], c);
}
// for (auto const & e : edges_)
// edge(e[0], e[1], gfx::black);
for (std::size_t i = 0; i < points_.size(); ++i)
{
edge(i, (i + 1) % points_.size(), gfx::black);
}
auto camera_transform = geom::orthographic_camera{view_bbox}.transform();
painter_.render(camera_transform);
if (auto m = mouse(); m)
{
geom::point<float, 2> m_world;
m_world[0] = geom::lerp(view_bbox[0], (*m)[0] * 1.f / width());
m_world[1] = geom::lerp(view_bbox[1], 1.f - (*m)[1] * 1.f / height());
auto compare = [&](auto i, auto j){
return geom::distance(points_[i], m_world) < geom::distance(points_[j], m_world);
};
std::size_t const n_closest = 8;
std::partial_sort(closest_points_.begin(), closest_points_.begin() + n_closest, closest_points_.end(), compare);
float max_distance = camera_size_ / 32.f;
for (std::size_t j = 0; j < n_closest; ++j)
{
auto i = closest_points_[j];
if (geom::distance(points_[i], m_world) > max_distance)
break;
gfx::painter::text_options opts;
opts.c = {0, 0, 0, 255};
opts.x = gfx::painter::x_align::left;
opts.y = gfx::painter::y_align::bottom;
opts.scale = 2.f;
auto p = geom::swizzle<0, 1>(geom::as_point(camera_transform * geom::homogeneous(geom::swizzle<0, 1, -1>(points_[i]))));
p[0] = std::round((p[0] * 0.5f + 0.5f) * width());
p[1] = std::round((0.5f - p[1] * 0.5f) * height());
painter_.text(p, util::to_string(i), opts);
}
}
painter_.render(geom::window_camera{width(), height()}.transform());
}
private:
@ -158,6 +201,7 @@ private:
float camera_size_tgt_;
std::optional<geom::point<int, 2>> drag_start_;
std::vector<std::uint16_t> closest_points_;
gfx::painter painter_;

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

@ -3,6 +3,7 @@
#include <psemek/geom/point.hpp>
#include <psemek/geom/orientation.hpp>
#include <psemek/geom/contains.hpp>
#include <psemek/cg/dcel.hpp>
#include <vector>
#include <set>
@ -10,39 +11,43 @@
namespace psemek::cg
{
template <typename IndexType = std::size_t, typename RobustTag, typename Iterator, typename OutputIterator>
OutputIterator monotone_triangulation(RobustTag robust_tag, Iterator begin, Iterator end, OutputIterator out)
template <typename IndexType = std::size_t, typename RobustTag, typename Iterator>
auto monotone_triangulation(RobustTag robust_tag, Iterator begin, Iterator end)
{
IndexType const count = std::distance(begin, end);
IndexType const null = IndexType(-1);
auto at = [begin](auto i){ return *(begin + i); };
auto at = [begin](IndexType 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);
using dcel_type = dcel<util::empty, util::empty, util::empty, IndexType>;
using point_handle = typename dcel_type::point_handle;
dcel_type result;
result.push_face();
result.faces[0].edge = count;
result.points.resize(count);
result.edges.resize(2 * count);
for (IndexType i = 0; i < count; ++i)
{
start[i] = i;
next[i] = (i + 1) % count;
prev[i] = (i + count - 1) % count;
twin[i] = null;
result.points[i].edge = i;
result.edges[i].face = result.null;
result.edges[i].origin = i;
result.edges[i].prev = (i + count - 1) % count;
result.edges[i].next = (i + 1) % count;
result.edges[i].twin = i + count;
result.edges[i + count].face = 0;
result.edges[i + count].origin = (i + 1) % count;
result.edges[i + count].prev = count + (i + 1) % count;
result.edges[i + count].next = count + (i + count - 1) % count;
result.edges[i + count].twin = i;
}
auto finish = [&](auto i)
{
return start[next[i]];
};
struct vertex_id
{
IndexType index;
};
struct sweepline_compare
{
struct is_transparent{};
@ -71,14 +76,14 @@ namespace psemek::cg
return geom::orientation(RobustTag{}, j0, j1, i0) == geom::sign_t::negative;
}
bool operator()(IndexType e, vertex_id const & v) const
bool operator()(IndexType e, point_handle v) const
{
return geom::orientation(RobustTag{}, 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
bool operator()(point_handle v, IndexType e) const
{
return geom::orientation(RobustTag{}, 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;
}
};
@ -129,11 +134,11 @@ namespace psemek::cg
auto const vi = at(i);
auto const vj = at(j);
auto ni = i;
auto ni = result.edge(i);
while (true)
{
auto f = at(finish(ni));
auto p = at(start[prev[ni]]);
auto f = at(ni.next().origin().index());
auto p = at(ni.prev().origin().index());
bool contains;
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;
@ -141,15 +146,15 @@ namespace psemek::cg
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]];
ni = ni.prev().twin();
}
auto pi = prev[ni];
auto pi = ni.prev();
auto nj = j;
auto nj = result.edge(j);
while (true)
{
auto f = at(finish(nj));
auto p = at(start[prev[nj]]);
auto f = at(nj.next().origin().index());
auto p = at(nj.prev().origin().index());
bool contains;
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;
@ -157,37 +162,22 @@ namespace psemek::cg
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]];
nj = nj.prev().twin();
}
auto pj = prev[nj];
auto pj = nj.prev();
IndexType eij = start.size();
IndexType eji = eij + 1;
auto eij = result.push_edge();
auto eji = result.push_edge();
start.emplace_back();
start.emplace_back();
next.emplace_back();
next.emplace_back();
prev.emplace_back();
prev.emplace_back();
twin.emplace_back();
twin.emplace_back();
eij.origin(result.point(i));
eij.next(nj); nj.prev(eij);
eij.prev(pi); pi.next(eij);
eij.twin(eji);
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;
eji.origin(result.point(j));
eji.next(ni); ni.prev(eji);
eji.prev(pj); pj.next(eji);
eji.twin(eij);
};
// Partition into monotone parts
@ -208,7 +198,7 @@ namespace psemek::cg
break;
case event_type::split:
{
auto l = *std::prev(sweepline.upper_bound(vertex_id{e}));
auto l = *std::prev(sweepline.upper_bound(result.point(e)));
connect(e, helper[l]);
helper[l] = e;
}
@ -220,7 +210,7 @@ namespace psemek::cg
connect(e, h);
sweepline.erase(p);
{
auto l = *std::prev(sweepline.upper_bound(vertex_id{e}));
auto l = *std::prev(sweepline.upper_bound(result.point(e)));
if (auto h = helper[l]; classify(h) == event_type::merge)
connect(e, h);
helper[l] = e;
@ -237,7 +227,7 @@ namespace psemek::cg
}
else
{
auto l = *std::prev(sweepline.upper_bound(vertex_id{e}));
auto l = *std::prev(sweepline.upper_bound(result.point(e)));
if (auto h = helper[l]; classify(h) == event_type::merge)
connect(e, h);
helper[l] = e;
@ -246,125 +236,156 @@ namespace psemek::cg
}
}
// 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;
fill_missing_faces(result);
IndexType face_count = result.faces.size();
// Triangulate monotone parts
std::vector<IndexType> stack;
for (IndexType f = 1; f < face_count; ++f)
{
auto e = result.face(f).edge();
// Find leftmost point
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])
for (auto ee = e.next(); ee != e; ee = ee.next())
{
used[ee] = true;
if (at(start[ee]) < at(start[le]))
if (at(ee.origin().index()) < at(le.origin().index()))
le = ee;
}
bool stack_left = true;
auto re = prev[le];
stack.push_back(start[le]);
le = next[le];
auto re = le.prev();
stack.push_back(le.index());
le = le.next();
auto emit_triangle = [&](auto e0_id, auto e1_id)
{
auto e0 = result.edge(e0_id);
auto e1 = result.edge(e1_id);
auto p = e0.prev();
auto n = e1.next();
auto face_out = e0.face();
auto face_in = result.push_face();
auto e_in = result.push_edge();
auto e_out = result.push_edge();
face_in.edge(e_in);
face_out.edge(e_out);
e_in.origin(n.origin());
e_out.origin(e0.origin());
e_in.twin(e_out);
e_out.twin(e_in);
e0.face(face_in);
e1.face(face_in);
e_in.face(face_in);
e_out.face(face_out);
e0.prev(e_in); e_in.next(e0);
e1.next(e_in); e_in.prev(e1);
p.next(e_out); e_out.prev(p);
n.prev(e_out); e_out.next(n);
return e_out.index();
};
while (true)
{
bool last = (le == re);
IndexType v;
auto ve = le;
bool left;
if (at(start[le]) < at(start[re]))
if (at(le.origin().index()) < at(re.origin().index()))
{
v = start[le];
le = next[le];
ve = le;
le = le.next();
left = true;
}
else
{
v = start[re];
re = prev[re];
ve = re;
re = re.prev();
left = false;
}
auto orientation = [&](auto e0, auto e1, auto e2)
{
auto point = [&](auto e){ return at(result.edge(e).origin().index()); };
return geom::orientation(robust_tag, point(e0), point(e1), point(e2));
};
if (left && stack_left)
{
while (stack.size() >= 2 && geom::orientation(robust_tag, at(v), at(stack.back()), at(stack[stack.size() - 2])) == geom::sign_t::negative)
while (stack.size() >= 2 && orientation(ve.index(), stack.back(), stack[stack.size() - 2]) == geom::sign_t::negative)
{
*out++ = v;
*out++ = stack[stack.size() - 2];
*out++ = stack.back();
auto new_edge = emit_triangle(stack[stack.size() - 2], stack.back());
stack.pop_back();
stack.pop_back();
stack.push_back(new_edge);
}
stack.push_back(v);
stack.push_back(ve.index());
}
else if (!left && !stack_left)
{
while (stack.size() >= 2 && geom::orientation(robust_tag, at(v), at(stack.back()), at(stack[stack.size() - 2])) == geom::sign_t::positive)
auto e = ve.index();
while (stack.size() >= 2 && orientation(ve.index(), stack.back(), stack[stack.size() - 2]) == geom::sign_t::positive)
{
*out++ = v;
*out++ = stack.back();
*out++ = stack[stack.size() - 2];
e = emit_triangle(e, stack.back());
stack.pop_back();
}
stack.push_back(v);
stack.push_back(e);
}
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);
auto e = stack.front();
for (auto it = std::next(stack.begin()); it != stack.end(); ++it)
e = emit_triangle(*it, e);
stack.clear();
stack.push_back(e);
stack.push_back(ve.index());
}
else if (!left && stack_left)
{
auto e = ve.index();
auto last = stack.back();
while (stack.size() >= 2)
{
*out++ = v;
*out++ = stack[stack.size() - 2];
*out++ = stack.back();
stack.pop_back();
}
stack.pop_back();
for (auto it = stack.begin(); it != std::prev(stack.end()); ++it)
e = emit_triangle(e, *it);
stack.clear();
stack.push_back(last);
stack.push_back(v);
stack.push_back(e);
}
stack_left = left;
if (last)
if (le == re)
break;
}
auto de = le.index();
if (stack_left)
{
for (IndexType i = 0; i + 2 < stack.size(); ++i)
de = emit_triangle(de, stack[i]);
}
else
{
for (IndexType i = stack.size() - 1; i > 1; --i)
de = emit_triangle(de, stack[i]);
}
stack.clear();
}
return out;
}
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<IndexType>(robust_tag, begin, end, std::back_inserter(result));
return result;
}
template <typename IndexType, typename Iterator>
auto monotone_triangulation(Iterator begin, Iterator end)
{