Interpolation

Interpolation is the process of constructing a function that takes on specified values at specified points. In engineering and science, one often has a number of data points, obtained by sampling or experimentation, which represent the values of a function for a limited number of values of the independent variable. It is often required to interpolate (i.e. estimate) the value of that function for an intermediate value of the independent variable.

Smile package smile.interpolation provides a variety of algorithms on 1d and 2d data. These algorithms implement the interface Interpolation (or Interpolation2D for 2d data), of which the method interpolate takes a value and return an interploated value.

Piecewise Linear Interpolation

LinearInterpolation is quick and easy, but it is not very precise. Another disadvantage is that the interpolant is not differentiable at the control points.


    val x = Array(0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
    val y = Array(0.0, 0.8415, 0.9093, 0.1411, -0.7568, -0.9589, -0.2794)
    val points = (0 until x.length).map(i => Array(x(i), y(i))).toArray

    val linear = new LinearInterpolation(x, y)

    val data = (0 to 60).map { i =>
      val x = i * 0.1
      val y = linear.interpolate(x)
      Array(x, y)
    }.toArray

    val canvas = plot(points, '@', Color.BLACK)
    canvas.add(LinePlot.of(data, Color.RED))
    show(canvas)
    

    double[] x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
    double[] y = {0.0, 0.8415, 0.9093, 0.1411, -0.7568, -0.9589, -0.2794};
    double[][] points = new double[x.length][2];
    for (int i = 0; i < x.length; i++) {
      points[i][0] = x[i];
      points[i][1] = y[i];
    }

    var linear = new LinearInterpolation(x, y);
    double[][] data = new double[61][2];
    for (int i = 0; i < data.length; i++) {
      data[i][0] = i * 0.1;
      data[i][1] = linear.interpolate(data[i][0]);
    }

    var canvas = ScatterPlot.of(points, '@', Color.BLACK).canvas();
    canvas.add(LinePlot.of(data, Color.RED));
    canvas.window();
          

BilinearInterpolation is an extension of linear interpolation for interpolating functions of two variables on a regular grid. The key idea is to perform linear interpolation first in one direction, and then again in the other direction.

Cubic Spline Interpolation

Spline interpolation uses low-degree polynomials in each of the intervals, and chooses the polynomial pieces such that they fit smoothly together. The resulting function is called a spline. The natural cubic spline is piecewise cubic and twice continuously differentiable. Furthermore, its second derivative is zero at the end points.

Like polynomial interpolation, spline interpolation incurs a smaller error than linear interpolation and the interpolant is smoother. However, the interpolant is easier to evaluate than the high-degree polynomials used in polynomial interpolation. It also does not suffer from Runge's phenomenon.


    val cubic = new CubicSplineInterpolation1D(x, y)

    val data = (0 to 60).map { i =>
      val x = i * 0.1
      val y = cubic.interpolate(x)
      Array(x, y)
    }.toArray

    val canvas = plot(points, '@', Color.BLACK)
    canvas.add(LinePlot.of(data, Color.BLUE))
    show(canvas)
          

    var cubic = new CubicSplineInterpolation1D(x, y);

    double[][] data = new double[61][2];
    for (int i = 0; i < data.length; i++) {
      data[i][0] = i * 0.1;
      data[i][1] = cubic.interpolate(data[i][0]);
    }

    var canvas = ScatterPlot.of(points, '@', Color.BLACK).canvas();
    canvas.add(LinePlot.of(data, Color.BLUE));
    canvas.window();
          

For 2d grid, we provide two implementations: CubicSplineInterpolation2D and BicubicInterpolation. CubicSplineInterpolation2D is similar to one-dimensional splines as it guarantees the continuity of the first and second function derivatives. In contrast, BicubicInterpolation guarantees continuity of only gradient and cross-derivative. Second derivatives could be discontinuous.

In image processing, bicubic interpolation is often chosen over bilinear interpolation or nearest neighbor in image resampling, when speed is not an issue. Images resampled with bicubic interpolation are smoother and have fewer interpolation artifacts.


    val x1 = Array(1950.0, 1960, 1970, 1980, 1990)
    val x2 = Array(10.0, 20, 30)
    val y = Array(
            Array(150.697, 199.592, 187.625),
            Array(179.323, 195.072, 250.287),
            Array(203.212, 179.092, 322.767),
            Array(226.505, 153.706, 426.730),
            Array(249.633, 120.281, 598.243)
        )

    val canvas = heatmap(y, Palette.jet(256))
    canvas.setTitle("Original")
    show(canvas)

    val cubic = new CubicSplineInterpolation2D(x1, x2, y)

    val data = Array.ofDim[Double](101, 101)
    for (i <- 0 to 100; j <- 0 to 100)
      data(i)(j) = cubic.interpolate(1950 + i*0.4, 10 + j*0.2)

    val cubicPlot = heatmap(data, Palette.jet(256))
    cubicPlot.setTitle("Cubic")
    show(cubicPlot)

    val bicubic = new BicubicInterpolation(x1, x2, y)

    for (i <- 0 to 100; j <- 0 to 100)
      data(i)(j) = bicubic.interpolate(1950 + i*0.4, 10 + j*0.2)

    val bicubicPlot = heatmap(data, Palette.jet(256))
    bicubicPlot.setTitle("Bicubic")
    show(bicubicPlot)
          

    double[] x1 = {1950.0, 1960, 1970, 1980, 1990};
    double[] x2 = {10.0, 20, 30};
    double[][] y = {
            {150.697, 199.592, 187.625},
            {179.323, 195.072, 250.287},
            {203.212, 179.092, 322.767},
            {226.505, 153.706, 426.730},
            {249.633, 120.281, 598.243}
        };

    var canvas = Heatmap.of(y, Palette.jet(256)).canvas();
    canvas.setTitle("Original");
    canvas.window();

    var cubic = new CubicSplineInterpolation2D(x1, x2, y);

    double[][] data = new double[101][101];
    for (int i = 0; i <= 100; i++) {
        for (int j = 0; j <= 100; j++) {
            data[i][j] = cubic.interpolate(1950 + i*0.4, 10 + j*0.2);
        }
    }

    var cubicPlot = Heatmap.of(data, Palette.jet(256)).canvas();
    cubicPlot.setTitle("Cubic");
    cubicPlot.window();

    var bicubic = new BicubicInterpolation(x1, x2, y);

    for (int i = 0; i <= 100; i++) {
        for (int j = 0; j <= 100; j++) {
            data[i][j] = cubic.interpolate(1950 + i*0.4, 10 + j*0.2);
        }
    }

    var bicubicPlot = Heatmap.of(data, Palette.jet(256)).canvas();
    bicubicPlot.setTitle("Bicubic");
    bicubicPlot.window();
          

Kriging Interpolation

Kriging interpolation can be used for the data points irregularly distributed in space. Kriging belongs to the family of linear least squares estimation algorithms, also known as Gauss-Markov estimation or Gaussian process regression. We implement ordinary kriging for interpolation with power variogram.


    val kriging = new KrigingInterpolation1D(x, y)

    val data = Array.ofDim[Double](61, 2)
    for (i <- 0 to 60) {
      data(i)(0) = i * 0.1
      data(i)(1) = kriging.interpolate(data(i)(0))
    }

    val canvas = plot(points, '@', Color.BLACK)
    canvas.add(LinePlot.of(data, Color.CYAN))
    show(canvas)
          

    var kriging = new KrigingInterpolation1D(x, y);

    double[][] data = new double[61][2];
    for (int i = 0; i < data.length; i++) {
      data[i][0] = i * 0.1;
      data[i][1] = kriging.interpolate(data[i][0]);
    }

    var canvas = ScatterPlot.of(points, '@', Color.BLACK).canvas();
    canvas.add(LinePlot.of(data, Color.CYAN));
    canvas.window();
          

RBF Interpolation

Radial basis function interpolation is a popular method for the data points are irregularly distributed in space. In its basic form, radial basis function interpolation is in the form y(x) = Σ wi φ(||x-ci||), where the approximating function y(x) is represented as a sum of N radial basis functions φ, each associated with a different center ci, and weighted by an appropriate coefficient wi. For distance, one usually chooses Euclidean distance. The weights wi can be estimated using the matrix methods of linear least squares, because the approximating function is linear in the weights.

The points ci often called the centers or collocation points of the RBF interpolant. Note also that the centers ci can be located at arbitrary points in the domain, and do not require a grid. For certain RBF exponential convergence has been shown. Radial basis functions were successfully applied to problems as diverse as computer graphics, neural networks, for the solution of differential equations via collocation methods and many other problems.

Other popular choices for φ comprise the Gaussian function and the so-called thin plate splines. Thin plate splines result from the solution of a variational problem. The advantage of the thin plate splines is that their conditioning is invariant under scalings. Gaussians, multi-quadrics and inverse multi-quadrics are infinitely smooth and involve a scale or shape parameter, r0 > 0. Decreasing r0 tends to flatten the basis function. For a given function, the quality of approximation may strongly depend on this parameter. In particular, increasing r0 has the effect of better conditioning (the separation distance of the scaled points increases).

A variant on RBF interpolation is normalized radial basis function (NRBF) interpolation, in which we require the sum of the basis functions to be unity. NRBF arises more naturally from a Bayesian statistical perspective. However, there is no evidence that either the NRBF method is consistently superior to the RBF method, or vice versa.


    val rbf = new RBFInterpolation1D(x, y, new smile.math.rbf.GaussianRadialBasis())

    val data = Array.ofDim[Double](61, 2)
    for (i <- 0 to 60) {
      data(i)(0) = i * 0.1
      data(i)(1) = rbf.interpolate(data(i)(0))
    }

    val canvas = plot(points, '@', Color.BLACK)
    canvas.add(LinePlot.of(data, Color.GREEN))
    show(canvas)
          

    var rbf = new RBFInterpolation1D(x, y, new smile.math.rbf.GaussianRadialBasis());

    double[][] data = new double[61][2];
    for (int i = 0; i < data.length; i++) {
      data[i][0] = i * 0.1;
      data[i][1] = rbf.interpolate(data[i][0]);
    }

    var canvas = ScatterPlot.of(points, '@', Color.BLACK).canvas();
    canvas.add(LinePlot.of(data, Color.GREEN));
    canvas.window();
          

Shepard Interpolation

Shepard's interpolation is a special case of normalized radial basis function interpolation if the function φ(r) goes to infinity as r → 0, and is finite for r > 0. In this case, the weights wi are just equal to the respective function values yi. So we need not solve linear equations, and thus it works for very large N. An example of such φ is φ(r) = r-p with (typically) 1 < p ≤ 3.


    val shepard = new ShepardInterpolation1D(x, y, 3)

    val data = Array.ofDim[Double](61, 2)
    for (i <- 0 to 60) {
      data(i)(0) = i * 0.1
      data(i)(1) = shepard.interpolate(data(i)(0))
    }

    val canvas = plot(points, '@', Color.BLACK)
    canvas.add(LinePlot.of(data, Color.PINK))
    show(canvas)
          

    var shepard = new ShepardInterpolation1D(x, y, 3);

    double[][] data = new double[61][2];
    for (int i = 0; i < data.length; i++) {
      data[i][0] = i * 0.1;
      data[i][1] = shepard.interpolate(data[i][0]);
    }

    var canvas = ScatterPlot.of(points, '@', Color.BLACK).canvas();
    canvas.add(LinePlot.of(data, Color.PINK));
    canvas.window();
          

Shepard's interpolation is rarely as accurate as the well-tuned application of other radial basis functions. However, it is simple, fast, and often sufficient for quick and dirty applications.

Laplace Interpolation

Laplace's interpolation can restore missing or unmeasured values on a 2-dimensional evenly spaced regular grid. In some sense, Laplace interpolation produces the smoothest possible interpolant, which are obtained by solving a very sparse linear equations with biconjugate gradient method.


    val zz = Array.ofDim[Double](101, 101)
    val ww = Array.ofDim[Double](101, 101)
    for (i <- 0 to 100; j <- 0 to 100) {
      zz(i)(j) = if (java.lang.Math.random() < 0.2) Double.NaN else data(i)(j)
      ww(i)(j) = zz(i)(j)
    }

    val canvas = heatmap(ww, Palette.jet(256))
    canvas.setTitle("Missing Values")
    show(canvas)

    LaplaceInterpolation.interpolate(zz)
    val canvas = heatmap(zz, Palette.jet(256))
    canvas.setTitle("Laplace")
    show(canvas)
          

    double[][] zz = new double[101][101];
    double[][] ww = new double[101][101];
    for (int i = 0; i <= 100; i++) {
        for (int j = 0; j <= 100; j++) {
            zz[i][j] = Math.random() < 0.2 ? Double.NaN : data[i][j];
            ww[i][j] = zz[i][j];
        }
    }

    var canvas = Heatmap.of(ww, Palette.jet(256)).canvas();
    canvas.setTitle("Missing Values");
    canvas.window();

    LaplaceInterpolation.interpolate(zz);
    var canvas = Heatmap.of(zz, Palette.jet(256)).canvas();
    canvas.setTitle("Laplace");
    canvas.window();
          
←   Wavelet Graph  →
Fork me on GitHub