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
7 // This implements the compactly supported cubic b spline algorithm described in
8 // 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 is
14 // 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 exactly
21
22 #ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_CUBIC_B_SPLINE_HPP
23 #define BOOST_MATH_INTERPOLATORS_CARDINAL_CUBIC_B_SPLINE_HPP
24
25 #include <boost/math/interpolators/detail/cardinal_cubic_b_spline_detail.hpp>
26
27 namespace boost{ namespace math{ namespace interpolators {
28
29 template <class Real>
30 class cardinal_cubic_b_spline
31 {
32 public:
33 // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
34 // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
35 template <class BidiIterator>
36 cardinal_cubic_b_spline(const BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
37 Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
38 Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
39 cardinal_cubic_b_spline(const Real* const f, size_t length, 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
43 cardinal_cubic_b_spline() = default;
44 Real operator()(Real x) const;
45
46 Real prime(Real x) const;
47
48 Real double_prime(Real x) const;
49
50 private:
51 std::shared_ptr<detail::cardinal_cubic_b_spline_imp<Real>> m_imp;
52 };
53
54 template<class Real>
cardinal_cubic_b_spline(const Real * const f,size_t length,Real left_endpoint,Real step_size,Real left_endpoint_derivative,Real right_endpoint_derivative)55 cardinal_cubic_b_spline<Real>::cardinal_cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size,
56 Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cardinal_cubic_b_spline_imp<Real>>(f, f + length, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
57 {
58 }
59
60 template <class Real>
61 template <class BidiIterator>
cardinal_cubic_b_spline(BidiIterator f,BidiIterator end_p,Real left_endpoint,Real step_size,Real left_endpoint_derivative,Real right_endpoint_derivative)62 cardinal_cubic_b_spline<Real>::cardinal_cubic_b_spline(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
63 Real left_endpoint_derivative, Real right_endpoint_derivative) : m_imp(std::make_shared<detail::cardinal_cubic_b_spline_imp<Real>>(f, end_p, left_endpoint, step_size, left_endpoint_derivative, right_endpoint_derivative))
64 {
65 }
66
67 template<class Real>
operator ()(Real x) const68 Real cardinal_cubic_b_spline<Real>::operator()(Real x) const
69 {
70 return m_imp->operator()(x);
71 }
72
73 template<class Real>
prime(Real x) const74 Real cardinal_cubic_b_spline<Real>::prime(Real x) const
75 {
76 return m_imp->prime(x);
77 }
78
79 template<class Real>
double_prime(Real x) const80 Real cardinal_cubic_b_spline<Real>::double_prime(Real x) const
81 {
82 return m_imp->double_prime(x);
83 }
84
85
86 }}}
87 #endif
88