From 42a858b9999384272457f175849e15f80b65882c Mon Sep 17 00:00:00 2001 From: bubnikv Date: Tue, 15 Oct 2019 09:40:40 +0200 Subject: [PATCH] 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 | 88 +++++++ src/libslic3r/Geometry.hpp | 9 + src/libslic3r/Line.cpp | 6 +- src/libslic3r/Point.hpp | 24 +- src/libslic3r/Polygon.cpp | 10 +- tests/CMakeLists.txt | 3 +- tests/fff_print/CMakeLists.txt | 7 + tests/fff_print/fff_print_tests.cpp | 15 ++ tests/libslic3r/CMakeLists.txt | 10 + tests/libslic3r/libslic3r_tests.cpp | 15 ++ tests/libslic3r/test_geometry.cpp | 375 ++++++++++++++++++++++++++++ 11 files changed, 548 insertions(+), 14 deletions(-) create mode 100644 tests/fff_print/CMakeLists.txt create mode 100644 tests/fff_print/fff_print_tests.cpp create mode 100644 tests/libslic3r/CMakeLists.txt create mode 100644 tests/libslic3r/libslic3r_tests.cpp create mode 100644 tests/libslic3r/test_geometry.cpp diff --git a/src/libslic3r/Geometry.cpp b/src/libslic3r/Geometry.cpp index 3adf8c670..e5dfacc09 100644 --- a/src/libslic3r/Geometry.cpp +++ b/src/libslic3r/Geometry.cpp @@ -16,6 +16,7 @@ #include #include +#include #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(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; diff --git a/src/libslic3r/Geometry.hpp b/src/libslic3r/Geometry.hpp index 32b66663e..44303711b 100644 --- a/src/libslic3r/Geometry.hpp +++ b/src/libslic3r/Geometry.hpp @@ -162,6 +162,15 @@ template 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); diff --git a/src/libslic3r/Line.cpp b/src/libslic3r/Line.cpp index baa04795a..e5f7b8fa9 100644 --- a/src/libslic3r/Line.cpp +++ b/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(); const Vec2d v2 = (l2.b - l2.a).cast(); - const Vec2d v12 = (l1.a - l2.a).cast(); 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 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) { diff --git a/src/libslic3r/Point.hpp b/src/libslic3r/Point.hpp index 994f45e59..5f5d86aa8 100644 --- a/src/libslic3r/Point.hpp +++ b/src/libslic3r/Point.hpp @@ -38,6 +38,7 @@ typedef std::vector PointPtrs; typedef std::vector PointConstPtrs; typedef std::vector Points3; typedef std::vector Pointfs; +typedef std::vector Vec2ds; typedef std::vector Pointf3s; typedef Eigen::Matrix 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 Point(const Eigen::MatrixBase &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. diff --git a/src/libslic3r/Polygon.cpp b/src/libslic3r/Polygon.cpp index 5ef5d0ceb..4c2db2848 100644 --- a/src/libslic3r/Polygon.cpp +++ b/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) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6cbfd5432..e957c0c20 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -21,6 +21,7 @@ endif() set_property(GLOBAL PROPERTY USE_FOLDERS ON) add_subdirectory(libnest2d) +add_subdirectory(libslic3r) add_subdirectory(timeutils) +add_subdirectory(fff_print) add_subdirectory(sla_print) - diff --git a/tests/fff_print/CMakeLists.txt b/tests/fff_print/CMakeLists.txt new file mode 100644 index 000000000..d0b51a01d --- /dev/null +++ b/tests/fff_print/CMakeLists.txt @@ -0,0 +1,7 @@ +get_filename_component(_TEST_NAME ${CMAKE_CURRENT_LIST_DIR} NAME) +add_executable(${_TEST_NAME}_tests ${_TEST_NAME}_tests.cpp) +target_link_libraries(${_TEST_NAME}_tests test_common libslic3r) +set_property(TARGET ${_TEST_NAME}_tests PROPERTY FOLDER "tests") + +# catch_discover_tests(${_TEST_NAME}_tests TEST_PREFIX "${_TEST_NAME}: ") +add_test(${_TEST_NAME}_tests ${_TEST_NAME}_tests "--durations yes") diff --git a/tests/fff_print/fff_print_tests.cpp b/tests/fff_print/fff_print_tests.cpp new file mode 100644 index 000000000..907304f57 --- /dev/null +++ b/tests/fff_print/fff_print_tests.cpp @@ -0,0 +1,15 @@ +#define CATCH_CONFIG_MAIN +#include + +#include "libslic3r/libslic3r.h" + +namespace { + +TEST_CASE("sort_remove_duplicates", "[utils]") { + std::vector data_src = { 3, 0, 2, 1, 15, 3, 5, 6, 3, 1, 0 }; + std::vector data_dst = { 0, 1, 2, 3, 5, 6, 15 }; + Slic3r::sort_remove_duplicates(data_src); + REQUIRE(data_src == data_dst); +} + +} diff --git a/tests/libslic3r/CMakeLists.txt b/tests/libslic3r/CMakeLists.txt new file mode 100644 index 000000000..ea83bf089 --- /dev/null +++ b/tests/libslic3r/CMakeLists.txt @@ -0,0 +1,10 @@ +get_filename_component(_TEST_NAME ${CMAKE_CURRENT_LIST_DIR} NAME) +add_executable(${_TEST_NAME}_tests + ${_TEST_NAME}_tests.cpp + test_geometry.cpp + ) +target_link_libraries(${_TEST_NAME}_tests test_common libslic3r) +set_property(TARGET ${_TEST_NAME}_tests PROPERTY FOLDER "tests") + +# catch_discover_tests(${_TEST_NAME}_tests TEST_PREFIX "${_TEST_NAME}: ") +add_test(${_TEST_NAME}_tests ${_TEST_NAME}_tests "--durations yes") diff --git a/tests/libslic3r/libslic3r_tests.cpp b/tests/libslic3r/libslic3r_tests.cpp new file mode 100644 index 000000000..907304f57 --- /dev/null +++ b/tests/libslic3r/libslic3r_tests.cpp @@ -0,0 +1,15 @@ +#define CATCH_CONFIG_MAIN +#include + +#include "libslic3r/libslic3r.h" + +namespace { + +TEST_CASE("sort_remove_duplicates", "[utils]") { + std::vector data_src = { 3, 0, 2, 1, 15, 3, 5, 6, 3, 1, 0 }; + std::vector data_dst = { 0, 1, 2, 3, 5, 6, 15 }; + Slic3r::sort_remove_duplicates(data_src); + REQUIRE(data_src == data_dst); +} + +} diff --git a/tests/libslic3r/test_geometry.cpp b/tests/libslic3r/test_geometry.cpp new file mode 100644 index 000000000..7a3b84c7d --- /dev/null +++ b/tests/libslic3r/test_geometry.cpp @@ -0,0 +1,375 @@ +#include + +#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" +#include "libslic3r/ShortestPath.hpp" + +using namespace Slic3r; + +TEST_CASE("Polygon::contains works properly", ""){ + // this test was failing on Windows (GH #1950) + auto polygon = Slic3r::Polygon(std::vector({ + Point(207802834,-57084522), + Point(196528149,-37556190), + Point(173626821,-25420928), + Point(171285751,-21366123), + Point(118673592,-21366123), + Point(116332562,-25420928), + Point(93431208,-37556191), + Point(82156517,-57084523), + Point(129714478,-84542120), + Point(160244873,-84542120) + })); + auto point = Point(95706562, -57294774); + REQUIRE(polygon.contains(point)); +} + +SCENARIO("Intersections of line segments"){ + GIVEN("Integer coordinates"){ + auto line1 = Line(Point(5,15),Point(30,15)); + auto line2 = Line(Point(10,20), Point(10,10)); + THEN("The intersection is valid"){ + Point point; + line1.intersection(line2,&point); + REQUIRE(Point(10,15) == point); + } + } + + GIVEN("Scaled coordinates"){ + auto line1 = Line(Point(73.6310778185108 / 0.00001, 371.74239268924 / 0.00001), Point(73.6310778185108 / 0.00001, 501.74239268924 / 0.00001)); + auto line2 = Line(Point(75/0.00001, 437.9853/0.00001), Point(62.7484/0.00001, 440.4223/0.00001)); + THEN("There is still an intersection"){ + Point point; + REQUIRE(line1.intersection(line2,&point)); + } + } +} + +/* +Tests for unused methods still written in perl +{ + my $polygon = Slic3r::Polygon->new( + [45919000, 515273900], [14726100, 461246400], [14726100, 348753500], [33988700, 315389800], + [43749700, 343843000], [45422300, 352251500], [52362100, 362637800], [62748400, 369577600], + [75000000, 372014700], [87251500, 369577600], [97637800, 362637800], [104577600, 352251500], + [107014700, 340000000], [104577600, 327748400], [97637800, 317362100], [87251500, 310422300], + [82789200, 309534700], [69846100, 294726100], [254081000, 294726100], [285273900, 348753500], + [285273900, 461246400], [254081000, 515273900], + ); + + # this points belongs to $polyline + # note: it's actually a vertex, while we should better check an intermediate point + my $point = Slic3r::Point->new(104577600, 327748400); + + local $Slic3r::Geometry::epsilon = 1E-5; + is_deeply Slic3r::Geometry::polygon_segment_having_point($polygon, $point)->pp, + [ [107014700, 340000000], [104577600, 327748400] ], + 'polygon_segment_having_point'; +} +{ + auto point = Point(736310778.185108, 5017423926.8924); + auto line = Line(Point((long int) 627484000, (long int) 3695776000), Point((long int) 750000000, (long int)3720147000)); + //is Slic3r::Geometry::point_in_segment($point, $line), 0, 'point_in_segment'; +} + +// Possible to delete +{ + //my $p1 = [10, 10]; + //my $p2 = [10, 20]; + //my $p3 = [10, 30]; + //my $p4 = [20, 20]; + //my $p5 = [0, 20]; + + THEN("Points in a line give the correct angles"){ + //is Slic3r::Geometry::angle3points($p2, $p3, $p1), PI(), 'angle3points'; + //is Slic3r::Geometry::angle3points($p2, $p1, $p3), PI(), 'angle3points'; + } + THEN("Left turns give the correct angle"){ + //is Slic3r::Geometry::angle3points($p2, $p4, $p3), PI()/2, 'angle3points'; + //is Slic3r::Geometry::angle3points($p2, $p1, $p4), PI()/2, 'angle3points'; + } + THEN("Right turns give the correct angle"){ + //is Slic3r::Geometry::angle3points($p2, $p3, $p4), PI()/2*3, 'angle3points'; + //is Slic3r::Geometry::angle3points($p2, $p1, $p5), PI()/2*3, 'angle3points'; + } + //my $p1 = [30, 30]; + //my $p2 = [20, 20]; + //my $p3 = [10, 10]; + //my $p4 = [30, 10]; + + //is Slic3r::Geometry::angle3points($p2, $p1, $p3), PI(), 'angle3points'; + //is Slic3r::Geometry::angle3points($p2, $p1, $p4), PI()/2*3, 'angle3points'; + //is Slic3r::Geometry::angle3points($p2, $p1, $p1), 2*PI(), 'angle3points'; +} + +SCENARIO("polygon_is_convex works"){ + GIVEN("A square of dimension 10"){ + //my $cw_square = [ [0,0], [0,10], [10,10], [10,0] ]; + THEN("It is not convex clockwise"){ + //is polygon_is_convex($cw_square), 0, 'cw square is not convex'; + } + THEN("It is convex counter-clockwise"){ + //is polygon_is_convex([ reverse @$cw_square ]), 1, 'ccw square is convex'; + } + + } + GIVEN("A concave polygon"){ + //my $convex1 = [ [0,0], [10,0], [10,10], [0,10], [0,6], [4,6], [4,4], [0,4] ]; + THEN("It is concave"){ + //is polygon_is_convex($convex1), 0, 'concave polygon'; + } + } +}*/ + + +TEST_CASE("Creating a polyline generates the obvious lines"){ + auto polyline = Slic3r::Polyline(); + polyline.points = std::vector({Point(0, 0), Point(10, 0), Point(20, 0)}); + REQUIRE(polyline.lines().at(0).a == Point(0,0)); + REQUIRE(polyline.lines().at(0).b == Point(10,0)); + REQUIRE(polyline.lines().at(1).a == Point(10,0)); + REQUIRE(polyline.lines().at(1).b == Point(20,0)); +} + +TEST_CASE("Splitting a Polygon generates a polyline correctly"){ + auto polygon = Slic3r::Polygon(std::vector({Point(0, 0), Point(10, 0), Point(5, 5)})); + auto split = polygon.split_at_index(1); + REQUIRE(split.points[0]==Point(10,0)); + REQUIRE(split.points[1]==Point(5,5)); + REQUIRE(split.points[2]==Point(0,0)); + REQUIRE(split.points[3]==Point(10,0)); +} + + +TEST_CASE("Bounding boxes are scaled appropriately"){ + auto bb = BoundingBox(std::vector({Point(0, 1), Point(10, 2), Point(20, 2)})); + bb.scale(2); + REQUIRE(bb.min == Point(0,2)); + REQUIRE(bb.max == Point(40,4)); +} + + +TEST_CASE("Offseting a line generates a polygon correctly"){ + Slic3r::Polyline tmp = { Point(10,10), Point(20,10) }; + Slic3r::Polygon area = offset(tmp,5).at(0); + REQUIRE(area.area() == Slic3r::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 Vec2ds arranged in a half-circle with approximately the same distance R from some point") { + Vec2d expected_center(-6, 0); + Vec2ds 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") { + Vec2d result_center(0,0); + result_center = Geometry::circle_taubin_newton(sample); + THEN("A center point of -6,0 is returned.") { + REQUIRE(is_approx(result_center, expected_center)); + } + } + WHEN("Circle fit is called on the first four points") { + 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(is_approx(result_center, expected_center)); + } + } + WHEN("Circle fit is called on the middle four points") { + 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(is_approx(result_center, expected_center)); + } + } + } + GIVEN("A vector of Vec2ds arranged in a half-circle with approximately the same distance R from some point") { + Vec2d expected_center(-3, 9); + Vec2ds 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") { + Vec2d result_center(0,0); + result_center = Geometry::circle_taubin_newton(sample); + THEN("A center point of 3,9 is returned.") { + REQUIRE(is_approx(result_center, expected_center)); + } + } + WHEN("Circle fit is called on the first four points") { + 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(is_approx(result_center, expected_center)); + } + } + WHEN("Circle fit is called on the middle four points") { + 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(is_approx(result_center, expected_center)); + } + } + } + GIVEN("A vector of Points arranged in a half-circle with approximately the same distance R from some point") { + Point expected_center { Point::new_scale(-3, 9)}; + Points sample {Point::new_scale(6.0, 0), Point::new_scale(5.1961524, 3), Point::new_scale(3 ,5.1961524), + Point::new_scale(0, 6.0), + Point::new_scale(3, 5.1961524), Point::new_scale(-5.1961524, 3), Point::new_scale(-6.0, 0)}; + + std::transform(sample.begin(), sample.end(), sample.begin(), [expected_center] (const Point& a) { return a + expected_center;}); + + + WHEN("Circle fit is called on the entire array") { + Point result_center(0,0); + result_center = Geometry::circle_taubin_newton(sample); + THEN("A center point of scaled 3,9 is returned.") { + REQUIRE(is_approx(result_center, expected_center)); + } + } + WHEN("Circle fit is called on the first four points") { + Point result_center(0,0); + result_center = Geometry::circle_taubin_newton(sample.cbegin(), sample.cbegin()+4); + THEN("A center point of scaled 3,9 is returned.") { + REQUIRE(is_approx(result_center, expected_center)); + } + } + WHEN("Circle fit is called on the middle four points") { + Point result_center(0,0); + result_center = Geometry::circle_taubin_newton(sample.cbegin()+2, sample.cbegin()+6); + THEN("A center point of scaled 3,9 is returned.") { + REQUIRE(is_approx(result_center, expected_center)); + } + } + } +} + +TEST_CASE("Chained path working correctly"){ + // if chained_path() works correctly, these points should be joined with no diagonal paths + // (thus 26 units long) + std::vector points = {Point(26,26),Point(52,26),Point(0,26),Point(26,52),Point(26,0),Point(0,52),Point(52,52),Point(52,0)}; + std::vector indices = chain_points(points); + for (Points::size_type i = 0; i + 1 < indices.size(); ++ i) { + double dist = (points.at(indices.at(i)).cast() - points.at(indices.at(i+1)).cast()).norm(); + REQUIRE(std::abs(dist-26) <= EPSILON); + } +} + +SCENARIO("Line distances"){ + GIVEN("A line"){ + auto line = Line(Point(0, 0), Point(20, 0)); + THEN("Points on the line segment have 0 distance"){ + REQUIRE(line.distance_to(Point(0, 0)) == 0); + REQUIRE(line.distance_to(Point(20, 0)) == 0); + REQUIRE(line.distance_to(Point(10, 0)) == 0); + + } + THEN("Points off the line have the appropriate distance"){ + REQUIRE(line.distance_to(Point(10, 10)) == 10); + REQUIRE(line.distance_to(Point(50, 0)) == 30); + } + } +} + +SCENARIO("Polygon convex/concave detection"){ + GIVEN(("A Square with dimension 100")){ + auto square = Slic3r::Polygon /*new_scale*/(std::vector({ + Point(100,100), + Point(200,100), + Point(200,200), + Point(100,200)})); + THEN("It has 4 convex points counterclockwise"){ + REQUIRE(square.concave_points(PI*4/3).size() == 0); + REQUIRE(square.convex_points(PI*2/3).size() == 4); + } + THEN("It has 4 concave points clockwise"){ + square.make_clockwise(); + REQUIRE(square.concave_points(PI*4/3).size() == 4); + REQUIRE(square.convex_points(PI*2/3).size() == 0); + } + } + GIVEN("A Square with an extra colinearvertex"){ + auto square = Slic3r::Polygon /*new_scale*/(std::vector({ + Point(150,100), + Point(200,100), + Point(200,200), + Point(100,200), + Point(100,100)})); + THEN("It has 4 convex points counterclockwise"){ + REQUIRE(square.concave_points(PI*4/3).size() == 0); + REQUIRE(square.convex_points(PI*2/3).size() == 4); + } + } + GIVEN("A Square with an extra collinear vertex in different order"){ + auto square = Slic3r::Polygon /*new_scale*/(std::vector({ + Point(200,200), + Point(100,200), + Point(100,100), + Point(150,100), + Point(200,100)})); + THEN("It has 4 convex points counterclockwise"){ + REQUIRE(square.concave_points(PI*4/3).size() == 0); + REQUIRE(square.convex_points(PI*2/3).size() == 4); + } + } + + GIVEN("A triangle"){ + auto triangle = Slic3r::Polygon(std::vector({ + Point(16000170,26257364), + Point(714223,461012), + Point(31286371,461008) + })); + THEN("it has three convex vertices"){ + REQUIRE(triangle.concave_points(PI*4/3).size() == 0); + REQUIRE(triangle.convex_points(PI*2/3).size() == 3); + } + } + + GIVEN("A triangle with an extra collinear point"){ + auto triangle = Slic3r::Polygon(std::vector({ + Point(16000170,26257364), + Point(714223,461012), + Point(20000000,461012), + Point(31286371,461012) + })); + THEN("it has three convex vertices"){ + REQUIRE(triangle.concave_points(PI*4/3).size() == 0); + REQUIRE(triangle.convex_points(PI*2/3).size() == 3); + } + } + GIVEN("A polygon with concave vertices with angles of specifically 4/3pi"){ + // Two concave vertices of this polygon have angle = PI*4/3, so this test fails + // if epsilon is not used. + auto polygon = Slic3r::Polygon(std::vector({ + Point(60246458,14802768),Point(64477191,12360001), + Point(63727343,11060995),Point(64086449,10853608), + Point(66393722,14850069),Point(66034704,15057334), + Point(65284646,13758387),Point(61053864,16200839), + Point(69200258,30310849),Point(62172547,42483120), + Point(61137680,41850279),Point(67799985,30310848), + Point(51399866,1905506),Point(38092663,1905506), + Point(38092663,692699),Point(52100125,692699) + })); + THEN("the correct number of points are detected"){ + REQUIRE(polygon.concave_points(PI*4/3).size() == 6); + REQUIRE(polygon.convex_points(PI*2/3).size() == 10); + } + } +} + +TEST_CASE("Triangle Simplification does not result in less than 3 points"){ + auto triangle = Slic3r::Polygon(std::vector({ + Point(16000170,26257364), Point(714223,461012), Point(31286371,461008) + })); + REQUIRE(triangle.simplify(250000).at(0).points.size() == 3); +} + +