PLplot
5.10.0
|
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