From 33e2b0e4fa292fc6978b1d49b8d9d3a66dc3e9da Mon Sep 17 00:00:00 2001 From: tlatorre Date: Fri, 19 Oct 2018 13:00:26 -0500 Subject: add interp2d() for fast bilinear 2D interpolation --- src/misc.c | 47 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) (limited to 'src/misc.c') diff --git a/src/misc.c b/src/misc.c index d20bc26..605f3ea 100644 --- a/src/misc.c +++ b/src/misc.c @@ -3,6 +3,7 @@ #include /* for size_t */ #include #include "vector.h" +#include static struct { int n; @@ -338,6 +339,52 @@ double kahan_sum(double *x, size_t n) return sum; } +double interp2d(double x, double y, double *xp, double *yp, double *zp, size_t n1, size_t n2) +{ + /* A fast bilinear interpolation routine which assumes that the values in + * `xp` and `yp` are evenly spaced. + * + * `zp` should be a 2D array indexed as zp[i,j] = zp[i*n2 + j]. + * + * If x < xp[0], x > xp[n-1], y < yp[0], or y > yp[n-1] prints an error and + * exits. + * + * See https://en.wikipedia.org/wiki/Bilinear_interpolation. */ + size_t i, j; + double q11, q12, q21, q22; + + if (x <= xp[0]) { + fprintf(stderr, "x < xp[0]!\n"); + exit(1); + } + + if (y <= yp[0]) { + fprintf(stderr, "y < xp[0]!\n"); + exit(1); + } + + i = (x-xp[0])/(xp[1]-xp[0]); + + if (i > n1-2) { + fprintf(stderr, "i > n1 - 2!\n"); + exit(1); + } + + j = (y-yp[0])/(yp[1]-yp[0]); + + if (j > n2-2) { + fprintf(stderr, "j > n2 - 2!\n"); + exit(1); + } + + q11 = zp[i*n2+j]; + q12 = zp[i*n2+j+1]; + q21 = zp[(i+1)*n2+j]; + q22 = zp[(i+1)*n2+j+1]; + + return (q11*(xp[i+1]-x)*(yp[j+1]-y) + q21*(x-xp[i])*(yp[j+1]-y) + q12*(xp[i+1]-x)*(y-yp[j]) + q22*(x-xp[i])*(y-yp[j]))/((xp[i+1]-xp[i])*(yp[i+1]-yp[i])); +} + double interp1d(double x, double *xp, double *yp, size_t n) { /* A fast interpolation routine which assumes that the values in `xp` are -- cgit