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 {
55 delaunay* d;
56 point * p;
57 double wmin;
58 //
59 // work variables
60 //
61 int nvertices;
62 int nallocated;
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 //
nnpi_create(delaunay * d)95 nnpi* nnpi_create( delaunay* d )
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 //
nnpi_destroy(nnpi * nn)115 void nnpi_destroy( nnpi* nn )
116 {
117 free( nn->weights );
118 free( nn->vertices );
119 free( nn );
120 }
121
nnpi_reset(nnpi * nn)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
nnpi_add_weight(nnpi * nn,int vertex,double w)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
triangle_scale_get(delaunay * d,triangle * t)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 //
nnpi_triangle_process(nnpi * nn,point * p,int i)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
nnpi_calculate_weights(nnpi * nn)259 void nnpi_calculate_weights( nnpi* nn )
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
nnpi_normalize_weights(nnpi * nn)279 void nnpi_normalize_weights( nnpi* nn )
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 //
nnpi_interpolate_point(nnpi * nn,point * p)297 void nnpi_interpolate_point( nnpi* nn, point* p )
298 {
299 delaunay* d = nn->d;
300 int i;
301
302 nnpi_reset( nn );
303 nn->p = p;
304 nnpi_calculate_weights( nn );
305 nnpi_normalize_weights( nn );
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 //
nnpi_interpolate_points(int nin,point pin[],double wmin,int nout,point pout[])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 //
nnpi_setwmin(nnpi * nn,double wmin)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 //
nnpi_set_point(nnpi * nn,point * p)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 //
nnpi_get_nvertices(nnpi * nn)429 int nnpi_get_nvertices( nnpi* nn )
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 //
nnpi_get_vertices(nnpi * nn)437 int* nnpi_get_vertices( nnpi* nn )
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 //
nnpi_get_weights(nnpi * nn)445 double* nnpi_get_weights( nnpi* nn )
446 {
447 return nn->weights;
448 }
449
450 //
451 // nnhpi
452 //
453
454 struct nnhpi
455 {
456 nnpi * nnpi;
457 hashtable* ht_data;
458 hashtable* ht_weights;
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 //
nnhpi_create(delaunay * d,int size)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
free_nn_weights(void * data)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 //
nnhpi_destroy(nnhpi * nn)505 void nnhpi_destroy( nnhpi* nn )
506 {
507 ht_destroy( nn->ht_data );
508 ht_process( nn->ht_weights, free_nn_weights );
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 //
nnhpi_interpolate(nnhpi * nnhp,point * p)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 //
nnhpi_modify_data(nnhpi * nnhp,point * p)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 //
nnhpi_setwmin(nnhpi * nn,double wmin)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
franke(double x,double y)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
usage()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
main(int argc,char * argv[])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++;
693 nn_rule = NON_SIBSONIAN;
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
915