PLplot
5.10.0
|
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