3DPrinting/PrusaSlicer-NonPlainar
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fix of infill connecting along perimeter lines,

Browse Source 
new 3-opt iterative improvement of infill path (currently disabled,
it is extremely slow)
This commit is contained in:
bubnikv 2019-12-02 15:01:52 +01:00
parent b2b5df7d43
commit afa72da9d1
5 changed files with 259 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/libslic3r/Fill/Fill3DHoneycomb.cpp
View File

@ -169,7 +169,7 @@ void Fill3DHoneycomb::_fill_surface_single(
if (params.dont_connect) if (params.dont_connect)
append(polylines_out, std::move(polylines_chained)); append(polylines_out, std::move(polylines_chained));
else 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);
} }
} }

149
src/libslic3r/Fill/FillBase.cpp
View File

@ -610,16 +610,15 @@ static inline SegmentPoint clip_start_segment_and_point(const Points &polyline,
// Initialized to "invalid". // Initialized to "invalid".
SegmentPoint out; SegmentPoint out;
if (polyline.size() >= 2) { if (polyline.size() >= 2) {
const double d2 = distance * distance;
Vec2d pt_prev = polyline.front().cast<double>(); Vec2d pt_prev = polyline.front().cast<double>();
for (int i = 1; i < polyline.size(); ++ i) { for (int i = 1; i < polyline.size(); ++ i) {
Vec2d pt = polyline[i].cast<double>(); Vec2d pt = polyline[i].cast<double>();
Vec2d v = pt - pt_prev; Vec2d v = pt - pt_prev;
double l2 = v.squaredNorm(); double l2 = v.squaredNorm();
if (l2 > d2) { if (l2 > distance * distance) {
out.idx_segment = i; out.idx_segment = i;
out.t = distance / sqrt(l2); out.t = distance / sqrt(l2);
out.point = pt + out.t * v; out.point = pt_prev + out.t * v;
break; break;
} }
distance -= sqrt(l2); distance -= sqrt(l2);
@ -636,16 +635,17 @@ static inline SegmentPoint clip_end_segment_and_point(const Points &polyline, do
// Initialized to "invalid". // Initialized to "invalid".
SegmentPoint out; SegmentPoint out;
if (polyline.size() >= 2) { if (polyline.size() >= 2) {
const double d2 = distance * distance;
Vec2d pt_next = polyline.back().cast<double>(); Vec2d pt_next = polyline.back().cast<double>();
for (int i = int(polyline.size()) - 2; i >= 0; -- i) { for (int i = int(polyline.size()) - 2; i >= 0; -- i) {
Vec2d pt = polyline[i].cast<double>(); Vec2d pt = polyline[i].cast<double>();
Vec2d v = pt - pt_next; Vec2d v = pt - pt_next;
double l2 = v.squaredNorm(); double l2 = v.squaredNorm();
if (l2 > d2) { if (l2 > distance * distance) {
out.idx_segment = i; out.idx_segment = i;
out.t = distance / sqrt(l2); 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; break;
} }
distance -= sqrt(l2); distance -= sqrt(l2);
@ -655,21 +655,26 @@ static inline SegmentPoint clip_end_segment_and_point(const Points &polyline, do
return out; 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) static inline double segment_point_distance_squared(const Vec2d &p1a, const Vec2d &p1b, const Vec2d &p2)
{ {
const Vec2d v = p1b - p1a; const Vec2d v = p1b - p1a;
const Vec2d va = p2 - p1a; const double l2 = v.squaredNorm();
const double l2 = v.squaredNorm();
if (l2 < EPSILON) if (l2 < EPSILON)
// p1a == p1b // p1a == p1b
return va.squaredNorm(); return (p2 - p1a).squaredNorm();
// Project p2 onto the (p1a, p1b) segment. return segment_point_distance_squared(p1a, p1b, v, v.squaredNorm(), p2);
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();
} }
// Distance to the closest point of line. // 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(); double l2_2 = v2.squaredNorm();
if (l2_2 < EPSILON) if (l2_2 < EPSILON)
// p2a == p2b: Return distance of p2a from the (p1a, p1b) segment. // p2a == p2b: Return distance of p2a from the (p1a, p1b) segment.
return segment_point_distance_squared(p1a, p1b, p2a); return segment_point_distance_squared(p1a, p1b, v1, l1_2, 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<double>::max();
};
return std::min( return std::min(
project_p2a_p2b_onto_seg_p1a_p1b(p1a, p1b, p2a, p2b, v1, l1_2), std::min(segment_point_distance_squared(p1a, p1b, v1, l1_2, p2a), segment_point_distance_squared(p1a, p1b, v1, l1_2, p2b)),
project_p2a_p2b_onto_seg_p1a_p1b(p2a, p2b, p1a, p1b, v2, l2_2)); 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. // 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<double>(); const Vec2d seg_pt2 = segment.second.cast<double>();
if (min_distance_of_segments(seg_pt1, seg_pt2, *this->pt1, *this->pt2) < this->dist2_max) { 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. // 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; bdp.segment_consumed = true;
// There is no need for checking seg_pt2 as it will be checked the next time. // 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; 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<coord_t>(), this->pt2->cast<coord_t>()), "magenta");
}
#endif
} }
} }
// Continue traversing the grid along the edge. // Continue traversing the grid along the edge.
@ -780,6 +768,7 @@ void mark_boundary_segments_touching_infill(
BoundingBoxf bboxf(boundary_bbox.min.cast<double>(), boundary_bbox.max.cast<double>()); BoundingBoxf bboxf(boundary_bbox.min.cast<double>(), boundary_bbox.max.cast<double>());
bboxf.offset(- SCALED_EPSILON); bboxf.offset(- SCALED_EPSILON);
for (const Polyline &polyline : infill) { for (const Polyline &polyline : infill) {
// Clip the infill polyline by the Eucledian distance along the polyline. // Clip the infill polyline by the Eucledian distance along the polyline.
SegmentPoint start_point = clip_start_segment_and_point(polyline.points, clip_distance); 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); visitor.init(pt1d, pt2d);
grid.visit_cells_intersecting_thick_line(pt1, pt2, distance_colliding, visitor); grid.visit_cells_intersecting_thick_line(pt1, pt2, distance_colliding, visitor);
#else #else
Vec2d pt1 = (point_idx == start_point.idx_segment) ? start_point.point : polyline.points[point_idx].cast<double>(); Vec2d pt1 = (point_idx == start_point.idx_segment) ? start_point.point : polyline.points[point_idx ].cast<double>();
Vec2d pt2 = (point_idx == end_point .idx_segment) ? end_point .point : polyline.points[point_idx].cast<double>(); Vec2d pt2 = (point_idx == end_point .idx_segment) ? end_point .point : polyline.points[point_idx + 1].cast<double>();
#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<coord_t>(), pt2.cast<coord_t>()), "magenta", scale_(0.1));
}
#endif
visitor.init(pt1, pt2); visitor.init(pt1, pt2);
// Simulate tracing of a thick line. This only works reliably if distance_colliding <= grid cell size. // Simulate tracing of a thick line. This only works reliably if distance_colliding <= grid cell size.
Vec2d v = (pt2 - pt1).normalized() * distance_colliding; 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 &params) void Fill::connect_infill(Polylines &&infill_ordered, const ExPolygon &boundary_src, Polylines &polylines_out, const double spacing, const FillParams &params)
{ {
assert(! infill_ordered.empty()); assert(! infill_ordered.empty());
assert(! boundary_src.contour.points.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. // 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); // @supermerill used 2. * scale_(spacing)
const double clip_distance = 3. * scale_(this->spacing); const double clip_distance = 3. * scale_(spacing);
const double distance_colliding = scale_(this->spacing); const double distance_colliding = 1.1 * scale_(spacing);
mark_boundary_segments_touching_infill(boundary, boundary_data, bbox, infill_ordered, clip_distance, distance_colliding); 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. // 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_(spacing);
// const float length_max = scale_((2. / params.density) * this->spacing); // const float length_max = scale_((2. / params.density) * spacing);
const float length_max = scale_((1000. / params.density) * this->spacing); const float length_max = scale_((1000. / params.density) * spacing);
std::vector<size_t> merged_with(infill_ordered.size()); std::vector<size_t> merged_with(infill_ordered.size());
for (size_t i = 0; i < merged_with.size(); ++ i) for (size_t i = 0; i < merged_with.size(); ++ i)
merged_with[i] = 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; size_t idx_chain_last = 0;
for (ConnectionCost &connection_cost : connections_sorted) { for (ConnectionCost &connection_cost : connections_sorted) {
const std::pair<size_t, size_t> *cp1 = &map_infill_end_point_to_boundary[connection_cost.idx_first * 2 + 1]; const std::pair<size_t, size_t> *cp1 = &map_infill_end_point_to_boundary[connection_cost.idx_first * 2 + 1];
const std::pair<size_t, size_t> *cp2 = &map_infill_end_point_to_boundary[(connection_cost.idx_first + 1) * 2]; const std::pair<size_t, size_t> *cp1prev = cp1 - 1;
const std::pair<size_t, size_t> *cp2 = &map_infill_end_point_to_boundary[(connection_cost.idx_first + 1) * 2];
const std::pair<size_t, size_t> *cp2next = cp2 + 1;
assert(cp1->first == cp2->first); assert(cp1->first == cp2->first);
std::vector<ContourPointData> &contour_data = boundary_data[cp1->first]; std::vector<ContourPointData> &contour_data = boundary_data[cp1->first];
if (connection_cost.reversed) if (connection_cost.reversed)
std::swap(cp1, cp2); 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)) { if (could_take(contour_data, cp1->second, cp2->second)) {
// Indices of the polygons to be connected. // Indices of the polygons to be connected.
size_t idx_first = connection_cost.idx_first; 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. // Mark the second polygon as merged with the first one.
merged_with[idx_second] = merged_with[idx_first]; 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(); })); 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) for (Polyline &pl : infill_ordered)

4
src/libslic3r/Fill/FillBase.hpp
View File

@ -111,9 +111,9 @@ protected:
virtual std::pair<float, Point> _infill_direction(const Surface *surface) const; virtual std::pair<float, Point> _infill_direction(const Surface *surface) const;
void connect_infill(Polylines &&infill_ordered, const ExPolygon &boundary, Polylines &polylines_out, const FillParams &params);
public: public:
static void connect_infill(Polylines &&infill_ordered, const ExPolygon &boundary, Polylines &polylines_out, double spacing, const FillParams &params);
static coord_t _adjust_solid_spacing(const coord_t width, const coord_t distance); 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, // Align a coordinate to a grid. The coordinate may be negative,

3
src/libslic3r/Fill/FillGyroid.cpp
View File

@ -185,6 +185,7 @@ void FillGyroid::_fill_surface_single(
if (! polylines.empty()) if (! polylines.empty())
// remove too small bits (larger than longer) // remove too small bits (larger than longer)
polylines.erase( 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); }), std::remove_if(polylines.begin(), polylines.end(), [this](const Polyline &pl) { return pl.length() < scale_(this->spacing * 3); }),
polylines.end()); polylines.end());
@ -195,7 +196,7 @@ void FillGyroid::_fill_surface_single(
if (params.dont_connect) if (params.dont_connect)
append(polylines_out, std::move(polylines)); append(polylines_out, std::move(polylines));
else 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 // new paths must be rotated back
if (abs(infill_angle) >= EPSILON) { if (abs(infill_angle) >= EPSILON) {
for (auto it = polylines_out.begin() + polylines_out_first_idx; it != polylines_out.end(); ++ it) for (auto it = polylines_out.begin() + polylines_out_first_idx; it != polylines_out.end(); ++ it)

173
src/libslic3r/ShortestPath.cpp
View File

@ -1065,7 +1065,7 @@ std::vector<size_t> chain_points(const Points &points, Point *start_near)
} }
#ifndef NDEBUG #ifndef NDEBUG
// #define DEBUG_SVG_OUTPUT // #define DEBUG_SVG_OUTPUT
#endif /* NDEBUG */ #endif /* NDEBUG */
#ifdef DEBUG_SVG_OUTPUT #ifdef DEBUG_SVG_OUTPUT
@ -1597,7 +1597,6 @@ static inline void reorder_by_two_exchanges_with_segment_flipping(std::vector<Fl
} }
} }
static inline void reorder_by_three_exchanges_with_segment_flipping(std::vector<FlipEdge> &edges) static inline void reorder_by_three_exchanges_with_segment_flipping(std::vector<FlipEdge> &edges)
{ {
if (edges.size() < 3) { if (edges.size() < 3) {
@ -1681,6 +1680,173 @@ static inline void reorder_by_three_exchanges_with_segment_flipping(std::vector<
} }
} }
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::DontAlign> 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<double, size_t> 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<FlipEdge> &edges)
{
if (edges.size() < 3) {
reorder_by_two_exchanges_with_segment_flipping(edges);
return;
}
std::vector<ConnectionCost> connections(edges.size());
std::vector<FlipEdge> edges_tmp(edges);
std::vector<std::pair<double, size_t>> connection_lengths(edges.size() - 1, std::pair<double, size_t>(0., 0));
std::vector<char> 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<double, size_t> &l, const std::pair<double, size_t> &r) { return l.first > r.first; });
std::fill(connection_tried.begin(), connection_tried.end(), false);
size_t crossover1_pos_final = std::numeric_limits<size_t>::max();
size_t crossover2_pos_final = std::numeric_limits<size_t>::max();
size_t crossover3_pos_final = std::numeric_limits<size_t>::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<double>::max());
#endif /* NDEBUG */
for (const std::pair<double, size_t> &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<size_t>::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<double, size_t> 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. // 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) 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 #if 1
reorder_by_two_exchanges_with_segment_flipping(edges); reorder_by_two_exchanges_with_segment_flipping(edges);
#else #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 #endif
Polylines out; Polylines out;
out.reserve(polylines.size()); out.reserve(polylines.size());