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