Meged with release_candidate_1_3

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