Removed various Point::ccw() and Point::ccw_angle() methods, they were

provided for Perl bindings and their semantic was confusing.
Implemented free function angle() to measure angle between two vectors.
Reworked Polygon::convex/concave_points(), changed the meaning of their
angle threshold parameter.
Removed some unused methods from Perl bindings and tests.
Reworked the "wipe inside at the external perimeter" function
after Point::ccw_angle() was removed.

src/libslic3r/GCode.cpp

@@ -2596,18 +2596,19 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
         // the side depends on the original winding order of the polygon (left for contours, right for holes)
         //FIXME improve the algorithm in case the loop is tiny.
         //FIXME improve the algorithm in case the loop is split into segments with a low number of points (see the Point b query).
-        Point a = paths.front().polyline.points[1];  // second point
-        Point b = *(paths.back().polyline.points.end()-3);       // second to last point
+        // Angle from the 2nd point to the last point.
+        double angle_inside = angle(paths.front().polyline.points[1]        - paths.front().first_point(),
+                                    *(paths.back().polyline.points.end()-3) - paths.front().first_point());
+        assert(angle_inside >= -M_PI && angle_inside <= M_PI);
+        // 3rd of this angle will be taken, thus make the angle monotonic before interpolation.
         if (was_clockwise) {
-            // swap points
-            Point c = a; a = b; b = c;
+            if (angle_inside > 0)
+                angle_inside -= 2.0 * M_PI;
+        } else {
+            if (angle_inside < 0)
+                angle_inside += 2.0 * M_PI;
         }
 
-        double angle = paths.front().first_point().ccw_angle(a, b) / 3;
-
-        // turn left if contour, turn right if hole
-        if (was_clockwise) angle *= -1;
-
         // create the destination point along the first segment and rotate it
         // we make sure we don't exceed the segment length because we don't know
         // the rotation of the second segment so we might cross the object boundary
@@ -2619,7 +2620,8 @@ std::string GCode::extrude_loop(ExtrusionLoop loop, std::string description, dou
         // Shift by no more than a nozzle diameter.
         //FIXME Hiding the seams will not work nicely for very densely discretized contours!
         Point  pt = ((nd * nd >= l2) ? p2 : (p1 + v * (nd / sqrt(l2)))).cast<coord_t>();
-        pt.rotate(angle, paths.front().polyline.points.front());
+        // Rotate pt inside around the seam point.
+        pt.rotate(angle_inside / 3., paths.front().polyline.points.front());
         // generate the travel move
         gcode += m_writer.travel_to_xy(this->point_to_gcode(pt), "move inwards before travel");
     }

src/libslic3r/Geometry.hpp

@@ -36,9 +36,9 @@ enum Orientation
 static inline Orientation orient(const Point &a, const Point &b, const Point &c)
 {
     static_assert(sizeof(coord_t) * 2 == sizeof(int64_t), "orient works with 32 bit coordinates");
-    int64_t u = int64_t(b(0)) * int64_t(c(1)) - int64_t(b(1)) * int64_t(c(0));
-    int64_t v = int64_t(a(0)) * int64_t(c(1)) - int64_t(a(1)) * int64_t(c(0));
-    int64_t w = int64_t(a(0)) * int64_t(b(1)) - int64_t(a(1)) * int64_t(b(0));
+    int64_t u = int64_t(b.x()) * int64_t(c.y()) - int64_t(b.y()) * int64_t(c.x());
+    int64_t v = int64_t(a.x()) * int64_t(c.y()) - int64_t(a.y()) * int64_t(c.x());
+    int64_t w = int64_t(a.x()) * int64_t(b.y()) - int64_t(a.y()) * int64_t(b.x());
     int64_t d = u - v + w;
     return (d > 0) ? ORIENTATION_CCW : ((d == 0) ? ORIENTATION_COLINEAR : ORIENTATION_CW);
 }

src/libslic3r/Geometry/ConvexHull.cpp

@@ -1,6 +1,7 @@
 #include "libslic3r.h"
 #include "ConvexHull.hpp"
 #include "BoundingBox.hpp"
+#include "../Geometry.hpp"
 
 #include <boost/multiprecision/integer.hpp>
 
@@ -19,13 +20,13 @@ Polygon convex_hull(Points pts)
         hull.points.resize(2 * n);
         // Build lower hull
         for (int i = 0; i < n; ++ i) {
-            while (k >= 2 && pts[i].ccw(hull[k-2], hull[k-1]) <= 0)
+            while (k >= 2 && Geometry::orient(pts[i], hull[k-2], hull[k-1]) != Geometry::ORIENTATION_CCW)
                 -- k;
             hull[k ++] = pts[i];
         }
         // Build upper hull
         for (int i = n-2, t = k+1; i >= 0; i--) {
-            while (k >= t && pts[i].ccw(hull[k-2], hull[k-1]) <= 0)
+            while (k >= t && Geometry::orient(pts[i], hull[k-2], hull[k-1]) != Geometry::ORIENTATION_CCW)
                 -- k;
             hull[k ++] = pts[i];
         }
@@ -58,7 +59,7 @@ Pointf3s convex_hull(Pointf3s points)
                 Point k1 = Point::new_scale(hull[k - 1](0), hull[k - 1](1));
                 Point k2 = Point::new_scale(hull[k - 2](0), hull[k - 2](1));
 
-                if (p.ccw(k2, k1) <= 0)
+                if (Geometry::orient(p, k2, k1) != Geometry::ORIENTATION_CCW)
                     --k;
                 else
                     break;
@@ -76,7 +77,7 @@ Pointf3s convex_hull(Pointf3s points)
                 Point k1 = Point::new_scale(hull[k - 1](0), hull[k - 1](1));
                 Point k2 = Point::new_scale(hull[k - 2](0), hull[k - 2](1));
 
-                if (p.ccw(k2, k1) <= 0)
+                if (Geometry::orient(p, k2, k1) != Geometry::ORIENTATION_CCW)
                     --k;
                 else
                     break;

src/libslic3r/Line.hpp

@@ -113,7 +113,6 @@ class Line
     Vector vector() const { return this->b - this->a; }
     Vector normal() const { return Vector((this->b(1) - this->a(1)), -(this->b(0) - this->a(0))); }
     bool   intersection(const Line& line, Point* intersection) const;
-    double ccw(const Point& point) const { return point.ccw(*this); }
     // Clip a line with a bounding box. Returns false if the line is completely outside of the bounding box.
 	bool   clip_with_bbox(const BoundingBox &bbox);
     // Extend the line from both sides by an offset.

src/libslic3r/Point.cpp

@@ -108,38 +108,6 @@ bool Point::nearest_point(const Points &points, Point* point) const
     return true;
 }
 
-/* Three points are a counter-clockwise turn if ccw > 0, clockwise if
- * ccw < 0, and collinear if ccw = 0 because ccw is a determinant that
- * gives the signed area of the triangle formed by p1, p2 and this point.
- * In other words it is the 2D cross product of p1-p2 and p1-this, i.e.
- * z-component of their 3D cross product.
- * We return double because it must be big enough to hold 2*max(|coordinate|)^2
- */
-double Point::ccw(const Point &p1, const Point &p2) const
-{
-    static_assert(sizeof(coord_t) == 4, "Point::ccw() requires a 32 bit coord_t");
-    return cross2((p2 - p1).cast<int64_t>(), (*this - p1).cast<int64_t>());
-//    return cross2((p2 - p1).cast<double>(), (*this - p1).cast<double>());
-}
-
-double Point::ccw(const Line &line) const
-{
-    return this->ccw(line.a, line.b);
-}
-
-// returns the CCW angle between this-p1 and this-p2
-// i.e. this assumes a CCW rotation from p1 to p2 around this
-double Point::ccw_angle(const Point &p1, const Point &p2) const
-{
-    const Point v1 = *this - p1;
-    const Point v2 = p2 - *this;
-    int64_t dot = int64_t(v1(0)) * int64_t(v2(0)) + int64_t(v1(1)) * int64_t(v2(1));
-    int64_t cross = int64_t(v1(0)) * int64_t(v2(1)) - int64_t(v1(1)) * int64_t(v2(0));
-    float angle = float(atan2(float(cross), float(dot)));
-    // we only want to return only positive angles
-    return angle <= 0 ? angle + 2*PI : angle;
-}
-
 Point Point::projection_onto(const MultiPoint &poly) const
 {
     Point running_projection = poly.first_point();

src/libslic3r/Point.hpp

@@ -79,24 +79,35 @@ inline const auto &identity3d = identity<3, double>;
 
 inline bool operator<(const Vec2d &lhs, const Vec2d &rhs) { return lhs.x() < rhs.x() || (lhs.x() == rhs.x() && lhs.y() < rhs.y()); }
 
-template<int Options>
-int32_t cross2(const Eigen::MatrixBase<Eigen::Matrix<int32_t, 2, 1, Options>> &v1, const Eigen::MatrixBase<Eigen::Matrix<int32_t, 2, 1, Options>> &v2) = delete;
-
-template<typename T, int Options>
-inline T cross2(const Eigen::MatrixBase<Eigen::Matrix<T, 2, 1, Options>> &v1, const Eigen::MatrixBase<Eigen::Matrix<T, 2, 1, Options>> &v2)
-{
-    return v1.x() * v2.y() - v1.y() * v2.x();
-}
-
+// Cross product of two 2D vectors.
+// None of the vectors may be of int32_t type as the result would overflow.
 template<typename Derived, typename Derived2>
 inline typename Derived::Scalar cross2(const Eigen::MatrixBase<Derived> &v1, const Eigen::MatrixBase<Derived2> &v2)
 {
+    static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "cross2(): first parameter is not a 2D vector");
+    static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "cross2(): first parameter is not a 2D vector");
+    static_assert(! std::is_same<typename Derived::Scalar, int32_t>::value, "cross2(): Scalar type must not be int32_t, otherwise the cross product would overflow.");
     static_assert(std::is_same<typename Derived::Scalar, typename Derived2::Scalar>::value, "cross2(): Scalar types of 1st and 2nd operand must be equal.");
     return v1.x() * v2.y() - v1.y() * v2.x();
 }
 
-template<typename T, int Options>
-inline Eigen::Matrix<T, 2, 1, Eigen::DontAlign> perp(const Eigen::MatrixBase<Eigen::Matrix<T, 2, 1, Options>> &v) { return Eigen::Matrix<T, 2, 1, Eigen::DontAlign>(- v.y(), v.x()); }
+// 2D vector perpendicular to the argument.
+template<typename Derived>
+inline Eigen::Matrix<typename Derived::Scalar, 2, 1, Eigen::DontAlign> perp(const Eigen::MatrixBase<Derived> &v)
+{ 
+    static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "perp(): parameter is not a 2D vector");
+    return { - v.y(), v.x() };
+}
+
+// Angle from v1 to v2, returning double atan2(y, x) normalized to <-PI, PI>.
+template<typename Derived, typename Derived2>
+inline double angle(const Eigen::MatrixBase<Derived> &v1, const Eigen::MatrixBase<Derived2> &v2) {
+    static_assert(Derived::IsVectorAtCompileTime && int(Derived::SizeAtCompileTime) == 2, "angle(): first parameter is not a 2D vector");
+    static_assert(Derived2::IsVectorAtCompileTime && int(Derived2::SizeAtCompileTime) == 2, "angle(): second parameter is not a 2D vector");
+    auto v1d = v1.cast<double>();
+    auto v2d = v2.cast<double>();
+    return atan2(cross2(v1d, v2d), v1d.dot(v2d));
+}
 
 template<class T, int N, int Options>
 Eigen::Matrix<T, 2, 1, Eigen::DontAlign> to_2d(const Eigen::MatrixBase<Eigen::Matrix<T, N, 1, Options>> &ptN) { return { ptN.x(), ptN.y() }; }
@@ -168,9 +179,6 @@ class Point : public Vec2crd
     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;
-    double ccw(const Point &p1, const Point &p2) const;
-    double ccw(const Line &line) const;
-    double ccw_angle(const Point &p1, const Point &p2) const;
     Point  projection_onto(const MultiPoint &poly) const;
     Point  projection_onto(const Line &line) const;
 };

src/libslic3r/Polygon.cpp

@@ -171,50 +171,56 @@ Point Polygon::centroid() const
     return Point(Vec2d(c / (3. * area_sum)));
 }
 
-// find all concave vertices (i.e. having an internal angle greater than the supplied angle)
-// (external = right side, thus we consider ccw orientation)
-Points Polygon::concave_points(double angle) const
-{
-    Points points;
-    angle = 2. * PI - angle + EPSILON;
-
-    // check whether first point forms a concave angle
-    if (this->points.front().ccw_angle(this->points.back(), *(this->points.begin()+1)) <= angle)
-        points.push_back(this->points.front());
-
-    // check whether points 1..(n-1) form concave angles
-    for (Points::const_iterator p = this->points.begin()+1; p != this->points.end()-1; ++ p)
-        if (p->ccw_angle(*(p-1), *(p+1)) <= angle)
-        	points.push_back(*p);
-
-    // check whether last point forms a concave angle
-    if (this->points.back().ccw_angle(*(this->points.end()-2), this->points.front()) <= angle)
-        points.push_back(this->points.back());
+// Filter points from poly to the output with the help of FilterFn.
+// filter function receives two vectors:
+// v1: this_point - previous_point
+// v2: next_point - this_point
+// and returns true if the point is to be copied to the output.
+template<typename FilterFn>
+Points filter_points_by_vectors(const Points &poly, FilterFn filter)
+{
+    // Last point is the first point visited.
+    Point p1 = poly.back();
+    // Previous vector to p1.
+    Vec2d v1 = (p1 - *(poly.end() - 2)).cast<double>();
+
+    Points out;
+    for (Point p2 : poly) {
+        // p2 is next point to the currently visited point p1.
+        Vec2d v2 = (p2 - p1).cast<double>();
+        if (filter(v1, v2))
+            out.emplace_back(p2);
+        v1 = v2;
+        p1 = p2;
+    }
 
-    return points;
+    return out;
 }
 
-// find all convex vertices (i.e. having an internal angle smaller than the supplied angle)
-// (external = right side, thus we consider ccw orientation)
-Points Polygon::convex_points(double angle) const
+template<typename ConvexConcaveFilterFn>
+Points filter_convex_concave_points_by_angle_threshold(const Points &poly, double angle_threshold, ConvexConcaveFilterFn convex_concave_filter)
 {
-    Points points;
-    angle = 2*PI - angle - EPSILON;
-
-    // check whether first point forms a convex angle
-    if (this->points.front().ccw_angle(this->points.back(), *(this->points.begin()+1)) >= angle)
-        points.push_back(this->points.front());
-
-    // check whether points 1..(n-1) form convex angles
-    for (Points::const_iterator p = this->points.begin()+1; p != this->points.end()-1; ++p) {
-        if (p->ccw_angle(*(p-1), *(p+1)) >= angle) points.push_back(*p);
+    assert(angle_threshold >= 0.);
+    if (angle_threshold < EPSILON) {
+        double cos_angle  = cos(angle_threshold);
+        return filter_points_by_vectors(poly, [convex_concave_filter, cos_angle](const Vec2d &v1, const Vec2d &v2){
+            return convex_concave_filter(v1, v2) && v1.normalized().dot(v2.normalized()) < cos_angle;
+        });
+    } else {
+        return filter_points_by_vectors(poly, [convex_concave_filter](const Vec2d &v1, const Vec2d &v2){
+            return convex_concave_filter(v1, v2);
+        });
     }
-
-    // check whether last point forms a convex angle
-    if (this->points.back().ccw_angle(*(this->points.end()-2), this->points.front()) >= angle)
-        points.push_back(this->points.back());
-
-    return points;
+}
+
+Points Polygon::convex_points(double angle_threshold) const
+{
+    return filter_convex_concave_points_by_angle_threshold(this->points, angle_threshold, [](const Vec2d &v1, const Vec2d &v2){ return cross2(v1, v2) > 0.; });
+}
+
+Points Polygon::concave_points(double angle_threshold) const
+{
+    return filter_convex_concave_points_by_angle_threshold(this->points, angle_threshold, [](const Vec2d &v1, const Vec2d &v2){ return cross2(v1, v2) < 0.; });
 }
 
 // Projection of a point onto the polygon.

src/libslic3r/Polygon.hpp

@@ -65,8 +65,11 @@ class Polygon : public MultiPoint
     void densify(float min_length, std::vector<float>* lengths = nullptr);
     void triangulate_convex(Polygons* polygons) const;
     Point centroid() const;
-    Points concave_points(double angle = PI) const;
-    Points convex_points(double angle = PI) 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.
+    Points convex_points(double angle_threshold = 0.) const;
+    Points concave_points(double angle_threshold = 0.) const;
     // Projection of a point onto the polygon.
     Point point_projection(const Point &point) const;
     std::vector<float> parameter_by_length() const;

t/geometry.t

@@ -2,7 +2,7 @@ use Test::More;
 use strict;
 use warnings;
 
-plan tests => 42;
+plan tests => 30;
 
 BEGIN {
     use FindBin;
@@ -189,63 +189,6 @@ my $polygons = [
     is +Slic3r::Point->new(10, 0)->distance_to_line($line), 0, 'distance_to';
 }
 
-#==========================================================
-
-{
-    my $square = Slic3r::Polygon->new_scale(
-        [100,100],
-        [200,100],
-        [200,200],
-        [100,200],
-    );
-    is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in ccw square';
-    is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in ccw square';
-
-    $square->make_clockwise;
-    is scalar(@{$square->concave_points(PI*4/3)}), 4, 'fuor concave vertices detected in cw square';
-    is scalar(@{$square->convex_points(PI*2/3)}), 0, 'no convex vertices detected in cw square';
-}
-
-{
-    my $square = Slic3r::Polygon->new_scale(
-        [150,100],
-        [200,100],
-        [200,200],
-        [100,200],
-        [100,100],
-    );
-    is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in convex polygon';
-    is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in square';
-}
-
-{
-    my $square = Slic3r::Polygon->new_scale(
-        [200,200],
-        [100,200],
-        [100,100],
-        [150,100],
-        [200,100],
-    );
-    is scalar(@{$square->concave_points(PI*4/3)}), 0, 'no concave vertices detected in convex polygon';
-    is scalar(@{$square->convex_points(PI*2/3)}), 4, 'four convex vertices detected in square';
-}
-
-{
-    my $triangle = Slic3r::Polygon->new(
-        [16000170,26257364], [714223,461012], [31286371,461008],
-    );
-    is scalar(@{$triangle->concave_points(PI*4/3)}), 0, 'no concave vertices detected in triangle';
-    is scalar(@{$triangle->convex_points(PI*2/3)}), 3, 'three convex vertices detected in triangle';
-}
-
-{
-    my $triangle = Slic3r::Polygon->new(
-        [16000170,26257364], [714223,461012], [20000000,461012], [31286371,461012],
-    );
-    is scalar(@{$triangle->concave_points(PI*4/3)}), 0, 'no concave vertices detected in triangle having collinear point';
-    is scalar(@{$triangle->convex_points(PI*2/3)}), 3, 'three convex vertices detected in triangle having collinear point';
-}
-
 {
     my $triangle = Slic3r::Polygon->new(
         [16000170,26257364], [714223,461012], [31286371,461008],