PLplot 5.15.0
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
static const char accounting for coordinate transforms n n n y1 to(\n\ x2,\n\ y2) . If a global coordinate transform is defined then the line is\n\ broken in to n segments to approximate the path. If no transform is\n\ defined then this simply acts like a call to pljoin.\n\ \n\ Redacted form reads the desired grid location from the input vectors n xg[nptsx] and yg[nptsy]
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
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