Fix of a typo in KDTreeIndirect.
Improvement of the infill path planning. Regression fix of Gyroid infill crashes. Some unit tests for elephant foot and path planning.
This commit is contained in:
bubnikv
parent
ae887d5833
commit
dd59945098
9 changed files with 443 additions and 145 deletions
|
|
@ -1205,7 +1205,7 @@ ExPolygons variable_offset_inner_ex(const ExPolygon &expoly, const std::vector<s
|
|||
{
|
||||
#ifndef NDEBUG
|
||||
// Verify that the deltas are all non positive.
|
||||
for (const std::vector<float>& ds : deltas)
|
||||
for (const std::vector<float>& ds : deltas)
|
||||
for (float delta : ds)
|
||||
assert(delta <= 0.);
|
||||
assert(expoly.holes.size() + 1 == deltas.size());
|
||||
|
|
|
|||
|
|
@ -60,9 +60,9 @@ std::vector<float> 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<float> contour_distance(const EdgeGrid::Grid &grid, const size_t idx
|
|||
dir *= t;
|
||||
this->pt_end = (this->pt + dir).cast<coord_t>();
|
||||
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<std::vector<float>> 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<ResampledPoint> resampled_point_parameters;
|
||||
poly.points = resample_polygon(poly.points, resample_interval, resampled_point_parameters);
|
||||
std::vector<float> 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<double>(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<std::vector<float>> 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<ResampledPoint> resampled_point_parameters;
|
||||
poly.points = resample_polygon(poly.points, resample_interval, resampled_point_parameters);
|
||||
std::vector<float> 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<double>(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)
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
};
|
||||
|
|
|
|||
|
|
@ -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<ContourPointData> &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<size_t> 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<ConnectionCost> 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<size_t, size_t> *cp1 = &map_infill_end_point_to_boundary[(idx_chain - 1) * 2 + 1];
|
||||
const std::pair<size_t, size_t> *cp2 = &map_infill_end_point_to_boundary[idx_chain * 2];
|
||||
const Points &contour = boundary[cp1->first];
|
||||
std::vector<ContourPointData> &contour_data = boundary_data[cp1->first];
|
||||
bool valid = false;
|
||||
bool reversed = false;
|
||||
const std::vector<ContourPointData> &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<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];
|
||||
assert(cp1->first == cp2->first);
|
||||
std::vector<ContourPointData> &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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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<size_t> &input, size_t node, int dimension, int left, int right)
|
||||
void build_recursive(std::vector<size_t> &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 <left, right> at k using QuickSelect method.
|
||||
// Partition the input m_nodes <left, right> at "k" and "dimension" using the QuickSelect method:
|
||||
// https://en.wikipedia.org/wiki/Quickselect
|
||||
void partition_input(std::vector<size_t> &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<size_t> &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 <left, right> 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) {
|
||||
|
|
|
|||
|
|
@ -237,11 +237,19 @@ std::vector<std::pair<size_t, bool>> 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<std::pair<size_t, bool>> 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<size_t> 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<double>().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<size_t>::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<Connection> &connections) const
|
||||
{ return ((this + 1)->start_point(polylines, connections) - this->end_point(polylines, connections)).squaredNorm(); }
|
||||
double norm(const Polylines &polylines, const std::vector<Connection> &connections) const
|
||||
{ return sqrt(this->squaredNorm(polylines, connections)); }
|
||||
double squaredNorm(const Polylines &polylines, const std::vector<Connection> &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<Connection> &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<Connection> &connections, bool flip = false) const
|
||||
{ const Polyline &pl = polylines[this - connections.data()]; return ((this->flipped == flip) ? pl.points.front() : pl.points.back()).cast<double>(); }
|
||||
Vec2d end_point(const Polylines &polylines, const std::vector<Connection> &connections, bool flip = false) const
|
||||
{ const Polyline &pl = polylines[this - connections.data()]; return ((this->flipped == flip) ? pl.points.back() : pl.points.front()).cast<double>(); }
|
||||
|
||||
bool in_queue() const { return this->heap_idx != std::numeric_limits<size_t>::max(); }
|
||||
void flip() { this->flipped = ! this->flipped; }
|
||||
};
|
||||
std::vector<Connection> 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<double>(connections[i - 1].norm(polylines, connections)),
|
||||
int(connections[i - 1].flipped), int(connections[i].flipped),
|
||||
unscale<double>(connections[i - 1].norm(polylines, connections, true, false)),
|
||||
unscale<double>(connections[i - 1].norm(polylines, connections, false, true)),
|
||||
unscale<double>(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 *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<size_t>::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<std::pair<size_t, bool>> ordered = chain_segments_greedy<Point, decltype(segment_end_point)>(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<std::pair<size_t, bool>> ordered = chain_segments_greedy<Point, decltype(segment_end_point)>(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<class T> static inline T chain_path_items(const Points &points, const T &items)
|
||||
|
|
|
|||
|
|
@ -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") {
|
||||
|
|
|
|||
|
|
@ -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<Point> 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<Points::size_type> indices = chain_points(points);
|
||||
for (Points::size_type i = 0; i + 1 < indices.size(); ++ i) {
|
||||
double dist = (points.at(indices.at(i)).cast<double>() - points.at(indices.at(i+1)).cast<double>()).norm();
|
||||
REQUIRE(std::abs(dist-26) <= EPSILON);
|
||||
}
|
||||
SCENARIO("Path chaining", "[Geometry]") {
|
||||
GIVEN("A path") {
|
||||
std::vector<Point> 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<Points::size_type> indices = chain_points(points);
|
||||
for (Points::size_type i = 0; i + 1 < indices.size(); ++ i) {
|
||||
double dist = (points.at(indices.at(i)).cast<double>() - points.at(indices.at(i+1)).cast<double>()).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]"){
|
||||
|
|
|
|||
Loading…
Reference in a new issue