From 7a57c1c1d837858831a9da6883b695436f88cdc4 Mon Sep 17 00:00:00 2001 From: lisyarus Date: Mon, 19 Sep 2022 23:35:33 +0300 Subject: [PATCH] Re-implement monotone chain triangulation to output a DCEL --- examples/triangulation.cpp | 56 +++- .../psemek/cg/triangulation/monotone.hpp | 283 ++++++++++-------- 2 files changed, 202 insertions(+), 137 deletions(-) diff --git a/examples/triangulation.cpp b/examples/triangulation.cpp index 86a7bfb9..1b26eb1b 100644 --- a/examples/triangulation.cpp +++ b/examples/triangulation.cpp @@ -64,14 +64,20 @@ struct triangulation_app { prof::profiler prof("triangulate"); - auto dcel = cg::ear_clipping(geom::fast, points_.begin(), points_.end()); +// auto dcel = cg::ear_clipping(geom::fast, points_.begin(), points_.end()); + auto dcel = cg::monotone_triangulation(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 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> drag_start_; + std::vector closest_points_; gfx::painter painter_; diff --git a/libs/cg/include/psemek/cg/triangulation/monotone.hpp b/libs/cg/include/psemek/cg/triangulation/monotone.hpp index aa52fec3..9d372424 100644 --- a/libs/cg/include/psemek/cg/triangulation/monotone.hpp +++ b/libs/cg/include/psemek/cg/triangulation/monotone.hpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -10,39 +11,43 @@ namespace psemek::cg { - template - OutputIterator monotone_triangulation(RobustTag robust_tag, Iterator begin, Iterator end, OutputIterator out) + template + 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 start(count, 0); - std::vector next(count, 0); - std::vector prev(count, 0); - std::vector twin(count, 0); + using dcel_type = dcel; + 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 used(start.size(), false); - std::vector 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 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 - 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(robust_tag, begin, end, std::back_inserter(result)); return result; } + template auto monotone_triangulation(Iterator begin, Iterator end) {