72 void nn_quit(
const char* format, ... );
83 #define EPS_SHIFT 1.0e-9
84 #define N_SEARCH_TURNON 20
85 #define BIGNUMBER 1.0e+100
87 #define min( x, y ) ( ( x ) < ( y ) ? ( x ) : ( y ) )
88 #define max( x, y ) ( ( x ) > ( y ) ? ( x ) : ( y ) )
97 nnpi* nn = malloc(
sizeof (
nnpi ) );
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 );
181 return xmax - xmin + ymax - ymin;
210 for ( j = 0; j < 3; ++j )
212 int j1 = ( j + 1 ) % 3;
213 int j2 = ( j + 2 ) % 3;
214 int v1 = t->
vids[j1];
215 int v2 = t->
vids[j2];
230 for ( j = 0; j < 3; ++j )
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 ) );
241 double d1 = c->
r - hypot( p->
x - c->
x, p->
y - c->
y );
243 for ( i = 0; i < 3; ++i )
245 int vid = t->
vids[i];
247 double d2 = hypot( p->
x - pp->
x, p->
y - pp->
y );
270 for ( i = 0; i < n; ++i )
274 for ( i = 0; i < n; ++i )
285 for ( i = 0; i < n; ++i )
288 for ( i = 0; i < n; ++i )
312 fprintf( stderr,
"weights:\n" );
313 fprintf( stderr,
" %d: {", nn->
n );
317 if ( i < nn->nvertices - 1 )
318 fprintf( stderr,
", " );
320 fprintf( stderr,
"}\n" );
336 fprintf( stderr,
"%15.7g %15.7g %15.7g\n", p->
x, p->
y, w );
351 double weight = nn->
weights[i];
353 if ( weight < nn->wmin )
381 fprintf( stderr,
"xytoi:\n" );
382 for ( i = 0; i < nout; ++i )
386 fprintf( stderr,
"(%.7g,%.7g) -> %d\n", p->
x, p->
y,
delaunay_xytoi( d, p, seed ) );
390 for ( i = 0; i < nout; ++i )
395 fprintf( stderr,
"output:\n" );
396 for ( i = 0; i < nout; ++i )
400 fprintf( stderr,
" %d:%15.7g %15.7g %15.7g\n", i, p->
x, p->
y, p->
z );
486 for ( i = 0; i < d->
npoints; ++i )
526 if (
ht_find( ht_weights, p ) != NULL )
528 weights =
ht_find( ht_weights, p );
530 fprintf( stderr,
" <hashtable>\n" );
541 weights->
weights = malloc(
sizeof (
double ) * (
size_t) ( nnp->
nvertices ) );
558 fprintf( stderr,
"weights:\n" );
559 fprintf( stderr,
" %d: {", nnp->
n );
565 if ( i < nnp->nvertices - 1 )
566 fprintf( stderr,
", " );
568 fprintf( stderr,
"}\n" );
585 fprintf( stderr,
"%15.7g %15.7g %15.7g\n", p->
x, p->
y, w );
601 for ( i = 0; i < weights->
nvertices; ++i )
624 assert( orig != NULL );
637 #if defined ( NNPHI_TEST )
639 #include <sys/time.h>
641 #define NPOINTSIN 10000
646 #define SQ( x ) ( ( x ) * ( x ) )
648 static double franke(
double x,
double y )
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 ) );
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" );
682 struct timeval tv0, tv1;
689 switch ( argv[i][1] )
698 nn_quit(
"no number of data points found after -n\n" );
699 nin = atoi( argv[i] );
702 nn_quit(
"no number of ouput points per side found after -i\n" );
703 nx = atoi( argv[i] );
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 );
732 printf(
" generating data:\n" );
734 pin = malloc( nin *
sizeof (
point ) );
735 for ( i = 0; i < nin; ++i )
739 p->
x = (double) random() / RAND_MAX;
740 p->
y = (double) random() / RAND_MAX;
741 p->
z = franke( p->
x, p->
y );
743 printf(
" (%f, %f, %f)\n", p->
x, p->
y, p->
z );
749 printf(
" triangulating:\n" );
757 cpi = ( nx / 2 ) * ( nx + 1 );
759 gettimeofday( &tv0, &tz );
764 printf(
" creating interpolator:\n" );
769 gettimeofday( &tv1, &tz );
771 long dt = 1000000 * ( tv1.tv_sec - tv0.tv_sec ) + tv1.tv_usec - tv0.tv_usec;
773 printf(
" interpolator creation time = %ld us (%.2f us / point)\n", dt, (
double) dt / nout );
779 printf(
" interpolating:\n" );
781 gettimeofday( &tv1, &tz );
782 for ( i = 0; i < nout; ++i )
788 printf(
" (%f, %f, %f)\n", p->
x, p->
y, p->
z );
792 gettimeofday( &tv0, &tz );
794 long dt = 1000000.0 * ( tv0.tv_sec - tv1.tv_sec ) + tv0.tv_usec - tv1.tv_usec;
796 printf(
" interpolation time = %ld us (%.2f us / point)\n", dt, (
double) dt / nout );
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 ) );
802 printf(
" interpolating one more time:\n" );
804 gettimeofday( &tv0, &tz );
805 for ( i = 0; i < nout; ++i )
811 printf(
" (%f, %f, %f)\n", p->
x, p->
y, p->
z );
815 gettimeofday( &tv1, &tz );
817 long dt = 1000000.0 * ( tv1.tv_sec - tv0.tv_sec ) + tv1.tv_usec - tv0.tv_usec;
819 printf(
" interpolation time = %ld us (%.2f us / point)\n", dt, (
double) dt / nout );
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 ) );
825 printf(
" entering new data:\n" );
827 for ( i = 0; i < nin; ++i )
831 p->
z = p->
x * p->
x - p->
y * p->
y;
834 printf(
" (%f, %f, %f)\n", p->
x, p->
y, p->
z );
837 printf(
" interpolating:\n" );
839 gettimeofday( &tv1, &tz );
840 for ( i = 0; i < nout; ++i )
846 printf(
" (%f, %f, %f)\n", p->
x, p->
y, p->
z );
850 gettimeofday( &tv0, &tz );
852 long dt = 1000000.0 * ( tv0.tv_sec - tv1.tv_sec ) + tv0.tv_usec - tv1.tv_usec;
854 printf(
" interpolation time = %ld us (%.2f us / point)\n", dt, (
double) dt / nout );
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 );
860 printf(
" restoring data:\n" );
862 for ( i = 0; i < nin; ++i )
866 p->
z = franke( p->
x, p->
y );
869 printf(
" (%f, %f, %f)\n", p->
x, p->
y, p->
z );
872 printf(
" interpolating:\n" );
874 gettimeofday( &tv0, &tz );
875 for ( i = 0; i < nout; ++i )
881 printf(
" (%f, %f, %f)\n", p->
x, p->
y, p->
z );
885 gettimeofday( &tv1, &tz );
887 long dt = 1000000.0 * ( tv1.tv_sec - tv0.tv_sec ) + tv1.tv_usec - tv0.tv_usec;
889 printf(
" interpolation time = %ld us (%.2f us / point)\n", dt, (
double) dt / nout );
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 ) );
895 printf(
" hashtable stats:\n" );
900 printf(
" input points: %d entries, %d table elements, %d filled elements\n", ht->
n, ht->
size, ht->
nhash );
902 printf(
" weights: %d entries, %d table elements, %d filled elements\n", ht->
n, ht->
size, ht->
nhash );
void * ht_find(hashtable *table, void *key)
static void nnpi_add_weight(nnpi *nn, int vertex, double w)
void nnhpi_modify_data(nnhpi *nnhp, point *p)
static PLCHAR_VECTOR usage
int circle_contains(circle *c, point *p)
double * nnpi_get_weights(nnpi *nn)
void delaunay_circles_find(delaunay *d, point *p, int *n, int **out)
void nnpi_interpolate_points(int nin, point pin[], double wmin, int nout, point pout[])
void ht_destroy(hashtable *table)
static void free_nn_weights(void *data)
void nnhpi_interpolate(nnhpi *nnhp, point *p)
void nnhpi_destroy(nnhpi *nn)
void nnpi_interpolate_point(nnpi *nn, point *p)
void nnhpi_setwmin(nnhpi *nn, double wmin)
void delaunay_destroy(delaunay *d)
int nnpi_get_nvertices(nnpi *nn)
nnpi * nnpi_create(delaunay *d)
nnhpi * nnhpi_create(delaunay *d, int size)
static void nnpi_triangle_process(nnpi *nn, point *p, int i)
delaunay * delaunay_build(int np, point points[], int ns, int segments[], int nh, double holes[])
int circle_build(circle *c, point *p0, point *p1, point *p2)
void ht_process(hashtable *table, void(*func)(void *))
void nnpi_destroy(nnpi *nn)
void nnpi_set_point(nnpi *nn, point *p)
void nnpi_reset(nnpi *nn)
void nn_quit(const char *format,...)
void nnpi_normalize_weights(nnpi *nn)
void nnpi_setwmin(nnpi *nn, double wmin)
void * ht_insert(hashtable *table, void *key, void *data)
void nnpi_calculate_weights(nnpi *nn)
int * nnpi_get_vertices(nnpi *nn)
int delaunay_xytoi(delaunay *d, point *p, int seed)
hashtable * ht_create_d2(int size)
void points_generate2(double xmin, double xmax, double ymin, double ymax, int nx, int ny, int *nout, point **pout)
static double triangle_scale_get(delaunay *d, triangle *t)