Raise C++ version to c++17

- update STL function to c++17 version
- fix root solving for axis parallel linear curves

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
   $<BUILD_INTERFACE:${Bezier_SOURCE_DIR}/include>
   $<INSTALL_INTERFACE:include>

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

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

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();
   }

example/qpolycurve.h

@@ -5,15 +5,16 @@
 
 #include "Bezier/declarations.h"
 #include "Bezier/polycurve.h"
+#include "Bezier/bezier.h"
 
 class qPolyCurve : public QGraphicsItem, public Bezier::PolyCurve
 {
 private:
   bool draw_control_points = false;
   bool draw_curvature_radious = false;
 public:
-  qPolyCurve(const std::vector<Bezier::Curve>& curve_list) : QGraphicsItem(), Bezier::PolyCurve(curve_list) {}
-  qPolyCurve(const Bezier::Curve& curve) : QGraphicsItem(), Bezier::PolyCurve(curve) {}
+  qPolyCurve(const std::deque<Bezier::Curve>& curve_list) : QGraphicsItem(), Bezier::PolyCurve(curve_list) {}
+  qPolyCurve(const Bezier::Curve& curve) : QGraphicsItem(), Bezier::PolyCurve(std::deque<Bezier::Curve>{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;

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'

include/Bezier/polycurve.h

@@ -19,7 +19,8 @@
 
 #include <deque>
 
-#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<double> 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<double> distance(const PointVector& point_vector) const;
+
 protected:
   /// Structure for holding underlying Bezier curves
   std::deque<Curve> curves_;

src/bezier.cpp

@@ -1,6 +1,8 @@
 #include "Bezier/bezier.h"
 #include "Bezier/legendre_gauss.h"
 
+#include <execution>
+#include <functional>
 #include <numeric>
 
 #include <unsupported/Eigen/MatrixFunctions>
@@ -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<Point, Point> 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<Curve*>(this))->cached_polyline_params_ = {smoothness, precision};
-    (const_cast<Curve*>(this))->cached_polyline_.reset(polyline);
+    const_cast<Curve*>(this)->cached_polyline_params_ = {smoothness, precision};
+    const_cast<Curve*>(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<double>& t_vector) const
@@ -200,7 +200,8 @@ PointVector Curve::valueAt(const std::vector<double>& t_vector) const
 
   auto t_matrix =
       Eigen::Map<const Eigen::VectorXd>(t_vector.data(), static_cast<int>(t_vector.size())).replicate(1, N_);
-  auto p_matrix = Eigen::ArrayXd::LinSpaced(N_, 0, N_ - 1).transpose().replicate(static_cast<int>(t_vector.size()), 1);
+  auto p_matrix =
+      Eigen::ArrayXd::LinSpaced(N_, 0, N_ - 1).transpose().replicate(static_cast<int>(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<double>& 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<Curve*>(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<Curve*>(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<double> Curve::roots() const
       std::vector<double> roots_X, roots_Y;
       Eigen::MatrixXd bezier_polynomial = bernsteinCoeffs(N_) * control_points_;
       Eigen::PolynomialSolver<double, Eigen::Dynamic> 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<Curve*>(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<Curve*>(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,21 +493,19 @@ double Curve::projectPoint(const Point& point) const
 
 std::vector<double> Curve::projectPoint(const PointVector& point_vector) const
 {
-  std::vector<double> t_vector;
-  t_vector.reserve(point_vector.size());
-  for (const auto& point : point_vector)
-    t_vector.emplace_back(projectPoint(point));
+  std::vector<double> 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 Curve::distance(const Point& point) const { return (point - valueAt(projectPoint(point))).norm(); }
 
 std::vector<double> Curve::distance(const PointVector& point_vector) const
 {
-  std::vector<double> dist_vector;
-  dist_vector.reserve(point_vector.size());
-  for (const auto& point : point_vector)
-    dist_vector.emplace_back(distance(point));
+  std::vector<double> 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;
 }