From afa72da9d18ac674bb062f2f8fbb9d15c4db1e17 Mon Sep 17 00:00:00 2001 From: bubnikv Date: Mon, 2 Dec 2019 15:01:52 +0100 Subject: [PATCH] Fix of infill connecting along perimeter lines, new 3-opt iterative improvement of infill path (currently disabled, it is extremely slow) --- src/libslic3r/Fill/Fill3DHoneycomb.cpp | 2 +- src/libslic3r/Fill/FillBase.cpp | 149 +++++++++++---------- src/libslic3r/Fill/FillBase.hpp | 4 +- src/libslic3r/Fill/FillGyroid.cpp | 3 +- src/libslic3r/ShortestPath.cpp | 173 ++++++++++++++++++++++++- 5 files changed, 259 insertions(+), 72 deletions(-) diff --git a/src/libslic3r/Fill/Fill3DHoneycomb.cpp b/src/libslic3r/Fill/Fill3DHoneycomb.cpp index 6c4b4d903..8aac6e49c 100644 --- a/src/libslic3r/Fill/Fill3DHoneycomb.cpp +++ b/src/libslic3r/Fill/Fill3DHoneycomb.cpp @@ -169,7 +169,7 @@ void Fill3DHoneycomb::_fill_surface_single( if (params.dont_connect) append(polylines_out, std::move(polylines_chained)); else - this->connect_infill(std::move(polylines_chained), expolygon, polylines_out, params); + this->connect_infill(std::move(polylines_chained), expolygon, polylines_out, this->spacing, params); } } diff --git a/src/libslic3r/Fill/FillBase.cpp b/src/libslic3r/Fill/FillBase.cpp index fa7dd2851..39679e3d6 100644 --- a/src/libslic3r/Fill/FillBase.cpp +++ b/src/libslic3r/Fill/FillBase.cpp @@ -610,16 +610,15 @@ static inline SegmentPoint clip_start_segment_and_point(const Points &polyline, // Initialized to "invalid". SegmentPoint out; if (polyline.size() >= 2) { - const double d2 = distance * distance; Vec2d pt_prev = polyline.front().cast(); for (int i = 1; i < polyline.size(); ++ i) { Vec2d pt = polyline[i].cast(); Vec2d v = pt - pt_prev; double l2 = v.squaredNorm(); - if (l2 > d2) { + if (l2 > distance * distance) { out.idx_segment = i; out.t = distance / sqrt(l2); - out.point = pt + out.t * v; + out.point = pt_prev + out.t * v; break; } distance -= sqrt(l2); @@ -636,16 +635,17 @@ static inline SegmentPoint clip_end_segment_and_point(const Points &polyline, do // Initialized to "invalid". SegmentPoint out; if (polyline.size() >= 2) { - const double d2 = distance * distance; Vec2d pt_next = polyline.back().cast(); for (int i = int(polyline.size()) - 2; i >= 0; -- i) { Vec2d pt = polyline[i].cast(); Vec2d v = pt - pt_next; double l2 = v.squaredNorm(); - if (l2 > d2) { + if (l2 > distance * distance) { out.idx_segment = i; out.t = distance / sqrt(l2); - out.point = pt + out.t * v; + out.point = pt_next + out.t * v; + // Store the parameter referenced to the starting point of a segment. + out.t = 1. - out.t; break; } distance -= sqrt(l2); @@ -655,21 +655,26 @@ static inline SegmentPoint clip_end_segment_and_point(const Points &polyline, do return out; } +// Optimized version with the precalculated v1 = p1b - p1a and l1_2 = v1.squaredNorm(). +// Assumption: l1_2 < EPSILON. +static inline double segment_point_distance_squared(const Vec2d &p1a, const Vec2d &p1b, const Vec2d &v1, const double l1_2, const Vec2d &p2) +{ + assert(l1_2 > EPSILON); + Vec2d v12 = p2 - p1a; + double t = v12.dot(v1); + return (t <= 0. ) ? v12.squaredNorm() : + (t >= l1_2) ? (p2 - p1a).squaredNorm() : + ((t / l1_2) * v1 - v12).squaredNorm(); +} + static inline double segment_point_distance_squared(const Vec2d &p1a, const Vec2d &p1b, const Vec2d &p2) { - const Vec2d v = p1b - p1a; - const Vec2d va = p2 - p1a; - const double l2 = v.squaredNorm(); + const Vec2d v = p1b - p1a; + const double l2 = v.squaredNorm(); if (l2 < EPSILON) // p1a == p1b - return va.squaredNorm(); - // Project p2 onto the (p1a, p1b) segment. - const double t = va.dot(v); - if (t < 0.) - return va.squaredNorm(); - else if (t > l2) - return (p2 - p1b).squaredNorm(); - return ((t / l2) * v - va).squaredNorm(); + return (p2 - p1a).squaredNorm(); + return segment_point_distance_squared(p1a, p1b, v, v.squaredNorm(), p2); } // Distance to the closest point of line. @@ -685,43 +690,11 @@ static inline double min_distance_of_segments(const Vec2d &p1a, const Vec2d &p1b double l2_2 = v2.squaredNorm(); if (l2_2 < EPSILON) // p2a == p2b: Return distance of p2a from the (p1a, p1b) segment. - return segment_point_distance_squared(p1a, p1b, p2a); - - // Project p2a, p2b onto the (p1a, p1b) segment. - auto project_p2a_p2b_onto_seg_p1a_p1b = [](const Vec2d& p1a, const Vec2d& p1b, const Vec2d& p2a, const Vec2d& p2b, const Vec2d& v1, const double l1_2) { - Vec2d v1a2a = p2a - p1a; - Vec2d v1a2b = p2b - p1a; - double t1 = v1a2a.dot(v1); - double t2 = v1a2b.dot(v1); - if (t1 <= 0.) { - if (t2 <= 0.) - // Both p2a and p2b are left of v1. - return (((t1 < t2) ? p2b : p2a) - p1a).squaredNorm(); - else if (t2 < l1_2) - // Project p2b onto the (p1a, p1b) segment. - return ((t2 / l1_2) * v1 - v1a2b).squaredNorm(); - } - else if (t1 >= l1_2) { - if (t2 >= l1_2) - // Both p2a and p2b are right of v1. - return (((t1 < t2) ? p2a : p2b) - p1b).squaredNorm(); - else if (t2 < l1_2) - // Project p2b onto the (p1a, p1b) segment. - return ((t2 / l1_2) * v1 - v1a2b).squaredNorm(); - } - else { - // Project p1b onto the (p1a, p1b) segment. - double dist_min = ((t2 / l1_2) * v1 - v1a2a).squaredNorm(); - if (t2 > 0. && t2 < l1_2) - dist_min = std::min(dist_min, ((t2 / l1_2) * v1 - v1a2b).squaredNorm()); - return dist_min; - } - return std::numeric_limits::max(); - }; + return segment_point_distance_squared(p1a, p1b, v1, l1_2, p2a); return std::min( - project_p2a_p2b_onto_seg_p1a_p1b(p1a, p1b, p2a, p2b, v1, l1_2), - project_p2a_p2b_onto_seg_p1a_p1b(p2a, p2b, p1a, p1b, v2, l2_2)); + std::min(segment_point_distance_squared(p1a, p1b, v1, l1_2, p2a), segment_point_distance_squared(p1a, p1b, v1, l1_2, p2b)), + std::min(segment_point_distance_squared(p2a, p2b, v2, l2_2, p1a), segment_point_distance_squared(p2a, p2b, v2, l2_2, p1b))); } // Mark the segments of split boundary as consumed if they are very close to some of the infill line. @@ -757,11 +730,26 @@ void mark_boundary_segments_touching_infill( const Vec2d seg_pt2 = segment.second.cast(); if (min_distance_of_segments(seg_pt1, seg_pt2, *this->pt1, *this->pt2) < this->dist2_max) { // Mark this boundary segment as touching the infill line. - ContourPointData&bdp = boundary_data[it_contour_and_segment->first][it_contour_and_segment->second]; + ContourPointData &bdp = boundary_data[it_contour_and_segment->first][it_contour_and_segment->second]; bdp.segment_consumed = true; // There is no need for checking seg_pt2 as it will be checked the next time. - if (segment_point_distance_squared(*this->pt1, *this->pt2, seg_pt1) < this->dist2_max) + bool point_touching = false; + if (segment_point_distance_squared(*this->pt1, *this->pt2, seg_pt1) < this->dist2_max) { + point_touching = true; bdp.point_consumed = true; + } +#if 0 + { + static size_t iRun = 0; + ExPolygon expoly(Polygon(*grid.contours().front())); + for (size_t i = 1; i < grid.contours().size(); ++i) + expoly.holes.emplace_back(Polygon(*grid.contours()[i])); + SVG svg(debug_out_path("%s-%d.svg", "FillBase-mark_boundary_segments_touching_infill", iRun ++).c_str(), get_extents(expoly)); + svg.draw(expoly, "green"); + svg.draw(Line(segment.first, segment.second), "red"); + svg.draw(Line(this->pt1->cast(), this->pt2->cast()), "magenta"); + } +#endif } } // Continue traversing the grid along the edge. @@ -780,6 +768,7 @@ void mark_boundary_segments_touching_infill( BoundingBoxf bboxf(boundary_bbox.min.cast(), boundary_bbox.max.cast()); bboxf.offset(- SCALED_EPSILON); + for (const Polyline &polyline : infill) { // Clip the infill polyline by the Eucledian distance along the polyline. SegmentPoint start_point = clip_start_segment_and_point(polyline.points, clip_distance); @@ -809,8 +798,20 @@ void mark_boundary_segments_touching_infill( visitor.init(pt1d, pt2d); grid.visit_cells_intersecting_thick_line(pt1, pt2, distance_colliding, visitor); #else - Vec2d pt1 = (point_idx == start_point.idx_segment) ? start_point.point : polyline.points[point_idx].cast(); - Vec2d pt2 = (point_idx == end_point .idx_segment) ? end_point .point : polyline.points[point_idx].cast(); + Vec2d pt1 = (point_idx == start_point.idx_segment) ? start_point.point : polyline.points[point_idx ].cast(); + Vec2d pt2 = (point_idx == end_point .idx_segment) ? end_point .point : polyline.points[point_idx + 1].cast(); +#if 0 + { + static size_t iRun = 0; + ExPolygon expoly(Polygon(*grid.contours().front())); + for (size_t i = 1; i < grid.contours().size(); ++i) + expoly.holes.emplace_back(Polygon(*grid.contours()[i])); + SVG svg(debug_out_path("%s-%d.svg", "FillBase-mark_boundary_segments_touching_infill0", iRun ++).c_str(), get_extents(expoly)); + svg.draw(expoly, "green"); + svg.draw(polyline, "blue"); + svg.draw(Line(pt1.cast(), pt2.cast()), "magenta", scale_(0.1)); + } +#endif visitor.init(pt1, pt2); // Simulate tracing of a thick line. This only works reliably if distance_colliding <= grid cell size. Vec2d v = (pt2 - pt1).normalized() * distance_colliding; @@ -829,7 +830,7 @@ void mark_boundary_segments_touching_infill( } } -void Fill::connect_infill(Polylines &&infill_ordered, const ExPolygon &boundary_src, Polylines &polylines_out, const FillParams ¶ms) +void Fill::connect_infill(Polylines &&infill_ordered, const ExPolygon &boundary_src, Polylines &polylines_out, const double spacing, const FillParams ¶ms) { assert(! infill_ordered.empty()); assert(! boundary_src.contour.points.empty()); @@ -905,16 +906,16 @@ void Fill::connect_infill(Polylines &&infill_ordered, const ExPolygon &boundary_ // Mark the points and segments of split boundary as consumed if they are very close to some of the infill line. { - //const double clip_distance = scale_(this->spacing); - const double clip_distance = 3. * scale_(this->spacing); - const double distance_colliding = scale_(this->spacing); + // @supermerill used 2. * scale_(spacing) + const double clip_distance = 3. * scale_(spacing); + const double distance_colliding = 1.1 * scale_(spacing); mark_boundary_segments_touching_infill(boundary, boundary_data, bbox, infill_ordered, clip_distance, distance_colliding); } // Connection from end of one infill line to the start of another infill line. - //const float length_max = scale_(this->spacing); -// const float length_max = scale_((2. / params.density) * this->spacing); - const float length_max = scale_((1000. / params.density) * this->spacing); + //const float length_max = scale_(spacing); +// const float length_max = scale_((2. / params.density) * spacing); + const float length_max = scale_((1000. / params.density) * spacing); std::vector merged_with(infill_ordered.size()); for (size_t i = 0; i < merged_with.size(); ++ i) merged_with[i] = i; @@ -956,12 +957,26 @@ void Fill::connect_infill(Polylines &&infill_ordered, const ExPolygon &boundary_ size_t idx_chain_last = 0; for (ConnectionCost &connection_cost : connections_sorted) { - const std::pair *cp1 = &map_infill_end_point_to_boundary[connection_cost.idx_first * 2 + 1]; - const std::pair *cp2 = &map_infill_end_point_to_boundary[(connection_cost.idx_first + 1) * 2]; + const std::pair *cp1 = &map_infill_end_point_to_boundary[connection_cost.idx_first * 2 + 1]; + const std::pair *cp1prev = cp1 - 1; + const std::pair *cp2 = &map_infill_end_point_to_boundary[(connection_cost.idx_first + 1) * 2]; + const std::pair *cp2next = cp2 + 1; assert(cp1->first == cp2->first); std::vector &contour_data = boundary_data[cp1->first]; if (connection_cost.reversed) std::swap(cp1, cp2); + // Mark the the other end points of the segments to be taken as consumed temporarily, so they will not be crossed + // by the new connection line. + bool prev_marked = false; + bool next_marked = false; + if (cp1prev->first == cp1->first && ! contour_data[cp1prev->second].point_consumed) { + contour_data[cp1prev->second].point_consumed = true; + prev_marked = true; + } + if (cp2next->first == cp1->first && ! contour_data[cp2next->second].point_consumed) { + contour_data[cp2next->second].point_consumed = true; + next_marked = true; + } if (could_take(contour_data, cp1->second, cp2->second)) { // Indices of the polygons to be connected. size_t idx_first = connection_cost.idx_first; @@ -980,6 +995,10 @@ void Fill::connect_infill(Polylines &&infill_ordered, const ExPolygon &boundary_ // Mark the second polygon as merged with the first one. merged_with[idx_second] = merged_with[idx_first]; } + if (prev_marked) + contour_data[cp1prev->second].point_consumed = false; + if (next_marked) + contour_data[cp2next->second].point_consumed = false; } polylines_out.reserve(polylines_out.size() + std::count_if(infill_ordered.begin(), infill_ordered.end(), [](const Polyline &pl) { return ! pl.empty(); })); for (Polyline &pl : infill_ordered) diff --git a/src/libslic3r/Fill/FillBase.hpp b/src/libslic3r/Fill/FillBase.hpp index d005f6e4a..517ce8383 100644 --- a/src/libslic3r/Fill/FillBase.hpp +++ b/src/libslic3r/Fill/FillBase.hpp @@ -111,9 +111,9 @@ protected: virtual std::pair _infill_direction(const Surface *surface) const; - void connect_infill(Polylines &&infill_ordered, const ExPolygon &boundary, Polylines &polylines_out, const FillParams ¶ms); - public: + static void connect_infill(Polylines &&infill_ordered, const ExPolygon &boundary, Polylines &polylines_out, double spacing, const FillParams ¶ms); + static coord_t _adjust_solid_spacing(const coord_t width, const coord_t distance); // Align a coordinate to a grid. The coordinate may be negative, diff --git a/src/libslic3r/Fill/FillGyroid.cpp b/src/libslic3r/Fill/FillGyroid.cpp index ed0586edf..913b0b0c0 100644 --- a/src/libslic3r/Fill/FillGyroid.cpp +++ b/src/libslic3r/Fill/FillGyroid.cpp @@ -185,6 +185,7 @@ void FillGyroid::_fill_surface_single( if (! polylines.empty()) // remove too small bits (larger than longer) polylines.erase( + //FIXME what is the small size? Removing tiny extrusions disconnects walls! std::remove_if(polylines.begin(), polylines.end(), [this](const Polyline &pl) { return pl.length() < scale_(this->spacing * 3); }), polylines.end()); @@ -195,7 +196,7 @@ void FillGyroid::_fill_surface_single( if (params.dont_connect) append(polylines_out, std::move(polylines)); else - this->connect_infill(std::move(polylines), expolygon, polylines_out, params); + this->connect_infill(std::move(polylines), expolygon, polylines_out, this->spacing, params); // new paths must be rotated back if (abs(infill_angle) >= EPSILON) { for (auto it = polylines_out.begin() + polylines_out_first_idx; it != polylines_out.end(); ++ it) diff --git a/src/libslic3r/ShortestPath.cpp b/src/libslic3r/ShortestPath.cpp index 9362e6043..8b5cbf483 100644 --- a/src/libslic3r/ShortestPath.cpp +++ b/src/libslic3r/ShortestPath.cpp @@ -1065,7 +1065,7 @@ std::vector chain_points(const Points &points, Point *start_near) } #ifndef NDEBUG -// #define DEBUG_SVG_OUTPUT + // #define DEBUG_SVG_OUTPUT #endif /* NDEBUG */ #ifdef DEBUG_SVG_OUTPUT @@ -1597,7 +1597,6 @@ static inline void reorder_by_two_exchanges_with_segment_flipping(std::vector &edges) { if (edges.size() < 3) { @@ -1681,6 +1680,173 @@ static inline void reorder_by_three_exchanges_with_segment_flipping(std::vector< } } +typedef Eigen::Matrix Matrixd; + +class FourOptCosts { +public: + FourOptCosts(const ConnectionCost &c1, const ConnectionCost &c2, const ConnectionCost &c3, const ConnectionCost &c4) : costs { &c1, &c2, &c3, &c4 } {} + + double operator()(size_t piece_idx, bool flipped) const { return flipped ? costs[piece_idx]->cost_flipped : costs[piece_idx]->cost; } + +private: + const ConnectionCost* costs[4]; +}; + +static inline std::pair minimum_crossover_cost( + const FourOptCosts &segment_costs, + const Matrixd &segment_end_point_distance_matrix, + const double cost_current) +{ + // Distance from the end of span1 to the start of span2. + auto end_point_distance = [&segment_end_point_distance_matrix](size_t span1, bool reversed1, bool flipped1, size_t span2, bool reversed2, bool flipped2) { + return segment_end_point_distance_matrix(span1 * 4 + (! reversed1) * 2 + flipped1, span2 * 4 + reversed2 * 2 + flipped2); + }; + auto connection_cost = [&segment_costs, end_point_distance]( + const size_t span1, bool reversed1, bool flipped1, + const size_t span2, bool reversed2, bool flipped2, + const size_t span3, bool reversed3, bool flipped3, + const size_t span4, bool reversed4, bool flipped4) { + // Calculate the cost of reverting chains and / or flipping segment orientations. + return segment_costs(span1, flipped1) + segment_costs(span2, flipped2) + segment_costs(span3, flipped3) + segment_costs(span4, flipped4) + + end_point_distance(span1, reversed1, flipped1, span2, reversed2, flipped2) + + end_point_distance(span2, reversed2, flipped2, span3, reversed3, flipped3) + + end_point_distance(span3, reversed3, flipped3, span4, reversed4, flipped4); + }; + +#ifndef NDEBUG + { + double c = connection_cost(0, false, false, 1, false, false, 2, false, false, 3, false, false); + assert(std::abs(c - cost_current) < SCALED_EPSILON); + } +#endif /* NDEBUG */ + + double cost_min = cost_current; + size_t flip_min = 0; // no flip, no improvement + for (size_t i = 0; i < (1 << 8); ++ i) { + // From the three combinations of 1,2,3 ordering, the other three are reversals of the first three. + size_t permutation = 0; + for (double c : { + (i == 0) ? cost_current : + connection_cost(0, (i & 1) != 0, (i & (1 << 1)) != 0, 1, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 2, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 3, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0), + connection_cost(0, (i & 1) != 0, (i & (1 << 1)) != 0, 1, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 3, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 2, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0), + connection_cost(0, (i & 1) != 0, (i & (1 << 1)) != 0, 2, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 1, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 3, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0), + connection_cost(0, (i & 1) != 0, (i & (1 << 1)) != 0, 2, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 3, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 1, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0), + connection_cost(0, (i & 1) != 0, (i & (1 << 1)) != 0, 3, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 1, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 2, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0), + connection_cost(0, (i & 1) != 0, (i & (1 << 1)) != 0, 3, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 2, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 1, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0), + connection_cost(1, (i & 1) != 0, (i & (1 << 1)) != 0, 0, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 2, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 3, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0), + connection_cost(1, (i & 1) != 0, (i & (1 << 1)) != 0, 0, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 3, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 2, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0), + connection_cost(1, (i & 1) != 0, (i & (1 << 1)) != 0, 2, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 0, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 3, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0), + connection_cost(1, (i & 1) != 0, (i & (1 << 1)) != 0, 3, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 0, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 2, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0), + connection_cost(2, (i & 1) != 0, (i & (1 << 1)) != 0, 0, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 1, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 3, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0), + connection_cost(2, (i & 1) != 0, (i & (1 << 1)) != 0, 1, (i & (1 << 2)) != 0, (i & (1 << 3)) != 0, 0, (i & (1 << 4)) != 0, (i & (1 << 5)) != 0, 3, (i & (1 << 6)) != 0, (i & (1 << 7)) != 0) + }) { + if (c < cost_min) { + cost_min = c; + flip_min = i + (permutation << 8); + } + ++ permutation; + } + } + return std::make_pair(cost_min, flip_min); +} + +static inline void reorder_by_three_exchanges_with_segment_flipping2(std::vector &edges) +{ + if (edges.size() < 3) { + reorder_by_two_exchanges_with_segment_flipping(edges); + return; + } + + std::vector connections(edges.size()); + std::vector edges_tmp(edges); + std::vector> connection_lengths(edges.size() - 1, std::pair(0., 0)); + std::vector connection_tried(edges.size(), false); + for (size_t iter = 0; iter < edges.size(); ++ iter) { + // Initialize connection costs and connection lengths. + for (size_t i = 1; i < edges.size(); ++ i) { + const FlipEdge &e1 = edges[i - 1]; + const FlipEdge &e2 = edges[i]; + ConnectionCost &c = connections[i]; + c = connections[i - 1]; + double l = (e2.p1 - e1.p2).norm(); + c.cost += l; + c.cost_flipped += (e2.p2 - e1.p1).norm(); + connection_lengths[i - 1] = std::make_pair(l, i); + } + std::sort(connection_lengths.begin(), connection_lengths.end(), [](const std::pair &l, const std::pair &r) { return l.first > r.first; }); + std::fill(connection_tried.begin(), connection_tried.end(), false); + size_t crossover1_pos_final = std::numeric_limits::max(); + size_t crossover2_pos_final = std::numeric_limits::max(); + size_t crossover3_pos_final = std::numeric_limits::max(); + size_t crossover_flip_final = 0; + // Distances between the end points of the four pieces of the current segment sequence. +#ifdef NDEBUG + Matrixd segment_end_point_distance_matrix(4 * 4, 4 * 4); +#else /* NDEBUG */ + Matrixd segment_end_point_distance_matrix = Matrixd::Constant(4 * 4, 4 * 4, std::numeric_limits::max()); +#endif /* NDEBUG */ + for (const std::pair &first_crossover_candidate : connection_lengths) { + double longest_connection_length = first_crossover_candidate.first; + size_t longest_connection_idx = first_crossover_candidate.second; + connection_tried[longest_connection_idx] = true; + // Find the second crossover connection with the lowest total chain cost. + size_t crossover_pos_min = std::numeric_limits::max(); + double crossover_cost_min = connections.back().cost; + for (size_t j = 1; j < connections.size(); ++ j) + if (! connection_tried[j]) { + for (size_t k = j + 1; k < connections.size(); ++ k) + if (! connection_tried[k]) { + size_t a = longest_connection_idx; + size_t b = j; + size_t c = k; + if (a > c) + std::swap(a, c); + if (a > b) + std::swap(a, b); + if (b > c) + std::swap(b, c); + const Vec2d* endpts[16] = { + &edges[0].p1, &edges[0].p2, &edges[a - 1].p2, &edges[a - 1].p1, + &edges[a].p1, &edges[a].p2, &edges[b - 1].p2, &edges[b - 1].p1, + &edges[b].p1, &edges[b].p2, &edges[c - 1].p2, &edges[c - 1].p1, + &edges[c].p1, &edges[c].p2, &edges.back().p2, &edges.back().p1 }; + for (size_t v = 0; v < 16; ++ v) { + const Vec2d &p1 = *endpts[v]; + for (size_t u = (v & (~3)) + 4; u < 16; ++ u) + segment_end_point_distance_matrix(u, v) = segment_end_point_distance_matrix(v, u) = (*endpts[u] - p1).norm(); + } + FourOptCosts segment_costs(connections[a - 1], connections[b - 1] - connections[a], connections[c - 1] - connections[b], connections.back() - connections[c]); + std::pair cost_and_flip = minimum_crossover_cost(segment_costs, segment_end_point_distance_matrix, connections.back().cost); + if (cost_and_flip.second > 0 && cost_and_flip.first < crossover_cost_min) { + crossover_cost_min = cost_and_flip.first; + crossover1_pos_final = a; + crossover2_pos_final = b; + crossover3_pos_final = c; + crossover_flip_final = cost_and_flip.second; + assert(crossover_cost_min < connections.back().cost + EPSILON); + } + } + } + if (crossover_flip_final > 0) { + // The cost of the chain with the proposed two crossovers has a lower total cost than the current chain. Apply the crossover. + break; + } else { + // Continue with another long candidate edge. + } + } + if (crossover_flip_final > 0) { + // Pair of cross over positions and flip / reverse constellation has been found, which improves the total cost of the connection. + // Perform a crossover. + do_crossover(edges, edges_tmp, std::make_pair(size_t(0), crossover1_pos_final), std::make_pair(crossover1_pos_final, crossover2_pos_final), + std::make_pair(crossover2_pos_final, crossover3_pos_final), std::make_pair(crossover3_pos_final, edges.size()), crossover_flip_final); + edges.swap(edges_tmp); + } else { + // No valid pair of cross over positions was found improving the total cost. Giving up. + break; + } + } +} + // Flip the sequences of polylines to lower the total length of connecting lines. static inline void improve_ordering_by_two_exchanges_with_segment_flipping(Polylines &polylines, bool fixed_start) { @@ -1707,7 +1873,8 @@ static inline void improve_ordering_by_two_exchanges_with_segment_flipping(Polyl #if 1 reorder_by_two_exchanges_with_segment_flipping(edges); #else - reorder_by_three_exchanges_with_segment_flipping(edges); + // reorder_by_three_exchanges_with_segment_flipping(edges); + reorder_by_three_exchanges_with_segment_flipping2(edges); #endif Polylines out; out.reserve(polylines.size());