diff --git a/src/libslic3r/ClipperUtils.cpp b/src/libslic3r/ClipperUtils.cpp index 25100b22fe..f053aea296 100644 --- a/src/libslic3r/ClipperUtils.cpp +++ b/src/libslic3r/ClipperUtils.cpp @@ -1205,7 +1205,7 @@ ExPolygons variable_offset_inner_ex(const ExPolygon &expoly, const std::vector& ds : deltas) + for (const std::vector& ds : deltas) for (float delta : ds) assert(delta <= 0.); assert(expoly.holes.size() + 1 == deltas.size()); diff --git a/src/libslic3r/ElephantFootCompensation.cpp b/src/libslic3r/ElephantFootCompensation.cpp index 0f4eb01350..828064c329 100644 --- a/src/libslic3r/ElephantFootCompensation.cpp +++ b/src/libslic3r/ElephantFootCompensation.cpp @@ -60,9 +60,9 @@ std::vector contour_distance(const EdgeGrid::Grid &grid, const size_t idx for (size_t axis = 0; axis < 2; ++ axis) { double dx = std::abs(dir(axis)); if (dx >= EPSILON) { - double tedge = (dir(axis) > 0) ? (double(bbox.max(axis)) - EPSILON - this->pt(axis)) : (this->pt(axis) - double(bbox.min(axis)) - EPSILON); + double tedge = (dir(axis) > 0) ? (double(bbox.max(axis)) - SCALED_EPSILON - this->pt(axis)) : (this->pt(axis) - double(bbox.min(axis)) - SCALED_EPSILON); if (tedge < dx) - t = tedge / dx; + t = std::min(t, tedge / dx); } } this->dir = dir; @@ -70,6 +70,7 @@ std::vector contour_distance(const EdgeGrid::Grid &grid, const size_t idx dir *= t; this->pt_end = (this->pt + dir).cast(); this->t_min = 1.; + assert(this->grid.bbox().contains(this->pt_start) && this->grid.bbox().contains(this->pt_end)); } bool operator()(coord_t iy, coord_t ix) { @@ -361,7 +362,7 @@ static inline void smooth_compensation_banded(const Points &contour, float band, } ExPolygon elephant_foot_compensation(const ExPolygon &input_expoly, const Flow &external_perimeter_flow, const double compensation) -{ +{ // The contour shall be wide enough to apply the external perimeter plus compensation on both sides. double min_contour_width = double(external_perimeter_flow.scaled_width() + external_perimeter_flow.scaled_spacing()); double scaled_compensation = scale_(compensation); @@ -369,39 +370,59 @@ ExPolygon elephant_foot_compensation(const ExPolygon &input_expoly, const Flow & // Make the search radius a bit larger for the averaging in contour_distance over a fan of rays to work. double search_radius = min_contour_width_compensated + min_contour_width * 0.5; - EdgeGrid::Grid grid; - ExPolygon simplified = input_expoly.simplify(SCALED_EPSILON).front(); - BoundingBox bbox = get_extents(simplified.contour); - bbox.offset(SCALED_EPSILON); - grid.set_bbox(bbox); - grid.create(simplified, coord_t(0.7 * search_radius)); - std::vector> deltas; - deltas.reserve(simplified.holes.size() + 1); - ExPolygon resampled(simplified); - double resample_interval = scale_(0.5); - for (size_t idx_contour = 0; idx_contour <= simplified.holes.size(); ++ idx_contour) { - Polygon &poly = (idx_contour == 0) ? resampled.contour : resampled.holes[idx_contour - 1]; - std::vector resampled_point_parameters; - poly.points = resample_polygon(poly.points, resample_interval, resampled_point_parameters); - std::vector dists = contour_distance(grid, idx_contour, poly.points, resampled_point_parameters, search_radius); - for (float &d : dists) { -// printf("Point %d, Distance: %lf\n", int(&d - dists.data()), unscale(d)); - // Convert contour width to available compensation distance. - if (d < min_contour_width) - d = 0.f; - else if (d > min_contour_width_compensated) - d = - float(scaled_compensation); - else - d = - (d - float(min_contour_width)) / 2.f; - assert(d >= - float(scaled_compensation) && d <= 0.f); + BoundingBox bbox = get_extents(input_expoly.contour); + Point bbox_size = bbox.size(); + ExPolygon out; + if (bbox_size.x() < min_contour_width_compensated + SCALED_EPSILON || + bbox_size.y() < min_contour_width_compensated + SCALED_EPSILON || + input_expoly.area() < min_contour_width_compensated * min_contour_width_compensated * 5.) + { + // The contour is tiny. Don't correct it. + out = input_expoly; + } + else + { + EdgeGrid::Grid grid; + ExPolygon simplified = input_expoly.simplify(SCALED_EPSILON).front(); + BoundingBox bbox = get_extents(simplified.contour); + bbox.offset(SCALED_EPSILON); + grid.set_bbox(bbox); + grid.create(simplified, coord_t(0.7 * search_radius)); + std::vector> deltas; + deltas.reserve(simplified.holes.size() + 1); + ExPolygon resampled(simplified); + double resample_interval = scale_(0.5); + for (size_t idx_contour = 0; idx_contour <= simplified.holes.size(); ++ idx_contour) { + Polygon &poly = (idx_contour == 0) ? resampled.contour : resampled.holes[idx_contour - 1]; + std::vector resampled_point_parameters; + poly.points = resample_polygon(poly.points, resample_interval, resampled_point_parameters); + std::vector dists = contour_distance(grid, idx_contour, poly.points, resampled_point_parameters, search_radius); + for (float &d : dists) { + // printf("Point %d, Distance: %lf\n", int(&d - dists.data()), unscale(d)); + // Convert contour width to available compensation distance. + if (d < min_contour_width) + d = 0.f; + else if (d > min_contour_width_compensated) + d = - float(scaled_compensation); + else + d = - (d - float(min_contour_width)) / 2.f; + assert(d >= - float(scaled_compensation) && d <= 0.f); + } + // smooth_compensation(dists, 0.4f, 10); + smooth_compensation_banded(poly.points, float(0.8 * resample_interval), dists, 0.3f, 3); + deltas.emplace_back(dists); } -// smooth_compensation(dists, 0.4f, 10); - smooth_compensation_banded(poly.points, float(0.8 * resample_interval), dists, 0.3f, 3); - deltas.emplace_back(dists); + + ExPolygons out_vec = variable_offset_inner_ex(resampled, deltas, 2.); + assert(out_vec.size() == 1); + if (out_vec.size() == 1) + out = std::move(out_vec.front()); + else + // Something went wrong, don't compensate. + out = input_expoly; } - ExPolygons out = variable_offset_inner_ex(resampled, deltas, 2.); - return out.front(); + return out; } ExPolygons elephant_foot_compensation(const ExPolygons &input, const Flow &external_perimeter_flow, const double compensation) diff --git a/src/libslic3r/ExtrusionEntity.hpp b/src/libslic3r/ExtrusionEntity.hpp index b22d85b657..b76991f1cc 100644 --- a/src/libslic3r/ExtrusionEntity.hpp +++ b/src/libslic3r/ExtrusionEntity.hpp @@ -267,6 +267,15 @@ public: //static inline std::string role_to_string(ExtrusionLoopRole role); +#ifndef NDEBUG + bool validate() const { + assert(this->first_point() == this->paths.back().polyline.points.back()); + for (size_t i = 1; i < paths.size(); ++ i) + assert(this->paths[i - 1].polyline.points.back() == this->paths[i].polyline.points.front()); + return true; + } +#endif /* NDEBUG */ + private: ExtrusionLoopRole m_loop_role; }; diff --git a/src/libslic3r/Fill/FillBase.cpp b/src/libslic3r/Fill/FillBase.cpp index 0ba75465f2..88eba9a51f 100644 --- a/src/libslic3r/Fill/FillBase.cpp +++ b/src/libslic3r/Fill/FillBase.cpp @@ -534,7 +534,8 @@ struct ContourPointData { // Verify whether the contour from point idx_start to point idx_end could be taken (whether all segments along the contour were not yet extruded). static bool could_take(const std::vector &contour_data, size_t idx_start, size_t idx_end) { - for (size_t i = idx_start; i < idx_end; ) { + assert(idx_start != idx_end); + for (size_t i = idx_start; i != idx_end; ) { if (contour_data[i].segment_consumed || contour_data[i].point_consumed) return false; if (++ i == contour_data.size()) @@ -899,63 +900,86 @@ 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 = scale_(this->spacing); + const double clip_distance = 3. * scale_(this->spacing); const double distance_colliding = scale_(this->spacing); mark_boundary_segments_touching_infill(boundary, boundary_data, bbox, infill_ordered, clip_distance, distance_colliding); } - // Chain infill_ordered. - //FIXME run the following loop through a heap sorted by the shortest perimeter edge that could be taken. - //length between two lines + // 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); - size_t idx_chain_last = 0; +// const float length_max = scale_((2. / params.density) * this->spacing); + const float length_max = scale_((1000. / params.density) * this->spacing); + std::vector merged_with(infill_ordered.size()); + for (size_t i = 0; i < merged_with.size(); ++ i) + merged_with[i] = i; + struct ConnectionCost { + ConnectionCost(size_t idx_first, double cost, bool reversed) : idx_first(idx_first), cost(cost), reversed(reversed) {} + size_t idx_first; + double cost; + bool reversed; + }; + std::vector connections_sorted; + connections_sorted.reserve(infill_ordered.size() * 2 - 2); for (size_t idx_chain = 1; idx_chain < infill_ordered.size(); ++ idx_chain) { - Polyline &pl1 = infill_ordered[idx_chain_last]; - Polyline &pl2 = infill_ordered[idx_chain]; + const Polyline &pl1 = infill_ordered[idx_chain - 1]; + const Polyline &pl2 = infill_ordered[idx_chain]; const std::pair *cp1 = &map_infill_end_point_to_boundary[(idx_chain - 1) * 2 + 1]; const std::pair *cp2 = &map_infill_end_point_to_boundary[idx_chain * 2]; - const Points &contour = boundary[cp1->first]; - std::vector &contour_data = boundary_data[cp1->first]; - bool valid = false; - bool reversed = false; + const std::vector &contour_data = boundary_data[cp1->first]; if (cp1->first == cp2->first) { // End points on the same contour. Try to connect them. - float param_lo = (cp1->second == 0) ? 0.f : contour_data[cp1->second].param; - float param_hi = (cp2->second == 0) ? 0.f : contour_data[cp2->second].param; + float param_lo = (cp1->second == 0) ? 0.f : contour_data[cp1->second].param; + float param_hi = (cp2->second == 0) ? 0.f : contour_data[cp2->second].param; float param_end = contour_data.front().param; + bool reversed = false; if (param_lo > param_hi) { std::swap(param_lo, param_hi); - std::swap(cp1, cp2); reversed = true; } assert(param_lo >= 0.f && param_lo <= param_end); assert(param_hi >= 0.f && param_hi <= param_end); - float dist1 = param_hi - param_lo; - float dist2 = param_lo + param_end - param_hi; - if (dist1 > dist2) { - std::swap(dist1, dist2); - std::swap(cp1, cp2); - reversed = ! reversed; - } - if (dist1 < length_max) { - // Try to connect the shorter path. - valid = could_take(contour_data, cp1->second, cp2->second); - // Try to connect the longer path. - if (! valid && dist2 < length_max) { - std::swap(cp1, cp2); - reversed = ! reversed; - valid = could_take(contour_data, cp1->second, cp2->second); - } - } + double len = param_hi - param_lo; + if (len < length_max) + connections_sorted.emplace_back(idx_chain - 1, len, reversed); + len = param_lo + param_end - param_hi; + if (len < length_max) + connections_sorted.emplace_back(idx_chain - 1, len, ! reversed); } - if (valid) - take(pl1, std::move(pl2), contour, contour_data, cp1->second, cp2->second, reversed); - else if (++ idx_chain_last < idx_chain) - infill_ordered[idx_chain_last] = std::move(pl2); } - infill_ordered.erase(infill_ordered.begin() + idx_chain_last + 1, infill_ordered.end()); - append(polylines_out, std::move(infill_ordered)); + std::sort(connections_sorted.begin(), connections_sorted.end(), [](const ConnectionCost& l, const ConnectionCost& r) { return l.cost < r.cost; }); + + 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]; + assert(cp1->first == cp2->first); + std::vector &contour_data = boundary_data[cp1->first]; + if (connection_cost.reversed) + std::swap(cp1, cp2); + if (could_take(contour_data, cp1->second, cp2->second)) { + // Indices of the polygons to be connected. + size_t idx_first = connection_cost.idx_first; + size_t idx_second = idx_first + 1; + for (size_t last = idx_first;;) { + size_t lower = merged_with[last]; + if (lower == last) { + merged_with[idx_first] = lower; + idx_first = lower; + break; + } + last = lower; + } + // Connect the two polygons using the boundary contour. + take(infill_ordered[idx_first], std::move(infill_ordered[idx_second]), boundary[cp1->first], contour_data, cp1->second, cp2->second, connection_cost.reversed); + // Mark the second polygon as merged with the first one. + merged_with[idx_second] = merged_with[idx_first]; + } + } + 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) + if (! pl.empty()) + polylines_out.emplace_back(std::move(pl)); } #endif diff --git a/src/libslic3r/Fill/FillGyroid.cpp b/src/libslic3r/Fill/FillGyroid.cpp index e7b4706acd..b12dfb2e74 100644 --- a/src/libslic3r/Fill/FillGyroid.cpp +++ b/src/libslic3r/Fill/FillGyroid.cpp @@ -166,7 +166,7 @@ void FillGyroid::_fill_surface_single( bb.merge(_align_to_grid(bb.min, Point(2*M_PI*distance, 2*M_PI*distance))); // generate pattern - Polylines polylines_square = make_gyroid_waves( + Polylines polylines = make_gyroid_waves( scale_(this->z), density_adjusted, this->spacing, @@ -174,22 +174,25 @@ void FillGyroid::_fill_surface_single( ceil(bb.size()(1) / distance) + 1.); // shift the polyline to the grid origin - for (Polyline &pl : polylines_square) + for (Polyline &pl : polylines) pl.translate(bb.min); - Polylines polylines_chained = chain_polylines(intersection_pl(polylines_square, to_polygons(expolygon))); + polylines = intersection_pl(polylines, to_polygons(expolygon)); - size_t polylines_out_first_idx = polylines_out.size(); - if (! polylines_chained.empty()) { - // connect lines + if (! polylines.empty()) + // remove too small bits (larger than longer) + polylines.erase( + std::remove_if(polylines.begin(), polylines.end(), [this](const Polyline &pl) { return pl.length() < scale_(this->spacing * 3); }), + polylines.end()); + + if (! polylines.empty()) { + polylines = chain_polylines(polylines); + // connect lines + size_t polylines_out_first_idx = polylines_out.size(); if (params.dont_connect) - append(polylines_out, std::move(polylines_chained)); + append(polylines_out, std::move(polylines)); else - this->connect_infill(std::move(polylines_chained), expolygon, polylines_out, params); - // remove too small bits (larger than longer) - polylines_out.erase( - std::remove_if(polylines_out.begin() + polylines_out_first_idx, polylines_out.end(), [this](const Polyline &pl){ return pl.length() < scale_(this->spacing * 3); }), - polylines_out.end()); + this->connect_infill(std::move(polylines), expolygon, polylines_out, 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/KDTreeIndirect.hpp b/src/libslic3r/KDTreeIndirect.hpp index 3cccfdafac..239008559c 100644 --- a/src/libslic3r/KDTreeIndirect.hpp +++ b/src/libslic3r/KDTreeIndirect.hpp @@ -46,9 +46,9 @@ public: if (indices.empty()) clear(); else { - // Allocate a next highest power of 2 nodes, because the incomplete binary tree will not have the leaves filled strictly from the left. + // Allocate enough memory for a full binary tree. m_nodes.assign(next_highest_power_of_2(indices.size() + 1), npos); - build_recursive(indices, 0, 0, 0, (int)(indices.size() - 1)); + build_recursive(indices, 0, 0, 0, indices.size() - 1); } indices.clear(); } @@ -81,7 +81,7 @@ public: private: // Build a balanced tree by splitting the input sequence by an axis aligned plane at a dimension. - void build_recursive(std::vector &input, size_t node, int dimension, int left, int right) + void build_recursive(std::vector &input, size_t node, const size_t dimension, const size_t left, const size_t right) { if (left > right) return; @@ -94,54 +94,56 @@ private: return; } - // Partition the input sequence to two equal halves. - int center = (left + right) >> 1; + // Partition the input to left / right pieces of the same length to produce a balanced tree. + size_t center = (left + right) / 2; partition_input(input, dimension, left, right, center); // Insert a node into the tree. m_nodes[node] = input[center]; - // Partition the left and right subtrees. - size_t next_dimension = (++ dimension == NumDimensions) ? 0 : dimension; - build_recursive(input, (node << 1) + 1, next_dimension, left, center - 1); - build_recursive(input, (node << 1) + 2, next_dimension, center + 1, right); + // Build up the left / right subtrees. + size_t next_dimension = dimension; + if (++ next_dimension == NumDimensions) + next_dimension = 0; + if (center > left) + build_recursive(input, node * 2 + 1, next_dimension, left, center - 1); + build_recursive(input, node * 2 + 2, next_dimension, center + 1, right); } - // Partition the input m_nodes at k using QuickSelect method. + // Partition the input m_nodes at "k" and "dimension" using the QuickSelect method: // https://en.wikipedia.org/wiki/Quickselect - void partition_input(std::vector &input, int dimension, int left, int right, int k) const + // Items left of the k'th item are lower than the k'th item in the "dimension", + // items right of the k'th item are higher than the k'th item in the "dimension", + void partition_input(std::vector &input, const size_t dimension, size_t left, size_t right, const size_t k) const { while (left < right) { - // Guess the k'th element. - // Pick the pivot as a median of first, center and last value. - // Sort first, center and last values. - int center = (left + right) >> 1; - auto left_value = this->coordinate(input[left], dimension); - auto center_value = this->coordinate(input[center], dimension); - auto right_value = this->coordinate(input[right], dimension); - if (center_value < left_value) { - std::swap(input[left], input[center]); - std::swap(left_value, center_value); + size_t center = (left + right) / 2; + CoordType pivot; + { + // Bubble sort the input[left], input[center], input[right], so that a median of the three values + // will end up in input[center]. + CoordType left_value = this->coordinate(input[left], dimension); + CoordType center_value = this->coordinate(input[center], dimension); + CoordType right_value = this->coordinate(input[right], dimension); + if (left_value > center_value) { + std::swap(input[left], input[center]); + std::swap(left_value, center_value); + } + if (left_value > right_value) { + std::swap(input[left], input[right]); + right_value = left_value; + } + if (center_value > right_value) { + std::swap(input[center], input[right]); + center_value = right_value; + } + pivot = center_value; } - if (right_value < left_value) { - std::swap(input[left], input[right]); - std::swap(left_value, right_value); - } - if (right_value < center_value) { - std::swap(input[center], input[right]); - // No need to do that, result is not used. - // std::swap(center_value, right_value); - } - // Only two or three values are left and those are sorted already. - if (left + 3 > right) + if (right <= left + 2) + // The interval is already sorted. break; - // left and right items are already at their correct positions. - // input[left].point[dimension] <= input[center].point[dimension] <= input[right].point[dimension] - // Move the pivot to the (right - 1) position. - std::swap(input[center], input[right - 1]); - // Pivot value. - double pivot = this->coordinate(input[right - 1], dimension); + size_t i = left; + size_t j = right - 1; + std::swap(input[center], input[j]); // Partition the set based on the pivot. - int i = left; - int j = right - 1; for (;;) { // Skip left points that are already at correct positions. // Search will certainly stop at position (right - 1), which stores the pivot. @@ -153,7 +155,7 @@ private: std::swap(input[i], input[j]); } // Restore pivot to the center of the sequence. - std::swap(input[i], input[right]); + std::swap(input[i], input[right - 1]); // Which side the kth element is in? if (k < i) right = i - 1; @@ -173,7 +175,7 @@ private: return; // Left / right child node index. - size_t left = (node << 1) + 1; + size_t left = node * 2 + 1; size_t right = left + 1; unsigned int mask = visitor(m_nodes[node], dimension); if ((mask & (unsigned int)VisitorReturnMask::STOP) == 0) { diff --git a/src/libslic3r/ShortestPath.cpp b/src/libslic3r/ShortestPath.cpp index a0fb050057..b38655e68a 100644 --- a/src/libslic3r/ShortestPath.cpp +++ b/src/libslic3r/ShortestPath.cpp @@ -237,11 +237,19 @@ std::vector> chain_segments_greedy_constrained_reversals // Chain the end points: find (num_segments - 1) shortest links not forming bifurcations or loops. assert(num_segments >= 2); +#ifndef NDEBUG + double distance_taken_last = 0.; +#endif /* NDEBUG */ for (int iter = int(num_segments) - 2;; -- iter) { assert(validate_graph_and_queue()); // Take the first end point, for which the link points to the currently closest valid neighbor. EndPoint &end_point1 = *queue.top(); - assert(end_point1.edge_out != nullptr); +#ifndef NDEBUG + // Each edge added shall be longer than the previous one taken. + assert(end_point1.distance_out > distance_taken_last - SCALED_EPSILON); + distance_taken_last = end_point1.distance_out; +#endif /* NDEBUG */ + assert(end_point1.edge_out != nullptr); // No point on the queue may be connected yet. assert(end_point1.chain_id == 0); // Take the closest end point to the first end point, @@ -313,6 +321,10 @@ std::vector> chain_segments_greedy_constrained_reversals assert(next_idx < end_points.size()); end_point1.edge_out = &end_points[next_idx]; end_point1.distance_out = (end_points[next_idx].pos - end_point1.pos).squaredNorm(); +#ifndef NDEBUG + // Each edge shall be longer than the last one removed from the queue. + assert(end_point1.distance_out > distance_taken_last - SCALED_EPSILON); +#endif /* NDEBUG */ // Update position of this end point in the queue based on the distance calculated at the line above. queue.update(end_point1.heap_idx); //FIXME Remove the other end point from the KD tree. @@ -460,18 +472,206 @@ std::vector chain_points(const Points &points, Point *start_near) return out; } +// Flip the sequences of polylines to lower the total length of connecting lines. +// #define DEBUG_SVG_OUTPUT +static inline void improve_ordering_by_segment_flipping(Polylines &polylines, bool fixed_start) +{ +#ifndef NDEBUG + auto cost = [&polylines]() { + double sum = 0.; + for (size_t i = 1; i < polylines.size(); ++i) + sum += (polylines[i].first_point() - polylines[i - 1].last_point()).cast().norm(); + return sum; + }; + double cost_initial = cost(); + + static int iRun = 0; + ++ iRun; + BoundingBox bbox = get_extents(polylines); +#ifdef DEBUG_SVG_OUTPUT + { + SVG svg(debug_out_path("improve_ordering_by_segment_flipping-initial-%d.svg", iRun).c_str(), bbox); + svg.draw(polylines); + for (size_t i = 1; i < polylines.size(); ++ i) + svg.draw(Line(polylines[i - 1].last_point(), polylines[i].first_point()), "red"); + } +#endif /* DEBUG_SVG_OUTPUT */ +#endif /* NDEBUG */ + + struct Connection { + Connection(size_t heap_idx = std::numeric_limits::max(), bool flipped = false) : heap_idx(heap_idx), flipped(flipped) {} + // Position of this object on MutablePriorityHeap. + size_t heap_idx; + // Is segment_idx flipped? + bool flipped; + + double squaredNorm(const Polylines &polylines, const std::vector &connections) const + { return ((this + 1)->start_point(polylines, connections) - this->end_point(polylines, connections)).squaredNorm(); } + double norm(const Polylines &polylines, const std::vector &connections) const + { return sqrt(this->squaredNorm(polylines, connections)); } + double squaredNorm(const Polylines &polylines, const std::vector &connections, bool try_flip1, bool try_flip2) const + { return ((this + 1)->start_point(polylines, connections, try_flip2) - this->end_point(polylines, connections, try_flip1)).squaredNorm(); } + double norm(const Polylines &polylines, const std::vector &connections, bool try_flip1, bool try_flip2) const + { return sqrt(this->squaredNorm(polylines, connections, try_flip1, try_flip2)); } + Vec2d start_point(const Polylines &polylines, const std::vector &connections, bool flip = false) const + { const Polyline &pl = polylines[this - connections.data()]; return ((this->flipped == flip) ? pl.points.front() : pl.points.back()).cast(); } + Vec2d end_point(const Polylines &polylines, const std::vector &connections, bool flip = false) const + { const Polyline &pl = polylines[this - connections.data()]; return ((this->flipped == flip) ? pl.points.back() : pl.points.front()).cast(); } + + bool in_queue() const { return this->heap_idx != std::numeric_limits::max(); } + void flip() { this->flipped = ! this->flipped; } + }; + std::vector connections(polylines.size()); + +#ifndef NDEBUG + auto cost_flipped = [fixed_start, &polylines, &connections]() { + assert(! fixed_start || ! connections.front().flipped); + double sum = 0.; + for (size_t i = 1; i < polylines.size(); ++ i) + sum += connections[i - 1].norm(polylines, connections); + return sum; + }; + double cost_prev = cost_flipped(); + assert(std::abs(cost_initial - cost_prev) < SCALED_EPSILON); + + auto print_statistics = [&polylines, &connections]() { +#if 0 + for (size_t i = 1; i < polylines.size(); ++ i) + printf("Connecting %d with %d: Current length %lf flip(%d, %d), left flipped: %lf, right flipped: %lf, both flipped: %lf, \n", + int(i - 1), int(i), + unscale(connections[i - 1].norm(polylines, connections)), + int(connections[i - 1].flipped), int(connections[i].flipped), + unscale(connections[i - 1].norm(polylines, connections, true, false)), + unscale(connections[i - 1].norm(polylines, connections, false, true)), + unscale(connections[i - 1].norm(polylines, connections, true, true))); +#endif + }; + print_statistics(); +#endif /* NDEBUG */ + + // Initialize a MutablePriorityHeap of connections between polylines. + auto queue = make_mutable_priority_queue( + [](Connection *connection, size_t idx){ connection->heap_idx = idx; }, + // Sort by decreasing connection distance. + [&polylines, &connections](Connection *l, Connection *r){ return l->squaredNorm(polylines, connections) > r->squaredNorm(polylines, connections); }); + queue.reserve(polylines.size() - 1); + for (size_t i = 0; i + 1 < polylines.size(); ++ i) + queue.push(&connections[i]); + + static constexpr size_t itercnt = 100; + size_t iter = 0; + for (; ! queue.empty() && iter < itercnt; ++ iter) { + Connection &connection = *queue.top(); + queue.pop(); + connection.heap_idx = std::numeric_limits::max(); + size_t idx_first = &connection - connections.data(); + // Try to flip segments starting with idx_first + 1 to the end. + // Calculate the last segment to be flipped to improve the total path length. + double length_current = connection.norm(polylines, connections); + double length_flipped = connection.norm(polylines, connections, false, true); + int best_idx_forward = int(idx_first); + double best_improvement_forward = 0.; + for (size_t i = idx_first + 1; i + 1 < connections.size(); ++ i) { + length_current += connections[i].norm(polylines, connections); + double this_improvement = length_current - length_flipped - connections[i].norm(polylines, connections, true, false); + length_flipped += connections[i].norm(polylines, connections, true, true); + if (this_improvement > best_improvement_forward) { + best_improvement_forward = this_improvement; + best_idx_forward = int(i); + } +// if (length_flipped > 1.5 * length_current) +// break; + } + if (length_current - length_flipped > best_improvement_forward) + // Best improvement by flipping up to the end. + best_idx_forward = int(connections.size()) - 1; + // Try to flip segments starting with idx_first - 1 to the start. + // Calculate the last segment to be flipped to improve the total path length. + length_current = connection.norm(polylines, connections); + length_flipped = connection.norm(polylines, connections, true, false); + int best_idx_backwards = int(idx_first); + double best_improvement_backwards = 0.; + for (int i = int(idx_first) - 1; i >= 0; -- i) { + length_current += connections[i].norm(polylines, connections); + double this_improvement = length_current - length_flipped - connections[i].norm(polylines, connections, false, true); + length_flipped += connections[i].norm(polylines, connections, true, true); + if (this_improvement > best_improvement_backwards) { + best_improvement_backwards = this_improvement; + best_idx_backwards = int(i); + } +// if (length_flipped > 1.5 * length_current) +// break; + } + if (! fixed_start && length_current - length_flipped > best_improvement_backwards) + // Best improvement by flipping up to the start including the first polyline. + best_idx_backwards = -1; + int update_begin = int(idx_first); + int update_end = best_idx_forward; + if (best_improvement_backwards > 0. && best_improvement_backwards > best_improvement_forward) { + // Flip the sequence of polylines from idx_first to best_improvement_forward + 1. + update_begin = best_idx_backwards; + update_end = int(idx_first); + } + assert(update_begin <= update_end); + if (update_begin == update_end) + continue; + for (int i = update_begin + 1; i <= update_end; ++ i) + connections[i].flip(); + +#ifndef NDEBUG + double cost = cost_flipped(); + assert(cost < cost_prev); + cost_prev = cost; + print_statistics(); +#endif /* NDEBUG */ + + update_end = std::min(update_end + 1, int(connections.size()) - 1); + for (int i = std::max(0, update_begin); i < update_end; ++ i) { + Connection &c = connections[i]; + if (c.in_queue()) + queue.update(c.heap_idx); + else + queue.push(&c); + } + } + + // Flip the segments based on the flip flag. + for (Polyline &pl : polylines) + if (connections[&pl - polylines.data()].flipped) + pl.reverse(); + +#ifndef NDEBUG + double cost_final = cost(); +#ifdef DEBUG_SVG_OUTPUT + { + SVG svg(debug_out_path("improve_ordering_by_segment_flipping-final-%d.svg", iRun).c_str(), bbox); + svg.draw(polylines); + for (size_t i = 1; i < polylines.size(); ++ i) + svg.draw(Line(polylines[i - 1].last_point(), polylines[i].first_point()), "red"); + } +#endif /* DEBUG_SVG_OUTPUT */ +#endif /* NDEBUG */ + + assert(cost_final <= cost_prev); + assert(cost_final <= cost_initial); +} + Polylines chain_polylines(Polylines &&polylines, const Point *start_near) { - auto segment_end_point = [&polylines](size_t idx, bool first_point) -> const Point& { return first_point ? polylines[idx].first_point() : polylines[idx].last_point(); }; - std::vector> ordered = chain_segments_greedy(segment_end_point, polylines.size(), start_near); Polylines out; - out.reserve(polylines.size()); - for (auto &segment_and_reversal : ordered) { - out.emplace_back(std::move(polylines[segment_and_reversal.first])); - if (segment_and_reversal.second) - out.back().reverse(); + if (! polylines.empty()) { + auto segment_end_point = [&polylines](size_t idx, bool first_point) -> const Point& { return first_point ? polylines[idx].first_point() : polylines[idx].last_point(); }; + std::vector> ordered = chain_segments_greedy(segment_end_point, polylines.size(), start_near); + out.reserve(polylines.size()); + for (auto &segment_and_reversal : ordered) { + out.emplace_back(std::move(polylines[segment_and_reversal.first])); + if (segment_and_reversal.second) + out.back().reverse(); + } + if (out.size() > 1) + improve_ordering_by_segment_flipping(out, start_near != nullptr); } - return out; + return out; } template static inline T chain_path_items(const Points &points, const T &items) diff --git a/tests/libslic3r/test_elephant_foot_compensation.cpp b/tests/libslic3r/test_elephant_foot_compensation.cpp index 9a7e652648..e5ea97e687 100644 --- a/tests/libslic3r/test_elephant_foot_compensation.cpp +++ b/tests/libslic3r/test_elephant_foot_compensation.cpp @@ -222,6 +222,21 @@ static ExPolygon vase_with_fins() SCENARIO("Elephant foot compensation", "[ElephantFoot]") { + GIVEN("Tiny contour") { + ExPolygon expoly({ { 133382606, 94912473 }, { 134232493, 95001115 }, { 133783926, 95159440 }, { 133441897, 95180666 }, { 133408242, 95191984 }, { 133339012, 95166830 }, { 132991642, 95011087 }, { 133206549, 94908304 } }); + WHEN("Compensated") { + ExPolygon expoly_compensated = elephant_foot_compensation(expoly, Flow(0.419999987f, 0.2f, 0.4f, false), 0.2f); +#ifdef TESTS_EXPORT_SVGS + SVG::export_expolygons(debug_out_path("elephant_foot_compensation_tiny.svg").c_str(), + { { { expoly }, { "gray", "black", "blue", coord_t(scale_(0.02)), 0.5f, "black", coord_t(scale_(0.05)) } }, + { { expoly_compensated }, { "gray", "black", "blue", coord_t(scale_(0.02)), 0.5f, "black", coord_t(scale_(0.05)) } } }); +#endif /* TESTS_EXPORT_SVGS */ + THEN("Tiny contour is not compensated") { + REQUIRE(expoly_compensated == expoly); + } + } + } + GIVEN("Large box") { ExPolygon expoly( { {50000000, 50000000 }, { 0, 50000000 }, { 0, 0 }, { 50000000, 0 } } ); WHEN("Compensated") { diff --git a/tests/libslic3r/test_geometry.cpp b/tests/libslic3r/test_geometry.cpp index 22755c2622..4a800b3d3f 100644 --- a/tests/libslic3r/test_geometry.cpp +++ b/tests/libslic3r/test_geometry.cpp @@ -252,15 +252,39 @@ SCENARIO("Circle Fit, TaubinFit with Newton's method", "[Geometry]") { } } -TEST_CASE("Chained path working correctly", "[Geometry]"){ - // if chained_path() works correctly, these points should be joined with no diagonal paths - // (thus 26 units long) - std::vector points = {Point(26,26),Point(52,26),Point(0,26),Point(26,52),Point(26,0),Point(0,52),Point(52,52),Point(52,0)}; - std::vector indices = chain_points(points); - for (Points::size_type i = 0; i + 1 < indices.size(); ++ i) { - double dist = (points.at(indices.at(i)).cast() - points.at(indices.at(i+1)).cast()).norm(); - REQUIRE(std::abs(dist-26) <= EPSILON); - } +SCENARIO("Path chaining", "[Geometry]") { + GIVEN("A path") { + std::vector points = { Point(26,26),Point(52,26),Point(0,26),Point(26,52),Point(26,0),Point(0,52),Point(52,52),Point(52,0) }; + THEN("Chained with no diagonals (thus 26 units long)") { + std::vector indices = chain_points(points); + for (Points::size_type i = 0; i + 1 < indices.size(); ++ i) { + double dist = (points.at(indices.at(i)).cast() - points.at(indices.at(i+1)).cast()).norm(); + REQUIRE(std::abs(dist-26) <= EPSILON); + } + } + } + GIVEN("Loop pieces") { + Point a { 2185796, 19058485 }; + Point b { 3957902, 18149382 }; + Point c { 2912841, 18790564 }; + Point d { 2831848, 18832390 }; + Point e { 3179601, 18627769 }; + Point f { 3137952, 18653370 }; + Polylines polylines = { { a, b }, + { c, d }, + { e, f }, + { d, a }, + { f, c }, + { b, e } }; + Polylines chained = chain_polylines(polylines, &a); + THEN("Connected without a gap") { + for (size_t i = 0; i < chained.size(); ++i) { + const Polyline &pl1 = (i == 0) ? chained.back() : chained[i - 1]; + const Polyline &pl2 = chained[i]; + REQUIRE(pl1.points.back() == pl2.points.front()); + } + } + } } SCENARIO("Line distances", "[Geometry]"){