From cad34ff7562abafe8bb30e4950a105772fcd4494 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Thu, 21 Dec 2017 16:12:24 -0700 Subject: [PATCH 01/20] First pass at a Catmull-Rom curve interpolator. --- doc/interpolators/catmull_rom.qbk | 99 ++++++++++ example/catmull_rom_example.cpp | 14 ++ .../boost/math/interpolators/catmull_rom.hpp | 179 ++++++++++++++++++ test/catmull_rom_test.cpp | 73 +++++++ 4 files changed, 365 insertions(+) create mode 100644 doc/interpolators/catmull_rom.qbk create mode 100644 example/catmull_rom_example.cpp create mode 100644 include/boost/math/interpolators/catmull_rom.hpp create mode 100644 test/catmull_rom_test.cpp diff --git a/doc/interpolators/catmull_rom.qbk b/doc/interpolators/catmull_rom.qbk new file mode 100644 index 000000000..00457986f --- /dev/null +++ b/doc/interpolators/catmull_rom.qbk @@ -0,0 +1,99 @@ +[/ + Copyright 2017 Nick Thompson + + 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:catmull_rom Catmull-Rom Splines] + +[heading Synopsis] + +`` +#include + +namespace boost{ namespace math{ + template + class catmull_rom + { + public: + catmull_rom(const Real* const points, size_t num_points, size_t dimension, bool closed = false, Real alpha = (Real) 1/ (Real) 2) + + Real operator()(Real s) const; + + Real max_parameter() const; + + Point prime(Real s) const; + }; + +}} +`` + +[heading Description] + +Catmull-Rom splines are a family of interpolating curves which are commonly used in computer graphics and animation. +Catmull-Rom splines enjoy the following nice properties: + +* Affine invariance: The interpolant commutes with affine transformations. +* Local support of the basis functions: This gives stability and fast evaluation. +* Smoothness +* Interpolation of control points: Many curves (such as Bezier) are *approximating*, they do not pass through their control points. This makes them more difficult to use than interpolating splines. + +The `catmull_rom` class provided by boost creates a cubic Catmull-Rom spline from an array of points in any dimension. +Since there are numerous ways to represent a point in /n/-dimensional space, +the class attempts to be flexible by templating on the point type. +The requirement of the point type is that it must have a dereference operator defined (e.g., `p[0]` is not a syntax error), +it must be able to be dereferenced up to `dimension -1`, and `p[i]` is of type `Real`. +The basic usage is shown here: + + std::vector> points(4); + points[0] = {0,0,0}; + points[1] = {1,0,0}; + points[2] = {0,1,0}; + points[3] = {0,0,1}; + catmull_rom, 3> cr(points.data(), points.size()); + // Interpolate at s = 0.1: + auto point = cr(0.1); + +The spline can be either open or *closed*, meaning that there is some /t/ such /P(t) = P(0)/. +The default is open, but this can be easily changed: + + // closed = true + catmull_rom, 3> cr(points.data(), points.size(), true); + +The Catmull-Rom curve admits an infinite number of parameterizations. +The default parameterization of the `catmull_rom` class is the so-called *centripedal* parameterization. +This parameterization has been shown to be the only parameterization that does not have cusps or self-intersections within segments. +However, for advanced users, other parameterizations can be chosen using the /alpha/ parameter: + + // alpha = 1 is the "chordal" parameterization. + catmull_rom, 3> cr(points.data(), points.size(), false, (Real) 1); + +Since the parameterization is chosen by the routine, +it is often of interest to query the maximum bound `s_max` such that `P(s_max) = P_f` where `P_f` is the final point on the list. +This is easily discovered: + + Real max_s = cr.max_parameter(); + +If you attempt to interpolate for `s > max_s`, an exception is thrown. + +Finally, the tangent vector to any point of the curve can be discovered via + + Real s = 0.1; + Point tangent = cr.prime(s); + +Since the magnitude of the tangent vector is dependent on the parameterization, it is not very meaningful. +However, its direction is meaningful, so the user may wish to normalize this result. + +[heading Examples] + +[import ../../example/catmull_rom_example.cpp] + +[section:catmull_rom_refs References] + +* Cem Yuksel, Scott Schaefer, and John Keyser, ['Parameterization and applications of Catmull–Rom curves], Computer-Aided Design 43 (2011) 747–755. +* Phillip J. Barry and Ronald N. Goldman, ['A Recursive Evaluation Algorithm for a Class of Catmull-Rom Splines], Computer Graphics, Volume 22, Number 4, August 1988 + +[endsect] +[endsect] [/section:catmull_rom Catmull-Rom Splines] diff --git a/example/catmull_rom_example.cpp b/example/catmull_rom_example.cpp new file mode 100644 index 000000000..40d0c5ad1 --- /dev/null +++ b/example/catmull_rom_example.cpp @@ -0,0 +1,14 @@ +// Copyright Nick Thompson, 2017 + +// 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). + +#include +#include + +int main() +{ + std::cout << "This shows how to use the Catmull-Rom class of Boost to create a logarithmic spiral.\n"; + return 0; +} diff --git a/include/boost/math/interpolators/catmull_rom.hpp b/include/boost/math/interpolators/catmull_rom.hpp new file mode 100644 index 000000000..1e0734f37 --- /dev/null +++ b/include/boost/math/interpolators/catmull_rom.hpp @@ -0,0 +1,179 @@ +// Copyright Nick Thompson, 2017 +// 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) + +// This computes the Catmull-Rom spline from a list of points. + +#ifndef BOOST_MATH_INTERPOLATORS_CATMULL_ROM +#define BOOST_MATH_INTERPOLATORS_CATMULL_ROM + +#include +#include + +namespace boost{ namespace math{ + + namespace detail + { + template + Real alpha_distance(Point const & p1, Point const & p2, Real alpha) + { + using std::pow; + Real dsq = 0; + for (size_t i = 0; i < dimension; ++i) + { + Real dx = p1[i] - p2[i]; + dsq += dx*dx; + } + return pow(dsq, alpha/2); + } + } + +template +class catmull_rom +{ +public: + catmull_rom(const Point* const points, size_t num_pnts, bool closed = false, Real alpha = (Real) 1/ (Real) 2) + { + if (num_pnts < 4) + { + throw std::domain_error("The Catmull-Rom curve requires at least 4 points.\n"); + } + using std::abs; + m_s.resize(num_pnts+3); + m_pnts.resize(num_pnts+3); + if (closed) + { + m_pnts[0] = points[num_pnts-1]; + for (size_t i = 0; i < num_pnts; ++i) { + m_pnts[i+1] = points[i]; + } + m_pnts[num_pnts+1] = points[0]; + m_pnts[num_pnts+2] = points[1]; + m_s[0] = -detail::alpha_distance(m_pnts[0], m_pnts[1], alpha); + if (abs(m_s[0]) < std::numeric_limits::epsilon()) + { + throw std::domain_error("The first and last point should not be the same for a closed curve-that is indicated by a boolean flag.\n"); + } + m_s[1] = 0; + for (size_t i = 2; i < m_s.size(); ++i) + { + Real d = detail::alpha_distance(m_pnts[i], m_pnts[i-1], alpha); + if (abs(d) < std::numeric_limits::epsilon()) + { + throw std::domain_error("The control points of the Catmull-Rom curve are too close together; this will lead to ill-conditioning.\n"); + } + m_s[i] = m_s[i-1] + d; + } + m_max_s = m_s[num_pnts]; + } + else + { + for(size_t i = 0; i < dimension; ++i) + { + m_pnts[0][i] = 0; + m_pnts[num_pnts+1][i] = 0; + m_pnts[num_pnts+2][i] = 0; + } + + for (size_t i = 0; i < num_pnts; ++i) + { + m_pnts[i+1] = points[i]; + } + m_pnts[num_pnts+1] = points[0]; + m_pnts[num_pnts+2] = points[1]; + m_s[0] = -1; + m_s[1] = 0; + for (size_t i = 2; i < num_pnts+1; ++i) + { + Real d = detail::alpha_distance(m_pnts[i-1], m_pnts[i-2], alpha); + if (abs(d) < std::numeric_limits::epsilon()) + { + throw std::domain_error("The control points of the Catmull-Rom curve are too close together; this will lead to ill-conditioning.\n"); + } + m_s[i] = m_s[i-1] + d; + } + m_max_s = m_s[num_pnts]; + m_s[num_pnts+2] = m_s[num_pnts+1] + 1; + } + } + + Real max_parameter() const + { + return m_max_s; + } + + Point operator()(Real s) const + { + if (s < 0 || s > m_max_s) + { + throw std::domain_error("Parameter outside bounds.\n"); + } + auto it = std::lower_bound(m_s.begin(), m_s.end(), s); + size_t i = std::distance(m_s.begin(), it); + //std::cout << "Index = " << i << std::endl; + assert(m_s[i] <= s && s < m_s[i+1]); + Point A1; + Real denom = 1/(m_s[i] - m_s[i-1]); + Real k1 = (m_s[i]-s)*denom; + Real k2 = (s - m_s[i-1])*denom; + for (size_t j = 0; j < dimension; ++j) + { + A1[j] = k1*m_pnts[i-1][j] + k2*m_pnts[i][j]; + } + + Point A2; + denom = 1/(m_s[i+1] - m_s[i]); + k1 = (m_s[i+1]-s)*denom; + k2 = (s - m_s[i])*denom; + for (size_t j = 0; j < dimension; ++j) + { + A2[j] = k1*m_pnts[i][j] + k2*m_pnts[i+1][j]; + } + + Point B1; + for (size_t j = 0; j < dimension; ++j) + { + B1[j] = k1*A1[j] + k2*A2[j]; + } + + Point A3; + denom = 1/(m_s[i+2] - m_s[i+1]); + k1 = (m_s[i+2]-s)*denom; + k2 = (s - m_s[i+1])*denom; + for (size_t j = 0; j < dimension; ++j) + { + A3[j] = k1*m_pnts[i+1][j] + k2*m_pnts[i+2][j]; + } + + Point B2; + denom = 1/(m_s[i+2] - m_s[i]); + k1 = (m_s[i+2]-s)*denom; + k2 = (s - m_s[i])*denom; + for (size_t j = 0; j < dimension; ++j) + { + B2[j] = k1*A2[j] + k2*A3[j]; + } + + Point C; + denom = 1/(m_s[i+1] - m_s[i]); + k1 = (m_s[i+1]-s)*denom; + k2 = (s - m_s[i])*denom; + for (size_t j = 0; j < dimension; ++j) + { + C[j] = k1*B1[j] + k2*B2[j]; + } + + + return C; + } + +private: + std::vector m_pnts; + std::vector m_s; + Real m_max_s; +}; + +}} +#endif diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp new file mode 100644 index 000000000..56f0e5db3 --- /dev/null +++ b/test/catmull_rom_test.cpp @@ -0,0 +1,73 @@ +/* + * Copyright Nick Thompson, 2017 + * 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) + */ +#define BOOST_TEST_MODULE catmull_rom_test + +#include +#include +#include +#include +#include +#include +#include +#include + + +using boost::math::catmull_rom; + +template +void test_linear() +{ + Real tol = std::numeric_limits::epsilon(); + std::vector> v(4); + v[0] = {0,0,0}; + v[1] = {1,0,0}; + v[2] = {2,0,0}; + v[3] = {3,0,0}; + catmull_rom, 3> cr(v.data(), v.size(), true); + + BOOST_CHECK_CLOSE_FRACTION(cr.max_parameter(), 3, tol); + auto p0 = cr(0.0); + BOOST_CHECK_SMALL(p0[0], tol); + BOOST_CHECK_SMALL(p0[1], tol); + BOOST_CHECK_SMALL(p0[2], tol); + auto p1 = cr(1.0); + BOOST_CHECK_CLOSE_FRACTION(p1[0], 1, tol); + BOOST_CHECK_SMALL(p1[1], tol); + BOOST_CHECK_SMALL(p1[2], tol); + + auto p2 = cr(2.0); + BOOST_CHECK_CLOSE_FRACTION(p2[0], 2, tol); + BOOST_CHECK_SMALL(p2[1], tol); + BOOST_CHECK_SMALL(p2[2], tol); + + + auto p3 = cr(3.0); + BOOST_CHECK_CLOSE_FRACTION(p3[0], 3, tol); + BOOST_CHECK_SMALL(p3[1], tol); + BOOST_CHECK_SMALL(p3[2], tol); + +} + +template +void test_affine_invariance() +{ + std::cout << "Testing that the Catmull-Rom spline is affine invariant.\n"; + +} + +template +void test_data_representations() +{ + std::cout << "Testing that the Catmull-Rom spline works with multiple data representations.\n"; +} + +BOOST_AUTO_TEST_CASE(catmull_rom_test) +{ + test_linear(); + test_affine_invariance(); + test_data_representations(); +} From 4cce07e15e85e6822cd044694e223f1b9e8a5323 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Fri, 22 Dec 2017 16:31:36 -0700 Subject: [PATCH 02/20] Implement tangent vector computation. Fix index lookup. Close all Catmull-Rom curves and document. --- doc/interpolators/catmull_rom.qbk | 33 +- example/catmull_rom_example.cpp | 33 +- .../boost/math/interpolators/catmull_rom.hpp | 344 +++++++++++------- test/catmull_rom_test.cpp | 74 +++- 4 files changed, 346 insertions(+), 138 deletions(-) diff --git a/doc/interpolators/catmull_rom.qbk b/doc/interpolators/catmull_rom.qbk index 00457986f..d93490d0e 100644 --- a/doc/interpolators/catmull_rom.qbk +++ b/doc/interpolators/catmull_rom.qbk @@ -62,6 +62,29 @@ The default is open, but this can be easily changed: // closed = true catmull_rom, 3> cr(points.data(), points.size(), true); +Inside `catmull_rom`, the Catmull-Rom curve is represented as closed. +This is because an open Catmull-Rom curve is *implicitly closed*, but the closing point is the zero vector. +There is no reason to suppose that the zero vector is a better closing point than the endpoint (or any other point, for that matter), +so traditionally Catmull-Rom splines leave the segment between the first and second point undefined, +as well as the segment between the second-to-last and last point. +We find this property of the traditional implementation of Catmull-Rom splines annoying and confusing to the user. +Hence internally, we close the curve so that the first and last segments are defined. +Of course, this causes the *tangent* vectors to the first and last points to be bizarre. +This is a "pick your poison" design decision-either the curve cannot interpolate in its first and last segments, +or the tangents along the first and last segments are meaningless. + +Since the routine internally represents the curve as closed, +a question arises: Why does the user have to specify if the curve is open or closed? +The answer is that the parameterization is chosen by the routine, +so it is of interest to the user to understand the values where a meaningful result is returned. + + Real max_s = cr.max_parameter(); + +If you attempt to interpolate for `s > max_s`, an exception is thrown. +If the curve is closed, then `cr(max_s) = p0`, where `p0` is the first point on the curve. +If the curve is open, then `cr(max_s) = pf`, where `pf` is the final point on the curve. + + The Catmull-Rom curve admits an infinite number of parameterizations. The default parameterization of the `catmull_rom` class is the so-called *centripedal* parameterization. This parameterization has been shown to be the only parameterization that does not have cusps or self-intersections within segments. @@ -70,20 +93,14 @@ However, for advanced users, other parameterizations can be chosen using the /al // alpha = 1 is the "chordal" parameterization. catmull_rom, 3> cr(points.data(), points.size(), false, (Real) 1); -Since the parameterization is chosen by the routine, -it is often of interest to query the maximum bound `s_max` such that `P(s_max) = P_f` where `P_f` is the final point on the list. -This is easily discovered: - - Real max_s = cr.max_parameter(); - -If you attempt to interpolate for `s > max_s`, an exception is thrown. Finally, the tangent vector to any point of the curve can be discovered via Real s = 0.1; Point tangent = cr.prime(s); -Since the magnitude of the tangent vector is dependent on the parameterization, it is not very meaningful. +Since the magnitude of the tangent vector is dependent on the parameterization, +it is not as meaningful as (say) arc-length parameterization. However, its direction is meaningful, so the user may wish to normalize this result. [heading Examples] diff --git a/example/catmull_rom_example.cpp b/example/catmull_rom_example.cpp index 40d0c5ad1..d124074a8 100644 --- a/example/catmull_rom_example.cpp +++ b/example/catmull_rom_example.cpp @@ -5,10 +5,41 @@ // copy at http://www.boost.org/LICENSE_1_0.txt). #include +#include +#include +#include #include +using boost::math::catmull_rom; + int main() { - std::cout << "This shows how to use the Catmull-Rom class of Boost to create a logarithmic spiral.\n"; + std::cout << "This shows how to use Boost's Catmull-Rom spline to create an Archimedean spiral.\n"; + + // The Archimedean spiral is given by r = a*theta. We have set a = 1. + std::vector> spiral_points(500); + double theta_max = M_PI; + for (size_t i = 0; i < spiral_points.size(); ++i) + { + double theta = i/theta_max; + spiral_points[i] = {theta*cos(theta), theta*sin(theta)}; + } + + auto archimedean = catmull_rom, 2>(spiral_points.data(), spiral_points.size()); + double max_s = archimedean.max_parameter(); + for (double s = 0; s < max_s; s += 0.1) + { + auto p = archimedean(s); + double x = p[0]; + double y = p[1]; + double r = sqrt(x*x + y*y); + double theta = atan2(y/r, x/r); + if (theta < 0) + { + theta += 2*M_PI; + } + std::cout << "r = " << r << ", theta = " << theta << std::endl; + } + return 0; } diff --git a/include/boost/math/interpolators/catmull_rom.hpp b/include/boost/math/interpolators/catmull_rom.hpp index 1e0734f37..e45601453 100644 --- a/include/boost/math/interpolators/catmull_rom.hpp +++ b/include/boost/math/interpolators/catmull_rom.hpp @@ -10,6 +10,7 @@ #define BOOST_MATH_INTERPOLATORS_CATMULL_ROM #include +#include #include namespace boost{ namespace math{ @@ -34,146 +35,239 @@ template class catmull_rom { public: - catmull_rom(const Point* const points, size_t num_pnts, bool closed = false, Real alpha = (Real) 1/ (Real) 2) - { - if (num_pnts < 4) - { - throw std::domain_error("The Catmull-Rom curve requires at least 4 points.\n"); - } - using std::abs; - m_s.resize(num_pnts+3); - m_pnts.resize(num_pnts+3); - if (closed) - { - m_pnts[0] = points[num_pnts-1]; - for (size_t i = 0; i < num_pnts; ++i) { - m_pnts[i+1] = points[i]; - } - m_pnts[num_pnts+1] = points[0]; - m_pnts[num_pnts+2] = points[1]; - m_s[0] = -detail::alpha_distance(m_pnts[0], m_pnts[1], alpha); - if (abs(m_s[0]) < std::numeric_limits::epsilon()) - { - throw std::domain_error("The first and last point should not be the same for a closed curve-that is indicated by a boolean flag.\n"); - } - m_s[1] = 0; - for (size_t i = 2; i < m_s.size(); ++i) - { - Real d = detail::alpha_distance(m_pnts[i], m_pnts[i-1], alpha); - if (abs(d) < std::numeric_limits::epsilon()) - { - throw std::domain_error("The control points of the Catmull-Rom curve are too close together; this will lead to ill-conditioning.\n"); - } - m_s[i] = m_s[i-1] + d; - } - m_max_s = m_s[num_pnts]; - } - else - { - for(size_t i = 0; i < dimension; ++i) - { - m_pnts[0][i] = 0; - m_pnts[num_pnts+1][i] = 0; - m_pnts[num_pnts+2][i] = 0; - } - for (size_t i = 0; i < num_pnts; ++i) - { - m_pnts[i+1] = points[i]; - } - m_pnts[num_pnts+1] = points[0]; - m_pnts[num_pnts+2] = points[1]; - m_s[0] = -1; - m_s[1] = 0; - for (size_t i = 2; i < num_pnts+1; ++i) - { - Real d = detail::alpha_distance(m_pnts[i-1], m_pnts[i-2], alpha); - if (abs(d) < std::numeric_limits::epsilon()) - { - throw std::domain_error("The control points of the Catmull-Rom curve are too close together; this will lead to ill-conditioning.\n"); - } - m_s[i] = m_s[i-1] + d; - } - m_max_s = m_s[num_pnts]; - m_s[num_pnts+2] = m_s[num_pnts+1] + 1; - } - } + catmull_rom(const Point* const points, size_t num_pnts, bool closed = false, Real alpha = (Real) 1/ (Real) 2); Real max_parameter() const { return m_max_s; } - Point operator()(Real s) const + Real parameter_at_point(size_t i) const { - if (s < 0 || s > m_max_s) - { - throw std::domain_error("Parameter outside bounds.\n"); - } - auto it = std::lower_bound(m_s.begin(), m_s.end(), s); - size_t i = std::distance(m_s.begin(), it); - //std::cout << "Index = " << i << std::endl; - assert(m_s[i] <= s && s < m_s[i+1]); - Point A1; - Real denom = 1/(m_s[i] - m_s[i-1]); - Real k1 = (m_s[i]-s)*denom; - Real k2 = (s - m_s[i-1])*denom; - for (size_t j = 0; j < dimension; ++j) - { - A1[j] = k1*m_pnts[i-1][j] + k2*m_pnts[i][j]; - } - - Point A2; - denom = 1/(m_s[i+1] - m_s[i]); - k1 = (m_s[i+1]-s)*denom; - k2 = (s - m_s[i])*denom; - for (size_t j = 0; j < dimension; ++j) - { - A2[j] = k1*m_pnts[i][j] + k2*m_pnts[i+1][j]; - } - - Point B1; - for (size_t j = 0; j < dimension; ++j) - { - B1[j] = k1*A1[j] + k2*A2[j]; - } - - Point A3; - denom = 1/(m_s[i+2] - m_s[i+1]); - k1 = (m_s[i+2]-s)*denom; - k2 = (s - m_s[i+1])*denom; - for (size_t j = 0; j < dimension; ++j) - { - A3[j] = k1*m_pnts[i+1][j] + k2*m_pnts[i+2][j]; - } - - Point B2; - denom = 1/(m_s[i+2] - m_s[i]); - k1 = (m_s[i+2]-s)*denom; - k2 = (s - m_s[i])*denom; - for (size_t j = 0; j < dimension; ++j) - { - B2[j] = k1*A2[j] + k2*A3[j]; - } - - Point C; - denom = 1/(m_s[i+1] - m_s[i]); - k1 = (m_s[i+1]-s)*denom; - k2 = (s - m_s[i])*denom; - for (size_t j = 0; j < dimension; ++j) - { - C[j] = k1*B1[j] + k2*B2[j]; - } - - - return C; + return m_s[i+1]; } + Point operator()(Real s) const; + + + Point prime(Real s) const; + private: std::vector m_pnts; std::vector m_s; Real m_max_s; }; +template +catmull_rom::catmull_rom(const Point* const points, size_t num_pnts, bool closed, Real alpha) +{ + if (num_pnts < 4) + { + throw std::domain_error("The Catmull-Rom curve requires at least 4 points.\n"); + } + using std::abs; + m_s.resize(num_pnts+3); + m_pnts.resize(num_pnts+3); + + m_pnts[0] = points[num_pnts-1]; + for (size_t i = 0; i < num_pnts; ++i) + { + m_pnts[i+1] = points[i]; + } + m_pnts[num_pnts+1] = points[0]; + m_pnts[num_pnts+2] = points[1]; + m_s[0] = -detail::alpha_distance(m_pnts[0], m_pnts[1], alpha); + if (abs(m_s[0]) < std::numeric_limits::epsilon()) + { + throw std::domain_error("The first and last point should not be the same.\n"); + } + m_s[1] = 0; + for (size_t i = 2; i < m_s.size(); ++i) + { + Real d = detail::alpha_distance(m_pnts[i], m_pnts[i-1], alpha); + if (abs(d) < std::numeric_limits::epsilon()) + { + throw std::domain_error("The control points of the Catmull-Rom curve are too close together; this will lead to ill-conditioning.\n"); + } + m_s[i] = m_s[i-1] + d; + } + if(closed) + { + m_max_s = m_s[num_pnts+1]; + } + else + { + m_max_s = m_s[num_pnts]; + } +} + + +template +Point catmull_rom::operator()(Real s) const +{ + if (s < 0 || s > m_max_s) + { + throw std::domain_error("Parameter outside bounds.\n"); + } + auto it = std::upper_bound(m_s.begin(), m_s.end(), s); + //Now *it >= s. We want the index such that m_s[i] <= s < m_s[i+1]: + size_t i = std::distance(m_s.begin(), it - 1); + // We'll keep the assert in here a while until we feel good that we've understood this algorithm. + assert(m_s[i] <= s && s < m_s[i+1]); + Point A1; + Real denom = 1/(m_s[i] - m_s[i-1]); + Real k1 = (m_s[i]-s)*denom; + Real k2 = (s - m_s[i-1])*denom; + for (size_t j = 0; j < dimension; ++j) + { + A1[j] = k1*m_pnts[i-1][j] + k2*m_pnts[i][j]; + } + + Point A2; + denom = 1/(m_s[i+1] - m_s[i]); + k1 = (m_s[i+1]-s)*denom; + k2 = (s - m_s[i])*denom; + for (size_t j = 0; j < dimension; ++j) + { + A2[j] = k1*m_pnts[i][j] + k2*m_pnts[i+1][j]; + } + + Point B1; + for (size_t j = 0; j < dimension; ++j) + { + B1[j] = k1*A1[j] + k2*A2[j]; + } + + Point A3; + denom = 1/(m_s[i+2] - m_s[i+1]); + k1 = (m_s[i+2]-s)*denom; + k2 = (s - m_s[i+1])*denom; + for (size_t j = 0; j < dimension; ++j) + { + A3[j] = k1*m_pnts[i+1][j] + k2*m_pnts[i+2][j]; + } + + Point B2; + denom = 1/(m_s[i+2] - m_s[i]); + k1 = (m_s[i+2]-s)*denom; + k2 = (s - m_s[i])*denom; + for (size_t j = 0; j < dimension; ++j) + { + B2[j] = k1*A2[j] + k2*A3[j]; + } + + Point C; + denom = 1/(m_s[i+1] - m_s[i]); + k1 = (m_s[i+1]-s)*denom; + k2 = (s - m_s[i])*denom; + for (size_t j = 0; j < dimension; ++j) + { + C[j] = k1*B1[j] + k2*B2[j]; + } + + + return C; +} + +template +Point catmull_rom::prime(Real s) const +{ + // https://math.stackexchange.com/questions/843595/how-can-i-calculate-the-derivative-of-a-catmull-rom-spline-with-nonuniform-param + // http://denkovacs.com/2016/02/catmull-rom-spline-derivatives/ + if (s < 0 || s > m_max_s) + { + throw std::domain_error("Parameter outside bounds.\n"); + } + auto it = std::upper_bound(m_s.begin(), m_s.end(), s); + //Now *it >= s. We want the index such that m_s[i] <= s < m_s[i+1]: + size_t i = std::distance(m_s.begin(), it - 1); + // We'll keep the assert in here a while until we feel good that we've understood this algorithm. + assert(m_s[i] <= s && s < m_s[i+1]); + Point A1; + Real denom = 1/(m_s[i] - m_s[i-1]); + Real k1 = (m_s[i]-s)*denom; + Real k2 = (s - m_s[i-1])*denom; + for (size_t j = 0; j < dimension; ++j) + { + A1[j] = k1*m_pnts[i-1][j] + k2*m_pnts[i][j]; + } + + Point A1p; + for (size_t j = 0; j < dimension; ++j) + { + A1p[j] = denom*(m_pnts[i][j] - m_pnts[i-1][j]); + } + + Point A2; + denom = 1/(m_s[i+1] - m_s[i]); + k1 = (m_s[i+1]-s)*denom; + k2 = (s - m_s[i])*denom; + for (size_t j = 0; j < dimension; ++j) + { + A2[j] = k1*m_pnts[i][j] + k2*m_pnts[i+1][j]; + } + + Point A2p; + for (size_t j = 0; j < dimension; ++j) + { + A2p[j] = denom*(m_pnts[i+1][j] - m_pnts[i][j]); + } + + + Point B1; + for (size_t j = 0; j < dimension; ++j) + { + B1[j] = k1*A1[j] + k2*A2[j]; + } + + Point A3; + denom = 1/(m_s[i+2] - m_s[i+1]); + k1 = (m_s[i+2]-s)*denom; + k2 = (s - m_s[i+1])*denom; + for (size_t j = 0; j < dimension; ++j) + { + A3[j] = k1*m_pnts[i+1][j] + k2*m_pnts[i+2][j]; + } + + Point A3p; + for (size_t j = 0; j < dimension; ++j) + { + A3p[j] = denom*(m_pnts[i+2][j] - m_pnts[i+1][j]); + } + + Point B2; + denom = 1/(m_s[i+2] - m_s[i]); + k1 = (m_s[i+2]-s)*denom; + k2 = (s - m_s[i])*denom; + for (size_t j = 0; j < dimension; ++j) + { + B2[j] = k1*A2[j] + k2*A3[j]; + } + + Point B1p; + denom = 1/(m_s[i+1] - m_s[i-1]); + for (size_t j = 0; j < dimension; ++j) + { + B1p[j] = denom*(A2[j] - A1[j] + (m_s[i+1]- s)*A1p[j] + (s-m_s[i-1])*A2p[j]); + } + + Point B2p; + denom = 1/(m_s[i+2] - m_s[i]); + for (size_t j = 0; j < dimension; ++j) + { + B2p[j] = denom*(A3[j] - A2[j] + (m_s[i+2] - s)*A2p[j] + (s - m_s[i])*A3p[j]); + } + + Point Cp; + denom = 1/(m_s[i+1] - m_s[i]); + for (size_t j = 0; j < dimension; ++j) + { + Cp[j] = denom*(B2[j] - B1[j] + (m_s[i+1] - s)*B1p[j] + (s - m_s[i])*B2p[j]); + } + return Cp; +} + + }} #endif diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index 56f0e5db3..b616bd46a 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -15,20 +15,50 @@ #include #include - +using boost::multiprecision::cpp_bin_float_50; using boost::math::catmull_rom; +template +void test_alpha_distance() +{ + Real tol = std::numeric_limits::epsilon(); + std::array v1 = {0,0,0}; + std::array v2 = {1,0,0}; + Real alpha = 0.5; + Real d = boost::math::detail::alpha_distance, 3>(v1, v2, alpha); + BOOST_CHECK_CLOSE_FRACTION(d, 1, tol); + + d = boost::math::detail::alpha_distance, 3>(v1, v2, 0); + BOOST_CHECK_CLOSE_FRACTION(d, 1, tol); + + d = boost::math::detail::alpha_distance, 3>(v1, v2, 1); + BOOST_CHECK_CLOSE_FRACTION(d, 1, tol); + + + v2[0] = 2; + d = boost::math::detail::alpha_distance, 3>(v1, v2, alpha); + BOOST_CHECK_CLOSE_FRACTION(d, pow(2, (Real)1/ (Real) 2), tol); + + d = boost::math::detail::alpha_distance, 3>(v1, v2, 0); + BOOST_CHECK_CLOSE_FRACTION(d, 1, tol); + + d = boost::math::detail::alpha_distance, 3>(v1, v2, 1); + BOOST_CHECK_CLOSE_FRACTION(d, 2, tol); +} + + template void test_linear() { - Real tol = std::numeric_limits::epsilon(); + Real tol = 10*std::numeric_limits::epsilon(); std::vector> v(4); v[0] = {0,0,0}; v[1] = {1,0,0}; v[2] = {2,0,0}; v[3] = {3,0,0}; - catmull_rom, 3> cr(v.data(), v.size(), true); + catmull_rom, 3> cr(v.data(), v.size()); + // Test that the interpolation condition is obeyed: BOOST_CHECK_CLOSE_FRACTION(cr.max_parameter(), 3, tol); auto p0 = cr(0.0); BOOST_CHECK_SMALL(p0[0], tol); @@ -50,6 +80,33 @@ void test_linear() BOOST_CHECK_SMALL(p3[1], tol); BOOST_CHECK_SMALL(p3[2], tol); + Real s = cr.parameter_at_point(0); + BOOST_CHECK_SMALL(s, tol); + + s = cr.parameter_at_point(1); + BOOST_CHECK_CLOSE_FRACTION(s, 1, tol); + + s = cr.parameter_at_point(2); + BOOST_CHECK_CLOSE_FRACTION(s, 2, tol); + + s = cr.parameter_at_point(3); + BOOST_CHECK_CLOSE_FRACTION(s, 3, tol); + + // Test that the function is linear on the interval [1,2]: + for (double s = 1; s < 2; s += 0.01) + { + auto p = cr(s); + BOOST_CHECK_CLOSE_FRACTION(p[0], s, tol); + BOOST_CHECK_SMALL(p[1], tol); + BOOST_CHECK_SMALL(p[2], tol); + + auto tangent = cr.prime(s); + // TODO: Fix the tangent vector! + //BOOST_CHECK_CLOSE_FRACTION(p[0], 1, tol); + BOOST_CHECK_SMALL(p[1], tol); + BOOST_CHECK_SMALL(p[2], tol); + } + } template @@ -67,7 +124,16 @@ void test_data_representations() BOOST_AUTO_TEST_CASE(catmull_rom_test) { + test_alpha_distance(); + test_alpha_distance(); + test_alpha_distance(); + test_alpha_distance(); + test_linear(); + test_linear(); + test_linear(); + test_linear(); + /* test_affine_invariance(); - test_data_representations(); + test_data_representations();*/ } From 1a29ce99c638248dc8f33eee906f29ac4d367470 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Fri, 22 Dec 2017 21:52:52 -0700 Subject: [PATCH 03/20] Demonstrate affine invariance of Catmull-Rom splines. --- .../boost/math/interpolators/catmull_rom.hpp | 1 + test/catmull_rom_test.cpp | 91 +++++++++++++++++-- 2 files changed, 82 insertions(+), 10 deletions(-) diff --git a/include/boost/math/interpolators/catmull_rom.hpp b/include/boost/math/interpolators/catmull_rom.hpp index e45601453..d827aef5a 100644 --- a/include/boost/math/interpolators/catmull_rom.hpp +++ b/include/boost/math/interpolators/catmull_rom.hpp @@ -62,6 +62,7 @@ private: template catmull_rom::catmull_rom(const Point* const points, size_t num_pnts, bool closed, Real alpha) { + static_assert(dimension > 0, "The dimension of the Catmull-Rom spline must be > 0\n"); if (num_pnts < 4) { throw std::domain_error("The Catmull-Rom curve requires at least 4 points.\n"); diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index b616bd46a..554b42c29 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -7,6 +7,7 @@ #define BOOST_TEST_MODULE catmull_rom_test #include +#include #include #include #include @@ -15,6 +16,7 @@ #include #include +using std::abs; using boost::multiprecision::cpp_bin_float_50; using boost::math::catmull_rom; @@ -34,7 +36,6 @@ void test_alpha_distance() d = boost::math::detail::alpha_distance, 3>(v1, v2, 1); BOOST_CHECK_CLOSE_FRACTION(d, 1, tol); - v2[0] = 2; d = boost::math::detail::alpha_distance, 3>(v1, v2, alpha); BOOST_CHECK_CLOSE_FRACTION(d, pow(2, (Real)1/ (Real) 2), tol); @@ -101,19 +102,73 @@ void test_linear() BOOST_CHECK_SMALL(p[2], tol); auto tangent = cr.prime(s); - // TODO: Fix the tangent vector! - //BOOST_CHECK_CLOSE_FRACTION(p[0], 1, tol); - BOOST_CHECK_SMALL(p[1], tol); - BOOST_CHECK_SMALL(p[2], tol); + BOOST_CHECK_CLOSE_FRACTION(tangent[0], 1, tol); + BOOST_CHECK_SMALL(tangent[1], tol); + BOOST_CHECK_SMALL(tangent[2], tol); } } -template +template void test_affine_invariance() { - std::cout << "Testing that the Catmull-Rom spline is affine invariant.\n"; + std::cout << "Testing that the Catmull-Rom spline is affine invariant in dimension " + << dimension << " on type " + << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = 1000*std::numeric_limits::epsilon(); + std::vector> v(100); + std::random_device rd; + std::mt19937_64 gen(rd()); + Real inv_denom = (Real) 100/( (Real) gen.max() + (Real) 2); + for(size_t j = 0; j < dimension; ++j) + { + v[0][j] = gen()*inv_denom; + } + + for (size_t i = 1; i < v.size(); ++i) + { + for(size_t j = 0; j < dimension; ++j) + { + v[i][j] = v[i-1][j] + gen()*inv_denom; + } + } + std::array affine_shift; + for (size_t j = 0; j < dimension; ++j) + { + affine_shift[j] = gen()*inv_denom; + } + + catmull_rom, dimension> cr1(v.data(), v.size()); + + for(size_t i = 0; i< v.size(); ++i) + { + for(size_t j = 0; j < dimension; ++j) + { + v[i][j] += affine_shift[j]; + } + } + + catmull_rom, dimension> cr2(v.data(), v.size()); + + BOOST_CHECK_CLOSE_FRACTION(cr1.max_parameter(), cr2.max_parameter(), tol); + + Real ds = cr1.max_parameter()/1024; + for (Real s = 0; s < cr1.max_parameter(); s += ds) + { + auto p0 = cr1(s); + auto p1 = cr2(s); + auto tangent0 = cr1.prime(s); + auto tangent1 = cr2.prime(s); + for (size_t j = 0; j < dimension; ++j) + { + BOOST_CHECK_CLOSE_FRACTION(p0[j] + affine_shift[j], p1[j], tol); + if (abs(tangent0[j]) > tol) + { + BOOST_CHECK_CLOSE_FRACTION(tangent0[j], tangent1[j], 20*tol); + } + } + } } template @@ -133,7 +188,23 @@ BOOST_AUTO_TEST_CASE(catmull_rom_test) test_linear(); test_linear(); test_linear(); - /* - test_affine_invariance(); - test_data_representations();*/ + + test_affine_invariance(); + test_affine_invariance(); + test_affine_invariance(); + test_affine_invariance(); + + test_affine_invariance(); + test_affine_invariance(); + test_affine_invariance(); + test_affine_invariance(); + + test_affine_invariance(); + test_affine_invariance(); + test_affine_invariance(); + test_affine_invariance(); + test_affine_invariance(); + test_affine_invariance(); + test_affine_invariance(); + test_affine_invariance(); } From d25677a7043c2f5765d094d7c3bb1c3bb587da13 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Fri, 22 Dec 2017 23:07:51 -0700 Subject: [PATCH 04/20] Fix Archmedean spiral example for Catmull-Rom Curve. --- example/catmull_rom_example.cpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/example/catmull_rom_example.cpp b/example/catmull_rom_example.cpp index d124074a8..acdea3360 100644 --- a/example/catmull_rom_example.cpp +++ b/example/catmull_rom_example.cpp @@ -10,6 +10,8 @@ #include #include +using std::sin; +using std::cos; using boost::math::catmull_rom; int main() @@ -21,24 +23,21 @@ int main() double theta_max = M_PI; for (size_t i = 0; i < spiral_points.size(); ++i) { - double theta = i/theta_max; + double theta = ((double) i/ (double) spiral_points.size())*theta_max; spiral_points[i] = {theta*cos(theta), theta*sin(theta)}; } auto archimedean = catmull_rom, 2>(spiral_points.data(), spiral_points.size()); double max_s = archimedean.max_parameter(); - for (double s = 0; s < max_s; s += 0.1) + std::cout << "Max s = " << max_s << std::endl; + for (double s = 0; s < max_s; s += 0.01) { auto p = archimedean(s); double x = p[0]; double y = p[1]; double r = sqrt(x*x + y*y); double theta = atan2(y/r, x/r); - if (theta < 0) - { - theta += 2*M_PI; - } - std::cout << "r = " << r << ", theta = " << theta << std::endl; + std::cout << "r = " << r << ", theta = " << theta << ", r - theta = " << r - theta << std::endl; } return 0; From 9dce6d8bb95210ccda1e2325e30ac52d51244dac Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Fri, 22 Dec 2017 23:30:31 -0700 Subject: [PATCH 05/20] Catmull-Rom: Add example and test to build, update docs. --- doc/interpolators/catmull_rom.qbk | 6 ++++-- doc/math.qbk | 1 + example/Jamfile.v2 | 1 + test/Jamfile.v2 | 1 + 4 files changed, 7 insertions(+), 2 deletions(-) diff --git a/doc/interpolators/catmull_rom.qbk b/doc/interpolators/catmull_rom.qbk index d93490d0e..561473596 100644 --- a/doc/interpolators/catmull_rom.qbk +++ b/doc/interpolators/catmull_rom.qbk @@ -14,16 +14,18 @@ #include namespace boost{ namespace math{ - template + template class catmull_rom { public: - catmull_rom(const Real* const points, size_t num_points, size_t dimension, bool closed = false, Real alpha = (Real) 1/ (Real) 2) + catmull_rom(const Real* const points, size_t num_points, bool closed = false, Real alpha = (Real) 1/ (Real) 2) Real operator()(Real s) const; Real max_parameter() const; + Real parameter_at_point(size_t i) const; + Point prime(Real s) const; }; diff --git a/doc/math.qbk b/doc/math.qbk index 95c7e73fe..75bb0679c 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -629,6 +629,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [mathpart interpolation Interpolation] [include interpolators/cubic_b_spline.qbk] [include interpolators/barycentric_rational_interpolation.qbk] +[include interpolators/catmull_rom.qbk] [endmathpart] [mathpart quadrature Quadrature] diff --git a/example/Jamfile.v2 b/example/Jamfile.v2 index 232bdbce5..5a9c0b29a 100644 --- a/example/Jamfile.v2 +++ b/example/Jamfile.v2 @@ -129,5 +129,6 @@ test-suite examples : [ run barycentric_interpolation_example.cpp : : : [ requires cxx11_smart_ptr ] ] [ run barycentric_interpolation_example_2.cpp : : : [ requires cxx11_smart_ptr cxx11_function_template_default_args ] ] [ run cubic_b_spline_example.cpp : : : [ requires cxx11_smart_ptr cxx11_hdr_random cxx11_defaulted_functions ] ] + [ run catmull_rom_example.cpp ] ; diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index e5e5f60d9..0cf9345ab 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -824,6 +824,7 @@ test-suite misc : [ run test_constant_generate.cpp : : : release USE_CPP_FLOAT=1 off:no ] [ run test_classify.cpp pch ../../test/build//boost_unit_test_framework ] [ run test_cubic_b_spline.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_smart_ptr cxx11_defaulted_functions ] ] + [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework ] [ run test_error_handling.cpp ../../test/build//boost_unit_test_framework ] [ run legendre_stieltjes_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : "-Bstatic -lquadmath -Bdynamic" ] ] [ run test_minima.cpp pch ../../test/build//boost_unit_test_framework ] From 0f007c23c6b65adc673ee3ce9c40c65b612f5e92 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sat, 23 Dec 2017 17:06:20 -0700 Subject: [PATCH 06/20] Reuse temporaries to increase performance. --- doc/interpolators/catmull_rom.qbk | 37 ++++++++++++ .../boost/math/interpolators/catmull_rom.hpp | 58 +++++++++---------- test/catmull_rom_test.cpp | 3 + 3 files changed, 66 insertions(+), 32 deletions(-) diff --git a/doc/interpolators/catmull_rom.qbk b/doc/interpolators/catmull_rom.qbk index 561473596..06179a2cb 100644 --- a/doc/interpolators/catmull_rom.qbk +++ b/doc/interpolators/catmull_rom.qbk @@ -109,6 +109,43 @@ However, its direction is meaningful, so the user may wish to normalize this res [import ../../example/catmull_rom_example.cpp] +[heading Performance] + +The following performance numbers were generated for a call to the Catmull-Rom interpolation method. +The number that follows the slash is the number of points passed to the interpolant. + + +Run on 2700 MHz CPU +CPU Caches: + L1 Data 32K (x2) + L1 Instruction 32K (x2) + L2 Unified 262K (x2) + L3 Unified 3145K (x1) +--------------------------------------------------------- +Benchmark Time CPU +--------------------------------------------------------- +BM_CatmullRom/4 20 ns 20 ns +BM_CatmullRom/8 21 ns 21 ns +BM_CatmullRom/16 23 ns 23 ns +BM_CatmullRom/32 24 ns 24 ns +BM_CatmullRom/64 27 ns 27 ns +BM_CatmullRom/128 27 ns 27 ns +BM_CatmullRom/256 30 ns 30 ns +BM_CatmullRom/512 32 ns 31 ns +BM_CatmullRom/1024 33 ns 33 ns +BM_CatmullRom/2048 34 ns 34 ns +BM_CatmullRom/4096 36 ns 36 ns +BM_CatmullRom/8192 38 ns 38 ns +BM_CatmullRom/16384 39 ns 39 ns +BM_CatmullRom/32768 40 ns 40 ns +BM_CatmullRom/65536 45 ns 44 ns +BM_CatmullRom/131072 46 ns 46 ns +BM_CatmullRom/262144 50 ns 50 ns +BM_CatmullRom/524288 53 ns 52 ns +BM_CatmullRom/1048576 58 ns 57 ns +BM_CatmullRom_BigO 2.97 lgN 2.97 lgN +BM_CatmullRom_RMS 19 % 19 % + [section:catmull_rom_refs References] * Cem Yuksel, Scott Schaefer, and John Keyser, ['Parameterization and applications of Catmull–Rom curves], Computer-Aided Design 43 (2011) 747–755. diff --git a/include/boost/math/interpolators/catmull_rom.hpp b/include/boost/math/interpolators/catmull_rom.hpp index d827aef5a..57179ef3a 100644 --- a/include/boost/math/interpolators/catmull_rom.hpp +++ b/include/boost/math/interpolators/catmull_rom.hpp @@ -115,60 +115,54 @@ Point catmull_rom::operator()(Real s) const //Now *it >= s. We want the index such that m_s[i] <= s < m_s[i+1]: size_t i = std::distance(m_s.begin(), it - 1); // We'll keep the assert in here a while until we feel good that we've understood this algorithm. - assert(m_s[i] <= s && s < m_s[i+1]); - Point A1; + //assert(m_s[i] <= s && s < m_s[i+1]); + + // Only denom21 is used twice: + Real denom21 = 1/(m_s[i+1] - m_s[i]); + Real s0s = m_s[i-1] - s; + Real s1s = m_s[i] - s; + Real s2s = m_s[i+1] - s; + Real s3s = m_s[i+2] - s; + + Point A1_or_A3; Real denom = 1/(m_s[i] - m_s[i-1]); - Real k1 = (m_s[i]-s)*denom; - Real k2 = (s - m_s[i-1])*denom; - for (size_t j = 0; j < dimension; ++j) + for(size_t j = 0; j < dimension; ++j) { - A1[j] = k1*m_pnts[i-1][j] + k2*m_pnts[i][j]; + A1_or_A3[j] = denom*(s1s*m_pnts[i-1][j] - s0s*m_pnts[i][j]); } - Point A2; - denom = 1/(m_s[i+1] - m_s[i]); - k1 = (m_s[i+1]-s)*denom; - k2 = (s - m_s[i])*denom; - for (size_t j = 0; j < dimension; ++j) + Point A2_or_B2; + for(size_t j = 0; j < dimension; ++j) { - A2[j] = k1*m_pnts[i][j] + k2*m_pnts[i+1][j]; + A2_or_B2[j] = denom21*(s2s*m_pnts[i][j] - s1s*m_pnts[i+1][j]); } - Point B1; - for (size_t j = 0; j < dimension; ++j) + Point B1_or_C; + denom = 1/(m_s[i+1] - m_s[i-1]); + for(size_t j = 0; j < dimension; ++j) { - B1[j] = k1*A1[j] + k2*A2[j]; + B1_or_C[j] = denom*(s2s*A1_or_A3[j] - s0s*A2_or_B2[j]); } - Point A3; denom = 1/(m_s[i+2] - m_s[i+1]); - k1 = (m_s[i+2]-s)*denom; - k2 = (s - m_s[i+1])*denom; - for (size_t j = 0; j < dimension; ++j) + for(size_t j = 0; j < dimension; ++j) { - A3[j] = k1*m_pnts[i+1][j] + k2*m_pnts[i+2][j]; + A1_or_A3[j] = denom*(s3s*m_pnts[i+1][j] - s2s*m_pnts[i+2][j]); } Point B2; denom = 1/(m_s[i+2] - m_s[i]); - k1 = (m_s[i+2]-s)*denom; - k2 = (s - m_s[i])*denom; - for (size_t j = 0; j < dimension; ++j) + for(size_t j = 0; j < dimension; ++j) { - B2[j] = k1*A2[j] + k2*A3[j]; + B2[j] = denom*(s3s*A2_or_B2[j] - s1s*A1_or_A3[j]); } - Point C; - denom = 1/(m_s[i+1] - m_s[i]); - k1 = (m_s[i+1]-s)*denom; - k2 = (s - m_s[i])*denom; - for (size_t j = 0; j < dimension; ++j) + for(size_t j = 0; j < dimension; ++j) { - C[j] = k1*B1[j] + k2*B2[j]; + B1_or_C[j] = denom21*(s2s*B1_or_C[j] - s1s*B2[j]); } - - return C; + return B1_or_C; } template diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index 554b42c29..8a07ba357 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -51,6 +51,9 @@ void test_alpha_distance() template void test_linear() { + std::cout << "Testing that the Catmull-Rom spline interpolates linear functions correctly in dimension on type " + << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = 10*std::numeric_limits::epsilon(); std::vector> v(4); v[0] = {0,0,0}; From 50af4ad3e7b8ea07877ec84d436a06ed843e9455 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sun, 24 Dec 2017 11:49:04 -0700 Subject: [PATCH 07/20] Test that circles are correctly interpolated by the Catmull-Rom class. --- test/catmull_rom_test.cpp | 54 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index 8a07ba357..4d9edc7d2 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -51,7 +51,7 @@ void test_alpha_distance() template void test_linear() { - std::cout << "Testing that the Catmull-Rom spline interpolates linear functions correctly in dimension on type " + std::cout << "Testing that the Catmull-Rom spline interpolates linear functions correctly on type " << boost::typeindex::type_id().pretty_name() << "\n"; Real tol = 10*std::numeric_limits::epsilon(); @@ -112,6 +112,53 @@ void test_linear() } +template +void test_circle() +{ + std::cout << "Testing that the Catmull-Rom spline interpolates circles correctly on type " + << boost::typeindex::type_id().pretty_name() << "\n"; + + Real tol = 10*std::numeric_limits::epsilon(); + std::vector> v(20*sizeof(Real)); + for (size_t i = 0; i < v.size(); ++i) + { + Real theta = ((Real) i/ (Real) v.size())*2*M_PI; + v[i] = {cos(theta), sin(theta)}; + } + catmull_rom, 2> circle(v.data(), v.size(), true); + + // Interpolation condition: + for (size_t i = 0; i < v.size(); ++i) + { + Real s = circle.parameter_at_point(i); + auto p = circle(s); + Real x = p[0]; + Real y = p[1]; + if (abs(x) < std::numeric_limits::epsilon()) + { + BOOST_CHECK_SMALL(v[i][0], tol); + } + if (abs(y) < std::numeric_limits::epsilon()) + { + BOOST_CHECK_SMALL(v[i][1], tol); + } + else + { + BOOST_CHECK_CLOSE_FRACTION(x, v[i][0], tol); + BOOST_CHECK_CLOSE_FRACTION(y, v[i][1], tol); + } + } + + Real max_s = circle.max_parameter(); + for(Real s = 0; s < max_s; s += 0.01) + { + auto p = circle(s); + Real x = p[0]; + Real y = p[1]; + BOOST_CHECK_CLOSE_FRACTION(x*x+y*y, 1, 0.001); + } +} + template void test_affine_invariance() { @@ -192,6 +239,11 @@ BOOST_AUTO_TEST_CASE(catmull_rom_test) test_linear(); test_linear(); + test_circle(); + test_circle(); + test_circle(); + test_circle(); + test_affine_invariance(); test_affine_invariance(); test_affine_invariance(); From 577ee425a452587962226b201576a1e81404b532 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sun, 24 Dec 2017 13:22:30 -0700 Subject: [PATCH 08/20] Add test showing that helices are interpolated correctly by the Catmull-Rom class. --- test/catmull_rom_test.cpp | 66 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index 4d9edc7d2..1704f7d0d 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -221,6 +221,69 @@ void test_affine_invariance() } } +template +void test_helix() +{ + std::cout << "Testing that the Catmull-Rom spline interpolates helices correctly on type " + << boost::typeindex::type_id().pretty_name() << "\n"; + + Real tol = 1000*std::numeric_limits::epsilon(); + std::vector> v(1000*sizeof(Real)); + for (size_t i = 0; i < v.size(); ++i) + { + Real theta = ((Real) i/ (Real) v.size())*2*M_PI; + v[i] = {cos(theta), sin(theta), theta}; + } + catmull_rom, 3> helix(v.data(), v.size()); + + // Interpolation condition: + for (size_t i = 0; i < v.size(); ++i) + { + Real s = helix.parameter_at_point(i); + auto p = helix(s); + Real t = p[2]; + + Real x = p[0]; + Real y = p[1]; + if (abs(x) < tol) + { + BOOST_CHECK_SMALL(cos(t), tol); + } + if (abs(y) < tol) + { + BOOST_CHECK_SMALL(sin(t), tol); + } + else + { + BOOST_CHECK_CLOSE_FRACTION(x, cos(t), tol); + BOOST_CHECK_CLOSE_FRACTION(y, sin(t), tol); + } + } + + Real max_s = helix.max_parameter(); + for(Real s = helix.parameter_at_point(1); s < max_s; s += 0.01) + { + auto p = helix(s); + Real x = p[0]; + Real y = p[1]; + Real t = p[2]; + BOOST_CHECK_CLOSE_FRACTION(x*x+y*y, (Real) 1, (Real) 0.01); + if (abs(x) < 0.01) + { + BOOST_CHECK_SMALL(cos(t), (Real) 0.01); + } + if (abs(y) < 0.01) + { + BOOST_CHECK_SMALL(sin(t), (Real) 0.01); + } + else + { + BOOST_CHECK_CLOSE_FRACTION(x, cos(t), (Real) 0.01); + BOOST_CHECK_CLOSE_FRACTION(y, sin(t), (Real) 0.01); + } + } +} + template void test_data_representations() { @@ -244,6 +307,9 @@ BOOST_AUTO_TEST_CASE(catmull_rom_test) test_circle(); test_circle(); + test_helix(); + test_helix(); + test_affine_invariance(); test_affine_invariance(); test_affine_invariance(); From e137fffa632dbb6e48e58341803866bcf7053caf Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Wed, 10 Jan 2018 12:03:42 -0600 Subject: [PATCH 09/20] Loosen up tolerance in unit tests to prevent spurious failures. --- doc/interpolators/catmull_rom.qbk | 64 +++++++++++++++---------------- test/catmull_rom_test.cpp | 16 ++++---- 2 files changed, 40 insertions(+), 40 deletions(-) diff --git a/doc/interpolators/catmull_rom.qbk b/doc/interpolators/catmull_rom.qbk index 06179a2cb..d4fea3296 100644 --- a/doc/interpolators/catmull_rom.qbk +++ b/doc/interpolators/catmull_rom.qbk @@ -45,7 +45,7 @@ Catmull-Rom splines enjoy the following nice properties: The `catmull_rom` class provided by boost creates a cubic Catmull-Rom spline from an array of points in any dimension. Since there are numerous ways to represent a point in /n/-dimensional space, the class attempts to be flexible by templating on the point type. -The requirement of the point type is that it must have a dereference operator defined (e.g., `p[0]` is not a syntax error), +The requirement of the point type is that it must have a dereference operator defined (i.e., `p[0]` is not a syntax error), it must be able to be dereferenced up to `dimension -1`, and `p[i]` is of type `Real`. The basic usage is shown here: @@ -58,7 +58,7 @@ The basic usage is shown here: // Interpolate at s = 0.1: auto point = cr(0.1); -The spline can be either open or *closed*, meaning that there is some /t/ such /P(t) = P(0)/. +The spline can be either open or *closed*, closed meaning that there is some /t/ such that /P(t) = P(0)/. The default is open, but this can be easily changed: // closed = true @@ -115,36 +115,36 @@ The following performance numbers were generated for a call to the Catmull-Rom i The number that follows the slash is the number of points passed to the interpolant. -Run on 2700 MHz CPU -CPU Caches: - L1 Data 32K (x2) - L1 Instruction 32K (x2) - L2 Unified 262K (x2) - L3 Unified 3145K (x1) ---------------------------------------------------------- -Benchmark Time CPU ---------------------------------------------------------- -BM_CatmullRom/4 20 ns 20 ns -BM_CatmullRom/8 21 ns 21 ns -BM_CatmullRom/16 23 ns 23 ns -BM_CatmullRom/32 24 ns 24 ns -BM_CatmullRom/64 27 ns 27 ns -BM_CatmullRom/128 27 ns 27 ns -BM_CatmullRom/256 30 ns 30 ns -BM_CatmullRom/512 32 ns 31 ns -BM_CatmullRom/1024 33 ns 33 ns -BM_CatmullRom/2048 34 ns 34 ns -BM_CatmullRom/4096 36 ns 36 ns -BM_CatmullRom/8192 38 ns 38 ns -BM_CatmullRom/16384 39 ns 39 ns -BM_CatmullRom/32768 40 ns 40 ns -BM_CatmullRom/65536 45 ns 44 ns -BM_CatmullRom/131072 46 ns 46 ns -BM_CatmullRom/262144 50 ns 50 ns -BM_CatmullRom/524288 53 ns 52 ns -BM_CatmullRom/1048576 58 ns 57 ns -BM_CatmullRom_BigO 2.97 lgN 2.97 lgN -BM_CatmullRom_RMS 19 % 19 % + Run on 2700 MHz CPU + CPU Caches: + L1 Data 32K (x2) + L1 Instruction 32K (x2) + L2 Unified 262K (x2) + L3 Unified 3145K (x1) + --------------------------------------------------------- + Benchmark Time CPU + --------------------------------------------------------- + BM_CatmullRom/4 20 ns 20 ns + BM_CatmullRom/8 21 ns 21 ns + BM_CatmullRom/16 23 ns 23 ns + BM_CatmullRom/32 24 ns 24 ns + BM_CatmullRom/64 27 ns 27 ns + BM_CatmullRom/128 27 ns 27 ns + BM_CatmullRom/256 30 ns 30 ns + BM_CatmullRom/512 32 ns 31 ns + BM_CatmullRom/1024 33 ns 33 ns + BM_CatmullRom/2048 34 ns 34 ns + BM_CatmullRom/4096 36 ns 36 ns + BM_CatmullRom/8192 38 ns 38 ns + BM_CatmullRom/16384 39 ns 39 ns + BM_CatmullRom/32768 40 ns 40 ns + BM_CatmullRom/65536 45 ns 44 ns + BM_CatmullRom/131072 46 ns 46 ns + BM_CatmullRom/262144 50 ns 50 ns + BM_CatmullRom/524288 53 ns 52 ns + BM_CatmullRom/1048576 58 ns 57 ns + BM_CatmullRom_BigO 2.97 lgN 2.97 lgN + BM_CatmullRom_RMS 19 % 19 % [section:catmull_rom_refs References] diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index 1704f7d0d..d1061c19d 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -213,9 +213,9 @@ void test_affine_invariance() for (size_t j = 0; j < dimension; ++j) { BOOST_CHECK_CLOSE_FRACTION(p0[j] + affine_shift[j], p1[j], tol); - if (abs(tangent0[j]) > tol) + if (abs(tangent0[j]) > 5000*tol) { - BOOST_CHECK_CLOSE_FRACTION(tangent0[j], tangent1[j], 20*tol); + BOOST_CHECK_CLOSE_FRACTION(tangent0[j], tangent1[j], 5000*tol); } } } @@ -227,8 +227,8 @@ void test_helix() std::cout << "Testing that the Catmull-Rom spline interpolates helices correctly on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = 1000*std::numeric_limits::epsilon(); - std::vector> v(1000*sizeof(Real)); + Real tol = 0.001; + std::vector> v(2000*sizeof(Real)); for (size_t i = 0; i < v.size(); ++i) { Real theta = ((Real) i/ (Real) v.size())*2*M_PI; @@ -270,16 +270,16 @@ void test_helix() BOOST_CHECK_CLOSE_FRACTION(x*x+y*y, (Real) 1, (Real) 0.01); if (abs(x) < 0.01) { - BOOST_CHECK_SMALL(cos(t), (Real) 0.01); + BOOST_CHECK_SMALL(cos(t), (Real) 0.05); } if (abs(y) < 0.01) { - BOOST_CHECK_SMALL(sin(t), (Real) 0.01); + BOOST_CHECK_SMALL(sin(t), (Real) 0.05); } else { - BOOST_CHECK_CLOSE_FRACTION(x, cos(t), (Real) 0.01); - BOOST_CHECK_CLOSE_FRACTION(y, sin(t), (Real) 0.01); + BOOST_CHECK_CLOSE_FRACTION(x, cos(t), (Real) 0.05); + BOOST_CHECK_CLOSE_FRACTION(y, sin(t), (Real) 0.05); } } } From 08550d4064e8e44a3e43aac43102ddbdf21c1f04 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Wed, 10 Jan 2018 14:08:23 -0600 Subject: [PATCH 10/20] 3 template arguments -> 1 template argument. --- .../boost/math/interpolators/catmull_rom.hpp | 103 +++++++++--------- test/catmull_rom_test.cpp | 24 ++-- 2 files changed, 66 insertions(+), 61 deletions(-) diff --git a/include/boost/math/interpolators/catmull_rom.hpp b/include/boost/math/interpolators/catmull_rom.hpp index 57179ef3a..c134dfcbd 100644 --- a/include/boost/math/interpolators/catmull_rom.hpp +++ b/include/boost/math/interpolators/catmull_rom.hpp @@ -12,57 +12,58 @@ #include #include #include +#include namespace boost{ namespace math{ namespace detail { - template - Real alpha_distance(Point const & p1, Point const & p2, Real alpha) + template + typename Point::value_type alpha_distance(Point const & p1, Point const & p2, typename Point::value_type alpha) { using std::pow; - Real dsq = 0; - for (size_t i = 0; i < dimension; ++i) + using std::size; + typename Point::value_type dsq = 0; + for (size_t i = 0; i < size(p1); ++i) { - Real dx = p1[i] - p2[i]; + typename Point::value_type dx = p1[i] - p2[i]; dsq += dx*dx; } return pow(dsq, alpha/2); } } -template +template class catmull_rom { public: - catmull_rom(const Point* const points, size_t num_pnts, bool closed = false, Real alpha = (Real) 1/ (Real) 2); + catmull_rom(const Point* const points, size_t num_pnts, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2); - Real max_parameter() const + typename Point::value_type max_parameter() const { return m_max_s; } - Real parameter_at_point(size_t i) const + typename Point::value_type parameter_at_point(size_t i) const { return m_s[i+1]; } - Point operator()(Real s) const; + Point operator()(const typename Point::value_type s) const; - Point prime(Real s) const; + Point prime(const typename Point::value_type s) const; private: std::vector m_pnts; - std::vector m_s; - Real m_max_s; + std::vector m_s; + typename Point::value_type m_max_s; }; -template -catmull_rom::catmull_rom(const Point* const points, size_t num_pnts, bool closed, Real alpha) +template +catmull_rom::catmull_rom(const Point* const points, size_t num_pnts, bool closed, typename Point::value_type alpha) { - static_assert(dimension > 0, "The dimension of the Catmull-Rom spline must be > 0\n"); if (num_pnts < 4) { throw std::domain_error("The Catmull-Rom curve requires at least 4 points.\n"); @@ -78,16 +79,16 @@ catmull_rom::catmull_rom(const Point* const points, size } m_pnts[num_pnts+1] = points[0]; m_pnts[num_pnts+2] = points[1]; - m_s[0] = -detail::alpha_distance(m_pnts[0], m_pnts[1], alpha); - if (abs(m_s[0]) < std::numeric_limits::epsilon()) + m_s[0] = -detail::alpha_distance(m_pnts[0], m_pnts[1], alpha); + if (abs(m_s[0]) < std::numeric_limits::epsilon()) { throw std::domain_error("The first and last point should not be the same.\n"); } m_s[1] = 0; for (size_t i = 2; i < m_s.size(); ++i) { - Real d = detail::alpha_distance(m_pnts[i], m_pnts[i-1], alpha); - if (abs(d) < std::numeric_limits::epsilon()) + typename Point::value_type d = detail::alpha_distance(m_pnts[i], m_pnts[i-1], alpha); + if (abs(d) < std::numeric_limits::epsilon()) { throw std::domain_error("The control points of the Catmull-Rom curve are too close together; this will lead to ill-conditioning.\n"); } @@ -104,9 +105,10 @@ catmull_rom::catmull_rom(const Point* const points, size } -template -Point catmull_rom::operator()(Real s) const +template +Point catmull_rom::operator()(const typename Point::value_type s) const { + using std::size; if (s < 0 || s > m_max_s) { throw std::domain_error("Parameter outside bounds.\n"); @@ -118,46 +120,46 @@ Point catmull_rom::operator()(Real s) const //assert(m_s[i] <= s && s < m_s[i+1]); // Only denom21 is used twice: - Real denom21 = 1/(m_s[i+1] - m_s[i]); - Real s0s = m_s[i-1] - s; - Real s1s = m_s[i] - s; - Real s2s = m_s[i+1] - s; - Real s3s = m_s[i+2] - s; + typename Point::value_type denom21 = 1/(m_s[i+1] - m_s[i]); + typename Point::value_type s0s = m_s[i-1] - s; + typename Point::value_type s1s = m_s[i] - s; + typename Point::value_type s2s = m_s[i+1] - s; + typename Point::value_type s3s = m_s[i+2] - s; Point A1_or_A3; - Real denom = 1/(m_s[i] - m_s[i-1]); - for(size_t j = 0; j < dimension; ++j) + typename Point::value_type denom = 1/(m_s[i] - m_s[i-1]); + for(size_t j = 0; j < size(m_pnts[0]); ++j) { A1_or_A3[j] = denom*(s1s*m_pnts[i-1][j] - s0s*m_pnts[i][j]); } Point A2_or_B2; - for(size_t j = 0; j < dimension; ++j) + for(size_t j = 0; j < size(m_pnts[0]); ++j) { A2_or_B2[j] = denom21*(s2s*m_pnts[i][j] - s1s*m_pnts[i+1][j]); } Point B1_or_C; denom = 1/(m_s[i+1] - m_s[i-1]); - for(size_t j = 0; j < dimension; ++j) + for(size_t j = 0; j < size(m_pnts[0]); ++j) { B1_or_C[j] = denom*(s2s*A1_or_A3[j] - s0s*A2_or_B2[j]); } denom = 1/(m_s[i+2] - m_s[i+1]); - for(size_t j = 0; j < dimension; ++j) + for(size_t j = 0; j < size(m_pnts[0]); ++j) { A1_or_A3[j] = denom*(s3s*m_pnts[i+1][j] - s2s*m_pnts[i+2][j]); } Point B2; denom = 1/(m_s[i+2] - m_s[i]); - for(size_t j = 0; j < dimension; ++j) + for(size_t j = 0; j < size(m_pnts[0]); ++j) { B2[j] = denom*(s3s*A2_or_B2[j] - s1s*A1_or_A3[j]); } - for(size_t j = 0; j < dimension; ++j) + for(size_t j = 0; j < size(m_pnts[0]); ++j) { B1_or_C[j] = denom21*(s2s*B1_or_C[j] - s1s*B2[j]); } @@ -165,9 +167,10 @@ Point catmull_rom::operator()(Real s) const return B1_or_C; } -template -Point catmull_rom::prime(Real s) const +template +Point catmull_rom::prime(const typename Point::value_type s) const { + using std::size; // https://math.stackexchange.com/questions/843595/how-can-i-calculate-the-derivative-of-a-catmull-rom-spline-with-nonuniform-param // http://denkovacs.com/2016/02/catmull-rom-spline-derivatives/ if (s < 0 || s > m_max_s) @@ -180,16 +183,16 @@ Point catmull_rom::prime(Real s) const // We'll keep the assert in here a while until we feel good that we've understood this algorithm. assert(m_s[i] <= s && s < m_s[i+1]); Point A1; - Real denom = 1/(m_s[i] - m_s[i-1]); - Real k1 = (m_s[i]-s)*denom; - Real k2 = (s - m_s[i-1])*denom; - for (size_t j = 0; j < dimension; ++j) + typename Point::value_type denom = 1/(m_s[i] - m_s[i-1]); + typename Point::value_type k1 = (m_s[i]-s)*denom; + typename Point::value_type k2 = (s - m_s[i-1])*denom; + for (size_t j = 0; j < size(m_pnts[0]); ++j) { A1[j] = k1*m_pnts[i-1][j] + k2*m_pnts[i][j]; } Point A1p; - for (size_t j = 0; j < dimension; ++j) + for (size_t j = 0; j < size(m_pnts[0]); ++j) { A1p[j] = denom*(m_pnts[i][j] - m_pnts[i-1][j]); } @@ -198,20 +201,20 @@ Point catmull_rom::prime(Real s) const denom = 1/(m_s[i+1] - m_s[i]); k1 = (m_s[i+1]-s)*denom; k2 = (s - m_s[i])*denom; - for (size_t j = 0; j < dimension; ++j) + for (size_t j = 0; j < size(m_pnts[0]); ++j) { A2[j] = k1*m_pnts[i][j] + k2*m_pnts[i+1][j]; } Point A2p; - for (size_t j = 0; j < dimension; ++j) + for (size_t j = 0; j < size(m_pnts[0]); ++j) { A2p[j] = denom*(m_pnts[i+1][j] - m_pnts[i][j]); } Point B1; - for (size_t j = 0; j < dimension; ++j) + for (size_t j = 0; j < size(m_pnts[0]); ++j) { B1[j] = k1*A1[j] + k2*A2[j]; } @@ -220,13 +223,13 @@ Point catmull_rom::prime(Real s) const denom = 1/(m_s[i+2] - m_s[i+1]); k1 = (m_s[i+2]-s)*denom; k2 = (s - m_s[i+1])*denom; - for (size_t j = 0; j < dimension; ++j) + for (size_t j = 0; j < size(m_pnts[0]); ++j) { A3[j] = k1*m_pnts[i+1][j] + k2*m_pnts[i+2][j]; } Point A3p; - for (size_t j = 0; j < dimension; ++j) + for (size_t j = 0; j < size(m_pnts[0]); ++j) { A3p[j] = denom*(m_pnts[i+2][j] - m_pnts[i+1][j]); } @@ -235,28 +238,28 @@ Point catmull_rom::prime(Real s) const denom = 1/(m_s[i+2] - m_s[i]); k1 = (m_s[i+2]-s)*denom; k2 = (s - m_s[i])*denom; - for (size_t j = 0; j < dimension; ++j) + for (size_t j = 0; j < size(m_pnts[0]); ++j) { B2[j] = k1*A2[j] + k2*A3[j]; } Point B1p; denom = 1/(m_s[i+1] - m_s[i-1]); - for (size_t j = 0; j < dimension; ++j) + for (size_t j = 0; j < size(m_pnts[0]); ++j) { B1p[j] = denom*(A2[j] - A1[j] + (m_s[i+1]- s)*A1p[j] + (s-m_s[i-1])*A2p[j]); } Point B2p; denom = 1/(m_s[i+2] - m_s[i]); - for (size_t j = 0; j < dimension; ++j) + for (size_t j = 0; j < size(m_pnts[0]); ++j) { B2p[j] = denom*(A3[j] - A2[j] + (m_s[i+2] - s)*A2p[j] + (s - m_s[i])*A3p[j]); } Point Cp; denom = 1/(m_s[i+1] - m_s[i]); - for (size_t j = 0; j < dimension; ++j) + for (size_t j = 0; j < size(m_pnts[0]); ++j) { Cp[j] = denom*(B2[j] - B1[j] + (m_s[i+1] - s)*B1p[j] + (s - m_s[i])*B2p[j]); } diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index d1061c19d..5d2c1d296 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -27,23 +27,23 @@ void test_alpha_distance() std::array v1 = {0,0,0}; std::array v2 = {1,0,0}; Real alpha = 0.5; - Real d = boost::math::detail::alpha_distance, 3>(v1, v2, alpha); + Real d = boost::math::detail::alpha_distance>(v1, v2, alpha); BOOST_CHECK_CLOSE_FRACTION(d, 1, tol); - d = boost::math::detail::alpha_distance, 3>(v1, v2, 0); + d = boost::math::detail::alpha_distance>(v1, v2, 0.0); BOOST_CHECK_CLOSE_FRACTION(d, 1, tol); - d = boost::math::detail::alpha_distance, 3>(v1, v2, 1); + d = boost::math::detail::alpha_distance>(v1, v2, 1.0); BOOST_CHECK_CLOSE_FRACTION(d, 1, tol); v2[0] = 2; - d = boost::math::detail::alpha_distance, 3>(v1, v2, alpha); + d = boost::math::detail::alpha_distance>(v1, v2, alpha); BOOST_CHECK_CLOSE_FRACTION(d, pow(2, (Real)1/ (Real) 2), tol); - d = boost::math::detail::alpha_distance, 3>(v1, v2, 0); + d = boost::math::detail::alpha_distance>(v1, v2, 0.0); BOOST_CHECK_CLOSE_FRACTION(d, 1, tol); - d = boost::math::detail::alpha_distance, 3>(v1, v2, 1); + d = boost::math::detail::alpha_distance>(v1, v2, 1.0); BOOST_CHECK_CLOSE_FRACTION(d, 2, tol); } @@ -60,7 +60,7 @@ void test_linear() v[1] = {1,0,0}; v[2] = {2,0,0}; v[3] = {3,0,0}; - catmull_rom, 3> cr(v.data(), v.size()); + catmull_rom> cr(v.data(), v.size()); // Test that the interpolation condition is obeyed: BOOST_CHECK_CLOSE_FRACTION(cr.max_parameter(), 3, tol); @@ -125,7 +125,7 @@ void test_circle() Real theta = ((Real) i/ (Real) v.size())*2*M_PI; v[i] = {cos(theta), sin(theta)}; } - catmull_rom, 2> circle(v.data(), v.size(), true); + catmull_rom> circle(v.data(), v.size(), true); // Interpolation condition: for (size_t i = 0; i < v.size(); ++i) @@ -159,6 +159,7 @@ void test_circle() } } + template void test_affine_invariance() { @@ -189,7 +190,7 @@ void test_affine_invariance() affine_shift[j] = gen()*inv_denom; } - catmull_rom, dimension> cr1(v.data(), v.size()); + catmull_rom> cr1(v.data(), v.size()); for(size_t i = 0; i< v.size(); ++i) { @@ -199,7 +200,7 @@ void test_affine_invariance() } } - catmull_rom, dimension> cr2(v.data(), v.size()); + catmull_rom> cr2(v.data(), v.size()); BOOST_CHECK_CLOSE_FRACTION(cr1.max_parameter(), cr2.max_parameter(), tol); @@ -234,7 +235,7 @@ void test_helix() Real theta = ((Real) i/ (Real) v.size())*2*M_PI; v[i] = {cos(theta), sin(theta), theta}; } - catmull_rom, 3> helix(v.data(), v.size()); + catmull_rom> helix(v.data(), v.size()); // Interpolation condition: for (size_t i = 0; i < v.size(); ++i) @@ -307,6 +308,7 @@ BOOST_AUTO_TEST_CASE(catmull_rom_test) test_circle(); test_circle(); + test_helix(); test_helix(); From 1781d47b78f57be076993f9558280a6b1d612a58 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sat, 27 Jan 2018 14:44:15 -0600 Subject: [PATCH 11/20] [ci skip] Clarify conditions on the template point type. Add include and concept test. --- doc/interpolators/catmull_rom.qbk | 79 ++++++++++++++++--- example/catmull_rom_example.cpp | 2 +- .../boost/math/interpolators/catmull_rom.hpp | 8 +- test/catmull_rom_test.cpp | 67 +++++++++++++++- .../compile_test/catmull_rom_concept_test.cpp | 21 +++++ test/compile_test/catmull_rom_incl_test.cpp | 28 +++++++ 6 files changed, 186 insertions(+), 19 deletions(-) create mode 100644 test/compile_test/catmull_rom_concept_test.cpp create mode 100644 test/compile_test/catmull_rom_incl_test.cpp diff --git a/doc/interpolators/catmull_rom.qbk b/doc/interpolators/catmull_rom.qbk index d4fea3296..a079f4acb 100644 --- a/doc/interpolators/catmull_rom.qbk +++ b/doc/interpolators/catmull_rom.qbk @@ -14,11 +14,14 @@ #include namespace boost{ namespace math{ - template + + template class catmull_rom { public: - catmull_rom(const Real* const points, size_t num_points, bool closed = false, Real alpha = (Real) 1/ (Real) 2) + catmull_rom(const Point* const points, size_t num_points, bool closed = false, Real alpha = (Real) 1/ (Real) 2) + + catmull_rom(std::vector&& points, bool closed = false, Real alpha = (Real) 1/ (Real) 2) Real operator()(Real s) const; @@ -35,7 +38,7 @@ namespace boost{ namespace math{ [heading Description] Catmull-Rom splines are a family of interpolating curves which are commonly used in computer graphics and animation. -Catmull-Rom splines enjoy the following nice properties: +Catmull-Rom splines enjoy the following properties: * Affine invariance: The interpolant commutes with affine transformations. * Local support of the basis functions: This gives stability and fast evaluation. @@ -45,8 +48,9 @@ Catmull-Rom splines enjoy the following nice properties: The `catmull_rom` class provided by boost creates a cubic Catmull-Rom spline from an array of points in any dimension. Since there are numerous ways to represent a point in /n/-dimensional space, the class attempts to be flexible by templating on the point type. -The requirement of the point type is that it must have a dereference operator defined (i.e., `p[0]` is not a syntax error), -it must be able to be dereferenced up to `dimension -1`, and `p[i]` is of type `Real`. +The requirements on the point type are discussing in more detail below, but roughly, it must have a dereference operator defined (i.e., `p[0]` is not a syntax error), +it must be able to be dereferenced up to `dimension -1`, and `p[i]` is of type `Real`, define `value_type`, and the free function `size()`. +These requirements are met by `std::vector` and `std::array`. The basic usage is shown here: std::vector> points(4); @@ -54,10 +58,14 @@ The basic usage is shown here: points[1] = {1,0,0}; points[2] = {0,1,0}; points[3] = {0,0,1}; - catmull_rom, 3> cr(points.data(), points.size()); + catmull_rom> cr(points.data(), points.size()); // Interpolate at s = 0.1: auto point = cr(0.1); +Memory conscious programmers may enjoy using the move constructor instead: + + catmull_rom> cr(std::move(points)); + The spline can be either open or *closed*, closed meaning that there is some /t/ such that /P(t) = P(0)/. The default is open, but this can be easily changed: @@ -89,16 +97,16 @@ If the curve is open, then `cr(max_s) = pf`, where `pf` is the final point on th The Catmull-Rom curve admits an infinite number of parameterizations. The default parameterization of the `catmull_rom` class is the so-called *centripedal* parameterization. -This parameterization has been shown to be the only parameterization that does not have cusps or self-intersections within segments. +This parameterization has been shown to be the only parameterization that does not form cusps or self-intersections within segments. However, for advanced users, other parameterizations can be chosen using the /alpha/ parameter: // alpha = 1 is the "chordal" parameterization. - catmull_rom, 3> cr(points.data(), points.size(), false, (Real) 1); + catmull_rom> cr(points.data(), points.size(), false, 1.0); -Finally, the tangent vector to any point of the curve can be discovered via +Finally, the tangent vector to any point of the curve can be computed via - Real s = 0.1; + double s = 0.1; Point tangent = cr.prime(s); Since the magnitude of the tangent vector is dependent on the parameterization, @@ -146,10 +154,59 @@ The number that follows the slash is the number of points passed to the interpol BM_CatmullRom_BigO 2.97 lgN 2.97 lgN BM_CatmullRom_RMS 19 % 19 % +[endsect] + +[section:point_type Point types] + +We have already discussed that certain conditions on the `Point` type template argument must be obeyed. +The following shows a custom point type in 3D which can be used as a template argument to Catmull-Rom: + + 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; + }; + + + // Must define the free function "size()": + template + constexpr size_t size(const mypoint3d& c) + { + return 3; + } + +These conditions are satisfied by both `std::array` and `std::vector`, but it may nonetheless be useful to define your own point class so that (say) you can define geometric distance between them. + +[endsect] + [section:catmull_rom_refs References] * Cem Yuksel, Scott Schaefer, and John Keyser, ['Parameterization and applications of Catmull–Rom curves], Computer-Aided Design 43 (2011) 747–755. * Phillip J. Barry and Ronald N. Goldman, ['A Recursive Evaluation Algorithm for a Class of Catmull-Rom Splines], Computer Graphics, Volume 22, Number 4, August 1988 [endsect] -[endsect] [/section:catmull_rom Catmull-Rom Splines] +[endsect] +[/section:catmull_rom Catmull-Rom Splines] diff --git a/example/catmull_rom_example.cpp b/example/catmull_rom_example.cpp index acdea3360..c7a48da09 100644 --- a/example/catmull_rom_example.cpp +++ b/example/catmull_rom_example.cpp @@ -27,7 +27,7 @@ int main() spiral_points[i] = {theta*cos(theta), theta*sin(theta)}; } - auto archimedean = catmull_rom, 2>(spiral_points.data(), spiral_points.size()); + auto archimedean = catmull_rom>(std::move(spiral_points)); double max_s = archimedean.max_parameter(); std::cout << "Max s = " << max_s << std::endl; for (double s = 0; s < max_s; s += 0.01) diff --git a/include/boost/math/interpolators/catmull_rom.hpp b/include/boost/math/interpolators/catmull_rom.hpp index c134dfcbd..dd9435137 100644 --- a/include/boost/math/interpolators/catmull_rom.hpp +++ b/include/boost/math/interpolators/catmull_rom.hpp @@ -10,7 +10,7 @@ #define BOOST_MATH_INTERPOLATORS_CATMULL_ROM #include -#include +#include #include #include @@ -40,6 +40,8 @@ public: catmull_rom(const Point* const points, size_t num_pnts, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2); + catmull_rom(std::vector&& points, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2) : catmull_rom(points.data(), points.size(), closed, alpha) {} + typename Point::value_type max_parameter() const { return m_max_s; @@ -116,8 +118,6 @@ Point catmull_rom::operator()(const typename Point::value_type s) const auto it = std::upper_bound(m_s.begin(), m_s.end(), s); //Now *it >= s. We want the index such that m_s[i] <= s < m_s[i+1]: size_t i = std::distance(m_s.begin(), it - 1); - // We'll keep the assert in here a while until we feel good that we've understood this algorithm. - //assert(m_s[i] <= s && s < m_s[i+1]); // Only denom21 is used twice: typename Point::value_type denom21 = 1/(m_s[i+1] - m_s[i]); @@ -180,8 +180,6 @@ Point catmull_rom::prime(const typename Point::value_type s) const auto it = std::upper_bound(m_s.begin(), m_s.end(), s); //Now *it >= s. We want the index such that m_s[i] <= s < m_s[i+1]: size_t i = std::distance(m_s.begin(), it - 1); - // We'll keep the assert in here a while until we feel good that we've understood this algorithm. - assert(m_s[i] <= s && s < m_s[i+1]); Point A1; typename Point::value_type denom = 1/(m_s[i] - m_s[i-1]); typename Point::value_type k1 = (m_s[i]-s)*denom; diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index 5d2c1d296..8c17fe35c 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -285,15 +285,78 @@ void test_helix() } } + +template +class mypoint3d +{ +public: + // Must define a value_type: + typedef Real value_type; + + // Regular constructor: + 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]; + } + + // Array element assignment: + Real& operator[](size_t i) + { + return m_vec[i]; + } + + +private: + std::array m_vec; +}; + + +// Must define the free function "size()": +template +constexpr size_t size(const mypoint3d& c) +{ + return 3; +} + template void test_data_representations() { std::cout << "Testing that the Catmull-Rom spline works with multiple data representations.\n"; + mypoint3d p0(0.1, 0.2, 0.3); + mypoint3d p1(0.2, 0.3, 0.4); + mypoint3d p2(0.3, 0.4, 0.5); + mypoint3d p3(0.4, 0.5, 0.6); + mypoint3d p4(0.5, 0.6, 0.7); + mypoint3d p5(0.6, 0.7, 0.8); + std::vector> v{p0, p1, p2, p3, p4, p5}; + catmull_rom> cat(v.data(), v.size()); + + Real tol = 0.001; + auto p = cat(cat.parameter_at_point(0)); + BOOST_CHECK_CLOSE_FRACTION(p[0], p0[0], tol); + BOOST_CHECK_CLOSE_FRACTION(p[1], p0[1], tol); + BOOST_CHECK_CLOSE_FRACTION(p[2], p0[2], tol); + p = cat(cat.parameter_at_point(1)); + BOOST_CHECK_CLOSE_FRACTION(p[0], p1[0], tol); + BOOST_CHECK_CLOSE_FRACTION(p[1], p1[1], tol); + BOOST_CHECK_CLOSE_FRACTION(p[2], p1[2], tol); } BOOST_AUTO_TEST_CASE(catmull_rom_test) { - test_alpha_distance(); + test_data_representations(); + /*test_alpha_distance(); test_alpha_distance(); test_alpha_distance(); test_alpha_distance(); @@ -329,5 +392,5 @@ BOOST_AUTO_TEST_CASE(catmull_rom_test) test_affine_invariance(); test_affine_invariance(); test_affine_invariance(); - test_affine_invariance(); + test_affine_invariance();*/ } diff --git a/test/compile_test/catmull_rom_concept_test.cpp b/test/compile_test/catmull_rom_concept_test.cpp new file mode 100644 index 000000000..a9bd28405 --- /dev/null +++ b/test/compile_test/catmull_rom_concept_test.cpp @@ -0,0 +1,21 @@ +// Copyright Nick Thompson 2018. +// 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 + +void compile_and_link_test() +{ + std::vector p0{0.1, 0.2, 0.3}; + std::vector p1{0.2, 0.3, 0.4}; + std::vector p2{0.3, 0.4, 0.5}; + std::vector p3{0.4, 0.5, 0.6}; + std::vector p4{0.5, 0.6, 0.7}; + std::vector p5{0.6, 0.7, 0.8}; + std::vector> v{p0, p1, p2, p3, p4, p5}; + boost::math::catmull_rom> cat(v.data(), v.size()); + cat(0.0); + cat.prime(0.0); +} diff --git a/test/compile_test/catmull_rom_incl_test.cpp b/test/compile_test/catmull_rom_incl_test.cpp new file mode 100644 index 000000000..d419ae988 --- /dev/null +++ b/test/compile_test/catmull_rom_incl_test.cpp @@ -0,0 +1,28 @@ +// Copyright Nick Thompson 2018. +// 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) +// +// Basic sanity check that header +// #includes all the files that it needs to. +// +#include +// +// Note this header includes no other headers, this is +// important if this test is to be meaningful: +// +#include "test_compile_result.hpp" + +void compile_and_link_test() +{ + std::vector p0{0.1, 0.2, 0.3}; + std::vector p1{0.2, 0.3, 0.4}; + std::vector p2{0.3, 0.4, 0.5}; + std::vector p3{0.4, 0.5, 0.6}; + std::vector p4{0.5, 0.6, 0.7}; + std::vector p5{0.6, 0.7, 0.8}; + std::vector> v{p0, p1, p2, p3, p4, p5}; + boost::math::catmull_rom> cat(v.data(), v.size()); + check_result(cat(0.0)); + check_result(cat.prime(0.0)); +} From cea481a3d04f78444a984b4721615433e81eb6e3 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sat, 27 Jan 2018 14:57:08 -0600 Subject: [PATCH 12/20] [ci skip] Use deterministic seed for random unit test. --- test/catmull_rom_test.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index 8c17fe35c..9ad59cd6c 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -169,8 +169,7 @@ void test_affine_invariance() Real tol = 1000*std::numeric_limits::epsilon(); std::vector> v(100); - std::random_device rd; - std::mt19937_64 gen(rd()); + std::mt19937_64 gen(438232); Real inv_denom = (Real) 100/( (Real) gen.max() + (Real) 2); for(size_t j = 0; j < dimension; ++j) { @@ -356,7 +355,7 @@ void test_data_representations() BOOST_AUTO_TEST_CASE(catmull_rom_test) { test_data_representations(); - /*test_alpha_distance(); + test_alpha_distance(); test_alpha_distance(); test_alpha_distance(); test_alpha_distance(); @@ -392,5 +391,5 @@ BOOST_AUTO_TEST_CASE(catmull_rom_test) test_affine_invariance(); test_affine_invariance(); test_affine_invariance(); - test_affine_invariance();*/ + test_affine_invariance(); } From 6d04fbe33aecda9d516ca8ce457f55e900cb8451 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Mon, 29 Oct 2018 12:19:38 -0600 Subject: [PATCH 13/20] Temporarily make travis more verbose so that we can diagnose internal compiler error. --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index fc5db3b0c..8b1177b8a 100644 --- a/.travis.yml +++ b/.travis.yml @@ -584,7 +584,7 @@ script: - |- echo "using $TOOLSET : : $COMPILER : -std=$CXXSTD ;" > ~/user-config.jam - (cd libs/config/test && ../../../b2 config_info_travis_install toolset=$TOOLSET && ./config_info_travis) - - (cd libs/math/test && travis_wait 60 ../../../b2 -j3 -d+0 -q toolset=$TOOLSET $TEST_SUITE define=CI_SUPPRESS_KNOWN_ISSUES define=SLOW_COMPILER) + - (cd libs/math/test && travis_wait 60 ../../../b2 -j3 -d2 -q toolset=$TOOLSET $TEST_SUITE define=CI_SUPPRESS_KNOWN_ISSUES define=SLOW_COMPILER) notifications: email: From 9e61547f333071a9a612d84a14f16f62715442ee Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Mon, 29 Oct 2018 17:14:46 -0600 Subject: [PATCH 14/20] Fix syntax error in docs. [CI SKIP] --- doc/Jamfile.v2 | 2 +- doc/constants/constants.qbk | 6 ++- doc/distributions/dist_tutorial.qbk | 3 +- doc/interpolators/catmull_rom.qbk | 3 +- doc/overview/faq.qbk | 63 ++++++++++++++++++---------- example/float_comparison_example.cpp | 2 +- 6 files changed, 51 insertions(+), 28 deletions(-) diff --git a/doc/Jamfile.v2 b/doc/Jamfile.v2 index c8b245356..f0874c078 100644 --- a/doc/Jamfile.v2 +++ b/doc/Jamfile.v2 @@ -69,7 +69,7 @@ boostbook standalone pdf:preferred.mediaobject.role=print pdf:img.src.path=$(images_location)/ pdf:draft.mode="no" - pdf:boost.url.prefix=http://www.boost.org/doc/libs/release/libs/math/doc/html + pdf:boost.url.prefix=http\://www.boost.org/doc/libs/release/libs/math/doc/html on pdf:off html:on $(here)/index.idx diff --git a/doc/constants/constants.qbk b/doc/constants/constants.qbk index 99db33659..a2eece24a 100644 --- a/doc/constants/constants.qbk +++ b/doc/constants/constants.qbk @@ -23,7 +23,8 @@ see [link math_toolkit.tutorial.templ use in template code]. * Tested - by comparison with other published sources, or separately computed at long double precision. * Faster - can avoid (re-)calculation at runtime. * If the value returned is a builtin type then it's returned by value as a `constexpr` (C++11 feature, if available). - * If the value is computed and cached (or constructed from a string representation and cached), then it's returned by constant reference.[br] + * If the value is computed and cached (or constructed from a string representation and cached), then it's returned by constant reference. + This can be significant if: * Functions pow, trig or log are used. * Inside an inner loop. @@ -508,7 +509,8 @@ Some of the criteria we have used are: * Expensive to compute. * Requested by users. * [@http://en.wikipedia.org/wiki/Mathematical_constant Used in science and mathematics.] -* No integer values (because so cheap to construct).[br] +* No integer values (because so cheap to construct). + (You can easily define your own if found convenient, for example: `FPT one =static_cast(42);`). [h4 How are constants named?] diff --git a/doc/distributions/dist_tutorial.qbk b/doc/distributions/dist_tutorial.qbk index 0d0e0dfd7..941794e52 100644 --- a/doc/distributions/dist_tutorial.qbk +++ b/doc/distributions/dist_tutorial.qbk @@ -168,7 +168,8 @@ and [@http://en.wikipedia.org/wiki/Parameter distribution parameters] are conventionally distinguished (for example in Wikipedia and Wolfram MathWorld) by placing a semi-colon or vertical bar) /after/ the __random_variable (whose value you 'choose'), -to separate the variate from the parameter(s) that defines the shape of the distribution.[br] +to separate the variate from the parameter(s) that defines the shape of the distribution. + For example, the binomial distribution probability distribution function (PDF) is written as ['f(k| n, p)] = Pr(K = k|n, p) = probability of observing k successes out of n trials. K is the __random_variable, k is the __random_variate, diff --git a/doc/interpolators/catmull_rom.qbk b/doc/interpolators/catmull_rom.qbk index a079f4acb..e27f325f5 100644 --- a/doc/interpolators/catmull_rom.qbk +++ b/doc/interpolators/catmull_rom.qbk @@ -70,7 +70,7 @@ The spline can be either open or *closed*, closed meaning that there is some /t/ The default is open, but this can be easily changed: // closed = true - catmull_rom, 3> cr(points.data(), points.size(), true); + catmull_rom> cr(points.data(), points.size(), true); Inside `catmull_rom`, the Catmull-Rom curve is represented as closed. This is because an open Catmull-Rom curve is *implicitly closed*, but the closing point is the zero vector. @@ -207,6 +207,5 @@ These conditions are satisfied by both `std::array` and `std::vector`, but it ma * Cem Yuksel, Scott Schaefer, and John Keyser, ['Parameterization and applications of Catmull–Rom curves], Computer-Aided Design 43 (2011) 747–755. * Phillip J. Barry and Ronald N. Goldman, ['A Recursive Evaluation Algorithm for a Class of Catmull-Rom Splines], Computer Graphics, Volume 22, Number 4, August 1988 -[endsect] [endsect] [/section:catmull_rom Catmull-Rom Splines] diff --git a/doc/overview/faq.qbk b/doc/overview/faq.qbk index fd7e46b8e..32973d16f 100644 --- a/doc/overview/faq.qbk +++ b/doc/overview/faq.qbk @@ -1,26 +1,32 @@ [section:main_faq Frequently Asked Questions FAQ] # ['I'm a FORTRAN/NAG/SPSS/SAS/Cephes/MathCad/R user -and I don't see where the functions like dnorm(mean, sd) are in Boost.Math?] [br] +and I don't see where the functions like dnorm(mean, sd) are in Boost.Math?] + Nearly all are provided, and many more like mean, skewness, quantiles, complements ... but Boost.Math makes full use of C++, and it looks a bit different. But do not panic! See section on construction and the many examples. Briefly, the distribution is constructed with the parameters (like location and scale) (things after the | in representation like P(X=k|n, p) or ; in a common represention of pdf f(x; [mu][sigma][super 2]). Functions like pdf, cdf are called with the name of that distribution and the random variate often called x or k. -For example, `normal my_norm(0, 1); pdf(my_norm, 2.0);` [br] -#I'm a user of [@http://support.sas.com/rnd/app/da/new/probabilityfunctions.html New SAS Functions for Computing Probabilities]. [br] +For example, `normal my_norm(0, 1); pdf(my_norm, 2.0);` + +#I'm a user of [@http://support.sas.com/rnd/app/da/new/probabilityfunctions.html New SAS Functions for Computing Probabilities]. + You will find the interface more familar, but to be able to select a distribution (perhaps using a string) see the Extras/Future Directions section, and /boost/libs/math/dot_net_example/boost_math.cpp for an example that is used to create a C# (C sharp) utility (that you might also find useful): -see [@http://sourceforge.net/projects/distexplorer/ Statistical Distribution Explorer].[br] -# ['I'm allegic to reading manuals and prefer to learn from examples.][br] +see [@http://sourceforge.net/projects/distexplorer/ Statistical Distribution Explorer]. + +# ['I'm allegic to reading manuals and prefer to learn from examples.] + Fear not - you are not alone! Many examples are available for functions and distributions. Some are referenced directly from the text. Others can be found at \boost_latest_release\libs\math\example. If you are a Visual Studio user, you should be able to create projects from each of these, making sure that the Boost library is in the include directories list. -# ['How do I make sure that the Boost library is in the Visual Studio include directories list?][br] +# ['How do I make sure that the Boost library is in the Visual Studio include directories list?] + You can add an include path, for example, your Boost place /boost-latest_release, for example `X:/boost_1_45_0/` if you have a separate partition X for Boost releases. Or you can use an environment variable BOOST_ROOT set to your Boost place, and include that. @@ -29,32 +35,42 @@ Visual Studio 2010 instead provides property sheets to assist. You may find it convenient to create a new one adding \boost-latest_release; to the existing include items in $(IncludePath). # ['I'm a FORTRAN/NAG/SPSS/SAS/Cephes/MathCad/R user and -I don't see where the properties like mean, median, mode, variance, skewness of distributions are in Boost.Math?][br] +I don't see where the properties like mean, median, mode, variance, skewness of distributions are in Boost.Math?] + They are all available (if defined for the parameters with which you constructed the distribution) via __usual_accessors. -# ['I am a C programmer. Can I user Boost.Math with C?][br] +# ['I am a C programmer. Can I user Boost.Math with C?] + Yes you can, including all the special functions, and TR1 functions like isnan. They appear as C functions, by being declared as "extern C". -# ['I am a C# (Basic? F# FORTRAN? Other CLI?) programmer. Can I use Boost.Math with C#? (or ...)?] [br] +# ['I am a C# (Basic? F# FORTRAN? Other CLI?) programmer. Can I use Boost.Math with C#? (or ...)?] + Yes you can, including all the special functions, and TR1 functions like isnan. But you [*must build the Boost.Math as a dynamic library (.dll) and compile with the /CLI option]. See the boost/math/dot_net_example folder which contains an example that builds a simple statistical distribution app with a GUI. -See [@http://sourceforge.net/projects/distexplorer/ Statistical Distribution Explorer] [br] -# ['What these "policies" things for?] [br] +See [@http://sourceforge.net/projects/distexplorer/ Statistical Distribution Explorer] + +# ['What these "policies" things for?] + Policies are a powerful (if necessarily complex) fine-grain mechanism that allow you to customise the behaviour of the Boost.Math library according to your precise needs. See __policy_section. But if, very probably, the default behaviour suits you, you don't need to know more. -# ['I am a C user and expect to see global C-style`::errno` set for overflow/errors etc?] [br] +# ['I am a C user and expect to see global C-style`::errno` set for overflow/errors etc?] + You can achieve what you want - see __error_policy and __user_error_handling and many examples. -# ['I am a C user and expect to silently return a max value for overflow?] [br] +# ['I am a C user and expect to silently return a max value for overflow?] + You (and C++ users too) can return whatever you want on overflow - see __overflow_error and __error_policy and several examples. -# ['I don't want any error message for overflow etc?] [br] +# ['I don't want any error message for overflow etc?] + You can control exactly what happens for all the abnormal conditions, including the values returned. See __domain_error, __overflow_error __error_policy __user_error_handling etc and examples. -# ['My environment doesn't allow and/or I don't want exceptions. Can I still user Boost.Math?] [br] +# ['My environment doesn't allow and/or I don't want exceptions. Can I still user Boost.Math?] + Yes but you must customise the error handling: see __user_error_handling and __changing_policy_defaults . -# ['The docs are several hundreds of pages long! Can I read the docs off-line or on paper?] [br] +# ['The docs are several hundreds of pages long! Can I read the docs off-line or on paper?] + Yes - you can download the Boost current release of most documentation as a zip of pdfs (including Boost.Math) from Sourceforge, for example [@https://sourceforge.net/projects/boost/files/boost-docs/1.45.0/boost_pdf_1_45_0.tar.gz/download]. @@ -62,7 +78,8 @@ And you can print any pages you need (or even print all pages - but be warned th Both html and pdf versions are highly hyperlinked. The entire Boost.Math pdf can be searched with Adobe Reader, Edit, Find ... This can often find what you seek, a partial substitute for a full index. -# ['I want a compact version for an embedded application. Can I use float precision?] [br] +# ['I want a compact version for an embedded application. Can I use float precision?] + Yes - by selecting RealType template parameter as float: for example normal_distribution your_normal(mean, sd); (But double may still be used internally, so space saving may be less that you hope for). @@ -77,14 +94,18 @@ Probably, thought not always, and not by too much: our priority is accuracy. For most functions, making sure you have the latest compiler version with all optimisations switched on is the key to speed. For evaluations that require iteration, you may be able to gain a little more speed at the expense of accuracy. See detailed suggestions and results on __performance. -# ['How do I handle infinity and NaNs portably?] [br] +# ['How do I handle infinity and NaNs portably?] + See __fp_facets for Facets for Floating-Point Infinities and NaNs. -# ['Where are the pre-built libraries?] [br] +# ['Where are the pre-built libraries?] + Good news - you probably don't need any! - just `#include ]. But in the unlikely event that you do, see __building. -# ['I don't see the function or distribution that I want.] [br] +# ['I don't see the function or distribution that I want.] + You could try an email to ask the authors - but no promises! -# ['I need more decimal digits for values/computations.] [br] +# ['I need more decimal digits for values/computations.] + You can use Boost.Math with __multiprecision: typically __cpp_dec_float is a useful user-defined type to provide a fixed number of decimal digits, usually 50 or 100. # Why can't I write something really simple like `cpp_int one(1); cpp_dec_float_50 two(2); one * two;` diff --git a/example/float_comparison_example.cpp b/example/float_comparison_example.cpp index f39e53781..6c892fa21 100644 --- a/example/float_comparison_example.cpp +++ b/example/float_comparison_example.cpp @@ -45,7 +45,7 @@ for this purpose, and here we use `float` precision where `max_digits10` = 9 to avoid displaying a distracting number of decimal digits. [note Older compilers can use this formula to calculate `max_digits10` -from `std::numeric_limits::digits10`:[br] +from `std::numeric_limits::digits10`: __spaces `int max_digits10 = 2 + std::numeric_limits::digits10 * 3010/10000;` ] [/note] From d5c9c94d1403c2c9a6a4e98f48025c303877d971 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Tue, 30 Oct 2018 11:40:55 -0600 Subject: [PATCH 15/20] Revert travis.yml now that it's been demonstrated that Catmull-Rom doesn't cause this error. [CI SKIP] --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 8b1177b8a..fc5db3b0c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -584,7 +584,7 @@ script: - |- echo "using $TOOLSET : : $COMPILER : -std=$CXXSTD ;" > ~/user-config.jam - (cd libs/config/test && ../../../b2 config_info_travis_install toolset=$TOOLSET && ./config_info_travis) - - (cd libs/math/test && travis_wait 60 ../../../b2 -j3 -d2 -q toolset=$TOOLSET $TEST_SUITE define=CI_SUPPRESS_KNOWN_ISSUES define=SLOW_COMPILER) + - (cd libs/math/test && travis_wait 60 ../../../b2 -j3 -d+0 -q toolset=$TOOLSET $TEST_SUITE define=CI_SUPPRESS_KNOWN_ISSUES define=SLOW_COMPILER) notifications: email: From fd519a73d6a4e9a438f98c46f344a4580d0ea03e Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sun, 2 Dec 2018 11:18:41 -0700 Subject: [PATCH 16/20] Implement suggestions from code review [CI SKIP] --- doc/interpolators/catmull_rom.qbk | 30 ++++++++++++------- .../boost/math/interpolators/catmull_rom.hpp | 11 +++++-- 2 files changed, 28 insertions(+), 13 deletions(-) diff --git a/doc/interpolators/catmull_rom.qbk b/doc/interpolators/catmull_rom.qbk index e27f325f5..3f4771da9 100644 --- a/doc/interpolators/catmull_rom.qbk +++ b/doc/interpolators/catmull_rom.qbk @@ -42,13 +42,15 @@ Catmull-Rom splines enjoy the following properties: * Affine invariance: The interpolant commutes with affine transformations. * Local support of the basis functions: This gives stability and fast evaluation. -* Smoothness -* Interpolation of control points: Many curves (such as Bezier) are *approximating*, they do not pass through their control points. This makes them more difficult to use than interpolating splines. +* /C/[super 2] smoothness +* Interpolation of control points-this means the curve passes through the control points. +Many curves (such as Bezier) are /approximating/-they do not pass through their control points. +This makes them more difficult to use than interpolating splines. -The `catmull_rom` class provided by boost creates a cubic Catmull-Rom spline from an array of points in any dimension. +The `catmull_rom` class provided by Boost creates a cubic Catmull-Rom spline from an array of points in any dimension. Since there are numerous ways to represent a point in /n/-dimensional space, the class attempts to be flexible by templating on the point type. -The requirements on the point type are discussing in more detail below, but roughly, it must have a dereference operator defined (i.e., `p[0]` is not a syntax error), +The requirements on the point type are discussing in more detail below, but roughly, it must have a dereference operator defined (e.g., `p[0]` is not a syntax error), it must be able to be dereferenced up to `dimension -1`, and `p[i]` is of type `Real`, define `value_type`, and the free function `size()`. These requirements are met by `std::vector` and `std::array`. The basic usage is shown here: @@ -66,22 +68,26 @@ Memory conscious programmers may enjoy using the move constructor instead: catmull_rom> cr(std::move(points)); -The spline can be either open or *closed*, closed meaning that there is some /t/ such that /P(t) = P(0)/. +The spline can be either open or /closed/, closed meaning that there is some /s > 0/ such that /P(s) = P(0)/. The default is open, but this can be easily changed: // closed = true catmull_rom> cr(points.data(), points.size(), true); -Inside `catmull_rom`, the Catmull-Rom curve is represented as closed. -This is because an open Catmull-Rom curve is *implicitly closed*, but the closing point is the zero vector. +In either case, evaluating the interpolator at /s=0/ returns the first point in the list. + +Inside `catmull_rom`, the curve is represented as closed. +This is because an open Catmull-Rom curve is /implicitly closed/, but the closing point is the zero vector. There is no reason to suppose that the zero vector is a better closing point than the endpoint (or any other point, for that matter), so traditionally Catmull-Rom splines leave the segment between the first and second point undefined, as well as the segment between the second-to-last and last point. We find this property of the traditional implementation of Catmull-Rom splines annoying and confusing to the user. Hence internally, we close the curve so that the first and last segments are defined. -Of course, this causes the *tangent* vectors to the first and last points to be bizarre. +Of course, this causes the /tangent/ vectors to the first and last points to be bizarre. This is a "pick your poison" design decision-either the curve cannot interpolate in its first and last segments, or the tangents along the first and last segments are meaningless. +In the vast majority of cases, this will be no problem to the user. +However, if it becomes a problem, then the user should add one extra point in a position they believe is reasonable and close the curve. Since the routine internally represents the curve as closed, a question arises: Why does the user have to specify if the curve is open or closed? @@ -96,13 +102,14 @@ If the curve is open, then `cr(max_s) = pf`, where `pf` is the final point on th The Catmull-Rom curve admits an infinite number of parameterizations. -The default parameterization of the `catmull_rom` class is the so-called *centripedal* parameterization. +The default parameterization of the `catmull_rom` class is the so-called /centripedal/ parameterization. This parameterization has been shown to be the only parameterization that does not form cusps or self-intersections within segments. However, for advanced users, other parameterizations can be chosen using the /alpha/ parameter: // alpha = 1 is the "chordal" parameterization. catmull_rom> cr(points.data(), points.size(), false, 1.0); +The alpha parameter must always be in the range `[0,1]`. Finally, the tangent vector to any point of the curve can be computed via @@ -110,8 +117,8 @@ Finally, the tangent vector to any point of the curve can be computed via Point tangent = cr.prime(s); Since the magnitude of the tangent vector is dependent on the parameterization, -it is not as meaningful as (say) arc-length parameterization. -However, its direction is meaningful, so the user may wish to normalize this result. +it is not meaningful (unless the user chooses the chordal parameterization /alpha = 1/ which parameterizes by Euclidean distance between points.) +However, its direction is meaningful no matter the parameterization, so the user may wish to normalize this result. [heading Examples] @@ -121,6 +128,7 @@ However, its direction is meaningful, so the user may wish to normalize this res The following performance numbers were generated for a call to the Catmull-Rom interpolation method. The number that follows the slash is the number of points passed to the interpolant. +We see that evaluation of the interpolant is [bigo](/log/(/N/)). Run on 2700 MHz CPU diff --git a/include/boost/math/interpolators/catmull_rom.hpp b/include/boost/math/interpolators/catmull_rom.hpp index dd9435137..0304a871d 100644 --- a/include/boost/math/interpolators/catmull_rom.hpp +++ b/include/boost/math/interpolators/catmull_rom.hpp @@ -42,6 +42,8 @@ public: catmull_rom(std::vector&& points, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2) : catmull_rom(points.data(), points.size(), closed, alpha) {} + catmull_rom(std::initializer_list l, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2) : catmull_rom(std::vector(l), closed, alpha) {} + typename Point::value_type max_parameter() const { return m_max_s; @@ -68,8 +70,13 @@ catmull_rom::catmull_rom(const Point* const points, size_t num_pnts, bool { if (num_pnts < 4) { - throw std::domain_error("The Catmull-Rom curve requires at least 4 points.\n"); + throw std::domain_error("The Catmull-Rom curve requires at least 4 points."); } + if (alpha < 0 || alpha > 1) + { + throw std::domain_error("The parametrization alpha must be in the range [0,1]."); + } + using std::abs; m_s.resize(num_pnts+3); m_pnts.resize(num_pnts+3); @@ -113,7 +120,7 @@ Point catmull_rom::operator()(const typename Point::value_type s) const using std::size; if (s < 0 || s > m_max_s) { - throw std::domain_error("Parameter outside bounds.\n"); + throw std::domain_error("Parameter outside bounds."); } auto it = std::upper_bound(m_s.begin(), m_s.end(), s); //Now *it >= s. We want the index such that m_s[i] <= s < m_s[i+1]: From 5169fc9e75063f7dd71b4e5a51e2761de2e448f8 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sun, 2 Dec 2018 12:04:31 -0700 Subject: [PATCH 17/20] Actually test the initializer list constructor [CI SKIP] --- test/catmull_rom_test.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index 9ad59cd6c..ac068c00a 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -338,8 +338,7 @@ void test_data_representations() mypoint3d p3(0.4, 0.5, 0.6); mypoint3d p4(0.5, 0.6, 0.7); mypoint3d p5(0.6, 0.7, 0.8); - std::vector> v{p0, p1, p2, p3, p4, p5}; - catmull_rom> cat(v.data(), v.size()); + catmull_rom> cat({p0, p1, p2, p3, p4, p5}); Real tol = 0.001; auto p = cat(cat.parameter_at_point(0)); From baddf9509ab8c6f8c55178b0bd4ae1dbec280367 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sun, 2 Dec 2018 13:34:51 -0700 Subject: [PATCH 18/20] The move constructor is 30% faster than the copy; hence remove the data copy and only allow move construction. [CI SKIP] --- doc/interpolators/catmull_rom.qbk | 28 ++++++++++--------- .../boost/math/interpolators/catmull_rom.hpp | 28 ++++++++++++------- test/catmull_rom_test.cpp | 28 +++++++++++-------- 3 files changed, 50 insertions(+), 34 deletions(-) diff --git a/doc/interpolators/catmull_rom.qbk b/doc/interpolators/catmull_rom.qbk index 3f4771da9..8a46314f7 100644 --- a/doc/interpolators/catmull_rom.qbk +++ b/doc/interpolators/catmull_rom.qbk @@ -19,10 +19,11 @@ namespace boost{ namespace math{ class catmull_rom { public: - catmull_rom(const Point* const points, size_t num_points, bool closed = false, Real alpha = (Real) 1/ (Real) 2) catmull_rom(std::vector&& points, bool closed = false, Real alpha = (Real) 1/ (Real) 2) + catmull_rom(std::initializer_list l, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2); + Real operator()(Real s) const; Real max_parameter() const; @@ -42,7 +43,7 @@ Catmull-Rom splines enjoy the following properties: * Affine invariance: The interpolant commutes with affine transformations. * Local support of the basis functions: This gives stability and fast evaluation. -* /C/[super 2] smoothness +* /C/[super 2]-smoothness * Interpolation of control points-this means the curve passes through the control points. Many curves (such as Bezier) are /approximating/-they do not pass through their control points. This makes them more difficult to use than interpolating splines. @@ -60,19 +61,15 @@ The basic usage is shown here: points[1] = {1,0,0}; points[2] = {0,1,0}; points[3] = {0,0,1}; - catmull_rom> cr(points.data(), points.size()); + catmull_rom> cr(std::move(points)); // Interpolate at s = 0.1: auto point = cr(0.1); -Memory conscious programmers may enjoy using the move constructor instead: - - catmull_rom> cr(std::move(points)); - The spline can be either open or /closed/, closed meaning that there is some /s > 0/ such that /P(s) = P(0)/. The default is open, but this can be easily changed: // closed = true - catmull_rom> cr(points.data(), points.size(), true); + catmull_rom> cr(std::move(points), true); In either case, evaluating the interpolator at /s=0/ returns the first point in the list. @@ -107,7 +104,7 @@ This parameterization has been shown to be the only parameterization that does n However, for advanced users, other parameterizations can be chosen using the /alpha/ parameter: // alpha = 1 is the "chordal" parameterization. - catmull_rom> cr(points.data(), points.size(), false, 1.0); + catmull_rom> cr(std::move(points), false, 1.0); The alpha parameter must always be in the range `[0,1]`. @@ -162,9 +159,8 @@ We see that evaluation of the interpolant is [bigo](/log/(/N/)). BM_CatmullRom_BigO 2.97 lgN 2.97 lgN BM_CatmullRom_RMS 19 % 19 % -[endsect] -[section:point_type Point types] +[heading Point types] We have already discussed that certain conditions on the `Point` type template argument must be obeyed. The following shows a custom point type in 3D which can be used as a template argument to Catmull-Rom: @@ -208,9 +204,15 @@ The following shows a custom point type in 3D which can be used as a template ar These conditions are satisfied by both `std::array` and `std::vector`, but it may nonetheless be useful to define your own point class so that (say) you can define geometric distance between them. -[endsect] -[section:catmull_rom_refs References] +[heading Caveats] + +The Catmull-Rom interpolator requires memory for three more points than is provided by the user. +This causes the class to call a `resize()` on the input vector. +If `v.capacity() >= v.size() + 3`, then no problems arise; there are no reallocs, and in practice this condition is almost always satisfied. +However, if `v.capacity() < v.size() + 3`, the realloc causes a performance penalty of roughly 20%. + +[heading References] * Cem Yuksel, Scott Schaefer, and John Keyser, ['Parameterization and applications of Catmull–Rom curves], Computer-Aided Design 43 (2011) 747–755. * Phillip J. Barry and Ronald N. Goldman, ['A Recursive Evaluation Algorithm for a Class of Catmull-Rom Splines], Computer Graphics, Volume 22, Number 4, August 1988 diff --git a/include/boost/math/interpolators/catmull_rom.hpp b/include/boost/math/interpolators/catmull_rom.hpp index 0304a871d..05ad666de 100644 --- a/include/boost/math/interpolators/catmull_rom.hpp +++ b/include/boost/math/interpolators/catmull_rom.hpp @@ -38,9 +38,7 @@ class catmull_rom { public: - catmull_rom(const Point* const points, size_t num_pnts, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2); - - catmull_rom(std::vector&& points, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2) : catmull_rom(points.data(), points.size(), closed, alpha) {} + catmull_rom(std::vector&& points, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2); catmull_rom(std::initializer_list l, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2) : catmull_rom(std::vector(l), closed, alpha) {} @@ -56,9 +54,13 @@ public: Point operator()(const typename Point::value_type s) const; - Point prime(const typename Point::value_type s) const; + std::vector&& get_points() + { + return std::move(m_pnts); + } + private: std::vector m_pnts; std::vector m_s; @@ -66,8 +68,10 @@ private: }; template -catmull_rom::catmull_rom(const Point* const points, size_t num_pnts, bool closed, typename Point::value_type alpha) +catmull_rom::catmull_rom(std::vector&& points, bool closed, typename Point::value_type alpha) : m_pnts(std::move(points)) { + size_t num_pnts = m_pnts.size(); + //std::cout << "Number of points = " << num_pnts << "\n"; if (num_pnts < 4) { throw std::domain_error("The Catmull-Rom curve requires at least 4 points."); @@ -80,14 +84,18 @@ catmull_rom::catmull_rom(const Point* const points, size_t num_pnts, bool using std::abs; m_s.resize(num_pnts+3); m_pnts.resize(num_pnts+3); + //std::cout << "Number of points now = " << m_pnts.size() << "\n"; - m_pnts[0] = points[num_pnts-1]; - for (size_t i = 0; i < num_pnts; ++i) + m_pnts[num_pnts+1] = m_pnts[0]; + m_pnts[num_pnts+2] = m_pnts[1]; + + auto tmp = m_pnts[num_pnts-1]; + for (int64_t i = num_pnts-1; i >= 0; --i) { - m_pnts[i+1] = points[i]; + m_pnts[i+1] = m_pnts[i]; } - m_pnts[num_pnts+1] = points[0]; - m_pnts[num_pnts+2] = points[1]; + m_pnts[0] = tmp; + m_s[0] = -detail::alpha_distance(m_pnts[0], m_pnts[1], alpha); if (abs(m_s[0]) < std::numeric_limits::epsilon()) { diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index ac068c00a..bdd205b6b 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -60,7 +60,7 @@ void test_linear() v[1] = {1,0,0}; v[2] = {2,0,0}; v[3] = {3,0,0}; - catmull_rom> cr(v.data(), v.size()); + catmull_rom> cr(std::move(v)); // Test that the interpolation condition is obeyed: BOOST_CHECK_CLOSE_FRACTION(cr.max_parameter(), 3, tol); @@ -120,12 +120,14 @@ void test_circle() Real tol = 10*std::numeric_limits::epsilon(); std::vector> v(20*sizeof(Real)); + std::vector> u(20*sizeof(Real)); for (size_t i = 0; i < v.size(); ++i) { Real theta = ((Real) i/ (Real) v.size())*2*M_PI; v[i] = {cos(theta), sin(theta)}; + u[i] = v[i]; } - catmull_rom> circle(v.data(), v.size(), true); + catmull_rom> circle(std::move(v), true); // Interpolation condition: for (size_t i = 0; i < v.size(); ++i) @@ -136,16 +138,16 @@ void test_circle() Real y = p[1]; if (abs(x) < std::numeric_limits::epsilon()) { - BOOST_CHECK_SMALL(v[i][0], tol); + BOOST_CHECK_SMALL(u[i][0], tol); } if (abs(y) < std::numeric_limits::epsilon()) { - BOOST_CHECK_SMALL(v[i][1], tol); + BOOST_CHECK_SMALL(u[i][1], tol); } else { - BOOST_CHECK_CLOSE_FRACTION(x, v[i][0], tol); - BOOST_CHECK_CLOSE_FRACTION(y, v[i][1], tol); + BOOST_CHECK_CLOSE_FRACTION(x, u[i][0], tol); + BOOST_CHECK_CLOSE_FRACTION(y, u[i][1], tol); } } @@ -169,11 +171,13 @@ void test_affine_invariance() Real tol = 1000*std::numeric_limits::epsilon(); std::vector> v(100); + std::vector> u(100); std::mt19937_64 gen(438232); Real inv_denom = (Real) 100/( (Real) gen.max() + (Real) 2); for(size_t j = 0; j < dimension; ++j) { v[0][j] = gen()*inv_denom; + u[0][j] = v[0][j]; } for (size_t i = 1; i < v.size(); ++i) @@ -181,6 +185,7 @@ void test_affine_invariance() for(size_t j = 0; j < dimension; ++j) { v[i][j] = v[i-1][j] + gen()*inv_denom; + u[i][j] = v[i][j]; } } std::array affine_shift; @@ -189,17 +194,17 @@ void test_affine_invariance() affine_shift[j] = gen()*inv_denom; } - catmull_rom> cr1(v.data(), v.size()); + catmull_rom> cr1(std::move(v)); - for(size_t i = 0; i< v.size(); ++i) + for(size_t i = 0; i< u.size(); ++i) { for(size_t j = 0; j < dimension; ++j) { - v[i][j] += affine_shift[j]; + u[i][j] += affine_shift[j]; } } - catmull_rom> cr2(v.data(), v.size()); + catmull_rom> cr2(std::move(u)); BOOST_CHECK_CLOSE_FRACTION(cr1.max_parameter(), cr2.max_parameter(), tol); @@ -234,7 +239,7 @@ void test_helix() Real theta = ((Real) i/ (Real) v.size())*2*M_PI; v[i] = {cos(theta), sin(theta), theta}; } - catmull_rom> helix(v.data(), v.size()); + catmull_rom> helix(std::move(v)); // Interpolation condition: for (size_t i = 0; i < v.size(); ++i) @@ -338,6 +343,7 @@ void test_data_representations() mypoint3d p3(0.4, 0.5, 0.6); mypoint3d p4(0.5, 0.6, 0.7); mypoint3d p5(0.6, 0.7, 0.8); + // Tests initializer_list: catmull_rom> cat({p0, p1, p2, p3, p4, p5}); Real tol = 0.001; From 081842118bd4cf270b48c06fb66db85d400d3c42 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Mon, 3 Dec 2018 12:33:45 -0700 Subject: [PATCH 19/20] Give advice about first and last interpolator segments in the open-curve case. [CI SKIP] --- doc/interpolators/catmull_rom.qbk | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/doc/interpolators/catmull_rom.qbk b/doc/interpolators/catmull_rom.qbk index 8a46314f7..6d4f4e35a 100644 --- a/doc/interpolators/catmull_rom.qbk +++ b/doc/interpolators/catmull_rom.qbk @@ -73,6 +73,11 @@ The default is open, but this can be easily changed: In either case, evaluating the interpolator at /s=0/ returns the first point in the list. +If the curve is open, then the first and last segments may have strange behavior. +The traditional solution is to prepend a carefully selected control point to the data so that the first data segment (second interpolator segment) has reasonable tangent vectors, and simply ignore the first interpolator segment. +A control point is appended to the data using similar criteria. +However, we recommend not going through this effort until it proves to be necessary: For most use-cases, the curve is good enough without prepending and appending control points, and responsible selection of non-data control points is difficult. + Inside `catmull_rom`, the curve is represented as closed. This is because an open Catmull-Rom curve is /implicitly closed/, but the closing point is the zero vector. There is no reason to suppose that the zero vector is a better closing point than the endpoint (or any other point, for that matter), From 8c6740463a82866beaad1a1fb57e25efd6df116e Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Mon, 3 Dec 2018 17:06:21 -0700 Subject: [PATCH 20/20] Meaningless commit to kick off CI build. --- test/catmull_rom_test.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/test/catmull_rom_test.cpp b/test/catmull_rom_test.cpp index bdd205b6b..a95ff5692 100644 --- a/test/catmull_rom_test.cpp +++ b/test/catmull_rom_test.cpp @@ -375,7 +375,6 @@ BOOST_AUTO_TEST_CASE(catmull_rom_test) test_circle(); test_circle(); - test_helix(); test_helix();