Back to index

plt-scheme  4.2.1
plgridd.c
Go to the documentation of this file.
00001 
00002 #include "plplotP.h"
00003 
00004 #ifdef WITH_CSA
00005 #include "csa.h"
00006 #endif
00007 #include "nan.h" /* this is handy */
00008 
00009 #ifdef HAVE_QHULL
00010 #include "../lib/nn/nn.h"
00011 #include <qhull/qhull_a.h>
00012 #endif
00013 
00014 #if !defined(HAVE_ISNAN)
00015 #define isnan(x) ((x) != (x))
00016 #endif
00017 
00018 /* forward declarations */
00019 static void
00020 grid_nnaidw (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00021             PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg);
00022 
00023 static void
00024 grid_nnli (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00025           PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg,
00026           PLFLT threshold);
00027 
00028 static void
00029 grid_nnidw (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00030            PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg,
00031            int knn_order);
00032 
00033 #ifdef WITH_CSA
00034 static void
00035 grid_csa (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00036          PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg);
00037 #endif
00038 
00039 #ifdef HAVE_QHULL
00040 static void
00041 grid_nni (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00042          PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg,
00043          PLFLT wmin);
00044 
00045 static void
00046 grid_dtli (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00047           PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg);
00048 #endif
00049 
00050 static void
00051 dist1(PLFLT gx, PLFLT gy, PLFLT *x, PLFLT *y, int npts, int knn_order);
00052 static void
00053 dist2(PLFLT gx, PLFLT gy, PLFLT *x, PLFLT *y, int npts);
00054 
00055 #define KNN_MAX_ORDER 100
00056 
00057 typedef struct pt{
00058   PLFLT dist;
00059   int   item;
00060 }PT;
00061 
00062 static PT items[KNN_MAX_ORDER];
00063 
00064 /*----------------------------------------------------------------------*\
00065  *
00066  * plgriddata(): grids data from irregularly sampled data.
00067  *
00068  *    Real world data is frequently irregularly sampled, but most 3D plots
00069  *    require regularly gridded data. This function does exactly this
00070  *    using several methods:
00071  *    Irregularly sampled data x[npts], y[npts], z[npts] is gridded into
00072  *    zg[nptsx, nptsy] according to methode 'type' and grid information
00073  *    xg[nptsx], yg[nptsy].
00074  *
00075  *    'type' can be:
00076  *
00077  *       GRID_CSA:    Bivariate Cubic Spline approximation (1)
00078  *       GRID_NNIDW:  Nearest Neighbors Inverse Distance Weighted
00079  *       GRID_NNLI:   Nearest Neighbors Linear Interpolation
00080  *       GRID_NNAIDW: Nearest Neighbors Around Inverse Distance Weighted
00081  *       GRID_DTLI:   Delaunay Triangulation Linear Interpolation (2)
00082  *       GRID_NNI:    Natural Neighbors interpolation (2)
00083  * 
00084  * (1): Copyright 2000-2002 CSIRO Marine Research, Pavel Sakov's csa library
00085  * (2): Copyright 2002 CSIRO Marine Research, Pavel Sakov's nn library
00086  *
00087 \*----------------------------------------------------------------------*/
00088 
00089 void
00090 c_plgriddata(PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00091           PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy,
00092           PLFLT **zg, int type, PLFLT data)
00093 {
00094   int i, j;
00095 
00096   if(npts < 1 || nptsx < 1 || nptsy < 1) {
00097     plabort("plgriddata: Bad array dimensions"); 
00098     return;
00099   }
00100 
00101   /* Check that points in xg and in yg are strictly increasing */
00102 
00103   for (i = 0; i < nptsx - 1; i++) {
00104     if (xg[i] >= xg[i + 1]) {
00105       plabort("plgriddata: xg array must be strictly increasing");
00106       return;
00107     }
00108   }
00109   for (i = 0; i < nptsy - 1; i++) {
00110     if (yg[i] >= yg[i + 1]) {
00111       plabort("plgriddata: yg array must be strictly increasing");
00112       return;
00113     }
00114   }
00115 
00116   /* clear array to return */
00117   for(i=0; i<nptsx; i++)
00118     for(j=0; j<nptsy; j++)
00119       zg[i][j] = 0.; /* NaN signals a not processed grid point */
00120 
00121   switch (type) {
00122     
00123   case (GRID_CSA):     /*  Bivariate Cubic Spline Approximation */
00124 #ifdef WITH_CSA
00125     grid_csa(x, y, z, npts, xg, nptsx, yg, nptsy, zg);
00126 #else
00127     plabort("plgriddata(): PLplot was configured to not use GRID_CSA.");
00128 #endif
00129     break;
00130 
00131   case (GRID_NNIDW): /* Nearest Neighbors Inverse Distance Weighted */
00132     grid_nnidw(x, y, z, npts, xg, nptsx, yg, nptsy, zg, (int) data);
00133     break;
00134 
00135   case (GRID_NNLI): /* Nearest Neighbors Linear Interpolation */
00136     grid_nnli(x, y, z, npts, xg, nptsx, yg, nptsy, zg, data);
00137     break;
00138 
00139   case (GRID_NNAIDW): /* Nearest Neighbors "Around" Inverse Distance Weighted */
00140     grid_nnaidw(x, y, z, npts, xg, nptsx, yg, nptsy, zg);
00141     break;
00142 
00143   case (GRID_DTLI): /* Delaunay Triangulation Linear Interpolation */
00144 #ifdef HAVE_QHULL
00145     grid_dtli(x, y, z, npts, xg, nptsx, yg, nptsy, zg);
00146 #else
00147     plabort("plgriddata(): you must have Qhull to use GRID_DTLI.");
00148 #endif
00149     break;
00150 
00151   case (GRID_NNI): /* Natural Neighbors */
00152 #ifdef HAVE_QHULL
00153     grid_nni(x, y, z, npts, xg, nptsx, yg, nptsy, zg, data);
00154 #else
00155     plabort("plgriddata(): you must have Qhull to use GRID_NNI.");
00156 #endif
00157     break;
00158 
00159   default:
00160     plabort("plgriddata: unknown algorithm type"); 
00161   }
00162 }
00163 
00164 #ifdef WITH_CSA
00165 /* 
00166  * Bivariate Cubic Spline Approximation using Pavel Sakov's csa package
00167  *
00168  * NaNs are returned where no interpolation can be done.
00169  */
00170 
00171 static void
00172 grid_csa (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00173          PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg)
00174 {
00175   PLFLT *xt, *yt, *zt;
00176   point *pin, *pgrid, *pt;
00177   csa* a = NULL;
00178   int i, j, nptsg;
00179 
00180   pin = (point *) malloc(npts * sizeof(point));
00181 
00182   xt = x; yt = y; zt = z; pt = pin;
00183   for(i=0; i<npts; i++) {
00184     pt->x = (double) *xt++;
00185     pt->y = (double) *yt++;
00186     pt->z = (double) *zt++;
00187     pt++;
00188   }
00189 
00190   nptsg = nptsx * nptsy;
00191   pgrid =  (point *) malloc(nptsg * sizeof(point));
00192 
00193   yt = yg; pt = pgrid;
00194   for(j=0; j<nptsy; j++) {
00195     xt = xg;
00196     for(i=0; i<nptsx; i++) {
00197       pt->x = (double) *xt++;
00198       pt->y = (double) *yt;
00199       pt++;
00200     }
00201     yt++;
00202   }
00203 
00204   a = csa_create();
00205   csa_addpoints(a, npts, pin);
00206   csa_calculatespline(a);
00207   csa_approximate_points(a, nptsg, pgrid);
00208 
00209   for(i=0; i<nptsx; i++) {
00210     for(j=0; j<nptsy; j++) {
00211       pt = &pgrid[j*nptsx + i];
00212       zg[i][j] = (PLFLT) pt->z;
00213     }
00214   }
00215 
00216   csa_destroy(a);
00217   free(pin);
00218   free(pgrid); 
00219 }
00220 #endif /* WITH_CSA */
00221 
00222 /* Nearest Neighbors Inverse Distance Weighted, brute force approach.
00223  *
00224  * The z value at the grid position will be the weighted average
00225  * of the z values of the KNN points found. The weigth is the
00226  * inverse squared distance between the grid point and each
00227  * neighbor.
00228  */
00229 
00230 static void
00231 grid_nnidw (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00232            PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg,
00233            int knn_order)
00234 {
00235   int i, j, k;
00236   PLFLT wi, nt;
00237 
00238   if (knn_order > KNN_MAX_ORDER) {
00239     plabort("plgriddata(): GRID_NNIDW: knn_order too big"); /* make sure it is smaller that KNN_MAX_ORDER */
00240     return;
00241   }
00242 
00243   if (knn_order == 0) {
00244     plwarn("plgriddata(): GRID_NNIDW: knn_order must be specified with 'data' arg. Using 15");
00245     knn_order = 15;;
00246   }
00247 
00248   for (i=0; i<nptsx; i++) {
00249     for (j=0; j<nptsy; j++) {
00250       dist1(xg[i], yg[j], x, y, npts, knn_order);
00251 
00252 #ifdef GMS /* alternative weight coeficients. I Don't like the results */
00253       /* find the maximum distance */
00254       md = items[0].dist;
00255       for (k=1; k<knn_order; k++)
00256        if (items[k].dist > md)
00257          md = items[k].dist;
00258 #endif
00259       zg[i][j] = 0.;
00260       nt = 0.;
00261 
00262       for (k=0; k<knn_order; k++) {
00263        if(items[k].item == -1) /* not enough neighbors found ?! */
00264          continue;
00265 #ifdef GMS
00266        wi = (md - items[k].dist)/(md * items[k].dist);
00267        wi = wi*wi;
00268 #else
00269        wi = 1./(items[k].dist*items[k].dist);
00270 #endif
00271        zg[i][j] += wi * z[items[k].item];
00272        nt += wi;
00273       }
00274       if (nt != 0.)
00275        zg[i][j] /= nt;
00276       else
00277        zg[i][j] = NaN;
00278     }
00279   }
00280 }
00281 
00282 /* Nearest Neighbors Linear Interpolation
00283  *
00284  * The z value at the grid position will be interpolated from the
00285  * plane passing through the 3 nearest neighbors.
00286  */
00287 
00288 static void
00289 grid_nnli (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00290           PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg,
00291           PLFLT threshold)
00292 {
00293   PLFLT xx[4], yy[4], zz[4], t, A, B, C, D, d1, d2, d3, max_thick;
00294   int i, j, ii, excl, cnt, excl_item;
00295       
00296   if (threshold == 0.) {
00297     plwarn("plgriddata(): GRID_NNLI: threshold must be specified with 'data' arg. Using 1.001");
00298     threshold = 1.001;
00299   } else if ( threshold > 2. ||  threshold < 1.) {
00300     plabort("plgriddata(): GRID_NNLI: 1. < threshold < 2.");
00301     return;
00302   }
00303 
00304   for (i=0; i<nptsx; i++) {
00305     for (j=0; j<nptsy; j++) {
00306       dist1(xg[i], yg[j], x, y, npts, 3);
00307 
00308       /* see if the triangle is a thin one */
00309       for (ii=0; ii<3; ii++) {
00310        xx[ii] = x[items[ii].item];
00311        yy[ii] = y[items[ii].item];
00312        zz[ii] = z[items[ii].item];
00313       }
00314 
00315       d1 = sqrt((xx[1]-xx[0])*(xx[1]-xx[0]) + (yy[1]-yy[0])*(yy[1]-yy[0]));
00316       d2 = sqrt((xx[2]-xx[1])*(xx[2]-xx[1]) + (yy[2]-yy[1])*(yy[2]-yy[1]));
00317       d3 = sqrt((xx[0]-xx[2])*(xx[0]-xx[2]) + (yy[0]-yy[2])*(yy[0]-yy[2]));
00318 
00319       if (d1 == 0. || d2 == 0. || d3 == 0.) { /* coincident points */
00320        zg[i][j] = NaN; 
00321        continue;
00322       }
00323 
00324       /* make d1 < d2 */
00325       if (d1 > d2) {
00326        t = d1; d1 = d2; d2 = t;
00327       }
00328 
00329       /* and d2 < d3 */
00330       if (d2 > d3) {
00331        t = d2; d2 = d3; d3 = t;
00332       }
00333 
00334       if ((d1+d2)/d3 < threshold) { /* thin triangle! */
00335        zg[i][j] = NaN; /* deal with it latter */
00336       } else {  /* calculate the plane passing through the three points */
00337                
00338        A = yy[0]*(zz[1]-zz[2]) + yy[1]*(zz[2]-zz[0]) + yy[2]*(zz[0]-zz[1]);
00339        B = zz[0]*(xx[1]-xx[2]) + zz[1]*(xx[2]-xx[0]) + zz[2]*(xx[0]-xx[1]);
00340        C = xx[0]*(yy[1]-yy[2]) + xx[1]*(yy[2]-yy[0]) + xx[2]*(yy[0]-yy[1]);
00341        D = - A*xx[0] - B*yy[0] - C*zz[0];
00342 
00343        /* and interpolate (or extrapolate...) */
00344        zg[i][j] = - xg[i]*A/C - yg[j]*B/C - D/C;
00345       }
00346     }
00347   }
00348 
00349   /* now deal with NaNs resulting from thin triangles. The idea is
00350    * to use the 4 KNN points and exclude one at a time, creating
00351    * four triangles, evaluating their thickness and choosing the
00352    * most thick as the final one from where the interpolating
00353    * plane will be build.  Now that I'm talking of interpolating,
00354    * one should really check that the target point is interior to
00355    * the candidate triangle... otherwise one is extrapolating
00356    */
00357 
00358   { 
00359     for (i=0; i<nptsx; i++) {
00360       for (j=0; j<nptsy; j++) {
00361 
00362        if (isnan(zg[i][j])) {
00363          dist1(xg[i], yg[j], x, y, npts, 4);
00364 
00365          /* sort by distances. Not really needed! 
00366             for (ii=3; ii>0; ii--) {
00367             for (jj=0; jj<ii; jj++) {
00368             if (items[jj].dist > items[jj+1].dist) {
00369             t = items[jj].dist;
00370             items[jj].dist = items[jj+1].dist;
00371             items[jj+1].dist = t;
00372             }
00373             }
00374             }             
00375          */
00376 
00377          max_thick = 0.; excl_item = -1;
00378          for (excl=0; excl<4; excl++) { /* the excluded point */
00379              
00380            cnt = 0;
00381            for (ii=0; ii<4; ii++) {
00382              if (ii != excl) {
00383               xx[cnt] = x[items[ii].item];
00384               yy[cnt] = y[items[ii].item];
00385               cnt++;
00386              }
00387            }
00388 
00389            d1 = sqrt((xx[1]-xx[0])*(xx[1]-xx[0]) + (yy[1]-yy[0])*(yy[1]-yy[0]));
00390            d2 = sqrt((xx[2]-xx[1])*(xx[2]-xx[1]) + (yy[2]-yy[1])*(yy[2]-yy[1]));
00391            d3 = sqrt((xx[0]-xx[2])*(xx[0]-xx[2]) + (yy[0]-yy[2])*(yy[0]-yy[2]));
00392            if (d1 == 0. || d2 == 0. || d3 == 0.) /* coincident points */                     continue;
00393 
00394            /* make d1 < d2 */
00395            if (d1 > d2) {
00396              t = d1; d1 = d2; d2 = t;
00397            }
00398            /* and d2 < d3 */
00399            if (d2 > d3) {
00400              t = d2; d2 = d3; d3 = t;
00401            }
00402 
00403            t = (d1+d2)/d3;
00404            if ( t > max_thick) {
00405              max_thick = t;
00406              excl_item = excl;
00407            }
00408          }
00409 
00410          if (excl_item == -1) /* all points are coincident? */
00411            continue;
00412 
00413          /* one has the thicker triangle constructed from the 4 KNN */
00414          cnt = 0;
00415          for (ii=0; ii<4; ii++) {
00416            if (ii != excl_item) {
00417              xx[cnt] = x[items[ii].item];
00418              yy[cnt] = y[items[ii].item];
00419              zz[cnt] = z[items[ii].item];
00420              cnt++;
00421            }
00422          }
00423         
00424          A = yy[0]*(zz[1]-zz[2]) + yy[1]*(zz[2]-zz[0]) + yy[2]*(zz[0]-zz[1]);
00425          B = zz[0]*(xx[1]-xx[2]) + zz[1]*(xx[2]-xx[0]) + zz[2]*(xx[0]-xx[1]);
00426          C = xx[0]*(yy[1]-yy[2]) + xx[1]*(yy[2]-yy[0]) + xx[2]*(yy[0]-yy[1]);
00427          D = - A*xx[0] - B*yy[0] - C*zz[0];
00428 
00429          /* and interpolate (or extrapolate...) */
00430          zg[i][j] = - xg[i]*A/C - yg[j]*B/C - D/C;
00431 
00432        }
00433       }
00434     }
00435   }
00436 }
00437 
00438 /* 
00439  * Nearest Neighbors "Around" Inverse Distance Weighted, brute force approach.
00440  * 
00441  * This uses the 1-KNN in each quadrant around the grid point, then
00442  * Inverse Distance Weighted is used as in GRID_NNIDW.
00443  */
00444 
00445 static void
00446 grid_nnaidw (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00447             PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg)
00448 {
00449   PLFLT d, nt;
00450   int i, j, k;
00451 
00452   for (i=0; i<nptsx; i++) {
00453     for (j=0; j<nptsy; j++) {
00454       dist2(xg[i], yg[j], x, y, npts);
00455       zg[i][j] = 0.;
00456       nt = 0.;
00457       for (k=0; k<4; k++) {
00458        if (items[k].item != -1) { /* was found */
00459          d = 1./(items[k].dist * items[k].dist); /* 1/square distance */
00460          zg[i][j] += d * z[items[k].item];
00461          nt += d;
00462        }
00463       }
00464       if (nt == 0.) /* no points found?! */
00465        zg[i][j] = NaN;
00466       else
00467        zg[i][j] /= nt;
00468     }
00469   }
00470 }
00471 
00472 #ifdef HAVE_QHULL
00473 /*
00474  * Delaunay Triangulation Linear Interpolation using Pavel Sakov's nn package
00475  *
00476  * The Delaunay Triangulation on the data points is build and for
00477  * each grid point the triangle where it is enclosed found and a
00478  * linear interpolation performed.
00479  *
00480  * Points exterior to the convex hull of the data points cannot
00481  * be interpolated and are set to NaN.
00482  */
00483 
00484 static void
00485 grid_dtli(PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00486          PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg)
00487 {
00488   point *pin, *pgrid, *pt;
00489   PLFLT *xt, *yt, *zt;
00490   int i, j, nptsg;
00491 
00492   if (sizeof(realT) != sizeof(double)) {
00493     plabort("plgridata: QHull was compiled for floats instead of doubles");
00494     return;
00495   }
00496 
00497   pin = (point *) malloc(npts * sizeof(point));
00498 
00499   xt = x; yt = y; zt = z; pt = pin;
00500   for(i=0; i<npts; i++) {
00501     pt->x = (double) *xt++;
00502     pt->y = (double) *yt++;
00503     pt->z = (double) *zt++;
00504     pt++;
00505   }
00506 
00507   nptsg = nptsx * nptsy;
00508   pgrid =  (point *) malloc(nptsg * sizeof(point));
00509 
00510   yt = yg; pt = pgrid;
00511   for(j=0; j<nptsy; j++) {
00512     xt = xg;
00513     for(i=0; i<nptsx; i++) {
00514       pt->x = (double) *xt++;
00515       pt->y = (double) *yt;
00516       pt++;
00517     }
00518     yt++;
00519   }
00520 
00521   lpi_interpolate_points(npts, pin, nptsg, pgrid);
00522   for(i=0; i<nptsx; i++) {
00523     for(j=0; j<nptsy; j++) {
00524       pt = &pgrid[j*nptsx + i];
00525       zg[i][j] = (PLFLT) pt->z;
00526     }
00527   }
00528 
00529   free(pin);
00530   free(pgrid); 
00531 }
00532 
00533 /*
00534  * Natural Neighbors using Pavel Sakov's nn package
00535  *
00536  * Points exterior to the convex hull of the data points cannot
00537  * be interpolated and are set to NaN.
00538  */
00539 
00540 static void
00541 grid_nni (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00542          PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg,
00543          PLFLT wmin)
00544 {
00545   PLFLT *xt, *yt, *zt;
00546   point *pin, *pgrid, *pt;
00547   int i, j, nptsg;
00548   nn_algorithm = NON_SIBSONIAN;
00549 
00550   if (sizeof(realT) != sizeof(double)) {
00551     plabort("plgridata: QHull was compiled for floats instead of doubles");
00552     return;
00553   }
00554 
00555   if (wmin == 0.) {/* only accept weights greater than wmin */
00556    plwarn("plgriddata(): GRID_NNI: wmin must be specified with 'data' arg. Using -PLFLT_MAX"); 
00557     wmin =  -PLFLT_MAX;
00558   }
00559 
00560   pin = (point *) malloc(npts * sizeof(point));
00561 
00562   xt = x; yt = y; zt = z; pt = pin;
00563   for(i=0; i<npts; i++) {
00564     pt->x = (double) *xt++;
00565     pt->y = (double) *yt++;
00566     pt->z = (double) *zt++;
00567     pt++;
00568   }
00569 
00570   nptsg = nptsx * nptsy;
00571   pgrid =  (point *) malloc(nptsg * sizeof(point));
00572 
00573   yt = yg; pt = pgrid;
00574   for(j=0; j<nptsy; j++) {
00575     xt = xg;
00576     for(i=0; i<nptsx; i++) {
00577       pt->x = (double) *xt++;
00578       pt->y = (double) *yt;
00579       pt++;
00580     }
00581     yt++;
00582   }
00583 
00584   nnpi_interpolate_points(npts, pin, wmin, nptsg, pgrid);
00585   for(i=0; i<nptsx; i++) {
00586     for(j=0; j<nptsy; j++) {
00587       pt = &pgrid[j*nptsx + i];
00588       zg[i][j] = (PLFLT) pt->z;
00589     }
00590   }
00591 
00592   free(pin);
00593   free(pgrid); 
00594 }
00595 #endif /* HAVE_QHULL*/
00596 
00597 /* 
00598  * this function just calculates the K Nearest Neighbors of grid point
00599  * [gx, gy].
00600  */
00601 
00602 static void
00603 dist1(PLFLT gx, PLFLT gy, PLFLT *x, PLFLT *y, int npts, int knn_order)
00604 {
00605 
00606   PLFLT d, max_dist;
00607   int   max_slot, i, j;
00608 
00609   max_dist = PLFLT_MAX;
00610   max_slot = 0;
00611 
00612   for (i=0; i<knn_order; i++) {
00613     items[i].dist = PLFLT_MAX;
00614     items[i].item = -1;
00615   }
00616 
00617   for (i=0; i<npts; i++) {
00618     d = ((gx - x[i])*(gx - x[i]) + (gy - y[i])*(gy - y[i])); /* save sqrt() time */
00619     
00620     if (d < max_dist) {
00621       /* found an item with a distance smaller than the
00622        * maximum distance found so far. Replace.
00623        */
00624 
00625       items[max_slot].dist = d;
00626       items[max_slot].item = i;
00627 
00628       /* find new maximum distance */
00629       max_dist = items[0].dist;
00630       max_slot = 0;
00631       for (j=1; j<knn_order; j++) {
00632        if (items[j].dist > max_dist) {
00633          max_dist = items[j].dist;
00634          max_slot = j;
00635        }
00636       }
00637     }
00638   }
00639   for (j=0; j<knn_order; j++)
00640     items[j].dist = sqrt(items[j].dist); /* now calculate the distance */
00641 }
00642 
00643 /*
00644  * This function searchs the 1-nearest neighbor in each quadrant around
00645  * the grid point.
00646  */
00647 
00648 static void
00649 dist2(PLFLT gx, PLFLT gy, PLFLT *x, PLFLT *y, int npts)
00650 {
00651 
00652   PLFLT d;
00653   int   i, quad;
00654 
00655   for (i=0; i<4; i++) {
00656     items[i].dist = PLFLT_MAX;
00657     items[i].item = -1;
00658   }
00659 
00660   for (i=0; i<npts; i++) {
00661     d = ((gx - x[i])*(gx - x[i]) + (gy - y[i])*(gy - y[i])); /* save sqrt() time */
00662 
00663     /* trick to quickly compute a quadrant. The determined quadrants will be
00664      * miss-assigned, i.e., 1->2, 2->0, 3->1, 4->3, but that is not important,
00665      * speed is. */
00666 
00667     quad = 2*(x[i] > gx) +  (y[i] < gy);
00668 
00669     /* try to use the octants around the grid point, as it will give smoother
00670      * (and slower) results.
00671      * Hint: use the quadrant info plus x[i]/y[i] to determine the octant */
00672 
00673     if (d < items[quad].dist) {
00674       items[quad].dist = d;
00675       items[quad].item = i;
00676     }
00677   }
00678 
00679   for (i=0; i<4; i++)
00680     if (items[i].item != -1)
00681       items[i].dist = sqrt(items[i].dist); /* now calculate the distance */
00682 }
00683 
00684 #ifdef NONN /* another DTLI, based only on QHULL, not nn */ 
00685 static void
00686 grid_adtli (PLFLT *x, PLFLT *y, PLFLT *z, int npts,
00687           PLFLT *xg, int nptsx, PLFLT *yg,  int nptsy, PLFLT **zg)
00688 {
00689   coordT *points;           /* array of coordinates for each point */
00690   boolT ismalloc = False;   /* True if qhull should free points */
00691   char flags[250];          /* option flags for qhull */
00692   facetT *facet;         /* set by FORALLfacets */
00693   vertexT *vertex, **vertexp;
00694   facetT *neighbor,**neighborp;
00695   int curlong, totlong;         /* memory remaining after qh_memfreeshort */
00696   FILE *outfile = NULL;
00697   FILE *errfile = stderr;   /* error messages from qhull code */
00698 
00699   int exitcode;
00700   int i, j, k, l;
00701   int dim = 2;
00702   PLFLT xt[3], yt[3], zt[3];
00703   PLFLT A, B, C, D;
00704   coordT point[3];
00705   boolT isoutside;
00706   realT bestdist;
00707   int totpart=0;
00708   int numfacets, numsimplicial, numridges;
00709   int totneighbors, numcoplanars, numtricoplanars;       
00710       
00711   plwarn("plgriddata: GRID_DTLI, If you have QHull knowledge, FIXME.");
00712 
00713   /* Could pass extra args to qhull through the 'data' argument of
00714      plgriddata() */
00715   sprintf(flags, "qhull d Qbb Qt"); 
00716   points = (coordT *) malloc(npts * (dim+1) * sizeof(coordT));
00717 
00718   for (i=0; i<npts; i++) {
00719     points[i*dim] = x[i];
00720     points[i*dim+1] = y[i];
00721   }
00722 
00723 #if 1 /* easy way */
00724   exitcode = qh_new_qhull (dim, npts, points, ismalloc,
00725                         flags, outfile, errfile); 
00726 #else
00727   qh_init_A (stdin, stdout, stderr, 0, NULL);
00728   exitcode = setjmp (qh errexit);
00729   if (!exitcode) {
00730     qh_initflags (flags);
00731     qh PROJECTdelaunay = True;
00732     qh_init_B (points, npts, dim, ismalloc);
00733     qh_qhull();
00734   }
00735 #endif
00736   if (!exitcode) {                  /* if no error */
00737 
00738 #if 0 /* print the triangles vertices */
00739     printf("Triangles\n");
00740     FORALLfacets {
00741       if (!facet->upperdelaunay) {
00742        FOREACHvertex_(facet->vertices)
00743          printf (" %d", qh_pointid (vertex->point)); /* vertices index */
00744        printf ("\n");
00745       }
00746     }
00747 #endif
00748 
00749 #if 0 /* print each triangle neighbors */
00750     printf("Neigbors\n");
00751 
00752     qh_findgood_all (qh facet_list);
00753     qh_countfacets (qh facet_list, NULL, !qh_ALL, &numfacets, &numsimplicial, 
00754                   &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
00755 
00756     FORALLfacets {
00757       if (!facet->upperdelaunay) {
00758        FOREACHneighbor_(facet) 
00759          printf (" %d", neighbor->visitid ? neighbor->visitid - 1: - neighbor->id);
00760        printf ("\n");
00761       }
00762     }
00763 #endif
00764 
00765     /* Without the setjmp(), Qhull will exit() after reporting an error */
00766     exitcode = setjmp (qh errexit);  
00767     if (!exitcode) {
00768       qh NOerrexit= False;
00769       for(i=0; i<nptsx; i++)
00770        for(j=0; j<nptsy; j++){
00771          l = 0;
00772          point[0] = xg[i];
00773          point[1] = yg[j];
00774          qh_setdelaunay (3, 1, point);
00775 
00776 
00777          /* several ways to find the triangle given a point follow.
00778             None but brute force works */
00779 #if 0
00780          facet = qh_findbestfacet (point, qh_ALL, &bestdist, &isoutside);
00781 #endif
00782 
00783 #if 0
00784          facet = qh_findbest (point, qh facet_list, qh_ALL,
00785                             !qh_ISnewfacets, /*qh_ALL*/ qh_NOupper,
00786                             &bestdist, &isoutside, &totpart);
00787 #endif
00788        
00789 #if 0
00790          vertex = qh_nearvertex (facet, point, &bestdist);
00791 #endif
00792 
00793          /* Until someone implements a working qh_findbestfacet(),
00794           * do an exautive search!
00795           *
00796           * As far as I understand it, qh_findbestfacet() and
00797           * qh_findbest() fails when 'point' does not belongs to
00798           * the convex hull, i.e., when the search becomes blocked
00799           * when a facet is upperdelaunay (although the error
00800           * message says that the facet may be upperdelaynay or
00801           * flipped, I never found a flipped one).
00802           *
00803           * Another possibility is to implement the 'walking
00804           * triangle algorithm */
00805                
00806          facet = qh_findfacet_all( point, &bestdist, &isoutside, &totpart );
00807 
00808          if (facet->upperdelaunay)
00809            zg[i][j] = NaN;
00810          else {            
00811            FOREACHvertex_(facet->vertices) {
00812              k = qh_pointid(vertex->point);
00813              xt[l] = x[k];
00814              yt[l] = y[k];
00815              zt[l] = z[k];
00816              l++;
00817            }
00818 
00819            /* calculate the plane passing through the three points */
00820 
00821            A = yt[0]*(zt[1]-zt[2]) + yt[1]*(zt[2]-zt[0]) + yt[2]*(zt[0]-zt[1]);
00822            B = zt[0]*(xt[1]-xt[2]) + zt[1]*(xt[2]-xt[0]) + zt[2]*(xt[0]-xt[1]);
00823            C = xt[0]*(yt[1]-yt[2]) + xt[1]*(yt[2]-yt[0]) + xt[2]*(yt[0]-yt[1]);
00824            D = - A*xt[0] - B*yt[0] - C*zt[0];
00825 
00826            /* and interpolate */
00827            zg[i][j] = - xg[i]*A/C - yg[j]*B/C - D/C;
00828 
00829          } 
00830        }
00831     }
00832     qh NOerrexit= True;
00833   }
00834 
00835   free(points);
00836   qh_freeqhull(!qh_ALL);                 /* free long memory */
00837   qh_memfreeshort (&curlong, &totlong);  /* free short memory and memory allocator */
00838   if (curlong || totlong) 
00839     fprintf (errfile,
00840             "qhull: did not free %d bytes of long memory (%d pieces)\n",
00841             totlong, curlong);
00842 }
00843 #endif /* NONN */