interpolators/detail/cubic_b_spline_detail.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 #ifndef CUBIC_B_SPLINE_DETAIL_HPP
0008 #define CUBIC_B_SPLINE_DETAIL_HPP
0009
0010 #include <limits>
0011 #include <cmath>
0012 #include <vector>
0013 #include <memory>
0014 #include <boost/math/constants/constants.hpp>
0015 #include <boost/math/special_functions/fpclassify.hpp>
0016
0017 namespace boost{ namespace math{ namespace detail{
0018
0019
0020 template <class Real>
0021 class cubic_b_spline_imp
0022 {
0023 public:
0024     // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
0025     // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
0026     template <class BidiIterator>
0027     cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
0028                        Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
0029                        Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
0030
0031     Real operator()(Real x) const;
0032
0033     Real prime(Real x) const;
0034
0035     Real double_prime(Real x) const;
0036
0037 private:
0038     std::vector<Real> m_beta;
0039     Real m_h_inv;
0040     Real m_a;
0041     Real m_avg;
0042 };
0043
0044
0045
0046 template <class Real>
0047 Real b3_spline(Real x)
0048 {
0049     using std::abs;
0050     Real absx = abs(x);
0051     if (absx < 1)
0052     {
0053         Real y = 2 - absx;
0054         Real z = 1 - absx;
0055         return boost::math::constants::sixth<Real>()*(y*y*y - 4*z*z*z);
0056     }
0057     if (absx < 2)
0058     {
0059         Real y = 2 - absx;
0060         return boost::math::constants::sixth<Real>()*y*y*y;
0061     }
0062     return static_cast<Real>(0);
0063 }
0064
0065 template<class Real>
0066 Real b3_spline_prime(Real x)
0067 {
0068     if (x < 0)
0069     {
0070         return -b3_spline_prime(-x);
0071     }
0072
0073     if (x < 1)
0074     {
0075         return x*(3*boost::math::constants::half<Real>()*x - 2);
0076     }
0077     if (x < 2)
0078     {
0079         return -boost::math::constants::half<Real>()*(2 - x)*(2 - x);
0080     }
0081     return static_cast<Real>(0);
0082 }
0083
0084 template<class Real>
0085 Real b3_spline_double_prime(Real x)
0086 {
0087     if (x < 0)
0088     {
0089         return b3_spline_double_prime(-x);
0090     }
0091
0092     if (x < 1)
0093     {
0094         return 3*x - 2;
0095     }
0096     if (x < 2)
0097     {
0098         return (2 - x);
0099     }
0100     return static_cast<Real>(0);
0101 }
0102
0103
0104 template <class Real>
0105 template <class BidiIterator>
0106 cubic_b_spline_imp<Real>::cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
0107                                              Real left_endpoint_derivative, Real right_endpoint_derivative) : m_a(left_endpoint), m_avg(0)
0108 {
0109     using boost::math::constants::third;
0110
0111     std::size_t length = end_p - f;
0112
0113     if (length < 5)
0114     {
0115         if (boost::math::isnan(left_endpoint_derivative) || boost::math::isnan(right_endpoint_derivative))
0116         {
0117             throw std::logic_error("Interpolation using a cubic b spline with derivatives estimated at the endpoints requires at least 5 points.\n");
0118         }
0119         if (length < 3)
0120         {
0121             throw std::logic_error("Interpolation using a cubic b spline requires at least 3 points.\n");
0122         }
0123     }
0124
0125     if (boost::math::isnan(left_endpoint))
0126     {
0127         throw std::logic_error("Left endpoint is NAN; this is disallowed.\n");
0128     }
0129     if (left_endpoint + length*step_size >= (std::numeric_limits<Real>::max)())
0130     {
0131         throw std::logic_error("Right endpoint overflows the maximum representable number of the specified precision.\n");
0132     }
0133     if (step_size <= 0)
0134     {
0135         throw std::logic_error("The step size must be strictly > 0.\n");
0136     }
0137     // Storing the inverse of the stepsize does provide a measurable speedup.
0138     // It's not huge, but nonetheless worthwhile.
0139     m_h_inv = 1/step_size;
0140
0141     // Following Kress's notation, s'(a) = a1, s'(b) = b1
0142     Real a1 = left_endpoint_derivative;
0143     // See the finite-difference table on Wikipedia for reference on how
0144     // to construct high-order estimates for one-sided derivatives:
0145     // https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_and_backward_finite_difference
0146     // Here, we estimate then to O(h^4), as that is the maximum accuracy we could obtain from this method.
0147     if (boost::math::isnan(a1))
0148     {
0149         // For simple functions (linear, quadratic, so on)
0150         // almost all the error comes from derivative estimation.
0151         // This does pairwise summation which gives us another digit of accuracy over naive summation.
0152         Real t0 = 4*(f[1] + third<Real>()*f[3]);
0153         Real t1 = -(25*third<Real>()*f[0] + f[4])/4  - 3*f[2];
0154         a1 = m_h_inv*(t0 + t1);
0155     }
0156
0157     Real b1 = right_endpoint_derivative;
0158     if (boost::math::isnan(b1))
0159     {
0160         size_t n = length - 1;
0161         Real t0 = 4*(f[n-3] + third<Real>()*f[n - 1]);
0162         Real t1 = -(25*third<Real>()*f[n - 4] + f[n])/4  - 3*f[n - 2];
0163
0164         b1 = m_h_inv*(t0 + t1);
0165     }
0166
0167     // s(x) = \sum \alpha_i B_{3}( (x- x_i - a)/h )
0168     // Of course we must reindex from Kress's notation, since he uses negative indices which make C++ unhappy.
0169     m_beta.resize(length + 2, std::numeric_limits<Real>::quiet_NaN());
0170
0171     // Since the splines have compact support, they decay to zero very fast outside the endpoints.
0172     // This is often very annoying; we'd like to evaluate the interpolant a little bit outside the
0173     // boundary [a,b] without massive error.
0174     // A simple way to deal with this is just to subtract the DC component off the signal, so we need the average.
0175     // This algorithm for computing the average is recommended in
0176     // http://www.heikohoffmann.de/htmlthesis/node134.html
0177     Real t = 1;
0178     for (size_t i = 0; i < length; ++i)
0179     {
0180         if (boost::math::isnan(f[i]))
0181         {
0182             std::string err = "This function you are trying to interpolate is a nan at index " + std::to_string(i) + "\n";
0183             throw std::logic_error(err);
0184         }
0185         m_avg += (f[i] - m_avg) / t;
0186         t += 1;
0187     }
0188
0189
0190     // Now we must solve an almost-tridiagonal system, which requires O(N) operations.
0191     // There are, in fact 5 diagonals, but they only differ from zero on the first and last row,
0192     // so we can patch up the tridiagonal row reduction algorithm to deal with two special rows.
0193     // See Kress, equations 8.41
0194     // The the "tridiagonal" matrix is:
0195     // 1  0 -1
0196     // 1  4  1
0197     //    1  4  1
0198     //       1  4  1
0199     //          ....
0200     //          1  4  1
0201     //          1  0 -1
0202     // Numerical estimate indicate that as N->Infinity, cond(A) -> 6.9, so this matrix is good.
0203     std::vector<Real> rhs(length + 2, std::numeric_limits<Real>::quiet_NaN());
0204     std::vector<Real> super_diagonal(length + 2, std::numeric_limits<Real>::quiet_NaN());
0205
0206     rhs[0] = -2*step_size*a1;
0207     rhs[rhs.size() - 1] = -2*step_size*b1;
0208
0209     super_diagonal[0] = 0;
0210
0211     for(size_t i = 1; i < rhs.size() - 1; ++i)
0212     {
0213         rhs[i] = 6*(f[i - 1] - m_avg);
0214         super_diagonal[i] = 1;
0215     }
0216
0217
0218     // One step of row reduction on the first row to patch up the 5-diagonal problem:
0219     // 1 0 -1 | r0
0220     // 1 4 1  | r1
0221     // mapsto:
0222     // 1 0 -1 | r0
0223     // 0 4 2  | r1 - r0
0224     // mapsto
0225     // 1 0 -1 | r0
0226     // 0 1 1/2| (r1 - r0)/4
0227     super_diagonal[1] = 0.5;
0228     rhs[1] = (rhs[1] - rhs[0])/4;
0229
0230     // Now do a tridiagonal row reduction the standard way, until just before the last row:
0231     for (size_t i = 2; i < rhs.size() - 1; ++i)
0232     {
0233         Real diagonal = 4 - super_diagonal[i - 1];
0234         rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
0235         super_diagonal[i] /= diagonal;
0236     }
0237
0238     // Now the last row, which is in the form
0239     // 1 sd[n-3] 0      | rhs[n-3]
0240     // 0  1     sd[n-2] | rhs[n-2]
0241     // 1  0     -1      | rhs[n-1]
0242     Real final_subdiag = -super_diagonal[rhs.size() - 3];
0243     rhs[rhs.size() - 1] = (rhs[rhs.size() - 1] - rhs[rhs.size() - 3])/final_subdiag;
0244     Real final_diag = -1/final_subdiag;
0245     // Now we're here:
0246     // 1 sd[n-3] 0         | rhs[n-3]
0247     // 0  1     sd[n-2]    | rhs[n-2]
0248     // 0  1     final_diag | (rhs[n-1] - rhs[n-3])/diag
0249
0250     final_diag = final_diag - super_diagonal[rhs.size() - 2];
0251     rhs[rhs.size() - 1] = rhs[rhs.size() - 1] - rhs[rhs.size() - 2];
0252
0253
0254     // Back substitutions:
0255     m_beta[rhs.size() - 1] = rhs[rhs.size() - 1]/final_diag;
0256     for(size_t i = rhs.size() - 2; i > 0; --i)
0257     {
0258         m_beta[i] = rhs[i] - super_diagonal[i]*m_beta[i + 1];
0259     }
0260     m_beta[0] = m_beta[2] + rhs[0];
0261 }
0262
0263 template<class Real>
0264 Real cubic_b_spline_imp<Real>::operator()(Real x) const
0265 {
0266     // See Kress, 8.40: Since B3 has compact support, we don't have to sum over all terms,
0267     // just the (at most 5) whose support overlaps the argument.
0268     Real z = m_avg;
0269     Real t = m_h_inv*(x - m_a) + 1;
0270
0271     using std::max;
0272     using std::min;
0273     using std::ceil;
0274     using std::floor;
0275
0276     size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2))));
0277     size_t k_max = static_cast<size_t>((max)((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))), 0l));
0278     for (size_t k = k_min; k <= k_max; ++k)
0279     {
0280         z += m_beta[k]*b3_spline(t - k);
0281     }
0282
0283     return z;
0284 }
0285
0286 template<class Real>
0287 Real cubic_b_spline_imp<Real>::prime(Real x) const
0288 {
0289     Real z = 0;
0290     Real t = m_h_inv*(x - m_a) + 1;
0291
0292     using std::max;
0293     using std::min;
0294     using std::ceil;
0295     using std::floor;
0296
0297     size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2))));
0298     size_t k_max = static_cast<size_t>((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))));
0299
0300     for (size_t k = k_min; k <= k_max; ++k)
0301     {
0302         z += m_beta[k]*b3_spline_prime(t - k);
0303     }
0304     return z*m_h_inv;
0305 }
0306
0307 template<class Real>
0308 Real cubic_b_spline_imp<Real>::double_prime(Real x) const
0309 {
0310     Real z = 0;
0311     Real t = m_h_inv*(x - m_a) + 1;
0312
0313     using std::max;
0314     using std::min;
0315     using std::ceil;
0316     using std::floor;
0317
0318     size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2))));
0319     size_t k_max = static_cast<size_t>((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))));
0320
0321     for (size_t k = k_min; k <= k_max; ++k)
0322     {
0323         z += m_beta[k]*b3_spline_double_prime(t - k);
0324     }
0325     return z*m_h_inv*m_h_inv;
0326 }
0327
0328
0329 }}}
0330 #endif