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