From 5e037c766f11060396903663e1e7d19d10a7391c Mon Sep 17 00:00:00 2001 From: supermerill Date: Mon, 1 Oct 2018 19:16:04 +0200 Subject: [PATCH] forgot medial_axis.hpp/cpp files --- xs/src/libslic3r/MedialAxis.cpp | 1413 +++++++++++++++++++++++++++++++ xs/src/libslic3r/MedialAxis.hpp | 63 ++ 2 files changed, 1476 insertions(+) create mode 100644 xs/src/libslic3r/MedialAxis.cpp create mode 100644 xs/src/libslic3r/MedialAxis.hpp diff --git a/xs/src/libslic3r/MedialAxis.cpp b/xs/src/libslic3r/MedialAxis.cpp new file mode 100644 index 00000000000..27b77e197a9 --- /dev/null +++ b/xs/src/libslic3r/MedialAxis.cpp @@ -0,0 +1,1413 @@ +#include "BoundingBox.hpp" +#include "ExPolygon.hpp" +#include "Geometry.hpp" +#include "Polygon.hpp" +#include "Line.hpp" +#include "ClipperUtils.hpp" +#include "SVG.hpp" +#include "polypartition.h" +#include "poly2tri/poly2tri.h" +#include +#include +#include + +namespace Slic3r { + +int MedialAxis::id = 0; + +void +MedialAxis::build(Polylines* polylines) +{ + ThickPolylines tp; + this->build(&tp); + polylines->insert(polylines->end(), tp.begin(), tp.end()); +} + +void +MedialAxis::polyline_from_voronoi(const Lines& voronoi_edges, ThickPolylines* polylines) +{ + this->lines = voronoi_edges; + construct_voronoi(lines.begin(), lines.end(), &this->vd); + + /* + // DEBUG: dump all Voronoi edges + { + for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) { + if (edge->is_infinite()) continue; + + ThickPolyline polyline; + polyline.points.push_back(Point( edge->vertex0()->x(), edge->vertex0()->y() )); + polyline.points.push_back(Point( edge->vertex1()->x(), edge->vertex1()->y() )); + polylines->push_back(polyline); + } + return; + } + */ + + typedef const VD::vertex_type vert_t; + typedef const VD::edge_type edge_t; + + // collect valid edges (i.e. prune those not belonging to MAT) + // note: this keeps twins, so it inserts twice the number of the valid edges + this->valid_edges.clear(); + { + std::set seen_edges; + for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) { + // if we only process segments representing closed loops, none if the + // infinite edges (if any) would be part of our MAT anyway + if (edge->is_secondary() || edge->is_infinite()) continue; + + // don't re-validate twins + if (seen_edges.find(&*edge) != seen_edges.end()) continue; // TODO: is this needed? + seen_edges.insert(&*edge); + seen_edges.insert(edge->twin()); + + if (!this->validate_edge(&*edge)) continue; + this->valid_edges.insert(&*edge); + this->valid_edges.insert(edge->twin()); + } + } + this->edges = this->valid_edges; + + // iterate through the valid edges to build polylines + while (!this->edges.empty()) { + const edge_t* edge = *this->edges.begin(); + + // start a polyline + ThickPolyline polyline; + polyline.points.push_back(Point( edge->vertex0()->x(), edge->vertex0()->y() )); + polyline.points.push_back(Point( edge->vertex1()->x(), edge->vertex1()->y() )); + polyline.width.push_back(this->thickness[edge].first); + polyline.width.push_back(this->thickness[edge].second); + + // remove this edge and its twin from the available edges + (void)this->edges.erase(edge); + (void)this->edges.erase(edge->twin()); + + // get next points + this->process_edge_neighbors(edge, &polyline); + + // get previous points + { + ThickPolyline rpolyline; + this->process_edge_neighbors(edge->twin(), &rpolyline); + polyline.points.insert(polyline.points.begin(), rpolyline.points.rbegin(), rpolyline.points.rend()); + polyline.width.insert(polyline.width.begin(), rpolyline.width.rbegin(), rpolyline.width.rend()); + polyline.endpoints.first = rpolyline.endpoints.second; + } + + assert(polyline.width.size() == polyline.points.size()); + + // prevent loop endpoints from being extended + if (polyline.first_point().coincides_with(polyline.last_point())) { + polyline.endpoints.first = false; + polyline.endpoints.second = false; + } + + // append polyline to result + polylines->push_back(polyline); + } + + #ifdef SLIC3R_DEBUG + { + static int iRun = 0; + dump_voronoi_to_svg(this->lines, this->vd, polylines, debug_out_path("MedialAxis-%d.svg", iRun ++).c_str()); + printf("Thick lines: "); + for (ThickPolylines::const_iterator it = polylines->begin(); it != polylines->end(); ++ it) { + ThickLines lines = it->thicklines(); + for (ThickLines::const_iterator it2 = lines.begin(); it2 != lines.end(); ++ it2) { + printf("%f,%f ", it2->a_width, it2->b_width); + } + } + printf("\n"); + } + #endif /* SLIC3R_DEBUG */ +} + +void +MedialAxis::process_edge_neighbors(const VD::edge_type* edge, ThickPolyline* polyline) +{ + while (true) { + // Since rot_next() works on the edge starting point but we want + // to find neighbors on the ending point, we just swap edge with + // its twin. + const VD::edge_type* twin = edge->twin(); + + // count neighbors for this edge + std::vector neighbors; + for (const VD::edge_type* neighbor = twin->rot_next(); neighbor != twin; + neighbor = neighbor->rot_next()) { + if (this->valid_edges.count(neighbor) > 0) neighbors.push_back(neighbor); + } + + // if we have a single neighbor then we can continue recursively + if (neighbors.size() == 1) { + const VD::edge_type* neighbor = neighbors.front(); + + // break if this is a closed loop + if (this->edges.count(neighbor) == 0) return; + + Point new_point(neighbor->vertex1()->x(), neighbor->vertex1()->y()); + polyline->points.push_back(new_point); + polyline->width.push_back(this->thickness[neighbor].second); + + (void)this->edges.erase(neighbor); + (void)this->edges.erase(neighbor->twin()); + edge = neighbor; + } else if (neighbors.size() == 0) { + polyline->endpoints.second = true; + return; + } else { + // T-shaped or star-shaped joint + return; + } + } +} + +bool +MedialAxis::validate_edge(const VD::edge_type* edge) +{ + // prevent overflows and detect almost-infinite edges + if (std::abs(edge->vertex0()->x()) > double(CLIPPER_MAX_COORD_UNSCALED) || + std::abs(edge->vertex0()->y()) > double(CLIPPER_MAX_COORD_UNSCALED) || + std::abs(edge->vertex1()->x()) > double(CLIPPER_MAX_COORD_UNSCALED) || + std::abs(edge->vertex1()->y()) > double(CLIPPER_MAX_COORD_UNSCALED)) + return false; + + // construct the line representing this edge of the Voronoi diagram + const Line line( + Point( edge->vertex0()->x(), edge->vertex0()->y() ), + Point( edge->vertex1()->x(), edge->vertex1()->y() ) + ); + + // discard edge if it lies outside the supplied shape + // this could maybe be optimized (checking inclusion of the endpoints + // might give false positives as they might belong to the contour itself) + if (line.a.coincides_with(line.b)) { + // in this case, contains(line) returns a false positive + if (!this->expolygon.contains(line.a)) return false; + } else { + if (!this->expolygon.contains(line)) return false; + } + + // retrieve the original line segments which generated the edge we're checking + const VD::cell_type* cell_l = edge->cell(); + const VD::cell_type* cell_r = edge->twin()->cell(); + const Line &segment_l = this->retrieve_segment(cell_l); + const Line &segment_r = this->retrieve_segment(cell_r); + + + //SVG svg("edge.svg"); + //svg.draw(this->expolygon.expolygon); + //svg.draw(line); + //svg.draw(segment_l, "red"); + //svg.draw(segment_r, "blue"); + //svg.Close(); + // + + /* Calculate thickness of the cross-section at both the endpoints of this edge. + Our Voronoi edge is part of a CCW sequence going around its Voronoi cell + located on the left side. (segment_l). + This edge's twin goes around segment_r. Thus, segment_r is + oriented in the same direction as our main edge, and segment_l is oriented + in the same direction as our twin edge. + We used to only consider the (half-)distances to segment_r, and that works + whenever segment_l and segment_r are almost specular and facing. However, + at curves they are staggered and they only face for a very little length + (our very short edge represents such visibility). + Both w0 and w1 can be calculated either towards cell_l or cell_r with equal + results by Voronoi definition. + When cell_l or cell_r don't refer to the segment but only to an endpoint, we + calculate the distance to that endpoint instead. */ + + coordf_t w0 = cell_r->contains_segment() + ? line.a.distance_to(segment_r)*2 + : line.a.distance_to(this->retrieve_endpoint(cell_r))*2; + + coordf_t w1 = cell_l->contains_segment() + ? line.b.distance_to(segment_l)*2 + : line.b.distance_to(this->retrieve_endpoint(cell_l))*2; + + //don't remove the line that goes to the intersection of the contour + // we use them to create nicer thin wall lines + //if (cell_l->contains_segment() && cell_r->contains_segment()) { + // // calculate the relative angle between the two boundary segments + // double angle = fabs(segment_r.orientation() - segment_l.orientation()); + // if (angle > PI) angle = 2*PI - angle; + // assert(angle >= 0 && angle <= PI); + // + // // fabs(angle) ranges from 0 (collinear, same direction) to PI (collinear, opposite direction) + // // we're interested only in segments close to the second case (facing segments) + // // so we allow some tolerance. + // // this filter ensures that we're dealing with a narrow/oriented area (longer than thick) + // // we don't run it on edges not generated by two segments (thus generated by one segment + // // and the endpoint of another segment), since their orientation would not be meaningful + // if (PI - angle > PI/8) { + // // angle is not narrow enough + // + // // only apply this filter to segments that are not too short otherwise their + // // angle could possibly be not meaningful + // if (w0 < SCALED_EPSILON || w1 < SCALED_EPSILON || line.length() >= this->min_width) + // return false; + // } + //} else { + // if (w0 < SCALED_EPSILON || w1 < SCALED_EPSILON) + // return false; + //} + + // don't do that before we try to fusion them + //if (w0 < this->min_width && w1 < this->min_width) + // return false; + // + + //shouldn't occur if perimeter_generator is well made + if (w0 > this->max_width && w1 > this->max_width) + return false; + + this->thickness[edge] = std::make_pair(w0, w1); + this->thickness[edge->twin()] = std::make_pair(w1, w0); + + return true; +} + +const Line& +MedialAxis::retrieve_segment(const VD::cell_type* cell) const +{ + return lines[cell->source_index()]; +} + +const Point& +MedialAxis::retrieve_endpoint(const VD::cell_type* cell) const +{ + const Line& line = this->retrieve_segment(cell); + if (cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) { + return line.a; + } else { + return line.b; + } +} + + +/// remove point that are at SCALED_EPSILON * 2 distance. +void +remove_point_too_near(ThickPolyline* to_reduce) +{ + const coord_t smallest = SCALED_EPSILON * 2; + size_t id = 1; + while (id < to_reduce->points.size() - 1) { + size_t newdist = min(to_reduce->points[id].distance_to(to_reduce->points[id - 1]) + , to_reduce->points[id].distance_to(to_reduce->points[id + 1])); + if (newdist < smallest) { + to_reduce->points.erase(to_reduce->points.begin() + id); + to_reduce->width.erase(to_reduce->width.begin() + id); + newdist = to_reduce->points[id].distance_to(to_reduce->points[id - 1]); + } + //go to next one + //if you removed a point, it check if the next one isn't too near from the previous one. + // if not, it byepass it. + if (newdist > smallest) { + ++id; + } + } +} + +/// add points from pattern to to_modify at the same % of the length +/// so not add if an other point is present at the correct position +void +add_point_same_percent(ThickPolyline* pattern, ThickPolyline* to_modify) +{ + const double to_modify_length = to_modify->length(); + const double percent_epsilon = SCALED_EPSILON / to_modify_length; + const double pattern_length = pattern->length(); + + double percent_length = 0; + for (size_t idx_point = 1; idx_point < pattern->points.size() - 1; ++idx_point) { + percent_length += pattern->points[idx_point-1].distance_to(pattern->points[idx_point]) / pattern_length; + //find position + size_t idx_other = 1; + double percent_length_other_before = 0; + double percent_length_other = 0; + while (idx_other < to_modify->points.size()) { + percent_length_other_before = percent_length_other; + percent_length_other += to_modify->points[idx_other-1].distance_to(to_modify->points[idx_other]) + / to_modify_length; + if (percent_length_other > percent_length - percent_epsilon) { + //if higher (we have gone over it) + break; + } + ++idx_other; + } + if (percent_length_other > percent_length + percent_epsilon) { + //insert a new point before the position + double percent_dist = (percent_length - percent_length_other_before) / (percent_length_other - percent_length_other_before); + coordf_t new_width = to_modify->width[idx_other - 1] * (1 - percent_dist); + new_width += to_modify->width[idx_other] * (percent_dist); + Point new_point; + new_point.x = (coord_t)((double)(to_modify->points[idx_other - 1].x) * (1 - percent_dist)); + new_point.x += (coord_t)((double)(to_modify->points[idx_other].x) * (percent_dist)); + new_point.y = (coord_t)((double)(to_modify->points[idx_other - 1].y) * (1 - percent_dist)); + new_point.y += (coord_t)((double)(to_modify->points[idx_other].y) * (percent_dist)); + to_modify->width.insert(to_modify->width.begin() + idx_other, new_width); + to_modify->points.insert(to_modify->points.begin() + idx_other, new_point); + } + } +} + +/// find the nearest angle in the contour (or 2 nearest if it's difficult to choose) +/// return 1 for an angle of 90° and 0 for an angle of 0° or 180° +/// find the nearest angle in the contour (or 2 nearest if it's difficult to choose) +/// return 1 for an angle of 90° and 0 for an angle of 0° or 180° +double +get_coeff_from_angle_countour(Point &point, const ExPolygon &contour, coord_t min_dist_between_point) { + double nearestDist = point.distance_to(contour.contour.points.front()); + Point nearest = contour.contour.points.front(); + size_t id_nearest = 0; + double nearDist = nearestDist; + Point near = nearest; + size_t id_near = 0; + for (size_t id_point = 1; id_point < contour.contour.points.size(); ++id_point) { + if (nearestDist > point.distance_to(contour.contour.points[id_point])) { + //update near + id_near = id_nearest; + near = nearest; + nearDist = nearestDist; + //update nearest + nearestDist = point.distance_to(contour.contour.points[id_point]); + nearest = contour.contour.points[id_point]; + id_nearest = id_point; + } + } + double angle = 0; + size_t id_before = id_nearest == 0 ? contour.contour.points.size() - 1 : id_nearest - 1; + Point point_before = id_nearest == 0 ? contour.contour.points.back() : contour.contour.points[id_nearest - 1]; + //Search one point far enough to be relevant + while (nearest.distance_to(point_before) < min_dist_between_point) { + point_before = id_before == 0 ? contour.contour.points.back() : contour.contour.points[id_before - 1]; + id_before = id_before == 0 ? contour.contour.points.size() - 1 : id_before - 1; + //don't loop + if (id_before == id_nearest) { + id_before = id_nearest == 0 ? contour.contour.points.size() - 1 : id_nearest - 1; + point_before = id_nearest == 0 ? contour.contour.points.back() : contour.contour.points[id_nearest - 1]; + break; + } + } + size_t id_after = id_nearest == contour.contour.points.size() - 1 ? 0 : id_nearest + 1; + Point point_after = id_nearest == contour.contour.points.size() - 1 ? contour.contour.points.front() : contour.contour.points[id_nearest + 1]; + //Search one point far enough to be relevant + while (nearest.distance_to(point_after) < min_dist_between_point) { + point_after = id_after == contour.contour.points.size() - 1 ? contour.contour.points.front() : contour.contour.points[id_after + 1]; + id_after = id_after == contour.contour.points.size() - 1 ? 0 : id_after + 1; + //don't loop + if (id_after == id_nearest) { + id_after = id_nearest == contour.contour.points.size() - 1 ? 0 : id_nearest + 1; + point_after = id_nearest == contour.contour.points.size() - 1 ? contour.contour.points.front() : contour.contour.points[id_nearest + 1]; + break; + } + } + //compute angle + angle = nearest.ccw_angle(point_before, point_after); + if (angle >= PI) angle = 2 * PI - angle; // smaller angle + //compute the diff from 90° + angle = abs(angle - PI / 2); + if (near.coincides_with(nearest) && max(nearestDist, nearDist) + SCALED_EPSILON < nearest.distance_to(near)) { + //not only nearest + Point point_before = id_near == 0 ? contour.contour.points.back() : contour.contour.points[id_near - 1]; + Point point_after = id_near == contour.contour.points.size() - 1 ? contour.contour.points.front() : contour.contour.points[id_near + 1]; + double angle2 = min(nearest.ccw_angle(point_before, point_after), nearest.ccw_angle(point_after, point_before)); + angle2 = abs(angle - PI / 2); + angle = (angle + angle2) / 2; + } + + return 1 - (angle / (PI / 2)); +} + +double +dot(Line l1, Line l2) +{ + Vectorf v_1 = normalize(Vectorf(l1.b.x - l1.a.x, l1.b.y - l1.a.y)); + Vectorf v_2 = normalize(Vectorf(l2.b.x - l2.a.x, l2.b.y - l2.a.y)); + return v_1.x*v_2.x + v_1.y*v_2.y; +} + +void +MedialAxis::fusion_curve(ThickPolylines &pp) +{ + //fusion Y with only 1 '0' value => the "0" branch "pull" the cross-point + bool changes = false; + for (size_t i = 0; i < pp.size(); ++i) { + ThickPolyline& polyline = pp[i]; + // only consider 2-point polyline with endpoint + if (polyline.points.size() != 2) continue; + if (polyline.endpoints.first) polyline.reverse(); + else if (!polyline.endpoints.second) continue; + if (polyline.width.back() > EPSILON) continue; + + //check my length is small + coord_t length = (coord_t)polyline.length(); + if (length > max_width) continue; + + size_t closest_point_idx = this->expolygon.contour.closest_point_index(polyline.points.back()); + + //check the 0-wodth point is on the contour. + if (closest_point_idx == (size_t)-1) continue; + + size_t prev_idx = closest_point_idx == 0 ? this->expolygon.contour.points.size() - 1 : closest_point_idx - 1; + size_t next_idx = closest_point_idx == this->expolygon.contour.points.size() - 1 ? 0 : closest_point_idx + 1; + double mindot = 1; + mindot = min(mindot, abs(dot(Line(polyline.points[polyline.points.size() - 1], polyline.points[polyline.points.size() - 2]), + (Line(this->expolygon.contour.points[closest_point_idx], this->expolygon.contour.points[prev_idx]))))); + mindot = min(mindot, abs(dot(Line(polyline.points[polyline.points.size() - 1], polyline.points[polyline.points.size() - 2]), + (Line(this->expolygon.contour.points[closest_point_idx], this->expolygon.contour.points[next_idx]))))); + + //compute angle + double coeff_contour_angle = this->expolygon.contour.points[closest_point_idx].ccw_angle(this->expolygon.contour.points[prev_idx], this->expolygon.contour.points[next_idx]); + if (coeff_contour_angle >= PI) coeff_contour_angle = 2 * PI - coeff_contour_angle; // smaller angle + //compute the diff from 90° + coeff_contour_angle = abs(coeff_contour_angle - PI / 2); + + + // look if other end is a cross point with almost 90° angle + double sum_dot = 0; + double min_dot = 0; + // look if other end is a cross point with multiple other branch + vector crosspoint; + for (size_t j = 0; j < pp.size(); ++j) { + if (j == i) continue; + ThickPolyline& other = pp[j]; + if (polyline.first_point().coincides_with(other.last_point())) { + other.reverse(); + crosspoint.push_back(j); + double dot_temp = dot(Line(polyline.points[0], polyline.points[1]), (Line(other.points[0], other.points[1]))); + min_dot = min(min_dot, abs(dot_temp)); + sum_dot += dot_temp; + } else if (polyline.first_point().coincides_with(other.first_point())) { + crosspoint.push_back(j); + double dot_temp = dot(Line(polyline.points[0], polyline.points[1]), (Line(other.points[0], other.points[1]))); + min_dot = min(min_dot, abs(dot_temp)); + sum_dot += dot_temp; + } + } + //only consider very shallow angle for contour + if (mindot > 0.15 && + (1 - (coeff_contour_angle / (PI / 2))) > 0.2) continue; + + //check if it's a line that we can pull + if (crosspoint.size() != 2) continue; + if (sum_dot > 0.2) continue; + if (min_dot > 0.5) continue; + + //don't pull, it distords the line if there are too many points. + //// pull it a bit, depends on my size, the dot?, and the coeff at my 0-end (~14% for a square, almost 0 for a gentle curve) + //coord_t length_pull = polyline.length(); + //length_pull *= 0.144 * get_coeff_from_angle_countour(polyline.points.back(), this->expolygon, min(min_width, polyline.length() / 2)); + + ////compute dir + //Vectorf pull_direction(polyline.points[1].x - polyline.points[0].x, polyline.points[1].y - polyline.points[0].y); + //pull_direction = normalize(pull_direction); + //pull_direction.x *= length_pull; + //pull_direction.y *= length_pull; + + ////pull the points + //Point &p1 = pp[crosspoint[0]].points[0]; + //p1.x = p1.x + (coord_t)pull_direction.x; + //p1.y = p1.y + (coord_t)pull_direction.y; + + //Point &p2 = pp[crosspoint[1]].points[0]; + //p2.x = p2.x + (coord_t)pull_direction.x; + //p2.y = p2.y + (coord_t)pull_direction.y; + + //delete the now unused polyline + pp.erase(pp.begin() + i); + --i; + changes = true; + } + if (changes) { + concatThickPolylines(pp); + ///reorder, in case of change + std::sort(pp.begin(), pp.end(), [](const ThickPolyline & a, const ThickPolyline & b) { return a.length() < b.length(); }); + } +} + +void +MedialAxis::fusion_corners(ThickPolylines &pp) +{ + + //fusion Y with only 1 '0' value => the "0" branch "pull" the cross-point + bool changes = false; + for (size_t i = 0; i < pp.size(); ++i) { + ThickPolyline& polyline = pp[i]; + // only consider polyline with 0-end + if (polyline.points.size() != 2) continue; + if (polyline.endpoints.first) polyline.reverse(); + else if (!polyline.endpoints.second) continue; + if (polyline.width.back() > 0) continue; + + //check my length is small + coord_t length = polyline.length(); + if (length > max_width) continue; + + // look if other end is a cross point with multiple other branch + vector crosspoint; + for (size_t j = 0; j < pp.size(); ++j) { + if (j == i) continue; + ThickPolyline& other = pp[j]; + if (polyline.first_point().coincides_with(other.last_point())) { + other.reverse(); + crosspoint.push_back(j); + } else if (polyline.first_point().coincides_with(other.first_point())) { + crosspoint.push_back(j); + } + } + //check if it's a line that we can pull + if (crosspoint.size() != 2) continue; + + // check if i am at the external side of a curve + double angle1 = polyline.points[0].ccw_angle(polyline.points[1], pp[crosspoint[0]].points[1]); + if (angle1 >= PI) angle1 = 2 * PI - angle1; // smaller angle + double angle2 = polyline.points[0].ccw_angle(polyline.points[1], pp[crosspoint[1]].points[1]); + if (angle2 >= PI) angle2 = 2 * PI - angle2; // smaller angle + if (angle1 + angle2 < PI) continue; + + //check if is smaller or the other ones are not endpoits + if (pp[crosspoint[0]].endpoints.second && length > pp[crosspoint[0]].length()) continue; + if (pp[crosspoint[1]].endpoints.second && length > pp[crosspoint[1]].length()) continue; + + //FIXME: also pull (a bit less) points that are near to this one. + // if true, pull it a bit, depends on my size, the dot?, and the coeff at my 0-end (~14% for a square, almost 0 for a gentle curve) + coord_t length_pull = polyline.length(); + length_pull *= 0.144 * get_coeff_from_angle_countour(polyline.points.back(), this->expolygon, min(min_width, polyline.length() / 2)); + + //compute dir + Vectorf pull_direction(polyline.points[1].x - polyline.points[0].x, polyline.points[1].y - polyline.points[0].y); + pull_direction = normalize(pull_direction); + pull_direction.x *= length_pull; + pull_direction.y *= length_pull; + + //pull the points + Point &p1 = pp[crosspoint[0]].points[0]; + p1.x = p1.x + pull_direction.x; + p1.y = p1.y + pull_direction.y; + + Point &p2 = pp[crosspoint[1]].points[0]; + p2.x = p2.x + pull_direction.x; + p2.y = p2.y + pull_direction.y; + + //delete the now unused polyline + pp.erase(pp.begin() + i); + --i; + changes = true; + } + if (changes) { + concatThickPolylines(pp); + ///reorder, in case of change + std::sort(pp.begin(), pp.end(), [](const ThickPolyline & a, const ThickPolyline & b) { return a.length() < b.length(); }); + } +} + + +void +MedialAxis::extends_line(ThickPolyline& polyline, const ExPolygons& anchors, const coord_t join_width) +{ + // extend initial and final segments of each polyline if they're actual endpoints + // We assign new endpoints to temporary variables because in case of a single-line + // polyline, after we extend the start point it will be caught by the intersection() + // call, so we keep the inner point until we perform the second intersection() as well + if (polyline.endpoints.second && !bounds.has_boundary_point(polyline.points.back())) { + Line line(*(polyline.points.end() - 2), polyline.points.back()); + + // prevent the line from touching on the other side, otherwise intersection() might return that solution + if (polyline.points.size() == 2) line.a = line.midpoint(); + + line.extend_end(max_width); + Point new_back; + if (this->expolygon.contour.has_boundary_point(polyline.points.back())) { + new_back = polyline.points.back(); + } else { + (void)this->expolygon.contour.first_intersection(line, &new_back); + polyline.points.push_back(new_back); + polyline.width.push_back(polyline.width.back()); + } + Point new_bound; + (void)bounds.contour.first_intersection(line, &new_bound); + /* if (new_bound.coincides_with_epsilon(new_back)) { + return; + }*/ + // find anchor + Point best_anchor; + double shortest_dist = max_width; + for (const ExPolygon& a : anchors) { + Point p_maybe_inside = a.contour.centroid(); + double test_dist = new_bound.distance_to(p_maybe_inside) + new_back.distance_to(p_maybe_inside); + //if (test_dist < max_width / 2 && (test_dist < shortest_dist || shortest_dist < 0)) { + double angle_test = new_back.ccw_angle(p_maybe_inside, line.a); + if (angle_test > PI) angle_test = 2 * PI - angle_test; + if (test_dist < max_width && test_dist PI / 2) { + shortest_dist = test_dist; + best_anchor = p_maybe_inside; + } + } + if (best_anchor.x != 0 && best_anchor.y != 0) { + Point p_obj = best_anchor + new_bound; + p_obj.x /= 2; + p_obj.y /= 2; + Line l2 = Line(new_back, p_obj); + l2.extend_end(max_width); + (void)bounds.contour.first_intersection(l2, &new_bound); + } + if (new_bound.coincides_with_epsilon(new_back)) { + return; + } + polyline.points.push_back(new_bound); + //polyline.width.push_back(join_width); + //it thickens the line a bit too early, imo + polyline.width.push_back(polyline.width.back()); + } +} + +void +MedialAxis::main_fusion(ThickPolylines& pp) +{ + //int idf = 0; + //reoder pp by length (ascending) It's really important to do that to avoid building the line from the width insteand of the length + std::sort(pp.begin(), pp.end(), [](const ThickPolyline & a, const ThickPolyline & b) { + bool ahas0 = a.width.front() == 0 || a.width.back() == 0; + bool bhas0 = b.width.front() == 0 || b.width.back() == 0; + if (ahas0 && !bhas0) return true; + if (!ahas0 && bhas0) return false; + return a.length() < b.length(); + }); + + bool changes = true; + map coeff_angle_cache; + while (changes) { + changes = false; + for (size_t i = 0; i < pp.size(); ++i) { + ThickPolyline& polyline = pp[i]; + + //simple check to see if i can be fusionned + if (!polyline.endpoints.first && !polyline.endpoints.second) continue; + + + ThickPolyline* best_candidate = nullptr; + float best_dot = -1; + size_t best_idx = 0; + double dot_poly_branch = 0; + double dot_candidate_branch = 0; + + // find another polyline starting here + for (size_t j = i + 1; j < pp.size(); ++j) { + ThickPolyline& other = pp[j]; + if (polyline.last_point().coincides_with(other.last_point())) { + polyline.reverse(); + other.reverse(); + } else if (polyline.first_point().coincides_with(other.last_point())) { + other.reverse(); + } else if (polyline.first_point().coincides_with(other.first_point())) { + } else if (polyline.last_point().coincides_with(other.first_point())) { + polyline.reverse(); + } else { + continue; + } + //std::cout << " try : " << i << ":" << j << " : " << + // (polyline.points.size() < 2 && other.points.size() < 2) << + // (!polyline.endpoints.second || !other.endpoints.second) << + // ((polyline.points.back().distance_to(other.points.back()) + // + (polyline.width.back() + other.width.back()) / 4) + // > max_width*1.05) << + // (abs(polyline.length() - other.length()) > max_width) << "\n"; + + //// mergeable tests + if (polyline.points.size() < 2 && other.points.size() < 2) continue; + if (!polyline.endpoints.second || !other.endpoints.second) continue; + // test if the new width will not be too big if a fusion occur + //note that this isn't the real calcul. It's just to avoid merging lines too far apart. + if ( + ((polyline.points.back().distance_to(other.points.back()) + + (polyline.width.back() + other.width.back()) / 4) + > max_width*1.05)) + continue; + // test if the lines are not too different in length. + if (abs(polyline.length() - other.length()) > max_width) continue; + + + //test if we don't merge with something too different and without any relevance. + double coeffSizePolyI = 1; + if (polyline.width.back() == 0) { + coeffSizePolyI = 0.1 + 0.9*get_coeff_from_angle_countour(polyline.points.back(), this->expolygon, min(min_width, polyline.length() / 2)); + } + double coeffSizeOtherJ = 1; + if (other.width.back() == 0) { + coeffSizeOtherJ = 0.1 + 0.9*get_coeff_from_angle_countour(other.points.back(), this->expolygon, min(min_width, polyline.length() / 2)); + } + //std::cout << " try2 : " << i << ":" << j << " : " + // << (abs(polyline.length()*coeffSizePolyI - other.length()*coeffSizeOtherJ) > max_width / 2) + // << (abs(polyline.length()*coeffSizePolyI - other.length()*coeffSizeOtherJ) > max_width) + // << "\n"; + if (abs(polyline.length()*coeffSizePolyI - other.length()*coeffSizeOtherJ) > max_width / 2) continue; + + + //compute angle to see if it's better than previous ones (straighter = better). + //we need to add how strait we are from our main. + float test_dot = dot(polyline.lines().front(), other.lines().front()); + + // Get the branch/line in wich we may merge, if possible + // with that, we can decide what is important, and how we can merge that. + // angle_poly - angle_candi =90° => one is useless + // both angle are equal => both are useful with same strength + // ex: Y => | both are useful to crete a nice line + // ex2: TTTTT => ----- these 90° useless lines should be discarded + bool find_main_branch = false; + size_t biggest_main_branch_id = 0; + coord_t biggest_main_branch_length = 0; + for (size_t k = 0; k < pp.size(); ++k) { + //std::cout << "try to find main : " << k << " ? " << i << " " << j << " "; + if (k == i | k == j) continue; + ThickPolyline& main = pp[k]; + if (polyline.first_point().coincides_with(main.last_point())) { + main.reverse(); + if (!main.endpoints.second) + find_main_branch = true; + else if (biggest_main_branch_length < main.length()) { + biggest_main_branch_id = k; + biggest_main_branch_length = main.length(); + } + } else if (polyline.first_point().coincides_with(main.first_point())) { + if (!main.endpoints.second) + find_main_branch = true; + else if (biggest_main_branch_length < main.length()) { + biggest_main_branch_id = k; + biggest_main_branch_length = main.length(); + } + } + if (find_main_branch) { + //use this variable to store the good index and break to compute it + biggest_main_branch_id = k; + break; + } + } + if (!find_main_branch && biggest_main_branch_length == 0) { + // nothing -> it's impossible! + dot_poly_branch = 0.707; + dot_candidate_branch = 0.707; + //std::cout << "no main branch... impossible!!\n"; + } else if (!find_main_branch && ( + (pp[biggest_main_branch_id].length() < polyline.length() && (polyline.width.back() != 0 || pp[biggest_main_branch_id].width.back() ==0)) + || (pp[biggest_main_branch_id].length() < other.length() && (other.width.back() != 0 || pp[biggest_main_branch_id].width.back() == 0)))) { + //the main branch should have no endpoint or be bigger! + //here, it have an endpoint, and is not the biggest -> bad! + continue; + } else { + //compute the dot (biggest_main_branch_id) + dot_poly_branch = -dot(Line(polyline.points[0], polyline.points[1]), Line(pp[biggest_main_branch_id].points[0], pp[biggest_main_branch_id].points[1])); + dot_candidate_branch = -dot(Line(other.points[0], other.points[1]), Line(pp[biggest_main_branch_id].points[0], pp[biggest_main_branch_id].points[1])); + if (dot_poly_branch < 0) dot_poly_branch = 0; + if (dot_candidate_branch < 0) dot_candidate_branch = 0; + if (pp[biggest_main_branch_id].width.back()>0) + test_dot += 2 * dot_poly_branch ; + } + //test if it's useful to merge or not + //ie, don't merge 'T' but ok for 'Y', merge only lines of not disproportionate different length (ratio max: 4) (or they are both with 0-width end) + if (dot_poly_branch < 0.1 || dot_candidate_branch < 0.1 || + ( + ((polyline.length()>other.length() ? polyline.length() / other.length() : other.length() / polyline.length()) > 4) + && !(polyline.width.back() == 0 && other.width.back()==0) + ) ){ + continue; + } + if (test_dot > best_dot) { + best_candidate = &other; + best_idx = j; + best_dot = test_dot; + } + } + if (best_candidate != nullptr) { + //idf++; + //std::cout << " == fusion " << id <<" : "<< idf << " ==\n"; + // delete very near points + remove_point_too_near(&polyline); + remove_point_too_near(best_candidate); + + // add point at the same pos than the other line to have a nicer fusion + add_point_same_percent(&polyline, best_candidate); + add_point_same_percent(best_candidate, &polyline); + + //get the angle of the nearest points of the contour to see : _| (good) \_ (average) __(bad) + //sqrt because the result are nicer this way: don't over-penalize /_ angles + //TODO: try if we can achieve a better result if we use a different algo if the angle is <90° + const double coeff_angle_poly = (coeff_angle_cache.find(polyline.points.back()) != coeff_angle_cache.end()) + ? coeff_angle_cache[polyline.points.back()] + : (get_coeff_from_angle_countour(polyline.points.back(), this->expolygon, min(min_width, polyline.length() / 2))); + const double coeff_angle_candi = (coeff_angle_cache.find(best_candidate->points.back()) != coeff_angle_cache.end()) + ? coeff_angle_cache[best_candidate->points.back()] + : (get_coeff_from_angle_countour(best_candidate->points.back(), this->expolygon, min(min_width, best_candidate->length() / 2))); + + //this will encourage to follow the curve, a little, because it's shorter near the center + //without that, it tends to go to the outter rim. + //std::cout << " max(polyline.length(), best_candidate->length())=" << max(polyline.length(), best_candidate->length()) + // << ", polyline.length()=" << polyline.length() + // << ", best_candidate->length()=" << best_candidate->length() + // << ", polyline.length() / max=" << (polyline.length() / max(polyline.length(), best_candidate->length())) + // << ", best_candidate->length() / max=" << (best_candidate->length() / max(polyline.length(), best_candidate->length())) + // << "\n"; + double weight_poly = 2 - (polyline.length() / max(polyline.length(), best_candidate->length())); + double weight_candi = 2 - (best_candidate->length() / max(polyline.length(), best_candidate->length())); + weight_poly *= coeff_angle_poly; + weight_candi *= coeff_angle_candi; + const double coeff_poly = (dot_poly_branch * weight_poly) / (dot_poly_branch * weight_poly + dot_candidate_branch * weight_candi); + const double coeff_candi = 1.0 - coeff_poly; + //std::cout << "coeff_angle_poly=" << coeff_angle_poly + // << ", coeff_angle_candi=" << coeff_angle_candi + // << ", weight_poly=" << (2 - (polyline.length() / max(polyline.length(), best_candidate->length()))) + // << ", weight_candi=" << (2 - (best_candidate->length() / max(polyline.length(), best_candidate->length()))) + // << ", sumpoly=" << weight_poly + // << ", sumcandi=" << weight_candi + // << ", dot_poly_branch=" << dot_poly_branch + // << ", dot_candidate_branch=" << dot_candidate_branch + // << ", coeff_poly=" << coeff_poly + // << ", coeff_candi=" << coeff_candi + // << "\n"; + //iterate the points + // as voronoi should create symetric thing, we can iterate synchonously + size_t idx_point = 1; + while (idx_point < min(polyline.points.size(), best_candidate->points.size())) { + //fusion + polyline.points[idx_point].x = polyline.points[idx_point].x * coeff_poly + best_candidate->points[idx_point].x * coeff_candi; + polyline.points[idx_point].y = polyline.points[idx_point].y * coeff_poly + best_candidate->points[idx_point].y * coeff_candi; + + // The width decrease with distance from the centerline. + // This formula is what works the best, even if it's not perfect (created empirically). 0->3% error on a gap fill on some tests. + //If someone find an other formula based on the properties of the voronoi algorithm used here, and it works better, please use it. + //or maybe just use the distance to nearest edge in bounds... + double value_from_current_width = 0.5*polyline.width[idx_point] * dot_poly_branch / max(dot_poly_branch, dot_candidate_branch); + value_from_current_width += 0.5*best_candidate->width[idx_point] * dot_candidate_branch / max(dot_poly_branch, dot_candidate_branch); + double value_from_dist = 2 * polyline.points[idx_point].distance_to(best_candidate->points[idx_point]); + value_from_dist *= sqrt(min(dot_poly_branch, dot_candidate_branch) / max(dot_poly_branch, dot_candidate_branch)); + polyline.width[idx_point] = value_from_current_width + value_from_dist; + //std::cout << "width:" << polyline.width[idx_point] << " = " << value_from_current_width << " + " << value_from_dist + // << " (<" << max_width << " && " << (bounds.contour.closest_point(polyline.points[idx_point])->distance_to(polyline.points[idx_point]) * 2.1)<<")\n"; + //failsafes + if (polyline.width[idx_point] > max_width) + polyline.width[idx_point] = max_width; + const coord_t max_width_contour = bounds.contour.closest_point(polyline.points[idx_point])->distance_to(polyline.points[idx_point]) * 2.1; + if (polyline.width[idx_point] > max_width_contour) + polyline.width[idx_point] = max_width_contour; + + ++idx_point; + } + if (idx_point < best_candidate->points.size()) { + if (idx_point + 1 < best_candidate->points.size()) { + //create a new polyline + pp.emplace_back(); + pp.back().endpoints.first = true; + pp.back().endpoints.second = best_candidate->endpoints.second; + for (size_t idx_point_new_line = idx_point; idx_point_new_line < best_candidate->points.size(); ++idx_point_new_line) { + pp.back().points.push_back(best_candidate->points[idx_point_new_line]); + pp.back().width.push_back(best_candidate->width[idx_point_new_line]); + } + } else { + //Add last point + polyline.points.push_back(best_candidate->points[idx_point]); + polyline.width.push_back(best_candidate->width[idx_point]); + //select if an end opccur + polyline.endpoints.second &= best_candidate->endpoints.second; + } + + } else { + //select if an end opccur + polyline.endpoints.second &= best_candidate->endpoints.second; + } + + //remove points that are the same or too close each other, ie simplify + for (size_t idx_point = 1; idx_point < polyline.points.size(); ++idx_point) { + if (polyline.points[idx_point - 1].distance_to(polyline.points[idx_point]) < SCALED_EPSILON) { + if (idx_point < polyline.points.size() - 1) { + polyline.points.erase(polyline.points.begin() + idx_point); + polyline.width.erase(polyline.width.begin() + idx_point); + } else { + polyline.points.erase(polyline.points.begin() + idx_point - 1); + polyline.width.erase(polyline.width.begin() + idx_point - 1); + } + --idx_point; + } + } + //remove points that are outside of the geometry + for (size_t idx_point = 0; idx_point < polyline.points.size(); ++idx_point) { + if (!bounds.contains_b(polyline.points[idx_point])) { + polyline.points.erase(polyline.points.begin() + idx_point); + polyline.width.erase(polyline.width.begin() + idx_point); + --idx_point; + } + } + + //update cache + coeff_angle_cache[polyline.points.back()] = coeff_angle_poly * coeff_poly + coeff_angle_candi * coeff_candi; + + + if (polyline.points.size() < 2) { + //remove self + pp.erase(pp.begin() + i); + --i; + --best_idx; + } + + pp.erase(pp.begin() + best_idx); + //{ + // stringstream stri; + // stri << "medial_axis_2.0_aft_fus_" << id << "_" << idf << ".svg"; + // SVG svg(stri.str()); + // svg.draw(bounds); + // svg.draw(this->expolygon); + // svg.draw(pp); + // svg.Close(); + //} + changes = true; + break; + } + } + if (changes) { + concatThickPolylines(pp); + ///reorder, in case of change + std::sort(pp.begin(), pp.end(), [](const ThickPolyline & a, const ThickPolyline & b) { + bool ahas0 = a.width.front() == 0 || a.width.back() == 0; + bool bhas0 = b.width.front() == 0 || b.width.back() == 0; + if (ahas0 && !bhas0) return true; + if (!ahas0 && bhas0) return false; + return a.length() < b.length(); + }); + } + } +} + +void +MedialAxis::remove_too_thin_extrusion(ThickPolylines& pp) +{ + // remove too thin extrusion at start & end of polylines + bool changes = false; + for (size_t i = 0; i < pp.size(); ++i) { + ThickPolyline& polyline = pp[i]; + // remove bits with too small extrusion + while (polyline.points.size() > 1 && polyline.width.front() < this->min_width && polyline.endpoints.first) { + //try to split if possible + if (polyline.width[1] > min_width) { + double percent_can_keep = (min_width - polyline.width[0]) / (polyline.width[1] - polyline.width[0]); + if (polyline.points.front().distance_to(polyline.points[1]) * percent_can_keep > this->max_width / 2 + && polyline.points.front().distance_to(polyline.points[1])* (1 - percent_can_keep) > this->max_width / 2) { + //Can split => move the first point and assign a new weight. + //the update of endpoints wil be performed in concatThickPolylines + polyline.points.front().x = polyline.points.front().x + + (coord_t)((polyline.points[1].x - polyline.points.front().x) * percent_can_keep); + polyline.points.front().y = polyline.points.front().y + + (coord_t)((polyline.points[1].y - polyline.points.front().y) * percent_can_keep); + polyline.width.front() = min_width; + changes = true; + break; + } + } + polyline.points.erase(polyline.points.begin()); + polyline.width.erase(polyline.width.begin()); + changes = true; + } + while (polyline.points.size() > 1 && polyline.width.back() < this->min_width && polyline.endpoints.second) { + //try to split if possible + if (polyline.width[polyline.points.size() - 2] > min_width) { + double percent_can_keep = (min_width - polyline.width.back()) / (polyline.width[polyline.points.size() - 2] - polyline.width.back()); + if (polyline.points.back().distance_to(polyline.points[polyline.points.size() - 2]) * percent_can_keep > this->max_width / 2 + && polyline.points.back().distance_to(polyline.points[polyline.points.size() - 2]) * (1 - percent_can_keep) > this->max_width / 2) { + //Can split => move the first point and assign a new weight. + //the update of endpoints wil be performed in concatThickPolylines + polyline.points.back().x = polyline.points.back().x + + (coord_t)((polyline.points[polyline.points.size() - 2].x - polyline.points.back().x) * percent_can_keep); + polyline.points.back().y = polyline.points.back().y + + (coord_t)((polyline.points[polyline.points.size() - 2].y - polyline.points.back().y) * percent_can_keep); + polyline.width.back() = min_width; + changes = true; + break; + } + } + polyline.points.erase(polyline.points.end() - 1); + polyline.width.erase(polyline.width.end() - 1); + changes = true; + } + if (polyline.points.size() < 2) { + //remove self if too small + pp.erase(pp.begin() + i); + --i; + } + } + if (changes) concatThickPolylines(pp); +} + +void +MedialAxis::concatenate_polylines_with_crossing(ThickPolylines& pp) +{ + + // concatenate, but even where multiple thickpolyline join, to create nice long strait polylines + /* If we removed any short polylines we now try to connect consecutive polylines + in order to allow loop detection. Note that this algorithm is greedier than + MedialAxis::process_edge_neighbors() as it will connect random pairs of + polylines even when more than two start from the same point. This has no + drawbacks since we optimize later using nearest-neighbor which would do the + same, but should we use a more sophisticated optimization algorithm we should + not connect polylines when more than two meet. + Optimisation of the old algorithm : now we select the most "strait line" choice + when we merge with an other line at a point with more than two meet. + */ + bool changes = false; + for (size_t i = 0; i < pp.size(); ++i) { + ThickPolyline& polyline = pp[i]; + if (polyline.endpoints.first && polyline.endpoints.second) continue; // optimization + + ThickPolyline* best_candidate = nullptr; + float best_dot = -1; + size_t best_idx = 0; + + // find another polyline starting here + for (size_t j = i + 1; j < pp.size(); ++j) { + ThickPolyline& other = pp[j]; + if (polyline.last_point().coincides_with(other.last_point())) { + other.reverse(); + } else if (polyline.first_point().coincides_with(other.last_point())) { + polyline.reverse(); + other.reverse(); + } else if (polyline.first_point().coincides_with(other.first_point())) { + polyline.reverse(); + } else if (!polyline.last_point().coincides_with(other.first_point())) { + continue; + } + + Pointf v_poly(polyline.lines().back().vector().x, polyline.lines().back().vector().y); + v_poly.scale(1 / std::sqrt(v_poly.x*v_poly.x + v_poly.y*v_poly.y)); + Pointf v_other(other.lines().front().vector().x, other.lines().front().vector().y); + v_other.scale(1 / std::sqrt(v_other.x*v_other.x + v_other.y*v_other.y)); + float other_dot = v_poly.x*v_other.x + v_poly.y*v_other.y; + if (other_dot > best_dot) { + best_candidate = &other; + best_idx = j; + best_dot = other_dot; + } + } + if (best_candidate != nullptr) { + + polyline.points.insert(polyline.points.end(), best_candidate->points.begin() + 1, best_candidate->points.end()); + polyline.width.insert(polyline.width.end(), best_candidate->width.begin() + 1, best_candidate->width.end()); + polyline.endpoints.second = best_candidate->endpoints.second; + assert(polyline.width.size() == polyline.points.size()); + changes = true; + pp.erase(pp.begin() + best_idx); + } + } + if (changes) concatThickPolylines(pp); +} + +void +MedialAxis::remove_too_thin_points(ThickPolylines& pp) +{ + //remove too thin polylines points (inside a polyline : split it) + for (size_t i = 0; i < pp.size(); ++i) { + ThickPolyline& polyline = pp[i]; + + // remove bits with too small extrusion + size_t idx_point = 0; + while (idx_point polyline.length()) { + shortest_size = polyline.length(); + shortest_idx = i; + } + + } + } + if (shortest_idx >= 0 && shortest_idx < pp.size()) { + pp.erase(pp.begin() + shortest_idx); + changes = true; + } + if (changes) concatThickPolylines(pp); + } +} + +void +MedialAxis::ensure_not_overextrude(ThickPolylines& pp) +{ + //ensure the volume extruded is correct for what we have been asked + // => don't over-extrude + double surface = 0; + double volume = 0; + for (ThickPolyline& polyline : pp) { + for (ThickLine l : polyline.thicklines()) { + surface += l.length() * (l.a_width + l.b_width) / 2; + double width_mean = (l.a_width + l.b_width) / 2; + volume += height * (width_mean - height * (1. - 0.25 * PI)) * l.length(); + } + } + + // compute bounds volume + double boundsVolume = 0; + boundsVolume += height*bounds.area(); + // add external "perimeter gap" + double perimeterRoundGap = bounds.contour.length() * height * (1 - 0.25*PI) * 0.5; + // add holes "perimeter gaps" + double holesGaps = 0; + for (auto hole = bounds.holes.begin(); hole != bounds.holes.end(); ++hole) { + holesGaps += hole->length() * height * (1 - 0.25*PI) * 0.5; + } + boundsVolume += perimeterRoundGap + holesGaps; + + if (boundsVolume < volume) { + //reduce width + double reduce_by = boundsVolume / volume; + for (ThickPolyline& polyline : pp) { + for (ThickLine l : polyline.thicklines()) { + l.a_width *= reduce_by; + l.b_width *= reduce_by; + } + } + } +} + +ExPolygon +MedialAxis::simplify_polygon_frontier() +{ + + //simplify the boundary between us and the bounds. + //int firstIdx = 0; + //while (firstIdx < contour.points.size() && bounds.contour.contains(contour.points[firstIdx])) firstIdx++; + ExPolygon simplified_poly = this->surface; + if (&this->surface != &this->bounds) { + bool need_intersect = false; + for (size_t i = 0; i < simplified_poly.contour.points.size(); i++) { + Point &p_check = simplified_poly.contour.points[i]; + //if (!find) { + if (!bounds.has_boundary_point(p_check)) { + //check if we put it at a bound point instead of delete it + size_t prev_i = i == 0 ? simplified_poly.contour.points.size() - 1 : (i - 1); + size_t next_i = i == simplified_poly.contour.points.size() - 1 ? 0 : (i + 1); + const Point* closest = bounds.contour.closest_point(p_check); + if (closest != nullptr && closest->distance_to(p_check) + SCALED_EPSILON + < min(p_check.distance_to(simplified_poly.contour.points[prev_i]), p_check.distance_to(simplified_poly.contour.points[next_i])) / 2) { + p_check.x = closest->x; + p_check.y = closest->y; + need_intersect = true; + } else { + simplified_poly.contour.points.erase(simplified_poly.contour.points.begin() + i); + i--; + } + } + } + if (need_intersect) { + ExPolygons simplified_polys = intersection_ex(simplified_poly, bounds); + if (simplified_polys.size() == 1) { + simplified_poly = simplified_polys[0]; + } else { + simplified_poly = this->surface; + } + } + } + + simplified_poly.remove_point_too_near(SCALED_RESOLUTION); + return simplified_poly; +} + +void +MedialAxis::build(ThickPolylines* polylines_out) +{ + this->id++; + + this->expolygon = simplify_polygon_frontier(); + + + // compute the Voronoi diagram and extract medial axis polylines + ThickPolylines pp; + this->polyline_from_voronoi(this->expolygon.lines(), &pp); + + + //{ + // stringstream stri; + // stri << "medial_axis_1_voronoi_" << id << ".svg"; + // SVG svg(stri.str()); + // svg.draw(bounds); + // svg.draw(this->expolygon); + // svg.draw(pp); + // svg.Close(); + //} + + /* Find the maximum width returned; we're going to use this for validating and + filtering the output segments. */ + double max_w = 0; + for (ThickPolylines::const_iterator it = pp.begin(); it != pp.end(); ++it) + max_w = fmaxf(max_w, *std::max_element(it->width.begin(), it->width.end())); + + + fusion_curve(pp); + //{ + // stringstream stri; + // stri << "medial_axis_2_curve_" << id << ".svg"; + // SVG svg(stri.str()); + // svg.draw(bounds); + // svg.draw(this->expolygon); + // svg.draw(pp); + // svg.Close(); + //} + + concatThickPolylines(pp); + + // Aligned fusion: Fusion the bits at the end of lines by "increasing thickness" + // For that, we have to find other lines, + // and with a next point no more distant than the max width. + // Then, we can merge the bit from the first point to the second by following the mean. + // + main_fusion(pp); + //{ + // stringstream stri; + // stri << "medial_axis_3_fusion_" << id << ".svg"; + // SVG svg(stri.str()); + // svg.draw(bounds); + // svg.draw(this->expolygon); + // svg.draw(pp); + // svg.Close(); + //} + + //fusion right-angle corners. + fusion_corners(pp); + + //reduce extrusion when it's too thin to be printable + remove_too_thin_extrusion(pp); + //{ + // stringstream stri; + // stri << "medial_axis_4_thinok_" << id << ".svg"; + // SVG svg(stri.str()); + // svg.draw(bounds); + // svg.draw(this->expolygon); + // svg.draw(pp); + // svg.Close(); + //} + + remove_too_thin_points(pp); + //{ + // stringstream stri; + // stri << "medial_axis_5.0_thuinner_" << id << ".svg"; + // SVG svg(stri.str()); + // svg.draw(bounds); + // svg.draw(this->expolygon); + // svg.draw(pp); + // svg.Close(); + //} + + // Loop through all returned polylines in order to extend their endpoints to the + // expolygon boundaries + const ExPolygons anchors = offset2_ex(diff_ex(this->bounds, this->expolygon), -SCALED_RESOLUTION, SCALED_RESOLUTION); + for (size_t i = 0; i < pp.size(); ++i) { + ThickPolyline& polyline = pp[i]; + extends_line(polyline, anchors, min_width); + polyline.reverse(); + extends_line(polyline, anchors, min_width); + } + //{ + // stringstream stri; + // stri << "medial_axis_5_expand_" << id << ".svg"; + // SVG svg(stri.str()); + // svg.draw(bounds); + // svg.draw(this->expolygon); + // svg.draw(pp); + // svg.Close(); + //} + + concatenate_polylines_with_crossing(pp); + //{ + // stringstream stri; + // stri << "medial_axis_6_concat_" << id << ".svg"; + // SVG svg(stri.str()); + // svg.draw(bounds); + // svg.draw(this->expolygon); + // svg.draw(pp); + // svg.Close(); + //} + + remove_too_short_polylines(pp, max_w * 2); + //{ + // stringstream stri; + // stri << "medial_axis_8_tooshort_" << id << ".svg"; + // SVG svg(stri.str()); + // svg.draw(bounds); + // svg.draw(this->expolygon); + // svg.draw(pp); + // svg.Close(); + //} + + //TODO: reduce the flow at the intersection ( + ) points ? + ensure_not_overextrude(pp); + //{ + // stringstream stri; + // stri << "medial_axis_9_endn_" << id << ".svg"; + // SVG svg(stri.str()); + // svg.draw(bounds); + // svg.draw(this->expolygon); + // svg.draw(pp); + // svg.Close(); + //} + + polylines_out->insert(polylines_out->end(), pp.begin(), pp.end()); + +} + +} // namespace Slic3r diff --git a/xs/src/libslic3r/MedialAxis.hpp b/xs/src/libslic3r/MedialAxis.hpp new file mode 100644 index 00000000000..c65cbdb84b1 --- /dev/null +++ b/xs/src/libslic3r/MedialAxis.hpp @@ -0,0 +1,63 @@ +#ifndef slic3r_MedialAxis_hpp_ +#define slic3r_MedialAxis_hpp_ + +#include "libslic3r.h" +#include "ExPolygon.hpp" +#include "Polyline.hpp" +#include "Geometry.hpp" +#include + +#include "boost/polygon/voronoi.hpp" +using boost::polygon::voronoi_builder; +using boost::polygon::voronoi_diagram; + +namespace Slic3r { + + class MedialAxis { + public: + Lines lines; //lines is here only to avoid appassing it in argument of amny method. Initialized in polyline_from_voronoi. + ExPolygon expolygon; + const ExPolygon& bounds; + const ExPolygon& surface; + const double max_width; + const double min_width; + const double height; + MedialAxis(const ExPolygon &_expolygon, const ExPolygon &_bounds, const double _max_width, const double _min_width, const double _height) + : surface(_expolygon), bounds(_bounds), max_width(_max_width), min_width(_min_width), height(_height) { + }; + void build(ThickPolylines* polylines_out); + void build(Polylines* polylines); + + private: + static int id; + class VD : public voronoi_diagram { + public: + typedef double coord_type; + typedef boost::polygon::point_data point_type; + typedef boost::polygon::segment_data segment_type; + typedef boost::polygon::rectangle_data rect_type; + }; + VD vd; + std::set edges, valid_edges; + std::map > thickness; + void process_edge_neighbors(const VD::edge_type* edge, ThickPolyline* polyline); + bool validate_edge(const VD::edge_type* edge); + const Line& retrieve_segment(const VD::cell_type* cell) const; + const Point& retrieve_endpoint(const VD::cell_type* cell) const; + void polyline_from_voronoi(const Lines& voronoi_edges, ThickPolylines* polylines_out); + + ExPolygon simplify_polygon_frontier(); + void fusion_curve(ThickPolylines &pp); + void main_fusion(ThickPolylines& pp); + void fusion_corners(ThickPolylines &pp); + void extends_line(ThickPolyline& polyline, const ExPolygons& anchors, const coord_t join_width); + void remove_too_thin_extrusion(ThickPolylines& pp); + void concatenate_polylines_with_crossing(ThickPolylines& pp); + void remove_too_thin_points(ThickPolylines& pp); + void remove_too_short_polylines(ThickPolylines& pp, const coord_t min_size); + void ensure_not_overextrude(ThickPolylines& pp); + }; +} + + +#endif