PLplot 5.15.0
Loading...
Searching...
No Matches
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
53struct nnpi
54{
57 double wmin;
58 //
59 // work variables
60 //
63 int * vertices; // vertex indices
64 double* weights;
65 int n; // number of points processed
66};
67
68int circle_build( circle* c, point* p0, point* p1, point* p2 );
69int circle_contains( circle* c, point* p );
70void delaunay_circles_find( delaunay* d, point* p, int* n, int** out );
71int delaunay_xytoi( delaunay* d, point* p, int seed );
72void nn_quit( const char* format, ... );
73void nnpi_reset( nnpi* nn );
76void nnpi_set_point( nnpi* nn, point* p );
77int nnpi_get_nvertices( nnpi* nn );
78int* nnpi_get_vertices( nnpi* nn );
79double* 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//
116{
117 free( nn->weights );
118 free( nn->vertices );
119 free( nn );
120}
121
122void 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
129static 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
168static 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//
189static 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//
370void 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//
412void 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//
421void 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//
445double* nnpi_get_weights( nnpi* nn )
446{
447 return nn->weights;
448}
449
450//
451// nnhpi
452//
453
454struct nnhpi
455{
459 int n; // number of points processed
460};
461
462typedef 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//
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
492static 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//
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//
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;
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//
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//
632void 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
648static 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
658static 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
672int 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
int main()
Definition cdexpert.c:35
#define NaN
Definition csa/nan.h:62
int delaunay_xytoi(delaunay *d, point *p, int id)
Definition delaunay.c:631
delaunay * delaunay_build(int np, point points[], int ns, int segments[], int nh, double holes[])
Definition delaunay.c:265
void delaunay_destroy(delaunay *d)
Definition delaunay.c:578
hashtable * ht_create_d2(int size)
Definition hash.c:401
void * ht_find(hashtable *table, void *key)
Definition hash.c:210
void ht_destroy(hashtable *table)
Definition hash.c:102
void ht_process(hashtable *table, void(*func)(void *))
Definition hash.c:290
void * ht_insert(hashtable *table, void *key, void *data)
Definition hash.c:135
void points_generate2(double xmin, double xmax, double ymin, double ymax, int nx, int ny, int *nout, point **pout)
Definition nncommon.c:334
@ NON_SIBSONIAN
Definition nn.h:23
@ SIBSON
Definition nn.h:23
int nn_test_vertice
Definition nncommon.c:44
int nn_verbose
Definition nncommon.c:43
NN_RULE nn_rule
Definition nncommon.c:45
static void nnpi_add_weight(nnpi *nn, int vertex, double w)
Definition nnpi.c:129
int circle_build(circle *c, point *p0, point *p1, point *p2)
Definition nncommon.c:68
static void free_nn_weights(void *data)
Definition nnpi.c:492
#define NSTART
Definition nnpi.c:81
void nnpi_destroy(nnpi *nn)
Definition nnpi.c:115
void nnpi_set_point(nnpi *nn, point *p)
Definition nnpi.c:421
void nnpi_normalize_weights(nnpi *nn)
Definition nnpi.c:279
int * nnpi_get_vertices(nnpi *nn)
Definition nnpi.c:437
void nnpi_interpolate_point(nnpi *nn, point *p)
Definition nnpi.c:297
void nnhpi_interpolate(nnhpi *nnhp, point *p)
Definition nnpi.c:518
void nnpi_setwmin(nnpi *nn, double wmin)
Definition nnpi.c:412
void nnhpi_modify_data(nnhpi *nnhp, point *p)
Definition nnpi.c:620
int circle_contains(circle *c, point *p)
Definition nncommon.c:98
#define BIGNUMBER
Definition nnpi.c:85
nnhpi * nnhpi_create(delaunay *d, int size)
Definition nnpi.c:475
int nnpi_get_nvertices(nnpi *nn)
Definition nnpi.c:429
void nn_quit(const char *format,...)
Definition nncommon.c:53
static void nnpi_triangle_process(nnpi *nn, point *p, int i)
Definition nnpi.c:189
void delaunay_circles_find(delaunay *d, point *p, int *n, int **out)
Definition delaunay.c:681
#define NINC
Definition nnpi.c:82
static double triangle_scale_get(delaunay *d, triangle *t)
Definition nnpi.c:168
#define min(x, y)
Definition nnpi.c:87
#define N_SEARCH_TURNON
Definition nnpi.c:84
#define max(x, y)
Definition nnpi.c:88
void nnpi_interpolate_points(int nin, point pin[], double wmin, int nout, point pout[])
Definition nnpi.c:370
int delaunay_xytoi(delaunay *d, point *p, int seed)
Definition delaunay.c:631
double * nnpi_get_weights(nnpi *nn)
Definition nnpi.c:445
void nnhpi_setwmin(nnhpi *nn, double wmin)
Definition nnpi.c:632
void nnpi_reset(nnpi *nn)
Definition nnpi.c:122
void nnpi_calculate_weights(nnpi *nn)
Definition nnpi.c:259
#define EPS_SHIFT
Definition nnpi.c:83
nnpi * nnpi_create(delaunay *d)
Definition nnpi.c:95
void nnhpi_destroy(nnhpi *nn)
Definition nnpi.c:505
static PLCHAR_VECTOR usage
Definition plargs.c:179
#define NX
static int argc
Definition qt.cpp:48
static char ** argv
Definition qt.cpp:49
double x
Definition delaunay.h:35
double r
Definition delaunay.h:37
double y
Definition delaunay.h:36
point * points
Definition delaunay.h:48
int npoints
Definition delaunay.h:47
int ntriangles
Definition delaunay.h:54
circle * circles
Definition delaunay.h:56
triangle * triangles
Definition delaunay.h:55
int * flags
Definition delaunay.h:73
int n
Definition hash.c:41
int nhash
Definition hash.c:43
int size
Definition hash.c:40
int nvertices
Definition nnai.c:34
int * vertices
Definition nnai.c:35
double * weights
Definition nnai.c:36
Definition nnpi.c:455
hashtable * ht_data
Definition nnpi.c:457
int n
Definition nnpi.c:459
nnpi * nnpi
Definition nnpi.c:456
hashtable * ht_weights
Definition nnpi.c:458
Definition nnpi.c:54
point * p
Definition nnpi.c:56
delaunay * d
Definition nnpi.c:55
int nvertices
Definition nnpi.c:61
double * weights
Definition nnpi.c:64
int * vertices
Definition nnpi.c:63
int n
Definition nnpi.c:65
int nallocated
Definition nnpi.c:62
double wmin
Definition nnpi.c:57
Definition csa.h:30
double y
Definition csa.h:32
double x
Definition csa.h:31
double z
Definition csa.h:33
Definition csa.c:56
int vids[3]
Definition delaunay.h:25