/* interpolation/bilinear.c
*
* Copyright 2012 David Zaslavsky
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_interp2d.h>
#define IDX2D(i, j, xsize, ysize) ((j) * (xsize) + (i))
static int
bilinear_init(void * state, const double xa[], const double ya[],
const double za[], size_t xsize, size_t ysize)
{
return GSL_SUCCESS;
}
static int
bilinear_eval(const void * state, const double xarr[], const double yarr[],
const double zarr[], size_t xsize, size_t ysize,
double x, double y, gsl_interp_accel * xa,
gsl_interp_accel * ya, double * z)
{
double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;
double dx, dy;
double t, u;
size_t xi, yi;
if (xa != NULL)
xi = gsl_interp_accel_find(xa, xarr, xsize, x);
else
xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);
if (ya != NULL)
yi = gsl_interp_accel_find(ya, yarr, ysize, y);
else
yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);
xmin = xarr[xi];
xmax = xarr[xi + 1];
ymin = yarr[yi];
ymax = yarr[yi + 1];
zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];
zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];
zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];
zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];
dx = xmax - xmin;
dy = ymax - ymin;
t = (x - xmin)/dx;
u = (y - ymin)/dy;
*z = (1.-t)*(1.-u)*zminmin + t*(1.-u)*zmaxmin + (1.-t)*u*zminmax + t*u*zmaxmax;
return GSL_SUCCESS;
}
static int
bilinear_deriv_x(const void * state, const double xarr[],
const double yarr[], const double zarr[],
size_t xsize, size_t ysize, double x, double y,
gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_p)
{
double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;
double dx, dy;
double dt, u;
size_t xi, yi;
if (xa != NULL)
xi = gsl_interp_accel_find(xa, xarr, xsize, x);
else
xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);
if (ya != NULL)
yi = gsl_interp_accel_find(ya, yarr, ysize, y);
else
yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);
xmin = xarr[xi];
xmax = xarr[xi + 1];
ymin = yarr[yi];
ymax = yarr[yi + 1];
zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];
zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];
zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];
zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];
dx = xmax - xmin;
dy = ymax - ymin;
dt = 1./dx; /* partial t / partial x */
u = (y - ymin)/dy;
*z_p = dt*(-(1.-u)*zminmin + (1.-u)*zmaxmin - u*zminmax + u*zmaxmax);
return GSL_SUCCESS;
}
static int
bilinear_deriv_y(const void * state, const double xarr[],
const double yarr[], const double zarr[],
size_t xsize, size_t ysize, double x, double y,
gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_p)
{
double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;
double dx, dy;
double t, du;
size_t xi, yi;
if (xa != NULL)
xi = gsl_interp_accel_find(xa, xarr, xsize, x);
else
xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);
if (ya != NULL)
yi = gsl_interp_accel_find(ya, yarr, ysize, y);
else
yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);
xmin = xarr[xi];
xmax = xarr[xi + 1];
ymin = yarr[yi];
ymax = yarr[yi + 1];
zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];
zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];
zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];
zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];
dx = xmax - xmin;
dy = ymax - ymin;
t = (x - xmin)/dx;
du = 1./dy; /* partial u / partial y */
*z_p = du*(-(1.-t)*zminmin - t*zmaxmin + (1.-t)*zminmax + t*zmaxmax);
return GSL_SUCCESS;
}
static int
bilinear_deriv2(const void * state, const double xarr[],
const double yarr[], const double zarr[],
size_t xsize, size_t ysize, double x, double y,
gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_pp)
{
*z_pp = 0.0;
return GSL_SUCCESS;
}
static int
bilinear_derivxy(const void * state, const double xarr[],
const double yarr[], const double zarr[],
size_t xsize, size_t ysize, double x, double y,
gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_pp)
{
double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;
double dx, dy;
double dt, du;
size_t xi, yi;
if (xa != NULL)
xi = gsl_interp_accel_find(xa, xarr, xsize, x);
else
xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);
if (ya != NULL)
yi = gsl_interp_accel_find(ya, yarr, ysize, y);
else
yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);
xmin = xarr[xi];
xmax = xarr[xi + 1];
ymin = yarr[yi];
ymax = yarr[yi + 1];
zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];
zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];
zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];
zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];
dx = xmax - xmin;
dy = ymax - ymin;
dt = 1./dx; /* partial t / partial x */
du = 1./dy; /* partial u / partial y */
*z_pp = dt*du*(zminmin-zmaxmin-zminmax+zmaxmax);
return GSL_SUCCESS;
}
static const gsl_interp2d_type bilinear_type = {
"bilinear",
2,
NULL,
&bilinear_init,
&bilinear_eval,
&bilinear_deriv_x,
&bilinear_deriv_y,
&bilinear_deriv2,
&bilinear_derivxy,
&bilinear_deriv2,
NULL
};
const gsl_interp2d_type * gsl_interp2d_bilinear = &bilinear_type;
#undef IDX2D