Introduction

Interpolation is a technique to approximate continuous data from sample points. If you already heard about vector graphics, the core idea there is to use interpolation to describe shapes. Another very common application of interpolation is for animations, where you change one property of an object gradually by interpolating between two values.

Interpolation is also useful for scientific purposes, for instance, numerical integration methods actually integrate an interpolation of a function when one doesn't know a closed form solution.

Here I will present an overview of what I have found for JavaScript.

Overview

All the interpolation methods considered are piecewise polynomial interpolations, nearest, previous and next with degree 0, linear with degree 1, and splines with degree 3.

![interp1_methods.svg](raw.githubusercontent.com/o-alexandre-felip..

The previous method consists of the nearest point before the given point

function to_previous(ts: number[], ys: number[], t: number): number{
  let y;
  for(let i = 0; i < ts.length; ++i){
    if(ts[i] <= t) y = ys[i];
    else return y;
  }
}

This method is supported by javascript interp1, and can more efficiently evaluated using binary-search-bounds.

const bsb = require('binary-search-bounds');
function to_previous(ts: number[], ys: number[], t: number): number {
  return ys[bsb.le(ts, t)];
}

function to_next(ts: number[], ys: number[], t: number): number {
  for(let i = 0; i < ts.length; ++i){
    if(ts[i] >= t) return ys[i];
  }
}

const bsb = require('binary-search-bounds');
function to_next(ts: number[], ys: number[], t: number): number {
  return ys[bsb.ge(ts, t)];
}

Nearest

function to_nearest(ts: number[], ys: number[], t: number): number {
  for(let i = 0; i < ts.length; ++i){
    if(ts[i] <= t && ts[i+1] >= t) {
      if(t - ts[i] < ts[i+1] - t) return ys[i];
      else return ys[i+1];
    }
  }
}

Linear

function piecewise_linear(ts: number[], ys: number[], t: number): number {
  for(let i = 0; i < ts.length; ++i){
    if(ts[i] <= t && ts[i+1] >= t) {
      // Use the equation for the line passing through 
      // (ts[i], ys[i]) and (ts[i+1], ys[i+1])
      return ys[i] + (t - ts[i]) * (ys[i+1] - ys[i]) / (ts[i+1] - ts[i])
    }
  }
}

Using polylinear-scale module

const PolyLinearScale = require('polylinear-scale');
function piecewise_linear(ts: number[], ys: number[], t: number): number {
  const curve = new PolyLinearScale(ts, ys, true);
  return curve(t);
}

Cubic Spline

This method of interpolation is more complicated and currently is not supported by interp1. You can get do it using cubic-spline module, but it uses a far from optimal spline construction algorithm.

function interp_spline(ts: number[], ys: number[], t: number): number {
  const curve = new CubicSpline(ts, ys);
  return curve.at(t);
}

Implementation analysis

I had a read in the source code of the packages listed above.

interp1@1.0.18

The nicer code among those(?) is well commented, has inline type definition (written in typescript). It incurs in some overhead by accepting unsorted data, the first thing it does is to make an array of pairs p[i] = [x[i], y[i]], and then sort by p[0]. This copy operation may require many memory allocations that could be slow and also sorting operation that is \( O(n \log n) \). Then for each element in qs, it finds an index using a binary search \( O(\log n) \) complexity, and compute the interpolation \( O(1) \).

polylinear-scale@1.1.6

This module uses a function that creates an interpolator, the initialization is very simple as it only copies references to the input arrays and check if they have the same lengh, and it create a scoped function, so the construction takes \( O(1) )) time. However the evaluation of a single point uses a linear search, taking \( O(n) )) time.

binary-search@1.3.6

Binary search is single function package for binary search on a sorted array in \( O(n) \), this is not properly an interpolation function but since binary search is used to find the interpolation interval the performance of this is very pertinent to the comparison. This package as is cannot be used for interpolation because it will return the index of a value only if it compares exactly to the needle.

binary-search-bounds@2.0.5

This code is an obscure and interesting code, it implement five variations of binary search, each in less than 10 short lines of code (< 60 cols, with indentation). Not only that, it makes comparison optional, and accept limits to the search range. This is by far the most efficient code, e.g. 2x faster than binary-search, while supporting the same functionality as a particular case.

cubic-spline@3.0.3

This package does a non-trivial interpolation, spline interpolation makes sure the first derivative is continuous and uses a remaining degree of freedom to make the second derivative continuous.

This package assumes sorted input, so it makes no unnecessary copy. It uses a binary search to get the interpolation interval, so evaluation takes \( O(n) \) time, the downside of it is that it solves the equations with a generic purpose method, that requires \( O(n^2) \) memory and time.

Comparison and more to come

Here is the benchmark of interpolation through 10 points evaluated at 10 random points. polylinear-scale and cubic-spline require some initialization that is considered separately. This initialization will dominate only a few points will be evaluated, and will be amortized if many points are computed.

There is room for improvement (sticking with pure javascript). I worked on the cubic-spline algorithm. My implementation of spline is faster for evaluation than the available methods available in NPM. The above plot replacing cubic-spline with my implementation looks like this

And the method scales well to large arrays. Excluding polylinear-scale, whose evaluation complexity is (( O(n) )) from the comparison, interpolating 1000 points doing 1000 evaluations it is still better than interp1, the only viable. These results encourage me to prepare a new package.

Soon I will write a post about the implementation details, numerical stability analysis, and covering the math behind this. I also plan to write about how to link the Bezier notation (using control points) and the curve from the polynomial (mathematical expression).

Alexandre Felipe's blog