### Archive

Posts Tagged ‘Linear Interpolation’

## 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