Refactored Point / MultiPoint / Polyline / Polygon:

1) Removed virtual methods. There was not really need for them.
2) Some of the virtual methods were using conversion to Lines, which
   was unnecessary and expensive.
3) Removed some nearest element search methods from Point.

src/libslic3r/Arachne/SkeletalTrapezoidation.cpp

@@ -274,7 +274,7 @@ std::vector<Point> SkeletalTrapezoidation::discretize(const vd_t::edge_type& vd_
         Point right_point = VoronoiUtils::getSourcePoint(*right_cell, segments);
         coord_t d = (right_point - left_point).cast<int64_t>().norm();
         Point middle = (left_point + right_point) / 2;
-        Point x_axis_dir = Point(right_point - left_point).rotate_90_degree_ccw();
+        Point x_axis_dir = perp(Point(right_point - left_point));
         coord_t x_axis_length = x_axis_dir.cast<int64_t>().norm();
 
         const auto projected_x = [x_axis_dir, x_axis_length, middle](Point from) //Project a point on the edge.

src/libslic3r/Arachne/utils/VoronoiUtils.cpp

@@ -165,7 +165,7 @@ std::vector<Point> VoronoiUtils::discretizeParabola(const Point& p, const Segmen
     Line(a, b).distance_to_infinite_squared(p, &pxx);
     const Point ppxx = pxx - p;
     const coord_t d = ppxx.cast<int64_t>().norm();
-    const PointMatrix rot = PointMatrix(ppxx.rotate_90_degree_ccw());
+    const PointMatrix rot = PointMatrix(perp(ppxx));
 
     if (d == 0)
     {

src/libslic3r/Arachne/utils/linearAlg2D.hpp

@@ -38,15 +38,11 @@ inline static bool isInsideCorner(const Point &a, const Point &b, const Point &c
         return (p0.cast<int64_t>() * int64_t(len) / _len).cast<coord_t>();
     };
 
-    auto rotate_90_degree_ccw = [](const Vec2d &p) -> Vec2d {
-        return {-p.y(), p.x()};
-    };
-
     constexpr coord_t normal_length = 10000; //Create a normal vector of reasonable length in order to reduce rounding error.
     const Point ba = normal(a - b, normal_length);
     const Point bc = normal(c - b, normal_length);
     const Vec2d bq = query_point.cast<double>() - b.cast<double>();
-    const Vec2d perpendicular = rotate_90_degree_ccw(bq); //The query projects to this perpendicular to coordinate 0.
+    const Vec2d perpendicular = perp(bq); //The query projects to this perpendicular to coordinate 0.
 
     const double project_a_perpendicular = ba.cast<double>().dot(perpendicular); //Project vertex A on the perpendicular line.
     const double project_c_perpendicular = bc.cast<double>().dot(perpendicular); //Project vertex C on the perpendicular line.

src/libslic3r/ExPolygon.cpp

@@ -98,30 +98,24 @@ bool ExPolygon::contains(const Polylines &polylines) const
     return pl_out.empty();
 }
 
-bool ExPolygon::contains(const Point &point) const
+bool ExPolygon::contains(const Point &point, bool border_result /* = true */) const
 {
-    if (! Slic3r::contains(contour, point, true))
+    if (! Slic3r::contains(contour, point, border_result))
         // Outside the outer contour, not on the contour boundary.
         return false;
     for (const Polygon &hole : this->holes)
-        if (Slic3r::contains(hole, point, false))
+        if (Slic3r::contains(hole, point, ! border_result))
             // Inside a hole, not on the hole boundary.
             return false;
     return true;
 }
 
-// inclusive version of contains() that also checks whether point is on boundaries
+bool ExPolygon::on_boundary(const Point &point, double eps) const
-bool ExPolygon::contains_b(const Point &point) const
 {
-    return this->contains(point) || this->has_boundary_point(point);
+    if (this->contour.on_boundary(point, eps))
-}
-
-bool ExPolygon::has_boundary_point(const Point &point) const
-{
-    if (this->contour.has_boundary_point(point))
         return true;
     for (const Polygon &hole : this->holes)
-        if (hole.has_boundary_point(point))
+        if (hole.on_boundary(point, eps))
             return true;
     return false;
 }
@@ -148,6 +142,9 @@ Point ExPolygon::point_projection(const Point &point) const
 
 bool ExPolygon::overlaps(const ExPolygon &other) const
 {
+    if (this->empty() || other.empty())
+        return false;
+
     #if 0
     BoundingBox bbox = get_extents(other);
     bbox.merge(get_extents(*this));
@@ -164,7 +161,7 @@ bool ExPolygon::overlaps(const ExPolygon &other) const
     if (! pl_out.empty())
         return true; 
     //FIXME ExPolygon::overlaps() shall be commutative, it is not!
-    return ! other.contour.points.empty() && this->contains_b(other.contour.points.front());
+    return this->contains(other.contour.points.front());
 }
 
 void ExPolygon::simplify_p(double tolerance, Polygons* polygons) const
@@ -241,7 +238,7 @@ ExPolygon::medial_axis(double max_width, double min_width, ThickPolylines* polyl
            call, so we keep the inner point until we perform the second intersection() as well */
         Point new_front = polyline.points.front();
         Point new_back  = polyline.points.back();
-        if (polyline.endpoints.first && !this->has_boundary_point(new_front)) {
+        if (polyline.endpoints.first && !this->on_boundary(new_front, SCALED_EPSILON)) {
             Vec2d p1 = polyline.points.front().cast<double>();
             Vec2d p2 = polyline.points[1].cast<double>();
             // prevent the line from touching on the other side, otherwise intersection() might return that solution
@@ -251,7 +248,7 @@ ExPolygon::medial_axis(double max_width, double min_width, ThickPolylines* polyl
             p1 -= (p2 - p1).normalized() * max_width;
             this->contour.intersection(Line(p1.cast<coord_t>(), p2.cast<coord_t>()), &new_front);
         }
-        if (polyline.endpoints.second && !this->has_boundary_point(new_back)) {
+        if (polyline.endpoints.second && !this->on_boundary(new_back, SCALED_EPSILON)) {
             Vec2d p1 = (polyline.points.end() - 2)->cast<double>();
             Vec2d p2 = polyline.points.back().cast<double>();
             // prevent the line from touching on the other side, otherwise intersection() might return that solution

src/libslic3r/ExPolygon.hpp

@@ -51,9 +51,9 @@ public:
     bool contains(const Line &line) const;
     bool contains(const Polyline &polyline) const;
     bool contains(const Polylines &polylines) const;
-    bool contains(const Point &point) const;
+    bool contains(const Point &point, bool border_result = true) const;
-    bool contains_b(const Point &point) const;
+    // Approximate on boundary test.
-    bool has_boundary_point(const Point &point) const;
+    bool on_boundary(const Point &point, double eps) const;
     // Projection of a point onto the polygon.
     Point point_projection(const Point &point) const;
 
@@ -377,10 +377,10 @@ inline void expolygons_rotate(ExPolygons &expolys, double angle)
         expoly.rotate(angle);
 }
 
-inline bool expolygons_contain(ExPolygons &expolys, const Point &pt)
+inline bool expolygons_contain(ExPolygons &expolys, const Point &pt, bool border_result = true)
 {
     for (const ExPolygon &expoly : expolys)
-        if (expoly.contains(pt))
+        if (expoly.contains(pt, border_result))
             return true;
     return false;
 }

src/libslic3r/Fill/FillConcentric.cpp

@@ -40,7 +40,7 @@ void FillConcentric::_fill_surface_single(
     size_t iPathFirst = polylines_out.size();
     Point last_pos(0, 0);
     for (const Polygon &loop : loops) {
-        polylines_out.emplace_back(loop.split_at_index(last_pos.nearest_point_index(loop.points)));
+        polylines_out.emplace_back(loop.split_at_index(nearest_point_index(loop.points, last_pos)));
         last_pos = polylines_out.back().last_point();
     }
 
@@ -100,7 +100,7 @@ void FillConcentric::_fill_surface_single(const FillParams              &params,
             if (extrusion->is_closed && thick_polyline.points.front() == thick_polyline.points.back() && thick_polyline.width.front() == thick_polyline.width.back()) {
                 thick_polyline.points.pop_back();
                 assert(thick_polyline.points.size() * 2 == thick_polyline.width.size());
-                int nearest_idx = last_pos.nearest_point_index(thick_polyline.points);
+                int nearest_idx = nearest_point_index(thick_polyline.points, last_pos);
                 std::rotate(thick_polyline.points.begin(), thick_polyline.points.begin() + nearest_idx, thick_polyline.points.end());
                 std::rotate(thick_polyline.width.begin(), thick_polyline.width.begin() + 2 * nearest_idx, thick_polyline.width.end());
                 thick_polyline.points.emplace_back(thick_polyline.points.front());

src/libslic3r/GCode.cpp

@@ -116,8 +116,7 @@ namespace Slic3r {
             Point pos = Point::new_scale(writer_pos(0), writer_pos(1));
 
             // find standby point
-            Point standby_point;
+            Point standby_point = nearest_point(this->standby_points, pos).first;
-            pos.nearest_point(this->standby_points, &standby_point);
 
             /*  We don't call gcodegen.travel_to() because we don't need retraction (it was already
                 triggered by the caller) nor avoid_crossing_perimeters and also because the coordinates

src/libslic3r/Line.hpp

@@ -54,11 +54,11 @@ double distance_to_squared(const L &line, const Vec<Dim<L>, Scalar<L>> &point, V
     // We find projection of this point onto the line.
     // It falls where t = [(this-a) . (b-a)] / |b-a|^2
     const double t = va.dot(v) / l2;
-    if (t < 0.0) {
+    if (t <= 0.0) {
         // beyond the 'a' end of the segment
         *nearest_point = get_a(line);
         return va.squaredNorm();
-    } else if (t > 1.0) {
+    } else if (t >= 1.0) {
         // beyond the 'b' end of the segment
         *nearest_point = get_b(line);
         return (point - get_b(line)).template cast<double>().squaredNorm();

src/libslic3r/MultiPoint.cpp

@@ -45,16 +45,6 @@ void MultiPoint::rotate(double angle, const Point &center)
     }
 }
 
-double MultiPoint::length() const
-{
-    const Lines& lines = this->lines();
-    double len = 0;
-    for (auto it = lines.cbegin(); it != lines.cend(); ++it) {
-        len += it->length();
-    }
-    return len;
-}
-
 int MultiPoint::find_point(const Point &point) const
 {
     for (const Point &pt : this->points)
@@ -81,12 +71,6 @@ int MultiPoint::find_point(const Point &point, double scaled_epsilon) const
     return dist2_min < eps2 ? idx_min : -1;
 }
 
-bool MultiPoint::has_boundary_point(const Point &point) const
-{
-    double dist = (point.projection_onto(*this) - point).cast<double>().norm();
-    return dist < SCALED_EPSILON;
-}
-
 BoundingBox MultiPoint::bounding_box() const
 {
     return BoundingBox(this->points);
@@ -119,49 +103,6 @@ bool MultiPoint::remove_duplicate_points()
     return false;
 }
 
-bool MultiPoint::intersection(const Line& line, Point* intersection) const
-{
-    Lines lines = this->lines();
-    for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) {
-        if (it->intersection(line, intersection)) return true;
-    }
-    return false;
-}
-
-bool MultiPoint::first_intersection(const Line& line, Point* intersection) const
-{
-    bool   found = false;
-    double dmin  = 0.;
-    for (const Line &l : this->lines()) {
-        Point ip;
-        if (l.intersection(line, &ip)) {
-            if (! found) {
-                found = true;
-                dmin = (line.a - ip).cast<double>().norm();
-                *intersection = ip;
-            } else {
-                double d = (line.a - ip).cast<double>().norm();
-                if (d < dmin) {
-                    dmin = d;
-                    *intersection = ip;
-                }
-            }
-        }
-    }
-    return found;
-}
-
-bool MultiPoint::intersections(const Line &line, Points *intersections) const
-{
-    size_t intersections_size = intersections->size();
-    for (const Line &polygon_line : this->lines()) {
-        Point intersection;
-        if (polygon_line.intersection(line, &intersection))
-            intersections->emplace_back(std::move(intersection));
-    }
-    return intersections->size() > intersections_size;
-}
-
 std::vector<Point> MultiPoint::_douglas_peucker(const std::vector<Point>& pts, const double tolerance)
 {
     std::vector<Point> result_pts;
@@ -363,14 +304,6 @@ void MultiPoint3::translate(const Point& vector)
     this->translate(vector(0), vector(1));
 }
 
-double MultiPoint3::length() const
-{
-    double len = 0.0;
-    for (const Line3& line : this->lines())
-        len += line.length();
-    return len;
-}
-
 BoundingBox3 MultiPoint3::bounding_box() const
 {
     return BoundingBox3(points);

src/libslic3r/MultiPoint.hpp

@@ -18,7 +18,6 @@ public:
     Points points;
     
     MultiPoint() = default;
-    virtual ~MultiPoint() = default;
     MultiPoint(const MultiPoint &other) : points(other.points) {}
     MultiPoint(MultiPoint &&other) : points(std::move(other.points)) {}
     MultiPoint(std::initializer_list<Point> list) : points(list) {}
@@ -37,11 +36,8 @@ public:
     const Point& front() const { return this->points.front(); }
     const Point& back() const { return this->points.back(); }
     const Point& first_point() const { return this->front(); }
-    virtual const Point& last_point() const = 0;
-    virtual Lines lines() const = 0;
     size_t size() const { return points.size(); }
     bool   empty() const { return points.empty(); }
-    double length() const;
     bool   is_valid() const { return this->points.size() >= 2; }
 
     // Return index of a polygon point exactly equal to point.
@@ -50,7 +46,6 @@ public:
     // Return index of the closest point to point closer than scaled_epsilon.
     // Return -1 if no such point exists.
     int  find_point(const Point &point, const double scaled_epsilon) const;
-    bool has_boundary_point(const Point &point) const;
     int  closest_point_index(const Point &point) const {
         int idx = -1;
         if (! this->points.empty()) {
@@ -86,10 +81,6 @@ public:
         }
     }
 
-    bool intersection(const Line& line, Point* intersection) const;
-    bool first_intersection(const Line& line, Point* intersection) const;
-    bool intersections(const Line &line, Points *intersections) const;
-
     static Points _douglas_peucker(const Points &points, const double tolerance);
     static Points visivalingam(const Points& pts, const double& tolerance);
 
@@ -110,8 +101,6 @@ public:
 
     void translate(double x, double y);
     void translate(const Point& vector);
-    virtual Lines3 lines() const = 0;
-    double length() const;
     bool is_valid() const { return this->points.size() >= 2; }
 
     BoundingBox3 bounding_box() const;

src/libslic3r/PerimeterGenerator.cpp

@@ -405,8 +405,8 @@ static ClipperLib_Z::Paths clip_extrusion(const ClipperLib_Z::Path &subject, con
                 Point       prev(subject.front().x(), subject.front().y());
                 for (auto it = std::next(subject.begin()); it != subject.end(); ++it) {
                     Point curr(it->x(), it->y());
-                    Point projected_pt = pt.projection_onto(Line(prev, curr));
+                    Point projected_pt;
-                    if (double dist_sqr = (projected_pt - pt).cast<double>().squaredNorm(); dist_sqr < dist_sqr_min) {
+                    if (double dist_sqr = line_alg::distance_to_squared(Line(prev, curr), pt, &projected_pt); dist_sqr < dist_sqr_min) {
                         dist_sqr_min     = dist_sqr;
                         projected_pt_min = projected_pt;
                         it_min           = std::prev(it);

src/libslic3r/Point.cpp

@@ -57,98 +57,6 @@ void Point::rotate(double angle, const Point &center)
     (*this)(1) = (coord_t)round( (double)center(1) + c * dy + s * dx );
 }
 
-int Point::nearest_point_index(const Points &points) const
-{
-    PointConstPtrs p;
-    p.reserve(points.size());
-    for (Points::const_iterator it = points.begin(); it != points.end(); ++it)
-        p.push_back(&*it);
-    return this->nearest_point_index(p);
-}
-
-int Point::nearest_point_index(const PointConstPtrs &points) const
-{
-    int idx = -1;
-    double distance = -1;  // double because long is limited to 2147483647 on some platforms and it's not enough
-    
-    for (PointConstPtrs::const_iterator it = points.begin(); it != points.end(); ++it) {
-        /* If the X distance of the candidate is > than the total distance of the
-           best previous candidate, we know we don't want it */
-        double d = sqr<double>((*this)(0) - (*it)->x());
-        if (distance != -1 && d > distance) continue;
-        
-        /* If the Y distance of the candidate is > than the total distance of the
-           best previous candidate, we know we don't want it */
-        d += sqr<double>((*this)(1) - (*it)->y());
-        if (distance != -1 && d > distance) continue;
-        
-        idx = it - points.begin();
-        distance = d;
-        
-        if (distance < EPSILON) break;
-    }
-    
-    return idx;
-}
-
-int Point::nearest_point_index(const PointPtrs &points) const
-{
-    PointConstPtrs p;
-    p.reserve(points.size());
-    for (PointPtrs::const_iterator it = points.begin(); it != points.end(); ++it)
-        p.push_back(*it);
-    return this->nearest_point_index(p);
-}
-
-bool Point::nearest_point(const Points &points, Point* point) const
-{
-    int idx = this->nearest_point_index(points);
-    if (idx == -1) return false;
-    *point = points.at(idx);
-    return true;
-}
-
-Point Point::projection_onto(const MultiPoint &poly) const
-{
-    Point running_projection = poly.first_point();
-    double running_min = (running_projection - *this).cast<double>().norm();
-    
-    Lines lines = poly.lines();
-    for (Lines::const_iterator line = lines.begin(); line != lines.end(); ++line) {
-        Point point_temp = this->projection_onto(*line);
-        if ((point_temp - *this).cast<double>().norm() < running_min) {
-	        running_projection = point_temp;
-	        running_min = (running_projection - *this).cast<double>().norm();
-        }
-    }
-    return running_projection;
-}
-
-Point Point::projection_onto(const Line &line) const
-{
-    if (line.a == line.b) return line.a;
-    
-    /*
-        (Ported from VisiLibity by Karl J. Obermeyer)
-        The projection of point_temp onto the line determined by
-        line_segment_temp can be represented as an affine combination
-        expressed in the form projection of
-        Point = theta*line_segment_temp.first + (1.0-theta)*line_segment_temp.second.
-        If theta is outside the interval [0,1], then one of the Line_Segment's endpoints
-        must be closest to calling Point.
-    */
-    double lx = (double)(line.b(0) - line.a(0));
-    double ly = (double)(line.b(1) - line.a(1));
-    double theta = ( (double)(line.b(0) - (*this)(0))*lx + (double)(line.b(1)- (*this)(1))*ly ) 
-          / ( sqr<double>(lx) + sqr<double>(ly) );
-    
-    if (0.0 <= theta && theta <= 1.0)
-        return (theta * line.a.cast<coordf_t>() + (1.0-theta) * line.b.cast<coordf_t>()).cast<coord_t>();
-    
-    // Else pick closest endpoint.
-    return ((line.a - *this).cast<double>().squaredNorm() < (line.b - *this).cast<double>().squaredNorm()) ? line.a : line.b;
-}
-
 bool has_duplicate_points(std::vector<Point> &&pts)
 {
     std::sort(pts.begin(), pts.end());
@@ -179,6 +87,29 @@ BoundingBoxf get_extents(const std::vector<Vec2d> &pts)
     return bbox;
 }
 
+int nearest_point_index(const Points &points, const Point &pt)
+{
+    int64_t distance = std::numeric_limits<int64_t>::max();
+    int     idx      = -1;
+
+    for (const Point &pt2 : points) {
+        // If the X distance of the candidate is > than the total distance of the
+        // best previous candidate, we know we don't want it.
+        int64_t d = sqr<int64_t>(pt2.x() - pt.x());
+        if (d < distance) {
+            // If the Y distance of the candidate is > than the total distance of the
+            // best previous candidate, we know we don't want it.
+            d += sqr<int64_t>(pt2.y() - pt.y());
+            if (d < distance) {
+                idx      = &pt2 - points.data();
+                distance = d;
+            }
+        }
+    }
+
+    return idx;
+}
+
 std::ostream& operator<<(std::ostream &stm, const Vec2d &pointf)
 {
     return stm << pointf(0) << "," << pointf(1);

src/libslic3r/Point.hpp

@@ -17,8 +17,6 @@ namespace Slic3r {
 
 class BoundingBox;
 class BoundingBoxf;
-class Line;
-class MultiPoint;
 class Point;
 using Vector = Point;
 
@@ -183,13 +181,6 @@ public:
     Point  rotated(double angle) const { Point res(*this); res.rotate(angle); return res; }
     Point  rotated(double cos_a, double sin_a) const { Point res(*this); res.rotate(cos_a, sin_a); return res; }
     Point  rotated(double angle, const Point &center) const { Point res(*this); res.rotate(angle, center); return res; }
-    Point  rotate_90_degree_ccw() const { return Point(-this->y(), this->x()); }
-    int    nearest_point_index(const Points &points) const;
-    int    nearest_point_index(const PointConstPtrs &points) const;
-    int    nearest_point_index(const PointPtrs &points) const;
-    bool   nearest_point(const Points &points, Point* point) const;
-    Point  projection_onto(const MultiPoint &poly) const;
-    Point  projection_onto(const Line &line) const;
 };
 
 inline bool operator<(const Point &l, const Point &r) 
@@ -242,6 +233,14 @@ BoundingBox get_extents(const Points &pts);
 BoundingBox get_extents(const std::vector<Points> &pts);
 BoundingBoxf get_extents(const std::vector<Vec2d> &pts);
 
+int nearest_point_index(const Points &points, const Point &pt);
+
+inline std::pair<Point, bool> nearest_point(const Points &points, const Point &pt)
+{
+    int idx = nearest_point_index(points, pt);
+    return idx == -1 ? std::make_pair(Point(), false) : std::make_pair(points[idx], true);
+}
+
 // Test for duplicate points in a vector of points.
 // The points are copied, sorted and checked for duplicates globally.
 bool        has_duplicate_points(std::vector<Point> &&pts);

src/libslic3r/Polygon.cpp

@@ -6,6 +6,17 @@
 
 namespace Slic3r {
 
+double Polygon::length() const
+{
+    double l = 0;
+    if (this->points.size() > 1) {
+        l = (this->points.back() - this->points.front()).cast<double>().norm();
+        for (size_t i = 1; i < this->points.size(); ++ i)
+            l += (this->points[i] - this->points[i - 1]).cast<double>().norm();
+    }
+    return l;
+}
+
 Lines Polygon::lines() const
 {
     return to_lines(*this);
@@ -88,12 +99,6 @@ void Polygon::douglas_peucker(double tolerance)
     this->points = std::move(p);
 }
 
-// Does an unoriented polygon contain a point?
-bool Polygon::contains(const Point &p) const
-{
-    return Slic3r::contains(*this, p, true);
-}
-
 // this only works on CCW polygons as CW will be ripped out by Clipper's simplify_polygons()
 Polygons Polygon::simplify(double tolerance) const
 {
@@ -150,6 +155,64 @@ Point Polygon::centroid() const
     return Point(Vec2d(c / (3. * area_sum)));
 }
 
+bool Polygon::intersection(const Line &line, Point *intersection) const
+{
+    if (this->points.size() < 2)
+        return false;
+    if (Line(this->points.front(), this->points.back()).intersection(line, intersection))
+        return true;
+    for (size_t i = 1; i < this->points.size(); ++ i)
+        if (Line(this->points[i - 1], this->points[i]).intersection(line, intersection))
+            return true;
+    return false;
+}
+
+bool Polygon::first_intersection(const Line& line, Point* intersection) const
+{
+    if (this->points.size() < 2)
+        return false;
+
+    bool   found = false;
+    double dmin  = 0.;
+    Line l(this->points.back(), this->points.front());
+    for (size_t i = 0; i < this->points.size(); ++ i) {
+        l.b = this->points[i];
+        Point ip;
+        if (l.intersection(line, &ip)) {
+            if (! found) {
+                found = true;
+                dmin = (line.a - ip).cast<double>().squaredNorm();
+                *intersection = ip;
+            } else {
+                double d = (line.a - ip).cast<double>().squaredNorm();
+                if (d < dmin) {
+                    dmin = d;
+                    *intersection = ip;
+                }
+            }
+        }
+        l.a = l.b;
+    }
+    return found;
+}
+
+bool Polygon::intersections(const Line &line, Points *intersections) const
+{
+    if (this->points.size() < 2)
+        return false;
+
+    size_t intersections_size = intersections->size();
+    Line l(this->points.back(), this->points.front());
+    for (size_t i = 0; i < this->points.size(); ++ i) {
+        l.b = this->points[i];
+        Point intersection;
+        if (l.intersection(line, &intersection))
+            intersections->emplace_back(std::move(intersection));
+        l.a = l.b;
+    }
+    return intersections->size() > intersections_size;
+}
+
 // Filter points from poly to the output with the help of FilterFn.
 // filter function receives two vectors:
 // v1: this_point - previous_point

src/libslic3r/Polygon.hpp

@@ -15,11 +15,14 @@ using Polygons          = std::vector<Polygon>;
 using PolygonPtrs       = std::vector<Polygon*>;
 using ConstPolygonPtrs  = std::vector<const Polygon*>;
 
+// Returns true if inside. Returns border_result if on boundary.
+bool contains(const Polygon& polygon, const Point& p, bool border_result = true);
+bool contains(const Polygons& polygons, const Point& p, bool border_result = true);
+
 class Polygon : public MultiPoint
 {
 public:
     Polygon() = default;
-    ~Polygon() override = default;
     explicit Polygon(const Points &points) : MultiPoint(points) {}
 	Polygon(std::initializer_list<Point> points) : MultiPoint(points) {}
     Polygon(const Polygon &other) : MultiPoint(other.points) {}
@@ -38,9 +41,10 @@ public:
     const Point& operator[](Points::size_type idx) const { return this->points[idx]; }
 
     // last point == first point for polygons
-    const Point& last_point() const override { return this->points.front(); }
+    const Point& last_point() const { return this->points.front(); }
 
-    Lines lines() const override;
+    double length() const;
+    Lines lines() const;
     Polyline split_at_vertex(const Point &point) const;
     // Split a closed polygon into an open polyline, with the split point duplicated at both ends.
     Polyline split_at_index(int index) const;
@@ -58,13 +62,21 @@ public:
     void douglas_peucker(double tolerance);
 
     // Does an unoriented polygon contain a point?
-    // Tested by counting intersections along a horizontal line.
+    bool contains(const Point &point) const { return Slic3r::contains(*this, point, true); }
-    bool contains(const Point &point) const;
+    // Approximate on boundary test.
+    bool on_boundary(const Point &point, double eps) const
+        { return (this->point_projection(point) - point).cast<double>().squaredNorm() < eps * eps; }
+
     Polygons simplify(double tolerance) const;
     void simplify(double tolerance, Polygons &polygons) const;
     void densify(float min_length, std::vector<float>* lengths = nullptr);
     void triangulate_convex(Polygons* polygons) const;
     Point centroid() const;
+
+    bool intersection(const Line& line, Point* intersection) const;
+    bool first_intersection(const Line& line, Point* intersection) const;
+    bool intersections(const Line &line, Points *intersections) const;
+
     // Considering CCW orientation of this polygon, find all convex resp. concave points
     // with the angle at the vertex larger than a threshold.
     // Zero angle_threshold means to accept all convex resp. concave points.
@@ -254,10 +266,6 @@ inline Polygons to_polygons(std::vector<Points> &&paths)
 // however their contours may be rotated.
 bool polygons_match(const Polygon &l, const Polygon &r);
 
-// Returns true if inside. Returns border_result if on boundary.
-bool contains(const Polygon& polygon, const Point& p, bool border_result = true);
-bool contains(const Polygons& polygons, const Point& p, bool border_result = true);
-
 Polygon make_circle(double radius, double error);
 Polygon make_circle_num_segments(double radius, size_t num_segments);
 

src/libslic3r/Polyline.cpp

@@ -19,6 +19,14 @@ const Point& Polyline::leftmost_point() const
     return *p;
 }
 
+double Polyline::length() const
+{
+    double l = 0;
+    for (size_t i = 1; i < this->points.size(); ++ i)
+        l += (this->points[i] - this->points[i - 1]).cast<double>().norm();
+    return l;
+}
+
 Lines Polyline::lines() const
 {
     Lines lines;
@@ -146,9 +154,8 @@ void Polyline::split_at(const Point &point, Polyline* p1, Polyline* p2) const
     auto  min_point_it = this->points.cbegin();
     Point prev         = this->points.front();
     for (auto it = this->points.cbegin() + 1; it != this->points.cend(); ++ it) {
-        Point proj = point.projection_onto(Line(prev, *it));
+        Point proj;
-        auto  d2 = (proj - point).cast<double>().squaredNorm();
+        if (double d2 = line_alg::distance_to_squared(Line(prev, *it), point, &proj); d2 < min_dist2) {
-        if (d2 < min_dist2) {
 	        min_dist2    = d2;
 	        min_point_it = it;
         }
@@ -235,9 +242,8 @@ std::pair<int, Point> foot_pt(const Points &polyline, const Point &pt)
     auto  it = polyline.begin();
     auto  it_proj = polyline.begin();
     for (++ it; it != polyline.end(); ++ it) {
-        Point  foot_pt = pt.projection_onto(Line(prev, *it));
+        Point foot_pt;
-        double d2 = (foot_pt - pt).cast<double>().squaredNorm();
+        if (double d2 = line_alg::distance_to_squared(Line(prev, *it), pt, &foot_pt); d2 < d2_min) {
-        if (d2 < d2_min) {
             d2_min      = d2;
             foot_pt_min = foot_pt;
             it_proj     = it;
@@ -286,6 +292,14 @@ void ThickPolyline::clip_end(double distance)
     assert(this->width.size() == (this->points.size() - 1) * 2);
 }
 
+double Polyline3::length() const
+{
+    double l = 0;
+    for (size_t i = 1; i < this->points.size(); ++ i)
+        l += (this->points[i] - this->points[i - 1]).cast<double>().norm();
+    return l;
+}
+
 Lines3 Polyline3::lines() const
 {
     Lines3 lines;

src/libslic3r/Polyline.hpp

@@ -17,7 +17,6 @@ typedef std::vector<ThickPolyline> ThickPolylines;
 class Polyline : public MultiPoint {
 public:
     Polyline() = default;
-    ~Polyline() override = default;
     Polyline(const Polyline &other) : MultiPoint(other.points) {}
     Polyline(Polyline &&other) : MultiPoint(std::move(other.points)) {}
     Polyline(std::initializer_list<Point> list) : MultiPoint(list) {}
@@ -64,11 +63,12 @@ public:
     Point& operator[](Points::size_type idx) { return this->points[idx]; }
     const Point& operator[](Points::size_type idx) const { return this->points[idx]; }
 
-    const Point& last_point() const override { return this->points.back(); }
+    double length() const;
+    const Point& last_point() const { return this->points.back(); }
     const Point& leftmost_point() const;
-    Lines lines() const override;
+    Lines lines() const;
 
-    virtual void clip_end(double distance);
+    void clip_end(double distance);
     void clip_start(double distance);
     void extend_end(double distance);
     void extend_start(double distance);
@@ -170,7 +170,7 @@ public:
         std::swap(this->endpoints.first, this->endpoints.second);
     }
 
-    void clip_end(double distance) override;
+    void clip_end(double distance);
 
     std::vector<coordf_t> width;
     std::pair<bool,bool>  endpoints;
@@ -191,7 +191,8 @@ inline ThickPolylines to_thick_polylines(Polylines &&polylines, const coordf_t w
 class Polyline3 : public MultiPoint3
 {
 public:
-    virtual Lines3 lines() const;
+    double length() const;
+    Lines3 lines() const;
 };
 
 typedef std::vector<Polyline3> Polylines3;

xs/t/03_point.t

@@ -4,7 +4,7 @@ use strict;
 use warnings;
 
 use Slic3r::XS;
-use Test::More tests => 21;
+use Test::More tests => 16;
 
 my $point = Slic3r::Point->new(10, 15);
 
@@ -59,22 +59,4 @@ is_deeply [ @$point2 ], [30, 15], 'translate';
     ok $p0->ccw($p1, $p2) < 0, 'ccw() does not overflow';
 }
 
-{
-    my $point = Slic3r::Point->new(15,15);
-    my $line = Slic3r::Line->new([10,10], [20,10]);
-    is_deeply $point->projection_onto_line($line)->pp, [15,10], 'project_onto_line';
-    
-    $point = Slic3r::Point->new(0, 15);
-    is_deeply $point->projection_onto_line($line)->pp, [10,10], 'project_onto_line';
-    
-    $point = Slic3r::Point->new(25, 15);
-    is_deeply $point->projection_onto_line($line)->pp, [20,10], 'project_onto_line';
-    
-    $point = Slic3r::Point->new(10,10);
-    is_deeply $point->projection_onto_line($line)->pp, [10,10], 'project_onto_line';
-    
-    $point = Slic3r::Point->new(12, 10);
-    is_deeply $point->projection_onto_line($line)->pp, [12,10], 'project_onto_line';
-}
-
 __END__

xs/xsp/Point.xsp

@@ -29,9 +29,8 @@
         %code{% (*THIS)(0) = val; %};
     void set_y(int val)
         %code{% (*THIS)(1) = val; %};
-    int nearest_point_index(Points points);
     Clone<Point> nearest_point(Points points)
-        %code{% Point p; THIS->nearest_point(points, &p); RETVAL = p; %};
+        %code{% RETVAL = nearest_point(points, *THIS).first; %};
     double distance_to(Point* point)
         %code{% RETVAL = (*point - *THIS).cast<double>().norm(); %};
     double distance_to_line(Line* line)
@@ -40,8 +39,6 @@
         %code{% RETVAL = line->perp_distance_to(*THIS); %};
     double ccw(Point* p1, Point* p2)
         %code{% RETVAL = cross2((*p1 - *THIS).cast<double>(), (*p2 - *p1).cast<double>()); %};
-    Point* projection_onto_line(Line* line)
-        %code{% RETVAL = new Point(THIS->projection_onto(*line)); %};
     Point* negative()
         %code{% RETVAL = new Point(- *THIS); %};
     std::string serialize() %code{% char buf[2048]; sprintf(buf, "%ld,%ld", (*THIS)(0), (*THIS)(1)); RETVAL = buf; %};