From 8fcccce0318e0ecc8761eade7ec3ee73fabbdf21 Mon Sep 17 00:00:00 2001 From: stribor14 Date: Tue, 29 Mar 2022 15:00:31 +0200 Subject: [PATCH] Raise C++ version to c++17 - update STL function to c++17 version - fix root solving for axis parallel linear curves --- CMakeLists.txt | 4 +- README.md | 4 +- example/bezier_example.pro | 4 +- example/customscene.cpp | 4 +- example/qpolycurve.h | 5 ++- include/Bezier/bezier.h | 2 +- include/Bezier/polycurve.h | 19 +++++++- src/bezier.cpp | 91 ++++++++++++++++++++------------------ src/polycurve.cpp | 35 ++++++++------- 9 files changed, 99 insertions(+), 69 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2629899..4f83dfc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ project(Bezier LANGUAGES CXX VERSION 0.2.0) -set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) @@ -41,6 +41,8 @@ else() add_library(bezier STATIC ${Bezier_SRC}) endif() +target_link_libraries(bezier PUBLIC tbb) + target_include_directories(bezier PUBLIC $ $ diff --git a/README.md b/README.md index cc29c2e..8103c26 100644 --- a/README.md +++ b/README.md @@ -39,8 +39,8 @@ target_link_libraries(target bezier) - [ ] More sophisticated example ## Dependencies - - c++11 - - Eigen3 + - c++17 + - Eigen3.3 ## Instalation ### System-wide installation diff --git a/example/bezier_example.pro b/example/bezier_example.pro index 2133222..b456aa9 100644 --- a/example/bezier_example.pro +++ b/example/bezier_example.pro @@ -22,7 +22,7 @@ DEFINES += QT_DEPRECATED_WARNINGS # You can also select to disable deprecated APIs only up to a certain version of Qt. #DEFINES += QT_DISABLE_DEPRECATED_BEFORE=0x060000 # disables all the APIs deprecated before Qt 6.0.0 -CONFIG += c++11 +CONFIG += c++17 INCLUDEPATH += /usr/include/eigen3 \ ../include @@ -51,6 +51,8 @@ HEADERS += \ FORMS += \ mainwindow.ui +LIBS+= -ltbb + # Default rules for deployment. qnx: target.path = /tmp/$${TARGET}/bin else: unix:!android: target.path = /opt/$${TARGET}/bin diff --git a/example/customscene.cpp b/example/customscene.cpp index ba323a7..febfbb5 100644 --- a/example/customscene.cpp +++ b/example/customscene.cpp @@ -226,12 +226,12 @@ void CustomScene::mouseMoveEvent(QGraphicsSceneMouseEvent* mouseEvent) if (is_curve) { c_curve->prepareGeometryChange(); - c_curve->manipulateControlPoint(cp_to_update.second, p); + c_curve->moveControlPoint(cp_to_update.second, p); } if (is_poly) { c_poly->prepareGeometryChange(); - c_poly->manipulateControlPoint(cp_to_update.second, p); + c_poly->moveControlPoint(cp_to_update.second, p); } update(); } diff --git a/example/qpolycurve.h b/example/qpolycurve.h index b946cee..0b0529c 100644 --- a/example/qpolycurve.h +++ b/example/qpolycurve.h @@ -5,6 +5,7 @@ #include "Bezier/declarations.h" #include "Bezier/polycurve.h" +#include "Bezier/bezier.h" class qPolyCurve : public QGraphicsItem, public Bezier::PolyCurve { @@ -12,8 +13,8 @@ class qPolyCurve : public QGraphicsItem, public Bezier::PolyCurve bool draw_control_points = false; bool draw_curvature_radious = false; public: - qPolyCurve(const std::vector& curve_list) : QGraphicsItem(), Bezier::PolyCurve(curve_list) {} - qPolyCurve(const Bezier::Curve& curve) : QGraphicsItem(), Bezier::PolyCurve(curve) {} + qPolyCurve(const std::deque& curve_list) : QGraphicsItem(), Bezier::PolyCurve(curve_list) {} + qPolyCurve(const Bezier::Curve& curve) : QGraphicsItem(), Bezier::PolyCurve(std::deque{curve}) {} int type() const Q_DECL_OVERRIDE; void paint(QPainter* painter, const QStyleOptionGraphicsItem* option, QWidget* widget) Q_DECL_OVERRIDE; QRectF boundingRect() const Q_DECL_OVERRIDE; diff --git a/include/Bezier/bezier.h b/include/Bezier/bezier.h index 84ed6f3..1eef992 100644 --- a/include/Bezier/bezier.h +++ b/include/Bezier/bezier.h @@ -129,7 +129,7 @@ class Curve * \param idx Index of chosen control point * \param point New control point */ - void manipulateControlPoint(uint idx, const Point& point); + void moveControlPoint(uint idx, const Point& point); /*! * \brief Manipulate the curve so that it passes through wanted point for given 't' diff --git a/include/Bezier/polycurve.h b/include/Bezier/polycurve.h index 856ad56..cbb70df 100644 --- a/include/Bezier/polycurve.h +++ b/include/Bezier/polycurve.h @@ -19,7 +19,8 @@ #include -#include "declarations.h" +#include "Bezier/declarations.h" +#include "Bezier/bezier.h" namespace Bezier { @@ -175,7 +176,7 @@ class PolyCurve * \param index Index of chosen control point * \param point New control point */ - void manipulateControlPoint(uint idx, const Point& point); + void moveControlPoint(uint idx, const Point& point); /*! * \brief Get the point on polycurve for a given t @@ -265,6 +266,20 @@ class PolyCurve */ std::vector projectPoint(const PointVector& point_vector) const; + /*! + * \brief Get distance of the point to the polycurve + * \param point Point to project on the polycurve + * \return Distance to the curve + */ + double distance(const Point& point) const; + + /*! + * \brief Get the distance vector of points to the polycurve + * \param point_vector Points to project on the polycurve + * \return Vector of distances + */ + std::vector distance(const PointVector& point_vector) const; + protected: /// Structure for holding underlying Bezier curves std::deque curves_; diff --git a/src/bezier.cpp b/src/bezier.cpp index 6b587d9..6d86d5c 100644 --- a/src/bezier.cpp +++ b/src/bezier.cpp @@ -1,6 +1,8 @@ #include "Bezier/bezier.h" #include "Bezier/legendre_gauss.h" +#include +#include #include #include @@ -13,7 +15,7 @@ inline double binomial(uint n, uint k) { return factorial(n) / (factorial(k) * f inline Eigen::VectorXd trimZeroes(const Eigen::VectorXd& vec) { auto idx = vec.size(); - while (vec(idx - 1) == 0.0) + while (idx && vec(idx - 1) == 0.0) --idx; return vec.head(idx); } @@ -48,10 +50,7 @@ PointVector Curve::controlPoints() const return points; } -Point Curve::controlPoint(uint idx) const -{ - return control_points_.row(idx); -} +Point Curve::controlPoint(uint idx) const { return control_points_.row(idx); } std::pair Curve::endPoints() const { return {control_points_.row(0), control_points_.row(N_ - 1)}; } @@ -84,8 +83,8 @@ PointVector Curve::polyline(double smoothness, double precision) const } } - (const_cast(this))->cached_polyline_params_ = {smoothness, precision}; - (const_cast(this))->cached_polyline_.reset(polyline); + const_cast(this)->cached_polyline_params_ = {smoothness, precision}; + const_cast(this)->cached_polyline_.reset(polyline); } return *cached_polyline_; } @@ -133,7 +132,7 @@ void Curve::reverse() resetCache(); } -void Curve::manipulateControlPoint(uint idx, const Point& point) +void Curve::moveControlPoint(uint idx, const Point& point) { control_points_.row(idx) = point; resetCache(); @@ -189,8 +188,9 @@ Point Curve::valueAt(double t) const { if (N_ == 0) return {0, 0}; - Eigen::VectorXd power_basis = Eigen::pow(t, Eigen::ArrayXd::LinSpaced(N_, 0, N_ - 1)); - return (power_basis.transpose() * bernsteinCoeffs(N_) * control_points_).transpose(); + return (Eigen::VectorXd(Eigen::pow(t, Eigen::ArrayXd::LinSpaced(N_, 0, N_ - 1))).transpose() * bernsteinCoeffs(N_) * + control_points_) + .transpose(); } PointVector Curve::valueAt(const std::vector& t_vector) const @@ -200,7 +200,8 @@ PointVector Curve::valueAt(const std::vector& t_vector) const auto t_matrix = Eigen::Map(t_vector.data(), static_cast(t_vector.size())).replicate(1, N_); - auto p_matrix = Eigen::ArrayXd::LinSpaced(N_, 0, N_ - 1).transpose().replicate(static_cast(t_vector.size()), 1); + auto p_matrix = + Eigen::ArrayXd::LinSpaced(N_, 0, N_ - 1).transpose().replicate(static_cast(t_vector.size()), 1); Eigen::MatrixXd power_basis = Eigen::pow(t_matrix.array(), p_matrix.array()); Eigen::MatrixXd points_eigen = (power_basis * bernsteinCoeffs(N_) * control_points_); @@ -212,17 +213,17 @@ PointVector Curve::valueAt(const std::vector& t_vector) const double Curve::curvatureAt(double t) const { - Point d1 = derivativeAt(t); - Point d2 = derivativeAt(2, t); + Vector d1 = derivativeAt(t); + Vector d2 = derivativeAt(2, t); return (d1.x() * d2.y() - d1.y() * d2.x()) / std::pow(d1.norm(), 3); } double Curve::curvatureDerivativeAt(double t) const { - auto d1 = derivativeAt(t); - auto d2 = derivativeAt(2, t); - auto d3 = derivativeAt(3, t); + Vector d1 = derivativeAt(t); + Vector d2 = derivativeAt(2, t); + Vector d3 = derivativeAt(3, t); return (d1.x() * d3.y() - d1.y() * d3.x()) / std::pow(d1.norm(), 3) - 3 * d1.dot(d2) * (d1.x() * d2.y() - d1.y() * d2.x()) / std::pow(d1.norm(), 5); @@ -246,11 +247,10 @@ const Curve& Curve::derivative() const { if (!cached_derivative_) { - (const_cast(this)) - ->cached_derivative_.reset( - N_ == 1 ? new Curve(PointVector{Point(0, 0)}) - : new Curve( - ((N_ - 1) * (control_points_.bottomRows(N_ - 1) - control_points_.topRows(N_ - 1))).eval())); + const_cast(this)->cached_derivative_.reset( + N_ == 1 + ? new Curve(PointVector{Point(0, 0)}) + : new Curve(((N_ - 1) * (control_points_.bottomRows(N_ - 1) - control_points_.topRows(N_ - 1))).eval())); } return *cached_derivative_; } @@ -279,15 +279,25 @@ std::vector Curve::roots() const std::vector roots_X, roots_Y; Eigen::MatrixXd bezier_polynomial = bernsteinCoeffs(N_) * control_points_; Eigen::PolynomialSolver poly_solver; - poly_solver.compute(trimZeroes(bezier_polynomial.col(0))); - poly_solver.realRoots(roots_X); - poly_solver.compute(trimZeroes(bezier_polynomial.col(1))); - poly_solver.realRoots(roots_Y); + auto trimmed = trimZeroes(bezier_polynomial.col(0)); + if (trimmed.size()) + { + poly_solver.compute(trimmed); + poly_solver.realRoots(roots_X); + } + trimmed = trimZeroes(bezier_polynomial.col(1)); + if (trimmed.size()) + { + poly_solver.compute(trimmed); + poly_solver.realRoots(roots_Y); + } roots->reserve(roots_X.size() + roots_Y.size()); - std::copy_if(std::make_move_iterator(roots_X.begin()), std::make_move_iterator(roots_X.end()), - std::back_inserter(*roots), [](double t) { return t >= 0 && t <= 1; }); - std::copy_if(std::make_move_iterator(roots_Y.begin()), std::make_move_iterator(roots_Y.end()), - std::back_inserter(*roots), [](double t) { return t >= 0 && t <= 1; }); + std::copy_if(std::execution::par_unseq, std::make_move_iterator(roots_X.begin()), + std::make_move_iterator(roots_X.end()), std::back_inserter(*roots), + [](double t) { return t >= 0 && t <= 1; }); + std::copy_if(std::execution::par_unseq, std::make_move_iterator(roots_Y.begin()), + std::make_move_iterator(roots_Y.end()), std::back_inserter(*roots), + [](double t) { return t >= 0 && t <= 1; }); } const_cast(this)->cached_roots_.reset(roots); } @@ -306,9 +316,9 @@ BoundingBox Curve::boundingBox() const extremes.emplace_back(control_points_.row(N_ - 1)); // find mininum and maximum along each axis - auto x_extremes = std::minmax_element(extremes.begin(), extremes.end(), + auto x_extremes = std::minmax_element(std::execution::par_unseq, extremes.begin(), extremes.end(), [](const Point& lhs, const Point& rhs) { return lhs.x() < rhs.x(); }); - auto y_extremes = std::minmax_element(extremes.begin(), extremes.end(), + auto y_extremes = std::minmax_element(std::execution::par_unseq, extremes.begin(), extremes.end(), [](const Point& lhs, const Point& rhs) { return lhs.y() < rhs.y(); }); const_cast(this)->cached_bounding_box_.reset(new BoundingBox( Point(x_extremes.first->x(), y_extremes.first->y()), Point(x_extremes.second->x(), y_extremes.second->y()))); @@ -373,8 +383,7 @@ PointVector Curve::intersections(const Curve& curve, double epsilon) const while (!subcurve_pairs.empty()) { - Eigen::MatrixX2d part_a = std::get<0>(subcurve_pairs.back()); - Eigen::MatrixX2d part_b = std::get<1>(subcurve_pairs.back()); + auto [part_a, part_b] = subcurve_pairs.back(); subcurve_pairs.pop_back(); BoundingBox bbox1 = bbox(part_a); @@ -390,7 +399,7 @@ PointVector Curve::intersections(const Curve& curve, double epsilon) const // segments converged, check if not already found and add new Point new_point = bbox1.center(); if (points_of_intersection.end() == - std::find_if(points_of_intersection.begin(), points_of_intersection.end(), + std::find_if(std::execution::par_unseq, points_of_intersection.begin(), points_of_intersection.end(), [new_point, epsilon](const Point& point) { return (point - new_point).norm() < epsilon; })) points_of_intersection.emplace_back(new_point); @@ -484,10 +493,9 @@ double Curve::projectPoint(const Point& point) const std::vector Curve::projectPoint(const PointVector& point_vector) const { - std::vector t_vector; - t_vector.reserve(point_vector.size()); - for (const auto& point : point_vector) - t_vector.emplace_back(projectPoint(point)); + std::vector t_vector(point_vector.size()); + std::transform(std::execution::par_unseq, point_vector.begin(), point_vector.end(), t_vector.begin(), + [this](const Point& point) { return projectPoint(point); }); return t_vector; } @@ -495,10 +503,9 @@ double Curve::distance(const Point& point) const { return (point - valueAt(proje std::vector Curve::distance(const PointVector& point_vector) const { - std::vector dist_vector; - dist_vector.reserve(point_vector.size()); - for (const auto& point : point_vector) - dist_vector.emplace_back(distance(point)); + std::vector dist_vector(point_vector.size()); + std::transform(std::execution::par_unseq, point_vector.begin(), point_vector.end(), dist_vector.begin(), + [this](const Point& point) { return distance(point); }); return dist_vector; } diff --git a/src/polycurve.cpp b/src/polycurve.cpp index 54843b5..351fb50 100644 --- a/src/polycurve.cpp +++ b/src/polycurve.cpp @@ -1,6 +1,7 @@ #include "Bezier/polycurve.h" -#include "Bezier/bezier.h" +#include +#include #include #include @@ -16,15 +17,7 @@ void PolyCurve::insertFront(Curve curve) { curves_.emplace_front(curve); } void PolyCurve::insertBack(Curve curve) { curves_.emplace_back(curve); } -void PolyCurve::removeAt(uint idx) -{ - if (idx == 0) - removeFirst(); - else if (idx == size() - 1) - removeBack(); - else - curves_.erase(curves_.begin() + idx); -} +void PolyCurve::removeAt(uint idx) { curves_.erase(curves_.begin() + idx); } void PolyCurve::removeFirst() { curves_.pop_front(); } @@ -39,6 +32,7 @@ uint PolyCurve::curveIdx(double t) const } Curve& PolyCurve::curve(uint idx) { return curves_[idx]; } + const Curve& PolyCurve::curve(uint idx) const { return curves_[idx]; } std::deque& PolyCurve::curves() { return curves_; } @@ -122,12 +116,12 @@ PointVector PolyCurve::controlPoints() const return cp; } -void PolyCurve::manipulateControlPoint(uint idx, const Point& point) +void PolyCurve::moveControlPoint(uint idx, const Point& point) { for (auto& curve : curves_) if (idx <= curve.order()) { - curve.manipulateControlPoint(idx, point); + curve.moveControlPoint(idx, point); break; } else @@ -237,9 +231,18 @@ double PolyCurve::projectPoint(const Point& point) const std::vector PolyCurve::projectPoint(const PointVector& point_vector) const { - std::vector t_vector; - t_vector.reserve(point_vector.size()); - for (const auto& point : point_vector) - t_vector.emplace_back(projectPoint(point)); + std::vector t_vector(point_vector.size()); + std::transform(std::execution::par_unseq, point_vector.begin(), point_vector.end(), t_vector.begin(), + [this](const Point& point) { return projectPoint(point); }); return t_vector; } + +double PolyCurve::distance(const Point& point) const { return (point - valueAt(projectPoint(point))).norm(); } + +std::vector PolyCurve::distance(const PointVector& point_vector) const +{ + std::vector dist_vector(point_vector.size()); + std::transform(std::execution::par_unseq, point_vector.begin(), point_vector.end(), dist_vector.begin(), + [this](const Point& point) { return distance(point); }); + return dist_vector; +}