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