boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Bezier polynomials. (#650)

Browse Source 
* 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.
This commit is contained in:
Nick
2021-07-01 19:31:51 -04:00
committed by GitHub
parent 8094e3e6bf
commit af14cdaf47
7 changed files with 706 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

175
doc/interpolators/bezier_polynomial.qbk Normal file
View File

@@ -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 <boost/math/interpolators/bezier_polynomials.hpp>
namespace boost::math::interpolators {
template<RandomAccessContainer>
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<RandomAccessContainer> 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<std::array<double, 3>> 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<double, 3> 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<double, 3> 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 /<</ is disappointing: The control points are simply printed.
Rendering the Bezier and its convex hull seems to be the best "print" statement for it, but this is essentially impossible in modern terminals.
[heading Caveats]
Do not confuse the Bezier polynomial with a Bezier spline.
A Bezier spline has a fixed polynomial order and subdivides the curve into low-order polynomial segments.
/This is not a spline!/
Passing /n/ control points to the `bezier_polynomial` class creates a polynomial of degree n-1, whereas a Bezier spline has a fixed order independent of the number of control points.
Requires C++17.
[heading Performance]
The following performance numbers were generated for evaluating the Bezier-polynomial.
The evaluation of the interpolant is [bigo](/N/^2), as expected from de Casteljau's algorithm.
Run on (16 X 2300 MHz CPU s)
CPU Caches:
L1 Data 32 KiB (x8)
L1 Instruction 32 KiB (x8)
L2 Unified 256 KiB (x8)
L3 Unified 16384 KiB (x1)
---------------------------------------------------------
Benchmark Time CPU
---------------------------------------------------------
BezierPolynomial<double>/2 9.07 ns 9.06 ns
BezierPolynomial<double>/3 13.2 ns 13.1 ns
BezierPolynomial<double>/4 17.5 ns 17.5 ns
BezierPolynomial<double>/5 21.7 ns 21.7 ns
BezierPolynomial<double>/6 27.4 ns 27.4 ns
BezierPolynomial<double>/7 32.4 ns 32.3 ns
BezierPolynomial<double>/8 40.4 ns 40.4 ns
BezierPolynomial<double>/9 51.9 ns 51.8 ns
BezierPolynomial<double>/10 65.9 ns 65.9 ns
BezierPolynomial<double>/11 79.1 ns 79.1 ns
BezierPolynomial<double>/12 83.0 ns 82.9 ns
BezierPolynomial<double>/13 108 ns 108 ns
BezierPolynomial<double>/14 119 ns 119 ns
BezierPolynomial<double>/15 140 ns 140 ns
BezierPolynomial<double>/16 137 ns 137 ns
BezierPolynomial<double>/17 151 ns 151 ns
BezierPolynomial<double>/18 171 ns 171 ns
BezierPolynomial<double>/19 194 ns 193 ns
BezierPolynomial<double>/20 213 ns 213 ns
BezierPolynomial<double>/21 220 ns 220 ns
BezierPolynomial<double>/22 260 ns 260 ns
BezierPolynomial<double>/23 266 ns 266 ns
BezierPolynomial<double>/24 293 ns 292 ns
BezierPolynomial<double>/25 319 ns 319 ns
BezierPolynomial<double>/26 336 ns 335 ns
BezierPolynomial<double>/27 370 ns 370 ns
BezierPolynomial<double>/28 429 ns 429 ns
BezierPolynomial<double>/29 443 ns 443 ns
BezierPolynomial<double>/30 421 ns 421 ns
BezierPolynomial<double>_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 Real>
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<Real, 3> 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]

1
doc/math.qbk
View File

@@ -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]

60
include/boost/math/interpolators/bezier_polynomial.hpp Normal file
View File

@@ -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 <memory>
#include <boost/math/interpolators/detail/bezier_polynomial_detail.hpp>
namespace boost::math::interpolators {
template <class RandomAccessContainer>
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<detail::bezier_polynomial_imp<RandomAccessContainer>>(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<RandomAccessContainer> indefinite_integral() const {
return std::move(m_imp->indefinite_integral());
}
friend std::ostream& operator<<(std::ostream& out, bezier_polynomial<RandomAccessContainer> const & bp) {
out << *bp.m_imp;
return out;
}
private:
std::shared_ptr<detail::bezier_polynomial_imp<RandomAccessContainer>> m_imp;
};
}
#endif

164
include/boost/math/interpolators/detail/bezier_polynomial_detail.hpp Normal file
View File

@@ -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 <stdexcept>
#include <iostream>
#include <string>
namespace boost::math::interpolators::detail {
template <class RandomAccessContainer>
static inline RandomAccessContainer& get_bezier_storage()
{
static thread_local RandomAccessContainer the_storage;
return the_storage;
}
template <class RandomAccessContainer>
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<RandomAccessContainer>();
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<Real>::quiet_NaN();
}
return p;
}
auto & scratch_space = get_bezier_storage<RandomAccessContainer>();
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<RandomAccessContainer>();
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<RandomAccessContainer> 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

42
reporting/performance/bezier_polynomial_performance.cpp Normal file
View File

@@ -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 <random>
#include <array>
#include <vector>
#include <benchmark/benchmark.h>
#include <boost/math/interpolators/bezier_polynomial.hpp>
using boost::math::interpolators::bezier_polynomial;
template<class Real>
void BezierPolynomial(benchmark::State& state)
{
std::random_device rd;
std::mt19937_64 mt(rd());
std::uniform_real_distribution<Real> unif(0, 10);
std::vector<std::array<Real, 3>> 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<Real>::epsilon();
}
state.SetComplexityN(state.range(0));
}
BENCHMARK_TEMPLATE(BezierPolynomial, double)->DenseRange(2, 30)->Complexity();
BENCHMARK_MAIN();

1
test/Jamfile.v2
View File

@@ -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" : <linkflags>-lquadmath ] ]
[ run cubic_hermite_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-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" : <linkflags>-lquadmath ] ]
[ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : <define>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 : : : <define>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 : : : <define>TEST=3 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_3 ]

263
test/bezier_polynomial_test.cpp Normal file
View File

@@ -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 <numeric>
#include <random>
#include <array>
#include <boost/core/demangle.hpp>
#include <boost/math/interpolators/bezier_polynomial.hpp>
#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/float128.hpp>
using boost::multiprecision::float128;
#endif
using boost::math::interpolators::bezier_polynomial;
template<typename Real>
void test_linear()
{
std::vector<std::array<Real, 2>> 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<Real, 2> 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<typename Real>
void test_quadratic()
{
std::vector<std::array<Real, 2>> 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<typename Real>
void test_convex_hull()
{
std::vector<std::array<Real, 2>> 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<typename Real>
void test_reversal_symmetry()
{
std::vector<std::array<Real, 3>> control_points(10);
std::uniform_real_distribution<Real> 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<typename Real>
void test_linear_precision()
{
std::vector<std::array<Real, 3>> control_points(10);
std::array<Real, 3> P0 = {1,1,1};
std::array<Real, 3> 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<Real, 3> 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<Real, 3> 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<typename Real>
void test_indefinite_integral()
{
std::vector<std::array<Real, 3>> control_points(2);
std::uniform_real_distribution<Real> 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<Real, 3> 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<float>();
test_linear<double>();
test_quadratic<float>();
test_quadratic<double>();
test_convex_hull<float>();
test_convex_hull<double>();
test_linear_precision<float>();
test_linear_precision<double>();
test_reversal_symmetry<float>();
test_reversal_symmetry<double>();
//test_indefinite_integral<float>();
//test_indefinite_integral<double>();
#ifdef BOOST_HAS_FLOAT128
test_linear<float128>();
test_quadratic<float128>();
test_convex_hull<float128>();
#endif
return boost::math::test::report_errors();
}
#endif