91 lines · plain
1// Copyright Nick Thompson, 20172// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0.4// (See accompanying file LICENSE_1_0.txt5// or copy at http://www.boost.org/LICENSE_1_0.txt)6 7// This implements the compactly supported cubic b spline algorithm described in8// Kress, Rainer. "Numerical analysis, volume 181 of Graduate Texts in Mathematics." (1998).9// Splines of compact support are faster to evaluate and are better conditioned than classical cubic splines.10 11// Let f be the function we are trying to interpolate, and s be the interpolating spline.12// The routine constructs the interpolant in O(N) time, and evaluating s at a point takes constant time.13// The order of accuracy depends on the regularity of the f, however, assuming f is14// four-times continuously differentiable, the error is of O(h^4).15// In addition, we can differentiate the spline and obtain a good interpolant for f'.16// The main restriction of this method is that the samples of f must be evenly spaced.17// Look for barycentric rational interpolation for non-evenly sampled data.18// Properties:19// - s(x_j) = f(x_j)20// - All cubic polynomials interpolated exactly21 22#ifndef BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP23#define BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP24 25#include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>26#include <boost/math/tools/header_deprecated.hpp>27 28BOOST_MATH_HEADER_DEPRECATED("<boost/math/interpolators/cardinal_cubic_b_spline.hpp>");29 30namespace boost{ namespace math{31 32template <class Real>33class cubic_b_spline34{35public:36 // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.37 // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).38 template <class BidiIterator>39 cubic_b_spline(const BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,40 Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),41 Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());42 cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,43 Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),44 Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());45 46 cubic_b_spline() = default;47 Real operator()(Real x) const;48 49 Real prime(Real x) const;50 51 Real double_prime(Real x) const;52 53private:54 std::shared_ptr<detail::cubic_b_spline_imp<Real>> m_imp;55};56 57template<class Real>58cubic_b_spline<Real>::cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,59 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))60{61}62 63template <class Real>64template <class BidiIterator>65cubic_b_spline<Real>::cubic_b_spline(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,66 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))67{68}69 70template<class Real>71Real cubic_b_spline<Real>::operator()(Real x) const72{73 return m_imp->operator()(x);74}75 76template<class Real>77Real cubic_b_spline<Real>::prime(Real x) const78{79 return m_imp->prime(x);80}81 82template<class Real>83Real cubic_b_spline<Real>::double_prime(Real x) const84{85 return m_imp->double_prime(x);86}87 88 89}}90#endif91