PLplot  5.10.0
nnai.c
Go to the documentation of this file.
00001 //--------------------------------------------------------------------------
00002 //
00003 // File:           nnai.c
00004 //
00005 // Created:        15/11/2002
00006 //
00007 // Author:         Pavel Sakov
00008 //                 CSIRO Marine Research
00009 //
00010 // Purpose:        Code for:
00011 //                 -- Natural Neighbours Array Interpolator
00012 //
00013 // Description:    `nnai' is a tructure for conducting
00014 //                 consequitive Natural Neighbours interpolations on a given
00015 //                 spatial data set in a given array of points. It allows to
00016 //                 modify Z coordinate of data in between interpolations.
00017 //                 `nnai' is the fastest of the three Natural
00018 //                 Neighbours interpolators in `nn' library.
00019 //
00020 // Revisions:      None
00021 //
00022 //--------------------------------------------------------------------------
00023 
00024 #include <stdlib.h>
00025 #include <stdio.h>
00026 #include <string.h>
00027 #include <math.h>
00028 #include "nn.h"
00029 #include "delaunay.h"
00030 #include "nan.h"
00031 
00032 typedef struct
00033 {
00034     int   nvertices;
00035     int   * vertices;           // vertex indices [nvertices]
00036     double* weights;            // vertex weights [nvertices]
00037 } nn_weights;
00038 
00039 struct nnai
00040 {
00041     delaunay  * d;
00042     double    wmin;
00043     double    n;                // number of output points
00044     double    * x;              // [n]
00045     double    * y;              // [n]
00046     nn_weights* weights;
00047 };
00048 
00049 void nn_quit( const char* format, ... );
00050 void nnpi_calculate_weights( nnpi* nn );
00051 int nnpi_get_nvertices( nnpi* nn );
00052 int* nnpi_get_vertices( nnpi* nn );
00053 double* nnpi_get_weights( nnpi* nn );
00054 void nnpi_normalize_weights( nnpi* nn );
00055 void nnpi_reset( nnpi* nn );
00056 void nnpi_set_point( nnpi* nn, point* p );
00057 
00058 // Builds Natural Neighbours array interpolator. This includes calculation of
00059 // weights used in nnai_interpolate().
00060 //
00061 // @param d Delaunay triangulation
00062 // @return Natural Neighbours interpolation
00063 //
00064 nnai* nnai_build( delaunay* d, int n, double* x, double* y )
00065 {
00066     nnai  * nn  = malloc( sizeof ( nnai ) );
00067     nnpi  * nnp = nnpi_create( d );
00068     int   * vertices;
00069     double* weights;
00070     int   i;
00071 
00072     if ( n <= 0 )
00073         nn_quit( "nnai_create(): n = %d\n", n );
00074 
00075     nn->d = d;
00076     nn->n = n;
00077     nn->x = malloc( (size_t) n * sizeof ( double ) );
00078     memcpy( nn->x, x, (size_t) n * sizeof ( double ) );
00079     nn->y = malloc( (size_t) n * sizeof ( double ) );
00080     memcpy( nn->y, y, (size_t) n * sizeof ( double ) );
00081     nn->weights = malloc( (size_t) n * sizeof ( nn_weights ) );
00082 
00083     for ( i = 0; i < n; ++i )
00084     {
00085         nn_weights* w = &nn->weights[i];
00086         point     p;
00087 
00088         p.x = x[i];
00089         p.y = y[i];
00090 
00091         nnpi_reset( nnp );
00092         nnpi_set_point( nnp, &p );
00093         nnpi_calculate_weights( nnp );
00094         nnpi_normalize_weights( nnp );
00095 
00096         vertices = nnpi_get_vertices( nnp );
00097         weights  = nnpi_get_weights( nnp );
00098 
00099         w->nvertices = nnpi_get_nvertices( nnp );
00100         w->vertices  = malloc( (size_t) ( w->nvertices ) * sizeof ( int ) );
00101         memcpy( w->vertices, vertices, (size_t) ( w->nvertices ) * sizeof ( int ) );
00102         w->weights = malloc( (size_t) ( w->nvertices ) * sizeof ( double ) );
00103         memcpy( w->weights, weights, (size_t) ( w->nvertices ) * sizeof ( double ) );
00104     }
00105 
00106     nnpi_destroy( nnp );
00107 
00108     return nn;
00109 }
00110 
00111 // Destroys Natural Neighbours array interpolator.
00112 //
00113 // @param nn Structure to be destroyed
00114 //
00115 void nnai_destroy( nnai* nn )
00116 {
00117     int i;
00118 
00119     for ( i = 0; i < nn->n; ++i )
00120     {
00121         nn_weights* w = &nn->weights[i];
00122 
00123         free( w->vertices );
00124         free( w->weights );
00125     }
00126 
00127     free( nn->x );
00128     free( nn->y );
00129     free( nn->weights );
00130     free( nn );
00131 }
00132 
00133 // Conducts NN interpolation in a fixed array of output points using
00134 // data specified for a fixed array of input points. Uses pre-calculated
00135 // weights.
00136 //
00137 // @param nn NN array interpolator
00138 // @param zin input data [nn->d->npoints]
00139 // @param zout output data [nn->n]. Must be pre-allocated!
00140 //
00141 void nnai_interpolate( nnai* nn, double* zin, double* zout )
00142 {
00143     int i;
00144 
00145     for ( i = 0; i < nn->n; ++i )
00146     {
00147         nn_weights* w = &nn->weights[i];
00148         double    z   = 0.0;
00149         int       j;
00150 
00151         for ( j = 0; j < w->nvertices; ++j )
00152         {
00153             double weight = w->weights[j];
00154 
00155             if ( weight < nn->wmin )
00156             {
00157                 z = NaN;
00158                 break;
00159             }
00160             z += weight * zin[w->vertices[j]];
00161         }
00162 
00163         zout[i] = z;
00164     }
00165 }
00166 
00167 //* Sets minimal allowed weight for Natural Neighbours interpolation.
00168 // @param nn Natural Neighbours array interpolator
00169 // @param wmin Minimal allowed weight
00170 //
00171 void nnai_setwmin( nnai* nn, double wmin )
00172 {
00173     nn->wmin = wmin;
00174 }
00175 
00176 // The rest of this file contains a number of test programs.
00177 //
00178 #if defined ( NNAI_TEST )
00179 
00180 #include <sys/time.h>
00181 
00182 #define NPOINTSIN    10000
00183 #define NMIN         10
00184 #define NX           101
00185 #define NXMIN        1
00186 
00187 #define SQ( x )    ( ( x ) * ( x ) )
00188 
00189 static double franke( double x, double y )
00190 {
00191     x *= 9.0;
00192     y *= 9.0;
00193     return 0.75 * exp( ( -SQ( x - 2.0 ) - SQ( y - 2.0 ) ) / 4.0 )
00194            + 0.75 * exp( -SQ( x - 2.0 ) / 49.0 - ( y - 2.0 ) / 10.0 )
00195            + 0.5 * exp( ( -SQ( x - 7.0 ) - SQ( y - 3.0 ) ) / 4.0 )
00196            - 0.2 * exp( -SQ( x - 4.0 ) - SQ( y - 7.0 ) );
00197 }
00198 
00199 static void usage()
00200 {
00201     printf(
00202         "Usage: nn_test [-v|-V] [-n <nin> <nxout>]\n"
00203         "Options:\n"
00204         "  -a              -- use non-Sibsonian interpolation rule\n"
00205         "  -n <nin> <nout>:\n"
00206         "            <nin> -- number of input points (default = 10000)\n"
00207         "           <nout> -- number of output points per side (default = 64)\n"
00208         "  -v              -- verbose\n"
00209         "  -V              -- very verbose\n"
00210         );
00211 }
00212 
00213 int main( int argc, char* argv[] )
00214 {
00215     int nin  = NPOINTSIN;
00216     int nx   = NX;
00217     int nout = 0;
00218     point           * pin  = NULL;
00219     delaunay        * d    = NULL;
00220     point           * pout = NULL;
00221     nnai            * nn   = NULL;
00222     double          * zin  = NULL;
00223     double          * xout = NULL;
00224     double          * yout = NULL;
00225     double          * zout = NULL;
00226     int             cpi    = -1; // control point index
00227     struct timeval  tv0, tv1, tv2;
00228     struct timezone tz;
00229     int             i;
00230 
00231     i = 1;
00232     while ( i < argc )
00233     {
00234         switch ( argv[i][1] )
00235         {
00236         case 'a':
00237             i++;
00238             nn_rule = NON_SIBSONIAN;
00239             break;
00240         case 'n':
00241             i++;
00242             if ( i >= argc )
00243                 nn_quit( "no number of data points found after -i\n" );
00244             nin = atoi( argv[i] );
00245             i++;
00246             if ( i >= argc )
00247                 nn_quit( "no number of ouput points per side found after -i\n" );
00248             nx = atoi( argv[i] );
00249             i++;
00250             break;
00251         case 'v':
00252             i++;
00253             nn_verbose = 1;
00254             break;
00255         case 'V':
00256             i++;
00257             nn_verbose = 2;
00258             break;
00259         default:
00260             usage();
00261             break;
00262         }
00263     }
00264 
00265     if ( nin < NMIN )
00266         nin = NMIN;
00267     if ( nx < NXMIN )
00268         nx = NXMIN;
00269 
00270     printf( "\nTest of Natural Neighbours array interpolator:\n\n" );
00271     printf( "  %d data points\n", nin );
00272     printf( "  %d output points\n", nx * nx );
00273 
00274     //
00275     // generate data
00276     //
00277     printf( "  generating data:\n" );
00278     fflush( stdout );
00279     pin = malloc( nin * sizeof ( point ) );
00280     zin = malloc( nin * sizeof ( double ) );
00281     for ( i = 0; i < nin; ++i )
00282     {
00283         point* p = &pin[i];
00284 
00285         p->x   = (double) random() / RAND_MAX;
00286         p->y   = (double) random() / RAND_MAX;
00287         p->z   = franke( p->x, p->y );
00288         zin[i] = p->z;
00289         if ( nn_verbose )
00290             printf( "    (%f, %f, %f)\n", p->x, p->y, p->z );
00291     }
00292 
00293     //
00294     // triangulate
00295     //
00296     printf( "  triangulating:\n" );
00297     fflush( stdout );
00298     d = delaunay_build( nin, pin, 0, NULL, 0, NULL );
00299 
00300     //
00301     // generate output points
00302     //
00303     points_generate2( -0.1, 1.1, -0.1, 1.1, nx, nx, &nout, &pout );
00304     xout = malloc( nout * sizeof ( double ) );
00305     yout = malloc( nout * sizeof ( double ) );
00306     zout = malloc( nout * sizeof ( double ) );
00307     for ( i = 0; i < nout; ++i )
00308     {
00309         point* p = &pout[i];
00310 
00311         xout[i] = p->x;
00312         yout[i] = p->y;
00313         zout[i] = NaN;
00314     }
00315     cpi = ( nx / 2 ) * ( nx + 1 );
00316 
00317     gettimeofday( &tv0, &tz );
00318 
00319     //
00320     // create interpolator
00321     //
00322     printf( "  creating interpolator:\n" );
00323     fflush( stdout );
00324     nn = nnai_build( d, nout, xout, yout );
00325 
00326     fflush( stdout );
00327     gettimeofday( &tv1, &tz );
00328     {
00329         long dt = 1000000 * ( tv1.tv_sec - tv0.tv_sec ) + tv1.tv_usec - tv0.tv_usec;
00330 
00331         printf( "    interpolator creation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout );
00332     }
00333 
00334     //
00335     // interpolate
00336     //
00337     printf( "  interpolating:\n" );
00338     fflush( stdout );
00339     nnai_interpolate( nn, zin, zout );
00340     if ( nn_verbose )
00341         for ( i = 0; i < nout; ++i )
00342             printf( "    (%f, %f, %f)\n", xout[i], yout[i], zout[i] );
00343 
00344     fflush( stdout );
00345     gettimeofday( &tv2, &tz );
00346     {
00347         long dt = 1000000.0 * ( tv2.tv_sec - tv1.tv_sec ) + tv2.tv_usec - tv1.tv_usec;
00348 
00349         printf( "    interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout );
00350     }
00351 
00352     if ( !nn_verbose )
00353         printf( "    control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], franke( xout[cpi], yout[cpi] ) );
00354 
00355     printf( "  interpolating one more time:\n" );
00356     fflush( stdout );
00357     nnai_interpolate( nn, zin, zout );
00358     if ( nn_verbose )
00359         for ( i = 0; i < nout; ++i )
00360             printf( "    (%f, %f, %f)\n", xout[i], yout[i], zout[i] );
00361 
00362     fflush( stdout );
00363     gettimeofday( &tv0, &tz );
00364     {
00365         long dt = 1000000.0 * ( tv0.tv_sec - tv2.tv_sec ) + tv0.tv_usec - tv2.tv_usec;
00366 
00367         printf( "    interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout );
00368     }
00369 
00370     if ( !nn_verbose )
00371         printf( "    control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], franke( xout[cpi], yout[cpi] ) );
00372 
00373     printf( "  entering new data:\n" );
00374     fflush( stdout );
00375     for ( i = 0; i < nin; ++i )
00376     {
00377         point* p = &pin[i];
00378 
00379         p->z   = p->x * p->x - p->y * p->y;
00380         zin[i] = p->z;
00381         if ( nn_verbose )
00382             printf( "    (%f, %f, %f)\n", p->x, p->y, p->z );
00383     }
00384 
00385     printf( "  interpolating:\n" );
00386     fflush( stdout );
00387     nnai_interpolate( nn, zin, zout );
00388     if ( nn_verbose )
00389         for ( i = 0; i < nout; ++i )
00390             printf( "    (%f, %f, %f)\n", xout[i], yout[i], zout[i] );
00391 
00392     if ( !nn_verbose )
00393         printf( "    control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], xout[cpi] * xout[cpi] - yout[cpi] * yout[cpi] );
00394 
00395     printf( "  restoring data:\n" );
00396     fflush( stdout );
00397     for ( i = 0; i < nin; ++i )
00398     {
00399         point* p = &pin[i];
00400 
00401         p->z   = franke( p->x, p->y );
00402         zin[i] = p->z;
00403         if ( nn_verbose )
00404             printf( "    (%f, %f, %f)\n", p->x, p->y, p->z );
00405     }
00406 
00407     printf( "  interpolating:\n" );
00408     fflush( stdout );
00409     nnai_interpolate( nn, zin, zout );
00410     if ( nn_verbose )
00411         for ( i = 0; i < nout; ++i )
00412             printf( "    (%f, %f, %f)\n", xout[i], yout[i], zout[i] );
00413 
00414     if ( !nn_verbose )
00415         printf( "    control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], franke( xout[cpi], yout[cpi] ) );
00416 
00417     printf( "\n" );
00418 
00419     nnai_destroy( nn );
00420     free( zin );
00421     free( xout );
00422     free( yout );
00423     free( zout );
00424     free( pout );
00425     delaunay_destroy( d );
00426     free( pin );
00427 
00428     return 0;
00429 }
00430 
00431 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines