From a116914fceb56bc6b8c8962db3589d021e31af52 Mon Sep 17 00:00:00 2001 From: Vojtech Bubnik Date: Fri, 29 Jan 2021 16:34:22 +0100 Subject: [PATCH] WIP VoronoiOffset: Squash merge of vb_voronoi_offset Working contour offsetting, skeleton_edges_rough() to detect "important" skeleton edges. Radius of an inscribed circle along the "important" skeleton edges changes slowly, therefore these "important" skeleton edges signify oblong regions possibly needing a gap fill. --- src/libslic3r/Point.hpp | 2 +- src/libslic3r/VoronoiOffset.cpp | 1705 +++++++++++++++++++------- src/libslic3r/VoronoiOffset.hpp | 130 +- src/libslic3r/VoronoiVisualUtils.hpp | 118 +- src/libslic3r/libslic3r.h | 10 + tests/libslic3r/test_voronoi.cpp | 601 ++++++++- 6 files changed, 2018 insertions(+), 548 deletions(-) diff --git a/src/libslic3r/Point.hpp b/src/libslic3r/Point.hpp index a6525af4e..d11af8b58 100644 --- a/src/libslic3r/Point.hpp +++ b/src/libslic3r/Point.hpp @@ -417,7 +417,7 @@ namespace boost { namespace polygon { typedef coord_t coordinate_type; static inline coordinate_type get(const Slic3r::Point& point, orientation_2d orient) { - return (coordinate_type)point((orient == HORIZONTAL) ? 0 : 1); + return static_cast(point((orient == HORIZONTAL) ? 0 : 1)); } }; diff --git a/src/libslic3r/VoronoiOffset.cpp b/src/libslic3r/VoronoiOffset.cpp index f328a8ab2..7ad0f7af4 100644 --- a/src/libslic3r/VoronoiOffset.cpp +++ b/src/libslic3r/VoronoiOffset.cpp @@ -1,18 +1,20 @@ // Polygon offsetting using Voronoi diagram prodiced by boost::polygon. #include "VoronoiOffset.hpp" +#include "libslic3r.h" #include // #define VORONOI_DEBUG_OUT +#include + #ifdef VORONOI_DEBUG_OUT #include #endif namespace Slic3r { - -using VD = Geometry::VoronoiDiagram; +namespace Voronoi { namespace detail { // Intersect a circle with a ray, return the two parameters. @@ -22,9 +24,11 @@ namespace detail { { const Vec2d d = pt - center; #ifndef NDEBUG + // Start point should be inside, end point should be outside the circle. double d0 = (pt - center).norm(); double d1 = (pt + v - center).norm(); - assert(r < std::max(d0, d1) + EPSILON); + assert(d0 < r + SCALED_EPSILON); + assert(d1 > r - SCALED_EPSILON); #endif /* NDEBUG */ const double a = v.squaredNorm(); const double b = 2. * d.dot(v); @@ -194,242 +198,1135 @@ namespace detail { return out; } + // Double vertex equal to a coord_t point after conversion to double. + template + inline bool vertex_equal_to_point(const VertexType &vertex, const Point &ipt) + { + // Convert ipt to doubles, force the 80bit FPU temporary to 64bit and then compare. + // This should work with any settings of math compiler switches and the C++ compiler + // shall understand the memcpies as type punning and it shall optimize them out. +#if 1 + using ulp_cmp_type = boost::polygon::detail::ulp_comparison; + ulp_cmp_type ulp_cmp; + static constexpr int ULPS = boost::polygon::voronoi_diagram_traits::vertex_equality_predicate_type::ULPS; + return ulp_cmp(vertex.x(), double(ipt.x()), ULPS) == ulp_cmp_type::EQUAL && + ulp_cmp(vertex.y(), double(ipt.y()), ULPS) == ulp_cmp_type::EQUAL; +#else + volatile double u = static_cast(ipt.x()); + volatile double v = vertex.x(); + if (u != v) + return false; + u = static_cast(ipt.y()); + v = vertex.y(); + return u == v; +#endif + }; + bool vertex_equal_to_point(const VD::vertex_type *vertex, const Point &ipt) + { return vertex_equal_to_point(*vertex, ipt); } + + double dist_to_site(const Lines &lines, const VD::cell_type &cell, const Vec2d &point) + { + const Line &line = lines[cell.source_index()]; + return cell.contains_point() ? + (((cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line.a : line.b).cast() - point).norm() : + (Geometry::foot_pt(line.a.cast(), (line.b - line.a).cast(), point) - point).norm(); + }; + + bool on_site(const Lines &lines, const VD::cell_type &cell, const Vec2d &pt) + { + const Line &line = lines[cell.source_index()]; + auto on_contour = [&pt](const Point &ipt) { return detail::vertex_equal_to_point(pt, ipt); }; + if (cell.contains_point()) { + return on_contour(contour_point(cell, line)); + } else { + assert(! (on_contour(line.a) && on_contour(line.b))); + return on_contour(line.a) || on_contour(line.b); + } + }; + + // For a Voronoi segment starting with voronoi_point1 and ending with voronoi_point2, + // defined by a bisector of Voronoi sites pt1_site and pt2_site (line) + // find two points on the Voronoi bisector, that delimit regions with dr/dl measure + // lower / higher than threshold_dr_dl. + // + // Linear segment from voronoi_point1 to return.first and + // linear segment from return.second to voronoi_point2 have dr/dl measure + // higher than threshold_dr_dl. + // If such respective segment does not exist, then return.first resp. return.second is nan. + std::pair point_point_dr_dl_thresholds( + // Two Voronoi sites + const Point &pt1_site, const Point &pt2_site, + // End points of a Voronoi segment + const Vec2d &voronoi_point1, const Vec2d &voronoi_point2, + // Threshold of the skeleton function, where alpha is an angle of a sharp convex corner with the same dr/dl. + const double threshold_tan_alpha_half) + { + // sympy code to calculate +-x + // of a linear bisector of pt1_site, pt2_site parametrized with pt + x * v, |v| = 1 + // where dr/dl = threshold_dr_dl + // equals d|pt1_site - pt + x * v| / dx = threshold_dr_dl + // + // y = sqrt(x^2 + d^2) + // dy = diff(y, x) + // solve(dy - c, x) + // + // Project voronoi_point1/2 to line_site. + Vec2d dir_y = (pt2_site - pt1_site).cast(); + Vec2d dir_x = Vec2d(- dir_y.y(), dir_y.x()).normalized(); + Vec2d cntr = 0.5 * (pt1_site.cast() + pt2_site.cast()); + double t1 = (voronoi_point1 - cntr).dot(dir_x); + double t2 = (voronoi_point2 - cntr).dot(dir_x); + if (t1 > t2) { + t1 = -t1; + t2 = -t2; + dir_x = - dir_x; + } + auto x = 0.5 * dir_y.norm() * threshold_tan_alpha_half; + static constexpr double nan = std::numeric_limits::quiet_NaN(); + auto out = std::make_pair(Vec2d(nan, nan), Vec2d(nan, nan)); + if (t2 > -x && t1 < x) { + // Intervals overlap. + dir_x *= x; + out.first = (t1 < -x) ? cntr - dir_x : voronoi_point1; + out.second = (t2 > +x) ? cntr + dir_x : voronoi_point2; + } + return out; + } + + // For a Voronoi segment starting with voronoi_point1 and ending with voronoi_point2, + // defined by a bisector of Voronoi sites pt_site and line site (parabolic arc) + // find two points on the Voronoi parabolic arc, that delimit regions with dr/dl measure + // lower / higher than threshold_dr_dl. + // + // Parabolic arc from voronoi_point1 to return.first and + // parabolic arc from return.second to voronoi_point2 have dr/dl measure + // higher than threshold_dr_dl. + // If such respective segment does not exist, then return.first resp. return.second is nan. + std::pair point_segment_dr_dl_thresholds( + // Two Voronoi sites + const Point &pt_site, const Line &line_site, + // End points of a Voronoi segment + const Vec2d &voronoi_point1, const Vec2d &voronoi_point2, + // Threshold of the skeleton function, where alpha is an angle of a sharp convex corner with the same dr/dl. + const double threshold_tan_alpha_half) + { + // sympy code to calculate +-x + // of a parabola y = ax^2 + b + // where dr/dl = threshold_dr_dl + // + // a = 1 / (4 * b) + // y = a*x**2 + b + // dy = diff(y, x) + // solve(dy / sqrt(1 + dy**2) - c, x) + // + // Foot point of the point site on the line site. + Vec2d ft = Geometry::foot_pt(line_site, pt_site); + // Minimum distance of the bisector (parabolic arc) from the two sites, squared. + Vec2d dir_pt_ft = pt_site.cast() - ft; + double b = 0.5 * dir_pt_ft.norm(); + static constexpr double nan = std::numeric_limits::quiet_NaN(); + auto out = std::make_pair(Vec2d(nan, nan), Vec2d(nan, nan)); + { + // +x, -x are the two parameters along the line_site, where threshold_tan_alpha_half is met. + double x = 2. * b * threshold_tan_alpha_half; + // Project voronoi_point1/2 to line_site. + Vec2d dir_x = (line_site.b - line_site.a).cast().normalized(); + double t1 = (voronoi_point1 - ft).dot(dir_x); + double t2 = (voronoi_point2 - ft).dot(dir_x); + if (t1 > t2) { + t1 = -t1; + t2 = -t2; + dir_x = - dir_x; + } + if (t2 > -x && t1 < x) { + // Intervals overlap. + bool t1_valid = t1 < -x; + bool t2_valid = t2 > +x; + // Direction of the Y axis of the parabola. + Vec2d dir_y(- dir_x.y(), dir_x.x()); + // Orient the Y axis towards the point site. + if (dir_y.dot(dir_pt_ft) < 0.) + dir_y = - dir_y; + // Equation of the parabola: y = b + a * x^2 + double a = 0.25 / b; + dir_x *= x; + dir_y *= b + a * x * x; + out.first = t1_valid ? ft - dir_x + dir_y : voronoi_point1; + out.second = t2_valid ? ft + dir_x + dir_y : voronoi_point2; + } + } + return out; + } + + std::pair point_point_skeleton_thresholds( + // Two Voronoi sites + const Point &pt1_site, const Point &pt2_site, + // End points of a Voronoi segment + const Vec2d &voronoi_point1, const Vec2d &voronoi_point2, + // Threshold of the skeleton function. + const double tan_alpha_half) + { + // Project voronoi_point1/2 to line_site. + Vec2d dir_y = (pt2_site - pt1_site).cast(); + Vec2d dir_x = Vec2d(- dir_y.y(), dir_y.x()).normalized(); + Vec2d cntr = 0.5 * (pt1_site.cast() + pt2_site.cast()); + double t1 = (voronoi_point1 - cntr).dot(dir_x); + double t2 = (voronoi_point2 - cntr).dot(dir_x); + if (t1 > t2) { + t1 = -t1; + t2 = -t2; + dir_x = - dir_x; + } + auto x = 0.5 * dir_y.norm() * tan_alpha_half; + static constexpr double nan = std::numeric_limits::quiet_NaN(); + auto out = std::make_pair(Vec2d(nan, nan), Vec2d(nan, nan)); + if (t2 > -x && t1 < x) { + // Intervals overlap. + dir_x *= x; + out.first = (t1 < -x) ? cntr - dir_x : voronoi_point1; + out.second = (t2 > +x) ? cntr + dir_x : voronoi_point2; + } + return out; + } + + std::pair point_segment_skeleton_thresholds( + // Two Voronoi sites + const Point &pt_site, const Line &line_site, + // End points of a Voronoi segment + const Vec2d &voronoi_point1, const Vec2d &voronoi_point2, + // Threshold of the skeleton function. + const double threshold_cos_alpha) + { + // Foot point of the point site on the line site. + Vec2d ft = Geometry::foot_pt(line_site, pt_site); + // Minimum distance of the bisector (parabolic arc) from the two sites, squared. + Vec2d dir_pt_ft = pt_site.cast() - ft; + // Distance of Voronoi point site from the Voronoi line site. + double l = dir_pt_ft.norm(); + static constexpr double nan = std::numeric_limits::quiet_NaN(); + auto out = std::make_pair(Vec2d(nan, nan), Vec2d(nan, nan)); + // +x, -x are the two parameters along the line_site, where threshold is met. + double r = l / (1. + threshold_cos_alpha); + double x2 = r * r - Slic3r::sqr(l - r); + double x = sqrt(x2); + // Project voronoi_point1/2 to line_site. + Vec2d dir_x = (line_site.b - line_site.a).cast().normalized(); + double t1 = (voronoi_point1 - ft).dot(dir_x); + double t2 = (voronoi_point2 - ft).dot(dir_x); + if (t1 > t2) { + t1 = -t1; + t2 = -t2; + dir_x = - dir_x; + } + if (t2 > -x && t1 < x) { + // Intervals overlap. + bool t1_valid = t1 < -x; + bool t2_valid = t2 > +x; + // Direction of the Y axis of the parabola. + Vec2d dir_y(- dir_x.y(), dir_x.x()); + // Orient the Y axis towards the point site. + if (dir_y.dot(dir_pt_ft) < 0.) + dir_y = - dir_y; + // Equation of the parabola: y = b + a * x^2 + double a = 0.5 / l; + dir_x *= x; + dir_y *= 0.5 * l + a * x2; + out.first = t1_valid ? ft - dir_x + dir_y : voronoi_point1; + out.second = t2_valid ? ft + dir_x + dir_y : voronoi_point2; + } + return out; + } + } // namespace detail -Polygons voronoi_offset( - const Geometry::VoronoiDiagram &vd, - const Lines &lines, - double offset_distance, - double discretization_error) -{ #ifndef NDEBUG +namespace debug +{ // Verify that twin halfedges are stored next to the other in vd. - for (size_t i = 0; i < vd.num_edges(); i += 2) { - const VD::edge_type &e = vd.edges()[i]; - const VD::edge_type &e2 = vd.edges()[i + 1]; - assert(e.twin() == &e2); - assert(e2.twin() == &e); - assert(e.is_secondary() == e2.is_secondary()); - if (e.is_secondary()) { - assert(e.cell()->contains_point() != e2.cell()->contains_point()); - const VD::edge_type &ex = (e.cell()->contains_point() ? e : e2); - // Verify that the Point defining the cell left of ex is an end point of a segment - // defining the cell right of ex. - const Line &line0 = lines[ex.cell()->source_index()]; - const Line &line1 = lines[ex.twin()->cell()->source_index()]; - const Point &pt = (ex.cell()->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b; - assert(pt == line1.a || pt == line1.b); + bool verify_twin_halfedges_successive(const VD &vd, const Lines &lines) + { + for (size_t i = 0; i < vd.num_edges(); i += 2) { + const VD::edge_type &e = vd.edges()[i]; + const VD::edge_type &e2 = vd.edges()[i + 1]; + assert(e.twin() == &e2); + assert(e2.twin() == &e); + assert(e.is_secondary() == e2.is_secondary()); + if (e.is_secondary()) { + assert(e.cell()->contains_point() != e2.cell()->contains_point()); + const VD::edge_type &ex = (e.cell()->contains_point() ? e : e2); + // Verify that the Point defining the cell left of ex is an end point of a segment + // defining the cell right of ex. + const Line &line0 = lines[ex.cell()->source_index()]; + const Line &line1 = lines[ex.twin()->cell()->source_index()]; + const Point &pt = (ex.cell()->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b; + assert(pt == line1.a || pt == line1.b); + } } + return true; } + + bool verify_inside_outside_annotations(const VD &vd) + { + // Verify that "Colors" are set at all Voronoi entities. + for (const VD::vertex_type &v : vd.vertices()) { + assert(! v.is_degenerate()); + assert(vertex_category(v) != VertexCategory::Unknown); + } + for (const VD::edge_type &e : vd.edges()) + assert(edge_category(e) != EdgeCategory::Unknown); + for (const VD::cell_type &c : vd.cells()) { + // Unfortunately denegerate cells could be created, which reference a null edge. + // https://github.com/boostorg/polygon/issues/47 + assert(c.is_degenerate() || cell_category(c) != CellCategory::Unknown); + } + + // Verify consistency between markings of Voronoi cells, edges and verticies. + for (const VD::cell_type &cell : vd.cells()) { + if (cell.is_degenerate()) { + // Unfortunately denegerate cells could be created, which reference a null edge. + // https://github.com/boostorg/polygon/issues/47 + continue; + } + const VD::edge_type *first_edge = cell.incident_edge(); + const VD::edge_type *edge = first_edge; + CellCategory cc = cell_category(cell); + size_t num_vertices_on_contour = 0; + size_t num_vertices_inside = 0; + size_t num_vertices_outside = 0; + size_t num_edges_point_to_contour = 0; + size_t num_edges_point_inside = 0; + size_t num_edges_point_outside = 0; + do { + { + EdgeCategory ec = edge_category(edge); + switch (ec) { + case EdgeCategory::PointsInside: + assert(edge->vertex0() != nullptr && edge->vertex1() != nullptr); + ++ num_edges_point_inside; break; + case EdgeCategory::PointsOutside: +// assert(edge->vertex0() != nullptr); + ++ num_edges_point_outside; break; + case EdgeCategory::PointsToContour: + assert(edge->vertex1() != nullptr); + ++ num_edges_point_to_contour; break; + default: + assert(false); + } + } + { + VertexCategory vc = (edge->vertex1() == nullptr) ? VertexCategory::Outside : vertex_category(edge->vertex1()); + switch (vc) { + case VertexCategory::Inside: + ++ num_vertices_inside; break; + case VertexCategory::Outside: + ++ num_vertices_outside; break; + case VertexCategory::OnContour: + ++ num_vertices_on_contour; break; + default: + assert(false); + } + } + { + const VD::cell_type *cell_other = edge->twin()->cell(); + const CellCategory cc_other = cell_category(cell_other); + assert(cc_other != CellCategory::Unknown); + switch (cc) { + case CellCategory::Boundary: + assert(cc_other != CellCategory::Boundary || cell_other->contains_segment()); + break; + case CellCategory::Inside: + assert(cc_other == CellCategory::Inside || cc_other ==CellCategory::Boundary); + break; + case CellCategory::Outside: + assert(cc_other == CellCategory::Outside || cc_other == CellCategory::Boundary); + break; + default: + assert(false); + break; + } + } + edge = edge->next(); + } while (edge != first_edge); + + switch (cc) { + case CellCategory::Boundary: + assert(cell.contains_segment()); + assert(num_edges_point_to_contour == 2); + assert(num_vertices_on_contour == 2); + assert(num_vertices_inside > 0); + assert(num_vertices_outside > 0); + assert(num_edges_point_inside > 0); + assert(num_edges_point_outside > 0); + break; + case CellCategory::Inside: + assert(num_vertices_on_contour <= 1); + assert(num_edges_point_to_contour <= 1); + assert(num_vertices_inside > 0); + assert(num_vertices_outside == 0); + assert(num_edges_point_inside > 0); + assert(num_edges_point_outside == 0); + break; + case CellCategory::Outside: + assert(num_vertices_on_contour <= 1); + assert(num_edges_point_to_contour <= 1); + assert(num_vertices_inside == 0); + assert(num_vertices_outside > 0); + assert(num_edges_point_inside == 0); + assert(num_edges_point_outside > 0); + break; + default: + assert(false); + break; + } + } + + return true; + } + + bool verify_vertices_on_contour(const VD &vd, const Lines &lines) + { + for (const VD::edge_type &edge : vd.edges()) { + const VD::vertex_type *v = edge.vertex0(); + if (v != nullptr) { + bool on_contour = vertex_category(v) == VertexCategory::OnContour; + assert(detail::on_site(lines, *edge.cell(), vertex_point(v)) == on_contour); + assert(detail::on_site(lines, *edge.twin()->cell(), vertex_point(v)) == on_contour); + } + } + return true; + } + + bool verify_signed_distances(const VD &vd, const Lines &lines, const std::vector &signed_distances) + { + for (const VD::edge_type &edge : vd.edges()) { + const VD::vertex_type *v = edge.vertex0(); + double d = (v == nullptr) ? std::numeric_limits::max() : signed_distances[v - &vd.vertices().front()]; + if (v == nullptr || vertex_category(v) == VertexCategory::Outside) + assert(d > 0.); + else if (vertex_category(v) == VertexCategory::OnContour) + assert(d == 0.); + else + assert(d < 0.); + if (v != nullptr) { + double err = std::abs(detail::dist_to_site(lines, *edge.cell(), vertex_point(v)) - std::abs(d)); + double err2 = std::abs(detail::dist_to_site(lines, *edge.twin()->cell(), vertex_point(v)) - std::abs(d)); + assert(err < SCALED_EPSILON); + assert(err2 < SCALED_EPSILON); + } + } + return true; + } + + bool verify_offset_intersection_points(const VD &vd, const Lines &lines, const double offset_distance, const std::vector &offset_intersection_points) + { + const VD::edge_type *front_edge = &vd.edges().front(); + const double d = std::abs(offset_distance); + for (const VD::edge_type &edge : vd.edges()) { + const Vec2d &p = offset_intersection_points[&edge - front_edge]; + if (edge_offset_has_intersection(p)) { + double err = std::abs(detail::dist_to_site(lines, *edge.cell(), p) - d); + double err2 = std::abs(detail::dist_to_site(lines, *edge.twin()->cell(), p) - d); + assert(err < SCALED_EPSILON); + assert(err2 < SCALED_EPSILON); + } + } + return true; + } + +} #endif // NDEBUG - enum class EdgeState : unsigned char { - // Initial state, don't know. - Unknown, - // This edge will certainly not be intersected by the offset curve. - Inactive, - // This edge will certainly be intersected by the offset curve. - Active, - // This edge will possibly be intersected by the offset curve. - Possible +void reset_inside_outside_annotations(VD &vd) +{ + for (const VD::vertex_type &v : vd.vertices()) + set_vertex_category(const_cast(v), VertexCategory::Unknown); + for (const VD::edge_type &e : vd.edges()) + set_edge_category(const_cast(e), EdgeCategory::Unknown); + for (const VD::cell_type &c : vd.cells()) + set_cell_category(const_cast(c), CellCategory::Unknown); +} + +void annotate_inside_outside(VD &vd, const Lines &lines) +{ + assert(debug::verify_twin_halfedges_successive(vd, lines)); + + reset_inside_outside_annotations(vd); + +#ifdef VORONOI_DEBUG_OUT + BoundingBox bbox; + { + bbox.merge(get_extents(lines)); + bbox.min -= (0.01 * bbox.size().cast()).cast(); + bbox.max += (0.01 * bbox.size().cast()).cast(); + } + static int irun = 0; + ++ irun; + dump_voronoi_to_svg(debug_out_path("voronoi-offset-initial-%d.svg", irun).c_str(), vd, Points(), lines); +#endif // VORONOI_DEBUG_OUT + + // Set a VertexCategory, verify validity of the operation. + auto annotate_vertex = [](const VD::vertex_type *vertex, VertexCategory new_vertex_category) { + VertexCategory vc = vertex_category(vertex); + assert(vc == VertexCategory::Unknown || vc == new_vertex_category); + assert(new_vertex_category == VertexCategory::Inside || + new_vertex_category == VertexCategory::Outside || + new_vertex_category == VertexCategory::OnContour); + set_vertex_category(const_cast(vertex), new_vertex_category); }; - enum class CellState : unsigned char { - // Initial state, don't know. - Unknown, - // Inactive cell is inside for outside curves and outside for inside curves. - Inactive, - // Active cell is outside for outside curves and inside for inside curves. - Active, - // Boundary cell is intersected by the input segment, part of it is active. - Boundary + // Set an EdgeCategory, verify validity of the operation. + auto annotate_edge = [](const VD::edge_type *edge, EdgeCategory new_edge_category) { + EdgeCategory ec = edge_category(edge); + assert(ec == EdgeCategory::Unknown || ec == new_edge_category); +#ifndef NDEBUG + switch (new_edge_category) { + case EdgeCategory::PointsInside: + assert(edge->vertex0() != nullptr); + assert(edge->vertex1() != nullptr); + break; + case EdgeCategory::PointsOutside: + // assert(edge->vertex0() != nullptr); + break; + case EdgeCategory::PointsToContour: + assert(edge->vertex1() != nullptr); + break; + default: + assert(false); + } + +#endif // NDEBUG + set_edge_category(const_cast(edge), new_edge_category); }; - // Mark edges with outward vertex pointing outside the polygons, thus there is a chance - // that such an edge will have an intersection with our desired offset curve. - bool outside = offset_distance > 0.; - std::vector edge_state(vd.num_edges(), EdgeState::Unknown); - std::vector cell_state(vd.num_cells(), CellState::Unknown); - const VD::edge_type *front_edge = &vd.edges().front(); - const VD::cell_type *front_cell = &vd.cells().front(); - auto set_edge_state_initial = [&edge_state, front_edge](const VD::edge_type *edge, EdgeState new_edge_type) { - EdgeState &edge_type = edge_state[edge - front_edge]; - assert(edge_type == EdgeState::Unknown || edge_type == new_edge_type); - assert(new_edge_type == EdgeState::Possible || new_edge_type == EdgeState::Inactive); - edge_type = new_edge_type; - }; - auto set_edge_state_final = [&edge_state, front_edge](const size_t edge_id, EdgeState new_edge_type) { - EdgeState &edge_type = edge_state[edge_id]; - assert(edge_type == EdgeState::Possible || edge_type == new_edge_type); - assert(new_edge_type == EdgeState::Active || new_edge_type == EdgeState::Inactive); - edge_type = new_edge_type; - }; - auto set_cell_state = [&cell_state, front_cell](const VD::cell_type *cell, CellState new_cell_type) -> bool { - CellState &cell_type = cell_state[cell - front_cell]; - assert(cell_type == CellState::Active || cell_type == CellState::Inactive || cell_type == CellState::Boundary || cell_type == CellState::Unknown); - assert(new_cell_type == CellState::Active || new_cell_type == CellState::Inactive || new_cell_type == CellState::Boundary); - switch (cell_type) { - case CellState::Unknown: + // Set a CellCategory, verify validity of the operation. + // Handle marking of boundary cells (first time the cell is marked as outside, the other time as inside). + // Returns true if the current cell category was modified. + auto annotate_cell = [](const VD::cell_type *cell, CellCategory new_cell_category) -> bool { + CellCategory cc = cell_category(cell); + assert(cc == CellCategory::Inside || cc == CellCategory::Outside || cc == CellCategory::Boundary || cc == CellCategory::Unknown); + assert(new_cell_category == CellCategory::Inside || new_cell_category == CellCategory::Outside || new_cell_category == CellCategory::Boundary); + switch (cc) { + case CellCategory::Unknown: + // Old category unknown, just write the new category. break; - case CellState::Active: - if (new_cell_type == CellState::Inactive) - new_cell_type = CellState::Boundary; + case CellCategory::Outside: + if (new_cell_category == CellCategory::Inside) + new_cell_category = CellCategory::Boundary; break; - case CellState::Inactive: - if (new_cell_type == CellState::Active) - new_cell_type = CellState::Boundary; + case CellCategory::Inside: + if (new_cell_category == CellCategory::Outside) + new_cell_category = CellCategory::Boundary; break; - case CellState::Boundary: + case CellCategory::Boundary: return false; } - if (cell_type != new_cell_type) { - cell_type = new_cell_type; + if (cc != new_cell_category) { + set_cell_category(const_cast(cell), new_cell_category); return true; } return false; }; + // The next loop is supposed to annotate the "on input contour" vertices, but due to + // merging very close Voronoi vertices together by boost::polygon Voronoi generator + // the condition may not always be met. It should be safe to mark all Voronoi very close + // to the input contour as on contour. + for (const VD::edge_type &edge : vd.edges()) { + const VD::vertex_type *v = edge.vertex0(); + if (v != nullptr) { + bool on_contour = detail::on_site(lines, *edge.cell(), vertex_point(v)); +#ifndef NDEBUG + bool on_contour2 = detail::on_site(lines, *edge.twin()->cell(), vertex_point(v)); + assert(on_contour == on_contour2); +#endif // NDEBUG + if (on_contour) + annotate_vertex(v, VertexCategory::OnContour); + } + } + + // One side of a secondary edge is always on the source contour. Mark these vertices as OnContour. + // See the comment at the loop before, the condition may not always be met. + for (const VD::edge_type &edge : vd.edges()) { + if (edge.is_secondary() && edge.vertex0() != nullptr) { + assert(edge.is_linear()); + assert(edge.cell()->contains_point() != edge.twin()->cell()->contains_point()); + // The point associated with the point site shall be equal with one vertex of this Voronoi edge. + const Point &pt_on_contour = edge.cell()->contains_point() ? contour_point(*edge.cell(), lines) : contour_point(*edge.twin()->cell(), lines); + auto on_contour = [&pt_on_contour](const VD::vertex_type *v) { return detail::vertex_equal_to_point(v, pt_on_contour); }; + if (edge.vertex1() == nullptr) { + assert(on_contour(edge.vertex0())); + annotate_vertex(edge.vertex0(), VertexCategory::OnContour); + } else { + // Only one of the two vertices may lie on input contour. + const VD::vertex_type *v0 = edge.vertex0(); + const VD::vertex_type *v1 = edge.vertex1(); + VertexCategory v0_category = vertex_category(v0); + VertexCategory v1_category = vertex_category(v1); + assert(v0_category != VertexCategory::OnContour || v1_category != VertexCategory::OnContour); + assert(! (on_contour(v0) && on_contour(v1))); + if (on_contour(v0)) + annotate_vertex(v0, VertexCategory::OnContour); + else { + assert(on_contour(v1)); + annotate_vertex(v1, VertexCategory::OnContour); + } + } + } + } + + assert(debug::verify_vertices_on_contour(vd, lines)); + for (const VD::edge_type &edge : vd.edges()) if (edge.vertex1() == nullptr) { // Infinite Voronoi edge separating two Point sites or a Point site and a Segment site. - // Infinite edge is always outside and it has at least one valid vertex. + // Infinite edge is always outside and it references at least one valid vertex. + assert(edge.is_infinite()); + assert(edge.is_linear()); assert(edge.vertex0() != nullptr); - set_edge_state_initial(&edge, outside ? EdgeState::Possible : EdgeState::Inactive); + const VD::cell_type *cell = edge.cell(); + const VD::cell_type *cell2 = edge.twin()->cell(); + // A Point-Segment secondary Voronoi edge touches the input contour, a Point-Point Voronoi + // edge does not. + assert(edge.is_secondary() ? (cell->contains_segment() != cell2->contains_segment()) : + (cell->contains_point() == cell2->contains_point())); + annotate_edge(&edge, EdgeCategory::PointsOutside); // Opposite edge of an infinite edge is certainly not active. - set_edge_state_initial(edge.twin(), EdgeState::Inactive); - if (edge.is_secondary()) { - // edge.vertex0() must lie on source contour. - const VD::cell_type *cell = edge.cell(); - const VD::cell_type *cell2 = edge.twin()->cell(); - if (cell->contains_segment()) - std::swap(cell, cell2); - // State of a cell containing a boundary point is known. - assert(cell->contains_point()); - set_cell_state(cell, outside ? CellState::Active : CellState::Inactive); - // State of a cell containing a boundary edge is Boundary. - assert(cell2->contains_segment()); - set_cell_state(cell2, CellState::Boundary); - } + annotate_edge(edge.twin(), edge.is_secondary() ? EdgeCategory::PointsToContour : EdgeCategory::PointsOutside); + annotate_vertex(edge.vertex0(), edge.is_secondary() ? VertexCategory::OnContour : VertexCategory::Outside); + // edge.vertex1() is null, it is implicitely outside. + if (cell->contains_segment()) + std::swap(cell, cell2); + // State of a cell containing a boundary point is certainly outside. + assert(cell->contains_point()); + annotate_cell(cell, CellCategory::Outside); + assert(edge.is_secondary() == cell2->contains_segment()); + annotate_cell(cell2, cell2->contains_point() ? CellCategory::Outside : CellCategory::Boundary); } else if (edge.vertex0() != nullptr) { - // Finite edge. + assert(edge.is_finite()); const VD::cell_type *cell = edge.cell(); const Line *line = cell->contains_segment() ? &lines[cell->source_index()] : nullptr; if (line == nullptr) { cell = edge.twin()->cell(); line = cell->contains_segment() ? &lines[cell->source_index()] : nullptr; } + // Only one of the two vertices may lie on input contour. + assert(! edge.is_linear() || vertex_category(edge.vertex0()) != VertexCategory::OnContour || vertex_category(edge.vertex1()) != VertexCategory::OnContour); + // Now classify the Voronoi vertices and edges as inside outside, if at least one Voronoi + // site is a Segment site. + // Inside / outside classification of Point - Point Voronoi edges will be done later + // by a propagation (seed fill). if (line) { const VD::vertex_type *v1 = edge.vertex1(); const VD::cell_type *cell2 = (cell == edge.cell()) ? edge.twin()->cell() : edge.cell(); - assert(v1); - const Point *pt_on_contour = nullptr; - if (cell == edge.cell() && edge.twin()->cell()->contains_segment()) { - // Constrained bisector of two segments. + assert(v1 != nullptr); + VertexCategory v0_category = vertex_category(edge.vertex0()); + VertexCategory v1_category = vertex_category(edge.vertex1()); + bool on_contour = v0_category == VertexCategory::OnContour || v1_category == VertexCategory::OnContour; +#ifndef NDEBUG + if (! on_contour && cell == edge.cell() && edge.twin()->cell()->contains_segment()) { + // Constrained bisector of two segments. Vojtech is not quite sure whether the Voronoi generator is robust enough + // to always connect at least one secondary edge to an input contour point. Catch it here. + assert(edge.is_linear()); + // OnContour state of this edge is not known yet. + const Point *pt_on_contour = nullptr; // If the two segments share a point, then one end of the current Voronoi edge shares this point as well. - // Find pt_on_contour if it exists. + // A bisector may not necessarily connect to the source contour. Find pt_on_contour if it exists. const Line &line2 = lines[cell2->source_index()]; if (line->a == line2.b) pt_on_contour = &line->a; else if (line->b == line2.a) pt_on_contour = &line->b; - } else if (edge.is_secondary()) { - assert(edge.is_linear()); - // One end of the current Voronoi edge shares a point of a contour. - assert(edge.cell()->contains_point() != edge.twin()->cell()->contains_point()); - const Line &line2 = lines[cell2->source_index()]; - pt_on_contour = &((cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line2.a : line2.b); - } - if (pt_on_contour) { - // One end of the current Voronoi edge shares a point of a contour. - // Find out which one it is. - const VD::vertex_type *v0 = edge.vertex0(); - Vec2d vec0(v0->x() - pt_on_contour->x(), v0->y() - pt_on_contour->y()); - Vec2d vec1(v1->x() - pt_on_contour->x(), v1->y() - pt_on_contour->y()); - double d0 = vec0.squaredNorm(); - double d1 = vec1.squaredNorm(); - assert(std::min(d0, d1) < SCALED_EPSILON * SCALED_EPSILON); - if (d0 < d1) { - // v0 is equal to pt. - } else { - // Skip secondary edge pointing to a contour point. - set_edge_state_initial(&edge, EdgeState::Inactive); - continue; + if (pt_on_contour) { + const VD::vertex_type *v0 = edge.vertex0(); + auto on_contour = [&pt_on_contour](const VD::vertex_type *v) { + return std::abs(v->x() - pt_on_contour->x()) < 0.5001 && + std::abs(v->y() - pt_on_contour->y()) < 0.5001; + }; + assert(! on_contour(v0) && ! on_contour(v1)); } } - Vec2d l0(line->a.cast()); - Vec2d lv((line->b - line->a).cast()); - double side = cross2(lv, Vec2d(v1->x(), v1->y()) - l0); - bool edge_active = outside ? (side < 0.) : (side > 0.); - set_edge_state_initial(&edge, edge_active ? EdgeState::Possible : EdgeState::Inactive); - assert(cell->contains_segment()); - set_cell_state(cell, - pt_on_contour ? CellState::Boundary : - edge_active ? CellState::Active : CellState::Inactive); - set_cell_state(cell2, - (pt_on_contour && cell2->contains_segment()) ? - CellState::Boundary : - edge_active ? CellState::Active : CellState::Inactive); +#endif // NDEBUG + if (on_contour && v1_category == VertexCategory::OnContour) { + // Skip secondary edge pointing to a contour point. + annotate_edge(&edge, EdgeCategory::PointsToContour); + } else { + // v0 is certainly not on the input polygons. + // Is v1 inside or outside the input polygons? + // The Voronoi vertex coordinate is in doubles, calculate orientation in doubles. + Vec2d l0(line->a.cast()); + Vec2d lv((line->b - line->a).cast()); + double side = cross2(Vec2d(v1->x(), v1->y()) - l0, lv); + // No Voronoi edge could connect two vertices of input polygons. + assert(side != 0.); + auto vc = side > 0. ? VertexCategory::Outside : VertexCategory::Inside; + annotate_vertex(v1, vc); + auto ec = vc == VertexCategory::Outside ? EdgeCategory::PointsOutside : EdgeCategory::PointsInside; + annotate_edge(&edge, ec); + // Annotate the twin edge and its vertex. As the Voronoi edge may never cross the input + // contour, the twin edge and its vertex will share the property of edge. + annotate_vertex(edge.vertex0(), on_contour ? VertexCategory::OnContour : vc); + annotate_edge(edge.twin(), on_contour ? EdgeCategory::PointsToContour : ec); + assert(cell->contains_segment()); + annotate_cell(cell, on_contour ? CellCategory::Boundary : + (vc == VertexCategory::Outside ? CellCategory::Outside : CellCategory::Inside)); + annotate_cell(cell2, (on_contour && cell2->contains_segment()) ? CellCategory::Boundary : + (vc == VertexCategory::Outside ? CellCategory::Outside : CellCategory::Inside)); + } } } - { - // Perform one round of expansion marking Voronoi edges and cells next to boundary cells as active / inactive. - std::vector cell_queue; - for (const VD::edge_type &edge : vd.edges()) - if (edge_state[&edge - front_edge] == EdgeState::Unknown) { - assert(edge.cell()->contains_point() && edge.twin()->cell()->contains_point()); - // Edge separating two point sources, not yet classified as inside / outside. - CellState cs = cell_state[edge.cell() - front_cell]; - CellState cs2 = cell_state[edge.twin()->cell() - front_cell]; - if (cs != CellState::Unknown || cs2 != CellState::Unknown) { - if (cs == CellState::Unknown) { - cs = cs2; - if (set_cell_state(edge.cell(), cs)) - cell_queue.emplace_back(edge.cell()); - } else if (set_cell_state(edge.twin()->cell(), cs)) - cell_queue.emplace_back(edge.twin()->cell()); - EdgeState es = (cs == CellState::Active) ? EdgeState::Possible : EdgeState::Inactive; - set_edge_state_initial(&edge, es); - set_edge_state_initial(edge.twin(), es); - } else { - const VD::edge_type *e = edge.twin()->rot_prev(); - do { - EdgeState es = edge_state[e->twin() - front_edge]; - if (es != EdgeState::Unknown) { - assert(es == EdgeState::Possible || es == EdgeState::Inactive); - set_edge_state_initial(&edge, es); - CellState cs = (es == EdgeState::Possible) ? CellState::Active : CellState::Inactive; - if (set_cell_state(edge.cell(), cs)) - cell_queue.emplace_back(edge.cell()); - if (set_cell_state(edge.twin()->cell(), cs)) - cell_queue.emplace_back(edge.twin()->cell()); - break; - } - e = e->rot_prev(); - } while (e != edge.twin()); + + assert(debug::verify_vertices_on_contour(vd, lines)); + + // Now most Voronoi vertices, edges and cells are annotated, with the exception of some + // edges separating two Point sites, their cells and vertices. + // Perform one round of expansion marking Voronoi edges and cells next to boundary cells. + std::vector cell_queue; + for (const VD::edge_type &edge : vd.edges()) { + assert((edge_category(edge) == EdgeCategory::Unknown) == (edge_category(edge.twin()) == EdgeCategory::Unknown)); + if (edge_category(edge) == EdgeCategory::Unknown) { + assert(edge.is_finite()); + const VD::cell_type &cell = *edge.cell(); + const VD::cell_type &cell2 = *edge.twin()->cell(); + assert(cell.contains_point() && cell2.contains_point()); + CellCategory cc = cell_category(cell); + CellCategory cc2 = cell_category(cell2); + assert(cc != CellCategory::Boundary && cc2 != CellCategory::Boundary); + CellCategory cc_new = cc; + if (cc_new == CellCategory::Unknown) + cc_new = cc2; + else + assert(cc2 == CellCategory::Unknown || cc == cc2); + if (cc_new == CellCategory::Unknown) { + VertexCategory vc = vertex_category(edge.vertex0()); + assert(vc != VertexCategory::OnContour); + if (vc != VertexCategory::Unknown) + cc_new = (vc == VertexCategory::Outside) ? CellCategory::Outside : CellCategory::Inside; + } + if (cc_new != CellCategory::Unknown) { + VertexCategory vc = (cc_new == CellCategory::Outside) ? VertexCategory::Outside : VertexCategory::Inside; + annotate_vertex(edge.vertex0(), vc); + annotate_vertex(edge.vertex1(), vc); + auto ec_new = (cc_new == CellCategory::Outside) ? EdgeCategory::PointsOutside : EdgeCategory::PointsInside; + annotate_edge(&edge, ec_new); + annotate_edge(edge.twin(), ec_new); + if (cc != cc_new) { + annotate_cell(&cell, cc_new); + cell_queue.emplace_back(&cell); + } + if (cc2 != cc_new) { + annotate_cell(&cell2, cc_new); + cell_queue.emplace_back(&cell2); } } - // Do a final seed fill over Voronoi cells and unmarked Voronoi edges. - while (! cell_queue.empty()) { - const VD::cell_type *cell = cell_queue.back(); - const CellState cs = cell_state[cell - front_cell]; - cell_queue.pop_back(); - const VD::edge_type *first_edge = cell->incident_edge(); - const VD::edge_type *edge = cell->incident_edge(); - EdgeState es = (cs == CellState::Active) ? EdgeState::Possible : EdgeState::Inactive; - do { - if (set_cell_state(edge->twin()->cell(), cs)) { - set_edge_state_initial(edge, es); - set_edge_state_initial(edge->twin(), es); - cell_queue.emplace_back(edge->twin()->cell()); - } - edge = edge->next(); - } while (edge != first_edge); } } + assert(debug::verify_vertices_on_contour(vd, lines)); + + // Do a final seed fill over Voronoi cells and unmarked Voronoi edges. + while (! cell_queue.empty()) { + const VD::cell_type *cell = cell_queue.back(); + const CellCategory cc = cell_category(cell); + assert(cc == CellCategory::Outside || cc == CellCategory::Inside); + cell_queue.pop_back(); + const VD::edge_type *first_edge = cell->incident_edge(); + const VD::edge_type *edge = first_edge; + const auto ec_new = (cc == CellCategory::Outside) ? EdgeCategory::PointsOutside : EdgeCategory::PointsInside; + do { + EdgeCategory ec = edge_category(edge); + if (ec == EdgeCategory::Unknown) { + assert(edge->cell()->contains_point() && edge->twin()->cell()->contains_point()); + annotate_edge(edge, ec_new); + annotate_edge(edge->twin(), ec_new); + const VD::cell_type *cell2 = edge->twin()->cell(); + CellCategory cc2 = cell_category(cell2); + assert(cc2 == CellCategory::Unknown || cc2 == cc); + if (cc2 != cc) { + annotate_cell(cell2, cc); + cell_queue.emplace_back(cell2); + } + } else { + assert(edge->vertex0() == nullptr || vertex_category(edge->vertex0()) != VertexCategory::Unknown); + assert(edge->vertex1() == nullptr || vertex_category(edge->vertex1()) != VertexCategory::Unknown); + assert(edge_category(edge->twin()) != EdgeCategory::Unknown); + assert(cell_category(edge->cell()) != CellCategory::Unknown); + assert(cell_category(edge->twin()->cell()) != CellCategory::Unknown); + } + edge = edge->next(); + } while (edge != first_edge); + } + + assert(debug::verify_vertices_on_contour(vd, lines)); + assert(debug::verify_inside_outside_annotations(vd)); +} + +std::vector signed_vertex_distances(const VD &vd, const Lines &lines) +{ + // vd shall be annotated. + assert(debug::verify_inside_outside_annotations(vd)); + + std::vector out(vd.vertices().size(), 0.); + const VD::vertex_type *first_vertex = &vd.vertices().front(); + for (const VD::vertex_type &vertex : vd.vertices()) { + const VertexCategory vc = vertex_category(vertex); + double dist; + if (vc == VertexCategory::OnContour) { + dist = 0.; + } else { + const VD::edge_type *first_edge = vertex.incident_edge(); + const VD::edge_type *edge = first_edge; + const VD::cell_type *point_cell = nullptr; + do { + if (edge->cell()->contains_point()) { + point_cell = edge->cell(); + break; + } + edge = edge->rot_next(); + } while (edge != first_edge); + if (point_cell == 0) { + // Project vertex onto a contour segment. + const Line &line = lines[edge->cell()->source_index()]; + dist = Geometry::ray_point_distance( + line.a.cast(), (line.b - line.a).cast(), vertex_point(vertex)); + } else { + // Distance to a contour point. + dist = (contour_point(*point_cell, lines).cast() - vertex_point(vertex)).norm(); + } + if (vc == VertexCategory::Inside) + dist = - dist; + } + out[&vertex - first_vertex] = dist; + } + + assert(debug::verify_signed_distances(vd, lines, out)); + + return out; +} + +std::vector edge_offset_contour_intersections( + const VD &vd, + const Lines &lines, + const std::vector &vertex_distances, + double offset_distance) +{ + // vd shall be annotated. + assert(debug::verify_inside_outside_annotations(vd)); + + bool outside = offset_distance > 0; if (! outside) offset_distance = - offset_distance; + assert(offset_distance > 0.); + const VD::vertex_type *first_vertex = &vd.vertices().front(); + const VD::edge_type *first_edge = &vd.edges().front(); + static constexpr double nan = std::numeric_limits::quiet_NaN(); + // By default none edge has an intersection with the offset curve. + std::vector out(vd.num_edges(), Vec2d(nan, 0.)); + + for (const VD::edge_type &edge : vd.edges()) { + size_t edge_idx = &edge - first_edge; + if (edge_offset_has_intersection(out[edge_idx]) || out[edge_idx].y() != 0.) + // This edge was already classified. + continue; + + const VD::vertex_type *v0 = edge.vertex0(); + const VD::vertex_type *v1 = edge.vertex1(); + if (v0 == nullptr) { + assert(vertex_category(v1) == VertexCategory::OnContour || vertex_category(v1) == VertexCategory::Outside); + continue; + } + + double d0 = (v0 == nullptr) ? std::numeric_limits::max() : vertex_distances[v0 - first_vertex]; + double d1 = (v1 == nullptr) ? std::numeric_limits::max() : vertex_distances[v1 - first_vertex]; + assert(d0 * d1 >= 0.); + if (! outside) { + d0 = - d0; + d1 = - d1; + } +#ifndef NDEBUG + { + double err = std::abs(detail::dist_to_site(lines, *edge.cell(), vertex_point(v0)) - std::abs(d0)); + double err2 = std::abs(detail::dist_to_site(lines, *edge.twin()->cell(), vertex_point(v0)) - std::abs(d0)); + assert(err < SCALED_EPSILON); + assert(err2 < SCALED_EPSILON); + if (v1 != nullptr) { + double err3 = std::abs(detail::dist_to_site(lines, *edge.cell(), vertex_point(v1)) - std::abs(d1)); + double err4 = std::abs(detail::dist_to_site(lines, *edge.twin()->cell(), vertex_point(v1)) - std::abs(d1)); + assert(err3 < SCALED_EPSILON); + assert(err4 < SCALED_EPSILON); + } + } +#endif // NDEBUG + double dmin, dmax; + if (d0 < d1) + dmin = d0, dmax = d1; + else + dmax = d0, dmin = d1; + // Offset distance may be lower than dmin, but never higher than dmax. + // Don't intersect an edge at dmax + // 1) To avoid zero edge length, zero area offset contours. + // 2) To ensure that the offset contours that cross a Voronoi vertex are traced consistently + // at one side of the offset curve only. + if (offset_distance >= dmax) + continue; + + // Edge candidate, intersection points were not calculated yet. + assert(v0 != nullptr); + const VD::cell_type *cell = edge.cell(); + const VD::cell_type *cell2 = edge.twin()->cell(); + const Line &line0 = lines[cell->source_index()]; + const Line &line1 = lines[cell2->source_index()]; + size_t edge_idx2 = edge.twin() - first_edge; + if (v1 == nullptr) { + assert(edge.is_infinite()); + assert(edge.is_linear()); + // Unconstrained edges have always montonous distance. + assert(d0 != d1); + if (offset_distance > dmin) { + // There is certainly an intersection with the offset curve. + if (cell->contains_point() && cell2->contains_point()) { + assert(! edge.is_secondary()); + const Point &pt0 = contour_point(*cell, line0); + const Point &pt1 = contour_point(*cell2, line1); + // pt is inside the circle (pt0, offset_distance), (pt + dir) is certainly outside the same circle. + Vec2d dir = Vec2d(double(pt0.y() - pt1.y()), double(pt1.x() - pt0.x())) * (2. * offset_distance); + Vec2d pt(v0->x(), v0->y()); + double t = detail::first_circle_segment_intersection_parameter(Vec2d(pt0.x(), pt0.y()), offset_distance, pt, dir); + assert(t > 0.); + out[edge_idx] = pt + t * dir; + } else { + // Infinite edges could not be created by two segment sites. + assert(cell->contains_point() != cell2->contains_point()); + // Linear edge goes through the endpoint of a segment. + assert(edge.is_secondary()); + const Point &ipt = cell->contains_segment() ? contour_point(*cell2, line1) : contour_point(*cell, line0); + #ifndef NDEBUG + if (cell->contains_segment()) { + const Point &pt1 = contour_point(*cell2, line1); + assert(pt1 == line0.a || pt1 == line0.b); + } else { + const Point &pt0 = contour_point(*cell, line0); + assert(pt0 == line1.a || pt0 == line1.b); + } + assert((vertex_point(v0) - ipt.cast()).norm() < SCALED_EPSILON); + #endif /* NDEBUG */ + // Infinite edge starts at an input contour, therefore there is always an intersection with an offset curve. + const Line &line = cell->contains_segment() ? line0 : line1; + assert(line.a == ipt || line.b == ipt); + out[edge_idx] = ipt.cast() + offset_distance * Vec2d(line.b.y() - line.a.y(), line.a.x() - line.b.x()).normalized(); + } + } else if (offset_distance == dmin) + out[edge_idx] = vertex_point(v0); + // The other edge of an unconstrained edge starting with null vertex shall never be intersected. Mark it as visited. + out[edge_idx2].y() = 1.; + } else { + assert(edge.is_finite()); + bool done = false; + // Bisector of two line segments, distance along the bisector is linear. + bool bisector = cell->contains_segment() && cell2->contains_segment(); + assert(edge.is_finite()); + // or a secondary line, again the distance along the secondary line is linear and starts at the contour (zero distance). + if (bisector || edge.is_secondary()) { + assert(edge.is_linear()); +#ifndef NDEBUG + if (edge.is_secondary()) { + assert(cell->contains_point() != cell2->contains_point()); + // One of the vertices is on the input contour. + assert((vertex_category(edge.vertex0()) == VertexCategory::OnContour) != + (vertex_category(edge.vertex1()) == VertexCategory::OnContour)); + assert(dmin == 0.); + } +#endif // NDEBUG + if (! bisector || (dmin != dmax && offset_distance >= dmin)) { + double t = (offset_distance - dmin) / (dmax - dmin); + t = clamp(0., 1., t); + if (d1 < d0) { + out[edge_idx2] = Slic3r::lerp(vertex_point(v1), vertex_point(v0), t); + // mark visited + out[edge_idx].y() = 1.; + } else { + out[edge_idx] = Slic3r::lerp(vertex_point(v0), vertex_point(v1), t); + // mark visited + out[edge_idx2].y() = 1.; + } + done = true; + } + } else { + // Point - Segment or Point - Point edge, distance along this Voronoi edge may not be monotonous: + // The distance function along the Voronoi edge may either have one maximum at one vertex and one minimum at the other vertex, + // or it may have two maxima at the vertices and a minimum somewhere along the Voronoi edge, and this Voronoi edge + // may be intersected twice by an offset curve. + // + // Tracing an offset curve accross Voronoi regions with linear edges of montonously increasing or decrasing distance + // to a Voronoi region is stable in a sense, that if the distance of Voronoi vertices is calculated correctly, there will + // be maximum one intersection of an offset curve found at each Voronoi edge and tracing these intersections shall + // produce a set of closed curves. + // + // Having a non-monotonous distance function between the Voronoi edge end points may lead to splitting of offset curves + // at these Voronoi edges. If a Voronoi edge is classified as having no intersection at all while it has some, + // the extracted offset curve will contain self intersections at this Voronoi edge. + // + // If on the other side the two intersection points are found by a numerical error even though none should be found, then + // it may happen that it would not be possible to connect these two points into a closed loop, which is likely worse + // than the issue above. + // + // While it is likely not possible to avoid all the numerical issues, one shall strive for the highest numerical robustness. + assert(cell->contains_point() || cell2->contains_point()); + size_t num_intersections = 0; + bool point_vs_segment = cell->contains_point() != cell2->contains_point(); + const Point &pt0 = cell->contains_point() ? contour_point(*cell, line0) : contour_point(*cell2, line1); + // Project p0 to line segment . + Vec2d p0(v0->x(), v0->y()); + Vec2d p1(v1->x(), v1->y()); + Vec2d px(pt0.x(), pt0.y()); + const Point *pt1 = nullptr; + Vec2d dir; + if (point_vs_segment) { + const Line &line = cell->contains_segment() ? line0 : line1; + dir = (line.b - line.a).cast(); + } else { + pt1 = &contour_point(*cell2, line1); + // Perpendicular to the (pt1 - pt0) direction. + dir = Vec2d(double(pt0.y() - pt1->y()), double(pt1->x() - pt0.x())); + } + double s0 = (p0 - px).dot(dir); + double s1 = (p1 - px).dot(dir); + if (offset_distance >= dmin) { + // This Voronoi edge is intersected by the offset curve just once. + // There may be numerical issues if dmin is close to the minimum of the non-monotonous distance function. + num_intersections = 1; + } else { + // This Voronoi edge may not be intersected by the offset curve, or it may split the offset curve + // into two loops. First classify this edge robustly with regard to the Point-Segment bisector or Point-Point bisector. + double dmin_new; + bool found = false; + if (point_vs_segment) { + if (s0 * s1 <= 0.) { + // One end of the Voronoi edge is on one side of the Point-Segment bisector, the other end of the Voronoi + // edge is on the other side of the bisector, therefore with a high probability we should find a minimum + // of the distance to a nearest site somewhere inside this Voronoi edge (at the intersection of the bisector + // and the Voronoi edge. + const Line &line = cell->contains_segment() ? line0 : line1; + dmin_new = 0.5 * (Geometry::foot_pt(line.a.cast(), dir, px) - px).norm(); + found = true; + } + } else { + // Point-Point Voronoi sites. + if (s0 * s1 <= 0.) { + // This Voronoi edge intersects connection line of the two Point sites. + dmin_new = 0.5 * (pt1->cast() - px).norm(); + found = true; + } + } + if (found) { + assert(dmin_new < dmax + SCALED_EPSILON); + assert(dmin_new < dmin + SCALED_EPSILON); + if (dmin_new <= offset_distance) { + // 1) offset_distance > dmin_new -> two new distinct intersection points are found. + // 2) offset_distance == dmin_new -> one duplicate point is found. + // If 2) is ignored, then two tangentially touching offset curves are created. + // If not ignored, then the two offset curves merge at this double point. + // We should merge the contours while pushing the the two copies of the tangent point away a bit. + dmin = dmin_new; + num_intersections = (offset_distance > dmin) + 1; + } + } + } + if (num_intersections > 0) { + detail::Intersections intersections; + if (point_vs_segment) { + assert(cell->contains_point() || cell2->contains_point()); + intersections = detail::line_point_equal_distance_points(cell->contains_segment() ? line0 : line1, pt0, offset_distance); + } else { + intersections = detail::point_point_equal_distance_points(pt0, *pt1, offset_distance); + } + // The functions above perform complex calculations in doubles, therefore the results may not be quite accurate and + // the number of intersections found may not be in accord to the number of intersections expected from evaluating simpler expressions. + // Adjust the result to the number of intersection points expected. + if (num_intersections == 2) { + switch (intersections.count) { + case 0: + // No intersection found even though one or two were expected to be found. + // Not trying to find the intersection means that we may produce offset curves, that intersect at this Voronoi edge. + //FIXME We are fine with that for now, but we may try to create artificial split points in further revisions. + break; + case 1: + // Tangential point found. + //FIXME We are fine with that for now, but we may try to create artificial split points in further revisions. + break; + default: + { + // Two intersection points found. Sort them. + assert(intersections.count == 2); + double q0 = (intersections.pts[0] - px).dot(dir); + double q1 = (intersections.pts[1] - px).dot(dir); + // Both Voronoi edge end points and offset contour intersection points should be separated by the bisector. + assert(q0 * q1 <= 0.); + assert(s0 * s1 <= 0.); + // Sort the intersection points by dir. + if ((q0 < q1) != (s0 < s1)) + std::swap(intersections.pts[0], intersections.pts[1]); + } + } + } else { + assert(num_intersections == 1); + switch (intersections.count) { + case 0: + // No intersection found. This should not happen. + // Create one artificial intersection point by repeating the dmin point, which is supposed to be + // close to the minimum. + intersections.pts[0] = (dmin == d0) ? p0 : p1; + intersections.count = 1; + break; + case 1: + // One intersection found. This is a tangential point. Use it. + break; + default: + // Two intersections found. + // Now decide which of the point fall on this Voronoi edge. + assert(intersections.count == 2); + double q0 = (intersections.pts[0] - px).dot(dir); + double q1 = (intersections.pts[1] - px).dot(dir); + // Offset contour intersection points should be separated by the bisector. + assert(q0 * q1 <= 0); + double s = (dmax == d0) ? s0 : s1; + bool take_2nd = (s > 0.) ? q1 > q0 : q1 < q0; + if (take_2nd) + intersections.pts[0] = intersections.pts[1]; + -- intersections.count; + } + } + assert(intersections.count > 0); + if (intersections.count == 2) { + out[edge_idx] = intersections.pts[1]; + out[edge_idx2] = intersections.pts[0]; + done = true; + } else if (intersections.count == 1) { + if (d1 < d0) + std::swap(edge_idx, edge_idx2); + out[edge_idx] = intersections.pts[0]; + out[edge_idx2].y() = 1.; + done = true; + } + } + } + if (! done) + out[edge_idx].y() = out[edge_idx2].y() = 1.; + } + } + + assert(debug::verify_offset_intersection_points(vd, lines, offset_distance, out)); + + return out; +} + +Polygons offset( + const Geometry::VoronoiDiagram &vd, + const Lines &lines, + const std::vector &signed_vertex_distances, + double offset_distance, + double discretization_error) +{ #ifdef VORONOI_DEBUG_OUT BoundingBox bbox; { @@ -442,10 +1339,10 @@ Polygons voronoi_offset( { Lines helper_lines; for (const VD::edge_type &edge : vd.edges()) - if (edge_state[&edge - front_edge] == EdgeState::Possible) { + if (edge_category(edge) == (offset_distance > 0 ? EdgeCategory::PointsOutside : EdgeCategory::PointsInside) && + edge.vertex0() != nullptr) { const VD::vertex_type *v0 = edge.vertex0(); const VD::vertex_type *v1 = edge.vertex1(); - assert(v0 != nullptr); Vec2d pt1(v0->x(), v0->y()); Vec2d pt2; if (v1 == nullptr) { @@ -490,294 +1387,39 @@ Polygons voronoi_offset( } #endif // VORONOI_DEBUG_OUT - std::vector edge_offset_point(vd.num_edges(), Vec2d()); - const double offset_distance2 = offset_distance * offset_distance; - for (const VD::edge_type &edge : vd.edges()) { - assert(edge_state[&edge - front_edge] != EdgeState::Unknown); - size_t edge_idx = &edge - front_edge; - if (edge_state[edge_idx] == EdgeState::Possible) { - // Edge candidate, intersection points were not calculated yet. - const VD::vertex_type *v0 = edge.vertex0(); - const VD::vertex_type *v1 = edge.vertex1(); - assert(v0 != nullptr); - const VD::cell_type *cell = edge.cell(); - const VD::cell_type *cell2 = edge.twin()->cell(); - const Line &line0 = lines[cell->source_index()]; - const Line &line1 = lines[cell2->source_index()]; - size_t edge_idx2 = edge.twin() - front_edge; - if (v1 == nullptr) { - assert(edge.is_infinite()); - assert(edge.is_linear()); - assert(edge_state[edge_idx2] == EdgeState::Inactive); - if (cell->contains_point() && cell2->contains_point()) { - assert(! edge.is_secondary()); - const Point &pt0 = (cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b; - const Point &pt1 = (cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b; - double dmin2 = (Vec2d(v0->x(), v0->y()) - pt0.cast()).squaredNorm(); - assert(dmin2 >= SCALED_EPSILON * SCALED_EPSILON); - if (dmin2 <= offset_distance2) { - // There shall be an intersection of this unconstrained edge with the offset curve. - // Direction vector of this unconstrained Voronoi edge. - Vec2d dir(double(pt0.y() - pt1.y()), double(pt1.x() - pt0.x())); - Vec2d pt(v0->x(), v0->y()); - double t = detail::first_circle_segment_intersection_parameter(Vec2d(pt0.x(), pt0.y()), offset_distance, pt, dir); - edge_offset_point[edge_idx] = pt + t * dir; - set_edge_state_final(edge_idx, EdgeState::Active); - } else - set_edge_state_final(edge_idx, EdgeState::Inactive); - } else { - // Infinite edges could not be created by two segment sites. - assert(cell->contains_point() != cell2->contains_point()); - // Linear edge goes through the endpoint of a segment. - assert(edge.is_secondary()); - const Point &ipt = cell->contains_segment() ? - ((cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b) : - ((cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b); - #ifndef NDEBUG - if (cell->contains_segment()) { - const Point &pt1 = (cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b; - assert((pt1.x() == line0.a.x() && pt1.y() == line0.a.y()) || - (pt1.x() == line0.b.x() && pt1.y() == line0.b.y())); - } else { - const Point &pt0 = (cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b; - assert((pt0.x() == line1.a.x() && pt0.y() == line1.a.y()) || - (pt0.x() == line1.b.x() && pt0.y() == line1.b.y())); - } - assert((Vec2d(v0->x(), v0->y()) - ipt.cast()).norm() < SCALED_EPSILON); - #endif /* NDEBUG */ - // Infinite edge starts at an input contour, therefore there is always an intersection with an offset curve. - const Line &line = cell->contains_segment() ? line0 : line1; - assert(line.a == ipt || line.b == ipt); - edge_offset_point[edge_idx] = ipt.cast() + offset_distance * Vec2d(line.b.y() - line.a.y(), line.a.x() - line.b.x()).normalized(); - set_edge_state_final(edge_idx, EdgeState::Active); - } - // The other edge of an unconstrained edge starting with null vertex shall never be intersected. - set_edge_state_final(edge_idx2, EdgeState::Inactive); - } else if (edge.is_secondary()) { - assert(edge.is_linear()); - assert(cell->contains_point() != cell2->contains_point()); - const Line &line0 = lines[edge.cell()->source_index()]; - const Line &line1 = lines[edge.twin()->cell()->source_index()]; - const Point &pt = cell->contains_point() ? - ((cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b) : - ((cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b); - #ifndef NDEBUG - const Line &line = cell->contains_segment() ? line0 : line1; - assert(pt == line.a || pt == line.b); - assert((pt.cast() - Vec2d(v0->x(), v0->y())).norm() < SCALED_EPSILON); - #endif // NDEBUG - Vec2d dir(v1->x() - v0->x(), v1->y() - v0->y()); - double l2 = dir.squaredNorm(); - if (offset_distance2 <= l2) { - edge_offset_point[edge_idx] = pt.cast() + (offset_distance / sqrt(l2)) * dir; - set_edge_state_final(edge_idx, EdgeState::Active); - } else { - set_edge_state_final(edge_idx, EdgeState::Inactive); - } - set_edge_state_final(edge_idx2, EdgeState::Inactive); - } else { - // Finite edge has valid points at both sides. - bool done = false; - if (cell->contains_segment() && cell2->contains_segment()) { - // This edge is a bisector of two line segments. Project v0, v1 onto one of the line segments. - Vec2d pt(line0.a.cast()); - Vec2d dir(line0.b.cast() - pt); - Vec2d vec0 = Vec2d(v0->x(), v0->y()) - pt; - Vec2d vec1 = Vec2d(v1->x(), v1->y()) - pt; - double l2 = dir.squaredNorm(); - assert(l2 > 0.); - double dmin = (dir * (vec0.dot(dir) / l2) - vec0).squaredNorm(); - double dmax = (dir * (vec1.dot(dir) / l2) - vec1).squaredNorm(); - bool flip = dmin > dmax; - if (flip) - std::swap(dmin, dmax); - if (offset_distance2 >= dmin && offset_distance2 <= dmax) { - // Intersect. Maximum one intersection will be found. - // This edge is a bisector of two line segments. Distance to the input polygon increases/decreases monotonically. - dmin = sqrt(dmin); - dmax = sqrt(dmax); - assert(offset_distance > dmin - EPSILON && offset_distance < dmax + EPSILON); - double ddif = dmax - dmin; - if (ddif == 0.) { - // line, line2 are exactly parallel. This is a singular case, the offset curve should miss it. - } else { - if (flip) { - std::swap(edge_idx, edge_idx2); - std::swap(v0, v1); - } - double t = clamp(0., 1., (offset_distance - dmin) / ddif); - edge_offset_point[edge_idx] = Vec2d(lerp(v0->x(), v1->x(), t), lerp(v0->y(), v1->y(), t)); - set_edge_state_final(edge_idx, EdgeState::Active); - set_edge_state_final(edge_idx2, EdgeState::Inactive); - done = true; - } - } - } else { - assert(cell->contains_point() || cell2->contains_point()); - bool point_vs_segment = cell->contains_point() != cell2->contains_point(); - const Point &pt0 = cell->contains_point() ? - ((cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b) : - ((cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b); - // Project p0 to line segment . - Vec2d p0(v0->x(), v0->y()); - Vec2d p1(v1->x(), v1->y()); - Vec2d px(pt0.x(), pt0.y()); - double d0 = (p0 - px).squaredNorm(); - double d1 = (p1 - px).squaredNorm(); - double dmin = std::min(d0, d1); - double dmax = std::max(d0, d1); - bool has_intersection = false; - bool possibly_two_points = false; - if (offset_distance2 <= dmax) { - if (offset_distance2 >= dmin) { - has_intersection = true; - } else { - double dmin_new = dmin; - if (point_vs_segment) { - // Project on the source segment. - const Line &line = cell->contains_segment() ? line0 : line1; - const Vec2d pt_line = line.a.cast(); - const Vec2d v_line = (line.b - line.a).cast(); - double t0 = (p0 - pt_line).dot(v_line); - double t1 = (p1 - pt_line).dot(v_line); - double tx = (px - pt_line).dot(v_line); - if ((tx >= t0 && tx <= t1) || (tx >= t1 && tx <= t0)) { - // Projection of the Point site falls between the projections of the Voronoi edge end points - // onto the Line site. - Vec2d ft = pt_line + (tx / v_line.squaredNorm()) * v_line; - dmin_new = (ft - px).squaredNorm() * 0.25; - } - } else { - // Point-Point Voronoi sites. Project point site onto the current Voronoi edge. - Vec2d v = p1 - p0; - auto l2 = v.squaredNorm(); - assert(l2 > 0); - auto t = v.dot(px - p0); - if (t >= 0. && t <= l2) { - // Projection falls onto the Voronoi edge. Calculate foot point and distance. - Vec2d ft = p0 + (t / l2) * v; - dmin_new = (ft - px).squaredNorm(); - } - } - assert(dmin_new < dmax + SCALED_EPSILON); - assert(dmin_new < dmin + SCALED_EPSILON); - if (dmin_new < dmin) { - dmin = dmin_new; - has_intersection = possibly_two_points = offset_distance2 >= dmin; - } - } - } - if (has_intersection) { - detail::Intersections intersections; - if (point_vs_segment) { - assert(cell->contains_point() || cell2->contains_point()); - intersections = detail::line_point_equal_distance_points(cell->contains_segment() ? line0 : line1, pt0, offset_distance); - } else { - const Point &pt1 = (cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b; - intersections = detail::point_point_equal_distance_points(pt0, pt1, offset_distance); - } - // If the span of distances of start / end point / foot point to the point site indicate an intersection, - // we should find one. - assert(intersections.count > 0); - if (intersections.count == 2) { - // Now decide which points fall on this Voronoi edge. - // Tangential points (single intersection) are ignored. - if (possibly_two_points) { - Vec2d v = p1 - p0; - double l2 = v.squaredNorm(); - double t0 = v.dot(intersections.pts[0] - p0); - double t1 = v.dot(intersections.pts[1] - p0); - if (t0 > t1) { - std::swap(t0, t1); - std::swap(intersections.pts[0], intersections.pts[1]); - } - // Remove points outside of the line range. - if (t0 < 0. || t0 > l2) { - if (t1 < 0. || t1 > l2) - intersections.count = 0; - else { - -- intersections.count; - t0 = t1; - intersections.pts[0] = intersections.pts[1]; - } - } else if (t1 < 0. || t1 > l2) - -- intersections.count; - } else { - // Take the point furthest from the end points of the Voronoi edge or a Voronoi parabolic arc. - double d0 = std::max((intersections.pts[0] - p0).squaredNorm(), (intersections.pts[0] - p1).squaredNorm()); - double d1 = std::max((intersections.pts[1] - p0).squaredNorm(), (intersections.pts[1] - p1).squaredNorm()); - if (d0 > d1) - intersections.pts[0] = intersections.pts[1]; - -- intersections.count; - } - assert(intersections.count > 0); - if (intersections.count == 2) { - set_edge_state_final(edge_idx, EdgeState::Active); - set_edge_state_final(edge_idx2, EdgeState::Active); - edge_offset_point[edge_idx] = intersections.pts[1]; - edge_offset_point[edge_idx2] = intersections.pts[0]; - done = true; - } else if (intersections.count == 1) { - if (d1 < d0) - std::swap(edge_idx, edge_idx2); - set_edge_state_final(edge_idx, EdgeState::Active); - set_edge_state_final(edge_idx2, EdgeState::Inactive); - edge_offset_point[edge_idx] = intersections.pts[0]; - done = true; - } - } - } - } - if (! done) { - set_edge_state_final(edge_idx, EdgeState::Inactive); - set_edge_state_final(edge_idx2, EdgeState::Inactive); - } - } - } - } - -#ifndef NDEBUG - for (const VD::edge_type &edge : vd.edges()) { - assert(edge_state[&edge - front_edge] == EdgeState::Inactive || edge_state[&edge - front_edge] == EdgeState::Active); - // None of a new edge candidate may start with null vertex. - assert(edge_state[&edge - front_edge] == EdgeState::Inactive || edge.vertex0() != nullptr); - assert(edge_state[edge.twin() - front_edge] == EdgeState::Inactive || edge.twin()->vertex0() != nullptr); - } -#endif // NDEBUG + std::vector edge_points = edge_offset_contour_intersections(vd, lines, signed_vertex_distances, offset_distance); + const VD::edge_type *front_edge = &vd.edges().front(); #ifdef VORONOI_DEBUG_OUT + Lines helper_lines; { - Lines helper_lines; for (const VD::edge_type &edge : vd.edges()) - if (edge_state[&edge - front_edge] == EdgeState::Active) - helper_lines.emplace_back(Line(Point(edge.vertex0()->x(), edge.vertex0()->y()), Point(edge_offset_point[&edge - front_edge].cast()))); + if (edge_offset_has_intersection(edge_points[&edge - front_edge])) + helper_lines.emplace_back(Line(Point(edge.vertex0()->x(), edge.vertex0()->y()), Point(edge_points[&edge - front_edge].cast()))); dump_voronoi_to_svg(debug_out_path("voronoi-offset-candidates2-%d.svg", irun).c_str(), vd, Points(), lines, Polygons(), helper_lines); } #endif // VORONOI_DEBUG_OUT - auto next_offset_edge = [&edge_state, front_edge](const VD::edge_type *start_edge) -> const VD::edge_type* { + auto next_offset_edge = [&edge_points, front_edge](const VD::edge_type *start_edge) -> const VD::edge_type* { for (const VD::edge_type *edge = start_edge->next(); edge != start_edge; edge = edge->next()) - if (edge_state[edge->twin() - front_edge] == EdgeState::Active) + if (edge_offset_has_intersection(edge_points[edge->twin() - front_edge])) return edge->twin(); // assert(false); return nullptr; }; -#ifndef NDEBUG - auto dist_to_site = [&lines](const VD::cell_type &cell, const Vec2d &point) { - const Line &line = lines[cell.source_index()]; - return cell.contains_point() ? - (((cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line.a : line.b).cast() - point).norm() : - (Geometry::foot_pt(line.a.cast(), (line.b - line.a).cast(), point) - point).norm(); - }; -#endif /* NDEBUG */ + const bool inside_offset = offset_distance < 0.; + if (inside_offset) + offset_distance = - offset_distance; - // Track the offset curves. - Polygons out; + // Track the offset curves. + Polygons out; double angle_step = 2. * acos((offset_distance - discretization_error) / offset_distance); double cos_threshold = cos(angle_step); - for (size_t seed_edge_idx = 0; seed_edge_idx < vd.num_edges(); ++ seed_edge_idx) - if (edge_state[seed_edge_idx] == EdgeState::Active) { + static constexpr double nan = std::numeric_limits::quiet_NaN(); + for (size_t seed_edge_idx = 0; seed_edge_idx < vd.num_edges(); ++ seed_edge_idx) { + Vec2d last_pt = edge_points[seed_edge_idx]; + if (edge_offset_has_intersection(last_pt)) { const VD::edge_type *start_edge = &vd.edges()[seed_edge_idx]; const VD::edge_type *edge = start_edge; Polygon poly; @@ -786,21 +1428,23 @@ Polygons voronoi_offset( const VD::edge_type *next_edge = next_offset_edge(edge); #ifdef VORONOI_DEBUG_OUT if (next_edge == nullptr) { - Lines helper_lines; - dump_voronoi_to_svg(debug_out_path("voronoi-offset-open-loop-%d.svg", irun).c_str(), vd, Points(), lines, Polygons(), to_lines(poly)); + Lines hl = helper_lines; + append(hl, to_lines(Polyline(poly.points))); + dump_voronoi_to_svg(debug_out_path("voronoi-offset-open-loop-%d.svg", irun).c_str(), vd, Points(), lines, Polygons(), hl); } #endif // VORONOI_DEBUG_OUT assert(next_edge); //std::cout << "offset-output: "; print_edge(edge); std::cout << " to "; print_edge(next_edge); std::cout << "\n"; // Interpolate a circular segment or insert a linear segment between edge and next_edge. const VD::cell_type *cell = edge->cell(); - edge_state[next_edge - front_edge] = EdgeState::Inactive; - Vec2d p1 = edge_offset_point[edge - front_edge]; - Vec2d p2 = edge_offset_point[next_edge - front_edge]; + // Mark the edge / offset curve intersection point as consumed. + Vec2d p1 = last_pt; + Vec2d p2 = edge_points[next_edge - front_edge]; + edge_points[next_edge - front_edge].x() = nan; #ifndef NDEBUG { - double err = dist_to_site(*cell, p1) - offset_distance; - double err2 = dist_to_site(*cell, p2) - offset_distance; + double err = detail::dist_to_site(lines, *cell, p1) - offset_distance; + double err2 = detail::dist_to_site(lines, *cell, p2) - offset_distance; #ifdef VORONOI_DEBUG_OUT if (std::max(err, err2) >= SCALED_EPSILON) { Lines helper_lines; @@ -847,11 +1491,142 @@ Polygons voronoi_offset( poly.points.emplace_back(pt_last); } edge = next_edge; + last_pt = p2; } while (edge != start_edge); - out.emplace_back(std::move(poly)); + + while (! poly.empty() && poly.points.front() == poly.points.back()) + poly.points.pop_back(); + if (poly.size() >= 3) + out.emplace_back(std::move(poly)); } + } return out; } +Polygons offset( + const VD &vd, + const Lines &lines, + double offset_distance, + double discretization_error) +{ + annotate_inside_outside(const_cast(vd), lines); + std::vector dist = signed_vertex_distances(vd, lines); + return offset(vd, lines, dist, offset_distance, discretization_error); +} + +// Produce a list of start positions of a skeleton segment at a halfedge. +// If the whole Voronoi edge is part of the skeleton, then zero start positions are assigned +// to both ends of the edge. Position "1" shall never be assigned to a halfedge. +// +// Skeleton edges must be inside a closed polygon, therefore these edges are finite. +// A Voronoi Edge-Edge bisector is either completely part of a skeleton, or not at all. +// An infinite Voronoi Edge-Point (parabola) or Point-Point (line) bisector is split into +// a center part close to the Voronoi sites (not skeleton) and the ends (skeleton), +// though either part could be clipped by the Voronoi segment. +// +// Further filtering of the skeleton may be necessary. +std::vector skeleton_edges_rough( + const VD &vd, + const Lines &lines, + // Angle threshold at a sharp convex corner, which is marked for a gap fill. + const double threshold_alpha) +{ + // vd shall be annotated. + assert(debug::verify_inside_outside_annotations(vd)); + + const VD::edge_type *first_edge = &vd.edges().front(); + static constexpr double nan = std::numeric_limits::quiet_NaN(); + // By default no edge is annotated as being part of the skeleton. + std::vector out(vd.num_edges(), Vec2d(nan, nan)); + // Threshold at a sharp corner, derived from a dot product of the sharp corner edges. + const double threshold_cos_alpha = cos(threshold_alpha); + // For sharp corners, dr/dl = sin(alpha/2). Substituting the dr/dl threshold with tan(alpha/2) threshold + // in Voronoi point - point and Voronoi point - line site functions. + const double threshold_tan_alpha_half = tan(0.5 * threshold_alpha); + + for (const VD::edge_type &edge : vd.edges()) { + size_t edge_idx = &edge - first_edge; + if ( + // Ignore secondary and unbounded edges, they shall never be part of the skeleton. + edge.is_secondary() || edge.is_infinite() || + // Skip the twin edge of an edge, that has already been processed. + &edge > edge.twin() || + // Ignore outer edges. + (edge_category(edge) != EdgeCategory::PointsInside && edge_category(edge.twin()) != EdgeCategory::PointsInside)) + continue; + const VD::vertex_type *v0 = edge.vertex0(); + const VD::vertex_type *v1 = edge.vertex1(); + const VD::cell_type *cell = edge.cell(); + const VD::cell_type *cell2 = edge.twin()->cell(); + const Line &line0 = lines[cell->source_index()]; + const Line &line1 = lines[cell2->source_index()]; + size_t edge_idx2 = edge.twin() - first_edge; + if (cell->contains_segment() && cell2->contains_segment()) { + // Bisector of two line segments, distance along the bisector is linear, + // dr/dl is constant. + // using trigonometric identity sin^2(a) = (1-cos(2a))/2 + Vec2d lv0 = (line0.b - line0.a).cast(); + Vec2d lv1 = (line1.b - line1.a).cast(); + double d = lv0.dot(lv1); + if (d < 0.) { + double cos_alpha = - d / (lv0.norm() * lv1.norm()); + if (cos_alpha > threshold_cos_alpha) { + // The whole bisector is a skeleton segment. + out[edge_idx] = vertex_point(v0); + out[edge_idx2] = vertex_point(v1); + } + } + } else { + // An infinite Voronoi Edge-Point (parabola) or Point-Point (line) bisector, clipped to a finite Voronoi segment. + // The infinite bisector has a distance (skeleton radius) minimum, which is also a minimum + // of the skeleton function dr / dt. + assert(cell->contains_point() || cell2->contains_point()); + if (cell->contains_point() != cell2->contains_point()) { + // Point - Segment + const Point &pt0 = cell->contains_point() ? contour_point(*cell, line0) : contour_point(*cell2, line1); + const Line &line = cell->contains_segment() ? line0 : line1; + std::tie(out[edge_idx], out[edge_idx2]) = detail::point_segment_dr_dl_thresholds( + pt0, line, vertex_point(v0), vertex_point(v1), threshold_tan_alpha_half); + } else { + // Point - Point + const Point &pt0 = contour_point(*cell, line0); + const Point &pt1 = contour_point(*cell2, line1); + std::tie(out[edge_idx], out[edge_idx2]) = detail::point_point_dr_dl_thresholds( + pt0, pt1, vertex_point(v0), vertex_point(v1), threshold_tan_alpha_half); + } + } + } + +#ifdef VORONOI_DEBUG_OUT + { + static int irun = 0; + ++ irun; + Lines helper_lines; + for (const VD::edge_type &edge : vd.edges()) + if (&edge < edge.twin() && edge.is_finite()) { + const Vec2d &skeleton_pt = out[&edge - first_edge]; + const Vec2d &skeleton_pt2 = out[edge.twin() - first_edge]; + bool has_skeleton_pt = ! std::isnan(skeleton_pt.x()); + bool has_skeleton_pt2 = ! std::isnan(skeleton_pt2.x()); + const Vec2d &vertex_pt = vertex_point(edge.vertex0()); + const Vec2d &vertex_pt2 = vertex_point(edge.vertex1()); + if (has_skeleton_pt && has_skeleton_pt2) { + // Complete edge is part of the skeleton. + helper_lines.emplace_back(Line(Point(vertex_pt.x(), vertex_pt.y()), Point(vertex_pt2.x(), vertex_pt2.y()))); + } else { + if (has_skeleton_pt) + helper_lines.emplace_back(Line(Point(vertex_pt2.x(), vertex_pt2.y()), Point(skeleton_pt.x(), skeleton_pt.y()))); + if (has_skeleton_pt2) + helper_lines.emplace_back(Line(Point(vertex_pt.x(), vertex_pt.y()), Point(skeleton_pt2.x(), skeleton_pt2.y()))); + } + } + dump_voronoi_to_svg(debug_out_path("voronoi-skeleton-edges-%d.svg", irun).c_str(), vd, Points(), lines, Polygons(), helper_lines); + } +#endif // VORONOI_DEBUG_OUT + + return out; +} + +} // namespace Voronoi } // namespace Slic3r diff --git a/src/libslic3r/VoronoiOffset.hpp b/src/libslic3r/VoronoiOffset.hpp index a21b44f93..17b590ed7 100644 --- a/src/libslic3r/VoronoiOffset.hpp +++ b/src/libslic3r/VoronoiOffset.hpp @@ -9,16 +9,136 @@ namespace Slic3r { +namespace Voronoi { + +using VD = Slic3r::Geometry::VoronoiDiagram; + +inline const Point& contour_point(const VD::cell_type &cell, const Line &line) + { return ((cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line.a : line.b); } +inline Point& contour_point(const VD::cell_type &cell, Line &line) + { return ((cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line.a : line.b); } + +inline const Point& contour_point(const VD::cell_type &cell, const Lines &lines) + { return contour_point(cell, lines[cell.source_index()]); } +inline Point& contour_point(const VD::cell_type &cell, Lines &lines) + { return contour_point(cell, lines[cell.source_index()]); } + +inline Vec2d vertex_point(const VD::vertex_type &v) { return Vec2d(v.x(), v.y()); } +inline Vec2d vertex_point(const VD::vertex_type *v) { return Vec2d(v->x(), v->y()); } + +// "Color" stored inside the boost::polygon Voronoi vertex. +enum class VertexCategory : unsigned char +{ + // Voronoi vertex is on the input contour. + // VD::vertex_type stores coordinates in double, though the coordinates shall match exactly + // with the coordinates of the input contour when converted to int32_t. + OnContour, + // Vertex is inside the CCW input contour, holes are respected. + Inside, + // Vertex is outside the CCW input contour, holes are respected. + Outside, + // Not known yet. + Unknown, +}; + +// "Color" stored inside the boost::polygon Voronoi edge. +// The Voronoi edge as represented by boost::polygon Voronoi module is really a half-edge, +// the half-edges are classified based on the target vertex (VD::vertex_type::vertex1()) +enum class EdgeCategory : unsigned char +{ + // This half-edge points onto the contour, this VD::edge_type::vertex1().color() is OnContour. + PointsToContour, + // This half-edge points inside, this VD::edge_type::vertex1().color() is Inside. + PointsInside, + // This half-edge points outside, this VD::edge_type::vertex1().color() is Outside. + PointsOutside, + // Not known yet. + Unknown +}; + +// "Color" stored inside the boost::polygon Voronoi cell. +enum class CellCategory : unsigned char +{ + // This Voronoi cell is split by an input segment to two halves, one is inside, the other is outside. + Boundary, + // This Voronoi cell is completely inside. + Inside, + // This Voronoi cell is completely outside. + Outside, + // Not known yet. + Unknown +}; + +inline VertexCategory vertex_category(const VD::vertex_type &v) + { return static_cast(v.color()); } +inline VertexCategory vertex_category(const VD::vertex_type *v) + { return static_cast(v->color()); } +inline void set_vertex_category(VD::vertex_type &v, VertexCategory c) + { v.color(static_cast(c)); } +inline void set_vertex_category(VD::vertex_type *v, VertexCategory c) + { v->color(static_cast(c)); } + +inline EdgeCategory edge_category(const VD::edge_type &e) + { return static_cast(e.color()); } +inline EdgeCategory edge_category(const VD::edge_type *e) + { return static_cast(e->color()); } +inline void set_edge_category(VD::edge_type &e, EdgeCategory c) + { e.color(static_cast(c)); } +inline void set_edge_category(VD::edge_type *e, EdgeCategory c) + { e->color(static_cast(c)); } + +inline CellCategory cell_category(const VD::cell_type &v) + { return static_cast(v.color()); } +inline CellCategory cell_category(const VD::cell_type *v) + { return static_cast(v->color()); } +inline void set_cell_category(const VD::cell_type &v, CellCategory c) + { v.color(static_cast(c)); } +inline void set_cell_category(const VD::cell_type *v, CellCategory c) + { v->color(static_cast(c)); } + +// Mark the "Color" of VD vertices, edges and cells as Unknown. +void reset_inside_outside_annotations(VD &vd); + +// Assign "Color" to VD vertices, edges and cells signifying whether the entity is inside or outside +// the input polygons defined by Lines. +void annotate_inside_outside(VD &vd, const Lines &lines); + +// Returns a signed distance to Voronoi vertices from the input polygons. +// (negative distances inside, positive distances outside). +std::vector signed_vertex_distances(const VD &vd, const Lines &lines); + +static inline bool edge_offset_no_intersection(const Vec2d &intersection_point) + { return std::isnan(intersection_point.x()); } +static inline bool edge_offset_has_intersection(const Vec2d &intersection_point) + { return ! edge_offset_no_intersection(intersection_point); } +std::vector edge_offset_contour_intersections( + const VD &vd, const Lines &lines, const std::vector &distances, + double offset_distance); + +std::vector skeleton_edges_rough( + const VD &vd, + const Lines &lines, + const double threshold_alpha); + +Polygons offset( + const Geometry::VoronoiDiagram &vd, + const Lines &lines, + const std::vector &signed_vertex_distances, + double offset_distance, + double discretization_error); + // Offset a polygon or a set of polygons possibly with holes by traversing a Voronoi diagram. // The input polygons are stored in lines and lines are referenced by vd. // Outer curve will be extracted for a positive offset_distance, // inner curve will be extracted for a negative offset_distance. // Circular arches will be discretized to achieve discretization_error. -Polygons voronoi_offset( - const Geometry::VoronoiDiagram &vd, - const Lines &lines, - double offset_distance, - double discretization_error); +Polygons offset( + const VD &vd, + const Lines &lines, + double offset_distance, + double discretization_error); + +} // namespace Voronoi } // namespace Slic3r diff --git a/src/libslic3r/VoronoiVisualUtils.hpp b/src/libslic3r/VoronoiVisualUtils.hpp index fa6a34241..177f27293 100644 --- a/src/libslic3r/VoronoiVisualUtils.hpp +++ b/src/libslic3r/VoronoiVisualUtils.hpp @@ -5,6 +5,8 @@ #include #include +#include "VoronoiOffset.hpp" + namespace boost { namespace polygon { // The following code for the visualization of the boost Voronoi diagram is based on: @@ -70,6 +72,7 @@ class voronoi_visual_utils { get_point_projection((*discretization)[0], segment); CT projection_end = sqr_segment_length * get_point_projection((*discretization)[1], segment); + assert(projection_start != projection_end); // Compute parabola parameters in the transformed space. // Parabola has next representation: @@ -99,13 +102,16 @@ class voronoi_visual_utils { // furthest from the current line segment. CT mid_x = (new_y - cur_y) / (new_x - cur_x) * rot_y + rot_x; CT mid_y = parabola_y(mid_x, rot_x, rot_y); + assert(mid_x != cur_x || mid_y != cur_y); + assert(mid_x != new_x || mid_y != new_y); // Compute maximum distance between the given parabolic arc // and line segment that discretize it. CT dist = (new_y - cur_y) * (mid_x - cur_x) - (new_x - cur_x) * (mid_y - cur_y); - dist = dist * dist / ((new_y - cur_y) * (new_y - cur_y) + - (new_x - cur_x) * (new_x - cur_x)); + CT div = (new_y - cur_y) * (new_y - cur_y) + (new_x - cur_x) * (new_x - cur_x); + assert(div != 0); + dist = dist * dist / div; if (dist <= max_dist_transformed) { // Distance between parabola and line segment is less than max_dist. point_stack.pop(); @@ -236,50 +242,39 @@ namespace Voronoi { namespace Internal { inline void clip_infinite_edge(const Points &points, const std::vector &segments, const edge_type& edge, coordinate_type bbox_max_size, std::vector* clipped_edge) { + assert(edge.is_infinite()); + assert((edge.vertex0() == nullptr) != (edge.vertex1() == nullptr)); + const cell_type& cell1 = *edge.cell(); const cell_type& cell2 = *edge.twin()->cell(); - point_type origin, direction; // Infinite edges could not be created by two segment sites. + assert(cell1.contains_point() || cell2.contains_point()); if (! cell1.contains_point() && ! cell2.contains_point()) { printf("Error! clip_infinite_edge - infinite edge separates two segment cells\n"); return; } + point_type direction; if (cell1.contains_point() && cell2.contains_point()) { + assert(! edge.is_secondary()); point_type p1 = retrieve_point(points, segments, cell1); point_type p2 = retrieve_point(points, segments, cell2); - origin.x((p1.x() + p2.x()) * 0.5); - origin.y((p1.y() + p2.y()) * 0.5); + if (edge.vertex0() == nullptr) + std::swap(p1, p2); direction.x(p1.y() - p2.y()); direction.y(p2.x() - p1.x()); } else { - origin = cell1.contains_segment() ? retrieve_point(points, segments, cell2) : retrieve_point(points, segments, cell1); + assert(edge.is_secondary()); segment_type segment = cell1.contains_segment() ? segments[cell1.source_index()] : segments[cell2.source_index()]; - coordinate_type dx = high(segment).x() - low(segment).x(); - coordinate_type dy = high(segment).y() - low(segment).y(); - if ((low(segment) == origin) ^ cell1.contains_point()) { - direction.x(dy); - direction.y(-dx); - } else { - direction.x(-dy); - direction.y(dx); - } + direction.x(high(segment).y() - low(segment).y()); + direction.y(low(segment).x() - high(segment).x()); } coordinate_type koef = bbox_max_size / (std::max)(fabs(direction.x()), fabs(direction.y())); - if (edge.vertex0() == NULL) { - clipped_edge->push_back(point_type( - origin.x() - direction.x() * koef, - origin.y() - direction.y() * koef)); + if (edge.vertex0() == nullptr) { + clipped_edge->push_back(point_type(edge.vertex1()->x() + direction.x() * koef, edge.vertex1()->y() + direction.y() * koef)); + clipped_edge->push_back(point_type(edge.vertex1()->x(), edge.vertex1()->y())); } else { - clipped_edge->push_back( - point_type(edge.vertex0()->x(), edge.vertex0()->y())); - } - if (edge.vertex1() == NULL) { - clipped_edge->push_back(point_type( - origin.x() + direction.x() * koef, - origin.y() + direction.y() * koef)); - } else { - clipped_edge->push_back( - point_type(edge.vertex1()->x(), edge.vertex1()->y())); + clipped_edge->push_back(point_type(edge.vertex0()->x(), edge.vertex0()->y())); + clipped_edge->push_back(point_type(edge.vertex0()->x() + direction.x() * koef, edge.vertex0()->y() + direction.y() * koef)); } } @@ -307,11 +302,16 @@ static inline void dump_voronoi_to_svg( const Lines &helper_lines = Lines(), double scale = 0) { + const bool internalEdgesOnly = false; + BoundingBox bbox; bbox.merge(get_extents(points)); bbox.merge(get_extents(lines)); bbox.merge(get_extents(offset_curves)); bbox.merge(get_extents(helper_lines)); + for (boost::polygon::voronoi_diagram::const_vertex_iterator it = vd.vertices().begin(); it != vd.vertices().end(); ++it) + if (! internalEdgesOnly || it->color() != Voronoi::Internal::EXTERNAL_COLOR) + bbox.merge(Point(it->x(), it->y())); bbox.min -= (0.01 * bbox.size().cast()).cast(); bbox.max += (0.01 * bbox.size().cast()).cast(); @@ -321,15 +321,17 @@ static inline void dump_voronoi_to_svg( 0.01 * std::min(bbox.size().x(), bbox.size().y()); else - scale /= SCALING_FACTOR; + scale *= SCALING_FACTOR; const std::string inputSegmentPointColor = "lightseagreen"; - const coord_t inputSegmentPointRadius = coord_t(0.09 * scale); + const coord_t inputSegmentPointRadius = std::max(1, coord_t(0.09 * scale)); const std::string inputSegmentColor = "lightseagreen"; const coord_t inputSegmentLineWidth = coord_t(0.03 * scale); const std::string voronoiPointColor = "black"; - const coord_t voronoiPointRadius = coord_t(0.06 * scale); + const std::string voronoiPointColorOutside = "red"; + const std::string voronoiPointColorInside = "blue"; + const coord_t voronoiPointRadius = std::max(1, coord_t(0.06 * scale)); const std::string voronoiLineColorPrimary = "black"; const std::string voronoiLineColorSecondary = "green"; const std::string voronoiArcColor = "red"; @@ -341,7 +343,6 @@ static inline void dump_voronoi_to_svg( const std::string helperLineColor = "orange"; const coord_t helperLineWidth = coord_t(0.04 * scale); - const bool internalEdgesOnly = false; const bool primaryEdgesOnly = false; ::Slic3r::SVG svg(path, bbox); @@ -360,9 +361,11 @@ static inline void dump_voronoi_to_svg( Voronoi::Internal::point_type(double(it->b(0)), double(it->b(1))))); // Color exterior edges. - for (boost::polygon::voronoi_diagram::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) - if (!it->is_finite()) - Voronoi::Internal::color_exterior(&(*it)); + if (internalEdgesOnly) { + for (boost::polygon::voronoi_diagram::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) + if (!it->is_finite()) + Voronoi::Internal::color_exterior(&(*it)); + } // Draw the end points of the input polygon. for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) { @@ -376,8 +379,19 @@ static inline void dump_voronoi_to_svg( #if 1 // Draw voronoi vertices. for (boost::polygon::voronoi_diagram::const_vertex_iterator it = vd.vertices().begin(); it != vd.vertices().end(); ++it) - if (! internalEdgesOnly || it->color() != Voronoi::Internal::EXTERNAL_COLOR) - svg.draw(Point(coord_t(it->x()), coord_t(it->y())), voronoiPointColor, voronoiPointRadius); + if (! internalEdgesOnly || it->color() != Voronoi::Internal::EXTERNAL_COLOR) { + const std::string *color = nullptr; + switch (Voronoi::vertex_category(*it)) { + case Voronoi::VertexCategory::OnContour: color = &voronoiPointColor; break; + case Voronoi::VertexCategory::Outside: color = &voronoiPointColorOutside; break; + case Voronoi::VertexCategory::Inside: color = &voronoiPointColorInside; break; + default: color = &voronoiPointColor; // assert(false); + } + Point pt(coord_t(it->x()), coord_t(it->y())); + if (it->x() * pt.x() >= 0. && it->y() * pt.y() >= 0.) + // Conversion to coord_t is valid. + svg.draw(Point(coord_t(it->x()), coord_t(it->y())), *color, voronoiPointRadius); + } for (boost::polygon::voronoi_diagram::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) { if (primaryEdgesOnly && !it->is_primary()) @@ -401,8 +415,32 @@ static inline void dump_voronoi_to_svg( } else if (! it->is_primary()) color = voronoiLineColorSecondary; } - for (std::size_t i = 0; i + 1 < samples.size(); ++i) - svg.draw(Line(Point(coord_t(samples[i].x()), coord_t(samples[i].y())), Point(coord_t(samples[i+1].x()), coord_t(samples[i+1].y()))), color, voronoiLineWidth); + for (std::size_t i = 0; i + 1 < samples.size(); ++ i) { + Vec2d a(samples[i].x(), samples[i].y()); + Vec2d b(samples[i+1].x(), samples[i+1].y()); + // Convert to coord_t. + Point ia = a.cast(); + Point ib = b.cast(); + // Is the conversion possible? Do the resulting points fit into int32_t? + auto in_range = [](const Point &ip, const Vec2d &p) { return p.x() * ip.x() >= 0. && p.y() * ip.y() >= 0.; }; + bool a_in_range = in_range(ia, a); + bool b_in_range = in_range(ib, b); + if (! a_in_range || ! b_in_range) { + if (! a_in_range && ! b_in_range) + // None fits, ignore. + continue; + // One fit, the other does not. Try to clip. + Vec2d v = b - a; + v.normalize(); + v *= bbox.size().cast().norm(); + auto p = a_in_range ? Vec2d(a + v) : Vec2d(b - v); + Point ip = p.cast(); + if (! in_range(ip, p)) + continue; + (a_in_range ? ib : ia) = ip; + } + svg.draw(Line(ia, ib), color, voronoiLineWidth); + } } #endif diff --git a/src/libslic3r/libslic3r.h b/src/libslic3r/libslic3r.h index 2ef258a4c..efdecbac8 100644 --- a/src/libslic3r/libslic3r.h +++ b/src/libslic3r/libslic3r.h @@ -228,6 +228,16 @@ ForwardIt binary_find_by_predicate(ForwardIt first, ForwardIt last, LowerThanKey return first != last && equal_to_key(*first) ? first : last; } +template inline bool contains(const ContainerType &c, const ValueType &v) + { return std::find(c.begin(), c.end(), v) != c.end(); } +template inline bool contains(const std::initializer_list &il, const T &v) + { return std::find(il.begin(), il.end(), v) != il.end(); } + +template inline bool one_of(const ValueType &v, const ContainerType &c) + { return contains(c, v); } +template inline bool one_of(const T& v, const std::initializer_list& il) + { return contains(il, v); } + template static inline T sqr(T x) { diff --git a/tests/libslic3r/test_voronoi.cpp b/tests/libslic3r/test_voronoi.cpp index b847a890b..189dde7ce 100644 --- a/tests/libslic3r/test_voronoi.cpp +++ b/tests/libslic3r/test_voronoi.cpp @@ -38,6 +38,26 @@ TEST_CASE("Voronoi missing edges - points 12067", "[Voronoi]") { -5, 0 } }; +#if 0 + for (Point &p : pts) { + Vec2d q = p.cast(); + p.x() = scale_(p.x()); + p.y() = scale_(p.y()); + } +#endif + +#if 0 + // Try to rotate, maybe the issue is caused by incorrect handling of vertical or horizontal edges? + double a = 0.18764587962597876897475f; + double c = cos(a); + double s = sin(a); + for (Point &p : poly.points) { + Vec2d q = p.cast(); + p.x() = coord_t(q.x() * c + q.y() * s + 0.5); + p.y() = coord_t(- q.x() * s + q.y() * c + 0.5); + } +#endif + // Construction of the Voronoi Diagram. VD vd; construct_voronoi(pts.begin(), pts.end(), &vd); @@ -133,33 +153,40 @@ TEST_CASE("Voronoi missing edges - Alessandro gapfill 12707", "[Voronoi]") { { 41427551, 0}, { 42127548, 699996 } } }; - Lines lines = to_lines(Polygon { - { 0, 10000000}, - { 700000, 1}, // it has to be 1, higher number, zero or -1 work. - { 700000, 9000000}, - { 9100000, 9000000}, - { 9100000, 0}, - {10000000, 10000000} - }); + Polygon poly { + { 0, 10000000}, + { 700000, 1}, // it has to be 1, higher number, zero or -1 work. + { 700000, 9000000}, + { 9100000, 9000000}, + { 9100000, 0}, + {10000000, 10000000} + }; + +#if 1 + // Try to rotate, maybe the issue is caused by incorrect handling of vertical or horizontal edges? + double a = 0.18764587962597876897475f; + double c = cos(a); + double s = sin(a); + for (Point &p : poly.points) { + Vec2d q = p.cast(); + p.x() = coord_t(q.x() * c + q.y() * s + 0.5); + p.y() = coord_t(- q.x() * s + q.y() * c + 0.5); + } +#endif - Polygon poly; std::mt19937 gen; std::uniform_int_distribution dist(-100, 100); - for (size_t i = 0; i < lines.size(); ++ i) { - Line &l1 = lines[i]; - Line &l2 = lines[(i + 1) % lines.size()]; - REQUIRE(l1.b.x() == l2.a.x()); - REQUIRE(l1.b.y() == l2.a.y()); + for (Point &p : poly.points) { #if 0 // Wiggle the points a bit to find out whether this fixes the voronoi diagram for this particular polygon. - l1.b.x() = (l2.a.x() += dist(gen)); - l1.b.y() = (l2.a.y() += dist(gen)); + p.x() = (p.x() += dist(gen)); + p.y() = (p.y() += dist(gen)); #endif - poly.points.emplace_back(l1.a); } REQUIRE(intersecting_edges({ poly }).empty()); + Lines lines = to_lines(poly); VD vd; construct_voronoi(lines.begin(), lines.end(), &vd); @@ -169,6 +196,114 @@ TEST_CASE("Voronoi missing edges - Alessandro gapfill 12707", "[Voronoi]") #endif } +TEST_CASE("Voronoi weirdness", "[Voronoi]") +{ + Polygon poly2 { + { 0, 0 }, + { 70000000, 0 }, + { 70000000, 1300000 }, +// { 70000001, 14000000 } + { 70700000, 14000000 } + }; + + Polygon poly5 { + { 35058884, -25732145 }, + { 35058884, -19586070 }, + { 32753739, -20246796 }, + { 28756244, -21010725 }, + { 24657532, -21657939 }, + { 20836260, -21960233 }, + { 16115145, -22070742 }, + { 11850152, -21839761 }, + { 7646240, -21470177 }, + { 3607605, -20786940 }, + { 1280947, -20329742 }, + { -292823, -19963790 }, + { -3844469, -18809741 }, + { -7237277, -17593723 }, + { -10225900, -16143761 }, + { -13030266, -14643721 }, + { -15404294, -12977561 }, + { -17601713, -11280712 }, + { -19241930, -9435607 }, + { -20714420, -7583739 }, + { -21726144, -5664355 }, + { -22579294, -3741947 }, + { -22966684, -1786321 }, + { -23200322, 170140 }, + { -22966684, 2126602 }, + { -22579296, 4082227 }, + { -21726148, 6004637 }, + { -20714424, 7924020 }, + { -19241932, 9775888 }, + { -17601717, 11620994 }, + { -15404423, 13317749 }, + { -13030276, 14984003 }, + { -10225910, 16484042 }, + { -7237288, 17934005 }, + { -3844482, 19150025 }, + { -292841, 20304074 }, + { 1280949, 20670031 }, + { 3607587, 21127226 }, + { 7646218, 21810465 }, + { 11850128, 22180055 }, + { 16115122, 22411036 }, + { 20836263, 22300531 }, + { 24657513, 21998239 }, + { 28756227, 21351025 }, + { 32753725, 20587092 }, + { 35058893, 19926309 }, + { 35058893, 35000000 }, + { -31657232, 35000000 }, + { -31657202, -35000000 }, + { 35058881, -35000000 } + }; + + Polygon poly7 { + { 35058884, -25732145 }, + { 35058884, -19586070 }, + { -31657202, -35000000 }, + { 35058881, -35000000 } + }; + +// coord_t shift = 35058881; + coord_t shift_ok = 17000000; + coord_t shift = 35058881; + Polygon poly { + // <-4, 0>: bug + // -5: ok + // 1 - ok + { 0 + shift, -35000000 }, + { 0 + shift, -25732145 }, + { 0 + shift, -19586070 }, + { -66716086 + shift, -35000000 } + }; + + REQUIRE(intersecting_edges({ poly }).empty()); + REQUIRE(poly.area() > 0.); + +#if 0 + // Try to rotate, maybe the issue is caused by incorrect handling of vertical or horizontal edges? + double a = 0.18764587962597876897475f; + double c = cos(a); + double s = sin(a); + for (Point &p : poly.points) { + Vec2d q = p.cast(); + p.x() = coord_t(q.x() * c + q.y() * s + 0.5); + p.y() = coord_t(- q.x() * s + q.y() * c + 0.5); + } +#endif + + VD vd; + Lines lines = to_lines(poly); + construct_voronoi(lines.begin(), lines.end(), &vd); + +#ifdef VORONOI_DEBUG_OUT + dump_voronoi_to_svg(debug_out_path("voronoi-weirdness.svg").c_str(), + vd, Points(), lines); +#endif +} + // https://svn.boost.org/trac10/ticket/12903 // division by zero reported, but this issue is most likely a non-issue, as it produces an infinity for the interval of validity // of the floating point calculation, therefore forcing a recalculation with extended accuracy. @@ -1228,7 +1363,7 @@ TEST_CASE("Voronoi offset", "[VoronoiOffset]") for (const OffsetTest &ot : { OffsetTest { scale_(0.2), 1, 1 }, OffsetTest { scale_(0.4), 1, 1 }, - OffsetTest { scale_(0.5), 1, 1 }, + OffsetTest { scale_(0.5), 1, 2 }, OffsetTest { scale_(0.505), 1, 2 }, OffsetTest { scale_(0.51), 1, 2 }, OffsetTest { scale_(0.52), 1, 1 }, @@ -1237,21 +1372,21 @@ TEST_CASE("Voronoi offset", "[VoronoiOffset]") OffsetTest { scale_(0.55), 1, 0 } }) { - Polygons offsetted_polygons_out = voronoi_offset(vd, lines, ot.distance, scale_(0.005)); - REQUIRE(offsetted_polygons_out.size() == ot.num_outer); - +#if 0 + Polygons offsetted_polygons_out = Slic3r::Voronoi::offset(vd, lines, ot.distance, scale_(0.005)); #ifdef VORONOI_DEBUG_OUT dump_voronoi_to_svg(debug_out_path("voronoi-offset-out-%lf.svg", ot.distance).c_str(), vd, Points(), lines, offsetted_polygons_out); #endif + REQUIRE(offsetted_polygons_out.size() == ot.num_outer); +#endif - Polygons offsetted_polygons_in = voronoi_offset(vd, lines, - ot.distance, scale_(0.005)); - REQUIRE(offsetted_polygons_in.size() == ot.num_inner); - + Polygons offsetted_polygons_in = Slic3r::Voronoi::offset(vd, lines, - ot.distance, scale_(0.005)); #ifdef VORONOI_DEBUG_OUT dump_voronoi_to_svg(debug_out_path("voronoi-offset-in-%lf.svg", ot.distance).c_str(), vd, Points(), lines, offsetted_polygons_in); #endif + REQUIRE(offsetted_polygons_in.size() == ot.num_inner); } } @@ -1296,21 +1431,20 @@ TEST_CASE("Voronoi offset 2", "[VoronoiOffset]") OffsetTest { scale_(0.4), 2, 2 }, OffsetTest { scale_(0.45), 2, 2 }, OffsetTest { scale_(0.48), 2, 2 }, -//FIXME Exact intersections of an Offset curve with any Voronoi vertex are not handled correctly yet. -// OffsetTest { scale_(0.5), 2, 2 }, + OffsetTest { scale_(0.5), 2, 4 }, OffsetTest { scale_(0.505), 2, 4 }, OffsetTest { scale_(0.7), 2, 0 }, OffsetTest { scale_(0.8), 1, 0 } }) { - Polygons offsetted_polygons_out = voronoi_offset(vd, lines, ot.distance, scale_(0.005)); + Polygons offsetted_polygons_out = Slic3r::Voronoi::offset(vd, lines, ot.distance, scale_(0.005)); #ifdef VORONOI_DEBUG_OUT dump_voronoi_to_svg(debug_out_path("voronoi-offset2-out-%lf.svg", ot.distance).c_str(), vd, Points(), lines, offsetted_polygons_out); #endif REQUIRE(offsetted_polygons_out.size() == ot.num_outer); - Polygons offsetted_polygons_in = voronoi_offset(vd, lines, - ot.distance, scale_(0.005)); + Polygons offsetted_polygons_in = Slic3r::Voronoi::offset(vd, lines, - ot.distance, scale_(0.005)); #ifdef VORONOI_DEBUG_OUT dump_voronoi_to_svg(debug_out_path("voronoi-offset2-in-%lf.svg", ot.distance).c_str(), vd, Points(), lines, offsetted_polygons_in); @@ -1360,14 +1494,13 @@ TEST_CASE("Voronoi offset 3", "[VoronoiOffset]") VD vd; Lines lines = to_lines(poly); - construct_voronoi(lines.begin(), lines.end(), &vd); + construct_voronoi(lines.begin(), lines.end(), &vd); for (const OffsetTest &ot : { OffsetTest { scale_(0.2), 2, 2 }, OffsetTest { scale_(0.4), 2, 2 }, OffsetTest { scale_(0.49), 2, 2 }, -//FIXME this fails -// OffsetTest { scale_(0.5), 2, 2 }, + OffsetTest { scale_(0.5), 2, 2 }, OffsetTest { scale_(0.51), 2, 2 }, OffsetTest { scale_(0.56), 2, 2 }, OffsetTest { scale_(0.6), 2, 2 }, @@ -1375,23 +1508,417 @@ TEST_CASE("Voronoi offset 3", "[VoronoiOffset]") OffsetTest { scale_(0.8), 1, 6 }, OffsetTest { scale_(0.9), 1, 6 }, OffsetTest { scale_(0.99), 1, 6 }, -//FIXME this fails -// OffsetTest { scale_(1.0), 1, 6 }, + OffsetTest { scale_(1.0), 1, 0 }, OffsetTest { scale_(1.01), 1, 0 }, }) { - Polygons offsetted_polygons_out = voronoi_offset(vd, lines, ot.distance, scale_(0.005)); + Polygons offsetted_polygons_out = Slic3r::Voronoi::offset(vd, lines, ot.distance, scale_(0.005)); #ifdef VORONOI_DEBUG_OUT - dump_voronoi_to_svg(debug_out_path("voronoi-offset2-out-%lf.svg", ot.distance).c_str(), + dump_voronoi_to_svg(debug_out_path("voronoi-offset3-out-%lf.svg", ot.distance).c_str(), vd, Points(), lines, offsetted_polygons_out); #endif REQUIRE(offsetted_polygons_out.size() == ot.num_outer); - Polygons offsetted_polygons_in = voronoi_offset(vd, lines, - ot.distance, scale_(0.005)); + Polygons offsetted_polygons_in = Slic3r::Voronoi::offset(vd, lines, - ot.distance, scale_(0.005)); #ifdef VORONOI_DEBUG_OUT - dump_voronoi_to_svg(debug_out_path("voronoi-offset2-in-%lf.svg", ot.distance).c_str(), + dump_voronoi_to_svg(debug_out_path("voronoi-offset3-in-%lf.svg", ot.distance).c_str(), vd, Points(), lines, offsetted_polygons_in); #endif REQUIRE(offsetted_polygons_in.size() == ot.num_inner); } } + +TEST_CASE("Voronoi offset with edge collapse", "[VoronoiOffset4]") +{ + Polygons poly + { + // Outer contour + { + { 434890, -64754575 }, + { 541060, -64669076 }, + { 585647, -64583538 }, + { 576999, -64484233 }, + { 485666, -64191173 }, + { 1029323, -63646261 }, + { 1416925, -63666243 }, + { 1541532, -63643075 }, + { 1541535, -63643075 }, + { 1541535, -63643074 }, + { 1552612, -63641015 }, + { 1631609, -63568270 }, + { 1678393, -63443639 }, + { 1725519, -63219942 }, + { 1770552, -62903321 }, + { 1742544, -62755934 }, + { 1758152, -62588378 }, + { 1702289, -62474826 }, + { 1638627, -62242058 }, + { 1671281, -61938568 }, + { 1792676, -61478332 }, + { 1934894, -60760597 }, + { 2160333, -60064089 }, + { 2638183, -58972890 }, + { 2933095, -58478536 }, + { 3173586, -58159527 }, + { 3557779, -57556727 }, + { 3707088, -57158174 }, + { 3758053, -56907139 }, + { 3935932, -56278649 }, + { 4077331, -56033426 }, + { 4151386, -55768739 }, + { 4293796, -55561683 }, + { 4383089, -55387391 }, + { 4614128, -55201443 }, + { 5268755, -54333831 }, + { 5547634, -53797296 }, + { 5573864, -53639250 }, + { 5736812, -53304157 }, + { 6090973, -52066302 }, + { 6148630, -51666664 }, + { 6189511, -50974991 }, + { 6244112, -50774489 }, + { 6302489, -50761556 }, + { 6511179, -50579861 }, + { 6778128, -50398667 }, + { 6896879, -50350264 }, + { 7388259, -49913262 }, + { 7630584, -49602449 }, + { 7886846, -48987881 }, + { 7871333, -48318666 }, + { 7770349, -47841179 }, + { 7570223, -47234733 }, + { 7283115, -46705503 }, + { 6842948, -46056539 }, + { 6612732, -45760004 }, + { 6150430, -44991494 }, + { 5564326, -44168114 }, + { 5193032, -43807544 }, + { 4932097, -43489047 }, + { 3857385, -42548846 }, + { 3764745, -42081468 }, + { 3539887, -41422273 }, + { 3283048, -40803856 }, + { 3005925, -40367981 }, + { 3136185, -39834533 }, + { 3333757, -38499972 }, + { 3360660, -38183895 }, + { 3353375, -38036005 }, + { 3398808, -37582504 }, + { 3604229, -37221781 }, + { 3698079, -36962934 }, + { 4000449, -36007804 }, + { 4256811, -35016707 }, + { 4405951, -34581428 }, + { 4364534, -34178439 }, + { 4124603, -33432250 }, + { 3806159, -32765167 }, + { 3359088, -32109020 }, + { 2880790, -31317192 }, + { 1548334, -30402139 }, + { -7, -30131023 }, + { -1548347, -30402139 }, + { -2880803, -31317189 }, + { -3359100, -32109018 }, + { -3806173, -32765166 }, + { -4124617, -33432248 }, + { -4364548, -34178436 }, + { -4405966, -34581423 }, + { -4256825, -35016707 }, + { -4000461, -36007800 }, + { -3698093, -36962933 }, + { -3604243, -37221781 }, + { -3398823, -37582501 }, + { -3353390, -38036003 }, + { -3360675, -38183892 }, + { -3333772, -38499972 }, + { -3136200, -39834530 }, + { -3005940, -40367980 }, + { -3283062, -40803852 }, + { -3539902, -41422273 }, + { -3764761, -42081468 }, + { -3857402, -42548846 }, + { -4932112, -43489047 }, + { -5193047, -43807544 }, + { -5564341, -44168114 }, + { -6150446, -44991494 }, + { -6612748, -45760000 }, + { -6842965, -46056536 }, + { -7283130, -46705503 }, + { -7570238, -47234733 }, + { -7770365, -47841179 }, + { -7871349, -48318666 }, + { -7886860, -48987873 }, + { -7630599, -49602451 }, + { -7388277, -49913260 }, + { -6896896, -50350264 }, + { -6778145, -50398667 }, + { -6511196, -50579862 }, + { -6302504, -50761557 }, + { -6244127, -50774488 }, + { -6189527, -50974989 }, + { -6148648, -51666664 }, + { -6090990, -52066302 }, + { -5736829, -53304157 }, + { -5573882, -53639250 }, + { -5547651, -53797296 }, + { -5268772, -54333831 }, + { -4614145, -55201443 }, + { -4383106, -55387389 }, + { -4293814, -55561683 }, + { -4151404, -55768739 }, + { -4077350, -56033423 }, + { -3935952, -56278647 }, + { -3758072, -56907137 }, + { -3707107, -57158170 }, + { -3557796, -57556727 }, + { -3173605, -58159527 }, + { -2933113, -58478536 }, + { -2638201, -58972890 }, + { -2295003, -59738435 }, + { -2160353, -60064089 }, + { -1978487, -60596626 }, + { -1792695, -61478332 }, + { -1671298, -61938574 }, + { -1638646, -62242058 }, + { -1702306, -62474826 }, + { -1758168, -62588380 }, + { -1742563, -62755934 }, + { -1770570, -62903317 }, + { -1725537, -63219946 }, + { -1678412, -63443639 }, + { -1631627, -63568270 }, + { -1553479, -63640201 }, + { -1416944, -63666243 }, + { -998854, -63645714 }, + { -485685, -64191173 }, + { -577019, -64484225 }, + { -585667, -64583538 }, + { -541079, -64669076 }, + { -434910, -64754575 }, + { -294192, -64810609 }, + { 294172, -64810609 }, + }, + // Hole 1 + { + { -842246, -45167428 }, + { -1473641, -45154392 }, + { -2181728, -45100161 }, + { -2804308, -44985842 }, + { -3100514, -44879269 }, + { -3836807, -44802155 }, + { -4035816, -44718913 }, + { -4167175, -44583299 }, + { -4496705, -44467674 }, + { -4486674, -44382045 }, + { -4343949, -44070577 }, + { -3740991, -43494686 }, + { -2701372, -43313358 }, + { -1493599, -43312838 }, + { -8, -43352868 }, + { 1414575, -43313336 }, + { 2701358, -43313358 }, + { 3095462, -43374735 }, + { 3740975, -43494686 }, + { 4343934, -44070583 }, + { 4486657, -44382044 }, + { 4496688, -44467670 }, + { 4167159, -44583300 }, + { 4035800, -44718913 }, + { 3836792, -44802155 }, + { 3100498, -44879269 }, + { 2804292, -44985842 }, + { 2181712, -45100161 }, + { 1473625, -45154389 }, + { 842231, -45167428 }, + { -8, -45232686 }, + }, + // Hole 2 + { + { 1657455, -63016679 }, + { 1723196, -63056286 }, + { 1541535, -63643074 }, + { 1541532, -63643075 }, + { 1030020, -63645562 }, + } + }; + + + VD vd; + Lines lines = to_lines(poly); + construct_voronoi(lines.begin(), lines.end(), &vd); + + for (const OffsetTest &ot : { + OffsetTest { scale_(0.2), 2, 2 }, + OffsetTest { scale_(0.4), 2, 2 }, + OffsetTest { scale_(0.49), 2, 3 }, + OffsetTest { scale_(0.51), 2, 2 }, + OffsetTest { scale_(0.56), 2, 2 }, + OffsetTest { scale_(0.6), 2, 2 }, + OffsetTest { scale_(0.7), 2, 2 }, + OffsetTest { scale_(0.8), 2, 2 }, + OffsetTest { scale_(0.9), 2, 2 }, + OffsetTest { scale_(0.99), 1, 2 }, + OffsetTest { scale_(1.0), 1, 2 }, + OffsetTest { scale_(1.01), 1, 2 }, + }) { + + Polygons offsetted_polygons_out = Slic3r::Voronoi::offset(vd, lines, ot.distance, scale_(0.005)); +#ifdef VORONOI_DEBUG_OUT + dump_voronoi_to_svg(debug_out_path("voronoi-offset3-out-%lf.svg", ot.distance).c_str(), + vd, Points(), lines, offsetted_polygons_out); +#endif + REQUIRE(offsetted_polygons_out.size() == ot.num_outer); + + Polygons offsetted_polygons_in = Slic3r::Voronoi::offset(vd, lines, - ot.distance, scale_(0.005)); +#ifdef VORONOI_DEBUG_OUT + dump_voronoi_to_svg(debug_out_path("voronoi-offset3-in-%lf.svg", ot.distance).c_str(), + vd, Points(), lines, offsetted_polygons_in); +#endif + REQUIRE(offsetted_polygons_in.size() == ot.num_inner); + } +} + +// A sample extracted from file medallion_printable_fixed-teeth.stl from https://www.thingiverse.com/thing:1347129 +// This test for offset scale_(2.9) and bigger +// triggers assert(r < std::max(d0, d1) + EPSILON) in function first_circle_segment_intersection_parameter. +TEST_CASE("Voronoi offset 5", "[VoronoiOffset5]") +{ + Polygons poly = { + Polygon { + { 0 , -16077501 }, + { 3015222 , -16142836 }, + { 3072642 , -16138163 }, + { 3094279 , -16105662 }, + { 3110660 , -15140728 }, + { 3157438 , -14105326 }, + { 3338367 , -11081394 }, + { 3443412 , -8381621 }, + { 3472489 , -6084497 }, + { 3445790 , -5962924 }, + { 3061725 , -6003484 }, + { 3030326 , -6030622 }, + { 2989343 , -6270378 }, + { 2903752 , -7368176 }, + { 2856704 , -7740619 }, + { 2795743 , -7978809 }, + { 2729231 , -8098866 }, + { 2666509 , -8131138 }, + { 2614169 , -8112308 }, + { 2561157 , -8032014 }, + { 2488290 , -7479351 }, + { 2453360 , -6911556 }, + { 2456148 , -6463146 }, + { 2546029 , -4580396 }, + { 2688181 , -2593262 }, + { 2717617 , -1700519 }, + { 2682232 , -1128411 }, + { 2631029 , -832886 }, + { 2535941 , -504483 }, + { 2399115 , -199303 }, + { 2290997 , -171213 }, + { 611824 , -132865 }, + { 6419 , -375849 }, + { -611825 , -132865 }, + { -2377355 , -185241 }, + { -2420740 , -231171 }, + { -2535942 , -504484 }, + { -2631030 , -832887 }, + { -2684831 , -1150821 }, + { -2714840 , -1586454 }, + { -2688181 , -2593262 }, + { -2546030 , -4580396 }, + { -2456149 , -6463145 }, + { -2453361 , -6911557 }, + { -2488291 , -7479352 }, + { -2561159 , -8032018 }, + { -2614171 , -8112309 }, + { -2666509 , -8131138 }, + { -2729233 , -8098868 }, + { -2795744 , -7978809 }, + { -2856706 , -7740619 }, + { -2903752 , -7368176 }, + { -2989345 , -6270378 }, + { -3030327 , -6030622 }, + { -3061726 , -6003484 }, + { -3445790 , -5962924 }, + { -3472490 , -6084498 }, + { -3468804 , -7244095 }, + { -3399287 , -9714025 }, + { -3338368 , -11081395 }, + { -3141015 , -14446051 }, + { -3094280 , -16105662 }, + { -3072643 , -16138163 }, + { -3008836 , -16143225 } + } + }; + double area = std::accumulate(poly.begin(), poly.end(), 0., [](double a, auto &poly){ return a + poly.area(); }); + REQUIRE(area > 0.); + + VD vd; + Lines lines = to_lines(poly); + construct_voronoi(lines.begin(), lines.end(), &vd); + + for (const OffsetTest &ot : { + OffsetTest { scale_(2.8), 1, 1 }, + OffsetTest { scale_(2.9), 1, 1 }, + OffsetTest { scale_(3.0), 1, 1 }, + }) { + + Polygons offsetted_polygons_out = Slic3r::Voronoi::offset(vd, lines, ot.distance, scale_(0.005)); +#ifdef VORONOI_DEBUG_OUT + dump_voronoi_to_svg(debug_out_path("voronoi-offset5-out-%lf.svg", ot.distance).c_str(), + vd, Points(), lines, offsetted_polygons_out); +#endif + REQUIRE(offsetted_polygons_out.size() == ot.num_outer); + + Polygons offsetted_polygons_in = Slic3r::Voronoi::offset(vd, lines, - ot.distance, scale_(0.005)); +#ifdef VORONOI_DEBUG_OUT + dump_voronoi_to_svg(debug_out_path("voronoi-offset5-in-%lf.svg", ot.distance).c_str(), + vd, Points(), lines, offsetted_polygons_in); +#endif + REQUIRE(offsetted_polygons_in.size() == ot.num_inner); + } +} + +TEST_CASE("Voronoi skeleton", "[VoronoiSkeleton]") +{ + coord_t mm = coord_t(scale_(1.)); + Polygons poly = { + Polygon { + { 0, 0 }, + { 1, 0 }, + { 1, 1 }, + { 2, 1 }, + { 2, 0 }, + { 3, 0 }, + { 3, 2 }, + { 0, 2 } + }, + Polygon { + { 0, - 1 - 2 }, + { 3, - 1 - 2 }, + { 3, - 1 - 0 }, + { 2, - 1 - 0 }, + { 2, - 1 - 1 }, + { 1, - 1 - 1 }, + { 1, - 1 - 0 }, + { 0, - 1 - 0 } + }, + }; + for (Polygon &p : poly) + for (Point &pt : p.points) + pt *= mm; + + double area = std::accumulate(poly.begin(), poly.end(), 0., [](double a, auto &poly){ return a + poly.area(); }); + REQUIRE(area > 0.); + + VD vd; + Lines lines = to_lines(poly); + construct_voronoi(lines.begin(), lines.end(), &vd); + Slic3r::Voronoi::annotate_inside_outside(vd, lines); + Slic3r::Voronoi::annotate_inside_outside(vd, lines); + static constexpr double threshold_alpha = M_PI / 12.; // 30 degrees + std::vector skeleton_edges = Slic3r::Voronoi::skeleton_edges_rough(vd, lines, threshold_alpha); + + REQUIRE(! skeleton_edges.empty()); +}