diff --git a/src/libslic3r/Geometry.cpp b/src/libslic3r/Geometry.cpp index 4f833f5f3..21abb8bd8 100644 --- a/src/libslic3r/Geometry.cpp +++ b/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(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; +} + } } diff --git a/src/libslic3r/Geometry.hpp b/src/libslic3r/Geometry.hpp index e08bcf063..72489d618 100644 --- a/src/libslic3r/Geometry.hpp +++ b/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 diff --git a/src/libslic3r/Polygon.cpp b/src/libslic3r/Polygon.cpp index c33983dfc..2fe3c4d18 100644 --- a/src/libslic3r/Polygon.cpp +++ b/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) diff --git a/src/libslic3r/libslic3r.h b/src/libslic3r/libslic3r.h index 1dbffbd26..08383b6bf 100644 --- a/src/libslic3r/libslic3r.h +++ b/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 }; diff --git a/src/test/CMakeLists.txt b/src/test/CMakeLists.txt index 888632a6d..a97c2d2f3 100644 --- a/src/test/CMakeLists.txt +++ b/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) diff --git a/src/test/libslic3r/test_geometry.cpp b/src/test/libslic3r/test_geometry.cpp index a65e85313..2c5bad3f8 100644 --- a/src/test/libslic3r/test_geometry.cpp +++ b/src/test/libslic3r/test_geometry.cpp @@ -1,19 +1,19 @@ #include -#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({ + 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(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); } }