diff options
author | bubnikv <bubnikv@gmail.com> | 2016-09-12 12:29:39 +0300 |
---|---|---|
committer | bubnikv <bubnikv@gmail.com> | 2016-09-12 12:29:39 +0300 |
commit | 9fcc8fe9aeacb8edf6206787574cf1a84c88a5cc (patch) | |
tree | 549776f3a63d4eb240d0caaf7a63fa6259dac031 /xs/src/libslic3r/Geometry.cpp | |
parent | ca98e2655a4a5875947598f404791ea6c0313609 (diff) | |
parent | d022493297a110a9f982786ed696528c1d79d0f1 (diff) |
Meged with release_candidate_1_3
Diffstat (limited to 'xs/src/libslic3r/Geometry.cpp')
-rw-r--r-- | xs/src/libslic3r/Geometry.cpp | 481 |
1 files changed, 481 insertions, 0 deletions
diff --git a/xs/src/libslic3r/Geometry.cpp b/xs/src/libslic3r/Geometry.cpp index 6d864f631..cdf913a27 100644 --- a/xs/src/libslic3r/Geometry.cpp +++ b/xs/src/libslic3r/Geometry.cpp @@ -11,12 +11,186 @@ #include <map> #include <set> #include <utility> +#include <stack> #include <vector> #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 <typename CT> +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 InCT1, class InCT2, + template<class> class Point, + template<class> class Segment> + static + typename enable_if< + typename gtl_and< + typename gtl_if< + typename is_point_concept< + typename geometry_concept< Point<InCT1> >::type + >::type + >::type, + typename gtl_if< + typename is_segment_concept< + typename geometry_concept< Segment<InCT2> >::type + >::type + >::type + >::type, + void + >::type discretize( + const Point<InCT1>& point, + const Segment<InCT2>& segment, + const CT max_dist, + std::vector< Point<CT> >* 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<CT> last_point = (*discretization)[1]; + discretization->pop_back(); + + // Use stack to avoid recursion. + std::stack<CT> 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<CT>(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 InCT, + template<class> class Point, + template<class> class Segment> + static + typename enable_if< + typename gtl_and< + typename gtl_if< + typename is_point_concept< + typename geometry_concept< Point<int> >::type + >::type + >::type, + typename gtl_if< + typename is_segment_concept< + typename geometry_concept< Segment<long> >::type + >::type + >::type + >::type, + CT + >::type get_point_projection( + const Point<CT>& point, const Segment<InCT>& 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 <typename InCT> + static CT cast(const InCT& value) { + return static_cast<CT>(value); + } +}; + +} } // namespace boost::polygon +#endif + using namespace boost::polygon; // provides also high() and low() namespace Slic3r { namespace Geometry { @@ -290,6 +464,294 @@ arrange(size_t total_parts, Pointf part, coordf_t dist, const BoundingBoxf* bb) return positions; } +#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<coordinate_type> point_type; + typedef boost::polygon::segment_data<coordinate_type> segment_type; + typedef boost::polygon::rectangle_data<coordinate_type> rect_type; +// typedef voronoi_builder<int> VB; + typedef boost::polygon::voronoi_diagram<coordinate_type> 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<segment_type> &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<segment_type> &segments, const edge_type& edge, coordinate_type bbox_max_size, std::vector<point_type>* 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<segment_type> &segments, const edge_type& edge, std::vector<point_type> &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<coordinate_type>::discretize(point, segment, max_dist, &sampled_edge); + } + +} /* namespace Internal */ } // namespace Voronoi + +static inline void dump_voronoi_to_svg(const Lines &lines, /* const */ voronoi_diagram<double> &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.x -= coord_t(1. / SCALING_FACTOR); + bbox.min.y -= coord_t(1. / SCALING_FACTOR); + bbox.max.x += coord_t(1. / SCALING_FACTOR); + bbox.max.y += 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.x - bbox.min.x) + double(bbox.max.y - bbox.min.y); + // 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<Voronoi::Internal::segment_type> 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.x), double(it->a.y)), + Voronoi::Internal::point_type(double(it->b.x), double(it->b.y)))); + + // Color exterior edges. + for (voronoi_diagram<double>::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.x), coord_t(it->a.y)), Point(coord_t(it->b.x), coord_t(it->b.y))), inputSegmentColor, inputSegmentLineWidth); + +#if 1 + // Draw voronoi vertices. + for (voronoi_diagram<double>::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 (voronoi_diagram<double>::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<Voronoi::Internal::point_type> 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<typename T> +T dist(const boost::polygon::point_data<T> &p1,const boost::polygon::point_data<T> &p2) +{ + T dx = p2.x() - p1.x(); + T dy = p2.y() - p1.y(); + return sqrt(dx*dx+dy*dy); +} + +// Find a foot point of "px" on a segment "seg". +template<typename segment_type, typename point_type> +inline point_type project_point_to_segment(segment_type &seg, point_type &px) +{ + typedef typename point_type::coordinate_type T; + const point_type &p0 = low(seg); + const point_type &p1 = high(seg); + const point_type dir(p1.x()-p0.x(), p1.y()-p0.y()); + const point_type dproj(px.x()-p0.x(), px.y()-p0.y()); + const T t = (dir.x()*dproj.x() + dir.y()*dproj.y()) / (dir.x()*dir.x() + dir.y()*dir.y()); + assert(t >= T(-1e-6) && t <= T(1. + 1e-6)); + return point_type(p0.x() + t*dir.x(), p0.y() + t*dir.y()); +} + +template<typename VD, typename SEGMENTS> +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<typename VD, typename SEGMENTS> +inline std::pair<typename VD::coord_type, typename VD::coord_type> +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, T>(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<VD>(cell2, segments); + return std::pair<T, T>(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<VD>(cell1, segments); + return std::pair<T, T>(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<VD>(cell1, segments); + return std::pair<T, T>(T(2.)*dist(pa, p1), T(2.)*dist(pb, p1)); + } +} + +// Converts the Line instances of Lines vector to VD::segment_type. +template<typename VD> +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.x), typename VD::coord_type(lines[idx].a.y)), + typename VD::point_type(typename VD::coord_type(lines[idx].b.x), typename VD::coord_type(lines[idx].b.y))); + } +private: + const Lines &lines; +}; + void MedialAxis::build(ThickPolylines* polylines) { @@ -373,6 +835,25 @@ MedialAxis::build(ThickPolylines* polylines) // append polyline to result polylines->push_back(polyline); } + + #ifdef SLIC3R_DEBUG + { + char path[2048]; + static int iRun = 0; + sprintf(path, "out/MedialAxis-%d.svg", iRun ++); + dump_voronoi_to_svg(this->lines, this->vd, polylines, path); + + + 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 |