Back to index

plt-scheme  4.2.1
Classes | Defines | Typedefs | Functions | Variables
plgridd.c File Reference
#include "plplotP.h"
#include "nan.h"

Go to the source code of this file.

Classes

struct  pt

Defines

#define isnan(x)   ((x) != (x))
#define KNN_MAX_ORDER   100

Typedefs

typedef struct pt PT

Functions

static void grid_nnaidw (PLFLT *x, PLFLT *y, PLFLT *z, int npts, PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg)
static void grid_nnli (PLFLT *x, PLFLT *y, PLFLT *z, int npts, PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg, PLFLT threshold)
static void grid_nnidw (PLFLT *x, PLFLT *y, PLFLT *z, int npts, PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg, int knn_order)
static void dist1 (PLFLT gx, PLFLT gy, PLFLT *x, PLFLT *y, int npts, int knn_order)
static void dist2 (PLFLT gx, PLFLT gy, PLFLT *x, PLFLT *y, int npts)
void c_plgriddata (PLFLT *x, PLFLT *y, PLFLT *z, int npts, PLFLT *xg, int nptsx, PLFLT *yg, int nptsy, PLFLT **zg, int type, PLFLT data)

Variables

static PT items [KNN_MAX_ORDER]

Class Documentation

struct pt

Definition at line 57 of file plgridd.c.

Class Members
PLFLT dist
int item

Define Documentation

#define isnan (   x)    ((x) != (x))

Definition at line 15 of file plgridd.c.

#define KNN_MAX_ORDER   100

Definition at line 55 of file plgridd.c.


Typedef Documentation

typedef struct pt PT

Function Documentation

void c_plgriddata ( PLFLT x,
PLFLT y,
PLFLT z,
int  npts,
PLFLT xg,
int  nptsx,
PLFLT yg,
int  nptsy,
PLFLT **  zg,
int  type,
PLFLT  data 
)

Definition at line 90 of file plgridd.c.

{
  int i, j;

  if(npts < 1 || nptsx < 1 || nptsy < 1) {
    plabort("plgriddata: Bad array dimensions"); 
    return;
  }

  /* Check that points in xg and in yg are strictly increasing */

  for (i = 0; i < nptsx - 1; i++) {
    if (xg[i] >= xg[i + 1]) {
      plabort("plgriddata: xg array must be strictly increasing");
      return;
    }
  }
  for (i = 0; i < nptsy - 1; i++) {
    if (yg[i] >= yg[i + 1]) {
      plabort("plgriddata: yg array must be strictly increasing");
      return;
    }
  }

  /* clear array to return */
  for(i=0; i<nptsx; i++)
    for(j=0; j<nptsy; j++)
      zg[i][j] = 0.; /* NaN signals a not processed grid point */

  switch (type) {
    
  case (GRID_CSA):     /*  Bivariate Cubic Spline Approximation */
#ifdef WITH_CSA
    grid_csa(x, y, z, npts, xg, nptsx, yg, nptsy, zg);
#else
    plabort("plgriddata(): PLplot was configured to not use GRID_CSA.");
#endif
    break;

  case (GRID_NNIDW): /* Nearest Neighbors Inverse Distance Weighted */
    grid_nnidw(x, y, z, npts, xg, nptsx, yg, nptsy, zg, (int) data);
    break;

  case (GRID_NNLI): /* Nearest Neighbors Linear Interpolation */
    grid_nnli(x, y, z, npts, xg, nptsx, yg, nptsy, zg, data);
    break;

  case (GRID_NNAIDW): /* Nearest Neighbors "Around" Inverse Distance Weighted */
    grid_nnaidw(x, y, z, npts, xg, nptsx, yg, nptsy, zg);
    break;

  case (GRID_DTLI): /* Delaunay Triangulation Linear Interpolation */
#ifdef HAVE_QHULL
    grid_dtli(x, y, z, npts, xg, nptsx, yg, nptsy, zg);
#else
    plabort("plgriddata(): you must have Qhull to use GRID_DTLI.");
#endif
    break;

  case (GRID_NNI): /* Natural Neighbors */
#ifdef HAVE_QHULL
    grid_nni(x, y, z, npts, xg, nptsx, yg, nptsy, zg, data);
#else
    plabort("plgriddata(): you must have Qhull to use GRID_NNI.");
#endif
    break;

  default:
    plabort("plgriddata: unknown algorithm type"); 
  }
}

Here is the call graph for this function:

static void dist1 ( PLFLT  gx,
PLFLT  gy,
PLFLT x,
PLFLT y,
int  npts,
int  knn_order 
) [static]

Definition at line 603 of file plgridd.c.

{

  PLFLT d, max_dist;
  int   max_slot, i, j;

  max_dist = PLFLT_MAX;
  max_slot = 0;

  for (i=0; i<knn_order; i++) {
    items[i].dist = PLFLT_MAX;
    items[i].item = -1;
  }

  for (i=0; i<npts; i++) {
    d = ((gx - x[i])*(gx - x[i]) + (gy - y[i])*(gy - y[i])); /* save sqrt() time */
    
    if (d < max_dist) {
      /* found an item with a distance smaller than the
       * maximum distance found so far. Replace.
       */

      items[max_slot].dist = d;
      items[max_slot].item = i;

      /* find new maximum distance */
      max_dist = items[0].dist;
      max_slot = 0;
      for (j=1; j<knn_order; j++) {
       if (items[j].dist > max_dist) {
         max_dist = items[j].dist;
         max_slot = j;
       }
      }
    }
  }
  for (j=0; j<knn_order; j++)
    items[j].dist = sqrt(items[j].dist); /* now calculate the distance */
}

Here is the caller graph for this function:

static void dist2 ( PLFLT  gx,
PLFLT  gy,
PLFLT x,
PLFLT y,
int  npts 
) [static]

Definition at line 649 of file plgridd.c.

{

  PLFLT d;
  int   i, quad;

  for (i=0; i<4; i++) {
    items[i].dist = PLFLT_MAX;
    items[i].item = -1;
  }

  for (i=0; i<npts; i++) {
    d = ((gx - x[i])*(gx - x[i]) + (gy - y[i])*(gy - y[i])); /* save sqrt() time */

    /* trick to quickly compute a quadrant. The determined quadrants will be
     * miss-assigned, i.e., 1->2, 2->0, 3->1, 4->3, but that is not important,
     * speed is. */

    quad = 2*(x[i] > gx) +  (y[i] < gy);

    /* try to use the octants around the grid point, as it will give smoother
     * (and slower) results.
     * Hint: use the quadrant info plus x[i]/y[i] to determine the octant */

    if (d < items[quad].dist) {
      items[quad].dist = d;
      items[quad].item = i;
    }
  }

  for (i=0; i<4; i++)
    if (items[i].item != -1)
      items[i].dist = sqrt(items[i].dist); /* now calculate the distance */
}

Here is the caller graph for this function:

static void grid_nnaidw ( PLFLT x,
PLFLT y,
PLFLT z,
int  npts,
PLFLT xg,
int  nptsx,
PLFLT yg,
int  nptsy,
PLFLT **  zg 
) [static]

Definition at line 446 of file plgridd.c.

{
  PLFLT d, nt;
  int i, j, k;

  for (i=0; i<nptsx; i++) {
    for (j=0; j<nptsy; j++) {
      dist2(xg[i], yg[j], x, y, npts);
      zg[i][j] = 0.;
      nt = 0.;
      for (k=0; k<4; k++) {
       if (items[k].item != -1) { /* was found */
         d = 1./(items[k].dist * items[k].dist); /* 1/square distance */
         zg[i][j] += d * z[items[k].item];
         nt += d;
       }
      }
      if (nt == 0.) /* no points found?! */
       zg[i][j] = NaN;
      else
       zg[i][j] /= nt;
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

static void grid_nnidw ( PLFLT x,
PLFLT y,
PLFLT z,
int  npts,
PLFLT xg,
int  nptsx,
PLFLT yg,
int  nptsy,
PLFLT **  zg,
int  knn_order 
) [static]

Definition at line 231 of file plgridd.c.

{
  int i, j, k;
  PLFLT wi, nt;

  if (knn_order > KNN_MAX_ORDER) {
    plabort("plgriddata(): GRID_NNIDW: knn_order too big"); /* make sure it is smaller that KNN_MAX_ORDER */
    return;
  }

  if (knn_order == 0) {
    plwarn("plgriddata(): GRID_NNIDW: knn_order must be specified with 'data' arg. Using 15");
    knn_order = 15;;
  }

  for (i=0; i<nptsx; i++) {
    for (j=0; j<nptsy; j++) {
      dist1(xg[i], yg[j], x, y, npts, knn_order);

#ifdef GMS /* alternative weight coeficients. I Don't like the results */
      /* find the maximum distance */
      md = items[0].dist;
      for (k=1; k<knn_order; k++)
       if (items[k].dist > md)
         md = items[k].dist;
#endif
      zg[i][j] = 0.;
      nt = 0.;

      for (k=0; k<knn_order; k++) {
       if(items[k].item == -1) /* not enough neighbors found ?! */
         continue;
#ifdef GMS
       wi = (md - items[k].dist)/(md * items[k].dist);
       wi = wi*wi;
#else
       wi = 1./(items[k].dist*items[k].dist);
#endif
       zg[i][j] += wi * z[items[k].item];
       nt += wi;
      }
      if (nt != 0.)
       zg[i][j] /= nt;
      else
       zg[i][j] = NaN;
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

static void grid_nnli ( PLFLT x,
PLFLT y,
PLFLT z,
int  npts,
PLFLT xg,
int  nptsx,
PLFLT yg,
int  nptsy,
PLFLT **  zg,
PLFLT  threshold 
) [static]

Definition at line 289 of file plgridd.c.

{
  PLFLT xx[4], yy[4], zz[4], t, A, B, C, D, d1, d2, d3, max_thick;
  int i, j, ii, excl, cnt, excl_item;
      
  if (threshold == 0.) {
    plwarn("plgriddata(): GRID_NNLI: threshold must be specified with 'data' arg. Using 1.001");
    threshold = 1.001;
  } else if ( threshold > 2. ||  threshold < 1.) {
    plabort("plgriddata(): GRID_NNLI: 1. < threshold < 2.");
    return;
  }

  for (i=0; i<nptsx; i++) {
    for (j=0; j<nptsy; j++) {
      dist1(xg[i], yg[j], x, y, npts, 3);

      /* see if the triangle is a thin one */
      for (ii=0; ii<3; ii++) {
       xx[ii] = x[items[ii].item];
       yy[ii] = y[items[ii].item];
       zz[ii] = z[items[ii].item];
      }

      d1 = sqrt((xx[1]-xx[0])*(xx[1]-xx[0]) + (yy[1]-yy[0])*(yy[1]-yy[0]));
      d2 = sqrt((xx[2]-xx[1])*(xx[2]-xx[1]) + (yy[2]-yy[1])*(yy[2]-yy[1]));
      d3 = sqrt((xx[0]-xx[2])*(xx[0]-xx[2]) + (yy[0]-yy[2])*(yy[0]-yy[2]));

      if (d1 == 0. || d2 == 0. || d3 == 0.) { /* coincident points */
       zg[i][j] = NaN; 
       continue;
      }

      /* make d1 < d2 */
      if (d1 > d2) {
       t = d1; d1 = d2; d2 = t;
      }

      /* and d2 < d3 */
      if (d2 > d3) {
       t = d2; d2 = d3; d3 = t;
      }

      if ((d1+d2)/d3 < threshold) { /* thin triangle! */
       zg[i][j] = NaN; /* deal with it latter */
      } else {  /* calculate the plane passing through the three points */
               
       A = yy[0]*(zz[1]-zz[2]) + yy[1]*(zz[2]-zz[0]) + yy[2]*(zz[0]-zz[1]);
       B = zz[0]*(xx[1]-xx[2]) + zz[1]*(xx[2]-xx[0]) + zz[2]*(xx[0]-xx[1]);
       C = xx[0]*(yy[1]-yy[2]) + xx[1]*(yy[2]-yy[0]) + xx[2]*(yy[0]-yy[1]);
       D = - A*xx[0] - B*yy[0] - C*zz[0];

       /* and interpolate (or extrapolate...) */
       zg[i][j] = - xg[i]*A/C - yg[j]*B/C - D/C;
      }
    }
  }

  /* now deal with NaNs resulting from thin triangles. The idea is
   * to use the 4 KNN points and exclude one at a time, creating
   * four triangles, evaluating their thickness and choosing the
   * most thick as the final one from where the interpolating
   * plane will be build.  Now that I'm talking of interpolating,
   * one should really check that the target point is interior to
   * the candidate triangle... otherwise one is extrapolating
   */

  { 
    for (i=0; i<nptsx; i++) {
      for (j=0; j<nptsy; j++) {

       if (isnan(zg[i][j])) {
         dist1(xg[i], yg[j], x, y, npts, 4);

         /* sort by distances. Not really needed! 
            for (ii=3; ii>0; ii--) {
            for (jj=0; jj<ii; jj++) {
            if (items[jj].dist > items[jj+1].dist) {
            t = items[jj].dist;
            items[jj].dist = items[jj+1].dist;
            items[jj+1].dist = t;
            }
            }
            }             
         */

         max_thick = 0.; excl_item = -1;
         for (excl=0; excl<4; excl++) { /* the excluded point */
             
           cnt = 0;
           for (ii=0; ii<4; ii++) {
             if (ii != excl) {
              xx[cnt] = x[items[ii].item];
              yy[cnt] = y[items[ii].item];
              cnt++;
             }
           }

           d1 = sqrt((xx[1]-xx[0])*(xx[1]-xx[0]) + (yy[1]-yy[0])*(yy[1]-yy[0]));
           d2 = sqrt((xx[2]-xx[1])*(xx[2]-xx[1]) + (yy[2]-yy[1])*(yy[2]-yy[1]));
           d3 = sqrt((xx[0]-xx[2])*(xx[0]-xx[2]) + (yy[0]-yy[2])*(yy[0]-yy[2]));
           if (d1 == 0. || d2 == 0. || d3 == 0.) /* coincident points */                     continue;

           /* make d1 < d2 */
           if (d1 > d2) {
             t = d1; d1 = d2; d2 = t;
           }
           /* and d2 < d3 */
           if (d2 > d3) {
             t = d2; d2 = d3; d3 = t;
           }

           t = (d1+d2)/d3;
           if ( t > max_thick) {
             max_thick = t;
             excl_item = excl;
           }
         }

         if (excl_item == -1) /* all points are coincident? */
           continue;

         /* one has the thicker triangle constructed from the 4 KNN */
         cnt = 0;
         for (ii=0; ii<4; ii++) {
           if (ii != excl_item) {
             xx[cnt] = x[items[ii].item];
             yy[cnt] = y[items[ii].item];
             zz[cnt] = z[items[ii].item];
             cnt++;
           }
         }
        
         A = yy[0]*(zz[1]-zz[2]) + yy[1]*(zz[2]-zz[0]) + yy[2]*(zz[0]-zz[1]);
         B = zz[0]*(xx[1]-xx[2]) + zz[1]*(xx[2]-xx[0]) + zz[2]*(xx[0]-xx[1]);
         C = xx[0]*(yy[1]-yy[2]) + xx[1]*(yy[2]-yy[0]) + xx[2]*(yy[0]-yy[1]);
         D = - A*xx[0] - B*yy[0] - C*zz[0];

         /* and interpolate (or extrapolate...) */
         zg[i][j] = - xg[i]*A/C - yg[j]*B/C - D/C;

       }
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:


Variable Documentation

PT items[KNN_MAX_ORDER] [static]

Definition at line 62 of file plgridd.c.