PLplot  5.10.0
nnpi.c
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines