From f9642786d3edce55698914856436f2e06957db04 Mon Sep 17 00:00:00 2001 From: Alessandro Ranellucci Date: Thu, 9 Jan 2014 19:56:12 +0100 Subject: [PATCH] More work to get Boost.Polygon to compile --- xs/src/BoundingBox.cpp | 18 +++- xs/src/BoundingBox.hpp | 9 +- xs/src/ExPolygon.cpp | 35 ++---- xs/src/Geometry.cpp | 122 +++++++++++++++++++++ xs/src/Geometry.hpp | 20 ++++ xs/src/Point.cpp | 5 + xs/src/Point.hpp | 1 + xs/src/voronoi_visual_utils.hpp | 186 ++++++++++++++++++++++++++++++++ 8 files changed, 360 insertions(+), 36 deletions(-) create mode 100644 xs/src/voronoi_visual_utils.hpp diff --git a/xs/src/BoundingBox.cpp b/xs/src/BoundingBox.cpp index f268ea1b9..70d688e50 100644 --- a/xs/src/BoundingBox.cpp +++ b/xs/src/BoundingBox.cpp @@ -4,7 +4,7 @@ namespace Slic3r { template -BoundingBoxBase::BoundingBoxBase(const std::vector points) +BoundingBoxBase::BoundingBoxBase(const std::vector &points) { typename std::vector::const_iterator it = points.begin(); this->min.x = this->max.x = it->x; @@ -16,10 +16,10 @@ BoundingBoxBase::BoundingBoxBase(const std::vector point this->max.y = std::max(it->y, this->max.y); } } -template BoundingBoxBase::BoundingBoxBase(const std::vector points); +template BoundingBoxBase::BoundingBoxBase(const std::vector &points); template -BoundingBox3Base::BoundingBox3Base(const std::vector points) +BoundingBox3Base::BoundingBox3Base(const std::vector &points) : BoundingBoxBase(points) { typename std::vector::const_iterator it = points.begin(); @@ -29,7 +29,17 @@ BoundingBox3Base::BoundingBox3Base(const std::vector poi this->max.z = std::max(it->z, this->max.z); } } -template BoundingBox3Base::BoundingBox3Base(const std::vector points); +template BoundingBox3Base::BoundingBox3Base(const std::vector &points); + +BoundingBox::BoundingBox(const Lines &lines) +{ + Points points; + for (Lines::const_iterator line = lines.begin(); line != lines.end(); ++line) { + points.push_back(line->a); + points.push_back(line->b); + } + *this = BoundingBox(points); +} void BoundingBox::polygon(Polygon* polygon) const diff --git a/xs/src/BoundingBox.hpp b/xs/src/BoundingBox.hpp index 8623b4909..92c2d063f 100644 --- a/xs/src/BoundingBox.hpp +++ b/xs/src/BoundingBox.hpp @@ -20,7 +20,7 @@ class BoundingBoxBase PointClass max; BoundingBoxBase() {}; - BoundingBoxBase(const std::vector points); + BoundingBoxBase(const std::vector &points); void merge(const PointClass &point); void merge(const BoundingBoxBase &bb); void scale(double factor); @@ -34,7 +34,7 @@ class BoundingBox3Base : public BoundingBoxBase { public: BoundingBox3Base() {}; - BoundingBox3Base(const std::vector points); + BoundingBox3Base(const std::vector &points); void merge(const PointClass &point); void merge(const BoundingBox3Base &bb); PointClass size() const; @@ -48,7 +48,8 @@ class BoundingBox : public BoundingBoxBase void polygon(Polygon* polygon) const; BoundingBox() {}; - BoundingBox(const Points points) : BoundingBoxBase(points) {}; + BoundingBox(const Points &points) : BoundingBoxBase(points) {}; + BoundingBox(const Lines &lines); }; /* @@ -59,7 +60,7 @@ class BoundingBox3 : public BoundingBox3Base {}; class BoundingBoxf3 : public BoundingBox3Base { public: BoundingBoxf3() {}; - BoundingBoxf3(const std::vector points) : BoundingBox3Base(points) {}; + BoundingBoxf3(const std::vector &points) : BoundingBox3Base(points) {}; }; } diff --git a/xs/src/ExPolygon.cpp b/xs/src/ExPolygon.cpp index 7191fc8e9..50ea5b094 100644 --- a/xs/src/ExPolygon.cpp +++ b/xs/src/ExPolygon.cpp @@ -1,11 +1,8 @@ #include "ExPolygon.hpp" +#include "Geometry.hpp" #include "Polygon.hpp" #include "Line.hpp" #include "ClipperUtils.hpp" -#include "boost/polygon/voronoi.hpp" - -using boost::polygon::voronoi_builder; -using boost::polygon::voronoi_diagram; namespace Slic3r { @@ -138,34 +135,16 @@ ExPolygon::simplify(double tolerance, ExPolygons &expolygons) const void ExPolygon::medial_axis(Polylines* polylines) const { + // init helper object + Slic3r::Geometry::MedialAxis ma; + // populate list of segments for the Voronoi diagram - Lines lines; - this->contour.lines(&lines); + this->contour.lines(&ma.lines); for (Polygons::const_iterator hole = this->holes.begin(); hole != this->holes.end(); ++hole) - hole->lines(&lines); + hole->lines(&ma.lines); // compute the Voronoi diagram - voronoi_diagram vd; - construct_voronoi(lines.begin(), lines.end(), &vd); - - // iterate through the diagram - int result = 0; - for (voronoi_diagram::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) { - if (it->is_primary()) ++result; - - Polyline p; - if (!it->is_finite()) { - clip_infinite_edge(*it, &p.points); - } else { - p.points.push_back(Point( it->vertex0()->x(), it->vertex0()->y() )); - p.points.push_back(Point( it->vertex1()->x(), it->vertex1()->y() )); - if (it->is_curved()) { - sample_curved_edge(*it, &p.points); - } - } - polylines->push_back(p); - } - printf("medial axis result = %d\n", result); + ma.build(polylines); // clip segments to our expolygon area intersection(*polylines, *this, *polylines); diff --git a/xs/src/Geometry.cpp b/xs/src/Geometry.cpp index be7848a05..a90e83b50 100644 --- a/xs/src/Geometry.cpp +++ b/xs/src/Geometry.cpp @@ -3,6 +3,9 @@ #include #include #include +#include "voronoi_visual_utils.hpp" + +using namespace boost::polygon; // provides also high() and low() namespace Slic3r { namespace Geometry { @@ -82,4 +85,123 @@ chained_path_items(Points &points, T &items, T &retval) } template void chained_path_items(Points &points, ClipperLib::PolyNodes &items, ClipperLib::PolyNodes &retval); +void +MedialAxis::build(Polylines* polylines) +{ + // build bounding box (we use it for clipping infinite segments) + this->bb = BoundingBox(this->lines); + + construct_voronoi(this->lines.begin(), this->lines.end(), &this->vd); + + // iterate through the diagram + int result = 0; + for (voronoi_diagram::const_edge_iterator it = this->vd.edges().begin(); it != this->vd.edges().end(); ++it) { + if (it->is_primary()) ++result; + + Polyline p; + if (!it->is_finite()) { + this->clip_infinite_edge(*it, &p.points); + } else { + p.points.push_back(Point( it->vertex0()->x(), it->vertex0()->y() )); + p.points.push_back(Point( it->vertex1()->x(), it->vertex1()->y() )); + if (it->is_curved()) { + this->sample_curved_edge(*it, &p.points); + } + } + polylines->push_back(p); + } + printf("medial axis result = %d\n", result); +} + +void +MedialAxis::clip_infinite_edge(const voronoi_diagram::edge_type& edge, Points* clipped_edge) +{ + const voronoi_diagram::cell_type& cell1 = *edge.cell(); + const voronoi_diagram::cell_type& cell2 = *edge.twin()->cell(); + Point origin, direction; + // Infinite edges could not be created by two segment sites. + if (cell1.contains_point() && cell2.contains_point()) { + Point p1 = retrieve_point(cell1); + Point p2 = retrieve_point(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(cell2) + : retrieve_point(cell1); + Line segment = cell1.contains_segment() + ? retrieve_segment(cell1) + : retrieve_segment(cell2); + coord_t dx = high(segment).x - low(segment).x; + coord_t 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; + } + } + coord_t side = this->bb.size().x; + coord_t koef = side / (std::max)(fabs(direction.x), fabs(direction.y)); + if (edge.vertex0() == NULL) { + clipped_edge->push_back(Point( + origin.x - direction.x * koef, + origin.y - direction.y * koef + )); + } else { + clipped_edge->push_back( + Point(edge.vertex0()->x(), edge.vertex0()->y())); + } + if (edge.vertex1() == NULL) { + clipped_edge->push_back(Point( + origin.x + direction.x * koef, + origin.y + direction.y * koef + )); + } else { + clipped_edge->push_back( + Point(edge.vertex1()->x(), edge.vertex1()->y())); + } +} + +void +MedialAxis::sample_curved_edge(const voronoi_diagram::edge_type& edge, Points* sampled_edge) +{ + Point point = edge.cell()->contains_point() + ? retrieve_point(*edge.cell()) + : retrieve_point(*edge.twin()->cell()); + + Line segment = edge.cell()->contains_point() + ? retrieve_segment(*edge.twin()->cell()) + : retrieve_segment(*edge.cell()); + + coord_t max_dist = 1E-3 * this->bb.size().x; + voronoi_visual_utils::discretize(point, segment, max_dist, sampled_edge); +} + +Point +MedialAxis::retrieve_point(const voronoi_diagram::cell_type& cell) +{ + voronoi_diagram::cell_type::source_index_type index = cell.source_index(); + voronoi_diagram::cell_type::source_category_type category = cell.source_category(); + if (category == SOURCE_CATEGORY_SINGLE_POINT) { + return this->points[index]; + } + index -= this->points.size(); + if (category == SOURCE_CATEGORY_SEGMENT_START_POINT) { + return low(this->lines[index]); + } else { + return high(this->lines[index]); + } +} + +Line +MedialAxis::retrieve_segment(const voronoi_diagram::cell_type& cell) +{ + voronoi_diagram::cell_type::source_index_type index = cell.source_index() - this->points.size(); + return this->lines[index]; +} + } } diff --git a/xs/src/Geometry.hpp b/xs/src/Geometry.hpp index ea6cb3242..fba942be2 100644 --- a/xs/src/Geometry.hpp +++ b/xs/src/Geometry.hpp @@ -1,8 +1,13 @@ #ifndef slic3r_Geometry_hpp_ #define slic3r_Geometry_hpp_ +#include "BoundingBox.hpp" #include "Polygon.hpp" +#include "boost/polygon/voronoi.hpp" +using boost::polygon::voronoi_builder; +using boost::polygon::voronoi_diagram; + namespace Slic3r { namespace Geometry { void convex_hull(Points &points, Polygon* hull); @@ -10,6 +15,21 @@ void chained_path(Points &points, std::vector &retval, Point void chained_path(Points &points, std::vector &retval); template void chained_path_items(Points &points, T &items, T &retval); +class MedialAxis { + public: + Points points; + Lines lines; + void build(Polylines* polylines); + void clip_infinite_edge(const voronoi_diagram::edge_type& edge, Points* clipped_edge); + void sample_curved_edge(const voronoi_diagram::edge_type& edge, Points* sampled_edge); + Point retrieve_point(const voronoi_diagram::cell_type& cell); + Line retrieve_segment(const voronoi_diagram::cell_type& cell); + + private: + voronoi_diagram vd; + BoundingBox bb; +}; + } } #endif diff --git a/xs/src/Point.cpp b/xs/src/Point.cpp index b6742a9d2..d951acb8a 100644 --- a/xs/src/Point.cpp +++ b/xs/src/Point.cpp @@ -4,6 +4,11 @@ namespace Slic3r { +inline bool +Point::operator==(const Point& rhs) const { + return this->coincides_with(rhs); +} + void Point::scale(double factor) { diff --git a/xs/src/Point.hpp b/xs/src/Point.hpp index f8b2514a9..20f1dcfb9 100644 --- a/xs/src/Point.hpp +++ b/xs/src/Point.hpp @@ -21,6 +21,7 @@ class Point coord_t x; coord_t y; explicit Point(coord_t _x = 0, coord_t _y = 0): x(_x), y(_y) {}; + bool operator==(const Point& rhs) const; void scale(double factor); void translate(double x, double y); void rotate(double angle, Point* center); diff --git a/xs/src/voronoi_visual_utils.hpp b/xs/src/voronoi_visual_utils.hpp new file mode 100644 index 000000000..4ef2e9dc2 --- /dev/null +++ b/xs/src/voronoi_visual_utils.hpp @@ -0,0 +1,186 @@ +// 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) + +// See http://www.boost.org for updates, documentation, and revision history. + +#ifndef BOOST_POLYGON_VORONOI_VISUAL_UTILS +#define BOOST_POLYGON_VORONOI_VISUAL_UTILS + +#include +#include + +#include +#include +#include +#include + +namespace boost { +namespace polygon { +// Utilities class, that contains set of routines handful for visualization. +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); + } +}; +} +} + +#endif // BOOST_POLYGON_VORONOI_VISUAL_UTILS