Simple interpolation functions for C++.

Interpolation.h

/*
 * Interpolation functions!
 * Author: Ian Dunn
 */

#pragma once

template <typename T>
T LinearInterpolate(T p1, T p2, T x)
{
	return p1*(1 - x) + p2*x;
}

template <typename T>
T LinearInterpolate(T p[2], T x)
{
	return p[0]*(1 - x) + p[1]*x;
}

template <typename T>
T BilinearInterpolate(T p[2][2], T x, T y)
{
	T arr[2];
	arr[0] = LinearInterpolate(p[0], y);
	arr[1] = LinearInterpolate(p[1], y);
	return LinearInterpolate(arr, x);
}

template <typename T>
T TrilinearInterpolate(T p[2][2][2], T x, T y, T z)
{
	T arr[2];
	arr[0] = BilinearInterpolate(p[0], y, z);
	arr[1] = BilinearInterpolate(p[1], y, z);
	return LinearInterpolate(arr, x);
}

template <typename T>
T CubicInterpolate(T p[4], T x)
{
	return p[1] + (T) 0.5 * x*(p[2] - p[0] + x*(2*p[0] - 5*p[1] + 4*p[2] - p[3] + x*(3*(p[1] - p[2]) + p[3] - p[0])));
}

template <typename T>
T BicubicInterpolate(T p[4][4], T x, T y)
{
	T arr[4];
	arr[0] = CubicInterpolate(p[0], y);
	arr[1] = CubicInterpolate(p[1], y);
	arr[2] = CubicInterpolate(p[2], y);
	arr[3] = CubicInterpolate(p[3], y);
	return CubicInterpolate(arr, x);
}

template <typename T>
T TricubicInterpolate(T p[4][4][4], T x, T y, T z)
{
	T arr[4];
	arr[0] = BicubicInterpolate(p[0], y, z);
	arr[1] = BicubicInterpolate(p[1], y, z);
	arr[2] = BicubicInterpolate(p[2], y, z);
	arr[3] = BicubicInterpolate(p[3], y, z);
	return CubicInterpolate(arr, x);
}

Tests.cpp

#define CATCH_CONFIG_MAIN
#include <catch.hpp>
#include "Interpolation.h"
#include <array>


TEST_CASE("Linear interpolation (3) works", "[LinearInterpolation]") {
	REQUIRE(LinearInterpolate(0.0, 1.0, 0.5) == 0.5);
	REQUIRE(LinearInterpolate(-1.0, 1.0, 0.5) == 0.0);
	REQUIRE(LinearInterpolate(-1.0, 1.0, 0.25) == -0.5);
	REQUIRE(LinearInterpolate(-1.0, 1.0, 0.75) == 0.5);
}

TEST_CASE("Linear interpolation (2) works", "[LinearInterpolation]") {
	std::array<double, 2> p;

	p = {{0, 1}};
	REQUIRE(LinearInterpolate(p.data(), 0.5) == 0.5);
	p = {{-1, 1}};
	REQUIRE(LinearInterpolate(p.data(), 0.5) == 0.0);
	REQUIRE(LinearInterpolate(p.data(), 0.25) == -0.5);
	REQUIRE(LinearInterpolate(p.data(), 0.75) == 0.5);
}

TEST_CASE("Bilinear interpolation works", "[BilinearInterpolate]") {
	double p[2][2];

	p[0][0] = 1;
	p[0][1] = 2;
	p[1][0] = 3;
	p[1][1] = 4;
	REQUIRE(BilinearInterpolate(p, 0.0, 0.0) == 1.0);
	REQUIRE(BilinearInterpolate(p, 0.0, 1.0) == 2.0);
	REQUIRE(BilinearInterpolate(p, 1.0, 0.0) == 3.0);
	REQUIRE(BilinearInterpolate(p, 1.0, 1.0) == 4.0);

	REQUIRE(BilinearInterpolate(p, 0.0, 0.5) == 1.5);
	REQUIRE(BilinearInterpolate(p, 0.5, 0.0) == 2.0);
	REQUIRE(BilinearInterpolate(p, 1.0, 0.5) == 3.5);
	REQUIRE(BilinearInterpolate(p, 0.5, 1.0) == 3.0);

	REQUIRE(BilinearInterpolate(p, 0.5, 0.5) == 2.5);
}

TEST_CASE("Trilinear interpolation works", "[TrilinearInterpolate]") {
	double p[2][2][2];

	p[0][0][0] = 1;
	p[0][0][1] = 2;
	p[0][1][0] = 3;
	p[0][1][1] = 4;
	p[1][0][0] = 5;
	p[1][0][1] = 6;
	p[1][1][0] = 7;
	p[1][1][1] = 8;

	REQUIRE(TrilinearInterpolate(p, 0.0, 0.0, 0.0) == 1.0);
	REQUIRE(TrilinearInterpolate(p, 0.0, 0.0, 1.0) == 2.0);
	REQUIRE(TrilinearInterpolate(p, 0.0, 1.0, 0.0) == 3.0);
	REQUIRE(TrilinearInterpolate(p, 0.0, 1.0, 1.0) == 4.0);
	REQUIRE(TrilinearInterpolate(p, 1.0, 0.0, 0.0) == 5.0);
	REQUIRE(TrilinearInterpolate(p, 1.0, 0.0, 1.0) == 6.0);
	REQUIRE(TrilinearInterpolate(p, 1.0, 1.0, 0.0) == 7.0);
	REQUIRE(TrilinearInterpolate(p, 1.0, 1.0, 1.0) == 8.0);

	REQUIRE(TrilinearInterpolate(p, 0.0, 0.0, 0.5) == 1.5);
	REQUIRE(TrilinearInterpolate(p, 0.0, 0.5, 0.0) == 2.0);
	REQUIRE(TrilinearInterpolate(p, 0.0, 1.0, 0.5) == 3.5);
	REQUIRE(TrilinearInterpolate(p, 0.0, 0.5, 1.0) == 3.0);
	REQUIRE(TrilinearInterpolate(p, 0.5, 0.0, 0.0) == 3.0);
	REQUIRE(TrilinearInterpolate(p, 0.5, 0.0, 1.0) == 4.0);
	REQUIRE(TrilinearInterpolate(p, 0.5, 1.0, 0.0) == 5.0);
	REQUIRE(TrilinearInterpolate(p, 0.5, 1.0, 1.0) == 6.0);
	REQUIRE(TrilinearInterpolate(p, 1.0, 0.0, 0.5) == 5.5);
	REQUIRE(TrilinearInterpolate(p, 1.0, 0.5, 0.0) == 6.0);
	REQUIRE(TrilinearInterpolate(p, 1.0, 1.0, 0.5) == 7.5);
	REQUIRE(TrilinearInterpolate(p, 1.0, 0.5, 1.0) == 7.0);

	REQUIRE(TrilinearInterpolate(p, 0.0, 0.5, 0.5) == 2.5);
	REQUIRE(TrilinearInterpolate(p, 0.5, 0.0, 0.5) == 3.5);
	REQUIRE(TrilinearInterpolate(p, 0.5, 0.5, 0.0) == 4.0);

	REQUIRE(TrilinearInterpolate(p, 0.5, 0.5, 0.5) == 4.5);
}

TEST_CASE("Cubic interpolation works", "[CubicInterpolate]") {
	std::array<double, 4> p;

	p = {{0, 1, 2, 3}};
	REQUIRE(CubicInterpolate(p.data(), 0.0) == 1.0);
	REQUIRE(CubicInterpolate(p.data(), 0.5) == 1.5);
	REQUIRE(CubicInterpolate(p.data(), 1.0) == 2.0);
	REQUIRE(CubicInterpolate(p.data(), 0.25) == 1.25);
	REQUIRE(CubicInterpolate(p.data(), 0.75) == 1.75);

	p = {{0, 1, 3, 4}};
	REQUIRE(CubicInterpolate(p.data(), 0.0) == 1.0);
	REQUIRE(CubicInterpolate(p.data(), 0.5) == 2.0);
	REQUIRE(CubicInterpolate(p.data(), 1.0) == 3.0);
}