PLplot 5.9.6
|
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 */