From 1c95ceaeaa49d4815e4a60f3f94b94bcb397aef0 Mon Sep 17 00:00:00 2001 From: Vojtech Bubnik Date: Thu, 11 Jun 2020 16:09:51 +0200 Subject: [PATCH] Reworked algorithm for Voronoi Offset curve extraction. Now the algorithm is very different from the OpenVoronoi implementation and hopefully it is now correct (save numerical issues, which will be a big PITA). --- src/libslic3r/CMakeLists.txt | 1 + src/libslic3r/Line.cpp | 2 +- src/libslic3r/Line.hpp | 4 +- src/libslic3r/VoronoiOffset.cpp | 643 +++++++++++++++++++-------- src/libslic3r/VoronoiOffset.hpp | 13 +- src/libslic3r/VoronoiVisualUtils.hpp | 407 +++++++++++++++++ tests/libslic3r/test_voronoi.cpp | 421 +----------------- 7 files changed, 908 insertions(+), 583 deletions(-) create mode 100644 src/libslic3r/VoronoiVisualUtils.hpp diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 1a58bdbbd..881466b39 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -190,6 +190,7 @@ add_library(libslic3r STATIC MTUtils.hpp VoronoiOffset.cpp VoronoiOffset.hpp + VoronoiVisualUtils.hpp Zipper.hpp Zipper.cpp MinAreaBoundingBox.hpp diff --git a/src/libslic3r/Line.cpp b/src/libslic3r/Line.cpp index 429cde372..974f585dc 100644 --- a/src/libslic3r/Line.cpp +++ b/src/libslic3r/Line.cpp @@ -135,4 +135,4 @@ BoundingBox get_extents(const Lines &lines) return bbox; } -} +} // namespace Slic3r diff --git a/src/libslic3r/Line.hpp b/src/libslic3r/Line.hpp index f255bee93..caab809f5 100644 --- a/src/libslic3r/Line.hpp +++ b/src/libslic3r/Line.hpp @@ -103,7 +103,7 @@ public: Vec3d b; }; -extern BoundingBox get_extents(const Lines &lines); +BoundingBox get_extents(const Lines &lines); } // namespace Slic3r @@ -125,4 +125,4 @@ namespace boost { namespace polygon { } } // end Boost -#endif +#endif // slic3r_Line_hpp_ diff --git a/src/libslic3r/VoronoiOffset.cpp b/src/libslic3r/VoronoiOffset.cpp index cd96e3cdc..7c736f19d 100644 --- a/src/libslic3r/VoronoiOffset.cpp +++ b/src/libslic3r/VoronoiOffset.cpp @@ -1,11 +1,15 @@ -// Polygon offsetting code inspired by OpenVoronoi by Anders Wallin -// https://github.com/aewallin/openvoronoi -// This offsetter uses results of boost::polygon Voronoi. +// Polygon offsetting using Voronoi diagram prodiced by boost::polygon. #include "VoronoiOffset.hpp" #include +// #define VORONOI_DEBUG_OUT + +#ifdef VORONOI_DEBUG_OUT +#include +#endif + namespace Slic3r { using VD = Geometry::VoronoiDiagram; @@ -48,6 +52,93 @@ namespace detail { } } + struct Intersections + { + int count; + Vec2d pts[2]; + }; + + // Return maximum two points, that are at distance "d" from both points + Intersections point_point_equal_distance_points(const Point &pt1, const Point &pt2, const double d) + { + // input points + const auto cx = double(pt1.x()); + const auto cy = double(pt1.y()); + const auto qx = double(pt2.x()); + const auto qy = double(pt2.y()); + + // Calculating determinant. + auto x0 = 2. * qy; + auto cx2 = cx * cx; + auto cy2 = cy * cy; + auto x5 = 2 * cx * qx; + auto x6 = cy * x0; + auto qx2 = qx * qx; + auto qy2 = qy * qy; + auto x9 = qx2 + qy2; + auto x10 = cx2 + cy2 - x5 - x6 + x9; + auto x11 = - cx2 - cy2; + auto discr = x10 * (4. * d + x11 + x5 + x6 - qx2 - qy2); + if (discr < 0.) + // No intersection point found, the two circles are too far away. + return Intersections { 0, { Vec2d(), Vec2d() } }; + + // Some intersections are found. + int npoints = (discr > 0) ? 2 : 1; + auto x1 = 2. * cy - x0; + auto x2 = cx - qx; + auto x12 = 0.5 * x2 * sqrt(discr) / x10; + auto x13 = 0.5 * (cy + qy); + auto x14 = - x12 + x13; + auto x15 = x11 + x9; + auto x16 = 0.5 / x2; + auto x17 = x12 + x13; + return Intersections { npoints, { Vec2d(- x16 * (x1 * x14 + x15), x14), + Vec2d(- x16 * (x1 * x17 + x15), x17) } }; + } + + // Return maximum two points, that are at distance "d" from both the line and point. + Intersections line_point_equal_distance_points(const Line &line, const Point &pt, const double d) + { + assert(line.a != pt && line.b != pt); + // Calculating two points of distance "d" to a ray and a point. + // Point. + auto x0 = double(pt.x()); + auto y0 = double(pt.y()); + // Ray equation. Vector (a, b) is perpendicular to line. + auto a = double(line.a.y() - line.b.y()); + auto b = double(line.b.x() - line.a.x()); + // pt shall not lie on line. + assert(std::abs((x0 - line.a.x()) * a + (y0 - line.a.y()) * b) < SCALED_EPSILON); + // Orient line so that the vector (a, b) points towards pt. + if (a * (x0 - line.a.x()) + b * (y0 - line.a.y()) < 0.) + std::swap(x0, y0); + double c = - a * double(line.a.x()) - b * double(line.a.y()); + // Calculate the two points. + double a2 = a * a; + double b2 = b * b; + double a2b2 = a2 + b2; + double d2 = d * d; + double s = a2*d2 - a2*sqr(x0) - 2*a*b*x0*y0 - 2*a*c*x0 + 2*a*d*x0 + b2*d2 - b2*sqr(y0) - 2*b*c*y0 + 2*b*d*y0 - sqr(c) + 2*c*d - d2; + if (s < 0.) + // Distance of pt from line is bigger than 2 * d. + return Intersections { 0 }; + double u; + int cnt; + if (s == 0.) { + // Distance of pt from line is 2 * d. + cnt = 1; + u = 0.; + } else { + // Distance of pt from line is smaller than 2 * d. + cnt = 2; + u = a*sqrt(s)/a2b2; + } + double v = (-a2*y0 + a*b*x0 + b*c - b*d)/a2b2; + return Intersections { cnt, { Vec2d((b * ( u + v) - c + d) / a, - u - v), + Vec2d((b * (- u + v) - c + d) / a, u - v) } }; + } + Vec2d voronoi_edge_offset_point( const VD &vd, const Lines &lines, @@ -131,174 +222,384 @@ namespace detail { } }; -Polygons voronoi_offset(const VD &vd, const Lines &lines, double offset_distance, double discretization_error) +static Vec2d foot_pt(const Line &iline, const Point &ipt) { - // Distance of a VD vertex to the closest site (input polygon edge or vertex). - std::vector vertex_dist(vd.num_vertices(), std::numeric_limits::max()); + Vec2d pt = iline.a.cast(); + Vec2d dir = (iline.b - iline.a).cast(); + Vec2d v = ipt.cast() - pt; + double l2 = dir.squaredNorm(); + double t = (l2 == 0.) ? 0. : v.dot(dir) / l2; + return pt + dir * t; +} - // Minium distance of a VD edge to the closest site (input polygon edge or vertex). - // For a parabolic segment the distance may be smaller than the distance of the two end points. - std::vector edge_dist(vd.num_edges(), std::numeric_limits::max()); - - // Calculate minimum distance of input polygons to voronoi vertices and voronoi edges. - for (const VD::edge_type &edge : vd.edges()) { - 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()]; - double d0, d1, dmin; - if (v0 == nullptr || v1 == nullptr) { - assert(edge.is_infinite()); - if (cell->contains_point() && cell2->contains_point()) { - 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; - d0 = d1 = std::numeric_limits::max(); - if (v0 == nullptr && v1 == nullptr) { - dmin = (pt1.cast() - pt0.cast()).norm(); - } else { - Vec2d pt((pt0 + pt1).cast() * 0.5); - Vec2d dir(double(pt0.y() - pt1.y()), double(pt1.x() - pt0.x())); - Vec2d pt0d(pt0.x(), pt0.y()); - if (v0) { - Vec2d a(v0->x(), v0->y()); - d0 = (a - pt0d).norm(); - dmin = ((a - pt).dot(dir) < 0.) ? (a - pt0d).norm() : d0; - vertex_dist[v0 - &vd.vertices().front()] = d0; - } else { - Vec2d a(v1->x(), v1->y()); - d1 = (a - pt0d).norm(); - dmin = ((a - pt).dot(dir) < 0.) ? (a - pt0d).norm() : d1; - vertex_dist[v1 - &vd.vertices().front()] = d1; - } - } - } 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_linear()); - assert(edge.is_secondary()); +Polygons voronoi_offset( + const Geometry::VoronoiDiagram &vd, + const Lines &lines, + double offset_distance, + double discretization_error) +{ #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())); + // 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); + } + } +#endif // NDEBUG + + // 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_candidate(vd.num_edges(), 2); // unknown state + const VD::edge_type *front_edge = &vd.edges().front(); + for (const VD::edge_type &edge : vd.edges()) + if (edge.vertex1() == nullptr) { + // Infinite Voronoi edge separating two Point sites. + // Infinite edge is always outside and it has at least one valid vertex. + assert(edge.vertex0() != nullptr); + edge_candidate[&edge - front_edge] = outside; + // Opposite edge of an infinite edge is certainly not active. + edge_candidate[edge.twin() - front_edge] = 0; + } else if (edge.vertex1() != nullptr) { + // Finite edge. + 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; + } + if (line) { + const VD::vertex_type *v1 = edge.vertex1(); + assert(v1); + Vec2d l0(line->a.cast()); + Vec2d lv((line->b - line->a).cast()); + double side = cross2(lv, Vec2d(v1->x(), v1->y()) - l0); + edge_candidate[&edge - front_edge] = outside ? (side < 0.) : (side > 0.); + } + } + for (const VD::edge_type &edge : vd.edges()) + if (edge_candidate[&edge - front_edge] == 2) { + assert(edge.cell()->contains_point() && edge.twin()->cell()->contains_point()); + // Edge separating two point sources, not yet classified as inside / outside. + const VD::edge_type *e = &edge; + char state; + do { + state = edge_candidate[e - front_edge]; + if (state != 2) + break; + e = e->next(); + } while (e != &edge); + e = &edge; + do { + char &s = edge_candidate[e - front_edge]; + if (s == 2) { + assert(e->cell()->contains_point() && e->twin()->cell()->contains_point()); + assert(edge_candidate[e->twin() - front_edge] == 2); + s = state; + edge_candidate[e->twin() - front_edge] = state; + } + e = e->next(); + } while (e != &edge); + } + if (! outside) + offset_distance = - offset_distance; + +#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(); + } + { + Lines helper_lines; + for (const VD::edge_type &edge : vd.edges()) + if (edge_candidate[&edge - front_edge]) { + 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) { + // Unconstrained edge. Calculate a trimmed position. + assert(edge.is_linear()); + 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()]; + if (cell->contains_point() && cell2->contains_point()) { + 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; + // Direction vector of this unconstrained Voronoi edge. + Vec2d dir(double(pt0.y() - pt1.y()), double(pt1.x() - pt0.x())); + pt2 = Vec2d(v0->x(), v0->y()) + dir.normalized() * scale_(10.); + } 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); + // 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); + // dir is perpendicular to line. + Vec2d dir(line.a.y() - line.b.y(), line.b.x() - line.a.x()); + assert(dir.norm() > 0.); + if (((line.a == ipt) == cell->contains_point()) == (v0 == nullptr)) + dir = - dir; + pt2 = ipt.cast() + dir.normalized() * scale_(10.); + } } else { + pt2 = Vec2d(v1->x(), v1->y()); + // Clip the line by the bounding box, so that the coloring of the line will be visible. + Geometry::liang_barsky_line_clipping(pt1, pt2, BoundingBoxf(bbox.min.cast(), bbox.max.cast())); + } + helper_lines.emplace_back(Line(Point(pt1.cast()), Point(((pt1 + pt2) * 0.5).cast()))); + } + dump_voronoi_to_svg(debug_out_path("voronoi-offset-candidates1.svg").c_str(), vd, Points(), lines, Polygons(), helper_lines); + } +#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_candidate[&edge - front_edge] != 2); + size_t edge_idx = &edge - front_edge; + if (edge_candidate[edge_idx] == 1) { + // 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_candidate[edge_idx2] == 0); + if (cell->contains_point() && cell2->contains_point()) { 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())); - } - const Point &pt = 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); -#endif /* NDEBUG */ - if (v0) { - assert((Point(v0->x(), v0->y()) - pt).cast().norm() < SCALED_EPSILON); - d0 = dmin = 0.; - vertex_dist[v0 - &vd.vertices().front()] = d0; + 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(); + 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; + edge_candidate[edge_idx] = 3; + } else + edge_candidate[edge_idx] = 0; } else { - assert((Point(v1->x(), v1->y()) - pt).cast().norm() < SCALED_EPSILON); - d1 = dmin = 0.; - vertex_dist[v1 - &vd.vertices().front()] = d1; + // 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_linear()); + 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); + Vec2d pt = ipt.cast(); + Vec2d dir(line.a.y() - line.b.y(), line.b.x() - line.a.x()); + assert(dir.norm() > 0.); + double t = offset_distance / dir.norm(); + if (((line.a == ipt) == cell->contains_point()) == (v0 == nullptr)) + t = - t; + edge_offset_point[edge_idx] = pt + t * dir; + edge_candidate[edge_idx] = 3; } - } - } else { - // Finite edge has valid points at both sides. - 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.); - d0 = (dir * (vec0.dot(dir) / l2) - vec0).norm(); - d1 = (dir * (vec1.dot(dir) / l2) - vec1).norm(); - dmin = std::min(d0, d1); - } else { - assert(cell->contains_point() || cell2->contains_point()); - const Point &pt0 = cell->contains_point() ? + // The other edge of an unconstrained edge starting with null vertex shall never be intersected. + edge_candidate[edge_idx2] = 0; + } else if (edge.is_secondary()) { + 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); - // Project p0 to line segment . - Vec2d p0(v0->x(), v0->y()); - Vec2d p1(v1->x(), v1->y()); - Vec2d px(pt0.x(), pt0.y()); - Vec2d v = p1 - p0; - d0 = (p0 - px).norm(); - d1 = (p1 - px).norm(); - double t = v.dot(px - p0); - double l2 = v.squaredNorm(); - if (t > 0. && t < l2) { - // Foot point on the line segment. - Vec2d foot = p0 + (t / l2) * v; - dmin = (foot - px).norm(); - } else - dmin = std::min(d0, d1); - } - vertex_dist[v0 - &vd.vertices().front()] = d0; - vertex_dist[v1 - &vd.vertices().front()] = d1; + 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); + 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; + edge_candidate[edge_idx] = 3; + } else { + edge_candidate[edge_idx] = 0; + } + edge_candidate[edge_idx2] = 0; + } 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)); + edge_candidate[edge_idx] = 3; + edge_candidate[edge_idx2] = 0; + 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; + if (offset_distance2 <= dmax) { + if (offset_distance2 >= dmin) { + has_intersection = true; + } else { + double dmin_new; + if (point_vs_segment) { + Vec2d ft = foot_pt(cell->contains_segment() ? line0 : line1, pt0); + dmin_new = (ft - px).squaredNorm() * 0.25; + } else { + // point vs. point + const Point &pt1 = (cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b; + dmin_new = (pt1.cast() - px).squaredNorm() * 0.25; + } + assert(dmin_new < dmax + SCALED_EPSILON); + assert(dmin_new < dmin + SCALED_EPSILON); + dmin = dmin_new; + has_intersection = 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 (intersections.count == 2) { + // Now decide which points fall on this Voronoi edge. + // Tangential points (single intersection) are ignored. + 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; + if (intersections.count == 2) { + edge_candidate[edge_idx] = edge_candidate[edge_idx2] = 3; + edge_offset_point[edge_idx] = intersections.pts[0]; + edge_offset_point[edge_idx2] = intersections.pts[1]; + done = true; + } else if (intersections.count == 1) { + if (d1 > d0) { + std::swap(edge_idx, edge_idx2); + edge_candidate[edge_idx] = 3; + edge_candidate[edge_idx2] = 0; + edge_offset_point[edge_idx] = intersections.pts[0]; + } + done = true; + } + } + if (! done) + edge_candidate[edge_idx] = edge_candidate[edge_idx2] = 0; + } + } + } } - edge_dist[&edge - &vd.edges().front()] = dmin; - } + } - // Mark cells intersected by the offset curve. - std::vector seed_cells(vd.num_cells(), false); - for (const VD::cell_type &cell : vd.cells()) { - const VD::edge_type *first_edge = cell.incident_edge(); - const VD::edge_type *edge = first_edge; - do { - double dmin = edge_dist[edge - &vd.edges().front()]; - double dmax = std::numeric_limits::max(); - const VD::vertex_type *v0 = edge->vertex0(); - const VD::vertex_type *v1 = edge->vertex1(); - if (v0 != nullptr) - dmax = vertex_dist[v0 - &vd.vertices().front()]; - if (v1 != nullptr) - dmax = std::max(dmax, vertex_dist[v1 - &vd.vertices().front()]); - if (offset_distance >= dmin && offset_distance <= dmax) { - // This cell is being intersected by the offset curve. - seed_cells[&cell - &vd.cells().front()] = true; - break; - } - edge = edge->next(); - } while (edge != first_edge); - } - auto edge_dir = [&vd, &vertex_dist, &edge_dist, offset_distance](const VD::edge_type *edge) { - const VD::vertex_type *v0 = edge->vertex0(); - const VD::vertex_type *v1 = edge->vertex1(); - double d0 = v0 ? vertex_dist[v0 - &vd.vertices().front()] : std::numeric_limits::max(); - double d1 = v1 ? vertex_dist[v1 - &vd.vertices().front()] : std::numeric_limits::max(); - if (d0 < offset_distance && offset_distance < d1) - return true; - else if (d1 < offset_distance && offset_distance < d0) - return false; - else { - assert(false); - return false; - } - }; +#ifdef VORONOI_DEBUG_OUT + { + Lines helper_lines; + for (const VD::edge_type &edge : vd.edges()) + if (edge_candidate[&edge - front_edge] == 3) + helper_lines.emplace_back(Line(Point(edge.vertex0()->x(), edge.vertex0()->y()), Point(edge_offset_point[&edge - front_edge].cast()))); + dump_voronoi_to_svg(debug_out_path("voronoi-offset-candidates2.svg").c_str(), vd, Points(), lines, Polygons(), helper_lines); + } +#endif // VORONOI_DEBUG_OUT - /// \brief starting at e, find the next edge on the face that brackets t - /// - /// we can be in one of two modes. - /// if direction==false then we are looking for an edge where src_t < t < trg_t - /// if direction==true we are looning for an edge where trg_t < t < src_t - auto next_offset_edge = - [&vd, &vertex_dist, &edge_dist, offset_distance] - (const VD::edge_type *start_edge, bool direction) -> const VD::edge_type* { - const VD::edge_type *edge = start_edge; - do { - const VD::vertex_type *v0 = edge->vertex0(); - const VD::vertex_type *v1 = edge->vertex1(); - double d0 = v0 ? vertex_dist[v0 - &vd.vertices().front()] : std::numeric_limits::max(); - double d1 = v1 ? vertex_dist[v1 - &vd.vertices().front()] : std::numeric_limits::max(); - if (direction ? (d1 < offset_distance && offset_distance < d0) : (d0 < offset_distance && offset_distance < d1)) - return edge; - edge = edge->next(); - } while (edge != start_edge); + auto next_offset_edge = [&edge_candidate, 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_candidate[edge->twin() - front_edge] == 3) + return edge->twin(); assert(false); return nullptr; }; @@ -316,28 +617,20 @@ Polygons voronoi_offset(const VD &vd, const Lines &lines, double offset_distance Polygons out; double angle_step = 2. * acos((offset_distance - discretization_error) / offset_distance); double sin_threshold = sin(angle_step) + EPSILON; - for (size_t seed_cell_idx = 0; seed_cell_idx < vd.num_cells(); ++ seed_cell_idx) - if (seed_cells[seed_cell_idx]) { - seed_cells[seed_cell_idx] = false; - // Initial direction should not matter, an offset curve shall intersect a cell at least at two points - // (if it is not just touching the cell at a single vertex), and such two intersection points shall have - // opposite direction. - bool direction = false; - // the first edge on the start-face - const VD::cell_type &cell = vd.cells()[seed_cell_idx]; - const VD::edge_type *start_edge = next_offset_edge(cell.incident_edge(), direction); - assert(start_edge->cell() == &cell); + for (size_t seed_edge_idx = 0; seed_edge_idx < vd.num_edges(); ++ seed_edge_idx) + if (edge_candidate[seed_edge_idx] == 3) { + const VD::edge_type *start_edge = &vd.edges()[seed_edge_idx]; const VD::edge_type *edge = start_edge; Polygon poly; do { - direction = edge_dir(edge); // find the next edge - const VD::edge_type *next_edge = next_offset_edge(edge->next(), direction); + const VD::edge_type *next_edge = next_offset_edge(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(); - Vec2d p1 = detail::voronoi_edge_offset_point(vd, lines, vertex_dist, edge_dist, *edge, offset_distance); - Vec2d p2 = detail::voronoi_edge_offset_point(vd, lines, vertex_dist, edge_dist, *next_edge, offset_distance); + edge_candidate[next_edge - front_edge] = 0; + Vec2d p1 = edge_offset_point[edge - front_edge]; + Vec2d p2 = edge_offset_point[next_edge - front_edge]; #ifndef NDEBUG { double err = dist_to_site(*cell, p1) - offset_distance; @@ -380,9 +673,7 @@ Polygons voronoi_offset(const VD &vd, const Lines &lines, double offset_distance } } poly.points.emplace_back(Point(coord_t(p2.x()), coord_t(p2.y()))); - // although we may revisit current_face (if it is non-convex), it seems safe to mark it "done" here. - seed_cells[cell - &vd.cells().front()] = false; - edge = next_edge->twin(); + edge = next_edge; } while (edge != start_edge); out.emplace_back(std::move(poly)); } diff --git a/src/libslic3r/VoronoiOffset.hpp b/src/libslic3r/VoronoiOffset.hpp index 9f5485c00..a21b44f93 100644 --- a/src/libslic3r/VoronoiOffset.hpp +++ b/src/libslic3r/VoronoiOffset.hpp @@ -1,3 +1,5 @@ +// Polygon offsetting using Voronoi diagram prodiced by boost::polygon. + #ifndef slic3r_VoronoiOffset_hpp_ #define slic3r_VoronoiOffset_hpp_ @@ -7,7 +9,16 @@ namespace Slic3r { -Polygons voronoi_offset(const Geometry::VoronoiDiagram &vd, const Lines &lines, 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); } // namespace Slic3r diff --git a/src/libslic3r/VoronoiVisualUtils.hpp b/src/libslic3r/VoronoiVisualUtils.hpp new file mode 100644 index 000000000..186bfb7ac --- /dev/null +++ b/src/libslic3r/VoronoiVisualUtils.hpp @@ -0,0 +1,407 @@ +#include + +#include +#include +#include +#include + +namespace boost { namespace polygon { + +// The following code for the visualization of the boost Voronoi diagram is based on: +// +// Boost.Polygon library voronoi_graphic_utils.hpp header file +// Copyright Andrii Sydorchuk 2010-2012. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +template +class voronoi_visual_utils { + public: + // Discretize parabolic Voronoi edge. + // Parabolic Voronoi edges are always formed by one point and one segment + // from the initial input set. + // + // Args: + // point: input point. + // segment: input segment. + // max_dist: maximum discretization distance. + // discretization: point discretization of the given Voronoi edge. + // + // Template arguments: + // InCT: coordinate type of the input geometries (usually integer). + // Point: point type, should model point concept. + // Segment: segment type, should model segment concept. + // + // Important: + // discretization should contain both edge endpoints initially. + template class Point, + template class Segment> + static + typename enable_if< + typename gtl_and< + typename gtl_if< + typename is_point_concept< + typename geometry_concept< Point >::type + >::type + >::type, + typename gtl_if< + typename is_segment_concept< + typename geometry_concept< Segment >::type + >::type + >::type + >::type, + void + >::type discretize( + const Point& point, + const Segment& segment, + const CT max_dist, + std::vector< Point >* discretization) { + // Apply the linear transformation to move start point of the segment to + // the point with coordinates (0, 0) and the direction of the segment to + // coincide the positive direction of the x-axis. + CT segm_vec_x = cast(x(high(segment))) - cast(x(low(segment))); + CT segm_vec_y = cast(y(high(segment))) - cast(y(low(segment))); + CT sqr_segment_length = segm_vec_x * segm_vec_x + segm_vec_y * segm_vec_y; + + // Compute x-coordinates of the endpoints of the edge + // in the transformed space. + CT projection_start = sqr_segment_length * + get_point_projection((*discretization)[0], segment); + CT projection_end = sqr_segment_length * + get_point_projection((*discretization)[1], segment); + + // Compute parabola parameters in the transformed space. + // Parabola has next representation: + // f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y). + CT point_vec_x = cast(x(point)) - cast(x(low(segment))); + CT point_vec_y = cast(y(point)) - cast(y(low(segment))); + CT rot_x = segm_vec_x * point_vec_x + segm_vec_y * point_vec_y; + CT rot_y = segm_vec_x * point_vec_y - segm_vec_y * point_vec_x; + + // Save the last point. + Point last_point = (*discretization)[1]; + discretization->pop_back(); + + // Use stack to avoid recursion. + std::stack point_stack; + point_stack.push(projection_end); + CT cur_x = projection_start; + CT cur_y = parabola_y(cur_x, rot_x, rot_y); + + // Adjust max_dist parameter in the transformed space. + const CT max_dist_transformed = max_dist * max_dist * sqr_segment_length; + while (!point_stack.empty()) { + CT new_x = point_stack.top(); + CT new_y = parabola_y(new_x, rot_x, rot_y); + + // Compute coordinates of the point of the parabola that is + // 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); + + // 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)); + if (dist <= max_dist_transformed) { + // Distance between parabola and line segment is less than max_dist. + point_stack.pop(); + CT inter_x = (segm_vec_x * new_x - segm_vec_y * new_y) / + sqr_segment_length + cast(x(low(segment))); + CT inter_y = (segm_vec_x * new_y + segm_vec_y * new_x) / + sqr_segment_length + cast(y(low(segment))); + discretization->push_back(Point(inter_x, inter_y)); + cur_x = new_x; + cur_y = new_y; + } else { + point_stack.push(mid_x); + } + } + + // Update last point. + discretization->back() = last_point; + } + + private: + // Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b). + static CT parabola_y(CT x, CT a, CT b) { + return ((x - a) * (x - a) + b * b) / (b + b); + } + + // Get normalized length of the distance between: + // 1) point projection onto the segment + // 2) start point of the segment + // Return this length divided by the segment length. This is made to avoid + // sqrt computation during transformation from the initial space to the + // transformed one and vice versa. The assumption is made that projection of + // the point lies between the start-point and endpoint of the segment. + template class Point, + template class Segment> + static + typename enable_if< + typename gtl_and< + typename gtl_if< + typename is_point_concept< + typename geometry_concept< Point >::type + >::type + >::type, + typename gtl_if< + typename is_segment_concept< + typename geometry_concept< Segment >::type + >::type + >::type + >::type, + CT + >::type get_point_projection( + const Point& point, const Segment& segment) { + CT segment_vec_x = cast(x(high(segment))) - cast(x(low(segment))); + CT segment_vec_y = cast(y(high(segment))) - cast(y(low(segment))); + CT point_vec_x = x(point) - cast(x(low(segment))); + CT point_vec_y = y(point) - cast(y(low(segment))); + CT sqr_segment_length = + segment_vec_x * segment_vec_x + segment_vec_y * segment_vec_y; + CT vec_dot = segment_vec_x * point_vec_x + segment_vec_y * point_vec_y; + return vec_dot / sqr_segment_length; + } + + template + static CT cast(const InCT& value) { + return static_cast(value); + } +}; + +} } // namespace boost::polygon + + +namespace Slic3r +{ + +// The following code for the visualization of the boost Voronoi diagram is based on: +// +// Boost.Polygon library voronoi_visualizer.cpp file +// Copyright Andrii Sydorchuk 2010-2012. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +namespace Voronoi { namespace Internal { + + using VD = Geometry::VoronoiDiagram; + typedef double coordinate_type; + typedef boost::polygon::point_data point_type; + typedef boost::polygon::segment_data segment_type; + typedef boost::polygon::rectangle_data rect_type; + typedef VD::cell_type cell_type; + typedef VD::cell_type::source_index_type source_index_type; + typedef VD::cell_type::source_category_type source_category_type; + typedef VD::edge_type edge_type; + typedef VD::cell_container_type cell_container_type; + typedef VD::cell_container_type vertex_container_type; + typedef VD::edge_container_type edge_container_type; + typedef VD::const_cell_iterator const_cell_iterator; + typedef VD::const_vertex_iterator const_vertex_iterator; + typedef VD::const_edge_iterator const_edge_iterator; + + static const std::size_t EXTERNAL_COLOR = 1; + + inline void color_exterior(const VD::edge_type* edge) + { + if (edge->color() == EXTERNAL_COLOR) + return; + edge->color(EXTERNAL_COLOR); + edge->twin()->color(EXTERNAL_COLOR); + const VD::vertex_type* v = edge->vertex1(); + if (v == NULL || !edge->is_primary()) + return; + v->color(EXTERNAL_COLOR); + const VD::edge_type* e = v->incident_edge(); + do { + color_exterior(e); + e = e->rot_next(); + } while (e != v->incident_edge()); + } + + inline point_type retrieve_point(const Points &points, const std::vector &segments, const cell_type& cell) + { + assert(cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT || cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_END_POINT || + cell.source_category() == boost::polygon::SOURCE_CATEGORY_SINGLE_POINT); + return cell.source_category() == boost::polygon::SOURCE_CATEGORY_SINGLE_POINT ? + Voronoi::Internal::point_type(double(points[cell.source_index()].x()), double(points[cell.source_index()].y())) : + (cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? + low(segments[cell.source_index()]) : high(segments[cell.source_index()]); + } + + 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) + { + 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. + if (! cell1.contains_point() && ! cell2.contains_point()) { + printf("Error! clip_infinite_edge - infinite edge separates two segment cells\n"); + return; + } + if (cell1.contains_point() && cell2.contains_point()) { + 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); + 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); + 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); + } + } + 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)); + } 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())); + } + } + + inline void sample_curved_edge(const Points &points, const std::vector &segments, const edge_type& edge, std::vector &sampled_edge, coordinate_type max_dist) + { + point_type point = edge.cell()->contains_point() ? + retrieve_point(points, segments, *edge.cell()) : + retrieve_point(points, segments, *edge.twin()->cell()); + segment_type segment = edge.cell()->contains_point() ? + segments[edge.twin()->cell()->source_index()] : + segments[edge.cell()->source_index()]; + ::boost::polygon::voronoi_visual_utils::discretize(point, segment, max_dist, &sampled_edge); + } + +} /* namespace Internal */ } // namespace Voronoi + +BoundingBox get_extents(const Lines &lines); + +static inline void dump_voronoi_to_svg( + const char *path, + const Geometry::VoronoiDiagram &vd, + const Points &points, + const Lines &lines, + const Polygons &offset_curves = Polygons(), + const Lines &helper_lines = Lines(), + const double scale = 0.7) // 0.2? +{ + const std::string inputSegmentPointColor = "lightseagreen"; + const coord_t inputSegmentPointRadius = coord_t(0.09 * scale / SCALING_FACTOR); + const std::string inputSegmentColor = "lightseagreen"; + const coord_t inputSegmentLineWidth = coord_t(0.03 * scale / SCALING_FACTOR); + + const std::string voronoiPointColor = "black"; + const coord_t voronoiPointRadius = coord_t(0.06 * scale / SCALING_FACTOR); + const std::string voronoiLineColorPrimary = "black"; + const std::string voronoiLineColorSecondary = "green"; + const std::string voronoiArcColor = "red"; + const coord_t voronoiLineWidth = coord_t(0.02 * scale / SCALING_FACTOR); + + const std::string offsetCurveColor = "magenta"; + const coord_t offsetCurveLineWidth = coord_t(0.09 * scale / SCALING_FACTOR); + + const std::string helperLineColor = "orange"; + const coord_t helperLineWidth = coord_t(0.09 * scale / SCALING_FACTOR); + + const bool internalEdgesOnly = false; + const bool primaryEdgesOnly = false; + + BoundingBox bbox; + bbox.merge(get_extents(points)); + bbox.merge(get_extents(lines)); + bbox.merge(get_extents(offset_curves)); + bbox.min -= (0.01 * bbox.size().cast()).cast(); + bbox.max += (0.01 * bbox.size().cast()).cast(); + + ::Slic3r::SVG svg(path, bbox); + +// bbox.scale(1.2); + // For clipping of half-lines to some reasonable value. + // The line will then be clipped by the SVG viewer anyway. + const double bbox_dim_max = double(std::max(bbox.size().x(), bbox.size().y())); + // For the discretization of the Voronoi parabolic segments. + const double discretization_step = 0.05 * bbox_dim_max; + + // Make a copy of the input segments with the double type. + std::vector segments; + for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++ it) + segments.push_back(Voronoi::Internal::segment_type( + Voronoi::Internal::point_type(double(it->a(0)), double(it->a(1))), + 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)); + + // Draw the end points of the input polygon. + for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) { + svg.draw(it->a, inputSegmentPointColor, inputSegmentPointRadius); + svg.draw(it->b, inputSegmentPointColor, inputSegmentPointRadius); + } + // Draw the input polygon. + for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) + svg.draw(Line(Point(coord_t(it->a(0)), coord_t(it->a(1))), Point(coord_t(it->b(0)), coord_t(it->b(1)))), inputSegmentColor, inputSegmentLineWidth); + +#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); + + for (boost::polygon::voronoi_diagram::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) { + if (primaryEdgesOnly && !it->is_primary()) + continue; + if (internalEdgesOnly && (it->color() == Voronoi::Internal::EXTERNAL_COLOR)) + continue; + std::vector samples; + std::string color = voronoiLineColorPrimary; + if (!it->is_finite()) { + Voronoi::Internal::clip_infinite_edge(points, segments, *it, bbox_dim_max, &samples); + if (! it->is_primary()) + color = voronoiLineColorSecondary; + } else { + // Store both points of the segment into samples. sample_curved_edge will split the initial line + // until the discretization_step is reached. + samples.push_back(Voronoi::Internal::point_type(it->vertex0()->x(), it->vertex0()->y())); + samples.push_back(Voronoi::Internal::point_type(it->vertex1()->x(), it->vertex1()->y())); + if (it->is_curved()) { + Voronoi::Internal::sample_curved_edge(points, segments, *it, samples, discretization_step); + color = voronoiArcColor; + } 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); + } +#endif + + svg.draw_outline(offset_curves, offsetCurveColor, offsetCurveLineWidth); + svg.draw(helper_lines, helperLineColor, helperLineWidth); + + svg.Close(); +} + +} // namespace Slic3r diff --git a/tests/libslic3r/test_voronoi.cpp b/tests/libslic3r/test_voronoi.cpp index ef05119ad..6d7211f37 100644 --- a/tests/libslic3r/test_voronoi.cpp +++ b/tests/libslic3r/test_voronoi.cpp @@ -1,16 +1,18 @@ #include #include -#include - #include #include #include #include + #include -#define BOOST_VORONOI_USE_GMP 1 -#include "boost/polygon/voronoi.hpp" +// #define VORONOI_DEBUG_OUT + +#ifdef VORONOI_DEBUG_OUT +#include +#endif using boost::polygon::voronoi_builder; using boost::polygon::voronoi_diagram; @@ -19,400 +21,6 @@ using namespace Slic3r; using VD = Geometry::VoronoiDiagram; -// #define VORONOI_DEBUG_OUT - -#ifdef VORONOI_DEBUG_OUT -#include -#endif - -#ifdef VORONOI_DEBUG_OUT -namespace boost { namespace polygon { - -// The following code for the visualization of the boost Voronoi diagram is based on: -// -// Boost.Polygon library voronoi_graphic_utils.hpp header file -// Copyright Andrii Sydorchuk 2010-2012. -// Distributed under the Boost Software License, Version 1.0. -// (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) -template -class voronoi_visual_utils { - public: - // Discretize parabolic Voronoi edge. - // Parabolic Voronoi edges are always formed by one point and one segment - // from the initial input set. - // - // Args: - // point: input point. - // segment: input segment. - // max_dist: maximum discretization distance. - // discretization: point discretization of the given Voronoi edge. - // - // Template arguments: - // InCT: coordinate type of the input geometries (usually integer). - // Point: point type, should model point concept. - // Segment: segment type, should model segment concept. - // - // Important: - // discretization should contain both edge endpoints initially. - template class Point, - template class Segment> - static - typename enable_if< - typename gtl_and< - typename gtl_if< - typename is_point_concept< - typename geometry_concept< Point >::type - >::type - >::type, - typename gtl_if< - typename is_segment_concept< - typename geometry_concept< Segment >::type - >::type - >::type - >::type, - void - >::type discretize( - const Point& point, - const Segment& segment, - const CT max_dist, - std::vector< Point >* discretization) { - // Apply the linear transformation to move start point of the segment to - // the point with coordinates (0, 0) and the direction of the segment to - // coincide the positive direction of the x-axis. - CT segm_vec_x = cast(x(high(segment))) - cast(x(low(segment))); - CT segm_vec_y = cast(y(high(segment))) - cast(y(low(segment))); - CT sqr_segment_length = segm_vec_x * segm_vec_x + segm_vec_y * segm_vec_y; - - // Compute x-coordinates of the endpoints of the edge - // in the transformed space. - CT projection_start = sqr_segment_length * - get_point_projection((*discretization)[0], segment); - CT projection_end = sqr_segment_length * - get_point_projection((*discretization)[1], segment); - - // Compute parabola parameters in the transformed space. - // Parabola has next representation: - // f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y). - CT point_vec_x = cast(x(point)) - cast(x(low(segment))); - CT point_vec_y = cast(y(point)) - cast(y(low(segment))); - CT rot_x = segm_vec_x * point_vec_x + segm_vec_y * point_vec_y; - CT rot_y = segm_vec_x * point_vec_y - segm_vec_y * point_vec_x; - - // Save the last point. - Point last_point = (*discretization)[1]; - discretization->pop_back(); - - // Use stack to avoid recursion. - std::stack point_stack; - point_stack.push(projection_end); - CT cur_x = projection_start; - CT cur_y = parabola_y(cur_x, rot_x, rot_y); - - // Adjust max_dist parameter in the transformed space. - const CT max_dist_transformed = max_dist * max_dist * sqr_segment_length; - while (!point_stack.empty()) { - CT new_x = point_stack.top(); - CT new_y = parabola_y(new_x, rot_x, rot_y); - - // Compute coordinates of the point of the parabola that is - // 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); - - // 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)); - if (dist <= max_dist_transformed) { - // Distance between parabola and line segment is less than max_dist. - point_stack.pop(); - CT inter_x = (segm_vec_x * new_x - segm_vec_y * new_y) / - sqr_segment_length + cast(x(low(segment))); - CT inter_y = (segm_vec_x * new_y + segm_vec_y * new_x) / - sqr_segment_length + cast(y(low(segment))); - discretization->push_back(Point(inter_x, inter_y)); - cur_x = new_x; - cur_y = new_y; - } else { - point_stack.push(mid_x); - } - } - - // Update last point. - discretization->back() = last_point; - } - - private: - // Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b). - static CT parabola_y(CT x, CT a, CT b) { - return ((x - a) * (x - a) + b * b) / (b + b); - } - - // Get normalized length of the distance between: - // 1) point projection onto the segment - // 2) start point of the segment - // Return this length divided by the segment length. This is made to avoid - // sqrt computation during transformation from the initial space to the - // transformed one and vice versa. The assumption is made that projection of - // the point lies between the start-point and endpoint of the segment. - template class Point, - template class Segment> - static - typename enable_if< - typename gtl_and< - typename gtl_if< - typename is_point_concept< - typename geometry_concept< Point >::type - >::type - >::type, - typename gtl_if< - typename is_segment_concept< - typename geometry_concept< Segment >::type - >::type - >::type - >::type, - CT - >::type get_point_projection( - const Point& point, const Segment& segment) { - CT segment_vec_x = cast(x(high(segment))) - cast(x(low(segment))); - CT segment_vec_y = cast(y(high(segment))) - cast(y(low(segment))); - CT point_vec_x = x(point) - cast(x(low(segment))); - CT point_vec_y = y(point) - cast(y(low(segment))); - CT sqr_segment_length = - segment_vec_x * segment_vec_x + segment_vec_y * segment_vec_y; - CT vec_dot = segment_vec_x * point_vec_x + segment_vec_y * point_vec_y; - return vec_dot / sqr_segment_length; - } - - template - static CT cast(const InCT& value) { - return static_cast(value); - } -}; - -} } // namespace boost::polygon - -// The following code for the visualization of the boost Voronoi diagram is based on: -// -// Boost.Polygon library voronoi_visualizer.cpp file -// Copyright Andrii Sydorchuk 2010-2012. -// Distributed under the Boost Software License, Version 1.0. -// (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) -namespace Voronoi { namespace Internal { - - typedef double coordinate_type; - typedef boost::polygon::point_data point_type; - typedef boost::polygon::segment_data segment_type; - typedef boost::polygon::rectangle_data rect_type; - typedef boost::polygon::voronoi_diagram VD; - typedef VD::cell_type cell_type; - typedef VD::cell_type::source_index_type source_index_type; - typedef VD::cell_type::source_category_type source_category_type; - typedef VD::edge_type edge_type; - typedef VD::cell_container_type cell_container_type; - typedef VD::cell_container_type vertex_container_type; - typedef VD::edge_container_type edge_container_type; - typedef VD::const_cell_iterator const_cell_iterator; - typedef VD::const_vertex_iterator const_vertex_iterator; - typedef VD::const_edge_iterator const_edge_iterator; - - static const std::size_t EXTERNAL_COLOR = 1; - - inline void color_exterior(const VD::edge_type* edge) - { - if (edge->color() == EXTERNAL_COLOR) - return; - edge->color(EXTERNAL_COLOR); - edge->twin()->color(EXTERNAL_COLOR); - const VD::vertex_type* v = edge->vertex1(); - if (v == NULL || !edge->is_primary()) - return; - v->color(EXTERNAL_COLOR); - const VD::edge_type* e = v->incident_edge(); - do { - color_exterior(e); - e = e->rot_next(); - } while (e != v->incident_edge()); - } - - inline point_type retrieve_point(const Points &points, const std::vector &segments, const cell_type& cell) - { - assert(cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT || cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_END_POINT || - cell.source_category() == boost::polygon::SOURCE_CATEGORY_SINGLE_POINT); - return cell.source_category() == boost::polygon::SOURCE_CATEGORY_SINGLE_POINT ? - Voronoi::Internal::point_type(double(points[cell.source_index()].x()), double(points[cell.source_index()].y())) : - (cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? - low(segments[cell.source_index()]) : high(segments[cell.source_index()]); - } - - 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) - { - 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. - if (! cell1.contains_point() && ! cell2.contains_point()) { - printf("Error! clip_infinite_edge - infinite edge separates two segment cells\n"); - return; - } - if (cell1.contains_point() && cell2.contains_point()) { - 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); - 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); - 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); - } - } - 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)); - } 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())); - } - } - - inline void sample_curved_edge(const Points &points, const std::vector &segments, const edge_type& edge, std::vector &sampled_edge, coordinate_type max_dist) - { - point_type point = edge.cell()->contains_point() ? - retrieve_point(points, segments, *edge.cell()) : - retrieve_point(points, segments, *edge.twin()->cell()); - segment_type segment = edge.cell()->contains_point() ? - segments[edge.twin()->cell()->source_index()] : - segments[edge.cell()->source_index()]; - ::boost::polygon::voronoi_visual_utils::discretize(point, segment, max_dist, &sampled_edge); - } - -} /* namespace Internal */ } // namespace Voronoi - -static inline void dump_voronoi_to_svg( - const char *path, - /* const */ VD &vd, - const Points &points, - const Lines &lines, - const Polygons &offset_curves = Polygons(), - const double scale = 0.7) // 0.2? -{ - const std::string inputSegmentPointColor = "lightseagreen"; - const coord_t inputSegmentPointRadius = coord_t(0.09 * scale / SCALING_FACTOR); - const std::string inputSegmentColor = "lightseagreen"; - const coord_t inputSegmentLineWidth = coord_t(0.03 * scale / SCALING_FACTOR); - - const std::string voronoiPointColor = "black"; - const coord_t voronoiPointRadius = coord_t(0.06 * scale / SCALING_FACTOR); - const std::string voronoiLineColorPrimary = "black"; - const std::string voronoiLineColorSecondary = "green"; - const std::string voronoiArcColor = "red"; - const coord_t voronoiLineWidth = coord_t(0.02 * scale / SCALING_FACTOR); - - const std::string offsetCurveColor = "magenta"; - const coord_t offsetCurveLineWidth = coord_t(0.09 * scale / SCALING_FACTOR); - - const bool internalEdgesOnly = false; - const bool primaryEdgesOnly = false; - - BoundingBox bbox; - bbox.merge(get_extents(points)); - bbox.merge(get_extents(lines)); - bbox.min -= (0.01 * bbox.size().cast()).cast(); - bbox.max += (0.01 * bbox.size().cast()).cast(); - - ::Slic3r::SVG svg(path, bbox); - -// bbox.scale(1.2); - // For clipping of half-lines to some reasonable value. - // The line will then be clipped by the SVG viewer anyway. - const double bbox_dim_max = double(std::max(bbox.size().x(), bbox.size().y())); - // For the discretization of the Voronoi parabolic segments. - const double discretization_step = 0.05 * bbox_dim_max; - - // Make a copy of the input segments with the double type. - std::vector segments; - for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++ it) - segments.push_back(Voronoi::Internal::segment_type( - Voronoi::Internal::point_type(double(it->a(0)), double(it->a(1))), - 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)); - - // Draw the end points of the input polygon. - for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) { - svg.draw(it->a, inputSegmentPointColor, inputSegmentPointRadius); - svg.draw(it->b, inputSegmentPointColor, inputSegmentPointRadius); - } - // Draw the input polygon. - for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) - svg.draw(Line(Point(coord_t(it->a(0)), coord_t(it->a(1))), Point(coord_t(it->b(0)), coord_t(it->b(1)))), inputSegmentColor, inputSegmentLineWidth); - -#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); - - for (boost::polygon::voronoi_diagram::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) { - if (primaryEdgesOnly && !it->is_primary()) - continue; - if (internalEdgesOnly && (it->color() == Voronoi::Internal::EXTERNAL_COLOR)) - continue; - std::vector samples; - std::string color = voronoiLineColorPrimary; - if (!it->is_finite()) { - Voronoi::Internal::clip_infinite_edge(points, segments, *it, bbox_dim_max, &samples); - if (! it->is_primary()) - color = voronoiLineColorSecondary; - } else { - // Store both points of the segment into samples. sample_curved_edge will split the initial line - // until the discretization_step is reached. - samples.push_back(Voronoi::Internal::point_type(it->vertex0()->x(), it->vertex0()->y())); - samples.push_back(Voronoi::Internal::point_type(it->vertex1()->x(), it->vertex1()->y())); - if (it->is_curved()) { - Voronoi::Internal::sample_curved_edge(points, segments, *it, samples, discretization_step); - color = voronoiArcColor; - } 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); - } -#endif - - svg.draw_outline(offset_curves, offsetCurveColor, offsetCurveLineWidth); - svg.Close(); -} -#endif - // https://svn.boost.org/trac10/ticket/12067 // This bug seems to be confirmed. // Vojtech supposes that there may be no Voronoi edges produced for @@ -1586,7 +1194,7 @@ TEST_CASE("Voronoi NaN coordinates 12139", "[Voronoi][!hide][!mayfail]") #ifdef VORONOI_DEBUG_OUT dump_voronoi_to_svg(debug_out_path("voronoi-NaNs.svg").c_str(), - vd, Points(), lines, Polygons(), 0.015); + vd, Points(), lines, Polygons(), Lines(), 0.015); #endif } @@ -1606,12 +1214,19 @@ TEST_CASE("Voronoi offset", "[VoronoiOffset]") Lines lines = to_lines(poly_with_hole); construct_voronoi(lines.begin(), lines.end(), &vd); - Polygons offsetted_polygons = voronoi_offset(vd, lines, scale_(0.2), scale_(0.005)); + Polygons offsetted_polygons_out = voronoi_offset(vd, lines, scale_(0.2), scale_(0.005)); + REQUIRE(offsetted_polygons_out.size() == 1); #ifdef VORONOI_DEBUG_OUT - dump_voronoi_to_svg(debug_out_path("voronoi-offset.svg").c_str(), - vd, Points(), lines, offsetted_polygons); + dump_voronoi_to_svg(debug_out_path("voronoi-offset-out.svg").c_str(), + vd, Points(), lines, offsetted_polygons_out); #endif - REQUIRE(offsetted_polygons.size() == 2); + Polygons offsetted_polygons_in = voronoi_offset(vd, lines, - scale_(0.2), scale_(0.005)); + REQUIRE(offsetted_polygons_in.size() == 1); + +#ifdef VORONOI_DEBUG_OUT + dump_voronoi_to_svg(debug_out_path("voronoi-offset-in.svg").c_str(), + vd, Points(), lines, offsetted_polygons_in); +#endif }