math/interpolators/cubic_b_spline.hpp

0001 // Copyright Nick Thompson, 2017
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 // This implements the compactly supported cubic b spline algorithm described in
0008 // Kress, Rainer. "Numerical analysis, volume 181 of Graduate Texts in Mathematics." (1998).
0009 // Splines of compact support are faster to evaluate and are better conditioned than classical cubic splines.
0010
0011 // Let f be the function we are trying to interpolate, and s be the interpolating spline.
0012 // The routine constructs the interpolant in O(N) time, and evaluating s at a point takes constant time.
0013 // The order of accuracy depends on the regularity of the f, however, assuming f is
0014 // four-times continuously differentiable, the error is of O(h^4).
0015 // In addition, we can differentiate the spline and obtain a good interpolant for f'.
0016 // The main restriction of this method is that the samples of f must be evenly spaced.
0017 // Look for barycentric rational interpolation for non-evenly sampled data.
0018 // Properties:
0019 // - s(x_j) = f(x_j)
0020 // - All cubic polynomials interpolated exactly
0021
0022 #ifndef BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
0023 #define BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP
0024
0025 #include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
0026 #include <boost/math/tools/header_deprecated.hpp>
0027
0028 BOOST_MATH_HEADER_DEPRECATED("<boost/math/interpolators/cardinal_cubic_b_spline.hpp>");
0029
0030 namespace boost{ namespace math{
0031
0032 template <class Real>
0033 class cubic_b_spline
0034 {
0035 public:
0036     // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
0037     // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
0038     template <class BidiIterator>
0039     cubic_b_spline(const BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
0040                    Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
0041                    Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
0042     cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
0043        Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
0044        Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
0045
0046     cubic_b_spline() = default;
0047     Real operator()(Real x) const;
0048
0049     Real prime(Real x) const;
0050
0051     Real double_prime(Real x) const;
0052
0053 private:
0054     std::shared_ptr<detail::cubic_b_spline_imp<Real>> m_imp;
0055 };
0056
0057 template<class Real>
0058 cubic_b_spline<Real>::cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
0059                                      Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, f + length, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
0060 {
0061 }
0062
0063 template <class Real>
0064 template <class BidiIterator>
0065 cubic_b_spline<Real>::cubic_b_spline(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
0066    Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cubic_b_spline_imp<Real>>(f, end_p, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
0067 {
0068 }
0069
0070 template<class Real>
0071 Real cubic_b_spline<Real>::operator()(Real x) const
0072 {
0073     return m_imp->operator()(x);
0074 }
0075
0076 template<class Real>
0077 Real cubic_b_spline<Real>::prime(Real x) const
0078 {
0079     return m_imp->prime(x);
0080 }
0081
0082 template<class Real>
0083 Real cubic_b_spline<Real>::double_prime(Real x) const
0084 {
0085     return m_imp->double_prime(x);
0086 }
0087
0088
0089 }}
0090 #endif