PLplot 5.15.0
Loading...
Searching...
No Matches
plgridd.c
Go to the documentation of this file.
1// plgridd.c: Plot grids data from irregularly sampled data.
2
3// Copyright (C) 2004 Joao Cardoso
4// Copyright (C) 2004-2015 Alan W. Irwin
5//
6// This file is part of PLplot.
7//
8// PLplot is free software; you can redistribute it and/or modify
9// it under the terms of the GNU Library General Public License as published
10// by the Free Software Foundation; either version 2 of the License, or
11// (at your option) any later version.
12//
13// PLplot is distributed in the hope that it will be useful,
14// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16// GNU Library General Public License for more details.
17//
18// You should have received a copy of the GNU Library General Public License
19// along with PLplot; if not, write to the Free Software
20// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21//
22//
23
24#include "plplotP.h"
25
26#ifdef WITH_CSA
27#include "../lib/csa/csa.h"
28#endif
29#include "../lib/csa/nan.h" // this is handy
30
31#ifdef PL_HAVE_QHULL
32#include "../lib/nn/nn.h"
33// PLPLOT_NONN not #defined or tested for more than a decade.
34#ifdef PLPLOT_NONN // another DTLI, based only on QHULL, not nn
35#ifdef HAS_LIBQHULL_INCLUDE
36#include <libqhull/qhull_a.h>
37#else //#ifdef HAS_LIBQHULL_INCLUDE
38#include <qhull/qhull_a.h>
39#endif //#ifdef HAS_LIBQHULL_INCLUDE
40#endif //#ifdef PLPLOT_NONN
41#endif //#ifdef PL_HAVE_QHUL
42
43// forward declarations
44static void
46 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy,
47 PLF2OPS zops, PLPointer zgp );
48
49static void
51 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy,
52 PLF2OPS zops, PLPointer zgp, PLFLT threshold );
53
54static void
56 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy,
57 PLF2OPS zops, PLPointer zgp, int knn_order );
58
59#ifdef WITH_CSA
60static void
61grid_csa( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts,
62 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy,
63 PLF2OPS zops, PLPointer zgp );
64#endif
65
66#ifdef PL_HAVE_QHULL
67static void
68grid_nni( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts,
69 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy,
70 PLF2OPS zops, PLPointer zgp, PLFLT wtmin );
71
72static void
73grid_dtli( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts,
74 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy,
75 PLF2OPS zops, PLPointer zgp );
76#endif
77
78static void
79dist1( PLFLT gx, PLFLT gy, PLFLT_VECTOR x, PLFLT_VECTOR y, int npts, int knn_order );
80static void
81dist2( PLFLT gx, PLFLT gy, PLFLT_VECTOR x, PLFLT_VECTOR y, int npts );
82
83#define KNN_MAX_ORDER 100
84
85typedef struct pt
86{
88 int item;
90
92
93//--------------------------------------------------------------------------
94//
95// plgriddata(): grids data from irregularly sampled data.
96//
97// Real world data is frequently irregularly sampled, but most 3D plots
98// require regularly gridded data. This function does exactly this
99// using several methods:
100// Irregularly sampled data x[npts], y[npts], z[npts] is gridded into
101// zg[nptsx, nptsy] according to methode 'type' and grid information
102// xg[nptsx], yg[nptsy].
103//
104// 'type' can be:
105//
106// GRID_CSA: Bivariate Cubic Spline approximation (1)
107// GRID_NNIDW: Nearest Neighbors Inverse Distance Weighted
108// GRID_NNLI: Nearest Neighbors Linear Interpolation
109// GRID_NNAIDW: Nearest Neighbors Around Inverse Distance Weighted
110// GRID_DTLI: Delaunay Triangulation Linear Interpolation (2)
111// GRID_NNI: Natural Neighbors interpolation (2)
112//
113// (1): Copyright 2000-2002 CSIRO Marine Research, Pavel Sakov's csa library
114// (2): Copyright 2002 CSIRO Marine Research, Pavel Sakov's nn library
115//
116//--------------------------------------------------------------------------
117
118void
120 PLFLT_VECTOR xg, PLINT nptsx, PLFLT_VECTOR yg, PLINT nptsy,
121 PLFLT **zg, PLINT type, PLFLT data )
122{
123 plfgriddata( x, y, z, npts, xg, nptsx, yg, nptsy, plf2ops_c(), (PLPointer) zg, type, data );
124}
125
126void
128 PLFLT_VECTOR xg, PLINT nptsx, PLFLT_VECTOR yg, PLINT nptsy,
129 PLF2OPS zops, PLPointer zgp, PLINT type, PLFLT data )
130{
131 int i, j;
132
133 if ( npts < 1 || nptsx < 1 || nptsy < 1 )
134 {
135 plabort( "plgriddata: Bad array dimensions" );
136 return;
137 }
138
139 // Check that points in xg and in yg are strictly increasing
140
141 for ( i = 0; i < nptsx - 1; i++ )
142 {
143 if ( xg[i] >= xg[i + 1] )
144 {
145 plabort( "plgriddata: xg array must be strictly increasing" );
146 return;
147 }
148 }
149 for ( i = 0; i < nptsy - 1; i++ )
150 {
151 if ( yg[i] >= yg[i + 1] )
152 {
153 plabort( "plgriddata: yg array must be strictly increasing" );
154 return;
155 }
156 }
157
158 // clear array to return
159 for ( i = 0; i < nptsx; i++ )
160 for ( j = 0; j < nptsy; j++ )
161 zops->set( zgp, i, j, 0.0 );
162 // NaN signals a not processed grid point
163
164 switch ( type )
165 {
166 case ( GRID_CSA ): // Bivariate Cubic Spline Approximation
167#ifdef WITH_CSA
168 grid_csa( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
169#else
170 plwarn( "plgriddata(): PLplot was configured to not use GRID_CSA.\n Reverting to GRID_NNAIDW." );
171 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
172#endif
173 break;
174
175 case ( GRID_NNIDW ): // Nearest Neighbors Inverse Distance Weighted
176 grid_nnidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, (int) data );
177 break;
178
179 case ( GRID_NNLI ): // Nearest Neighbors Linear Interpolation
180 grid_nnli( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, data );
181 break;
182
183 case ( GRID_NNAIDW ): // Nearest Neighbors "Around" Inverse Distance Weighted
184 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
185 break;
186
187 case ( GRID_DTLI ): // Delaunay Triangulation Linear Interpolation
188#ifdef PL_HAVE_QHULL
189 grid_dtli( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
190#else
191 plwarn( "plgriddata(): you must have the Qhull library installed to use GRID_DTLI.\n Reverting to GRID_NNAIDW." );
192 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
193#endif
194 break;
195
196 case ( GRID_NNI ): // Natural Neighbors
197#ifdef PL_HAVE_QHULL
198 grid_nni( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp, data );
199#else
200 plwarn( "plgriddata(): you must have the Qhull library installed to use GRID_NNI.\n Reverting to GRID_NNAIDW." );
201 grid_nnaidw( x, y, z, npts, xg, nptsx, yg, nptsy, zops, zgp );
202#endif
203 break;
204
205 default:
206 plabort( "plgriddata: unknown algorithm type" );
207 }
208}
209
210#ifdef WITH_CSA
211//
212// Bivariate Cubic Spline Approximation using Pavel Sakov's csa package
213//
214// NaNs are returned where no interpolation can be done.
215//
216
217static void
218grid_csa( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts,
219 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy,
220 PLF2OPS zops, PLPointer zgp )
221{
222 PLFLT_VECTOR xt, yt, zt;
223 point *pin, *pgrid, *pt;
224 csa * a = NULL;
225 int i, j, nptsg;
226
227 if ( ( pin = (point *) malloc( (size_t) npts * sizeof ( point ) ) ) == NULL )
228 {
229 plexit( "grid_csa: Insufficient memory" );
230 }
231
232 xt = x;
233 yt = y;
234 zt = z;
235 pt = pin;
236 for ( i = 0; i < npts; i++ )
237 {
238 pt->x = (double) *xt++;
239 pt->y = (double) *yt++;
240 pt->z = (double) *zt++;
241 pt++;
242 }
243
244 nptsg = nptsx * nptsy;
245 if ( ( pgrid = (point *) malloc( (size_t) nptsg * sizeof ( point ) ) ) == NULL )
246 {
247 plexit( "grid_csa: Insufficient memory" );
248 }
249
250 yt = yg;
251 pt = pgrid;
252 for ( j = 0; j < nptsy; j++ )
253 {
254 xt = xg;
255 for ( i = 0; i < nptsx; i++ )
256 {
257 pt->x = (double) *xt++;
258 pt->y = (double) *yt;
259 pt++;
260 }
261 yt++;
262 }
263
264 a = csa_create();
265 csa_addpoints( a, npts, pin );
267 csa_approximate_points( a, nptsg, pgrid );
268
269 for ( i = 0; i < nptsx; i++ )
270 {
271 for ( j = 0; j < nptsy; j++ )
272 {
273 pt = &pgrid[j * nptsx + i];
274 zops->set( zgp, i, j, (PLFLT) pt->z );
275 }
276 }
277
278 csa_destroy( a );
279 free( pin );
280 free( pgrid );
281}
282#endif // WITH_CSA
283
284// Nearest Neighbors Inverse Distance Weighted, brute force approach.
285//
286// The z value at the grid position will be the weighted average
287// of the z values of the KNN points found. The weigth is the
288// inverse squared distance between the grid point and each
289// neighbor.
290//
291
292static void
294 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy,
295 PLF2OPS zops, PLPointer zgp, int knn_order )
296{
297 int i, j, k;
298 PLFLT wi, nt;
299
300 if ( knn_order > KNN_MAX_ORDER )
301 {
302 plabort( "plgriddata(): GRID_NNIDW: knn_order too big" ); // make sure it is smaller that KNN_MAX_ORDER
303 return;
304 }
305
306 if ( knn_order == 0 )
307 {
308 plwarn( "plgriddata(): GRID_NNIDW: knn_order must be specified with 'data' arg. Using 15" );
309 knn_order = 15;;
310 }
311
312 for ( i = 0; i < nptsx; i++ )
313 {
314 for ( j = 0; j < nptsy; j++ )
315 {
316 dist1( xg[i], yg[j], x, y, npts, knn_order );
317
318#ifdef GMS // alternative weight coeficients. I Don't like the results
319 // find the maximum distance
320 md = items[0].dist;
321 for ( k = 1; k < knn_order; k++ )
322 if ( items[k].dist > md )
323 md = items[k].dist;
324#endif
325 zops->set( zgp, i, j, 0.0 );
326 nt = 0.;
327
328 for ( k = 0; k < knn_order; k++ )
329 {
330 if ( items[k].item == -1 ) // not enough neighbors found ?!
331 continue;
332#ifdef GMS
333 wi = ( md - items[k].dist ) / ( md * items[k].dist );
334 wi = wi * wi;
335#else
336 wi = 1. / ( items[k].dist * items[k].dist );
337#endif
338 zops->add( zgp, i, j, wi * z[items[k].item] );
339 nt += wi;
340 }
341 if ( nt != 0. )
342 zops->div( zgp, i, j, nt );
343 else
344 zops->set( zgp, i, j, NaN );
345 }
346 }
347}
348
349// Nearest Neighbors Linear Interpolation
350//
351// The z value at the grid position will be interpolated from the
352// plane passing through the 3 nearest neighbors.
353//
354
355static void
357 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy,
358 PLF2OPS zops, PLPointer zgp, PLFLT threshold )
359{
360 PLFLT xx[4], yy[4], zz[4], t, A, B, C, D, d1, d2, d3, max_thick;
361 int i, j, ii, excl, cnt, excl_item;
362
363 if ( threshold == 0. )
364 {
365 plwarn( "plgriddata(): GRID_NNLI: threshold must be specified with 'data' arg. Using 1.001" );
366 threshold = 1.001;
367 }
368 else if ( threshold > 2. || threshold < 1. )
369 {
370 plabort( "plgriddata(): GRID_NNLI: 1. < threshold < 2." );
371 return;
372 }
373
374 for ( i = 0; i < nptsx; i++ )
375 {
376 for ( j = 0; j < nptsy; j++ )
377 {
378 dist1( xg[i], yg[j], x, y, npts, 3 );
379
380 // see if the triangle is a thin one
381 for ( ii = 0; ii < 3; ii++ )
382 {
383 xx[ii] = x[items[ii].item];
384 yy[ii] = y[items[ii].item];
385 zz[ii] = z[items[ii].item];
386 }
387
388 d1 = sqrt( ( xx[1] - xx[0] ) * ( xx[1] - xx[0] ) + ( yy[1] - yy[0] ) * ( yy[1] - yy[0] ) );
389 d2 = sqrt( ( xx[2] - xx[1] ) * ( xx[2] - xx[1] ) + ( yy[2] - yy[1] ) * ( yy[2] - yy[1] ) );
390 d3 = sqrt( ( xx[0] - xx[2] ) * ( xx[0] - xx[2] ) + ( yy[0] - yy[2] ) * ( yy[0] - yy[2] ) );
391
392 if ( d1 == 0. || d2 == 0. || d3 == 0. ) // coincident points
393 {
394 zops->set( zgp, i, j, NaN );
395 continue;
396 }
397
398 // make d1 < d2
399 if ( d1 > d2 )
400 {
401 t = d1; d1 = d2; d2 = t;
402 }
403
404 // and d2 < d3
405 if ( d2 > d3 )
406 {
407 t = d2; d2 = d3; d3 = t;
408 }
409
410 if ( ( d1 + d2 ) / d3 < threshold ) // thin triangle!
411 {
412 zops->set( zgp, i, j, NaN ); // deal with it later
413 }
414 else // calculate the plane passing through the three points
415
416 {
417 A = yy[0] * ( zz[1] - zz[2] ) + yy[1] * ( zz[2] - zz[0] ) + yy[2] * ( zz[0] - zz[1] );
418 B = zz[0] * ( xx[1] - xx[2] ) + zz[1] * ( xx[2] - xx[0] ) + zz[2] * ( xx[0] - xx[1] );
419 C = xx[0] * ( yy[1] - yy[2] ) + xx[1] * ( yy[2] - yy[0] ) + xx[2] * ( yy[0] - yy[1] );
420 D = -A * xx[0] - B * yy[0] - C * zz[0];
421
422 // and interpolate (or extrapolate...)
423 zops->set( zgp, i, j, -xg[i] * A / C - yg[j] * B / C - D / C );
424 }
425 }
426 }
427
428 // now deal with NaNs resulting from thin triangles. The idea is
429 // to use the 4 KNN points and exclude one at a time, creating
430 // four triangles, evaluating their thickness and choosing the
431 // most thick as the final one from where the interpolating
432 // plane will be build. Now that I'm talking of interpolating,
433 // one should really check that the target point is interior to
434 // the candidate triangle... otherwise one is extrapolating
435 //
436
437 {
438 for ( i = 0; i < nptsx; i++ )
439 {
440 for ( j = 0; j < nptsy; j++ )
441 {
442 if ( zops->is_nan( zgp, i, j ) )
443 {
444 dist1( xg[i], yg[j], x, y, npts, 4 );
445
446 // sort by distances. Not really needed!
447 // for (ii=3; ii>0; ii--) {
448 // for (jj=0; jj<ii; jj++) {
449 // if (items[jj].dist > items[jj+1].dist) {
450 // t = items[jj].dist;
451 // items[jj].dist = items[jj+1].dist;
452 // items[jj+1].dist = t;
453 // }
454 // }
455 // }
456 //
457
458 max_thick = 0.; excl_item = -1;
459 for ( excl = 0; excl < 4; excl++ ) // the excluded point
460
461 {
462 cnt = 0;
463 for ( ii = 0; ii < 4; ii++ )
464 {
465 if ( ii != excl )
466 {
467 xx[cnt] = x[items[ii].item];
468 yy[cnt] = y[items[ii].item];
469 cnt++;
470 }
471 }
472
473 d1 = sqrt( ( xx[1] - xx[0] ) * ( xx[1] - xx[0] ) + ( yy[1] - yy[0] ) * ( yy[1] - yy[0] ) );
474 d2 = sqrt( ( xx[2] - xx[1] ) * ( xx[2] - xx[1] ) + ( yy[2] - yy[1] ) * ( yy[2] - yy[1] ) );
475 d3 = sqrt( ( xx[0] - xx[2] ) * ( xx[0] - xx[2] ) + ( yy[0] - yy[2] ) * ( yy[0] - yy[2] ) );
476 if ( d1 == 0. || d2 == 0. || d3 == 0. ) // coincident points
477 continue;
478
479 // make d1 < d2
480 if ( d1 > d2 )
481 {
482 t = d1; d1 = d2; d2 = t;
483 }
484 // and d2 < d3
485 if ( d2 > d3 )
486 {
487 t = d2; d2 = d3; d3 = t;
488 }
489
490 t = ( d1 + d2 ) / d3;
491 if ( t > max_thick )
492 {
493 max_thick = t;
494 excl_item = excl;
495 }
496 }
497
498 if ( excl_item == -1 ) // all points are coincident?
499 continue;
500
501 // one has the thicker triangle constructed from the 4 KNN
502 cnt = 0;
503 for ( ii = 0; ii < 4; ii++ )
504 {
505 if ( ii != excl_item )
506 {
507 xx[cnt] = x[items[ii].item];
508 yy[cnt] = y[items[ii].item];
509 zz[cnt] = z[items[ii].item];
510 cnt++;
511 }
512 }
513
514 A = yy[0] * ( zz[1] - zz[2] ) + yy[1] * ( zz[2] - zz[0] ) + yy[2] * ( zz[0] - zz[1] );
515 B = zz[0] * ( xx[1] - xx[2] ) + zz[1] * ( xx[2] - xx[0] ) + zz[2] * ( xx[0] - xx[1] );
516 C = xx[0] * ( yy[1] - yy[2] ) + xx[1] * ( yy[2] - yy[0] ) + xx[2] * ( yy[0] - yy[1] );
517 D = -A * xx[0] - B * yy[0] - C * zz[0];
518
519 // and interpolate (or extrapolate...)
520 zops->set( zgp, i, j, -xg[i] * A / C - yg[j] * B / C - D / C );
521 }
522 }
523 }
524 }
525}
526
527//
528// Nearest Neighbors "Around" Inverse Distance Weighted, brute force approach.
529//
530// This uses the 1-KNN in each quadrant around the grid point, then
531// Inverse Distance Weighted is used as in GRID_NNIDW.
532//
533
534static void
536 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy, PLF2OPS zops, PLPointer zgp )
537{
538 PLFLT d, nt;
539 int i, j, k;
540
541 for ( i = 0; i < nptsx; i++ )
542 {
543 for ( j = 0; j < nptsy; j++ )
544 {
545 dist2( xg[i], yg[j], x, y, npts );
546 zops->set( zgp, i, j, 0. );
547 nt = 0.;
548 for ( k = 0; k < 4; k++ )
549 {
550 if ( items[k].item != -1 ) // was found
551 {
552 d = 1. / ( items[k].dist * items[k].dist ); // 1/square distance
553 zops->add( zgp, i, j, d * z[items[k].item] );
554 nt += d;
555 }
556 }
557 if ( nt == 0. ) // no points found?!
558 zops->set( zgp, i, j, NaN );
559 else
560 zops->div( zgp, i, j, nt );
561 }
562 }
563}
564
565#ifdef PL_HAVE_QHULL
566//
567// Delaunay Triangulation Linear Interpolation using Pavel Sakov's nn package
568//
569// The Delaunay Triangulation on the data points is build and for
570// each grid point the triangle where it is enclosed found and a
571// linear interpolation performed.
572//
573// Points exterior to the convex hull of the data points cannot
574// be interpolated and are set to NaN.
575//
576
577static void
578grid_dtli( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts,
579 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy, PLF2OPS zops, PLPointer zgp )
580{
581 point *pin, *pgrid, *pt;
582 PLFLT_VECTOR xt, yt, zt;
583 int i, j, nptsg;
584
585 // Build system has already checked that sizeof ( realT ) == sizeof ( double )
586
587 if ( ( pin = (point *) malloc( (size_t) npts * sizeof ( point ) ) ) == NULL )
588 {
589 plexit( "grid_dtli: Insufficient memory" );
590 }
591
592 xt = x;
593 yt = y;
594 zt = z;
595 pt = pin;
596 for ( i = 0; i < npts; i++ )
597 {
598 pt->x = (double) *xt++;
599 pt->y = (double) *yt++;
600 pt->z = (double) *zt++;
601 pt++;
602 }
603
604 nptsg = nptsx * nptsy;
605
606 if ( ( pgrid = (point *) malloc( (size_t) nptsg * sizeof ( point ) ) ) == NULL )
607 {
608 plexit( "grid_dtli: Insufficient memory" );
609 }
610
611 yt = yg;
612 pt = pgrid;
613 for ( j = 0; j < nptsy; j++ )
614 {
615 xt = xg;
616 for ( i = 0; i < nptsx; i++ )
617 {
618 pt->x = (double) *xt++;
619 pt->y = (double) *yt;
620 pt++;
621 }
622 yt++;
623 }
624
625 lpi_interpolate_points( npts, pin, nptsg, pgrid );
626 for ( i = 0; i < nptsx; i++ )
627 {
628 for ( j = 0; j < nptsy; j++ )
629 {
630 pt = &pgrid[j * nptsx + i];
631 zops->set( zgp, i, j, (PLFLT) pt->z );
632 }
633 }
634
635 free( pin );
636 free( pgrid );
637}
638
639//
640// Natural Neighbors using Pavel Sakov's nn package
641//
642// Points exterior to the convex hull of the data points cannot
643// be interpolated and are set to NaN.
644//
645
646static void
647grid_nni( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts,
648 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy, PLF2OPS zops, PLPointer zgp,
649 PLFLT wtmin )
650{
651 PLFLT_VECTOR xt, yt, zt;
652 point *pin, *pgrid, *pt;
653 int i, j, nptsg;
655
656 // Build system has already checked that sizeof ( realT ) == sizeof ( double )
657
658 if ( wtmin == 0. ) // only accept weights greater than wtmin
659 {
660 plwarn( "plgriddata(): GRID_NNI: wtmin must be specified with 'data' arg. Using -PLFLT_MAX" );
661 wtmin = -PLFLT_MAX;
662 }
663
664 if ( ( pin = (point *) malloc( (size_t) npts * sizeof ( point ) ) ) == NULL )
665 {
666 plexit( "plgridata: Insufficient memory" );
667 }
668
669 xt = x;
670 yt = y;
671 zt = z;
672 pt = pin;
673 for ( i = 0; i < npts; i++ )
674 {
675 pt->x = (double) *xt++;
676 pt->y = (double) *yt++;
677 pt->z = (double) *zt++;
678 pt++;
679 }
680
681 nptsg = nptsx * nptsy;
682
683 if ( ( pgrid = (point *) malloc( (size_t) nptsg * sizeof ( point ) ) ) == NULL )
684 {
685 plexit( "plgridata: Insufficient memory" );
686 }
687
688 yt = yg;
689 pt = pgrid;
690 for ( j = 0; j < nptsy; j++ )
691 {
692 xt = xg;
693 for ( i = 0; i < nptsx; i++ )
694 {
695 pt->x = (double) *xt++;
696 pt->y = (double) *yt;
697 pt++;
698 }
699 yt++;
700 }
701
702 nnpi_interpolate_points( npts, pin, wtmin, nptsg, pgrid );
703 for ( i = 0; i < nptsx; i++ )
704 {
705 for ( j = 0; j < nptsy; j++ )
706 {
707 pt = &pgrid[j * nptsx + i];
708 zops->set( zgp, i, j, (PLFLT) pt->z );
709 }
710 }
711
712 free( pin );
713 free( pgrid );
714}
715#endif // PL_HAVE_QHULL
716
717//
718// this function just calculates the K Nearest Neighbors of grid point
719// [gx, gy].
720//
721
722static void
723dist1( PLFLT gx, PLFLT gy, PLFLT_VECTOR x, PLFLT_VECTOR y, int npts, int knn_order )
724{
725 PLFLT d, max_dist;
726 int max_slot, i, j;
727
728 max_dist = PLFLT_MAX;
729 max_slot = 0;
730
731 for ( i = 0; i < knn_order; i++ )
732 {
733 items[i].dist = PLFLT_MAX;
734 items[i].item = -1;
735 }
736
737 for ( i = 0; i < npts; i++ )
738 {
739 d = ( ( gx - x[i] ) * ( gx - x[i] ) + ( gy - y[i] ) * ( gy - y[i] ) ); // save sqrt() time
740
741 if ( d < max_dist )
742 {
743 // found an item with a distance smaller than the
744 // maximum distance found so far. Replace.
745 //
746
747 items[max_slot].dist = d;
748 items[max_slot].item = i;
749
750 // find new maximum distance
751 max_dist = items[0].dist;
752 max_slot = 0;
753 for ( j = 1; j < knn_order; j++ )
754 {
755 if ( items[j].dist > max_dist )
756 {
757 max_dist = items[j].dist;
758 max_slot = j;
759 }
760 }
761 }
762 }
763 for ( j = 0; j < knn_order; j++ )
764 items[j].dist = sqrt( items[j].dist ); // now calculate the distance
765}
766
767//
768// This function searchs the 1-nearest neighbor in each quadrant around
769// the grid point.
770//
771
772static void
773dist2( PLFLT gx, PLFLT gy, PLFLT_VECTOR x, PLFLT_VECTOR y, int npts )
774{
775 PLFLT d;
776 int i, quad;
777
778 for ( i = 0; i < 4; i++ )
779 {
780 items[i].dist = PLFLT_MAX;
781 items[i].item = -1;
782 }
783
784 for ( i = 0; i < npts; i++ )
785 {
786 d = ( ( gx - x[i] ) * ( gx - x[i] ) + ( gy - y[i] ) * ( gy - y[i] ) ); // save sqrt() time
787
788 // trick to quickly compute a quadrant. The determined quadrants will be
789 // miss-assigned, i.e., 1->2, 2->0, 3->1, 4->3, but that is not important,
790 // speed is.
791
792 quad = 2 * ( x[i] > gx ) + ( y[i] < gy );
793
794 // try to use the octants around the grid point, as it will give smoother
795 // (and slower) results.
796 // Hint: use the quadrant info plus x[i]/y[i] to determine the octant
797
798 if ( d < items[quad].dist )
799 {
800 items[quad].dist = d;
801 items[quad].item = i;
802 }
803 }
804
805 for ( i = 0; i < 4; i++ )
806 if ( items[i].item != -1 )
807 items[i].dist = sqrt( items[i].dist );
808 // now calculate the distance
809}
810
811#ifdef PLPLOT_NONN // another DTLI, based only on QHULL, not nn
812static void
813grid_adtli( PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts,
814 PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy, PLF2OPS zops, PLPointer zgp )
815{
816 coordT *points; // array of coordinates for each point
817 boolT ismalloc = False; // True if qhull should free points
818 char flags[250]; // option flags for qhull
819 facetT *facet; // set by FORALLfacets
820 vertexT *vertex, **vertexp;
821 facetT *neighbor, **neighborp;
822 int curlong, totlong; // memory remaining after qh_memfreeshort
823 FILE *outfile = NULL;
824 FILE *errfile = stderr; // error messages from qhull code
825
826 int exitcode;
827 int i, j, k, l;
828 int dim = 2;
829 PLFLT xt[3], yt[3], zt[3];
830 PLFLT A, B, C, D;
831 coordT point[3];
832 boolT isoutside;
833 realT bestdist;
834 int totpart = 0;
835 int numfacets, numsimplicial, numridges;
836 int totneighbors, numcoplanars, numtricoplanars;
837
838 plwarn( "plgriddata: GRID_DTLI, If you have QHull knowledge, FIXME." );
839
840 // Could pass extra args to qhull through the 'data' argument of
841 // plgriddata()
842 strcpy( flags, "qhull d Qbb Qt", 250 );
843
844 if ( ( points = (coordT *) malloc( npts * ( dim + 1 ) * sizeof ( coordT ) ) ) == NULL )
845 {
846 plexit( "grid_adtli: Insufficient memory" );
847 }
848
849 for ( i = 0; i < npts; i++ )
850 {
851 points[i * dim] = x[i];
852 points[i * dim + 1] = y[i];
853 }
854
855#if 1 // easy way
856 exitcode = qh_new_qhull( dim, npts, points, ismalloc,
857 flags, outfile, errfile );
858#else
859 qh_init_A( stdin, stdout, stderr, 0, NULL );
860 exitcode = setjmp( qh errexit );
861 if ( !exitcode )
862 {
863 qh_initflags( flags );
864 qh PROJECTdelaunay = True;
865 qh_init_B( points, npts, dim, ismalloc );
866 qh_qhull();
867 }
868#endif
869 if ( !exitcode ) // if no error
870
871 {
872#if 0 // print the triangles vertices
873 printf( "Triangles\n" );
874 FORALLfacets {
875 if ( !facet->upperdelaunay )
876 {
877 FOREACHvertex_( facet->vertices )
878 printf( " %d", qh_pointid( vertex->point ) ); // vertices index
879 printf( "\n" );
880 }
881 }
882#endif
883
884#if 0 // print each triangle neighbors
885 printf( "Neigbors\n" );
886
887 qh_findgood_all( qh facet_list );
888 qh_countfacets( qh facet_list, NULL, !qh_ALL, &numfacets, &numsimplicial,
889 &totneighbors, &numridges, &numcoplanars, &numtricoplanars );
890
891 FORALLfacets {
892 if ( !facet->upperdelaunay )
893 {
894 FOREACHneighbor_( facet )
895 printf( " %d", neighbor->visitid ? neighbor->visitid - 1 : -neighbor->id );
896 printf( "\n" );
897 }
898 }
899#endif
900
901 // Without the setjmp(), Qhull will exit() after reporting an error
902 exitcode = setjmp( qh errexit );
903 if ( !exitcode )
904 {
905 qh NOerrexit = False;
906 for ( i = 0; i < nptsx; i++ )
907 for ( j = 0; j < nptsy; j++ )
908 {
909 l = 0;
910 point[0] = xg[i];
911 point[1] = yg[j];
912 qh_setdelaunay( 3, 1, point );
913
914
915 // several ways to find the triangle given a point follow.
916 // None but brute force works
917#if 0
918 facet = qh_findbestfacet( point, qh_ALL, &bestdist, &isoutside );
919#endif
920
921#if 0
922 facet = qh_findbest( point, qh facet_list, qh_ALL,
923 !qh_ISnewfacets, //qh_ALL
924 qh_NOupper,
925 &bestdist, &isoutside, &totpart );
926#endif
927
928#if 0
929 vertex = qh_nearvertex( facet, point, &bestdist );
930#endif
931
932 // Until someone implements a working qh_findbestfacet(),
933 // do an exautive search!
934 //
935 // As far as I understand it, qh_findbestfacet() and
936 // qh_findbest() fails when 'point' does not belongs to
937 // the convex hull, i.e., when the search becomes blocked
938 // when a facet is upperdelaunay (although the error
939 // message says that the facet may be upperdelaynay or
940 // flipped, I never found a flipped one).
941 //
942 // Another possibility is to implement the 'walking
943 // triangle algorithm
944
945 facet = qh_findfacet_all( point, &bestdist, &isoutside, &totpart );
946
947 if ( facet->upperdelaunay )
948 zops->set( zgp, i, j, NaN );
949 else
950 {
951 FOREACHvertex_( facet->vertices )
952 {
953 k = qh_pointid( vertex->point );
954 xt[l] = x[k];
955 yt[l] = y[k];
956 zt[l] = z[k];
957 l++;
958 }
959
960 // calculate the plane passing through the three points
961
962 A = yt[0] * ( zt[1] - zt[2] ) + yt[1] * ( zt[2] - zt[0] ) + yt[2] * ( zt[0] - zt[1] );
963 B = zt[0] * ( xt[1] - xt[2] ) + zt[1] * ( xt[2] - xt[0] ) + zt[2] * ( xt[0] - xt[1] );
964 C = xt[0] * ( yt[1] - yt[2] ) + xt[1] * ( yt[2] - yt[0] ) + xt[2] * ( yt[0] - yt[1] );
965 D = -A * xt[0] - B * yt[0] - C * zt[0];
966
967 // and interpolate
968 zops->set( zgp, i, j, -xg[i] * A / C - yg[j] * B / C - D / C );
969 }
970 }
971 }
972 qh NOerrexit = True;
973 }
974
975 free( points );
976 qh_freeqhull( !qh_ALL ); // free long memory
977 qh_memfreeshort( &curlong, &totlong ); // free short memory and memory allocator
978 if ( curlong || totlong )
979 fprintf( errfile,
980 "qhull: did not free %d bytes of long memory (%d pieces)\n",
981 totlong, curlong );
982}
983#endif // PLPLOT_NONN
#define NaN
Definition csa/nan.h:62
void csa_addpoints(csa *a, int n, point points[])
Definition csa.c:398
void csa_calculatespline(csa *a)
Definition csa.c:1665
void csa_approximate_points(csa *a, int n, point *points)
Definition csa.c:1747
csa * csa_create()
Definition csa.c:351
void csa_destroy(csa *a)
Definition csa.c:380
void lpi_interpolate_points(int nin, point pin[], int nout, point pout[])
Definition lpi.c:135
@ NON_SIBSONIAN
Definition nn.h:23
NNDLLIMPEXP void nnpi_interpolate_points(int nin, point pin[], double wmin, int nout, point pout[])
Definition nnpi.c:370
NN_RULE nn_rule
Definition nncommon.c:45
void plwarn(PLCHAR_VECTOR errormsg)
Definition plctrl.c:1863
void plexit(PLCHAR_VECTOR errormsg)
Definition plctrl.c:1958
void plabort(PLCHAR_VECTOR errormsg)
Definition plctrl.c:1894
PLF2OPS plf2ops_c()
Definition plf2ops.c:126
struct pt PT
void plfgriddata(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, PLINT npts, PLFLT_VECTOR xg, PLINT nptsx, PLFLT_VECTOR yg, PLINT nptsy, PLF2OPS zops, PLPointer zgp, PLINT type, PLFLT data)
Definition plgridd.c:127
static void dist2(PLFLT gx, PLFLT gy, PLFLT_VECTOR x, PLFLT_VECTOR y, int npts)
Definition plgridd.c:773
static void grid_nnli(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts, PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy, PLF2OPS zops, PLPointer zgp, PLFLT threshold)
Definition plgridd.c:356
static PT items[KNN_MAX_ORDER]
Definition plgridd.c:91
void c_plgriddata(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, PLINT npts, PLFLT_VECTOR xg, PLINT nptsx, PLFLT_VECTOR yg, PLINT nptsy, PLFLT **zg, PLINT type, PLFLT data)
Definition plgridd.c:119
static void grid_nnidw(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts, PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy, PLF2OPS zops, PLPointer zgp, int knn_order)
Definition plgridd.c:293
static void grid_nnaidw(PLFLT_VECTOR x, PLFLT_VECTOR y, PLFLT_VECTOR z, int npts, PLFLT_VECTOR xg, int nptsx, PLFLT_VECTOR yg, int nptsy, PLF2OPS zops, PLPointer zgp)
Definition plgridd.c:535
#define KNN_MAX_ORDER
Definition plgridd.c:83
static void dist1(PLFLT gx, PLFLT gy, PLFLT_VECTOR x, PLFLT_VECTOR y, int npts, int knn_order)
Definition plgridd.c:723
float PLFLT
Definition plplot.h:163
#define GRID_NNI
Definition plplot.h:1196
const PLFLT * PLFLT_VECTOR
Definition plplot.h:244
#define GRID_CSA
Definition plplot.h:1194
#define PLFLT_MAX
Definition plplot.h:164
#define GRID_DTLI
Definition plplot.h:1195
#define GRID_NNAIDW
Definition plplot.h:1199
int PLINT
Definition plplot.h:181
#define GRID_NNLI
Definition plplot.h:1198
void * PLPointer
Definition plplot.h:209
#define GRID_NNIDW
Definition plplot.h:1197
plgriddata(x, y, z, xg, yg, type, data)\n\ \n\ \n\ This function is used in example 21.\n\ \n\ \n\ \n\ SYNOPSIS:\n\ \n\ plgriddata(x, y, z, npts, xg, nptsx, yg, nptsy, zg, type, data)\n\ \n\ ARGUMENTS:\n\ \n\ x(PLFLT_VECTOR, input) : The input x vector.\n\ \n\ y(PLFLT_VECTOR, input) : The input y vector.\n\ \n\ z(PLFLT_VECTOR, input) : The input z vector. Each triple x[i],\n\ y[i], z[i] represents one data sample coordinate.\n\ \n\ npts(PLINT, input) : The number of data samples in the x, y and z\n\ vectors.\n\ \n\ xg(PLFLT_VECTOR, input) : A vector that specifies the grid spacing\n\ in the x direction. Usually xg has nptsx equally spaced values\n\ from the minimum to the maximum values of the x input vector.\n\ \n\ nptsx(PLINT, input) : The number of points in the xg vector.\n\ \n\ yg(PLFLT_VECTOR, input) : A vector that specifies the grid spacing\n\ in the y direction. Similar to the xg parameter.\n\ \n\ nptsy(PLINT, input) : The number of points in the yg vector.\n\ \n\ zg(PLFLT_NC_MATRIX, output) : The matrix of interpolated results\n\ where data lies in the grid specified by xg and yg. Therefore the\n\ zg matrix must be dimensioned\n\ nptsx by\n\ nptsy.\n\ \n\ type(PLINT, input) : The type of grid interpolation algorithm to\n\ use, which can be:GRID_CSA:Bivariate Cubic Spline approximation\n\ GRID_DTLI:Delaunay Triangulation Linear Interpolation\n\ GRID_NNI:Natural Neighbors Interpolation\n\ GRID_NNIDW:Nearest Neighbors Inverse Distance Weighted\n\ GRID_NNLI:Nearest Neighbors Linear Interpolation\n\ GRID_NNAIDW: Nearest Neighbors Around Inverse Distance\n\ Weighted\n\ For details of the algorithms read the source file plgridd.c.\n\ \n\ data(PLFLT, input) : Some gridding algorithms require extra data,\n\ which can be specified through this argument. Currently, for\n\ algorithm:GRID_NNIDW, data specifies the number of neighbors to\n\ use, the lower the value, the noisier(more local) the\n\ approximation is.\n\ GRID_NNLI, data specifies what a thin triangle is, in the\n\ range[1. .. 2.]. High values enable the usage of very thin\n\ triangles for interpolation, possibly resulting in error in\n\ the approximation.\n\ GRID_NNI, only weights greater than data will be accepted. If\n\ 0, all weights will be accepted.\n\ " zg
static const char shade or gradient plots n n or n gradient plots(See pllegend for similar functionality for creating\n\ legends with discrete elements). The arguments of plcolorbar provide\n\ control over the location and size of the color bar as well as the\n\ location and characteristics of the elements(most of which are\n\ optional) within that color bar. The resulting color bar is clipped\n\ at the boundaries of the current subpage.(N.B. the adopted coordinate\n\ system used for some of the parameters is defined in the documentation\n\ of the position parameter.)\n\ \n\ Redacted form reads the desired grid location from the input vectors n xg[nptsx] and yg[nptsy]
Definition csa.c:99
PLINT(* is_nan)(PLPointer p, PLINT ix, PLINT iy)
Definition plplot.h:609
PLFLT(* add)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)
Definition plplot.h:605
PLFLT(* div)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)
Definition plplot.h:608
PLFLT(* set)(PLPointer p, PLINT ix, PLINT iy, PLFLT z)
Definition plplot.h:604
Definition csa.h:30
Definition plgridd.c:86
PLFLT dist
Definition plgridd.c:87
int item
Definition plgridd.c:88