Bulldozer00's Blog

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?