PLplot
5.10.0
|
00001 //-------------------------------------------------------------------------- 00002 // 00003 // File: nnpi.c 00004 // 00005 // Created: 15/11/2002 00006 // 00007 // Author: Pavel Sakov 00008 // CSIRO Marine Research 00009 // 00010 // Purpose: Code for: 00011 // -- Natural Neighbours Point Interpolator 00012 // -- Natural Neighbours Point Hashing Interpolator 00013 // 00014 // Description: `nnpi' -- "Natural Neighbours Point 00015 // Interpolator" -- is a structure for conducting Natural 00016 // Neighbours interpolation on a given data on a 00017 // "point-to-point" basis. Because it involves weight 00018 // calculation for each next output point, it is not 00019 // particularly suitable for consequitive interpolations on 00020 // the same set of observation points -- use 00021 // `nnhpi' or `nnai' 00022 // in these cases. 00023 // 00024 // `nnhpi' is a structure for 00025 // conducting consequitive Natural Neighbours interpolations 00026 // on a given spatial data set in a random sequence of points 00027 // from a set of finite size, taking advantage of repeated 00028 // interpolations in the same point. It allows to modify Z 00029 // coordinate of data in between interpolations. 00030 // 00031 // 00032 // Revisions: 01/04/2003 PS: modified nnpi_triangle_process(): for 00033 // Sibson interpolation, if circle_build fails(), now a 00034 // local copy of a point is moved slightly rather than the 00035 // data point itself. The later approach have found leading 00036 // to inconsistencies of the new point position with the 00037 // earlier built triangulation. 00038 // 00039 //-------------------------------------------------------------------------- 00040 00041 #include <stdlib.h> 00042 #include <stdio.h> 00043 #include <limits.h> 00044 #include <float.h> 00045 #include <string.h> 00046 #include <assert.h> 00047 #include <math.h> 00048 #include "nn.h" 00049 #include "delaunay.h" 00050 #include "nan.h" 00051 #include "hash.h" 00052 00053 struct nnpi 00054 { 00055 delaunay* d; 00056 point * p; 00057 double wmin; 00058 // 00059 // work variables 00060 // 00061 int nvertices; 00062 int nallocated; 00063 int * vertices; // vertex indices 00064 double* weights; 00065 int n; // number of points processed 00066 }; 00067 00068 int circle_build( circle* c, point* p0, point* p1, point* p2 ); 00069 int circle_contains( circle* c, point* p ); 00070 void delaunay_circles_find( delaunay* d, point* p, int* n, int** out ); 00071 int delaunay_xytoi( delaunay* d, point* p, int seed ); 00072 void nn_quit( const char* format, ... ); 00073 void nnpi_reset( nnpi* nn ); 00074 void nnpi_calculate_weights( nnpi* nn ); 00075 void nnpi_normalize_weights( nnpi* nn ); 00076 void nnpi_set_point( nnpi* nn, point* p ); 00077 int nnpi_get_nvertices( nnpi* nn ); 00078 int* nnpi_get_vertices( nnpi* nn ); 00079 double* nnpi_get_weights( nnpi* nn ); 00080 00081 #define NSTART 10 00082 #define NINC 10 00083 #define EPS_SHIFT 1.0e-9 00084 #define N_SEARCH_TURNON 20 00085 #define BIGNUMBER 1.0e+100 00086 00087 #define min( x, y ) ( ( x ) < ( y ) ? ( x ) : ( y ) ) 00088 #define max( x, y ) ( ( x ) > ( y ) ? ( x ) : ( y ) ) 00089 00090 // Creates Natural Neighbours point interpolator. 00091 // 00092 // @param d Delaunay triangulation 00093 // @return Natural Neighbours interpolation 00094 // 00095 nnpi* nnpi_create( delaunay* d ) 00096 { 00097 nnpi* nn = malloc( sizeof ( nnpi ) ); 00098 00099 nn->d = d; 00100 nn->wmin = -DBL_MAX; 00101 nn->vertices = calloc( NSTART, sizeof ( int ) ); 00102 nn->weights = calloc( NSTART, sizeof ( double ) ); 00103 nn->nvertices = 0; 00104 nn->nallocated = NSTART; 00105 nn->p = NULL; 00106 nn->n = 0; 00107 00108 return nn; 00109 } 00110 00111 // Destroys Natural Neighbours point interpolator. 00112 // 00113 // @param nn Structure to be destroyed 00114 // 00115 void nnpi_destroy( nnpi* nn ) 00116 { 00117 free( nn->weights ); 00118 free( nn->vertices ); 00119 free( nn ); 00120 } 00121 00122 void nnpi_reset( nnpi* nn ) 00123 { 00124 nn->nvertices = 0; 00125 nn->p = NULL; 00126 memset( nn->d->flags, 0, (size_t) ( nn->d->ntriangles ) * sizeof ( int ) ); 00127 } 00128 00129 static void nnpi_add_weight( nnpi* nn, int vertex, double w ) 00130 { 00131 int i; 00132 00133 // 00134 // find whether the vertex is already in the list 00135 // 00136 for ( i = 0; i < nn->nvertices; ++i ) 00137 if ( nn->vertices[i] == vertex ) 00138 break; 00139 00140 if ( i == nn->nvertices ) // not in the list 00141 { // 00142 // get more memory if necessary 00143 // 00144 if ( nn->nvertices == nn->nallocated ) 00145 { 00146 nn->vertices = realloc( nn->vertices, (size_t) ( nn->nallocated + NINC ) * sizeof ( int ) ); 00147 nn->weights = realloc( nn->weights, (size_t) ( nn->nallocated + NINC ) * sizeof ( double ) ); 00148 nn->nallocated += NINC; 00149 } 00150 00151 // 00152 // add the vertex to the list 00153 // 00154 nn->vertices[i] = vertex; 00155 nn->weights[i] = w; 00156 nn->nvertices++; 00157 } 00158 else // in the list 00159 00160 { 00161 if ( nn_rule == SIBSON ) 00162 nn->weights[i] += w; 00163 else if ( w > nn->weights[i] ) 00164 nn->weights[i] = w; 00165 } 00166 } 00167 00168 static double triangle_scale_get( delaunay* d, triangle* t ) 00169 { 00170 double x0 = d->points[t->vids[0]].x; 00171 double x1 = d->points[t->vids[1]].x; 00172 double x2 = d->points[t->vids[2]].x; 00173 double y0 = d->points[t->vids[0]].y; 00174 double y1 = d->points[t->vids[1]].y; 00175 double y2 = d->points[t->vids[2]].y; 00176 double xmin = min( min( x0, x1 ), x2 ); 00177 double xmax = max( max( x0, x1 ), x2 ); 00178 double ymin = min( min( y0, y1 ), y2 ); 00179 double ymax = max( max( y0, y1 ), y2 ); 00180 00181 return xmax - xmin + ymax - ymin; 00182 } 00183 00184 // This is a central procedure for the Natural Neighbours interpolation. It 00185 // uses the Watson's algorithm for the required areas calculation and implies 00186 // that the vertices of the delaunay triangulation are listed in uniform 00187 // (clockwise or counterclockwise) order. 00188 // 00189 static void nnpi_triangle_process( nnpi* nn, point* p, int i ) 00190 { 00191 delaunay* d = nn->d; 00192 triangle* t = &d->triangles[i]; 00193 circle * c = &d->circles[i]; 00194 circle cs[3]; 00195 int j; 00196 00197 assert( circle_contains( c, p ) ); 00198 00199 if ( nn_rule == SIBSON ) 00200 { 00201 point pp; 00202 00203 pp.x = p->x; 00204 pp.y = p->y; 00205 // 00206 // Sibson interpolation by using Watson's algorithm 00207 // 00208 do 00209 { 00210 for ( j = 0; j < 3; ++j ) 00211 { 00212 int j1 = ( j + 1 ) % 3; 00213 int j2 = ( j + 2 ) % 3; 00214 int v1 = t->vids[j1]; 00215 int v2 = t->vids[j2]; 00216 00217 if ( !circle_build( &cs[j], &d->points[v1], &d->points[v2], &pp ) ) 00218 { 00219 double scale = triangle_scale_get( d, t ); 00220 00221 if ( d->points[v1].y == d->points[v2].y ) 00222 pp.y += EPS_SHIFT * scale; 00223 else 00224 pp.x += EPS_SHIFT * scale; 00225 break; 00226 } 00227 } 00228 } while ( j != 3 ); 00229 00230 for ( j = 0; j < 3; ++j ) 00231 { 00232 int j1 = ( j + 1 ) % 3; 00233 int j2 = ( j + 2 ) % 3; 00234 double det = ( ( cs[j1].x - c->x ) * ( cs[j2].y - c->y ) - ( cs[j2].x - c->x ) * ( cs[j1].y - c->y ) ); 00235 00236 nnpi_add_weight( nn, t->vids[j], det ); 00237 } 00238 } 00239 else if ( nn_rule == NON_SIBSONIAN ) 00240 { 00241 double d1 = c->r - hypot( p->x - c->x, p->y - c->y ); 00242 00243 for ( i = 0; i < 3; ++i ) 00244 { 00245 int vid = t->vids[i]; 00246 point * pp = &d->points[vid]; 00247 double d2 = hypot( p->x - pp->x, p->y - pp->y ); 00248 00249 if ( d2 == 0.0 ) 00250 nnpi_add_weight( nn, vid, BIGNUMBER ); 00251 else 00252 nnpi_add_weight( nn, vid, d1 / d2 ); 00253 } 00254 } 00255 else 00256 nn_quit( "unknown rule\n" ); 00257 } 00258 00259 void nnpi_calculate_weights( nnpi* nn ) 00260 { 00261 point* p = nn->p; 00262 int n = nn->d->ntriangles; 00263 int i; 00264 00265 if ( n > N_SEARCH_TURNON ) 00266 { 00267 int* tids; 00268 00269 delaunay_circles_find( nn->d, p, &n, &tids ); 00270 for ( i = 0; i < n; ++i ) 00271 nnpi_triangle_process( nn, p, tids[i] ); 00272 } 00273 else 00274 for ( i = 0; i < n; ++i ) 00275 if ( circle_contains( &nn->d->circles[i], p ) ) 00276 nnpi_triangle_process( nn, p, i ); 00277 } 00278 00279 void nnpi_normalize_weights( nnpi* nn ) 00280 { 00281 int n = nn->nvertices; 00282 double sum = 0.0; 00283 int i; 00284 00285 for ( i = 0; i < n; ++i ) 00286 sum += nn->weights[i]; 00287 00288 for ( i = 0; i < n; ++i ) 00289 nn->weights[i] /= sum; 00290 } 00291 00292 // Finds Natural Neighbours-interpolated value for a point. 00293 // 00294 // @param nn NN interpolation 00295 // @param p Point to be interpolated (p->x, p->y -- input; p->z -- output) 00296 // 00297 void nnpi_interpolate_point( nnpi* nn, point* p ) 00298 { 00299 delaunay* d = nn->d; 00300 int i; 00301 00302 nnpi_reset( nn ); 00303 nn->p = p; 00304 nnpi_calculate_weights( nn ); 00305 nnpi_normalize_weights( nn ); 00306 00307 if ( nn_verbose ) 00308 { 00309 if ( nn_test_vertice == -1 ) 00310 { 00311 if ( nn->n == 0 ) 00312 fprintf( stderr, "weights:\n" ); 00313 fprintf( stderr, " %d: {", nn->n ); 00314 for ( i = 0; i < nn->nvertices; ++i ) 00315 { 00316 fprintf( stderr, "(%d,%.5g)", nn->vertices[i], nn->weights[i] ); 00317 if ( i < nn->nvertices - 1 ) 00318 fprintf( stderr, ", " ); 00319 } 00320 fprintf( stderr, "}\n" ); 00321 } 00322 else 00323 { 00324 double w = 0.0; 00325 00326 if ( nn->n == 0 ) 00327 fprintf( stderr, "weights for vertex %d:\n", nn_test_vertice ); 00328 for ( i = 0; i < nn->nvertices; ++i ) 00329 { 00330 if ( nn->vertices[i] == nn_test_vertice ) 00331 { 00332 w = nn->weights[i]; 00333 break; 00334 } 00335 } 00336 fprintf( stderr, "%15.7g %15.7g %15.7g\n", p->x, p->y, w ); 00337 } 00338 } 00339 00340 nn->n++; 00341 00342 if ( nn->nvertices == 0 ) 00343 { 00344 p->z = NaN; 00345 return; 00346 } 00347 00348 p->z = 0.0; 00349 for ( i = 0; i < nn->nvertices; ++i ) 00350 { 00351 double weight = nn->weights[i]; 00352 00353 if ( weight < nn->wmin ) 00354 { 00355 p->z = NaN; 00356 return; 00357 } 00358 p->z += d->points[nn->vertices[i]].z * weight; 00359 } 00360 } 00361 00362 // Performs Natural Neighbours interpolation for an array of points. 00363 // 00364 // @param nin Number of input points 00365 // @param pin Array of input points [pin] 00366 // @param wmin Minimal allowed weight 00367 // @param nout Number of output points 00368 // @param pout Array of output points [nout] 00369 // 00370 void nnpi_interpolate_points( int nin, point pin[], double wmin, int nout, point pout[] ) 00371 { 00372 delaunay* d = delaunay_build( nin, pin, 0, NULL, 0, NULL ); 00373 nnpi * nn = nnpi_create( d ); 00374 int seed = 0; 00375 int i; 00376 00377 nn->wmin = wmin; 00378 00379 if ( nn_verbose ) 00380 { 00381 fprintf( stderr, "xytoi:\n" ); 00382 for ( i = 0; i < nout; ++i ) 00383 { 00384 point* p = &pout[i]; 00385 00386 fprintf( stderr, "(%.7g,%.7g) -> %d\n", p->x, p->y, delaunay_xytoi( d, p, seed ) ); 00387 } 00388 } 00389 00390 for ( i = 0; i < nout; ++i ) 00391 nnpi_interpolate_point( nn, &pout[i] ); 00392 00393 if ( nn_verbose ) 00394 { 00395 fprintf( stderr, "output:\n" ); 00396 for ( i = 0; i < nout; ++i ) 00397 { 00398 point* p = &pout[i]; 00399 00400 fprintf( stderr, " %d:%15.7g %15.7g %15.7g\n", i, p->x, p->y, p->z ); 00401 } 00402 } 00403 00404 nnpi_destroy( nn ); 00405 delaunay_destroy( d ); 00406 } 00407 00408 // Sets minimal allowed weight for Natural Neighbours interpolation. 00409 // @param nn Natural Neighbours point interpolator 00410 // @param wmin Minimal allowed weight 00411 // 00412 void nnpi_setwmin( nnpi* nn, double wmin ) 00413 { 00414 nn->wmin = wmin; 00415 } 00416 00417 // Sets point to interpolate in. 00418 // @param nn Natural Neighbours point interpolator 00419 // @param p Point to interpolate in 00420 // 00421 void nnpi_set_point( nnpi* nn, point* p ) 00422 { 00423 nn->p = p; 00424 } 00425 00426 // Gets number of data points involved in current interpolation. 00427 // @return Number of data points involved in current interpolation 00428 // 00429 int nnpi_get_nvertices( nnpi* nn ) 00430 { 00431 return nn->nvertices; 00432 } 00433 00434 // Gets indices of data points involved in current interpolation. 00435 // @return indices of data points involved in current interpolation 00436 // 00437 int* nnpi_get_vertices( nnpi* nn ) 00438 { 00439 return nn->vertices; 00440 } 00441 00442 // Gets weights of data points involved in current interpolation. 00443 // @return weights of data points involved in current interpolation 00444 // 00445 double* nnpi_get_weights( nnpi* nn ) 00446 { 00447 return nn->weights; 00448 } 00449 00450 // 00451 // nnhpi 00452 // 00453 00454 struct nnhpi 00455 { 00456 nnpi * nnpi; 00457 hashtable* ht_data; 00458 hashtable* ht_weights; 00459 int n; // number of points processed 00460 }; 00461 00462 typedef struct 00463 { 00464 int nvertices; 00465 int * vertices; // vertex indices [nvertices] 00466 double* weights; // vertex weights [nvertices] 00467 } nn_weights; 00468 00469 // Creates Natural Neighbours hashing point interpolator. 00470 // 00471 // @param d Delaunay triangulation 00472 // @param size Hash table size (should be of order of number of output points) 00473 // @return Natural Neighbours interpolation 00474 // 00475 nnhpi* nnhpi_create( delaunay* d, int size ) 00476 { 00477 nnhpi* nn = malloc( sizeof ( nnhpi ) ); 00478 int i; 00479 00480 nn->nnpi = nnpi_create( d ); 00481 00482 nn->ht_data = ht_create_d2( d->npoints ); 00483 nn->ht_weights = ht_create_d2( size ); 00484 nn->n = 0; 00485 00486 for ( i = 0; i < d->npoints; ++i ) 00487 ht_insert( nn->ht_data, &d->points[i], &d->points[i] ); 00488 00489 return nn; 00490 } 00491 00492 static void free_nn_weights( void* data ) 00493 { 00494 nn_weights* weights = (nn_weights *) data; 00495 00496 free( weights->vertices ); 00497 free( weights->weights ); 00498 free( weights ); 00499 } 00500 00501 // Destroys Natural Neighbours hashing point interpolation. 00502 // 00503 // @param nn Structure to be destroyed 00504 // 00505 void nnhpi_destroy( nnhpi* nn ) 00506 { 00507 ht_destroy( nn->ht_data ); 00508 ht_process( nn->ht_weights, free_nn_weights ); 00509 ht_destroy( nn->ht_weights ); 00510 nnpi_destroy( nn->nnpi ); 00511 } 00512 00513 // Finds Natural Neighbours-interpolated value in a point. 00514 // 00515 // @param nnhp NN point hashing interpolator 00516 // @param p Point to be interpolated (p->x, p->y -- input; p->z -- output) 00517 // 00518 void nnhpi_interpolate( nnhpi* nnhp, point* p ) 00519 { 00520 nnpi * nnp = nnhp->nnpi; 00521 delaunay * d = nnp->d; 00522 hashtable * ht_weights = nnhp->ht_weights; 00523 nn_weights* weights; 00524 int i; 00525 00526 if ( ht_find( ht_weights, p ) != NULL ) 00527 { 00528 weights = ht_find( ht_weights, p ); 00529 if ( nn_verbose ) 00530 fprintf( stderr, " <hashtable>\n" ); 00531 } 00532 else 00533 { 00534 nnpi_reset( nnp ); 00535 nnp->p = p; 00536 nnpi_calculate_weights( nnp ); 00537 nnpi_normalize_weights( nnp ); 00538 00539 weights = malloc( sizeof ( nn_weights ) ); 00540 weights->vertices = malloc( sizeof ( int ) * (size_t) ( nnp->nvertices ) ); 00541 weights->weights = malloc( sizeof ( double ) * (size_t) ( nnp->nvertices ) ); 00542 00543 weights->nvertices = nnp->nvertices; 00544 00545 for ( i = 0; i < nnp->nvertices; ++i ) 00546 { 00547 weights->vertices[i] = nnp->vertices[i]; 00548 weights->weights[i] = nnp->weights[i]; 00549 } 00550 00551 ht_insert( ht_weights, p, weights ); 00552 00553 if ( nn_verbose ) 00554 { 00555 if ( nn_test_vertice == -1 ) 00556 { 00557 if ( nnp->n == 0 ) 00558 fprintf( stderr, "weights:\n" ); 00559 fprintf( stderr, " %d: {", nnp->n ); 00560 00561 for ( i = 0; i < nnp->nvertices; ++i ) 00562 { 00563 fprintf( stderr, "(%d,%.5g)", nnp->vertices[i], nnp->weights[i] ); 00564 00565 if ( i < nnp->nvertices - 1 ) 00566 fprintf( stderr, ", " ); 00567 } 00568 fprintf( stderr, "}\n" ); 00569 } 00570 else 00571 { 00572 double w = 0.0; 00573 00574 if ( nnp->n == 0 ) 00575 fprintf( stderr, "weights for vertex %d:\n", nn_test_vertice ); 00576 for ( i = 0; i < nnp->nvertices; ++i ) 00577 { 00578 if ( nnp->vertices[i] == nn_test_vertice ) 00579 { 00580 w = nnp->weights[i]; 00581 00582 break; 00583 } 00584 } 00585 fprintf( stderr, "%15.7g %15.7g %15.7g\n", p->x, p->y, w ); 00586 } 00587 } 00588 00589 nnp->n++; 00590 } 00591 00592 nnhp->n++; 00593 00594 if ( weights->nvertices == 0 ) 00595 { 00596 p->z = NaN; 00597 return; 00598 } 00599 00600 p->z = 0.0; 00601 for ( i = 0; i < weights->nvertices; ++i ) 00602 { 00603 if ( weights->weights[i] < nnp->wmin ) 00604 { 00605 p->z = NaN; 00606 return; 00607 } 00608 p->z += d->points[weights->vertices[i]].z * weights->weights[i]; 00609 } 00610 } 00611 00612 // Modifies interpolated data. 00613 // Finds point* pd in the underlying Delaunay triangulation such that 00614 // pd->x = p->x and pd->y = p->y, and copies p->z to pd->z. Exits with error 00615 // if the point is not found. 00616 // 00617 // @param nnhp Natural Neighbours hashing point interpolator 00618 // @param p New data 00619 // 00620 void nnhpi_modify_data( nnhpi* nnhp, point* p ) 00621 { 00622 point* orig = ht_find( nnhp->ht_data, p ); 00623 00624 assert( orig != NULL ); 00625 orig->z = p->z; 00626 } 00627 00628 // Sets minimal allowed weight for Natural Neighbours interpolation. 00629 // @param nn Natural Neighbours point hashing interpolator 00630 // @param wmin Minimal allowed weight 00631 // 00632 void nnhpi_setwmin( nnhpi* nn, double wmin ) 00633 { 00634 nn->nnpi->wmin = wmin; 00635 } 00636 00637 #if defined ( NNPHI_TEST ) 00638 00639 #include <sys/time.h> 00640 00641 #define NPOINTSIN 10000 00642 #define NMIN 10 00643 #define NX 101 00644 #define NXMIN 1 00645 00646 #define SQ( x ) ( ( x ) * ( x ) ) 00647 00648 static double franke( double x, double y ) 00649 { 00650 x *= 9.0; 00651 y *= 9.0; 00652 return 0.75 * exp( ( -SQ( x - 2.0 ) - SQ( y - 2.0 ) ) / 4.0 ) 00653 + 0.75 * exp( -SQ( x - 2.0 ) / 49.0 - ( y - 2.0 ) / 10.0 ) 00654 + 0.5 * exp( ( -SQ( x - 7.0 ) - SQ( y - 3.0 ) ) / 4.0 ) 00655 - 0.2 * exp( -SQ( x - 4.0 ) - SQ( y - 7.0 ) ); 00656 } 00657 00658 static void usage() 00659 { 00660 printf( "Usage: nnhpi_test [-a] [-n <nin> <nxout>] [-v|-V]\n" ); 00661 printf( "Options:\n" ); 00662 printf( " -a -- use non-Sibsonian interpolation rule\n" ); 00663 printf( " -n <nin> <nout>:\n" ); 00664 printf( " <nin> -- number of input points (default = 10000)\n" ); 00665 printf( " <nout> -- number of output points per side (default = 64)\n" ); 00666 printf( " -v -- verbose\n" ); 00667 printf( " -V -- very verbose\n" ); 00668 00669 exit( 0 ); 00670 } 00671 00672 int main( int argc, char* argv[] ) 00673 { 00674 int nin = NPOINTSIN; 00675 int nx = NX; 00676 int nout = 0; 00677 point * pin = NULL; 00678 delaunay * d = NULL; 00679 point * pout = NULL; 00680 nnhpi * nn = NULL; 00681 int cpi = -1; // control point index 00682 struct timeval tv0, tv1; 00683 struct timezone tz; 00684 int i; 00685 00686 i = 1; 00687 while ( i < argc ) 00688 { 00689 switch ( argv[i][1] ) 00690 { 00691 case 'a': 00692 i++; 00693 nn_rule = NON_SIBSONIAN; 00694 break; 00695 case 'n': 00696 i++; 00697 if ( i >= argc ) 00698 nn_quit( "no number of data points found after -n\n" ); 00699 nin = atoi( argv[i] ); 00700 i++; 00701 if ( i >= argc ) 00702 nn_quit( "no number of ouput points per side found after -i\n" ); 00703 nx = atoi( argv[i] ); 00704 i++; 00705 break; 00706 case 'v': 00707 i++; 00708 nn_verbose = 1; 00709 break; 00710 case 'V': 00711 i++; 00712 nn_verbose = 2; 00713 break; 00714 default: 00715 usage(); 00716 break; 00717 } 00718 } 00719 00720 if ( nin < NMIN ) 00721 nin = NMIN; 00722 if ( nx < NXMIN ) 00723 nx = NXMIN; 00724 00725 printf( "\nTest of Natural Neighbours hashing point interpolator:\n\n" ); 00726 printf( " %d data points\n", nin ); 00727 printf( " %d output points\n", nx * nx ); 00728 00729 // 00730 // generate data 00731 // 00732 printf( " generating data:\n" ); 00733 fflush( stdout ); 00734 pin = malloc( nin * sizeof ( point ) ); 00735 for ( i = 0; i < nin; ++i ) 00736 { 00737 point* p = &pin[i]; 00738 00739 p->x = (double) random() / RAND_MAX; 00740 p->y = (double) random() / RAND_MAX; 00741 p->z = franke( p->x, p->y ); 00742 if ( nn_verbose ) 00743 printf( " (%f, %f, %f)\n", p->x, p->y, p->z ); 00744 } 00745 00746 // 00747 // triangulate 00748 // 00749 printf( " triangulating:\n" ); 00750 fflush( stdout ); 00751 d = delaunay_build( nin, pin, 0, NULL, 0, NULL ); 00752 00753 // 00754 // generate output points 00755 // 00756 points_generate2( -0.1, 1.1, -0.1, 1.1, nx, nx, &nout, &pout ); 00757 cpi = ( nx / 2 ) * ( nx + 1 ); 00758 00759 gettimeofday( &tv0, &tz ); 00760 00761 // 00762 // create interpolator 00763 // 00764 printf( " creating interpolator:\n" ); 00765 fflush( stdout ); 00766 nn = nnhpi_create( d, nout ); 00767 00768 fflush( stdout ); 00769 gettimeofday( &tv1, &tz ); 00770 { 00771 long dt = 1000000 * ( tv1.tv_sec - tv0.tv_sec ) + tv1.tv_usec - tv0.tv_usec; 00772 00773 printf( " interpolator creation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout ); 00774 } 00775 00776 // 00777 // interpolate 00778 // 00779 printf( " interpolating:\n" ); 00780 fflush( stdout ); 00781 gettimeofday( &tv1, &tz ); 00782 for ( i = 0; i < nout; ++i ) 00783 { 00784 point* p = &pout[i]; 00785 00786 nnhpi_interpolate( nn, p ); 00787 if ( nn_verbose ) 00788 printf( " (%f, %f, %f)\n", p->x, p->y, p->z ); 00789 } 00790 00791 fflush( stdout ); 00792 gettimeofday( &tv0, &tz ); 00793 { 00794 long dt = 1000000.0 * ( tv0.tv_sec - tv1.tv_sec ) + tv0.tv_usec - tv1.tv_usec; 00795 00796 printf( " interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout ); 00797 } 00798 00799 if ( !nn_verbose ) 00800 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 ) ); 00801 00802 printf( " interpolating one more time:\n" ); 00803 fflush( stdout ); 00804 gettimeofday( &tv0, &tz ); 00805 for ( i = 0; i < nout; ++i ) 00806 { 00807 point* p = &pout[i]; 00808 00809 nnhpi_interpolate( nn, p ); 00810 if ( nn_verbose ) 00811 printf( " (%f, %f, %f)\n", p->x, p->y, p->z ); 00812 } 00813 00814 fflush( stdout ); 00815 gettimeofday( &tv1, &tz ); 00816 { 00817 long dt = 1000000.0 * ( tv1.tv_sec - tv0.tv_sec ) + tv1.tv_usec - tv0.tv_usec; 00818 00819 printf( " interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout ); 00820 } 00821 00822 if ( !nn_verbose ) 00823 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 ) ); 00824 00825 printf( " entering new data:\n" ); 00826 fflush( stdout ); 00827 for ( i = 0; i < nin; ++i ) 00828 { 00829 point* p = &pin[i]; 00830 00831 p->z = p->x * p->x - p->y * p->y; 00832 nnhpi_modify_data( nn, p ); 00833 if ( nn_verbose ) 00834 printf( " (%f, %f, %f)\n", p->x, p->y, p->z ); 00835 } 00836 00837 printf( " interpolating:\n" ); 00838 fflush( stdout ); 00839 gettimeofday( &tv1, &tz ); 00840 for ( i = 0; i < nout; ++i ) 00841 { 00842 point* p = &pout[i]; 00843 00844 nnhpi_interpolate( nn, p ); 00845 if ( nn_verbose ) 00846 printf( " (%f, %f, %f)\n", p->x, p->y, p->z ); 00847 } 00848 00849 fflush( stdout ); 00850 gettimeofday( &tv0, &tz ); 00851 { 00852 long dt = 1000000.0 * ( tv0.tv_sec - tv1.tv_sec ) + tv0.tv_usec - tv1.tv_usec; 00853 00854 printf( " interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout ); 00855 } 00856 00857 if ( !nn_verbose ) 00858 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 ); 00859 00860 printf( " restoring data:\n" ); 00861 fflush( stdout ); 00862 for ( i = 0; i < nin; ++i ) 00863 { 00864 point* p = &pin[i]; 00865 00866 p->z = franke( p->x, p->y ); 00867 nnhpi_modify_data( nn, p ); 00868 if ( nn_verbose ) 00869 printf( " (%f, %f, %f)\n", p->x, p->y, p->z ); 00870 } 00871 00872 printf( " interpolating:\n" ); 00873 fflush( stdout ); 00874 gettimeofday( &tv0, &tz ); 00875 for ( i = 0; i < nout; ++i ) 00876 { 00877 point* p = &pout[i]; 00878 00879 nnhpi_interpolate( nn, p ); 00880 if ( nn_verbose ) 00881 printf( " (%f, %f, %f)\n", p->x, p->y, p->z ); 00882 } 00883 00884 fflush( stdout ); 00885 gettimeofday( &tv1, &tz ); 00886 { 00887 long dt = 1000000.0 * ( tv1.tv_sec - tv0.tv_sec ) + tv1.tv_usec - tv0.tv_usec; 00888 00889 printf( " interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout ); 00890 } 00891 00892 if ( !nn_verbose ) 00893 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 ) ); 00894 00895 printf( " hashtable stats:\n" ); 00896 fflush( stdout ); 00897 { 00898 hashtable* ht = nn->ht_data; 00899 00900 printf( " input points: %d entries, %d table elements, %d filled elements\n", ht->n, ht->size, ht->nhash ); 00901 ht = nn->ht_weights; 00902 printf( " weights: %d entries, %d table elements, %d filled elements\n", ht->n, ht->size, ht->nhash ); 00903 } 00904 printf( "\n" ); 00905 00906 nnhpi_destroy( nn ); 00907 free( pout ); 00908 delaunay_destroy( d ); 00909 free( pin ); 00910 00911 return 0; 00912 } 00913 00914 #endif