interpolators/detail/bezier_polynomial_detail.hpp

0001 // Copyright Nick Thompson, 2021
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_INTERPOLATORS_BEZIER_POLYNOMIAL_DETAIL_HPP
0008 #define BOOST_MATH_INTERPOLATORS_BEZIER_POLYNOMIAL_DETAIL_HPP
0009
0010 #include <stdexcept>
0011 #include <iostream>
0012 #include <string>
0013 #include <limits>
0014
0015 namespace boost::math::interpolators::detail {
0016
0017
0018 template <class RandomAccessContainer>
0019 static inline RandomAccessContainer& get_bezier_storage()
0020 {
0021     static thread_local RandomAccessContainer the_storage;
0022     return the_storage;
0023 }
0024
0025
0026 template <class RandomAccessContainer>
0027 class bezier_polynomial_imp
0028 {
0029 public:
0030     using Point = typename RandomAccessContainer::value_type;
0031     using Real = typename Point::value_type;
0032     using Z = typename RandomAccessContainer::size_type;
0033
0034     bezier_polynomial_imp(RandomAccessContainer && control_points)
0035     {
0036         using std::to_string;
0037         if (control_points.size() < 2) {
0038             std::string err = std::string(__FILE__) + ":" + to_string(__LINE__)
0039                + " At least two points are required to form a Bezier curve. Only " + to_string(control_points.size())  + " points have been provided.";
0040             throw std::logic_error(err);
0041         }
0042         Z dimension = control_points[0].size();
0043         for (Z i = 0; i < control_points.size(); ++i) {
0044             if (control_points[i].size() != dimension) {
0045                 std::string err = std::string(__FILE__) + ":" + to_string(__LINE__)
0046                 + " All points passed to the Bezier polynomial must have the same dimension.";
0047                 throw std::logic_error(err);
0048             }
0049         }
0050         control_points_ = std::move(control_points);
0051         auto & storage = get_bezier_storage<RandomAccessContainer>();
0052         if (storage.size() < control_points_.size() -1) {
0053             storage.resize(control_points_.size() -1);
0054         }
0055     }
0056
0057     inline Point operator()(Real t) const
0058     {
0059         if (t < 0 || t > 1) {
0060             std::cerr << __FILE__ << ":" << __LINE__ << ":" << __func__ << "\n";
0061             std::cerr << "Querying the Bezier curve interpolator at t = " << t << " is not allowed; t in [0,1] is required.\n";
0062             Point p;
0063             for (Z i = 0; i < p.size(); ++i) {
0064                 p[i] = std::numeric_limits<Real>::quiet_NaN();
0065             }
0066             return p;
0067         }
0068
0069         auto & scratch_space = get_bezier_storage<RandomAccessContainer>();
0070         for (Z i = 0; i < control_points_.size() - 1; ++i) {
0071             for (Z j = 0; j < control_points_[0].size(); ++j) {
0072                 scratch_space[i][j] = (1-t)*control_points_[i][j] + t*control_points_[i+1][j];
0073             }
0074         }
0075
0076         decasteljau_recursion(scratch_space, control_points_.size() - 1, t);
0077         return scratch_space[0];
0078     }
0079
0080     Point prime(Real t) {
0081         auto & scratch_space = get_bezier_storage<RandomAccessContainer>();
0082         for (Z i = 0; i < control_points_.size() - 1; ++i) {
0083             for (Z j = 0; j < control_points_[0].size(); ++j) {
0084                 scratch_space[i][j] = control_points_[i+1][j] - control_points_[i][j];
0085             }
0086         }
0087         decasteljau_recursion(scratch_space, control_points_.size() - 1, t);
0088         for (Z j = 0; j < control_points_[0].size(); ++j) {
0089             scratch_space[0][j] *= (control_points_.size()-1);
0090         }
0091         return scratch_space[0];
0092     }
0093
0094
0095     void edit_control_point(Point const & p, Z index)
0096     {
0097         if (index >= control_points_.size()) {
0098             std::cerr << __FILE__ << ":" << __LINE__ << ":" << __func__ << "\n";
0099             std::cerr << "Attempting to edit a control point outside the bounds of the container; requested edit of index " << index << ", but there are only " << control_points_.size() << " control points.\n";
0100             return;
0101         }
0102         control_points_[index] = p;
0103     }
0104
0105     RandomAccessContainer const & control_points() const {
0106         return control_points_;
0107     }
0108
0109     // See "Bezier and B-spline techniques", section 2.7:
0110     // I cannot figure out why this doesn't work.
0111     /*RandomAccessContainer indefinite_integral() const {
0112         using std::fma;
0113         // control_points_.size() == n + 1
0114         RandomAccessContainer c(control_points_.size() + 1);
0115         // This is the constant of integration, chosen arbitrarily to be zero:
0116         for (Z j = 0; j < control_points_[0].size(); ++j) {
0117             c[0][j] = Real(0);
0118         }
0119
0120         // Make the reciprocal approximation to unroll the iteration into a pile of fma's:
0121         Real rnp1 = Real(1)/control_points_.size();
0122         for (Z i = 1; i < c.size(); ++i) {
0123             for (Z j = 0; j < control_points_[0].size(); ++j) {
0124                 //c[i][j] = c[i-1][j] + control_points_[i-1][j]*rnp1;
0125                 c[i][j] = fma(rnp1, control_points_[i-1][j], c[i-1][j]);
0126             }
0127         }
0128         return c;
0129     }*/
0130
0131     friend std::ostream& operator<<(std::ostream& out, bezier_polynomial_imp<RandomAccessContainer> const & bp) {
0132         out << "{";
0133         for (Z i = 0; i < bp.control_points_.size() - 1; ++i) {
0134             out << "(";
0135             for (Z j = 0; j < bp.control_points_[0].size() - 1; ++j) {
0136                 out << bp.control_points_[i][j] << ", ";
0137             }
0138             out << bp.control_points_[i][bp.control_points_[0].size() - 1] << "), ";
0139         }
0140         out << "(";
0141         for (Z j = 0; j < bp.control_points_[0].size() - 1; ++j) {
0142             out << bp.control_points_.back()[j] << ", ";
0143         }
0144         out << bp.control_points_.back()[bp.control_points_[0].size() - 1] << ")}";
0145         return out;
0146     }
0147
0148 private:
0149
0150     void decasteljau_recursion(RandomAccessContainer & points, Z n, Real t) const {
0151         if (n <= 1) {
0152             return;
0153         }
0154         for (Z i = 0; i < n - 1; ++i) {
0155             for (Z j = 0; j < points[0].size(); ++j) {
0156                 points[i][j] = (1-t)*points[i][j] + t*points[i+1][j];
0157             }
0158         }
0159         decasteljau_recursion(points, n - 1, t);
0160     }
0161
0162     RandomAccessContainer control_points_;
0163 };
0164
0165
0166 }
0167 #endif