PLplot  5.10.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
nnpi.c
Go to the documentation of this file.
1 //--------------------------------------------------------------------------
2 //
3 // File: nnpi.c
4 //
5 // Created: 15/11/2002
6 //
7 // Author: Pavel Sakov
8 // CSIRO Marine Research
9 //
10 // Purpose: Code for:
11 // -- Natural Neighbours Point Interpolator
12 // -- Natural Neighbours Point Hashing Interpolator
13 //
14 // Description: `nnpi' -- "Natural Neighbours Point
15 // Interpolator" -- is a structure for conducting Natural
16 // Neighbours interpolation on a given data on a
17 // "point-to-point" basis. Because it involves weight
18 // calculation for each next output point, it is not
19 // particularly suitable for consequitive interpolations on
20 // the same set of observation points -- use
21 // `nnhpi' or `nnai'
22 // in these cases.
23 //
24 // `nnhpi' is a structure for
25 // conducting consequitive Natural Neighbours interpolations
26 // on a given spatial data set in a random sequence of points
27 // from a set of finite size, taking advantage of repeated
28 // interpolations in the same point. It allows to modify Z
29 // coordinate of data in between interpolations.
30 //
31 //
32 // Revisions: 01/04/2003 PS: modified nnpi_triangle_process(): for
33 // Sibson interpolation, if circle_build fails(), now a
34 // local copy of a point is moved slightly rather than the
35 // data point itself. The later approach have found leading
36 // to inconsistencies of the new point position with the
37 // earlier built triangulation.
38 //
39 //--------------------------------------------------------------------------
40 
41 #include <stdlib.h>
42 #include <stdio.h>
43 #include <limits.h>
44 #include <float.h>
45 #include <string.h>
46 #include <assert.h>
47 #include <math.h>
48 #include "nn.h"
49 #include "delaunay.h"
50 #include "nan.h"
51 #include "hash.h"
52 
53 struct nnpi
54 {
56  point * p;
57  double wmin;
58  //
59  // work variables
60  //
61  int nvertices;
63  int * vertices; // vertex indices
64  double* weights;
65  int n; // number of points processed
66 };
67 
68 int circle_build( circle* c, point* p0, point* p1, point* p2 );
69 int circle_contains( circle* c, point* p );
70 void delaunay_circles_find( delaunay* d, point* p, int* n, int** out );
71 int delaunay_xytoi( delaunay* d, point* p, int seed );
72 void nn_quit( const char* format, ... );
73 void nnpi_reset( nnpi* nn );
74 void nnpi_calculate_weights( nnpi* nn );
75 void nnpi_normalize_weights( nnpi* nn );
76 void nnpi_set_point( nnpi* nn, point* p );
77 int nnpi_get_nvertices( nnpi* nn );
78 int* nnpi_get_vertices( nnpi* nn );
79 double* nnpi_get_weights( nnpi* nn );
80 
81 #define NSTART 10
82 #define NINC 10
83 #define EPS_SHIFT 1.0e-9
84 #define N_SEARCH_TURNON 20
85 #define BIGNUMBER 1.0e+100
86 
87 #define min( x, y ) ( ( x ) < ( y ) ? ( x ) : ( y ) )
88 #define max( x, y ) ( ( x ) > ( y ) ? ( x ) : ( y ) )
89 
90 // Creates Natural Neighbours point interpolator.
91 //
92 // @param d Delaunay triangulation
93 // @return Natural Neighbours interpolation
94 //
96 {
97  nnpi* nn = malloc( sizeof ( nnpi ) );
98 
99  nn->d = d;
100  nn->wmin = -DBL_MAX;
101  nn->vertices = calloc( NSTART, sizeof ( int ) );
102  nn->weights = calloc( NSTART, sizeof ( double ) );
103  nn->nvertices = 0;
104  nn->nallocated = NSTART;
105  nn->p = NULL;
106  nn->n = 0;
107 
108  return nn;
109 }
110 
111 // Destroys Natural Neighbours point interpolator.
112 //
113 // @param nn Structure to be destroyed
114 //
115 void nnpi_destroy( nnpi* nn )
116 {
117  free( nn->weights );
118  free( nn->vertices );
119  free( nn );
120 }
121 
122 void nnpi_reset( nnpi* nn )
123 {
124  nn->nvertices = 0;
125  nn->p = NULL;
126  memset( nn->d->flags, 0, (size_t) ( nn->d->ntriangles ) * sizeof ( int ) );
127 }
128 
129 static void nnpi_add_weight( nnpi* nn, int vertex, double w )
130 {
131  int i;
132 
133  //
134  // find whether the vertex is already in the list
135  //
136  for ( i = 0; i < nn->nvertices; ++i )
137  if ( nn->vertices[i] == vertex )
138  break;
139 
140  if ( i == nn->nvertices ) // not in the list
141  { //
142  // get more memory if necessary
143  //
144  if ( nn->nvertices == nn->nallocated )
145  {
146  nn->vertices = realloc( nn->vertices, (size_t) ( nn->nallocated + NINC ) * sizeof ( int ) );
147  nn->weights = realloc( nn->weights, (size_t) ( nn->nallocated + NINC ) * sizeof ( double ) );
148  nn->nallocated += NINC;
149  }
150 
151  //
152  // add the vertex to the list
153  //
154  nn->vertices[i] = vertex;
155  nn->weights[i] = w;
156  nn->nvertices++;
157  }
158  else // in the list
159 
160  {
161  if ( nn_rule == SIBSON )
162  nn->weights[i] += w;
163  else if ( w > nn->weights[i] )
164  nn->weights[i] = w;
165  }
166 }
167 
168 static double triangle_scale_get( delaunay* d, triangle* t )
169 {
170  double x0 = d->points[t->vids[0]].x;
171  double x1 = d->points[t->vids[1]].x;
172  double x2 = d->points[t->vids[2]].x;
173  double y0 = d->points[t->vids[0]].y;
174  double y1 = d->points[t->vids[1]].y;
175  double y2 = d->points[t->vids[2]].y;
176  double xmin = min( min( x0, x1 ), x2 );
177  double xmax = max( max( x0, x1 ), x2 );
178  double ymin = min( min( y0, y1 ), y2 );
179  double ymax = max( max( y0, y1 ), y2 );
180 
181  return xmax - xmin + ymax - ymin;
182 }
183 
184 // This is a central procedure for the Natural Neighbours interpolation. It
185 // uses the Watson's algorithm for the required areas calculation and implies
186 // that the vertices of the delaunay triangulation are listed in uniform
187 // (clockwise or counterclockwise) order.
188 //
189 static void nnpi_triangle_process( nnpi* nn, point* p, int i )
190 {
191  delaunay* d = nn->d;
192  triangle* t = &d->triangles[i];
193  circle * c = &d->circles[i];
194  circle cs[3];
195  int j;
196 
197  assert( circle_contains( c, p ) );
198 
199  if ( nn_rule == SIBSON )
200  {
201  point pp;
202 
203  pp.x = p->x;
204  pp.y = p->y;
205  //
206  // Sibson interpolation by using Watson's algorithm
207  //
208  do
209  {
210  for ( j = 0; j < 3; ++j )
211  {
212  int j1 = ( j + 1 ) % 3;
213  int j2 = ( j + 2 ) % 3;
214  int v1 = t->vids[j1];
215  int v2 = t->vids[j2];
216 
217  if ( !circle_build( &cs[j], &d->points[v1], &d->points[v2], &pp ) )
218  {
219  double scale = triangle_scale_get( d, t );
220 
221  if ( d->points[v1].y == d->points[v2].y )
222  pp.y += EPS_SHIFT * scale;
223  else
224  pp.x += EPS_SHIFT * scale;
225  break;
226  }
227  }
228  } while ( j != 3 );
229 
230  for ( j = 0; j < 3; ++j )
231  {
232  int j1 = ( j + 1 ) % 3;
233  int j2 = ( j + 2 ) % 3;
234  double det = ( ( cs[j1].x - c->x ) * ( cs[j2].y - c->y ) - ( cs[j2].x - c->x ) * ( cs[j1].y - c->y ) );
235 
236  nnpi_add_weight( nn, t->vids[j], det );
237  }
238  }
239  else if ( nn_rule == NON_SIBSONIAN )
240  {
241  double d1 = c->r - hypot( p->x - c->x, p->y - c->y );
242 
243  for ( i = 0; i < 3; ++i )
244  {
245  int vid = t->vids[i];
246  point * pp = &d->points[vid];
247  double d2 = hypot( p->x - pp->x, p->y - pp->y );
248 
249  if ( d2 == 0.0 )
250  nnpi_add_weight( nn, vid, BIGNUMBER );
251  else
252  nnpi_add_weight( nn, vid, d1 / d2 );
253  }
254  }
255  else
256  nn_quit( "unknown rule\n" );
257 }
258 
260 {
261  point* p = nn->p;
262  int n = nn->d->ntriangles;
263  int i;
264 
265  if ( n > N_SEARCH_TURNON )
266  {
267  int* tids;
268 
269  delaunay_circles_find( nn->d, p, &n, &tids );
270  for ( i = 0; i < n; ++i )
271  nnpi_triangle_process( nn, p, tids[i] );
272  }
273  else
274  for ( i = 0; i < n; ++i )
275  if ( circle_contains( &nn->d->circles[i], p ) )
276  nnpi_triangle_process( nn, p, i );
277 }
278 
280 {
281  int n = nn->nvertices;
282  double sum = 0.0;
283  int i;
284 
285  for ( i = 0; i < n; ++i )
286  sum += nn->weights[i];
287 
288  for ( i = 0; i < n; ++i )
289  nn->weights[i] /= sum;
290 }
291 
292 // Finds Natural Neighbours-interpolated value for a point.
293 //
294 // @param nn NN interpolation
295 // @param p Point to be interpolated (p->x, p->y -- input; p->z -- output)
296 //
298 {
299  delaunay* d = nn->d;
300  int i;
301 
302  nnpi_reset( nn );
303  nn->p = p;
306 
307  if ( nn_verbose )
308  {
309  if ( nn_test_vertice == -1 )
310  {
311  if ( nn->n == 0 )
312  fprintf( stderr, "weights:\n" );
313  fprintf( stderr, " %d: {", nn->n );
314  for ( i = 0; i < nn->nvertices; ++i )
315  {
316  fprintf( stderr, "(%d,%.5g)", nn->vertices[i], nn->weights[i] );
317  if ( i < nn->nvertices - 1 )
318  fprintf( stderr, ", " );
319  }
320  fprintf( stderr, "}\n" );
321  }
322  else
323  {
324  double w = 0.0;
325 
326  if ( nn->n == 0 )
327  fprintf( stderr, "weights for vertex %d:\n", nn_test_vertice );
328  for ( i = 0; i < nn->nvertices; ++i )
329  {
330  if ( nn->vertices[i] == nn_test_vertice )
331  {
332  w = nn->weights[i];
333  break;
334  }
335  }
336  fprintf( stderr, "%15.7g %15.7g %15.7g\n", p->x, p->y, w );
337  }
338  }
339 
340  nn->n++;
341 
342  if ( nn->nvertices == 0 )
343  {
344  p->z = NaN;
345  return;
346  }
347 
348  p->z = 0.0;
349  for ( i = 0; i < nn->nvertices; ++i )
350  {
351  double weight = nn->weights[i];
352 
353  if ( weight < nn->wmin )
354  {
355  p->z = NaN;
356  return;
357  }
358  p->z += d->points[nn->vertices[i]].z * weight;
359  }
360 }
361 
362 // Performs Natural Neighbours interpolation for an array of points.
363 //
364 // @param nin Number of input points
365 // @param pin Array of input points [pin]
366 // @param wmin Minimal allowed weight
367 // @param nout Number of output points
368 // @param pout Array of output points [nout]
369 //
370 void nnpi_interpolate_points( int nin, point pin[], double wmin, int nout, point pout[] )
371 {
372  delaunay* d = delaunay_build( nin, pin, 0, NULL, 0, NULL );
373  nnpi * nn = nnpi_create( d );
374  int seed = 0;
375  int i;
376 
377  nn->wmin = wmin;
378 
379  if ( nn_verbose )
380  {
381  fprintf( stderr, "xytoi:\n" );
382  for ( i = 0; i < nout; ++i )
383  {
384  point* p = &pout[i];
385 
386  fprintf( stderr, "(%.7g,%.7g) -> %d\n", p->x, p->y, delaunay_xytoi( d, p, seed ) );
387  }
388  }
389 
390  for ( i = 0; i < nout; ++i )
391  nnpi_interpolate_point( nn, &pout[i] );
392 
393  if ( nn_verbose )
394  {
395  fprintf( stderr, "output:\n" );
396  for ( i = 0; i < nout; ++i )
397  {
398  point* p = &pout[i];
399 
400  fprintf( stderr, " %d:%15.7g %15.7g %15.7g\n", i, p->x, p->y, p->z );
401  }
402  }
403 
404  nnpi_destroy( nn );
405  delaunay_destroy( d );
406 }
407 
408 // Sets minimal allowed weight for Natural Neighbours interpolation.
409 // @param nn Natural Neighbours point interpolator
410 // @param wmin Minimal allowed weight
411 //
412 void nnpi_setwmin( nnpi* nn, double wmin )
413 {
414  nn->wmin = wmin;
415 }
416 
417 // Sets point to interpolate in.
418 // @param nn Natural Neighbours point interpolator
419 // @param p Point to interpolate in
420 //
421 void nnpi_set_point( nnpi* nn, point* p )
422 {
423  nn->p = p;
424 }
425 
426 // Gets number of data points involved in current interpolation.
427 // @return Number of data points involved in current interpolation
428 //
430 {
431  return nn->nvertices;
432 }
433 
434 // Gets indices of data points involved in current interpolation.
435 // @return indices of data points involved in current interpolation
436 //
438 {
439  return nn->vertices;
440 }
441 
442 // Gets weights of data points involved in current interpolation.
443 // @return weights of data points involved in current interpolation
444 //
445 double* nnpi_get_weights( nnpi* nn )
446 {
447  return nn->weights;
448 }
449 
450 //
451 // nnhpi
452 //
453 
454 struct nnhpi
455 {
459  int n; // number of points processed
460 };
461 
462 typedef struct
463 {
464  int nvertices;
465  int * vertices; // vertex indices [nvertices]
466  double* weights; // vertex weights [nvertices]
467 } nn_weights;
468 
469 // Creates Natural Neighbours hashing point interpolator.
470 //
471 // @param d Delaunay triangulation
472 // @param size Hash table size (should be of order of number of output points)
473 // @return Natural Neighbours interpolation
474 //
475 nnhpi* nnhpi_create( delaunay* d, int size )
476 {
477  nnhpi* nn = malloc( sizeof ( nnhpi ) );
478  int i;
479 
480  nn->nnpi = nnpi_create( d );
481 
482  nn->ht_data = ht_create_d2( d->npoints );
483  nn->ht_weights = ht_create_d2( size );
484  nn->n = 0;
485 
486  for ( i = 0; i < d->npoints; ++i )
487  ht_insert( nn->ht_data, &d->points[i], &d->points[i] );
488 
489  return nn;
490 }
491 
492 static void free_nn_weights( void* data )
493 {
494  nn_weights* weights = (nn_weights *) data;
495 
496  free( weights->vertices );
497  free( weights->weights );
498  free( weights );
499 }
500 
501 // Destroys Natural Neighbours hashing point interpolation.
502 //
503 // @param nn Structure to be destroyed
504 //
505 void nnhpi_destroy( nnhpi* nn )
506 {
507  ht_destroy( nn->ht_data );
509  ht_destroy( nn->ht_weights );
510  nnpi_destroy( nn->nnpi );
511 }
512 
513 // Finds Natural Neighbours-interpolated value in a point.
514 //
515 // @param nnhp NN point hashing interpolator
516 // @param p Point to be interpolated (p->x, p->y -- input; p->z -- output)
517 //
518 void nnhpi_interpolate( nnhpi* nnhp, point* p )
519 {
520  nnpi * nnp = nnhp->nnpi;
521  delaunay * d = nnp->d;
522  hashtable * ht_weights = nnhp->ht_weights;
523  nn_weights* weights;
524  int i;
525 
526  if ( ht_find( ht_weights, p ) != NULL )
527  {
528  weights = ht_find( ht_weights, p );
529  if ( nn_verbose )
530  fprintf( stderr, " <hashtable>\n" );
531  }
532  else
533  {
534  nnpi_reset( nnp );
535  nnp->p = p;
536  nnpi_calculate_weights( nnp );
537  nnpi_normalize_weights( nnp );
538 
539  weights = malloc( sizeof ( nn_weights ) );
540  weights->vertices = malloc( sizeof ( int ) * (size_t) ( nnp->nvertices ) );
541  weights->weights = malloc( sizeof ( double ) * (size_t) ( nnp->nvertices ) );
542 
543  weights->nvertices = nnp->nvertices;
544 
545  for ( i = 0; i < nnp->nvertices; ++i )
546  {
547  weights->vertices[i] = nnp->vertices[i];
548  weights->weights[i] = nnp->weights[i];
549  }
550 
551  ht_insert( ht_weights, p, weights );
552 
553  if ( nn_verbose )
554  {
555  if ( nn_test_vertice == -1 )
556  {
557  if ( nnp->n == 0 )
558  fprintf( stderr, "weights:\n" );
559  fprintf( stderr, " %d: {", nnp->n );
560 
561  for ( i = 0; i < nnp->nvertices; ++i )
562  {
563  fprintf( stderr, "(%d,%.5g)", nnp->vertices[i], nnp->weights[i] );
564 
565  if ( i < nnp->nvertices - 1 )
566  fprintf( stderr, ", " );
567  }
568  fprintf( stderr, "}\n" );
569  }
570  else
571  {
572  double w = 0.0;
573 
574  if ( nnp->n == 0 )
575  fprintf( stderr, "weights for vertex %d:\n", nn_test_vertice );
576  for ( i = 0; i < nnp->nvertices; ++i )
577  {
578  if ( nnp->vertices[i] == nn_test_vertice )
579  {
580  w = nnp->weights[i];
581 
582  break;
583  }
584  }
585  fprintf( stderr, "%15.7g %15.7g %15.7g\n", p->x, p->y, w );
586  }
587  }
588 
589  nnp->n++;
590  }
591 
592  nnhp->n++;
593 
594  if ( weights->nvertices == 0 )
595  {
596  p->z = NaN;
597  return;
598  }
599 
600  p->z = 0.0;
601  for ( i = 0; i < weights->nvertices; ++i )
602  {
603  if ( weights->weights[i] < nnp->wmin )
604  {
605  p->z = NaN;
606  return;
607  }
608  p->z += d->points[weights->vertices[i]].z * weights->weights[i];
609  }
610 }
611 
612 // Modifies interpolated data.
613 // Finds point* pd in the underlying Delaunay triangulation such that
614 // pd->x = p->x and pd->y = p->y, and copies p->z to pd->z. Exits with error
615 // if the point is not found.
616 //
617 // @param nnhp Natural Neighbours hashing point interpolator
618 // @param p New data
619 //
620 void nnhpi_modify_data( nnhpi* nnhp, point* p )
621 {
622  point* orig = ht_find( nnhp->ht_data, p );
623 
624  assert( orig != NULL );
625  orig->z = p->z;
626 }
627 
628 // Sets minimal allowed weight for Natural Neighbours interpolation.
629 // @param nn Natural Neighbours point hashing interpolator
630 // @param wmin Minimal allowed weight
631 //
632 void nnhpi_setwmin( nnhpi* nn, double wmin )
633 {
634  nn->nnpi->wmin = wmin;
635 }
636 
637 #if defined ( NNPHI_TEST )
638 
639 #include <sys/time.h>
640 
641 #define NPOINTSIN 10000
642 #define NMIN 10
643 #define NX 101
644 #define NXMIN 1
645 
646 #define SQ( x ) ( ( x ) * ( x ) )
647 
648 static double franke( double x, double y )
649 {
650  x *= 9.0;
651  y *= 9.0;
652  return 0.75 * exp( ( -SQ( x - 2.0 ) - SQ( y - 2.0 ) ) / 4.0 )
653  + 0.75 * exp( -SQ( x - 2.0 ) / 49.0 - ( y - 2.0 ) / 10.0 )
654  + 0.5 * exp( ( -SQ( x - 7.0 ) - SQ( y - 3.0 ) ) / 4.0 )
655  - 0.2 * exp( -SQ( x - 4.0 ) - SQ( y - 7.0 ) );
656 }
657 
658 static void usage()
659 {
660  printf( "Usage: nnhpi_test [-a] [-n <nin> <nxout>] [-v|-V]\n" );
661  printf( "Options:\n" );
662  printf( " -a -- use non-Sibsonian interpolation rule\n" );
663  printf( " -n <nin> <nout>:\n" );
664  printf( " <nin> -- number of input points (default = 10000)\n" );
665  printf( " <nout> -- number of output points per side (default = 64)\n" );
666  printf( " -v -- verbose\n" );
667  printf( " -V -- very verbose\n" );
668 
669  exit( 0 );
670 }
671 
672 int main( int argc, char* argv[] )
673 {
674  int nin = NPOINTSIN;
675  int nx = NX;
676  int nout = 0;
677  point * pin = NULL;
678  delaunay * d = NULL;
679  point * pout = NULL;
680  nnhpi * nn = NULL;
681  int cpi = -1; // control point index
682  struct timeval tv0, tv1;
683  struct timezone tz;
684  int i;
685 
686  i = 1;
687  while ( i < argc )
688  {
689  switch ( argv[i][1] )
690  {
691  case 'a':
692  i++;
694  break;
695  case 'n':
696  i++;
697  if ( i >= argc )
698  nn_quit( "no number of data points found after -n\n" );
699  nin = atoi( argv[i] );
700  i++;
701  if ( i >= argc )
702  nn_quit( "no number of ouput points per side found after -i\n" );
703  nx = atoi( argv[i] );
704  i++;
705  break;
706  case 'v':
707  i++;
708  nn_verbose = 1;
709  break;
710  case 'V':
711  i++;
712  nn_verbose = 2;
713  break;
714  default:
715  usage();
716  break;
717  }
718  }
719 
720  if ( nin < NMIN )
721  nin = NMIN;
722  if ( nx < NXMIN )
723  nx = NXMIN;
724 
725  printf( "\nTest of Natural Neighbours hashing point interpolator:\n\n" );
726  printf( " %d data points\n", nin );
727  printf( " %d output points\n", nx * nx );
728 
729  //
730  // generate data
731  //
732  printf( " generating data:\n" );
733  fflush( stdout );
734  pin = malloc( nin * sizeof ( point ) );
735  for ( i = 0; i < nin; ++i )
736  {
737  point* p = &pin[i];
738 
739  p->x = (double) random() / RAND_MAX;
740  p->y = (double) random() / RAND_MAX;
741  p->z = franke( p->x, p->y );
742  if ( nn_verbose )
743  printf( " (%f, %f, %f)\n", p->x, p->y, p->z );
744  }
745 
746  //
747  // triangulate
748  //
749  printf( " triangulating:\n" );
750  fflush( stdout );
751  d = delaunay_build( nin, pin, 0, NULL, 0, NULL );
752 
753  //
754  // generate output points
755  //
756  points_generate2( -0.1, 1.1, -0.1, 1.1, nx, nx, &nout, &pout );
757  cpi = ( nx / 2 ) * ( nx + 1 );
758 
759  gettimeofday( &tv0, &tz );
760 
761  //
762  // create interpolator
763  //
764  printf( " creating interpolator:\n" );
765  fflush( stdout );
766  nn = nnhpi_create( d, nout );
767 
768  fflush( stdout );
769  gettimeofday( &tv1, &tz );
770  {
771  long dt = 1000000 * ( tv1.tv_sec - tv0.tv_sec ) + tv1.tv_usec - tv0.tv_usec;
772 
773  printf( " interpolator creation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout );
774  }
775 
776  //
777  // interpolate
778  //
779  printf( " interpolating:\n" );
780  fflush( stdout );
781  gettimeofday( &tv1, &tz );
782  for ( i = 0; i < nout; ++i )
783  {
784  point* p = &pout[i];
785 
786  nnhpi_interpolate( nn, p );
787  if ( nn_verbose )
788  printf( " (%f, %f, %f)\n", p->x, p->y, p->z );
789  }
790 
791  fflush( stdout );
792  gettimeofday( &tv0, &tz );
793  {
794  long dt = 1000000.0 * ( tv0.tv_sec - tv1.tv_sec ) + tv0.tv_usec - tv1.tv_usec;
795 
796  printf( " interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout );
797  }
798 
799  if ( !nn_verbose )
800  printf( " control point: (%f, %f, %f) (expected z = %f)\n", pout[cpi].x, pout[cpi].y, pout[cpi].z, franke( pout[cpi].x, pout[cpi].y ) );
801 
802  printf( " interpolating one more time:\n" );
803  fflush( stdout );
804  gettimeofday( &tv0, &tz );
805  for ( i = 0; i < nout; ++i )
806  {
807  point* p = &pout[i];
808 
809  nnhpi_interpolate( nn, p );
810  if ( nn_verbose )
811  printf( " (%f, %f, %f)\n", p->x, p->y, p->z );
812  }
813 
814  fflush( stdout );
815  gettimeofday( &tv1, &tz );
816  {
817  long dt = 1000000.0 * ( tv1.tv_sec - tv0.tv_sec ) + tv1.tv_usec - tv0.tv_usec;
818 
819  printf( " interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout );
820  }
821 
822  if ( !nn_verbose )
823  printf( " control point: (%f, %f, %f) (expected z = %f)\n", pout[cpi].x, pout[cpi].y, pout[cpi].z, franke( pout[cpi].x, pout[cpi].y ) );
824 
825  printf( " entering new data:\n" );
826  fflush( stdout );
827  for ( i = 0; i < nin; ++i )
828  {
829  point* p = &pin[i];
830 
831  p->z = p->x * p->x - p->y * p->y;
832  nnhpi_modify_data( nn, p );
833  if ( nn_verbose )
834  printf( " (%f, %f, %f)\n", p->x, p->y, p->z );
835  }
836 
837  printf( " interpolating:\n" );
838  fflush( stdout );
839  gettimeofday( &tv1, &tz );
840  for ( i = 0; i < nout; ++i )
841  {
842  point* p = &pout[i];
843 
844  nnhpi_interpolate( nn, p );
845  if ( nn_verbose )
846  printf( " (%f, %f, %f)\n", p->x, p->y, p->z );
847  }
848 
849  fflush( stdout );
850  gettimeofday( &tv0, &tz );
851  {
852  long dt = 1000000.0 * ( tv0.tv_sec - tv1.tv_sec ) + tv0.tv_usec - tv1.tv_usec;
853 
854  printf( " interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout );
855  }
856 
857  if ( !nn_verbose )
858  printf( " control point: (%f, %f, %f) (expected z = %f)\n", pout[cpi].x, pout[cpi].y, pout[cpi].z, pout[cpi].x * pout[cpi].x - pout[cpi].y * pout[cpi].y );
859 
860  printf( " restoring data:\n" );
861  fflush( stdout );
862  for ( i = 0; i < nin; ++i )
863  {
864  point* p = &pin[i];
865 
866  p->z = franke( p->x, p->y );
867  nnhpi_modify_data( nn, p );
868  if ( nn_verbose )
869  printf( " (%f, %f, %f)\n", p->x, p->y, p->z );
870  }
871 
872  printf( " interpolating:\n" );
873  fflush( stdout );
874  gettimeofday( &tv0, &tz );
875  for ( i = 0; i < nout; ++i )
876  {
877  point* p = &pout[i];
878 
879  nnhpi_interpolate( nn, p );
880  if ( nn_verbose )
881  printf( " (%f, %f, %f)\n", p->x, p->y, p->z );
882  }
883 
884  fflush( stdout );
885  gettimeofday( &tv1, &tz );
886  {
887  long dt = 1000000.0 * ( tv1.tv_sec - tv0.tv_sec ) + tv1.tv_usec - tv0.tv_usec;
888 
889  printf( " interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout );
890  }
891 
892  if ( !nn_verbose )
893  printf( " control point: (%f, %f, %f) (expected z = %f)\n", pout[cpi].x, pout[cpi].y, pout[cpi].z, franke( pout[cpi].x, pout[cpi].y ) );
894 
895  printf( " hashtable stats:\n" );
896  fflush( stdout );
897  {
898  hashtable* ht = nn->ht_data;
899 
900  printf( " input points: %d entries, %d table elements, %d filled elements\n", ht->n, ht->size, ht->nhash );
901  ht = nn->ht_weights;
902  printf( " weights: %d entries, %d table elements, %d filled elements\n", ht->n, ht->size, ht->nhash );
903  }
904  printf( "\n" );
905 
906  nnhpi_destroy( nn );
907  free( pout );
908  delaunay_destroy( d );
909  free( pin );
910 
911  return 0;
912 }
913 
914 #endif