First pass at a Catmull-Rom curve interpolator.

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 <boost/math/interpolators/catmull_rom.hpp>
+
+namespace boost{ namespace math{
+    template<class Real, class Point>
+    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<std::array<Real, 3>> points(4);
+    points[0] = {0,0,0};
+    points[1] = {1,0,0};
+    points[2] = {0,1,0};
+    points[3] = {0,0,1};
+    catmull_rom<Real, std::array<Real, 3>, 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<Real, std::array<Real, 3>, 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<Real, std::array<Real, 3>, 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]

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 <iostream>
+#include <boost/math/interpolators/catmull_rom.hpp>
+
+int main()
+{
+    std::cout << "This shows how to use the Catmull-Rom class of Boost to create a logarithmic spiral.\n";
+    return 0;
+}

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 <cmath>
+#include <algorithm>
+
+namespace boost{ namespace math{
+
+    namespace detail
+    {
+        template<class Real, class Point, size_t dimension>
+        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 Real, class Point, size_t dimension>
+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<Real, Point, dimension>(m_pnts[0], m_pnts[1], alpha);
+            if (abs(m_s[0]) < std::numeric_limits<Real>::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<Real, Point, dimension>(m_pnts[i], m_pnts[i-1], alpha);
+                if (abs(d) < std::numeric_limits<Real>::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<Real, Point, dimension>(m_pnts[i-1], m_pnts[i-2], alpha);
+                if (abs(d) < std::numeric_limits<Real>::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<Point> m_pnts;
+    std::vector<Real> m_s;
+    Real m_max_s;
+};
+
+}}
+#endif

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 <array>
+#include <boost/cstdfloat.hpp>
+#include <boost/type_index.hpp>
+#include <boost/test/included/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+#include <boost/math/interpolators/catmull_rom.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+
+using boost::math::catmull_rom;
+
+template<class Real>
+void test_linear()
+{
+    Real tol = std::numeric_limits<Real>::epsilon();
+    std::vector<std::array<Real, 3>> v(4);
+    v[0] = {0,0,0};
+    v[1] = {1,0,0};
+    v[2] = {2,0,0};
+    v[3] = {3,0,0};
+    catmull_rom<Real, std::array<Real, 3>, 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<class Real>
+void test_affine_invariance()
+{
+    std::cout << "Testing that the Catmull-Rom spline is affine invariant.\n";
+
+}
+
+template<class Real>
+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<float>();
+    test_affine_invariance<float>();
+    test_data_representations<float>();
+}