Home > C++ > Linear Interpolation In C++

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:

linear interpolation

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?

Categories: C++ Tags: ,
  1. Vladi
    May 10, 2016 at 2:56 am

    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 (:

    • May 10, 2016 at 4:23 am

      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.

  2. Arup Debnath
    July 22, 2019 at 4:15 am

    When the plot size upto 200 Sqm. Than permissible ground coverage equal to 65% of plot size
    When the plot size above 500 Sqm. Than permissible ground coverage equal to 50% of plot size.

    Now how to calculate the permissible ground coverage of Plot size 397.96 ?

    • July 22, 2019 at 7:10 am

      Given:

      x = 397.96
      x2 = 200, y2 = 65
      x3 = 500, y3 = 50

      We have:

      y = 65 + (397.96 – 200) * (50 – 65) / (500 – 200) = 55.10

  1. No trackbacks yet.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.

%d bloggers like this: