home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / gmt_bcr.h < prev    next >
C/C++ Source or Header  |  1995-02-04  |  4KB  |  105 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)gmt_bcr.h    3.2  2/4/95
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /* bcr.h  --  stuff for implementing bicubic rectangle grid interpolation.
  8.    Has bilinear case as a subset.
  9.    The bilinear interpolation will give continuity of the function z only.
  10.    The bicubic interpolation will give continuity of z, dz/dx, dz/dy, and
  11.    d2z/dxdy.
  12.    Designed to operate with at least one spare row/column around each edge,
  13.    to allow derivative computations.
  14.    Follows the outline in Lancaster and Salkauskas, section 9.3
  15.    Meant to replace Taylor routines in GMT in version 3.0 of GMT.
  16.  
  17.    Author:    Walter H F Smith
  18.    Date:    23 September, 1993
  19.  
  20.    This include file defines structures and functions used.
  21. */
  22.  
  23. struct BCR {
  24.     double    nodal_value[4][4];    /* z, dz/dx, dz/dy, d2z/dxdy at 4 corners  */
  25.     double    bcr_basis[4][4];    /* multiply on nodal vals, yields z at point */
  26.     double    bl_basis[4];        /* bilinear basis functions  */
  27.     double    rx_inc;            /* 1.0 / grd.x_inc  */
  28.     double    ry_inc;            /* 1.0 / grd.y_inc  */
  29.     double    offset;            /* 0 or 0.5 for grid or pixel registration  */
  30. /* If we later want to estimate of dz/dx or dz/dy, we will need [4][4] basis for these  */
  31.     int    ij_move[4];        /* add to ij of zero vertex to get other vertex ij  */
  32.     int    i;            /* Location of current nodal_values  */
  33.     int    j;            /* Ditto.   */
  34.     int    bilinear;        /* T/F use bilinear instead of bicubic  */
  35.     int    nan_condition;        /* T/F we cannot evaluate; return z = NaN  */
  36.     int    ioff;            /* Padding on west side of array  */
  37.     int    joff;            /* Padding on north side of array  */
  38.     int    mx;            /* Padded array dimension  */
  39.     int    my;            /* Ditto  */
  40. };
  41.  
  42. /* extern struct BCR bcr;  */
  43.  
  44. struct BCR bcr;
  45.  
  46. void    bcr_init();        /* Initialize bcr structure.  */
  47. void    get_bcr_cardinals();    /* Load bcr.bcr_basis[][] and bcr.bl_basis[] given x,y */
  48. void    get_bcr_ij();        /* Find bcr.i, bcr.j given xx, yy, grdheader  */
  49. void    get_bcr_xy();        /* Find x,y given xx, yy, grdheader, bcr  */
  50. void    get_bcr_nodal_values();    /* Load bcr.nodal_values[][] and check for NaN, etc.  */
  51. double    get_bcr_z();        /* Compute z(x,y) from bcr structure  */
  52.  
  53. /* Here are some more remarks:
  54.  
  55.     Define x,y on the bcr so that they are in the range [0,1).  For pixel grids
  56. with points close to edges, x,y will have to be in [-0.5,1) or [0, 1.5].  Now the
  57. rectangle has 4 vertices, which we number thus:
  58.     vertex 0:    x=0, y=0;
  59.     vertex 1:    x=1, y=0;
  60.     vertex 2:    x=0, y=1;
  61.     vertex 3:    x=1, y=1.
  62.  
  63. i,j in struct BCR refers to the location of vertex 0.  The i,j will be referred to
  64. the struct GRD header values, not to the location in the padded array, so that i=0,
  65. j=0 is the point h.x_min, h.y_max.
  66.  
  67. The nodal values are specified at each vertex as nodal_value[vertex][value], where
  68. value stands for the following:
  69.     value 0:    z;
  70.     value 1:    dz/dx;
  71.     value 2:    dzdy;
  72.     value 3:    d2zdxdy.
  73.  
  74. If we want a bilinear estimate of z(x,y), only nodal_value[vertex][0] is needed
  75. and the estimate is obtained thus:
  76.     z = 0.0;
  77.     for (vertex = 0; vertex < 4; vertex++)
  78.         z += bl_basis[vertex] * nodal_value[vertex][0];
  79.  
  80. This means that bl_basis[i] is the function that gives 1 at vertex i and 0 at
  81. the other three vertices:
  82.     bl_basis[0] = (1.0 - x)*(1.0 - y);
  83.     bl_basis[1] = x*(1.0 - y);
  84.     bl_basis[2] = y*(1.0 - x);
  85.     bl_basis[3] = x*y;
  86.  
  87.  
  88.     If we want a bicubic surface, then there are sixteen multiply-add operations,
  89. looping over vertex and over value.  There are 16 basis functions, which are more
  90. complicated.  See the table in Lancaster and Salkauskas or the source code for
  91. load_bcr_cardinals().
  92.  
  93.     Given a point xx,yy in user's units, we need to be able to find i,j for
  94. that point and x,y.  This will depend on struct GRD_HEADER information.
  95.  
  96.     If i,j for x,y is not equal to the last i,j we need to recompute the
  97. nodal_values.  If the current i,j is within +/- 1 of the last i,j we can save
  98. some old nodal_values.  If we discover a NaN during the computation of these,
  99. we cannot continue.
  100.  
  101.     If nan_condition is false, we can evaluate the cardinals and do the
  102. multiply-adds.  If x=0 or 1 or y = 0 or 1, there may be a trivial solution.
  103.  
  104. */
  105.