diff --git a/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp b/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp
index 04dede546..04a1042d8 100644
--- a/src/libslic3r/Arachne/SkeletalTrapezoidation.cpp
+++ b/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.
diff --git a/src/libslic3r/Arachne/utils/VoronoiUtils.cpp b/src/libslic3r/Arachne/utils/VoronoiUtils.cpp
index 82bd79523..069e1f5ad 100644
--- a/src/libslic3r/Arachne/utils/VoronoiUtils.cpp
+++ b/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)
     {
diff --git a/src/libslic3r/Arachne/utils/linearAlg2D.hpp b/src/libslic3r/Arachne/utils/linearAlg2D.hpp
index 797bae0b9..304984b1c 100644
--- a/src/libslic3r/Arachne/utils/linearAlg2D.hpp
+++ b/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.
diff --git a/src/libslic3r/ExPolygon.cpp b/src/libslic3r/ExPolygon.cpp
index 311194102..61c204044 100644
--- a/src/libslic3r/ExPolygon.cpp
+++ b/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::contains_b(const Point &point) const
+bool ExPolygon::on_boundary(const Point &point, double eps) const
 {
-    return this->contains(point) || this->has_boundary_point(point);
-}
-
-bool ExPolygon::has_boundary_point(const Point &point) const
-{
-    if (this->contour.has_boundary_point(point))
+    if (this->contour.on_boundary(point, eps))
         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
diff --git a/src/libslic3r/ExPolygon.hpp b/src/libslic3r/ExPolygon.hpp
index c0f3d95a4..58f0d33e9 100644
--- a/src/libslic3r/ExPolygon.hpp
+++ b/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_b(const Point &point) const;
-    bool has_boundary_point(const Point &point) const;
+    bool contains(const Point &point, bool border_result = true) const;
+    // Approximate on boundary test.
+    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;
 }
diff --git a/src/libslic3r/Fill/FillConcentric.cpp b/src/libslic3r/Fill/FillConcentric.cpp
index 69f530720..7b005ee35 100644
--- a/src/libslic3r/Fill/FillConcentric.cpp
+++ b/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());
diff --git a/src/libslic3r/GCode.cpp b/src/libslic3r/GCode.cpp
index f42ae4f83..12342ddca 100644
--- a/src/libslic3r/GCode.cpp
+++ b/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;
-            pos.nearest_point(this->standby_points, &standby_point);
+            Point standby_point = nearest_point(this->standby_points, pos).first;
 
             /*  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
diff --git a/src/libslic3r/Line.hpp b/src/libslic3r/Line.hpp
index f9a8977d5..f7a14c251 100644
--- a/src/libslic3r/Line.hpp
+++ b/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();
diff --git a/src/libslic3r/MultiPoint.cpp b/src/libslic3r/MultiPoint.cpp
index 5ed9eb23c..f18720bd6 100644
--- a/src/libslic3r/MultiPoint.cpp
+++ b/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);
diff --git a/src/libslic3r/MultiPoint.hpp b/src/libslic3r/MultiPoint.hpp
index bc9cf761d..778b30c57 100644
--- a/src/libslic3r/MultiPoint.hpp
+++ b/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;
diff --git a/src/libslic3r/PerimeterGenerator.cpp b/src/libslic3r/PerimeterGenerator.cpp
index 92759238e..e864a37ad 100644
--- a/src/libslic3r/PerimeterGenerator.cpp
+++ b/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));
-                    if (double dist_sqr = (projected_pt - pt).cast<double>().squaredNorm(); dist_sqr < dist_sqr_min) {
+                    Point projected_pt;
+                    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);
diff --git a/src/libslic3r/Point.cpp b/src/libslic3r/Point.cpp
index 9f56f96d6..a1145a952 100644
--- a/src/libslic3r/Point.cpp
+++ b/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);
diff --git a/src/libslic3r/Point.hpp b/src/libslic3r/Point.hpp
index 949ddbad1..3d0bae2c3 100644
--- a/src/libslic3r/Point.hpp
+++ b/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);
diff --git a/src/libslic3r/Polygon.cpp b/src/libslic3r/Polygon.cpp
index 26b326b1c..717e25aed 100644
--- a/src/libslic3r/Polygon.cpp
+++ b/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
diff --git a/src/libslic3r/Polygon.hpp b/src/libslic3r/Polygon.hpp
index ccfe213e8..8a4a2ba0a 100644
--- a/src/libslic3r/Polygon.hpp
+++ b/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;
+    bool contains(const Point &point) const { return Slic3r::contains(*this, point, true); }
+    // 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);
 
diff --git a/src/libslic3r/Polyline.cpp b/src/libslic3r/Polyline.cpp
index 42d89c882..e2dcabfe1 100644
--- a/src/libslic3r/Polyline.cpp
+++ b/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));
-        auto  d2 = (proj - point).cast<double>().squaredNorm();
-        if (d2 < min_dist2) {
+        Point proj;
+        if (double d2 = line_alg::distance_to_squared(Line(prev, *it), point, &proj); 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));
-        double d2 = (foot_pt - pt).cast<double>().squaredNorm();
-        if (d2 < d2_min) {
+        Point foot_pt;
+        if (double d2 = line_alg::distance_to_squared(Line(prev, *it), pt, &foot_pt); 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;
diff --git a/src/libslic3r/Polyline.hpp b/src/libslic3r/Polyline.hpp
index 547664d4f..bee5e51ba 100644
--- a/src/libslic3r/Polyline.hpp
+++ b/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;
diff --git a/xs/t/03_point.t b/xs/t/03_point.t
index f888349b3..389e120c7 100644
--- a/xs/t/03_point.t
+++ b/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__
diff --git a/xs/xsp/Point.xsp b/xs/xsp/Point.xsp
index 2d70e0203..0d44ea364 100644
--- a/xs/xsp/Point.xsp
+++ b/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; %};