add test_geometry

add epsilon for concave_points / convecx_points

src/libslic3r/Geometry.cpp

@@ -1189,4 +1189,108 @@ Transformation Transformation::operator * (const Transformation& other) const
     return Transformation(get_matrix() * other.get_matrix());
 }
 
+//-------------- for tests ----------------------//
+
+Point
+circle_taubin_newton(const Points& input, size_t cycles)
+{
+    return circle_taubin_newton(input.cbegin(), input.cend(), cycles);
+}
+
+Point
+circle_taubin_newton(const Points::const_iterator& input_begin, const Points::const_iterator& input_end, size_t cycles)
+{
+    Pointfs 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); });
+    return Point::new_scale(circle_taubin_newton(tmp.cbegin(), tmp.end(), cycles));
+}
+
+Vec2d
+circle_taubin_newton(const Pointfs& input, size_t cycles)
+{
+    return circle_taubin_newton(input.cbegin(), input.cend(), cycles);
+}
+
+
+/// 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 Pointfs::const_iterator& input_begin, const Pointfs::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 = center / (DET / 2);
+
+    return center + centroid;
+}
+
 } }

src/libslic3r/Geometry.hpp

@@ -241,6 +241,17 @@ public:
     Transformation operator * (const Transformation& other) const;
 };
 
+//------------for tests------------//
+ 
+/// Find the center of the circle corresponding to the vector of Points as an arc.
+Point circle_taubin_newton(const Points& input, size_t cycles = 20);
+Point circle_taubin_newton(const Points::const_iterator& input_start, const Points::const_iterator& input_end, size_t cycles = 20);
+
+/// Find the center of the circle corresponding to the vector of Pointfs as an arc.
+Vec2d circle_taubin_newton(const Pointfs& input, size_t cycles = 20);
+Vec2d circle_taubin_newton(const Pointfs::const_iterator& input_start, const Pointfs::const_iterator& input_end, size_t cycles = 20);
+
+
 } }
 
 #endif

src/libslic3r/Polygon.cpp

@@ -236,7 +236,7 @@ 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)
@@ -260,7 +260,7 @@ 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)

src/libslic3r/libslic3r.h

@@ -87,6 +87,7 @@ inline T unscale(Q v) { return T(v) * T(SCALING_FACTOR); }
 
 inline double unscaled(double v) { return v * SCALING_FACTOR; }
 inline coordf_t unscale_(coord_t v) { return v * SCALING_FACTOR; }
+inline coord_t scale_t(coordf_t v) { return (coord_t)(v / SCALING_FACTOR); }
 
 enum Axis { X=0, Y, Z, E, F, NUM_AXES };
 

src/test/CMakeLists.txt

@@ -14,6 +14,7 @@ set(SLIC3R_TEST_SOURCES
     libslic3r/test_fill.cpp
     libslic3r/test_flow.cpp
     libslic3r/test_gcodewriter.cpp
+    libslic3r/test_geometry.cpp
 )
 
 if (NOT TARGET Catch)

src/test/libslic3r/test_geometry.cpp

@@ -1,19 +1,19 @@
 
 #include <catch.hpp>
 
-#include "Point.hpp"
-#include "BoundingBox.hpp"
-#include "Polygon.hpp"
-#include "Polyline.hpp"
-#include "Line.hpp"
-#include "Geometry.hpp"
-#include "ClipperUtils.hpp"
+#include "../../libslic3r/Point.hpp"
+#include "../../libslic3r/BoundingBox.hpp"
+#include "../../libslic3r/Polygon.hpp"
+#include "../../libslic3r/Polyline.hpp"
+#include "../../libslic3r/Line.hpp"
+#include "../../libslic3r/Geometry.hpp"
+#include "../../libslic3r/ClipperUtils.hpp"
 
 using namespace Slic3r;
 
 TEST_CASE("Polygon::contains works properly", ""){
    // this test was failing on Windows (GH #1950)
-    auto polygon = Polygon(std::vector<Point>({
+    auto polygon = Polygon(Points({
         Point(207802834,-57084522),
         Point(196528149,-37556190),
         Point(173626821,-25420928),
@@ -155,34 +155,33 @@ TEST_CASE("Bounding boxes are scaled appropriately"){
 
 
 TEST_CASE("Offseting a line generates a polygon correctly"){
-    auto line = Line(Point(10,10), Point(20,10));
-    Polyline tmp(line);
+    Polyline tmp({ Point(10,10), Point(20,10) });
     Polygon area = offset(tmp,5).at(0);
     REQUIRE(area.area() == Polygon(std::vector<Point>({Point(10,5),Point(20,5),Point(20,15),Point(10,15)})).area());
 }
 
 SCENARIO("Circle Fit, TaubinFit with Newton's method") {
     GIVEN("A vector of Pointfs arranged in a half-circle with approximately the same distance R from some point") {
-        Pointf expected_center(-6, 0);
-        Pointfs sample {Pointf(6.0, 0), Pointf(5.1961524, 3), Pointf(3 ,5.1961524), Pointf(0, 6.0), Pointf(3, 5.1961524), Pointf(-5.1961524, 3), Pointf(-6.0, 0)};
-        std::transform(sample.begin(), sample.end(), sample.begin(), [expected_center] (const Pointf& a) { return a + expected_center;});
+        Vec2d expected_center(-6, 0);
+        Pointfs sample {Vec2d(6.0, 0), Vec2d(5.1961524, 3), Vec2d(3 ,5.1961524), Vec2d(0, 6.0), Vec2d(-3, 5.1961524), Vec2d(-5.1961524, 3), Vec2d(-6.0, 0)};
+        std::transform(sample.begin(), sample.end(), sample.begin(), [expected_center] (const Vec2d& a) { return a + expected_center;});
 
         WHEN("Circle fit is called on the entire array") {
-            Pointf result_center(0,0);
+            Vec2d result_center(0,0);
             result_center = Geometry::circle_taubin_newton(sample);
             THEN("A center point of -6,0 is returned.") {
                 REQUIRE(result_center == expected_center);
             }
         }
         WHEN("Circle fit is called on the first four points") {
-            Pointf result_center(0,0);
+            Vec2d result_center(0,0);
             result_center = Geometry::circle_taubin_newton(sample.cbegin(), sample.cbegin()+4);
             THEN("A center point of -6,0 is returned.") {
                 REQUIRE(result_center == expected_center);
             }
         }
         WHEN("Circle fit is called on the middle four points") {
-            Pointf result_center(0,0);
+            Vec2d result_center(0,0);
             result_center = Geometry::circle_taubin_newton(sample.cbegin()+2, sample.cbegin()+6);
             THEN("A center point of -6,0 is returned.") {
                 REQUIRE(result_center == expected_center);
@@ -190,30 +189,30 @@ SCENARIO("Circle Fit, TaubinFit with Newton's method") {
         }
     }
     GIVEN("A vector of Pointfs arranged in a half-circle with approximately the same distance R from some point") {
-        Pointf expected_center(-3, 9);
-        Pointfs sample {Pointf(6.0, 0), Pointf(5.1961524, 3), Pointf(3 ,5.1961524), 
-                        Pointf(0, 6.0), 
-                        Pointf(3, 5.1961524), Pointf(-5.1961524, 3), Pointf(-6.0, 0)};
+        Vec2d expected_center(-3, 9);
+        Pointfs sample {Vec2d(6.0, 0), Vec2d(5.1961524, 3), Vec2d(3 ,5.1961524), 
+                        Vec2d(0, 6.0), 
+                        Vec2d(3, 5.1961524), Vec2d(-5.1961524, 3), Vec2d(-6.0, 0)};
 
-        std::transform(sample.begin(), sample.end(), sample.begin(), [expected_center] (const Pointf& a) { return a + expected_center;});
+        std::transform(sample.begin(), sample.end(), sample.begin(), [expected_center] (const Vec2d& a) { return a + expected_center;});
 
 
         WHEN("Circle fit is called on the entire array") {
-            Pointf result_center(0,0);
+            Vec2d result_center(0,0);
             result_center = Geometry::circle_taubin_newton(sample);
             THEN("A center point of 3,9 is returned.") {
                 REQUIRE(result_center == expected_center);
             }
         }
         WHEN("Circle fit is called on the first four points") {
-            Pointf result_center(0,0);
+            Vec2d result_center(0,0);
             result_center = Geometry::circle_taubin_newton(sample.cbegin(), sample.cbegin()+4);
             THEN("A center point of 3,9 is returned.") {
                 REQUIRE(result_center == expected_center);
             }
         }
         WHEN("Circle fit is called on the middle four points") {
-            Pointf result_center(0,0);
+            Vec2d result_center(0,0);
             result_center = Geometry::circle_taubin_newton(sample.cbegin()+2, sample.cbegin()+6);
             THEN("A center point of 3,9 is returned.") {
                 REQUIRE(result_center == expected_center);
@@ -262,7 +261,7 @@ TEST_CASE("Chained path working correctly"){
     Geometry::chained_path(points,indices);
     for(Points::size_type i = 0; i < indices.size()-1;i++){
         double dist = points.at(indices.at(i)).distance_to(points.at(indices.at(i+1)));
-        REQUIRE(abs(dist-26) <= Geometry::epsilon);
+        REQUIRE(abs(dist-26) <= EPSILON);
     }
 }