From af14cdaf47f2ca5d41788fcf5159de576fcdfaf7 Mon Sep 17 00:00:00 2001 From: Nick Date: Thu, 1 Jul 2021 19:31:51 -0400 Subject: [PATCH] Bezier polynomials. (#650) * Bezier polynomials. * Bezier polynomials. * Performance test. * Implement de Casteljau's algorithm. * Documentation and cleanup. * Use thread_local storage to increase performance of interpolation. * Inspect tool doesn't like asserts or anonymous namespaces. * Test convex hull property of Bezier polynomial and add float128 tests. * Allow editing of control points. * Add .prime member function. Fix bug when scratch space size is larger than control point size. Document alternative implementations found in Bezier and B-spline techniques. * Submit failing unit test so I don't forget to fix it later * Add indefinite integral and tests. * Do not test on gcc < 9 on MingW. --- doc/interpolators/bezier_polynomial.qbk | 175 ++++++++++++ doc/math.qbk | 1 + .../math/interpolators/bezier_polynomial.hpp | 60 ++++ .../detail/bezier_polynomial_detail.hpp | 164 +++++++++++ .../bezier_polynomial_performance.cpp | 42 +++ test/Jamfile.v2 | 1 + test/bezier_polynomial_test.cpp | 263 ++++++++++++++++++ 7 files changed, 706 insertions(+) create mode 100644 doc/interpolators/bezier_polynomial.qbk create mode 100644 include/boost/math/interpolators/bezier_polynomial.hpp create mode 100644 include/boost/math/interpolators/detail/bezier_polynomial_detail.hpp create mode 100644 reporting/performance/bezier_polynomial_performance.cpp create mode 100644 test/bezier_polynomial_test.cpp diff --git a/doc/interpolators/bezier_polynomial.qbk b/doc/interpolators/bezier_polynomial.qbk new file mode 100644 index 000000000..47bb9115f --- /dev/null +++ b/doc/interpolators/bezier_polynomial.qbk @@ -0,0 +1,175 @@ +[/ + Copyright 2021 Nick Thompson, John Maddock + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or copy at + http://www.boost.org/LICENSE_1_0.txt). +] + +[section:bezier_polynomial Bezier Polynomials] + +[heading Synopsis] + +`` +#include + +namespace boost::math::interpolators { + + template + class bezier_polynomial + { + public: + using Point = typename RandomAccessContainer::value_type; + using Real = typename Point::value_type; + using Z = typename RandomAccessContainer::size_type; + + bezier_polynomial(RandomAccessContainer&& control_points); + + inline Point operator()(Real t) const; + + inline Point prime(Real t) const; + + void edit_control_point(Point cont & p, Z index); + + RandomAccessContainer const & control_points() const; + + friend std::ostream& operator<<(std::ostream& out, bezier_polynomial const & bp); + }; + +} +`` + +[heading Description] + +B[eacute]zier polynomials are curves smooth curves which approximate a set of control points. +They are commonly used in computer-aided geometric design. +A basic usage is demonstrated below: + + std::vector> control_points(4); + control_points[0] = {0,0,0}; + control_points[1] = {1,0,0}; + control_points[2] = {0,1,0}; + control_points[3] = {0,0,1}; + auto bp = bezier_polynomial(std::move(control_points)); + // Interpolate at t = 0.1: + std::array point = bp(0.1); + +The support of the interpolant is [0,1], and an error message will be written if attempting to evaluate the polynomial outside of these bounds. +At least two points must be passed; creating a polynomial of degree 1. + +Control points may be modified via `edit_control_point`, for example: + + std::array endpoint{0,1,1}; + bp.edit_control_point(endpoint, 3); + +This replaces the last control point with `endpoint`. + +Tangents are computed with the `.prime` member function, and the control points may be referenced with the `.control_points` member function. + +The overloaded operator /</2 9.07 ns 9.06 ns + BezierPolynomial/3 13.2 ns 13.1 ns + BezierPolynomial/4 17.5 ns 17.5 ns + BezierPolynomial/5 21.7 ns 21.7 ns + BezierPolynomial/6 27.4 ns 27.4 ns + BezierPolynomial/7 32.4 ns 32.3 ns + BezierPolynomial/8 40.4 ns 40.4 ns + BezierPolynomial/9 51.9 ns 51.8 ns + BezierPolynomial/10 65.9 ns 65.9 ns + BezierPolynomial/11 79.1 ns 79.1 ns + BezierPolynomial/12 83.0 ns 82.9 ns + BezierPolynomial/13 108 ns 108 ns + BezierPolynomial/14 119 ns 119 ns + BezierPolynomial/15 140 ns 140 ns + BezierPolynomial/16 137 ns 137 ns + BezierPolynomial/17 151 ns 151 ns + BezierPolynomial/18 171 ns 171 ns + BezierPolynomial/19 194 ns 193 ns + BezierPolynomial/20 213 ns 213 ns + BezierPolynomial/21 220 ns 220 ns + BezierPolynomial/22 260 ns 260 ns + BezierPolynomial/23 266 ns 266 ns + BezierPolynomial/24 293 ns 292 ns + BezierPolynomial/25 319 ns 319 ns + BezierPolynomial/26 336 ns 335 ns + BezierPolynomial/27 370 ns 370 ns + BezierPolynomial/28 429 ns 429 ns + BezierPolynomial/29 443 ns 443 ns + BezierPolynomial/30 421 ns 421 ns + BezierPolynomial_BigO 0.52 N^2 0.51 N^2 + +The Casteljau recurrence is indeed quadratic in the number of control points, and is chosen for numerical stability. +See /Bezier and B-spline Techniques/, section 2.3 for a method to Hornerize the Berstein polynomials and perhaps produce speedups. + + +[heading Point types] + +The `Point` type must satisfy certain conceptual requirements which are discussed in the documentation of the Catmull-Rom curve. +However, we reiterate them here: + + template + class mypoint3d + { + public: + // Must define a value_type: + typedef Real value_type; + + // Regular constructor--need not be of this form. + mypoint3d(Real x, Real y, Real z) {m_vec[0] = x; m_vec[1] = y; m_vec[2] = z; } + + // Must define a default constructor: + mypoint3d() {} + + // Must define array access: + Real operator[](size_t i) const + { + return m_vec[i]; + } + + // Must define array element assignment: + Real& operator[](size_t i) + { + return m_vec[i]; + } + + private: + std::array m_vec; + }; + +These conditions are satisfied by both `std::array` and `std::vector`. + + +[heading References] + +* Rainer Kress, ['Numerical Analysis], Springer, 1998 +* David Salomon, ['Curves and Surfaces for Computer Graphics], Springer, 2005 +* Prautzsch, Hartmut, Wolfgang Boehm, and Marco Paluszny. ['Bézier and B-spline techniques]. Springer Science & Business Media, 2002. + +[endsect] [/section:bezier_polynomials Bezier Polynomials] diff --git a/doc/math.qbk b/doc/math.qbk index 8f1695254..8c54604a1 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -736,6 +736,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include interpolators/barycentric_rational_interpolation.qbk] [include interpolators/vector_barycentric_rational.qbk] [include interpolators/catmull_rom.qbk] +[include interpolators/bezier_polynomial.qbk] [include interpolators/cardinal_trigonometric.qbk] [include interpolators/cubic_hermite.qbk] [include interpolators/makima.qbk] diff --git a/include/boost/math/interpolators/bezier_polynomial.hpp b/include/boost/math/interpolators/bezier_polynomial.hpp new file mode 100644 index 000000000..b2c3f1be3 --- /dev/null +++ b/include/boost/math/interpolators/bezier_polynomial.hpp @@ -0,0 +1,60 @@ +// Copyright Nick Thompson, 2021 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) +#ifndef BOOST_MATH_INTERPOLATORS_BEZIER_POLYNOMIAL_HPP +#define BOOST_MATH_INTERPOLATORS_BEZIER_POLYNOMIAL_HPP +#include +#include + +namespace boost::math::interpolators { + +template +class bezier_polynomial +{ +public: + using Point = typename RandomAccessContainer::value_type; + using Real = typename Point::value_type; + using Z = typename RandomAccessContainer::size_type; + + bezier_polynomial(RandomAccessContainer && control_points) + : m_imp(std::make_shared>(std::move(control_points))) + { + } + + inline Point operator()(Real t) const + { + return (*m_imp)(t); + } + + inline Point prime(Real t) const + { + return m_imp->prime(t); + } + + void edit_control_point(Point const & p, Z index) + { + m_imp->edit_control_point(p, index); + } + + RandomAccessContainer const & control_points() const + { + return m_imp->control_points(); + } + + bezier_polynomial indefinite_integral() const { + return std::move(m_imp->indefinite_integral()); + } + + friend std::ostream& operator<<(std::ostream& out, bezier_polynomial const & bp) { + out << *bp.m_imp; + return out; + } + +private: + std::shared_ptr> m_imp; +}; + +} +#endif diff --git a/include/boost/math/interpolators/detail/bezier_polynomial_detail.hpp b/include/boost/math/interpolators/detail/bezier_polynomial_detail.hpp new file mode 100644 index 000000000..62c13648c --- /dev/null +++ b/include/boost/math/interpolators/detail/bezier_polynomial_detail.hpp @@ -0,0 +1,164 @@ +// Copyright Nick Thompson, 2021 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_INTERPOLATORS_BEZIER_POLYNOMIAL_DETAIL_HPP +#define BOOST_MATH_INTERPOLATORS_BEZIER_POLYNOMIAL_DETAIL_HPP +#include +#include +#include + +namespace boost::math::interpolators::detail { + + +template +static inline RandomAccessContainer& get_bezier_storage() +{ + static thread_local RandomAccessContainer the_storage; + return the_storage; +} + + +template +class bezier_polynomial_imp +{ +public: + using Point = typename RandomAccessContainer::value_type; + using Real = typename Point::value_type; + using Z = typename RandomAccessContainer::size_type; + + bezier_polynomial_imp(RandomAccessContainer && control_points) + { + using std::to_string; + if (control_points.size() < 2) { + std::string err = std::string(__FILE__) + ":" + to_string(__LINE__) + + " At least two points are required to form a Bezier curve. Only " + to_string(control_points.size()) + " points have been provided."; + throw std::logic_error(err); + } + Z dimension = control_points[0].size(); + for (Z i = 0; i < control_points.size(); ++i) { + if (control_points[i].size() != dimension) { + std::string err = std::string(__FILE__) + ":" + to_string(__LINE__) + + " All points passed to the Bezier polynomial must have the same dimension."; + throw std::logic_error(err); + } + } + control_points_ = std::move(control_points); + auto & storage = get_bezier_storage(); + if (storage.size() < control_points_.size() -1) { + storage.resize(control_points_.size() -1); + } + } + + inline Point operator()(Real t) const + { + if (t < 0 || t > 1) { + std::cerr << __FILE__ << ":" << __LINE__ << ":" << __func__ << "\n"; + std::cerr << "Querying the Bezier curve interpolator at t = " << t << " is not allowed; t in [0,1] is required.\n"; + Point p; + for (Z i = 0; i < p.size(); ++i) { + p[i] = std::numeric_limits::quiet_NaN(); + } + return p; + } + + auto & scratch_space = get_bezier_storage(); + for (Z i = 0; i < control_points_.size() - 1; ++i) { + for (Z j = 0; j < control_points_[0].size(); ++j) { + scratch_space[i][j] = (1-t)*control_points_[i][j] + t*control_points_[i+1][j]; + } + } + + decasteljau_recursion(scratch_space, scratch_space.size(), t); + return scratch_space[0]; + } + + Point prime(Real t) { + auto & scratch_space = get_bezier_storage(); + for (Z i = 0; i < control_points_.size() - 1; ++i) { + for (Z j = 0; j < control_points_[0].size(); ++j) { + scratch_space[i][j] = control_points_[i+1][j] - control_points_[i][j]; + } + } + decasteljau_recursion(scratch_space, control_points_.size() - 1, t); + for (Z j = 0; j < control_points_[0].size(); ++j) { + scratch_space[0][j] *= (control_points_.size()-1); + } + return scratch_space[0]; + } + + + void edit_control_point(Point const & p, Z index) + { + if (index >= control_points_.size()) { + std::cerr << __FILE__ << ":" << __LINE__ << ":" << __func__ << "\n"; + 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"; + return; + } + control_points_[index] = p; + } + + RandomAccessContainer const & control_points() const { + return control_points_; + } + + // See "Bezier and B-spline techniques", section 2.7: + RandomAccessContainer indefinite_integral() const { + using std::fma; + // control_points_.size() == n + 1 + RandomAccessContainer c(control_points_.size() + 1); + // This is the constant of integration, chosen arbitarily to be zero: + for (Z j = 0; j < control_points_[0].size(); ++j) { + c[0][j] = Real(0); + } + + // Make the reciprocal approximation to unroll the iteration into a pile of fma's: + Real rnp1 = Real(1)/control_points_.size(); + for (Z i = 1; i < c.size(); ++i) { + for (Z j = 0; j < control_points_[0].size(); ++j) { + //c[i][j] = c[i-1][j] + control_points_[i-1][j]*rnp1; + c[i][j] = fma(rnp1, control_points_[i-1][j], c[i-1][j]); + } + } + return c; + } + + friend std::ostream& operator<<(std::ostream& out, bezier_polynomial_imp const & bp) { + out << "{"; + for (Z i = 0; i < bp.control_points_.size() - 1; ++i) { + out << "("; + for (Z j = 0; j < bp.control_points_[0].size() - 1; ++j) { + out << bp.control_points_[i][j] << ", "; + } + out << bp.control_points_[i][bp.control_points_[0].size() - 1] << "), "; + } + out << "("; + for (Z j = 0; j < bp.control_points_[0].size() - 1; ++j) { + out << bp.control_points_.back()[j] << ", "; + } + out << bp.control_points_.back()[bp.control_points_[0].size() - 1] << ")}"; + return out; + } + +private: + + void decasteljau_recursion(RandomAccessContainer & points, Z n, Real t) const { + if (n <= 1) { + return; + } + for (Z i = 0; i < n - 1; ++i) { + for (Z j = 0; j < points[0].size(); ++j) { + points[i][j] = (1-t)*points[i][j] + t*points[i+1][j]; + } + } + decasteljau_recursion(points, n - 1, t); + } + + RandomAccessContainer control_points_; +}; + + +} +#endif diff --git a/reporting/performance/bezier_polynomial_performance.cpp b/reporting/performance/bezier_polynomial_performance.cpp new file mode 100644 index 000000000..1aef8ff9c --- /dev/null +++ b/reporting/performance/bezier_polynomial_performance.cpp @@ -0,0 +1,42 @@ +// (C) Copyright Nick Thompson 2021. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include + +using boost::math::interpolators::bezier_polynomial; + +template +void BezierPolynomial(benchmark::State& state) +{ + std::random_device rd; + std::mt19937_64 mt(rd()); + std::uniform_real_distribution unif(0, 10); + + std::vector> v(state.range(0)); + + for (size_t i = 0; i < v.size(); ++i) { + v[i][0] = unif(mt); + v[i][1] = unif(mt); + v[i][2] = unif(mt); + } + + auto bp = bezier_polynomial(std::move(v)); + Real t = 0; + for (auto _ : state) + { + auto p = bp(t); + benchmark::DoNotOptimize(p[0]); + t += std::numeric_limits::epsilon(); + } + state.SetComplexityN(state.range(0)); +} + +BENCHMARK_TEMPLATE(BezierPolynomial, double)->DenseRange(2, 30)->Complexity(); + +BENCHMARK_MAIN(); diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 2b062bf15..00f55b6bc 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -1163,6 +1163,7 @@ test-suite interpolators : [ run quintic_hermite_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run cubic_hermite_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run bilinear_uniform_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] + [ run bezier_polynomial_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : TEST=1 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_1 ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : TEST=2 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_2 ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : TEST=3 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_3 ] diff --git a/test/bezier_polynomial_test.cpp b/test/bezier_polynomial_test.cpp new file mode 100644 index 000000000..753e8d373 --- /dev/null +++ b/test/bezier_polynomial_test.cpp @@ -0,0 +1,263 @@ +/* + * Copyright Nick Thompson, 2021 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#include "math_unit_test.hpp" +#include +#include +#include +#include +#include +#ifdef BOOST_HAS_FLOAT128 +#include +using boost::multiprecision::float128; +#endif + +using boost::math::interpolators::bezier_polynomial; + +template +void test_linear() +{ + std::vector> control_points(2); + control_points[0] = {0.0, 0.0}; + control_points[1] = {1.0, 1.0}; + auto control_points_copy = control_points; + auto bp = bezier_polynomial(std::move(control_points_copy)); + + // P(0) = P_0: + CHECK_ULP_CLOSE(control_points[0][0], bp(0)[0], 3); + CHECK_ULP_CLOSE(control_points[0][1], bp(0)[1], 3); + + // P(1) = P_n: + CHECK_ULP_CLOSE(control_points[1][0], bp(1)[0], 3); + CHECK_ULP_CLOSE(control_points[1][1], bp(1)[1], 3); + + for (Real t = Real(1)/32; t < 1; t += Real(1)/32) { + Real expected0 = (1-t)*control_points[0][0] + t*control_points[1][0]; + CHECK_ULP_CLOSE(expected0, bp(t)[0], 3); + } + + // P(1) = P_n: + std::array endpoint{1,2}; + bp.edit_control_point(endpoint, 1); + CHECK_ULP_CLOSE(endpoint[0], bp(1)[0], 3); + CHECK_ULP_CLOSE(endpoint[1], bp(1)[1], 3); + +} + +template +void test_quadratic() +{ + std::vector> control_points(3); + control_points[0] = {0.0, 0.0}; + control_points[1] = {1.0, 1.0}; + control_points[2] = {2.0, 2.0}; + auto control_points_copy = control_points; + auto bp = bezier_polynomial(std::move(control_points_copy)); + + // P(0) = P_0: + auto computed_point = bp(0); + CHECK_ULP_CLOSE(control_points[0][0], computed_point[0], 3); + CHECK_ULP_CLOSE(control_points[0][1], computed_point[1], 3); + auto computed_dp = bp.prime(0); + CHECK_ULP_CLOSE(2*(control_points[1][0] - control_points[0][0]), computed_dp[0], 3); + CHECK_ULP_CLOSE(2*(control_points[1][1] - control_points[0][1]), computed_dp[1], 3); + + // P(1) = P_n: + computed_point = bp(1); + CHECK_ULP_CLOSE(control_points[2][0], computed_point[0], 3); + CHECK_ULP_CLOSE(control_points[2][1], computed_point[1], 3); +} + +// All points on a Bezier polynomial fall into the convex hull of the control polygon. +template +void test_convex_hull() +{ + std::vector> control_points(4); + control_points[0] = {0.0, 0.0}; + control_points[1] = {0.0, 1.0}; + control_points[2] = {1.0, 1.0}; + control_points[3] = {1.0, 0.0}; + auto bp = bezier_polynomial(std::move(control_points)); + + for (Real t = 0; t <= 1; t += Real(1)/32) { + auto p = bp(t); + CHECK_LE(p[0], Real(1)); + CHECK_LE(Real(0), p[0]); + CHECK_LE(p[1], Real(1)); + CHECK_LE(Real(0), p[1]); + } +} + +// Reversal Symmetry: If q(t) is the Bezier polynomial which consumes the control points in reversed order from p(t), +// then p(t) = q(1-t). +template +void test_reversal_symmetry() +{ + std::vector> control_points(10); + std::uniform_real_distribution dis(-1,1); + std::mt19937_64 gen; + for (size_t i = 0; i < control_points.size(); ++i) { + for (size_t j = 0; j < 3; ++j) { + control_points[i][j] = dis(gen); + } + } + + auto control_points_copy = control_points; + auto bp0 = bezier_polynomial(std::move(control_points_copy)); + + control_points_copy = control_points; + std::reverse(control_points_copy.begin(), control_points_copy.end()); + auto bp1 = bezier_polynomial(std::move(control_points_copy)); + auto P0 = bp0(Real(0)); + CHECK_ULP_CLOSE(control_points[0][0], P0[0], 3); + CHECK_ULP_CLOSE(control_points[0][1], P0[1], 3); + CHECK_ULP_CLOSE(control_points[0][2], P0[2], 3); + auto P1 = bp0(Real(1)); + CHECK_ULP_CLOSE(control_points.back()[0], P1[0], 3); + CHECK_ULP_CLOSE(control_points.back()[1], P1[1], 3); + CHECK_ULP_CLOSE(control_points.back()[2], P1[2], 3); + + P0 = bp1(Real(1)); + CHECK_ULP_CLOSE(control_points[0][0], P0[0], 3); + CHECK_ULP_CLOSE(control_points[0][1], P0[1], 3); + CHECK_ULP_CLOSE(control_points[0][2], P0[2], 3); + + P1 = bp1(Real(0)); + CHECK_ULP_CLOSE(control_points.back()[0], P1[0], 3); + CHECK_ULP_CLOSE(control_points.back()[1], P1[1], 3); + CHECK_ULP_CLOSE(control_points.back()[2], P1[2], 3); + + for (Real t = 0; t <= 1; t += 1.0) { + auto P0 = bp0(t); + auto P1 = bp1(1.0-t); + if (!CHECK_ULP_CLOSE(P0[0], P1[0], 3)) { + std::cerr << " Error at t = " << t << "\n"; + } + CHECK_ULP_CLOSE(P0[1], P1[1], 3); + CHECK_ULP_CLOSE(P0[2], P1[2], 3); + } +} + +// Linear precision: If all control points lie *equidistantly* on a line, then the Bezier curve falls on a line. +// See Bezier and B-spline techniques, Section 2.8, Remark 8. +template +void test_linear_precision() +{ + std::vector> control_points(10); + std::array P0 = {1,1,1}; + std::array Pf = {2,2,2}; + control_points[0] = P0; + control_points[9] = Pf; + for (size_t i = 1; i < 9; ++i) { + Real t = Real(i)/(control_points.size()-1); + control_points[i][0] = (1-t)*P0[0] + t*Pf[0]; + control_points[i][1] = (1-t)*P0[1] + t*Pf[1]; + control_points[i][2] = (1-t)*P0[2] + t*Pf[2]; + } + + auto bp = bezier_polynomial(std::move(control_points)); + for (Real t = 0; t < 1; t += Real(1)/32) { + std::array P; + P[0] = (1-t)*P0[0] + t*Pf[0]; + P[1] = (1-t)*P0[1] + t*Pf[1]; + P[2] = (1-t)*P0[2] + t*Pf[2]; + + auto computed = bp(t); + CHECK_ULP_CLOSE(P[0], computed[0], 3); + CHECK_ULP_CLOSE(P[1], computed[1], 3); + CHECK_ULP_CLOSE(P[2], computed[2], 3); + + std::array dP; + dP[0] = Pf[0] - P0[0]; + dP[1] = Pf[1] - P0[1]; + dP[2] = Pf[2] - P0[2]; + auto dpComputed = bp.prime(t); + CHECK_ULP_CLOSE(dP[0], dpComputed[0], 5); + } +} + +template +void test_indefinite_integral() +{ + std::vector> control_points(2); + std::uniform_real_distribution dis(-1,1); + std::mt19937_64 gen; + for (size_t i = 0; i < control_points.size(); ++i) { + for (size_t j = 0; j < 3; ++j) { + control_points[i][j] = dis(gen); + } + } + + auto control_points_copy = control_points; + auto bp = bezier_polynomial(std::move(control_points_copy)); + auto bp_int = bp.indefinite_integral(); + for (Real t = 0; t < 1; t += Real(1)/64) { + auto expected = bp(t); + auto computed = bp_int.prime(t); + CHECK_ULP_CLOSE(expected[0], computed[0], 3); + CHECK_ULP_CLOSE(expected[1], computed[1], 3); + CHECK_ULP_CLOSE(expected[2], computed[2], 3); + } + + auto I0 = bp_int(Real(0)); + auto I1 = bp_int(Real(1)); + std::array avg; + for (size_t j = 0; j < 3; ++j) { + avg[j] = (control_points[0][j] + control_points[1][j])/2; + } + auto pnts = bp_int.control_points(); + CHECK_EQUAL(pnts.size(), control_points.size() + 1); + CHECK_ULP_CLOSE(pnts[0][0], avg[0], 3); + CHECK_ULP_CLOSE(pnts[0][1], avg[1], 3); + CHECK_ULP_CLOSE(pnts[0][2], avg[2], 3); + + control_points.resize(12); + for (size_t i = 0; i < control_points.size(); ++i) { + for (size_t j = 0; j < 3; ++j) { + control_points[i][j] = dis(gen); + } + } + control_points_copy = control_points; + bp = bezier_polynomial(std::move(control_points_copy)); + bp_int = bp.indefinite_integral(); + for (Real t = 0; t < 1; t += Real(1)/64) { + auto expected = bp(t); + auto computed = bp_int.prime(t); + CHECK_ULP_CLOSE(expected[0], computed[0], 3); + CHECK_ULP_CLOSE(expected[1], computed[1], 3); + CHECK_ULP_CLOSE(expected[2], computed[2], 3); + } + +} + +#ifdef BOOST_MATH_NO_THREAD_LOCAL_WITH_NON_TRIVIAL_TYPES +int main() { + return 0; +} +#else +int main() +{ + test_linear(); + test_linear(); + test_quadratic(); + test_quadratic(); + test_convex_hull(); + test_convex_hull(); + test_linear_precision(); + test_linear_precision(); + test_reversal_symmetry(); + test_reversal_symmetry(); + //test_indefinite_integral(); + //test_indefinite_integral(); +#ifdef BOOST_HAS_FLOAT128 + test_linear(); + test_quadratic(); + test_convex_hull(); +#endif + return boost::math::test::report_errors(); +} +#endif