cubic_b_spline.hpp 3.6 KB

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