Semi-working Boost.Polygon medial axis. Hangs after perimeters, though. Needs pruning and chaining

lib/Slic3r/Layer/Region.pm

@@ -224,6 +224,16 @@ sub make_perimeters {
     # process thin walls by collapsing slices to single passes
     if (@thin_walls) {
         my @p = map @{$_->medial_axis($pspacing)}, @thin_walls;
+        if (0) {
+            use Slic3r::SVG;
+            Slic3r::SVG::output(
+                "medial_axis.svg",
+                no_arrows => 1,
+                expolygons => \@thin_walls,
+                polylines => \@p,
+            );
+            exit;
+        }
         my @paths = ();
         for my $p (@p) {
             next if $p->length <= $pspacing * 2;

xs/src/Geometry.cpp

@@ -2,6 +2,7 @@
 #include "clipper.hpp"
 #include <algorithm>
 #include <map>
+#include <set>
 #include <vector>
 #include "voronoi_visual_utils.hpp"
 
@@ -93,10 +94,15 @@ MedialAxis::build(Polylines* polylines)
     
     construct_voronoi(this->lines.begin(), this->lines.end(), &this->vd);
     
+    // prepare a cache of twin edges to prevent getting the same edge twice
+    // (Boost.Polygon returns it duplicated in both directions)
+    std::set<const voronoi_diagram<double>::edge_type*> edge_cache;
+    
     // 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;
+        (void)edge_cache.insert(it->twin());
+        if (edge_cache.count(&*it) > 0) continue;
+        if (!it->is_primary()) continue;
         
         Polyline p;
         if (!it->is_finite()) {
@@ -110,7 +116,6 @@ MedialAxis::build(Polylines* polylines)
         }
         polylines->push_back(p);
     }
-    printf("medial axis result = %d\n", result);
 }
 
 void
@@ -177,8 +182,8 @@ MedialAxis::sample_curved_edge(const voronoi_diagram<double>::edge_type& edge, P
         ? 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);
+    double max_dist = 1E-3 * this->bb.size().x;
+    voronoi_visual_utils<double>::discretize<coord_t,coord_t,Point,Line>(point, segment, max_dist, sampled_edge);
 }
 
 Point

xs/src/Point.cpp

@@ -4,7 +4,7 @@
 
 namespace Slic3r {
 
-inline bool
+bool
 Point::operator==(const Point& rhs) const {
     return this->coincides_with(rhs);
 }

xs/src/voronoi_visual_utils.hpp

@@ -42,28 +42,28 @@ class voronoi_visual_utils {
   // Important:
   //   discretization should contain both edge endpoints initially.
   template <class InCT1, class InCT2,
-            template<class> class Point,
-            template<class> class Segment>
+            class Point,
+            class Segment>
   static
   typename enable_if<
     typename gtl_and<
       typename gtl_if<
         typename is_point_concept<
-          typename geometry_concept< Point<InCT1> >::type
+          typename geometry_concept< Point >::type
         >::type
       >::type,
       typename gtl_if<
         typename is_segment_concept<
-          typename geometry_concept< Segment<InCT2> >::type
+          typename geometry_concept< Segment >::type
         >::type
       >::type
     >::type,
     void
   >::type discretize(
-      const Point<InCT1>& point,
-      const Segment<InCT2>& segment,
+      const Point& point,
+      const Segment& segment,
       const CT max_dist,
-      std::vector< Point<CT> >* discretization) {
+      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.
@@ -74,9 +74,9 @@ class voronoi_visual_utils {
     // 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);
+        get_point_projection<InCT1>((*discretization)[0], segment);
     CT projection_end = sqr_segment_length *
-        get_point_projection((*discretization)[1], segment);
+        get_point_projection<InCT1>((*discretization)[1], segment);
 
     // Compute parabola parameters in the transformed space.
     // Parabola has next representation:
@@ -87,7 +87,7 @@ class voronoi_visual_utils {
     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];
+    Point last_point = (*discretization)[1];
     discretization->pop_back();
 
     // Use stack to avoid recursion.
@@ -120,7 +120,7 @@ class voronoi_visual_utils {
             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));
+        discretization->push_back(Point(inter_x, inter_y));
         cur_x = new_x;
         cur_y = new_y;
       } else {
@@ -146,25 +146,25 @@ class voronoi_visual_utils {
   // 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>
+            class Point,
+            class Segment>
   static
   typename enable_if<
     typename gtl_and<
       typename gtl_if<
         typename is_point_concept<
-          typename geometry_concept< Point<int> >::type
+          typename geometry_concept< Point >::type
         >::type
       >::type,
       typename gtl_if<
         typename is_segment_concept<
-          typename geometry_concept< Segment<long> >::type
+          typename geometry_concept< Segment >::type
         >::type
       >::type
     >::type,
     CT
   >::type get_point_projection(
-      const Point<CT>& point, const Segment<InCT>& segment) {
+      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)));