More work to get Boost.Polygon to compile

xs/src/BoundingBox.cpp

@@ -4,7 +4,7 @@
 namespace Slic3r {
 
 template <class PointClass>
-BoundingBoxBase<PointClass>::BoundingBoxBase(const std::vector<PointClass> points)
+BoundingBoxBase<PointClass>::BoundingBoxBase(const std::vector<PointClass> &points)
 {
     typename std::vector<PointClass>::const_iterator it = points.begin();
     this->min.x = this->max.x = it->x;
@@ -16,10 +16,10 @@ BoundingBoxBase<PointClass>::BoundingBoxBase(const std::vector<PointClass> point
         this->max.y = std::max(it->y, this->max.y);
     }
 }
-template BoundingBoxBase<Point>::BoundingBoxBase(const std::vector<Point> points);
+template BoundingBoxBase<Point>::BoundingBoxBase(const std::vector<Point> &points);
 
 template <class PointClass>
-BoundingBox3Base<PointClass>::BoundingBox3Base(const std::vector<PointClass> points)
+BoundingBox3Base<PointClass>::BoundingBox3Base(const std::vector<PointClass> &points)
     : BoundingBoxBase<PointClass>(points)
 {
     typename std::vector<PointClass>::const_iterator it = points.begin();
@@ -29,7 +29,17 @@ BoundingBox3Base<PointClass>::BoundingBox3Base(const std::vector<PointClass> poi
         this->max.z = std::max(it->z, this->max.z);
     }
 }
-template BoundingBox3Base<Pointf3>::BoundingBox3Base(const std::vector<Pointf3> points);
+template BoundingBox3Base<Pointf3>::BoundingBox3Base(const std::vector<Pointf3> &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

xs/src/BoundingBox.hpp

@@ -20,7 +20,7 @@ class BoundingBoxBase
     PointClass max;
     
     BoundingBoxBase() {};
-    BoundingBoxBase(const std::vector<PointClass> points);
+    BoundingBoxBase(const std::vector<PointClass> &points);
     void merge(const PointClass &point);
     void merge(const BoundingBoxBase<PointClass> &bb);
     void scale(double factor);
@@ -34,7 +34,7 @@ class BoundingBox3Base : public BoundingBoxBase<PointClass>
 {
     public:
     BoundingBox3Base() {};
-    BoundingBox3Base(const std::vector<PointClass> points);
+    BoundingBox3Base(const std::vector<PointClass> &points);
     void merge(const PointClass &point);
     void merge(const BoundingBox3Base<PointClass> &bb);
     PointClass size() const;
@@ -48,7 +48,8 @@ class BoundingBox : public BoundingBoxBase<Point>
     void polygon(Polygon* polygon) const;
     
     BoundingBox() {};
-    BoundingBox(const Points points) : BoundingBoxBase<Point>(points) {};
+    BoundingBox(const Points &points) : BoundingBoxBase<Point>(points) {};
+    BoundingBox(const Lines &lines);
 };
 
 /*
@@ -59,7 +60,7 @@ class BoundingBox3  : public BoundingBox3Base<Point3> {};
 class BoundingBoxf3 : public BoundingBox3Base<Pointf3> {
     public:
     BoundingBoxf3() {};
-    BoundingBoxf3(const std::vector<Pointf3> points) : BoundingBox3Base<Pointf3>(points) {};
+    BoundingBoxf3(const std::vector<Pointf3> &points) : BoundingBox3Base<Pointf3>(points) {};
 };
 
 }

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(&ma.lines);
-    this->contour.lines(&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<double> vd;
+    ma.build(polylines);
-    construct_voronoi(lines.begin(), lines.end(), &vd);
-    
-    // iterate through the diagram
-    int result = 0;
-    for (voronoi_diagram<double>::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);
     
     // clip segments to our expolygon area
     intersection(*polylines, *this, *polylines);

xs/src/Geometry.cpp

@@ -3,6 +3,9 @@
 #include <algorithm>
 #include <map>
 #include <vector>
+#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<double>::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<double>::edge_type& edge, Points* clipped_edge)
+{
+    const voronoi_diagram<double>::cell_type& cell1 = *edge.cell();
+    const voronoi_diagram<double>::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<double>::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<coord_t>::discretize(point, segment, max_dist, sampled_edge);
+}
+
+Point
+MedialAxis::retrieve_point(const voronoi_diagram<double>::cell_type& cell)
+{
+    voronoi_diagram<double>::cell_type::source_index_type index = cell.source_index();
+    voronoi_diagram<double>::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<double>::cell_type& cell)
+{
+    voronoi_diagram<double>::cell_type::source_index_type index = cell.source_index() - this->points.size();
+    return this->lines[index];
+}
+
 } }

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<Points::size_type> &retval, Point
 void chained_path(Points &points, std::vector<Points::size_type> &retval);
 template<class T> 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<double>::edge_type& edge, Points* clipped_edge);
+    void sample_curved_edge(const voronoi_diagram<double>::edge_type& edge, Points* sampled_edge);
+    Point retrieve_point(const voronoi_diagram<double>::cell_type& cell);
+    Line retrieve_segment(const voronoi_diagram<double>::cell_type& cell);
+    
+    private:
+    voronoi_diagram<double> vd;
+    BoundingBox bb;
+};
+
 } }
 
 #endif

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)
 {

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);

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 <stack>
+#include <vector>
+
+#include <boost/polygon/isotropy.hpp>
+#include <boost/polygon/point_concept.hpp>
+#include <boost/polygon/segment_concept.hpp>
+#include <boost/polygon/rectangle_concept.hpp>
+
+namespace boost {
+namespace polygon {
+// Utilities class, that contains set of routines handful for visualization.
+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);
+  }
+};
+}
+}
+
+#endif  // BOOST_POLYGON_VORONOI_VISUAL_UTILS