Promising approach to medial axis pruning

lib/Slic3r/Geometry.pm

@@ -5,7 +5,7 @@ use warnings;
 require Exporter;
 our @ISA = qw(Exporter);
 our @EXPORT_OK = qw(
-    PI X Y Z A B X1 Y1 X2 Y2 Z1 Z2 MIN MAX epsilon slope line_atan lines_parallel 
+    PI X Y Z A B X1 Y1 X2 Y2 Z1 Z2 MIN MAX epsilon slope 
     line_point_belongs_to_segment points_coincide distance_between_points 
     normalize tan move_points_3D
     point_in_polygon point_in_segment segment_in_segment
@@ -15,7 +15,7 @@ our @EXPORT_OK = qw(
     rotate_points move_points
     dot perp polygon_points_visibility
     line_intersection bounding_box bounding_box_intersect
-    angle3points three_points_aligned line_direction
+    angle3points
     chained_path chained_path_from collinear scale unscale
     rad2deg_dir bounding_box_center line_intersects_any douglas_peucker
     polyline_remove_short_segments normal triangle_normal polygon_is_convex
@@ -57,30 +57,6 @@ sub slope {
     return ($line->[B][Y] - $line->[A][Y]) / ($line->[B][X] - $line->[A][X]);
 }
 
-sub line_atan {
-    my ($line) = @_;
-    return atan2($line->[B][Y] - $line->[A][Y], $line->[B][X] - $line->[A][X]);
-}
-
-sub line_direction {
-    my ($line) = @_;
-    my $atan2 = line_atan($line);
-    return ($atan2 == PI) ? 0
-        : ($atan2 < 0) ? ($atan2 + PI)
-        : $atan2;
-}
-
-sub lines_parallel {
-    my ($line1, $line2) = @_;
-    
-    return abs(line_direction($line1) - line_direction($line2)) < $parallel_degrees_limit;
-}
-
-sub three_points_aligned {
-    my ($p1, $p2, $p3) = @_;
-    return lines_parallel([$p1, $p2], [$p2, $p3]);
-}
-
 # this subroutine checks whether a given point may belong to a given
 # segment given the hypothesis that it belongs to the line containing
 # the segment

lib/Slic3r/Layer/Region.pm

@@ -224,7 +224,7 @@ sub make_perimeters {
     # process thin walls by collapsing slices to single passes
     my $min_thin_wall_width = $pwidth/3;
     my $min_thin_wall_length = 2*$pwidth;
-    @thin_walls = @{offset2_ex([ map @$_, @thin_walls ], -0.5*$min_thin_wall_width, +0.5*$min_thin_wall_width)};
+    #@thin_walls = @{offset2_ex([ map @$_, @thin_walls ], -0.5*$min_thin_wall_width, +0.5*$min_thin_wall_width)};
     if (@thin_walls) {
         my @p = map @{$_->medial_axis($pspacing)}, @thin_walls;
         

lib/Slic3r/Line.pm

@@ -7,16 +7,6 @@ use parent 'Slic3r::Polyline';
 
 use Slic3r::Geometry qw(A B X Y);
 
-sub atan {
-    my $self = shift;
-    return Slic3r::Geometry::line_atan($self);
-}
-
-sub direction {
-    my $self = shift;
-    return Slic3r::Geometry::line_direction($self);
-}
-
 sub intersection {
     my $self = shift;
     my ($line, $require_crossing) = @_;

t/angles.t

@@ -10,7 +10,7 @@ BEGIN {
 }
 
 use Slic3r;
-use Slic3r::Geometry qw(line_atan line_direction rad2deg_dir angle3points PI);
+use Slic3r::Geometry qw(rad2deg_dir angle3points PI);
 
 #==========================================================
 
@@ -61,3 +61,13 @@ use Slic3r::Geometry qw(line_atan line_direction rad2deg_dir angle3points PI);
 }
 
 #==========================================================
+
+sub line_atan {
+    my ($l) = @_;
+    return Slic3r::Line->new(@$l)->atan2_;
+}
+
+sub line_direction {
+    my ($l) = @_;
+    return Slic3r::Line->new(@$l)->direction;
+}

xs/src/ExPolygon.cpp

@@ -135,10 +135,10 @@ ExPolygon::simplify(double tolerance, ExPolygons &expolygons) const
 }
 
 void
-ExPolygon::medial_axis(Polylines* polylines) const
+ExPolygon::medial_axis(double width, Polylines* polylines) const
 {
     // init helper object
-    Slic3r::Geometry::MedialAxis ma;
+    Slic3r::Geometry::MedialAxis ma(width);
     
     // populate list of segments for the Voronoi diagram
     this->contour.lines(&ma.lines);

xs/src/ExPolygon.hpp

@@ -26,7 +26,7 @@ class ExPolygon
     Polygons simplify_p(double tolerance) const;
     ExPolygons simplify(double tolerance) const;
     void simplify(double tolerance, ExPolygons &expolygons) const;
-    void medial_axis(Polylines* polylines) const;
+    void medial_axis(double width, Polylines* polylines) const;
     
     #ifdef SLIC3RXS
     void from_SV(SV* poly_sv);

xs/src/Geometry.cpp

@@ -3,10 +3,10 @@
 #include "PolylineCollection.hpp"
 #include "clipper.hpp"
 #include <algorithm>
+#include <list>
 #include <map>
 #include <set>
 #include <vector>
-//#include "voronoi_visual_utils.hpp"
 
 #ifdef SLIC3R_DEBUG
 #include "SVG.hpp"
@@ -92,6 +92,16 @@ chained_path_items(Points &points, T &items, T &retval)
 }
 template void chained_path_items(Points &points, ClipperLib::PolyNodes &items, ClipperLib::PolyNodes &retval);
 
+Line
+MedialAxis::edge_to_line(const VD::edge_type &edge) {
+    Line line;
+    line.a.x = edge.vertex0()->x();
+    line.a.y = edge.vertex0()->y();
+    line.b.x = edge.vertex1()->x();
+    line.b.y = edge.vertex1()->y();
+    return line;
+}
+
 void
 MedialAxis::build(Polylines* polylines)
 {
@@ -104,11 +114,73 @@ MedialAxis::build(Polylines* polylines)
     construct_voronoi(this->lines.begin(), this->lines.end(), &this->vd);
     
     // collect valid edges (i.e. prune those not belonging to MAT)
+    // note: this keeps twins, so it contains twice the number of the valid edges
     this->edges.clear();
     for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) {
         if (this->is_valid_edge(*edge)) this->edges.insert(&*edge);
     }
     
+    // count valid segments for each vertex
+    std::map< const VD::vertex_type*,std::list<const VD::edge_type*> > vertex_edges;
+    std::list<const VD::vertex_type*> entry_nodes;
+    for (VD::const_vertex_iterator vertex = this->vd.vertices().begin(); vertex != this->vd.vertices().end(); ++vertex) {
+        // get a reference to the list of valid edges originating from this vertex
+        std::list<const VD::edge_type*>& edges = vertex_edges[&*vertex];
+        
+        // get one random edge originating from this vertex
+        const VD::edge_type* edge = vertex->incident_edge();
+        do {
+            if (this->edges.count(edge) > 0)    // only count valid edges
+                edges.push_back(edge);
+            edge = edge->rot_next();            // next edge originating from this vertex
+        } while (edge != vertex->incident_edge());
+        
+        // if there's only one edge starting at this vertex then it's a leaf
+        if (edges.size() == 1) entry_nodes.push_back(&*vertex);
+    }
+    
+    // iterate through the leafs to prune short branches
+    for (std::list<const VD::vertex_type*>::const_iterator vertex = entry_nodes.begin(); vertex != entry_nodes.end(); ++vertex) {
+        const VD::vertex_type* v = *vertex;
+        
+        // start a polyline from this vertex
+        Polyline polyline;
+        polyline.points.push_back(Point(v->x(), v->y()));
+        
+        // keep track of visited edges to prevent infinite loops
+        std::set<const VD::edge_type*> visited_edges;
+        
+        do {
+            // get edge starting from v
+            const VD::edge_type* edge = vertex_edges[v].front();
+            
+            // if we picked the edge going backwards (thus the twin of the previous edge)
+            if (visited_edges.count(edge->twin()) > 0) {
+                edge = vertex_edges[v].back();
+            }
+            
+            // avoid getting twice on the same edge
+            if (visited_edges.count(edge) > 0) break;
+            visited_edges.insert(edge);
+            
+            // get ending vertex for this edge and append it to the polyline
+            v = edge->vertex1();
+            polyline.points.push_back(Point( v->x(), v->y() ));
+            
+            // if two edges start at this vertex (one forward one backwards) then
+            // it's not branching and we can go on
+        } while (vertex_edges[v].size() == 2);
+        
+        // if this branch is too short, invalidate all of its edges so that 
+        // they will be ignored when building actual polylines in the loop below
+        if (polyline.length() < this->width) {
+            for (std::set<const VD::edge_type*>::const_iterator edge = visited_edges.begin(); edge != visited_edges.end(); ++edge) {
+                (void)this->edges.erase(*edge);
+                (void)this->edges.erase((*edge)->twin());
+            }
+        }
+    }
+    
     // iterate through the valid edges to build polylines
     while (!this->edges.empty()) {
         const VD::edge_type& edge = **this->edges.begin();
@@ -176,107 +248,26 @@ MedialAxis::is_valid_edge(const VD::edge_type& edge) const
        but I don't know how to do it. Maybe we could check the relative angle of
        the two segments (we are only interested in facing segments). */
     
-    const voronoi_diagram<double>::cell_type &cell1 = *edge.cell();
+    const VD::cell_type &cell1 = *edge.cell();
-    const voronoi_diagram<double>::cell_type &cell2 = *edge.twin()->cell();
+    const VD::cell_type &cell2 = *edge.twin()->cell();
     if (cell1.contains_segment() && cell2.contains_segment()) {
         Line segment1 = this->retrieve_segment(cell1);
         Line segment2 = this->retrieve_segment(cell2);
         if (segment1.a == segment2.b || segment1.b == segment2.a) return false;
+        if (fabs(segment1.atan2_() - segment2.atan2_()) < PI/3) {
+            //printf("segment1 atan2 = %f, segment2 atan2 = %f\n", segment1.atan2_(), segment2.atan2_());
+            //printf("  => SAME ATAN2\n");
+            return false;
+        }
     }
     
     return true;
 }
 
-/*
-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());
-    
-    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
-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) const
+MedialAxis::retrieve_segment(const VD::cell_type& cell) const
 {
-    voronoi_diagram<double>::cell_type::source_index_type index = cell.source_index() - this->points.size();
+    VD::cell_type::source_index_type index = cell.source_index() - this->points.size();
     return this->lines[index];
 }
 

xs/src/Geometry.hpp

@@ -20,19 +20,18 @@ class MedialAxis {
     public:
     Points points;
     Lines lines;
+    double width;
+    MedialAxis(double _width) : width(_width) {};
     void build(Polylines* polylines);
-    void process_edge_neighbors(const voronoi_diagram<double>::edge_type& edge, Points* points);
-    bool is_valid_edge(const voronoi_diagram<double>::edge_type& edge) const;
-    //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) const;
     
     private:
     typedef voronoi_diagram<double> VD;
     VD vd;
-    //BoundingBox bb;
     std::set<const VD::edge_type*> edges;
+    Line edge_to_line(const VD::edge_type &edge);
+    void process_edge_neighbors(const voronoi_diagram<double>::edge_type& edge, Points* points);
+    bool is_valid_edge(const voronoi_diagram<double>::edge_type& edge) const;
+    Line retrieve_segment(const voronoi_diagram<double>::cell_type& cell) const;
 };
 
 } }

xs/src/Line.cpp

@@ -1,6 +1,7 @@
 #include "Line.hpp"
 #include "Polyline.hpp"
 #include <algorithm>
+#include <cmath>
 #include <sstream>
 
 namespace Slic3r {
@@ -85,6 +86,21 @@ Line::distance_to(const Point* point) const
     return point->distance_to(this);
 }
 
+double
+Line::atan2_() const
+{
+    return atan2(this->b.y - this->a.y, this->b.x - this->a.x);
+}
+
+double
+Line::direction() const
+{
+    double atan2 = this->atan2_();
+    return (atan2 == PI) ? 0
+        : (atan2 < 0) ? (atan2 + PI)
+        : atan2;
+}
+
 #ifdef SLIC3RXS
 void
 Line::from_SV(SV* line_sv)

xs/src/Line.hpp

@@ -28,6 +28,8 @@ class Line
     Point* point_at(double distance) const;
     bool coincides_with(const Line* line) const;
     double distance_to(const Point* point) const;
+    double atan2_() const;
+    double direction() const;
     
     #ifdef SLIC3RXS
     void from_SV(SV* line_sv);

xs/src/voronoi_visual_utils.hpp

@@ -1,186 +0,0 @@
-// 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,
-            class Point,
-            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<InCT1>((*discretization)[0], segment);
-    CT projection_end = sqr_segment_length *
-        get_point_projection<InCT1>((*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<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(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,
-            class Point,
-            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 <typename InCT>
-  static CT cast(const InCT& value) {
-    return static_cast<CT>(value);
-  }
-};
-}
-}
-
-#endif  // BOOST_POLYGON_VORONOI_VISUAL_UTILS

xs/xsp/ExPolygon.xsp

@@ -26,7 +26,7 @@
     ExPolygons simplify(double tolerance);
     Polygons simplify_p(double tolerance);
     Polylines medial_axis(double width)
-        %code{% THIS->medial_axis(&RETVAL); %};
+        %code{% THIS->medial_axis(width, &RETVAL); %};
 %{
 
 ExPolygon*

xs/xsp/Line.xsp

@@ -21,6 +21,8 @@
     void scale(double factor);
     void translate(double x, double y);
     double length();
+    double atan2_();
+    double direction();
     Point* midpoint()
          %code{% const char* CLASS = "Slic3r::Point"; RETVAL = THIS->midpoint(); %};
     Point* point_at(double distance)