...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#include <boost/math/interpolators/pchip.hpp> namespace boost::math::interpolators { template <class RandomAccessContainer> class pchip { public: using Real = RandomAccessContainer::value_type; pchip(RandomAccessContainer&& abscissas, RandomAccessContainer&& ordinates, Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(), Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN()); Real operator()(Real x) const; Real prime(Real x) const; void push_back(Real x, Real y); friend std::ostream& operator<<(std::ostream & os, const pchip & m); }; } // namespaces
The PCHIP interpolant takes non-equispaced data and interpolates between them via cubic Hermite polynomials whose slopes are chosen so that the resulting interpolant is monotonic; see Fritsch and Carlson for details. The interpolant is C1 and evaluation has 𝑶(log(N)) complexity. An example usage is as follows:
std::vector<double> x{1, 5, 9 , 12}; std::vector<double> y{8,17, 4, -3}; using boost::math::interpolators::pchip; auto spline = pchip(std::move(x), std::move(y)); // evaluate at a point: double z = spline(3.4); // evaluate derivative at a point: double zprime = spline.prime(3.4);
Periodically, it is helpful to see what data the interpolator has, and the slopes it has chosen. This can be achieved via
std::cout << spline << "\n";
Note that the interpolator is pimpl'd, so that copying the class is cheap,
and hence it can be shared between threads. (The call operator and .prime()
are threadsafe; push_back
is
not.)
This interpolant can be updated in constant time. Hence we can use boost::circular_buffer
to do real-time interpolation:
#include <boost/circular_buffer.hpp> ... boost::circular_buffer<double> initial_x{1,2,3,4}; boost::circular_buffer<double> initial_y{4,5,6,7}; auto circular_pchip = pchip(std::move(initial_x), std::move(initial_y)); // interpolate via call operation: double y = circular_pchip(3.5); // add new data: circular_pchip.push_back(5, 8); // interpolate at 4.5: y = circular_pchip(4.5);
This interpolator chooses the slopes and forwards data to the cubic Hermite
interpolator, so the performance is stated in the documentation for cubic_hermite.hpp
.