/*--------------------------------------------------------------------
* The GMT-system: @(#)gmt_bcr.h 3.13 06/17/99
*
* Copyright (c) 1991-1999 by P. Wessel and W. H. F. Smith
* See COPYING file for copying and redistribution conditions.
*
* 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; version 2 of the License.
*
* 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.
*
* Contact info: www.soest.hawaii.edu/gmt
*--------------------------------------------------------------------*/
/* bcr.h -- stuff for implementing bicubic rectangle grid interpolation.
Has bilinear case as a subset.
The bilinear interpolation will give continuity of the function z only.
The bicubic interpolation will give continuity of z, dz/dx, dz/dy, and
d2z/dxdy.
Designed to operate with at least one spare row/column around each edge,
to allow derivative computations.
Follows the outline in Lancaster and Salkauskas, section 9.3
Meant to replace Taylor routines in GMT in version 3.0 of GMT.
Author: Walter H F Smith
Date: 23 September, 1993
Revised: 20 September, 1998
This include file defines structures and functions used.
*/
EXTERN_MSC struct BCR bcr;
EXTERN_MSC void GMT_bcr_init (struct GRD_HEADER *grd, int *pad, int bilinear);
EXTERN_MSC double GMT_get_bcr_z (struct GRD_HEADER *grd, double xx, double yy, float *data, struct GMT_EDGEINFO *edgeinfo); /* Compute z(x,y) from bcr structure */
/*----------------------------------------------------------------
Here are some more remarks:
Define x,y on the bcr so that they are in the range [0,1). For pixel grids
with points close to edges, x,y will have to be in [-0.5,1) or [0, 1.5]. Now the
rectangle has 4 vertices, which we number thus:
vertex 0: x=0, y=0;
vertex 1: x=1, y=0;
vertex 2: x=0, y=1;
vertex 3: x=1, y=1.
i,j in struct BCR refers to the location of vertex 0. The i,j will be referred to
the struct GRD header values, not to the location in the padded array, so that i=0,
j=0 is the point h.x_min, h.y_max.
The nodal values are specified at each vertex as nodal_value[vertex][value], where
value stands for the following:
value 0: z;
value 1: dz/dx;
value 2: dzdy;
value 3: d2zdxdy.
If we want a bilinear estimate of z(x,y), only nodal_value[vertex][0] is needed
and the estimate is obtained thus:
z = 0.0;
for (vertex = 0; vertex < 4; vertex++)
z += bl_basis[vertex] * nodal_value[vertex][0];
This means that bl_basis[i] is the function that gives 1 at vertex i and 0 at
the other three vertices:
bl_basis[0] = (1.0 - x)*(1.0 - y);
bl_basis[1] = x*(1.0 - y);
bl_basis[2] = y*(1.0 - x);
bl_basis[3] = x*y;
If we want a bicubic surface, then there are sixteen multiply-add operations,
looping over vertex and over value. There are 16 basis functions, which are more
complicated. See the table in Lancaster and Salkauskas or the source code for
GMT_get_bcr_cardinals().
Given a point xx,yy in user's units, we need to be able to find i,j for
that point and x,y. This will depend on struct GRD_HEADER information.
If i,j for x,y is not equal to the last i,j we need to recompute the
nodal_values. If the current i,j is within +/- 1 of the last i,j we can save
some old nodal_values. If we discover a NaN during the computation of these,
we cannot continue.
If nan_condition is false, we can evaluate the cardinals and do the
multiply-adds. If x=0 or 1 or y = 0 or 1, there may be a trivial solution.
*/