diff --git a/src/libslic3r/CMakeLists.txt b/src/libslic3r/CMakeLists.txt index 646a64d51..b5e576d72 100644 --- a/src/libslic3r/CMakeLists.txt +++ b/src/libslic3r/CMakeLists.txt @@ -116,6 +116,12 @@ add_library(libslic3r STATIC Geometry.hpp Geometry/Circle.cpp Geometry/Circle.hpp + Geometry/MedialAxis.Cpp + Geometry/MedialAxis.hpp + Geometry/Voronoi.hpp + Geometry/VoronoiOffset.cpp + Geometry/VoronoiOffset.hpp + Geometry/VoronoiVisualUtils.hpp Int128.hpp KDTreeIndirect.hpp Layer.cpp @@ -220,9 +226,6 @@ add_library(libslic3r STATIC TriangleSelector.cpp TriangleSelector.hpp MTUtils.hpp - VoronoiOffset.cpp - VoronoiOffset.hpp - VoronoiVisualUtils.hpp Zipper.hpp Zipper.cpp MinAreaBoundingBox.hpp diff --git a/src/libslic3r/ExPolygon.cpp b/src/libslic3r/ExPolygon.cpp index bfe0479fe..069c6e602 100644 --- a/src/libslic3r/ExPolygon.cpp +++ b/src/libslic3r/ExPolygon.cpp @@ -1,7 +1,7 @@ #include "BoundingBox.hpp" #include "ExPolygon.hpp" #include "Exception.hpp" -#include "Geometry.hpp" +#include "Geometry/MedialAxis.hpp" #include "Polygon.hpp" #include "Line.hpp" #include "ClipperUtils.hpp" diff --git a/src/libslic3r/Geometry.cpp b/src/libslic3r/Geometry.cpp index a999f248b..4c24116e4 100644 --- a/src/libslic3r/Geometry.cpp +++ b/src/libslic3r/Geometry.cpp @@ -26,185 +26,6 @@ #include -#ifdef SLIC3R_DEBUG -#include "SVG.hpp" -#endif - -#ifdef SLIC3R_DEBUG -namespace boost { namespace polygon { - -// The following code for the visualization of the boost Voronoi diagram is based on: -// -// Boost.Polygon library voronoi_graphic_utils.hpp header file -// Copyright Andrii Sydorchuk 2010-2012. -// Distributed under the Boost Software License, Version 1.0. -// (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) -template -class voronoi_visual_utils { - public: - // Discretize parabolic Voronoi edge. - // Parabolic Voronoi edges are always formed by one point and one segment - // from the initial input set. - // - // Args: - // point: input point. - // segment: input segment. - // max_dist: maximum discretization distance. - // discretization: point discretization of the given Voronoi edge. - // - // Template arguments: - // InCT: coordinate type of the input geometries (usually integer). - // Point: point type, should model point concept. - // Segment: segment type, should model segment concept. - // - // Important: - // discretization should contain both edge endpoints initially. - template class Point, - template class Segment> - static - typename enable_if< - typename gtl_and< - typename gtl_if< - typename is_point_concept< - typename geometry_concept< Point >::type - >::type - >::type, - typename gtl_if< - typename is_segment_concept< - typename geometry_concept< Segment >::type - >::type - >::type - >::type, - void - >::type discretize( - const Point& point, - const Segment& segment, - const CT max_dist, - std::vector< Point >* discretization) { - // Apply the linear transformation to move start point of the segment to - // the point with coordinates (0, 0) and the direction of the segment to - // coincide the positive direction of the x-axis. - CT segm_vec_x = cast(x(high(segment))) - cast(x(low(segment))); - CT segm_vec_y = cast(y(high(segment))) - cast(y(low(segment))); - CT sqr_segment_length = segm_vec_x * segm_vec_x + segm_vec_y * segm_vec_y; - - // Compute x-coordinates of the endpoints of the edge - // in the transformed space. - CT projection_start = sqr_segment_length * - get_point_projection((*discretization)[0], segment); - CT projection_end = sqr_segment_length * - get_point_projection((*discretization)[1], segment); - - // Compute parabola parameters in the transformed space. - // Parabola has next representation: - // f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y). - CT point_vec_x = cast(x(point)) - cast(x(low(segment))); - CT point_vec_y = cast(y(point)) - cast(y(low(segment))); - CT rot_x = segm_vec_x * point_vec_x + segm_vec_y * point_vec_y; - CT rot_y = segm_vec_x * point_vec_y - segm_vec_y * point_vec_x; - - // Save the last point. - Point last_point = (*discretization)[1]; - discretization->pop_back(); - - // Use stack to avoid recursion. - std::stack point_stack; - point_stack.push(projection_end); - CT cur_x = projection_start; - CT cur_y = parabola_y(cur_x, rot_x, rot_y); - - // Adjust max_dist parameter in the transformed space. - const CT max_dist_transformed = max_dist * max_dist * sqr_segment_length; - while (!point_stack.empty()) { - CT new_x = point_stack.top(); - CT new_y = parabola_y(new_x, rot_x, rot_y); - - // Compute coordinates of the point of the parabola that is - // furthest from the current line segment. - CT mid_x = (new_y - cur_y) / (new_x - cur_x) * rot_y + rot_x; - CT mid_y = parabola_y(mid_x, rot_x, rot_y); - - // Compute maximum distance between the given parabolic arc - // and line segment that discretize it. - CT dist = (new_y - cur_y) * (mid_x - cur_x) - - (new_x - cur_x) * (mid_y - cur_y); - dist = dist * dist / ((new_y - cur_y) * (new_y - cur_y) + - (new_x - cur_x) * (new_x - cur_x)); - if (dist <= max_dist_transformed) { - // Distance between parabola and line segment is less than max_dist. - point_stack.pop(); - CT inter_x = (segm_vec_x * new_x - segm_vec_y * new_y) / - sqr_segment_length + cast(x(low(segment))); - CT inter_y = (segm_vec_x * new_y + segm_vec_y * new_x) / - sqr_segment_length + cast(y(low(segment))); - discretization->push_back(Point(inter_x, inter_y)); - cur_x = new_x; - cur_y = new_y; - } else { - point_stack.push(mid_x); - } - } - - // Update last point. - discretization->back() = last_point; - } - - private: - // Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b). - static CT parabola_y(CT x, CT a, CT b) { - return ((x - a) * (x - a) + b * b) / (b + b); - } - - // Get normalized length of the distance between: - // 1) point projection onto the segment - // 2) start point of the segment - // Return this length divided by the segment length. This is made to avoid - // sqrt computation during transformation from the initial space to the - // transformed one and vice versa. The assumption is made that projection of - // the point lies between the start-point and endpoint of the segment. - template class Point, - template class Segment> - static - typename enable_if< - typename gtl_and< - typename gtl_if< - typename is_point_concept< - typename geometry_concept< Point >::type - >::type - >::type, - typename gtl_if< - typename is_segment_concept< - typename geometry_concept< Segment >::type - >::type - >::type - >::type, - CT - >::type get_point_projection( - const Point& point, const Segment& segment) { - CT segment_vec_x = cast(x(high(segment))) - cast(x(low(segment))); - CT segment_vec_y = cast(y(high(segment))) - cast(y(low(segment))); - CT point_vec_x = x(point) - cast(x(low(segment))); - CT point_vec_y = y(point) - cast(y(low(segment))); - CT sqr_segment_length = - segment_vec_x * segment_vec_x + segment_vec_y * segment_vec_y; - CT vec_dot = segment_vec_x * point_vec_x + segment_vec_y * point_vec_y; - return vec_dot / sqr_segment_length; - } - - template - static CT cast(const InCT& value) { - return static_cast(value); - } -}; - -} } // namespace boost::polygon -#endif - -using namespace boost::polygon; // provides also high() and low() - namespace Slic3r { namespace Geometry { // This implementation is based on Andrew's monotone chain 2D convex hull algorithm @@ -562,212 +383,6 @@ arrange(size_t total_parts, const Vec2d &part_size, coordf_t dist, const Boundin } #endif -#ifdef SLIC3R_DEBUG -// The following code for the visualization of the boost Voronoi diagram is based on: -// -// Boost.Polygon library voronoi_visualizer.cpp file -// Copyright Andrii Sydorchuk 2010-2012. -// Distributed under the Boost Software License, Version 1.0. -// (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) -namespace Voronoi { namespace Internal { - - typedef double coordinate_type; - typedef boost::polygon::point_data point_type; - typedef boost::polygon::segment_data segment_type; - typedef boost::polygon::rectangle_data rect_type; - typedef boost::polygon::voronoi_diagram VD; - typedef VD::cell_type cell_type; - typedef VD::cell_type::source_index_type source_index_type; - typedef VD::cell_type::source_category_type source_category_type; - typedef VD::edge_type edge_type; - typedef VD::cell_container_type cell_container_type; - typedef VD::cell_container_type vertex_container_type; - typedef VD::edge_container_type edge_container_type; - typedef VD::const_cell_iterator const_cell_iterator; - typedef VD::const_vertex_iterator const_vertex_iterator; - typedef VD::const_edge_iterator const_edge_iterator; - - static const std::size_t EXTERNAL_COLOR = 1; - - inline void color_exterior(const VD::edge_type* edge) - { - if (edge->color() == EXTERNAL_COLOR) - return; - edge->color(EXTERNAL_COLOR); - edge->twin()->color(EXTERNAL_COLOR); - const VD::vertex_type* v = edge->vertex1(); - if (v == NULL || !edge->is_primary()) - return; - v->color(EXTERNAL_COLOR); - const VD::edge_type* e = v->incident_edge(); - do { - color_exterior(e); - e = e->rot_next(); - } while (e != v->incident_edge()); - } - - inline point_type retrieve_point(const std::vector &segments, const cell_type& cell) - { - assert(cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT || cell.source_category() == SOURCE_CATEGORY_SEGMENT_END_POINT); - return (cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT) ? low(segments[cell.source_index()]) : high(segments[cell.source_index()]); - } - - inline void clip_infinite_edge(const std::vector &segments, const edge_type& edge, coordinate_type bbox_max_size, std::vector* clipped_edge) - { - const cell_type& cell1 = *edge.cell(); - const cell_type& cell2 = *edge.twin()->cell(); - point_type origin, direction; - // Infinite edges could not be created by two segment sites. - if (cell1.contains_point() && cell2.contains_point()) { - point_type p1 = retrieve_point(segments, cell1); - point_type p2 = retrieve_point(segments, cell2); - origin.x((p1.x() + p2.x()) * 0.5); - origin.y((p1.y() + p2.y()) * 0.5); - direction.x(p1.y() - p2.y()); - direction.y(p2.x() - p1.x()); - } else { - origin = cell1.contains_segment() ? retrieve_point(segments, cell2) : retrieve_point(segments, cell1); - segment_type segment = cell1.contains_segment() ? segments[cell1.source_index()] : segments[cell2.source_index()]; - coordinate_type dx = high(segment).x() - low(segment).x(); - coordinate_type dy = high(segment).y() - low(segment).y(); - if ((low(segment) == origin) ^ cell1.contains_point()) { - direction.x(dy); - direction.y(-dx); - } else { - direction.x(-dy); - direction.y(dx); - } - } - coordinate_type koef = bbox_max_size / (std::max)(fabs(direction.x()), fabs(direction.y())); - if (edge.vertex0() == NULL) { - clipped_edge->push_back(point_type( - origin.x() - direction.x() * koef, - origin.y() - direction.y() * koef)); - } else { - clipped_edge->push_back( - point_type(edge.vertex0()->x(), edge.vertex0()->y())); - } - if (edge.vertex1() == NULL) { - clipped_edge->push_back(point_type( - origin.x() + direction.x() * koef, - origin.y() + direction.y() * koef)); - } else { - clipped_edge->push_back( - point_type(edge.vertex1()->x(), edge.vertex1()->y())); - } - } - - inline void sample_curved_edge(const std::vector &segments, const edge_type& edge, std::vector &sampled_edge, coordinate_type max_dist) - { - point_type point = edge.cell()->contains_point() ? - retrieve_point(segments, *edge.cell()) : - retrieve_point(segments, *edge.twin()->cell()); - segment_type segment = edge.cell()->contains_point() ? - segments[edge.twin()->cell()->source_index()] : - segments[edge.cell()->source_index()]; - ::boost::polygon::voronoi_visual_utils::discretize(point, segment, max_dist, &sampled_edge); - } - -} /* namespace Internal */ } // namespace Voronoi - -static inline void dump_voronoi_to_svg(const Lines &lines, /* const */ boost::polygon::voronoi_diagram &vd, const ThickPolylines *polylines, const char *path) -{ - const double scale = 0.2; - const std::string inputSegmentPointColor = "lightseagreen"; - const coord_t inputSegmentPointRadius = coord_t(0.09 * scale / SCALING_FACTOR); - const std::string inputSegmentColor = "lightseagreen"; - const coord_t inputSegmentLineWidth = coord_t(0.03 * scale / SCALING_FACTOR); - - const std::string voronoiPointColor = "black"; - const coord_t voronoiPointRadius = coord_t(0.06 * scale / SCALING_FACTOR); - const std::string voronoiLineColorPrimary = "black"; - const std::string voronoiLineColorSecondary = "green"; - const std::string voronoiArcColor = "red"; - const coord_t voronoiLineWidth = coord_t(0.02 * scale / SCALING_FACTOR); - - const bool internalEdgesOnly = false; - const bool primaryEdgesOnly = false; - - BoundingBox bbox = BoundingBox(lines); - bbox.min(0) -= coord_t(1. / SCALING_FACTOR); - bbox.min(1) -= coord_t(1. / SCALING_FACTOR); - bbox.max(0) += coord_t(1. / SCALING_FACTOR); - bbox.max(1) += coord_t(1. / SCALING_FACTOR); - - ::Slic3r::SVG svg(path, bbox); - - if (polylines != NULL) - svg.draw(*polylines, "lime", "lime", voronoiLineWidth); - -// bbox.scale(1.2); - // For clipping of half-lines to some reasonable value. - // The line will then be clipped by the SVG viewer anyway. - const double bbox_dim_max = double(bbox.max(0) - bbox.min(0)) + double(bbox.max(1) - bbox.min(1)); - // For the discretization of the Voronoi parabolic segments. - const double discretization_step = 0.0005 * bbox_dim_max; - - // Make a copy of the input segments with the double type. - std::vector segments; - for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++ it) - segments.push_back(Voronoi::Internal::segment_type( - Voronoi::Internal::point_type(double(it->a(0)), double(it->a(1))), - Voronoi::Internal::point_type(double(it->b(0)), double(it->b(1))))); - - // Color exterior edges. - for (boost::polygon::voronoi_diagram::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) - if (!it->is_finite()) - Voronoi::Internal::color_exterior(&(*it)); - - // Draw the end points of the input polygon. - for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) { - svg.draw(it->a, inputSegmentPointColor, inputSegmentPointRadius); - svg.draw(it->b, inputSegmentPointColor, inputSegmentPointRadius); - } - // Draw the input polygon. - for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) - svg.draw(Line(Point(coord_t(it->a(0)), coord_t(it->a(1))), Point(coord_t(it->b(0)), coord_t(it->b(1)))), inputSegmentColor, inputSegmentLineWidth); - -#if 1 - // Draw voronoi vertices. - for (boost::polygon::voronoi_diagram::const_vertex_iterator it = vd.vertices().begin(); it != vd.vertices().end(); ++it) - if (! internalEdgesOnly || it->color() != Voronoi::Internal::EXTERNAL_COLOR) - svg.draw(Point(coord_t(it->x()), coord_t(it->y())), voronoiPointColor, voronoiPointRadius); - - for (boost::polygon::voronoi_diagram::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) { - if (primaryEdgesOnly && !it->is_primary()) - continue; - if (internalEdgesOnly && (it->color() == Voronoi::Internal::EXTERNAL_COLOR)) - continue; - std::vector samples; - std::string color = voronoiLineColorPrimary; - if (!it->is_finite()) { - Voronoi::Internal::clip_infinite_edge(segments, *it, bbox_dim_max, &samples); - if (! it->is_primary()) - color = voronoiLineColorSecondary; - } else { - // Store both points of the segment into samples. sample_curved_edge will split the initial line - // until the discretization_step is reached. - samples.push_back(Voronoi::Internal::point_type(it->vertex0()->x(), it->vertex0()->y())); - samples.push_back(Voronoi::Internal::point_type(it->vertex1()->x(), it->vertex1()->y())); - if (it->is_curved()) { - Voronoi::Internal::sample_curved_edge(segments, *it, samples, discretization_step); - color = voronoiArcColor; - } else if (! it->is_primary()) - color = voronoiLineColorSecondary; - } - for (std::size_t i = 0; i + 1 < samples.size(); ++i) - svg.draw(Line(Point(coord_t(samples[i].x()), coord_t(samples[i].y())), Point(coord_t(samples[i+1].x()), coord_t(samples[i+1].y()))), color, voronoiLineWidth); - } -#endif - - if (polylines != NULL) - svg.draw(*polylines, "blue", voronoiLineWidth); - - svg.Close(); -} -#endif /* SLIC3R_DEBUG */ - // Euclidian distance of two boost::polygon points. template T dist(const boost::polygon::point_data &p1,const boost::polygon::point_data &p2) @@ -791,331 +406,6 @@ inline point_type project_point_to_segment(segment_type &seg, point_type &px) return point_type(p0(0) + t*dir(0), p0(1) + t*dir(1)); } -template -inline const typename VD::point_type retrieve_cell_point(const typename VD::cell_type& cell, const SEGMENTS &segments) -{ - assert(cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT || cell.source_category() == SOURCE_CATEGORY_SEGMENT_END_POINT); - return (cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT) ? low(segments[cell.source_index()]) : high(segments[cell.source_index()]); -} - -template -inline std::pair -measure_edge_thickness(const VD &vd, const typename VD::edge_type& edge, const SEGMENTS &segments) -{ - typedef typename VD::coord_type T; - const typename VD::point_type pa(edge.vertex0()->x(), edge.vertex0()->y()); - const typename VD::point_type pb(edge.vertex1()->x(), edge.vertex1()->y()); - const typename VD::cell_type &cell1 = *edge.cell(); - const typename VD::cell_type &cell2 = *edge.twin()->cell(); - if (cell1.contains_segment()) { - if (cell2.contains_segment()) { - // Both cells contain a linear segment, the left / right cells are symmetric. - // Project pa, pb to the left segment. - const typename VD::segment_type segment1 = segments[cell1.source_index()]; - const typename VD::point_type p1a = project_point_to_segment(segment1, pa); - const typename VD::point_type p1b = project_point_to_segment(segment1, pb); - return std::pair(T(2.)*dist(pa, p1a), T(2.)*dist(pb, p1b)); - } else { - // 1st cell contains a linear segment, 2nd cell contains a point. - // The medial axis between the cells is a parabolic arc. - // Project pa, pb to the left segment. - const typename VD::point_type p2 = retrieve_cell_point(cell2, segments); - return std::pair(T(2.)*dist(pa, p2), T(2.)*dist(pb, p2)); - } - } else if (cell2.contains_segment()) { - // 1st cell contains a point, 2nd cell contains a linear segment. - // The medial axis between the cells is a parabolic arc. - const typename VD::point_type p1 = retrieve_cell_point(cell1, segments); - return std::pair(T(2.)*dist(pa, p1), T(2.)*dist(pb, p1)); - } else { - // Both cells contain a point. The left / right regions are triangular and symmetric. - const typename VD::point_type p1 = retrieve_cell_point(cell1, segments); - return std::pair(T(2.)*dist(pa, p1), T(2.)*dist(pb, p1)); - } -} - -// Converts the Line instances of Lines vector to VD::segment_type. -template -class Lines2VDSegments -{ -public: - Lines2VDSegments(const Lines &alines) : lines(alines) {} - typename VD::segment_type operator[](size_t idx) const { - return typename VD::segment_type( - typename VD::point_type(typename VD::coord_type(lines[idx].a(0)), typename VD::coord_type(lines[idx].a(1))), - typename VD::point_type(typename VD::coord_type(lines[idx].b(0)), typename VD::coord_type(lines[idx].b(1)))); - } -private: - const Lines &lines; -}; - -void -MedialAxis::build(ThickPolylines* polylines) -{ - construct_voronoi(this->lines.begin(), this->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()*2 - 2); - - // prevent loop endpoints from being extended - if (polyline.first_point() == 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::build(Polylines* polylines) -{ - ThickPolylines tp; - this->build(&tp); - polylines->insert(polylines->end(), tp.begin(), tp.end()); -} - -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].first); - 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 -#ifndef CLIPPERLIB_INT32 - 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; -#endif // CLIPPERLIB_INT32 - - // 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 (this->expolygon != NULL) { - if (line.a == 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); - 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() - ? segment_r.distance_to(line.a)*2 - : (this->retrieve_endpoint(cell_r) - line.a).cast().norm()*2; - - coordf_t w1 = cell_l->contains_segment() - ? segment_l.distance_to(line.b)*2 - : (this->retrieve_endpoint(cell_l) - line.b).cast().norm()*2; - - 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; - } - - if (w0 < this->min_width && w1 < this->min_width) - return false; - - 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 this->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() == SOURCE_CATEGORY_SEGMENT_START_POINT) { - return line.a; - } else { - return line.b; - } -} - void assemble_transform(Transform3d& transform, const Vec3d& translation, const Vec3d& rotation, const Vec3d& scale, const Vec3d& mirror) { transform = Transform3d::Identity(); diff --git a/src/libslic3r/Geometry.hpp b/src/libslic3r/Geometry.hpp index 5ecfcaeff..179f062ca 100644 --- a/src/libslic3r/Geometry.hpp +++ b/src/libslic3r/Geometry.hpp @@ -10,18 +10,6 @@ // Serialization through the Cereal library #include -#define BOOST_VORONOI_USE_GMP 1 - -#ifdef _MSC_VER -// Suppress warning C4146 in OpenVDB: unary minus operator applied to unsigned type, result still unsigned -#pragma warning(push) -#pragma warning(disable : 4146) -#endif // _MSC_VER -#include "boost/polygon/voronoi.hpp" -#ifdef _MSC_VER -#pragma warning(pop) -#endif // _MSC_VER - namespace Slic3r { namespace ClipperLib { @@ -47,7 +35,7 @@ enum Orientation // and d is limited to 63 bits + signum and we are good. static inline Orientation orient(const Point &a, const Point &b, const Point &c) { - // BOOST_STATIC_ASSERT(sizeof(coord_t) * 2 == sizeof(int64_t)); + static_assert(sizeof(coord_t) * 2 == sizeof(int64_t), "orient works with 32 bit coordinates"); int64_t u = int64_t(b(0)) * int64_t(c(1)) - int64_t(b(1)) * int64_t(c(0)); int64_t v = int64_t(a(0)) * int64_t(c(1)) - int64_t(a(1)) * int64_t(c(0)); int64_t w = int64_t(a(0)) * int64_t(b(1)) - int64_t(a(1)) * int64_t(b(0)); @@ -361,36 +349,6 @@ bool arrange( // output Pointfs &positions); -class VoronoiDiagram : public boost::polygon::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; -}; - -class MedialAxis { -public: - Lines lines; - const ExPolygon* expolygon; - double max_width; - double min_width; - MedialAxis(double _max_width, double _min_width, const ExPolygon* _expolygon = NULL) - : expolygon(_expolygon), max_width(_max_width), min_width(_min_width) {}; - void build(ThickPolylines* polylines); - void build(Polylines* polylines); - -private: - using VD = VoronoiDiagram; - 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; -}; - // Sets the given transform by assembling the given transformations in the following order: // 1) mirror // 2) scale diff --git a/src/libslic3r/Geometry/MedialAxis.cpp b/src/libslic3r/Geometry/MedialAxis.cpp new file mode 100644 index 000000000..60af8fa47 --- /dev/null +++ b/src/libslic3r/Geometry/MedialAxis.cpp @@ -0,0 +1,710 @@ +#include "MedialAxis.hpp" + +#ifdef SLIC3R_DEBUG +namespace boost { namespace polygon { + +// The following code for the visualization of the boost Voronoi diagram is based on: +// +// Boost.Polygon library voronoi_graphic_utils.hpp header file +// Copyright Andrii Sydorchuk 2010-2012. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +template +class voronoi_visual_utils { + public: + // Discretize parabolic Voronoi edge. + // Parabolic Voronoi edges are always formed by one point and one segment + // from the initial input set. + // + // Args: + // point: input point. + // segment: input segment. + // max_dist: maximum discretization distance. + // discretization: point discretization of the given Voronoi edge. + // + // Template arguments: + // InCT: coordinate type of the input geometries (usually integer). + // Point: point type, should model point concept. + // Segment: segment type, should model segment concept. + // + // Important: + // discretization should contain both edge endpoints initially. + template class Point, + template class Segment> + static + typename enable_if< + typename gtl_and< + typename gtl_if< + typename is_point_concept< + typename geometry_concept< Point >::type + >::type + >::type, + typename gtl_if< + typename is_segment_concept< + typename geometry_concept< Segment >::type + >::type + >::type + >::type, + void + >::type discretize( + const Point& point, + const Segment& segment, + const CT max_dist, + std::vector< Point >* discretization) { + // Apply the linear transformation to move start point of the segment to + // the point with coordinates (0, 0) and the direction of the segment to + // coincide the positive direction of the x-axis. + CT segm_vec_x = cast(x(high(segment))) - cast(x(low(segment))); + CT segm_vec_y = cast(y(high(segment))) - cast(y(low(segment))); + CT sqr_segment_length = segm_vec_x * segm_vec_x + segm_vec_y * segm_vec_y; + + // Compute x-coordinates of the endpoints of the edge + // in the transformed space. + CT projection_start = sqr_segment_length * + get_point_projection((*discretization)[0], segment); + CT projection_end = sqr_segment_length * + get_point_projection((*discretization)[1], segment); + + // Compute parabola parameters in the transformed space. + // Parabola has next representation: + // f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y). + CT point_vec_x = cast(x(point)) - cast(x(low(segment))); + CT point_vec_y = cast(y(point)) - cast(y(low(segment))); + CT rot_x = segm_vec_x * point_vec_x + segm_vec_y * point_vec_y; + CT rot_y = segm_vec_x * point_vec_y - segm_vec_y * point_vec_x; + + // Save the last point. + Point last_point = (*discretization)[1]; + discretization->pop_back(); + + // Use stack to avoid recursion. + std::stack point_stack; + point_stack.push(projection_end); + CT cur_x = projection_start; + CT cur_y = parabola_y(cur_x, rot_x, rot_y); + + // Adjust max_dist parameter in the transformed space. + const CT max_dist_transformed = max_dist * max_dist * sqr_segment_length; + while (!point_stack.empty()) { + CT new_x = point_stack.top(); + CT new_y = parabola_y(new_x, rot_x, rot_y); + + // Compute coordinates of the point of the parabola that is + // furthest from the current line segment. + CT mid_x = (new_y - cur_y) / (new_x - cur_x) * rot_y + rot_x; + CT mid_y = parabola_y(mid_x, rot_x, rot_y); + + // Compute maximum distance between the given parabolic arc + // and line segment that discretize it. + CT dist = (new_y - cur_y) * (mid_x - cur_x) - + (new_x - cur_x) * (mid_y - cur_y); + dist = dist * dist / ((new_y - cur_y) * (new_y - cur_y) + + (new_x - cur_x) * (new_x - cur_x)); + if (dist <= max_dist_transformed) { + // Distance between parabola and line segment is less than max_dist. + point_stack.pop(); + CT inter_x = (segm_vec_x * new_x - segm_vec_y * new_y) / + sqr_segment_length + cast(x(low(segment))); + CT inter_y = (segm_vec_x * new_y + segm_vec_y * new_x) / + sqr_segment_length + cast(y(low(segment))); + discretization->push_back(Point(inter_x, inter_y)); + cur_x = new_x; + cur_y = new_y; + } else { + point_stack.push(mid_x); + } + } + + // Update last point. + discretization->back() = last_point; + } + + private: + // Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b). + static CT parabola_y(CT x, CT a, CT b) { + return ((x - a) * (x - a) + b * b) / (b + b); + } + + // Get normalized length of the distance between: + // 1) point projection onto the segment + // 2) start point of the segment + // Return this length divided by the segment length. This is made to avoid + // sqrt computation during transformation from the initial space to the + // transformed one and vice versa. The assumption is made that projection of + // the point lies between the start-point and endpoint of the segment. + template class Point, + template class Segment> + static + typename enable_if< + typename gtl_and< + typename gtl_if< + typename is_point_concept< + typename geometry_concept< Point >::type + >::type + >::type, + typename gtl_if< + typename is_segment_concept< + typename geometry_concept< Segment >::type + >::type + >::type + >::type, + CT + >::type get_point_projection( + const Point& point, const Segment& segment) { + CT segment_vec_x = cast(x(high(segment))) - cast(x(low(segment))); + CT segment_vec_y = cast(y(high(segment))) - cast(y(low(segment))); + CT point_vec_x = x(point) - cast(x(low(segment))); + CT point_vec_y = y(point) - cast(y(low(segment))); + CT sqr_segment_length = + segment_vec_x * segment_vec_x + segment_vec_y * segment_vec_y; + CT vec_dot = segment_vec_x * point_vec_x + segment_vec_y * point_vec_y; + return vec_dot / sqr_segment_length; + } + + template + static CT cast(const InCT& value) { + return static_cast(value); + } +}; + +} } // namespace boost::polygon +#endif // SLIC3R_DEBUG + +namespace Slic3r { namespace Geometry { + + +#ifdef SLIC3R_DEBUG +// The following code for the visualization of the boost Voronoi diagram is based on: +// +// Boost.Polygon library voronoi_visualizer.cpp file +// Copyright Andrii Sydorchuk 2010-2012. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +namespace Voronoi { namespace Internal { + + typedef double coordinate_type; + typedef boost::polygon::point_data point_type; + typedef boost::polygon::segment_data segment_type; + typedef boost::polygon::rectangle_data rect_type; + typedef boost::polygon::voronoi_diagram VD; + typedef VD::cell_type cell_type; + typedef VD::cell_type::source_index_type source_index_type; + typedef VD::cell_type::source_category_type source_category_type; + typedef VD::edge_type edge_type; + typedef VD::cell_container_type cell_container_type; + typedef VD::cell_container_type vertex_container_type; + typedef VD::edge_container_type edge_container_type; + typedef VD::const_cell_iterator const_cell_iterator; + typedef VD::const_vertex_iterator const_vertex_iterator; + typedef VD::const_edge_iterator const_edge_iterator; + + static const std::size_t EXTERNAL_COLOR = 1; + + inline void color_exterior(const VD::edge_type* edge) + { + if (edge->color() == EXTERNAL_COLOR) + return; + edge->color(EXTERNAL_COLOR); + edge->twin()->color(EXTERNAL_COLOR); + const VD::vertex_type* v = edge->vertex1(); + if (v == NULL || !edge->is_primary()) + return; + v->color(EXTERNAL_COLOR); + const VD::edge_type* e = v->incident_edge(); + do { + color_exterior(e); + e = e->rot_next(); + } while (e != v->incident_edge()); + } + + inline point_type retrieve_point(const std::vector &segments, const cell_type& cell) + { + assert(cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT || cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_END_POINT); + return (cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? low(segments[cell.source_index()]) : high(segments[cell.source_index()]); + } + + inline void clip_infinite_edge(const std::vector &segments, const edge_type& edge, coordinate_type bbox_max_size, std::vector* clipped_edge) + { + const cell_type& cell1 = *edge.cell(); + const cell_type& cell2 = *edge.twin()->cell(); + point_type origin, direction; + // Infinite edges could not be created by two segment sites. + if (cell1.contains_point() && cell2.contains_point()) { + point_type p1 = retrieve_point(segments, cell1); + point_type p2 = retrieve_point(segments, cell2); + origin.x((p1.x() + p2.x()) * 0.5); + origin.y((p1.y() + p2.y()) * 0.5); + direction.x(p1.y() - p2.y()); + direction.y(p2.x() - p1.x()); + } else { + origin = cell1.contains_segment() ? retrieve_point(segments, cell2) : retrieve_point(segments, cell1); + segment_type segment = cell1.contains_segment() ? segments[cell1.source_index()] : segments[cell2.source_index()]; + coordinate_type dx = high(segment).x() - low(segment).x(); + coordinate_type dy = high(segment).y() - low(segment).y(); + if ((low(segment) == origin) ^ cell1.contains_point()) { + direction.x(dy); + direction.y(-dx); + } else { + direction.x(-dy); + direction.y(dx); + } + } + coordinate_type koef = bbox_max_size / (std::max)(fabs(direction.x()), fabs(direction.y())); + if (edge.vertex0() == NULL) { + clipped_edge->push_back(point_type( + origin.x() - direction.x() * koef, + origin.y() - direction.y() * koef)); + } else { + clipped_edge->push_back( + point_type(edge.vertex0()->x(), edge.vertex0()->y())); + } + if (edge.vertex1() == NULL) { + clipped_edge->push_back(point_type( + origin.x() + direction.x() * koef, + origin.y() + direction.y() * koef)); + } else { + clipped_edge->push_back( + point_type(edge.vertex1()->x(), edge.vertex1()->y())); + } + } + + inline void sample_curved_edge(const std::vector &segments, const edge_type& edge, std::vector &sampled_edge, coordinate_type max_dist) + { + point_type point = edge.cell()->contains_point() ? + retrieve_point(segments, *edge.cell()) : + retrieve_point(segments, *edge.twin()->cell()); + segment_type segment = edge.cell()->contains_point() ? + segments[edge.twin()->cell()->source_index()] : + segments[edge.cell()->source_index()]; + ::boost::polygon::voronoi_visual_utils::discretize(point, segment, max_dist, &sampled_edge); + } + +} /* namespace Internal */ } // namespace Voronoi + +void dump_voronoi_to_svg(const Lines &lines, /* const */ boost::polygon::voronoi_diagram &vd, const ThickPolylines *polylines, const char *path) +{ + const double scale = 0.2; + const std::string inputSegmentPointColor = "lightseagreen"; + const coord_t inputSegmentPointRadius = coord_t(0.09 * scale / SCALING_FACTOR); + const std::string inputSegmentColor = "lightseagreen"; + const coord_t inputSegmentLineWidth = coord_t(0.03 * scale / SCALING_FACTOR); + + const std::string voronoiPointColor = "black"; + const coord_t voronoiPointRadius = coord_t(0.06 * scale / SCALING_FACTOR); + const std::string voronoiLineColorPrimary = "black"; + const std::string voronoiLineColorSecondary = "green"; + const std::string voronoiArcColor = "red"; + const coord_t voronoiLineWidth = coord_t(0.02 * scale / SCALING_FACTOR); + + const bool internalEdgesOnly = false; + const bool primaryEdgesOnly = false; + + BoundingBox bbox = BoundingBox(lines); + bbox.min(0) -= coord_t(1. / SCALING_FACTOR); + bbox.min(1) -= coord_t(1. / SCALING_FACTOR); + bbox.max(0) += coord_t(1. / SCALING_FACTOR); + bbox.max(1) += coord_t(1. / SCALING_FACTOR); + + ::Slic3r::SVG svg(path, bbox); + + if (polylines != NULL) + svg.draw(*polylines, "lime", "lime", voronoiLineWidth); + +// bbox.scale(1.2); + // For clipping of half-lines to some reasonable value. + // The line will then be clipped by the SVG viewer anyway. + const double bbox_dim_max = double(bbox.max(0) - bbox.min(0)) + double(bbox.max(1) - bbox.min(1)); + // For the discretization of the Voronoi parabolic segments. + const double discretization_step = 0.0005 * bbox_dim_max; + + // Make a copy of the input segments with the double type. + std::vector segments; + for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++ it) + segments.push_back(Voronoi::Internal::segment_type( + Voronoi::Internal::point_type(double(it->a(0)), double(it->a(1))), + Voronoi::Internal::point_type(double(it->b(0)), double(it->b(1))))); + + // Color exterior edges. + for (boost::polygon::voronoi_diagram::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) + if (!it->is_finite()) + Voronoi::Internal::color_exterior(&(*it)); + + // Draw the end points of the input polygon. + for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) { + svg.draw(it->a, inputSegmentPointColor, inputSegmentPointRadius); + svg.draw(it->b, inputSegmentPointColor, inputSegmentPointRadius); + } + // Draw the input polygon. + for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) + svg.draw(Line(Point(coord_t(it->a(0)), coord_t(it->a(1))), Point(coord_t(it->b(0)), coord_t(it->b(1)))), inputSegmentColor, inputSegmentLineWidth); + +#if 1 + // Draw voronoi vertices. + for (boost::polygon::voronoi_diagram::const_vertex_iterator it = vd.vertices().begin(); it != vd.vertices().end(); ++it) + if (! internalEdgesOnly || it->color() != Voronoi::Internal::EXTERNAL_COLOR) + svg.draw(Point(coord_t(it->x()), coord_t(it->y())), voronoiPointColor, voronoiPointRadius); + + for (boost::polygon::voronoi_diagram::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) { + if (primaryEdgesOnly && !it->is_primary()) + continue; + if (internalEdgesOnly && (it->color() == Voronoi::Internal::EXTERNAL_COLOR)) + continue; + std::vector samples; + std::string color = voronoiLineColorPrimary; + if (!it->is_finite()) { + Voronoi::Internal::clip_infinite_edge(segments, *it, bbox_dim_max, &samples); + if (! it->is_primary()) + color = voronoiLineColorSecondary; + } else { + // Store both points of the segment into samples. sample_curved_edge will split the initial line + // until the discretization_step is reached. + samples.push_back(Voronoi::Internal::point_type(it->vertex0()->x(), it->vertex0()->y())); + samples.push_back(Voronoi::Internal::point_type(it->vertex1()->x(), it->vertex1()->y())); + if (it->is_curved()) { + Voronoi::Internal::sample_curved_edge(segments, *it, samples, discretization_step); + color = voronoiArcColor; + } else if (! it->is_primary()) + color = voronoiLineColorSecondary; + } + for (std::size_t i = 0; i + 1 < samples.size(); ++i) + svg.draw(Line(Point(coord_t(samples[i].x()), coord_t(samples[i].y())), Point(coord_t(samples[i+1].x()), coord_t(samples[i+1].y()))), color, voronoiLineWidth); + } +#endif + + if (polylines != NULL) + svg.draw(*polylines, "blue", voronoiLineWidth); + + svg.Close(); +} +#endif // SLIC3R_DEBUG + +template +inline const typename VD::point_type retrieve_cell_point(const typename VD::cell_type& cell, const SEGMENTS &segments) +{ + assert(cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT || cell.source_category() == SOURCE_CATEGORY_SEGMENT_END_POINT); + return (cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT) ? low(segments[cell.source_index()]) : high(segments[cell.source_index()]); +} + +template +inline std::pair +measure_edge_thickness(const VD &vd, const typename VD::edge_type& edge, const SEGMENTS &segments) +{ + typedef typename VD::coord_type T; + const typename VD::point_type pa(edge.vertex0()->x(), edge.vertex0()->y()); + const typename VD::point_type pb(edge.vertex1()->x(), edge.vertex1()->y()); + const typename VD::cell_type &cell1 = *edge.cell(); + const typename VD::cell_type &cell2 = *edge.twin()->cell(); + if (cell1.contains_segment()) { + if (cell2.contains_segment()) { + // Both cells contain a linear segment, the left / right cells are symmetric. + // Project pa, pb to the left segment. + const typename VD::segment_type segment1 = segments[cell1.source_index()]; + const typename VD::point_type p1a = project_point_to_segment(segment1, pa); + const typename VD::point_type p1b = project_point_to_segment(segment1, pb); + return std::pair(T(2.)*dist(pa, p1a), T(2.)*dist(pb, p1b)); + } else { + // 1st cell contains a linear segment, 2nd cell contains a point. + // The medial axis between the cells is a parabolic arc. + // Project pa, pb to the left segment. + const typename VD::point_type p2 = retrieve_cell_point(cell2, segments); + return std::pair(T(2.)*dist(pa, p2), T(2.)*dist(pb, p2)); + } + } else if (cell2.contains_segment()) { + // 1st cell contains a point, 2nd cell contains a linear segment. + // The medial axis between the cells is a parabolic arc. + const typename VD::point_type p1 = retrieve_cell_point(cell1, segments); + return std::pair(T(2.)*dist(pa, p1), T(2.)*dist(pb, p1)); + } else { + // Both cells contain a point. The left / right regions are triangular and symmetric. + const typename VD::point_type p1 = retrieve_cell_point(cell1, segments); + return std::pair(T(2.)*dist(pa, p1), T(2.)*dist(pb, p1)); + } +} + +// Converts the Line instances of Lines vector to VD::segment_type. +template +class Lines2VDSegments +{ +public: + Lines2VDSegments(const Lines &alines) : lines(alines) {} + typename VD::segment_type operator[](size_t idx) const { + return typename VD::segment_type( + typename VD::point_type(typename VD::coord_type(lines[idx].a(0)), typename VD::coord_type(lines[idx].a(1))), + typename VD::point_type(typename VD::coord_type(lines[idx].b(0)), typename VD::coord_type(lines[idx].b(1)))); + } +private: + const Lines &lines; +}; + +void +MedialAxis::build(ThickPolylines* polylines) +{ + construct_voronoi(this->lines.begin(), this->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()*2 - 2); + + // prevent loop endpoints from being extended + if (polyline.first_point() == 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::build(Polylines* polylines) +{ + ThickPolylines tp; + this->build(&tp); + polylines->insert(polylines->end(), tp.begin(), tp.end()); +} + +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].first); + 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 +#ifndef CLIPPERLIB_INT32 + 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; +#endif // CLIPPERLIB_INT32 + + // 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 (this->expolygon != NULL) { + if (line.a == 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); + 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() + ? segment_r.distance_to(line.a)*2 + : (this->retrieve_endpoint(cell_r) - line.a).cast().norm()*2; + + coordf_t w1 = cell_l->contains_segment() + ? segment_l.distance_to(line.b)*2 + : (this->retrieve_endpoint(cell_l) - line.b).cast().norm()*2; + + 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; + } + + if (w0 < this->min_width && w1 < this->min_width) + return false; + + 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 this->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; + } +} + +} } // namespace Slicer::Geometry diff --git a/src/libslic3r/Geometry/MedialAxis.hpp b/src/libslic3r/Geometry/MedialAxis.hpp new file mode 100644 index 000000000..4189b5d4e --- /dev/null +++ b/src/libslic3r/Geometry/MedialAxis.hpp @@ -0,0 +1,32 @@ +#ifndef slic3r_Geometry_MedialAxis_hpp_ +#define slic3r_Geometry_MedialAxis_hpp_ + +#include "Voronoi.hpp" + +namespace Slic3r { namespace Geometry { + +class MedialAxis { +public: + Lines lines; + const ExPolygon* expolygon; + double max_width; + double min_width; + MedialAxis(double _max_width, double _min_width, const ExPolygon* _expolygon = NULL) + : expolygon(_expolygon), max_width(_max_width), min_width(_min_width) {}; + void build(ThickPolylines* polylines); + void build(Polylines* polylines); + +private: + using VD = VoronoiDiagram; + 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; +}; + +} } // namespace Slicer::Geometry + +#endif // slic3r_Geometry_MedialAxis_hpp_ diff --git a/src/libslic3r/Geometry/Voronoi.hpp b/src/libslic3r/Geometry/Voronoi.hpp new file mode 100644 index 000000000..5529750f3 --- /dev/null +++ b/src/libslic3r/Geometry/Voronoi.hpp @@ -0,0 +1,33 @@ +#ifndef slic3r_Geometry_Voronoi_hpp_ +#define slic3r_Geometry_Voronoi_hpp_ + +#include "../Line.hpp" +#include "../Polyline.hpp" + +#define BOOST_VORONOI_USE_GMP 1 + +#ifdef _MSC_VER +// Suppress warning C4146 in OpenVDB: unary minus operator applied to unsigned type, result still unsigned +#pragma warning(push) +#pragma warning(disable : 4146) +#endif // _MSC_VER +#include "boost/polygon/voronoi.hpp" +#ifdef _MSC_VER +#pragma warning(pop) +#endif // _MSC_VER + +namespace Slic3r { + +namespace Geometry { + +class VoronoiDiagram : public boost::polygon::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; +}; + +} } // namespace Slicer::Geometry + +#endif // slic3r_Geometry_Voronoi_hpp_ diff --git a/src/libslic3r/VoronoiOffset.cpp b/src/libslic3r/Geometry/VoronoiOffset.cpp similarity index 99% rename from src/libslic3r/VoronoiOffset.cpp rename to src/libslic3r/Geometry/VoronoiOffset.cpp index e8d13a6ad..46105220a 100644 --- a/src/libslic3r/VoronoiOffset.cpp +++ b/src/libslic3r/Geometry/VoronoiOffset.cpp @@ -1,5 +1,6 @@ // Polygon offsetting using Voronoi diagram prodiced by boost::polygon. +#include "Geometry.hpp" #include "VoronoiOffset.hpp" #include "libslic3r.h" diff --git a/src/libslic3r/VoronoiOffset.hpp b/src/libslic3r/Geometry/VoronoiOffset.hpp similarity index 99% rename from src/libslic3r/VoronoiOffset.hpp rename to src/libslic3r/Geometry/VoronoiOffset.hpp index 17b590ed7..359fe010c 100644 --- a/src/libslic3r/VoronoiOffset.hpp +++ b/src/libslic3r/Geometry/VoronoiOffset.hpp @@ -3,9 +3,9 @@ #ifndef slic3r_VoronoiOffset_hpp_ #define slic3r_VoronoiOffset_hpp_ -#include "libslic3r.h" +#include "../libslic3r.h" -#include "Geometry.hpp" +#include "Voronoi.hpp" namespace Slic3r { diff --git a/src/libslic3r/VoronoiVisualUtils.hpp b/src/libslic3r/Geometry/VoronoiVisualUtils.hpp similarity index 100% rename from src/libslic3r/VoronoiVisualUtils.hpp rename to src/libslic3r/Geometry/VoronoiVisualUtils.hpp diff --git a/src/libslic3r/MultiMaterialSegmentation.cpp b/src/libslic3r/MultiMaterialSegmentation.cpp index 3220a1e02..900073d3a 100644 --- a/src/libslic3r/MultiMaterialSegmentation.cpp +++ b/src/libslic3r/MultiMaterialSegmentation.cpp @@ -3,7 +3,7 @@ #include "EdgeGrid.hpp" #include "Layer.hpp" #include "Print.hpp" -#include "VoronoiVisualUtils.hpp" +#include "Geometry/VoronoiVisualUtils.hpp" #include "MutablePolygon.hpp" #include "format.hpp" diff --git a/tests/libslic3r/test_voronoi.cpp b/tests/libslic3r/test_voronoi.cpp index c78849c01..db12e2fec 100644 --- a/tests/libslic3r/test_voronoi.cpp +++ b/tests/libslic3r/test_voronoi.cpp @@ -6,8 +6,8 @@ #include #include -#include -#include +#include +#include #include