## Linear Interpolation In C++

In scientific programming and embedded sensor systems applications, linear interpolation is often used to estimate a value from a series of discrete data points. The problem is stated and the solution is given as follows:

The solution *assumes* that any two points in a set of given data points represents a straight line. Hence, it takes the form of the classic textbook equation for a line, y = b + mx, where b is the y intercept and m is the slope of the line.

If the set of data points does not represent a linear underlying phenomenon, more sophisticated polynomial interpolation techniques that use additional data points around the point of interest can be utilized to get a more accurate estimate.

The code snippets below give the definition and implementation of a C++ functor class that performs linear interpolation. I chose to use a vector of pairs to represent the set of data points. What would’ve you used?

**Interpolator.h**

#ifndef INTERPOLATOR_H_ #define INTERPOLATOR_H_ #include <utility> #include <vector> class Interpolator { public: //On construction, we take in a vector of data point pairs //that represent the line we will use to interpolate Interpolator(const std::vector<std::pair<double, double>>& points); //Computes the corresponding Y value //for X using linear interpolation double findValue(double x) const; private: //Our container of (x,y) data points //std::pair::<double, double>.first = x value //std::pair::<double, double>.second = y value std::vector<std::pair<double, double>> _points; }; #endif /* INTERPOLATOR_H_ */

**Interpolator.cpp**

#include "Interpolator.h" #include <algorithm> #include <stdexcept> Interpolator::Interpolator(const std::vector<std::pair<double, double>>& points) : _points(points) { //Defensive programming. Assume the caller has not sorted the table in //in ascending order std::sort(_points.begin(), _points.end()); //Ensure that no 2 adjacent x values are equal, //lest we try to divide by zero when we interpolate. const double EPSILON{1.0E-8}; for(std::size_t i=1; i<_points.size(); ++i) { double deltaX{std::abs(_points[i].first - _points[i-1].first)}; if(deltaX < EPSILON ) { std::string err{"Potential Divide By Zero: Points " + std::to_string(i-1) + " And " + std::to_string(i) + " Are Too Close In Value"}; throw std::range_error(err); } } } double Interpolator::findValue(double x) const { //Define a lambda that returns true if the x value //of a point pair is < the caller's x value auto lessThan = [](const std::pair<double, double>& point, double x) {return point.first < x;}; //Find the first table entry whose value is >= caller's x value auto iter = std::lower_bound(_points.cbegin(), _points.cend(), x, lessThan); //If the caller's X value is greater than the largest //X value in the table, we can't interpolate. if(iter == _points.cend()) { return (_points.cend() - 1)->second; } //If the caller's X value is less than the smallest X value in the table, //we can't interpolate. if(iter == _points.cbegin() and x <= _points.cbegin()->first) { return _points.cbegin()->second; } //We can interpolate! double upperX{iter->first}; double upperY{iter->second}; double lowerX{(iter - 1)->first}; double lowerY{(iter - 1)->second}; double deltaY{upperY - lowerY}; double deltaX{upperX - lowerX}; return lowerY + ((x - lowerX)/ deltaX) * deltaY; }

In the constructor, the code attempts to establish the invariant conditions required before any post-construction interpolation can be attempted:

- It sorts the data points in ascending X order – just in case the caller “forgot” to do it.
- It ensures that no two adjacent X values have the same value – which could cause a divide-by-zero during the interpolation computation. If this constraint is not satisfied, an exception is thrown.

Here are the unit tests I ran on the implementation:

#include "catch.hpp" #include "Interpolator.h" TEST_CASE( "Test", "[Interpolator]" ) { //Construct with an unsorted set of data points Interpolator interp1{ { //{X(i),Y(i) {7.5, 32.0}, {1.5, 20.0}, {0.5, 10.0}, {3.5, 28.0}, } }; //X value too low REQUIRE(10.0 == interp1.findValue(.2)); //X value too high REQUIRE(32.0 == interp1.findValue(8.5)); //X value in 1st sorted slot REQUIRE(15.0 == interp1.findValue(1.0)); //X value in last sorted slot REQUIRE(30.0 == interp1.findValue(5.5)); //X value in second sorted slot REQUIRE(24.0 == interp1.findValue(2.5)); //Two points have the same X value std::vector<std::pair<double, double>> badData{ {.5, 32.0}, {1.5, 20.0}, {3.5, 28.0}, {1.5, 10.0}, }; REQUIRE_THROWS(Interpolator interp2{badData}); }

How can this code be improved?

Total nitpicks most of these, but here goes:

– Precompute the deltaY/deltaX for each segment. This seems to be the most expensive part of the calculation.

– rename the findValue function to be operator() so this object can be used as a functor in different algorithms

– if out of bounds values are likely, you can store the smallest and largest x in the object and compare to them to save the indirection of access to the vector data, no need to run the whole search.

– the lessThen lambda is created on the stack in every function call. The compiler will probably optimize this, but there’s no reason the can’t be a static inline class member function (or just a free function)

– template on the data type to be able to use floats for your values

– You’re assuming linear behaviour between data point, not along the whole thing, right? Otherwise just calculate the slope and be done with it (:

They’re not nits. They’re terrific suggestions. Thanks.

On your last sentence, the answer is yes. Although the intro picture implies otherwise, I’m only assuming linearity between any two points.