From 8ac60ccc7ac245adb87356bf9fbef5cae720c563 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Luk=C3=A1=C5=A1=20Hejl?= Date: Sat, 17 Sep 2022 13:21:10 +0200 Subject: [PATCH] Extend the algorithm for detecting non-planar Voronoi diagrams to include testing orientation between a line and a parabola and testing orientation between two parabolas. It fixed most of the issues reported in #8846. Co-authored-by: Vojtech Bubnik --- .../Arachne/SkeletalTrapezoidation.cpp | 8 +- src/libslic3r/Geometry/VoronoiUtilsCgal.cpp | 216 ++++++++++++++++-- src/libslic3r/Geometry/VoronoiUtilsCgal.hpp | 3 +- 3 files changed, 204 insertions(+), 23 deletions(-) diff --git a/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp b/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp index 04a1042d8c..34caee3111 100644 --- a/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp +++ b/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp @@ -619,8 +619,8 @@ void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys) // When any Voronoi vertex is missing, or the Voronoi diagram is not // planar, rotate the input polygon and try again. const bool has_missing_voronoi_vertex = detect_missing_voronoi_vertex(voronoi_diagram, segments); - // Detection of non-planar Voronoi diagram detects at least GH issues #8474, #8514 and #8446. - const bool is_voronoi_diagram_planar = Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram); + // Detection of non-planar Voronoi diagram detects at least GH issues #8474, #8514, #8446 and #8846. + const bool is_voronoi_diagram_planar = Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram, segments); const double fix_angle = PI / 6; std::unordered_map vertex_mapping; @@ -635,10 +635,10 @@ void SkeletalTrapezoidation::constructFromPolygons(const Polygons& polys) vertex_mapping = try_to_fix_degenerated_voronoi_diagram_by_rotation(voronoi_diagram, polys, polys_copy, segments, fix_angle); assert(!detect_missing_voronoi_vertex(voronoi_diagram, segments)); - assert(Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram)); + assert(Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram, segments)); if (detect_missing_voronoi_vertex(voronoi_diagram, segments)) BOOST_LOG_TRIVIAL(error) << "Detected missing Voronoi vertex even after the rotation of input."; - else if (!Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram)) + else if (!Geometry::VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(voronoi_diagram, segments)) BOOST_LOG_TRIVIAL(error) << "Detected non-planar Voronoi diagram even after the rotation of input."; } diff --git a/src/libslic3r/Geometry/VoronoiUtilsCgal.cpp b/src/libslic3r/Geometry/VoronoiUtilsCgal.cpp index 062a3b3979..c3348110b1 100644 --- a/src/libslic3r/Geometry/VoronoiUtilsCgal.cpp +++ b/src/libslic3r/Geometry/VoronoiUtilsCgal.cpp @@ -8,17 +8,135 @@ #include "VoronoiUtilsCgal.hpp" using VD = Slic3r::Geometry::VoronoiDiagram; +using namespace Slic3r::Arachne; namespace Slic3r::Geometry { -using CGAL_Point = CGAL::Exact_predicates_exact_constructions_kernel::Point_2; -using CGAL_Segment = CGAL::Arr_segment_traits_2::Curve_2; +// The tangent vector of the parabola is computed based on the Proof of the reflective property. +// https://en.wikipedia.org/wiki/Parabola#Proof_of_the_reflective_property +// https://math.stackexchange.com/q/2439647/2439663#comment5039739_2439663 +namespace impl { + using K = CGAL::Simple_cartesian; + using FK = CGAL::Simple_cartesian; + using EK = CGAL::Simple_cartesian; + using C2E = CGAL::Cartesian_converter; + using C2F = CGAL::Cartesian_converter; + class Epick : public CGAL::Filtered_kernel_adaptor::Type, Epick>, true> {}; -inline static CGAL_Point to_cgal_point(const VD::vertex_type &pt) { return {pt.x(), pt.y()}; } + template + inline typename K::Vector_2 calculate_parabolic_tangent_vector( + // Test point on the parabola, where the tangent will be calculated. + const typename K::Point_2 &p, + // Focus point of the parabola. + const typename K::Point_2 &f, + // Points of a directrix of the parabola. + const typename K::Point_2 &u, + const typename K::Point_2 &v, + // On which side of the parabolic segment endpoints the focus point is, which determines the orientation of the tangent. + const typename K::Orientation &tangent_orientation) + { + using RT = typename K::RT; + using Vector_2 = typename K::Vector_2; + + const Vector_2 directrix_vec = v - u; + const RT directrix_vec_sqr_length = CGAL::scalar_product(directrix_vec, directrix_vec); + Vector_2 focus_vec = (f - u) * directrix_vec_sqr_length - directrix_vec * CGAL::scalar_product(directrix_vec, p - u); + Vector_2 tangent_vec = focus_vec.perpendicular(tangent_orientation); + return tangent_vec; + } + + template struct ParabolicTangentToSegmentOrientationPredicate + { + using Point_2 = typename K::Point_2; + using Vector_2 = typename K::Vector_2; + using Orientation = typename K::Orientation; + using result_type = typename K::Orientation; + + result_type operator()( + // Test point on the parabola, where the tangent will be calculated. + const Point_2 &p, + // End of the linear segment (p, q), for which orientation towards the tangent to parabola will be evaluated. + const Point_2 &q, + // Focus point of the parabola. + const Point_2 &f, + // Points of a directrix of the parabola. + const Point_2 &u, + const Point_2 &v, + // On which side of the parabolic segment endpoints the focus point is, which determines the orientation of the tangent. + const Orientation &tangent_orientation) const + { + assert(tangent_orientation == CGAL::Orientation::LEFT_TURN || tangent_orientation == CGAL::Orientation::RIGHT_TURN); + + Vector_2 tangent_vec = calculate_parabolic_tangent_vector(p, f, u, v, tangent_orientation); + Vector_2 linear_vec = q - p; + + return CGAL::sign(tangent_vec.x() * linear_vec.y() - tangent_vec.y() * linear_vec.x()); + } + }; + + template struct ParabolicTangentToParabolicTangentOrientationPredicate + { + using Point_2 = typename K::Point_2; + using Vector_2 = typename K::Vector_2; + using Orientation = typename K::Orientation; + using result_type = typename K::Orientation; + + result_type operator()( + // Common point on both parabolas, where the tangent will be calculated. + const Point_2 &p, + // Focus point of the first parabola. + const Point_2 &f_0, + // Points of a directrix of the first parabola. + const Point_2 &u_0, + const Point_2 &v_0, + // On which side of the parabolic segment endpoints the focus point is, which determines the orientation of the tangent. + const Orientation &tangent_orientation_0, + // Focus point of the second parabola. + const Point_2 &f_1, + // Points of a directrix of the second parabola. + const Point_2 &u_1, + const Point_2 &v_1, + // On which side of the parabolic segment endpoints the focus point is, which determines the orientation of the tangent. + const Orientation &tangent_orientation_1) const + { + assert(tangent_orientation_0 == CGAL::Orientation::LEFT_TURN || tangent_orientation_0 == CGAL::Orientation::RIGHT_TURN); + assert(tangent_orientation_1 == CGAL::Orientation::LEFT_TURN || tangent_orientation_1 == CGAL::Orientation::RIGHT_TURN); + + Vector_2 tangent_vec_0 = calculate_parabolic_tangent_vector(p, f_0, u_0, v_0, tangent_orientation_0); + Vector_2 tangent_vec_1 = calculate_parabolic_tangent_vector(p, f_1, u_1, v_1, tangent_orientation_1); + + return CGAL::sign(tangent_vec_0.x() * tangent_vec_1.y() - tangent_vec_0.y() * tangent_vec_1.x()); + } + }; + + using ParabolicTangentToSegmentOrientationPredicateFiltered = CGAL::Filtered_predicate, ParabolicTangentToSegmentOrientationPredicate, C2E, C2F>; + using ParabolicTangentToParabolicTangentOrientationPredicateFiltered = CGAL::Filtered_predicate, ParabolicTangentToParabolicTangentOrientationPredicate, C2E, C2F>; +} // namespace impl + +using ParabolicTangentToSegmentOrientation = impl::ParabolicTangentToSegmentOrientationPredicateFiltered; +using ParabolicTangentToParabolicTangentOrientation = impl::ParabolicTangentToParabolicTangentOrientationPredicateFiltered; +using CGAL_Point = impl::K::Point_2; + +inline static CGAL_Point to_cgal_point(const VD::vertex_type *pt) { return {pt->x(), pt->y()}; } +inline static CGAL_Point to_cgal_point(const Point &pt) { return {pt.x(), pt.y()}; } +inline static CGAL_Point to_cgal_point(const Vec2d &pt) { return {pt.x(), pt.y()}; } + +inline static Linef make_linef(const VD::edge_type &edge) +{ + const VD::vertex_type *v0 = edge.vertex0(); + const VD::vertex_type *v1 = edge.vertex1(); + return {Vec2d(v0->x(), v0->y()), Vec2d(v1->x(), v1->y())}; +} + +inline static bool is_equal(const VD::vertex_type &first, const VD::vertex_type &second) { return first.x() == second.x() && first.y() == second.y(); } // FIXME Lukas H.: Also includes parabolic segments. bool VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(const VD &voronoi_diagram) { + using CGAL_Point = CGAL::Exact_predicates_exact_constructions_kernel::Point_2; + using CGAL_Segment = CGAL::Arr_segment_traits_2::Curve_2; + auto to_cgal_point = [](const VD::vertex_type &pt) -> CGAL_Point { return {pt.x(), pt.y()}; }; + assert(std::all_of(voronoi_diagram.edges().cbegin(), voronoi_diagram.edges().cend(), [](const VD::edge_type &edge) { return edge.color() == 0; })); @@ -30,7 +148,7 @@ bool VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(const VD &voronoi_ continue; if (edge.is_finite() && edge.is_linear() && edge.vertex0() != nullptr && edge.vertex1() != nullptr && - Arachne::VoronoiUtils::is_finite(*edge.vertex0()) && Arachne::VoronoiUtils::is_finite(*edge.vertex1())) { + VoronoiUtils::is_finite(*edge.vertex0()) && VoronoiUtils::is_finite(*edge.vertex1())) { segments.emplace_back(to_cgal_point(*edge.vertex0()), to_cgal_point(*edge.vertex1())); edge.color(1); assert(edge.twin() != nullptr); @@ -46,37 +164,101 @@ bool VoronoiUtilsCgal::is_voronoi_diagram_planar_intersection(const VD &voronoi_ return intersections_pt.empty(); } -static bool check_if_three_vectors_are_ccw(const CGAL_Point &common_pt, const CGAL_Point &pt_1, const CGAL_Point &pt_2, const CGAL_Point &test_pt) { - CGAL::Orientation orientation = CGAL::orientation(common_pt, pt_1, pt_2); +struct ParabolicSegment +{ + const Point focus; + const Line directrix; + // Two points on the parabola; + const Linef segment; + // Indicate if focus point is on the left side or right side relative to parabolic segment endpoints. + const CGAL::Orientation is_focus_on_left; +}; + +inline static ParabolicSegment get_parabolic_segment(const VD::edge_type &edge, const std::vector &segments) +{ + assert(edge.is_curved()); + + const VD::cell_type *left_cell = edge.cell(); + const VD::cell_type *right_cell = edge.twin()->cell(); + + const Point focus_pt = VoronoiUtils::getSourcePoint(*(left_cell->contains_point() ? left_cell : right_cell), segments); + const VoronoiUtils::Segment &directrix = VoronoiUtils::getSourceSegment(*(left_cell->contains_point() ? right_cell : left_cell), segments); + CGAL::Orientation focus_side = CGAL::opposite(CGAL::orientation(to_cgal_point(edge.vertex0()), to_cgal_point(edge.vertex1()), to_cgal_point(focus_pt))); + + assert(focus_side == CGAL::Orientation::LEFT_TURN || focus_side == CGAL::Orientation::RIGHT_TURN); + return {focus_pt, Line(directrix.from(), directrix.to()), make_linef(edge), focus_side}; +} + +inline static CGAL::Orientation orientation_of_two_edges(const VD::edge_type &edge_a, const VD::edge_type &edge_b, const std::vector &segments) { + assert(is_equal(*edge_a.vertex0(), *edge_b.vertex0())); + CGAL::Orientation orientation; + if (edge_a.is_linear() && edge_b.is_linear()) { + orientation = CGAL::orientation(to_cgal_point(edge_a.vertex0()), to_cgal_point(edge_a.vertex1()), to_cgal_point(edge_b.vertex1())); + } else if (edge_a.is_curved() && edge_b.is_curved()) { + const ParabolicSegment parabolic_a = get_parabolic_segment(edge_a, segments); + const ParabolicSegment parabolic_b = get_parabolic_segment(edge_b, segments); + orientation = ParabolicTangentToParabolicTangentOrientation{}(to_cgal_point(parabolic_a.segment.a), + to_cgal_point(parabolic_a.focus), + to_cgal_point(parabolic_a.directrix.a), + to_cgal_point(parabolic_a.directrix.b), + parabolic_a.is_focus_on_left, + to_cgal_point(parabolic_b.focus), + to_cgal_point(parabolic_b.directrix.a), + to_cgal_point(parabolic_b.directrix.b), + parabolic_b.is_focus_on_left); + return orientation; + } else { + assert(edge_a.is_curved() != edge_b.is_curved()); + + const VD::edge_type &linear_edge = edge_a.is_curved() ? edge_b : edge_a; + const VD::edge_type ¶bolic_edge = edge_a.is_curved() ? edge_a : edge_b; + const ParabolicSegment parabolic = get_parabolic_segment(parabolic_edge, segments); + orientation = ParabolicTangentToSegmentOrientation{}(to_cgal_point(parabolic.segment.a), to_cgal_point(linear_edge.vertex1()), + to_cgal_point(parabolic.focus), + to_cgal_point(parabolic.directrix.a), + to_cgal_point(parabolic.directrix.b), + parabolic.is_focus_on_left); + + if (edge_b.is_curved()) + orientation = CGAL::opposite(orientation); + } + + return orientation; +} + +static bool check_if_three_edges_are_ccw(const VD::edge_type &first, const VD::edge_type &second, const VD::edge_type &third, const std::vector &segments) +{ + assert(is_equal(*first.vertex0(), *second.vertex0()) && is_equal(*second.vertex0(), *third.vertex0())); + + CGAL::Orientation orientation = orientation_of_two_edges(first, second, segments); if (orientation == CGAL::Orientation::COLLINEAR) { // The first two edges are collinear, so the third edge must be on the right side on the first of them. - return CGAL::orientation(common_pt, pt_1, test_pt) == CGAL::Orientation::RIGHT_TURN; + return orientation_of_two_edges(first, third, segments) == CGAL::Orientation::RIGHT_TURN; } else if (orientation == CGAL::Orientation::LEFT_TURN) { // CCW oriented angle between vectors (common_pt, pt1) and (common_pt, pt2) is bellow PI. // So we need to check if test_pt isn't between them. - CGAL::Orientation orientation1 = CGAL::orientation(common_pt, pt_1, test_pt); - CGAL::Orientation orientation2 = CGAL::orientation(common_pt, pt_2, test_pt); + CGAL::Orientation orientation1 = orientation_of_two_edges(first, third, segments); + CGAL::Orientation orientation2 = orientation_of_two_edges(second, third, segments); return (orientation1 != CGAL::Orientation::LEFT_TURN || orientation2 != CGAL::Orientation::RIGHT_TURN); } else { assert(orientation == CGAL::Orientation::RIGHT_TURN); // CCW oriented angle between vectors (common_pt, pt1) and (common_pt, pt2) is upper PI. // So we need to check if test_pt is between them. - CGAL::Orientation orientation1 = CGAL::orientation(common_pt, pt_1, test_pt); - CGAL::Orientation orientation2 = CGAL::orientation(common_pt, pt_2, test_pt); + CGAL::Orientation orientation1 = orientation_of_two_edges(first, third, segments); + CGAL::Orientation orientation2 = orientation_of_two_edges(second, third, segments); return (orientation1 == CGAL::Orientation::RIGHT_TURN || orientation2 == CGAL::Orientation::LEFT_TURN); } } -bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VoronoiDiagram &voronoi_diagram) +bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VoronoiDiagram &voronoi_diagram, const std::vector &segments) { for (const VD::vertex_type &vertex : voronoi_diagram.vertices()) { std::vector edges; const VD::edge_type *edge = vertex.incident_edge(); do { - // FIXME Lukas H.: Also process parabolic segments. - if (edge->is_finite() && edge->is_linear() && edge->vertex0() != nullptr && edge->vertex1() != nullptr && - Arachne::VoronoiUtils::is_finite(*edge->vertex0()) && Arachne::VoronoiUtils::is_finite(*edge->vertex1())) + if (edge->is_finite() && edge->vertex0() != nullptr && edge->vertex1() != nullptr && + VoronoiUtils::is_finite(*edge->vertex0()) && VoronoiUtils::is_finite(*edge->vertex1())) edges.emplace_back(edge); edge = edge->rot_next(); @@ -89,8 +271,7 @@ bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VoronoiDiagram &vor const Geometry::VoronoiDiagram::edge_type *curr_edge = *edge_it; const Geometry::VoronoiDiagram::edge_type *next_edge = std::next(edge_it) == edges.end() ? edges.front() : *std::next(edge_it); - if (!check_if_three_vectors_are_ccw(to_cgal_point(*prev_edge->vertex0()), to_cgal_point(*prev_edge->vertex1()), - to_cgal_point(*curr_edge->vertex1()), to_cgal_point(*next_edge->vertex1()))) + if (!check_if_three_edges_are_ccw(*prev_edge, *curr_edge, *next_edge, segments)) return false; } } @@ -99,5 +280,4 @@ bool VoronoiUtilsCgal::is_voronoi_diagram_planar_angle(const VoronoiDiagram &vor return true; } - } // namespace Slic3r::Geometry \ No newline at end of file diff --git a/src/libslic3r/Geometry/VoronoiUtilsCgal.hpp b/src/libslic3r/Geometry/VoronoiUtilsCgal.hpp index 897891bd93..cad54615bf 100644 --- a/src/libslic3r/Geometry/VoronoiUtilsCgal.hpp +++ b/src/libslic3r/Geometry/VoronoiUtilsCgal.hpp @@ -2,6 +2,7 @@ #define slic3r_VoronoiUtilsCgal_hpp_ #include "Voronoi.hpp" +#include "../Arachne/utils/VoronoiUtils.hpp" namespace Slic3r::Geometry { class VoronoiDiagram; @@ -13,7 +14,7 @@ public: static bool is_voronoi_diagram_planar_intersection(const VoronoiDiagram &voronoi_diagram); // Check if the Voronoi diagram is planar using verification that all neighboring edges are ordered CCW for each vertex. - static bool is_voronoi_diagram_planar_angle(const VoronoiDiagram &voronoi_diagram); + static bool is_voronoi_diagram_planar_angle(const VoronoiDiagram &voronoi_diagram, const std::vector &segments); }; } // namespace Slic3r::Geometry