1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92

/*
* The GMTsystem: @(#)gmt_bcr.h 3.13 06/17/99
*
* Copyright (c) 19911999 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 multiplyadd 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
multiplyadds. If x=0 or 1 or y = 0 or 1, there may be a trivial solution.
*/
