Added test projects for libslic3r and fff_print.

Added test_geometry.cpp from upstream slic3r, thanks @lordofhyphens
Added circle_taubin_newton() for circle center calculation, thanks @lordofhyphens

src/libslic3r/Geometry.cpp

@@ -16,6 +16,7 @@
 
 #include <boost/algorithm/string/classification.hpp>
 #include <boost/algorithm/string/split.hpp>
+#include <boost/log/trivial.hpp>
 
 #ifdef SLIC3R_DEBUG
 #include "SVG.hpp"
@@ -335,6 +336,93 @@ double rad2deg_dir(double angle)
     return rad2deg(angle);
 }
 
+Point circle_taubin_newton(const Points::const_iterator& input_begin, const Points::const_iterator& input_end, size_t cycles)
+{
+    Vec2ds tmp;
+    tmp.reserve(std::distance(input_begin, input_end));
+    std::transform(input_begin, input_end, std::back_inserter(tmp), [] (const Point& in) { return unscale(in); } );
+    Vec2d center = circle_taubin_newton(tmp.cbegin(), tmp.end(), cycles);
+	return Point::new_scale(center.x(), center.y());
+}
+
+/// Adapted from work in "Circular and Linear Regression: Fitting circles and lines by least squares", pg 126
+/// Returns a point corresponding to the center of a circle for which all of the points from input_begin to input_end
+/// lie on.
+Vec2d circle_taubin_newton(const Vec2ds::const_iterator& input_begin, const Vec2ds::const_iterator& input_end, size_t cycles)
+{
+    // calculate the centroid of the data set
+    const Vec2d sum = std::accumulate(input_begin, input_end, Vec2d(0,0));
+    const size_t n = std::distance(input_begin, input_end);
+    const double n_flt = static_cast<double>(n);
+    const Vec2d centroid { sum / n_flt };
+
+    // Compute the normalized moments of the data set.
+    double Mxx = 0, Myy = 0, Mxy = 0, Mxz = 0, Myz = 0, Mzz = 0;
+    for (auto it = input_begin; it < input_end; ++it) {
+        // center/normalize the data.
+        double Xi {it->x() - centroid.x()};
+        double Yi {it->y() - centroid.y()};
+        double Zi {Xi*Xi + Yi*Yi};
+        Mxy += (Xi*Yi);
+        Mxx += (Xi*Xi);
+        Myy += (Yi*Yi);
+        Mxz += (Xi*Zi);
+        Myz += (Yi*Zi);
+        Mzz += (Zi*Zi);
+    }
+
+    // divide by number of points to get the moments
+    Mxx /= n_flt;
+    Myy /= n_flt;
+    Mxy /= n_flt;
+    Mxz /= n_flt;
+    Myz /= n_flt;
+    Mzz /= n_flt;
+
+    // Compute the coefficients of the characteristic polynomial for the circle
+    // eq 5.60
+    const double Mz {Mxx + Myy}; // xx + yy = z
+    const double Cov_xy {Mxx*Myy - Mxy*Mxy}; // this shows up a couple times so cache it here.
+    const double C3 {4.0*Mz};
+    const double C2 {-3.0*(Mz*Mz) - Mzz};
+    const double C1 {Mz*(Mzz - (Mz*Mz)) + 4.0*Mz*Cov_xy - (Mxz*Mxz) - (Myz*Myz)};
+    const double C0 {(Mxz*Mxz)*Myy + (Myz*Myz)*Mxx - 2.0*Mxz*Myz*Mxy - Cov_xy*(Mzz - (Mz*Mz))};
+
+    const double C22 = {C2 + C2};
+    const double C33 = {C3 + C3 + C3};
+
+    // solve the characteristic polynomial with Newton's method.
+    double xnew = 0.0;
+    double ynew = 1e20;
+
+    for (size_t i = 0; i < cycles; ++i) {
+        const double yold {ynew};
+        ynew = C0 + xnew * (C1 + xnew*(C2 + xnew * C3));
+        if (std::abs(ynew) > std::abs(yold)) {
+			BOOST_LOG_TRIVIAL(error) << "Geometry: Fit is going in the wrong direction.\n";
+            return Vec2d(std::nan(""), std::nan(""));
+        }
+        const double Dy {C1 + xnew*(C22 + xnew*C33)};
+
+        const double xold {xnew};
+        xnew = xold - (ynew / Dy);
+
+        if (std::abs((xnew-xold) / xnew) < 1e-12) i = cycles; // converged, we're done here
+
+        if (xnew < 0) {
+            // reset, we went negative
+            xnew = 0.0;
+        }
+    }
+    
+    // compute the determinant and the circle's parameters now that we've solved.
+    double DET = xnew*xnew - xnew*Mz + Cov_xy;
+
+    Vec2d center(Mxz * (Myy - xnew) - Myz * Mxy, Myz * (Mxx - xnew) - Mxz*Mxy);
+    center /= (DET * 2.);
+    return center + centroid;
+}
+
 void simplify_polygons(const Polygons &polygons, double tolerance, Polygons* retval)
 {
     Polygons pp;

src/libslic3r/Geometry.hpp

@@ -162,6 +162,15 @@ template<typename T> T angle_to_0_2PI(T angle)
 
     return angle;
 }
+
+/// Find the center of the circle corresponding to the vector of Points as an arc.
+Point circle_taubin_newton(const Points::const_iterator& input_start, const Points::const_iterator& input_end, size_t cycles = 20);
+inline Point circle_taubin_newton(const Points& input, size_t cycles = 20) { return circle_taubin_newton(input.cbegin(), input.cend(), cycles); }
+
+/// Find the center of the circle corresponding to the vector of Pointfs as an arc.
+Vec2d circle_taubin_newton(const Vec2ds::const_iterator& input_start, const Vec2ds::const_iterator& input_end, size_t cycles = 20);
+inline Vec2d circle_taubin_newton(const Vec2ds& input, size_t cycles = 20) { return circle_taubin_newton(input.cbegin(), input.cend(), cycles); }
+
 void simplify_polygons(const Polygons &polygons, double tolerance, Polygons* retval);
 
 double linint(double value, double oldmin, double oldmax, double newmin, double newmax);

src/libslic3r/Line.cpp

@@ -86,10 +86,7 @@ bool Line::intersection(const Line &l2, Point *intersection) const
     const Line  &l1  = *this;
     const Vec2d  v1  = (l1.b - l1.a).cast<double>();
     const Vec2d  v2  = (l2.b - l2.a).cast<double>();
-    const Vec2d  v12 = (l1.a - l2.a).cast<double>();
     double       denom  = cross2(v1, v2);
-    double       nume_a = cross2(v2, v12);
-    double       nume_b = cross2(v1, v12);
     if (fabs(denom) < EPSILON)
 #if 0
         // Lines are collinear. Return true if they are coincident (overlappign).
@@ -97,6 +94,9 @@ bool Line::intersection(const Line &l2, Point *intersection) const
 #else
         return false;
 #endif
+    const Vec2d v12 = (l1.a - l2.a).cast<double>();
+    double nume_a = cross2(v2, v12);
+    double nume_b = cross2(v1, v12);
     double t1 = nume_a / denom;
     double t2 = nume_b / denom;
     if (t1 >= 0 && t1 <= 1.0f && t2 >= 0 && t2 <= 1.0f) {

src/libslic3r/Point.hpp

@@ -38,6 +38,7 @@ typedef std::vector<Point*>                             PointPtrs;
 typedef std::vector<const Point*>                       PointConstPtrs;
 typedef std::vector<Vec3crd>                            Points3;
 typedef std::vector<Vec2d>                              Pointfs;
+typedef std::vector<Vec2d>                              Vec2ds;
 typedef std::vector<Vec3d>                              Pointf3s;
 
 typedef Eigen::Matrix<float,  2, 2, Eigen::DontAlign> Matrix2f;
@@ -87,12 +88,13 @@ class Point : public Vec2crd
 public:
     typedef coord_t coord_type;
 
-    Point() : Vec2crd() { (*this)(0) = 0; (*this)(1) = 0; }
-    Point(coord_t x, coord_t y) { (*this)(0) = x; (*this)(1) = y; }
-    Point(int64_t x, int64_t y) { (*this)(0) = coord_t(x); (*this)(1) = coord_t(y); } // for Clipper
-    Point(double x, double y) { (*this)(0) = coord_t(lrint(x)); (*this)(1) = coord_t(lrint(y)); }
+    Point() : Vec2crd(0, 0) {}
+    Point(coord_t x, coord_t y) : Vec2crd(x, y) {}
+    Point(int64_t x, int64_t y) : Vec2crd(coord_t(x), coord_t(y)) {} // for Clipper
+    Point(double x, double y) : Vec2crd(coord_t(lrint(x)), coord_t(lrint(y))) {}
     Point(const Point &rhs) { *this = rhs; }
-    // This constructor allows you to construct Point from Eigen expressions
+	explicit Point(const Vec2d& rhs) : Vec2crd(coord_t(lrint(rhs.x())), coord_t(lrint(rhs.y()))) {}
+	// This constructor allows you to construct Point from Eigen expressions
     template<typename OtherDerived>
     Point(const Eigen::MatrixBase<OtherDerived> &other) : Vec2crd(other) {}
     static Point new_scale(coordf_t x, coordf_t y) { return Point(coord_t(scale_(x)), coord_t(scale_(y))); }
@@ -126,6 +128,18 @@ public:
     Point  projection_onto(const Line &line) const;
 };
 
+inline bool is_approx(const Point &p1, const Point &p2, coord_t epsilon = coord_t(SCALED_EPSILON))
+{
+	Point d = (p2 - p1).cwiseAbs();
+	return d.x() < epsilon && d.y() < epsilon;
+}
+
+inline bool is_approx(const Vec2d &p1, const Vec2d &p2, double epsilon = EPSILON)
+{
+	Vec2d d = (p2 - p1).cwiseAbs();
+	return d.x() < epsilon && d.y() < epsilon;
+}
+
 namespace int128 {
     // Exact orientation predicate,
     // returns +1: CCW, 0: collinear, -1: CW.

src/libslic3r/Polygon.cpp

@@ -175,16 +175,16 @@ Point Polygon::centroid() const
 Points Polygon::concave_points(double angle) const
 {
     Points points;
-    angle = 2*PI - angle;
+    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);
-    }
+    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)
@@ -198,7 +198,7 @@ Points Polygon::concave_points(double angle) const
 Points Polygon::convex_points(double angle) const
 {
     Points points;
-    angle = 2*PI - angle;
+    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)