busstoptaktik / geodesy

Rust geodesy
Apache License 2.0
66 stars 6 forks source link

Alternative interpolation algorithms #51

Open busstoptaktik opened 1 year ago

busstoptaktik commented 1 year ago

For grid interpolation, this article and its code companion shows Constrained Bicubic, which is "bicubic without overshoots"

// Constrained bicubic interpolation
template<typename T>
void cbi (const Matrix<T>& in, Matrix<T>& out)
{
  //scaling factors
  double xr = (double)out.cols () / (in.cols () - 1);
  double yr = (double)out.rows () / (in.rows () - 1);

  assert (xr >= 1 && yr >= 1);

  //weighting function
  auto w = [](double x, double y) -> double {
    return x * x * y * y * (9 - 6 * x - 6 * y + 4 * x * y);
  };

  for (int i = 0; i < out.rows (); i++)
  {
    int ii = (int)floor (i / xr);
    for (int j = 0; j < out.cols (); j++)
    {
      int jj = (int)floor (j / yr);
      T v00 = in (ii, jj), v01 = in (ii, jj + 1),
        v10 = in (ii + 1, jj), v11 = in (ii + 1, jj + 1);
      double fi = i / xr - ii, fj = j / yr - jj;
      out (i, j) = w(1 - fi, 1 - fj) * v00 + w(1 - fi, fj) * v01 +
        w(fi, 1 - fj) * v10 + w(fi, fj) * v11;
    }
  }
}