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