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