PLplot 5.15.0
nnai.c
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// File: nnai.c
4//
5// Created: 15/11/2002
6//
7// Author: Pavel Sakov
8// CSIRO Marine Research
9//
10// Purpose: Code for:
11// -- Natural Neighbours Array Interpolator
12//
13// Description: `nnai' is a tructure for conducting
14// consequitive Natural Neighbours interpolations on a given
15// spatial data set in a given array of points. It allows to
16// modify Z coordinate of data in between interpolations.
17// `nnai' is the fastest of the three Natural
18// Neighbours interpolators in `nn' library.
19//
20// Revisions: None
21//
22//--------------------------------------------------------------------------
23
24#include <stdlib.h>
25#include <stdio.h>
26#include <string.h>
27#include <math.h>
28#include "nn.h"
29#include "delaunay.h"
30#include "nan.h"
31
32typedef struct
33{
35 int * vertices; // vertex indices [nvertices]
36 double* weights; // vertex weights [nvertices]
38
39struct nnai
40{
42 double wmin;
43 double n; // number of output points
44 double * x; // [n]
45 double * y; // [n]
47};
48
49void nn_quit( const char* format, ... );
51int nnpi_get_nvertices( nnpi* nn );
52int* nnpi_get_vertices( nnpi* nn );
53double* nnpi_get_weights( nnpi* nn );
55void nnpi_reset( nnpi* nn );
56void nnpi_set_point( nnpi* nn, point* p );
57
58// Builds Natural Neighbours array interpolator. This includes calculation of
59// weights used in nnai_interpolate().
60//
61// @param d Delaunay triangulation
62// @return Natural Neighbours interpolation
63//
64nnai* nnai_build( delaunay* d, int n, double* x, double* y )
65{
66 nnai * nn = malloc( sizeof ( nnai ) );
67 nnpi * nnp = nnpi_create( d );
68 int * vertices;
69 double* weights;
70 int i;
71
72 if ( n <= 0 )
73 nn_quit( "nnai_create(): n = %d\n", n );
74
75 nn->d = d;
76 nn->n = n;
77 nn->x = malloc( (size_t) n * sizeof ( double ) );
78 memcpy( nn->x, x, (size_t) n * sizeof ( double ) );
79 nn->y = malloc( (size_t) n * sizeof ( double ) );
80 memcpy( nn->y, y, (size_t) n * sizeof ( double ) );
81 nn->weights = malloc( (size_t) n * sizeof ( nn_weights ) );
82
83 for ( i = 0; i < n; ++i )
84 {
85 nn_weights* w = &nn->weights[i];
86 point p;
87
88 p.x = x[i];
89 p.y = y[i];
90
91 nnpi_reset( nnp );
92 nnpi_set_point( nnp, &p );
95
96 vertices = nnpi_get_vertices( nnp );
97 weights = nnpi_get_weights( nnp );
98
99 w->nvertices = nnpi_get_nvertices( nnp );
100 w->vertices = malloc( (size_t) ( w->nvertices ) * sizeof ( int ) );
101 memcpy( w->vertices, vertices, (size_t) ( w->nvertices ) * sizeof ( int ) );
102 w->weights = malloc( (size_t) ( w->nvertices ) * sizeof ( double ) );
103 memcpy( w->weights, weights, (size_t) ( w->nvertices ) * sizeof ( double ) );
104 }
105
106 nnpi_destroy( nnp );
107
108 return nn;
109}
110
111// Destroys Natural Neighbours array interpolator.
112//
113// @param nn Structure to be destroyed
114//
116{
117 int i;
118
119 for ( i = 0; i < nn->n; ++i )
120 {
121 nn_weights* w = &nn->weights[i];
122
123 free( w->vertices );
124 free( w->weights );
125 }
126
127 free( nn->x );
128 free( nn->y );
129 free( nn->weights );
130 free( nn );
131}
132
133// Conducts NN interpolation in a fixed array of output points using
134// data specified for a fixed array of input points. Uses pre-calculated
135// weights.
136//
137// @param nn NN array interpolator
138// @param zin input data [nn->d->npoints]
139// @param zout output data [nn->n]. Must be pre-allocated!
140//
141void nnai_interpolate( nnai* nn, double* zin, double* zout )
142{
143 int i;
144
145 for ( i = 0; i < nn->n; ++i )
146 {
147 nn_weights* w = &nn->weights[i];
148 double z = 0.0;
149 int j;
150
151 for ( j = 0; j < w->nvertices; ++j )
152 {
153 double weight = w->weights[j];
154
155 if ( weight < nn->wmin )
156 {
157 z = NaN;
158 break;
159 }
160 z += weight * zin[w->vertices[j]];
161 }
162
163 zout[i] = z;
164 }
165}
166
167//* Sets minimal allowed weight for Natural Neighbours interpolation.
168// @param nn Natural Neighbours array interpolator
169// @param wmin Minimal allowed weight
170//
171void nnai_setwmin( nnai* nn, double wmin )
172{
173 nn->wmin = wmin;
174}
175
176// The rest of this file contains a number of test programs.
177//
178#if defined ( NNAI_TEST )
179
180#include <sys/time.h>
181
182#define NPOINTSIN 10000
183#define NMIN 10
184#define NX 101
185#define NXMIN 1
186
187#define SQ( x ) ( ( x ) * ( x ) )
188
189static double franke( double x, double y )
190{
191 x *= 9.0;
192 y *= 9.0;
193 return 0.75 * exp( ( -SQ( x - 2.0 ) - SQ( y - 2.0 ) ) / 4.0 )
194 + 0.75 * exp( -SQ( x - 2.0 ) / 49.0 - ( y - 2.0 ) / 10.0 )
195 + 0.5 * exp( ( -SQ( x - 7.0 ) - SQ( y - 3.0 ) ) / 4.0 )
196 - 0.2 * exp( -SQ( x - 4.0 ) - SQ( y - 7.0 ) );
197}
198
199static void usage()
200{
201 printf(
202 "Usage: nn_test [-v|-V] [-n <nin> <nxout>]\n"
203 "Options:\n"
204 " -a -- use non-Sibsonian interpolation rule\n"
205 " -n <nin> <nout>:\n"
206 " <nin> -- number of input points (default = 10000)\n"
207 " <nout> -- number of output points per side (default = 64)\n"
208 " -v -- verbose\n"
209 " -V -- very verbose\n"
210 );
211}
212
213int main( int argc, char* argv[] )
214{
215 int nin = NPOINTSIN;
216 int nx = NX;
217 int nout = 0;
218 point * pin = NULL;
219 delaunay * d = NULL;
220 point * pout = NULL;
221 nnai * nn = NULL;
222 double * zin = NULL;
223 double * xout = NULL;
224 double * yout = NULL;
225 double * zout = NULL;
226 int cpi = -1; // control point index
227 struct timeval tv0, tv1, tv2;
228 struct timezone tz;
229 int i;
230
231 i = 1;
232 while ( i < argc )
233 {
234 switch ( argv[i][1] )
235 {
236 case 'a':
237 i++;
239 break;
240 case 'n':
241 i++;
242 if ( i >= argc )
243 nn_quit( "no number of data points found after -i\n" );
244 nin = atoi( argv[i] );
245 i++;
246 if ( i >= argc )
247 nn_quit( "no number of ouput points per side found after -i\n" );
248 nx = atoi( argv[i] );
249 i++;
250 break;
251 case 'v':
252 i++;
253 nn_verbose = 1;
254 break;
255 case 'V':
256 i++;
257 nn_verbose = 2;
258 break;
259 default:
260 usage();
261 break;
262 }
263 }
264
265 if ( nin < NMIN )
266 nin = NMIN;
267 if ( nx < NXMIN )
268 nx = NXMIN;
269
270 printf( "\nTest of Natural Neighbours array interpolator:\n\n" );
271 printf( " %d data points\n", nin );
272 printf( " %d output points\n", nx * nx );
273
274 //
275 // generate data
276 //
277 printf( " generating data:\n" );
278 fflush( stdout );
279 pin = malloc( nin * sizeof ( point ) );
280 zin = malloc( nin * sizeof ( double ) );
281 for ( i = 0; i < nin; ++i )
282 {
283 point* p = &pin[i];
284
285 p->x = (double) random() / RAND_MAX;
286 p->y = (double) random() / RAND_MAX;
287 p->z = franke( p->x, p->y );
288 zin[i] = p->z;
289 if ( nn_verbose )
290 printf( " (%f, %f, %f)\n", p->x, p->y, p->z );
291 }
292
293 //
294 // triangulate
295 //
296 printf( " triangulating:\n" );
297 fflush( stdout );
298 d = delaunay_build( nin, pin, 0, NULL, 0, NULL );
299
300 //
301 // generate output points
302 //
303 points_generate2( -0.1, 1.1, -0.1, 1.1, nx, nx, &nout, &pout );
304 xout = malloc( nout * sizeof ( double ) );
305 yout = malloc( nout * sizeof ( double ) );
306 zout = malloc( nout * sizeof ( double ) );
307 for ( i = 0; i < nout; ++i )
308 {
309 point* p = &pout[i];
310
311 xout[i] = p->x;
312 yout[i] = p->y;
313 zout[i] = NaN;
314 }
315 cpi = ( nx / 2 ) * ( nx + 1 );
316
317 gettimeofday( &tv0, &tz );
318
319 //
320 // create interpolator
321 //
322 printf( " creating interpolator:\n" );
323 fflush( stdout );
324 nn = nnai_build( d, nout, xout, yout );
325
326 fflush( stdout );
327 gettimeofday( &tv1, &tz );
328 {
329 long dt = 1000000 * ( tv1.tv_sec - tv0.tv_sec ) + tv1.tv_usec - tv0.tv_usec;
330
331 printf( " interpolator creation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout );
332 }
333
334 //
335 // interpolate
336 //
337 printf( " interpolating:\n" );
338 fflush( stdout );
339 nnai_interpolate( nn, zin, zout );
340 if ( nn_verbose )
341 for ( i = 0; i < nout; ++i )
342 printf( " (%f, %f, %f)\n", xout[i], yout[i], zout[i] );
343
344 fflush( stdout );
345 gettimeofday( &tv2, &tz );
346 {
347 long dt = 1000000.0 * ( tv2.tv_sec - tv1.tv_sec ) + tv2.tv_usec - tv1.tv_usec;
348
349 printf( " interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout );
350 }
351
352 if ( !nn_verbose )
353 printf( " control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], franke( xout[cpi], yout[cpi] ) );
354
355 printf( " interpolating one more time:\n" );
356 fflush( stdout );
357 nnai_interpolate( nn, zin, zout );
358 if ( nn_verbose )
359 for ( i = 0; i < nout; ++i )
360 printf( " (%f, %f, %f)\n", xout[i], yout[i], zout[i] );
361
362 fflush( stdout );
363 gettimeofday( &tv0, &tz );
364 {
365 long dt = 1000000.0 * ( tv0.tv_sec - tv2.tv_sec ) + tv0.tv_usec - tv2.tv_usec;
366
367 printf( " interpolation time = %ld us (%.2f us / point)\n", dt, (double) dt / nout );
368 }
369
370 if ( !nn_verbose )
371 printf( " control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], franke( xout[cpi], yout[cpi] ) );
372
373 printf( " entering new data:\n" );
374 fflush( stdout );
375 for ( i = 0; i < nin; ++i )
376 {
377 point* p = &pin[i];
378
379 p->z = p->x * p->x - p->y * p->y;
380 zin[i] = p->z;
381 if ( nn_verbose )
382 printf( " (%f, %f, %f)\n", p->x, p->y, p->z );
383 }
384
385 printf( " interpolating:\n" );
386 fflush( stdout );
387 nnai_interpolate( nn, zin, zout );
388 if ( nn_verbose )
389 for ( i = 0; i < nout; ++i )
390 printf( " (%f, %f, %f)\n", xout[i], yout[i], zout[i] );
391
392 if ( !nn_verbose )
393 printf( " control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], xout[cpi] * xout[cpi] - yout[cpi] * yout[cpi] );
394
395 printf( " restoring data:\n" );
396 fflush( stdout );
397 for ( i = 0; i < nin; ++i )
398 {
399 point* p = &pin[i];
400
401 p->z = franke( p->x, p->y );
402 zin[i] = p->z;
403 if ( nn_verbose )
404 printf( " (%f, %f, %f)\n", p->x, p->y, p->z );
405 }
406
407 printf( " interpolating:\n" );
408 fflush( stdout );
409 nnai_interpolate( nn, zin, zout );
410 if ( nn_verbose )
411 for ( i = 0; i < nout; ++i )
412 printf( " (%f, %f, %f)\n", xout[i], yout[i], zout[i] );
413
414 if ( !nn_verbose )
415 printf( " control point: (%f, %f, %f) (expected z = %f)\n", xout[cpi], yout[cpi], zout[cpi], franke( xout[cpi], yout[cpi] ) );
416
417 printf( "\n" );
418
419 nnai_destroy( nn );
420 free( zin );
421 free( xout );
422 free( yout );
423 free( zout );
424 free( pout );
425 delaunay_destroy( d );
426 free( pin );
427
428 return 0;
429}
430
431#endif
int main()
Definition: cdexpert.c:35
#define NaN
Definition: csa/nan.h:62
delaunay * delaunay_build(int np, point points[], int ns, int segments[], int nh, double holes[])
Definition: delaunay.c:265
void delaunay_destroy(delaunay *d)
Definition: delaunay.c:578
void points_generate2(double xmin, double xmax, double ymin, double ymax, int nx, int ny, int *nout, point **pout)
Definition: nncommon.c:334
@ NON_SIBSONIAN
Definition: nn.h:23
int nn_verbose
Definition: nncommon.c:43
void nnpi_destroy(nnpi *nn)
Definition: nnpi.c:115
nnpi * nnpi_create(delaunay *d)
Definition: nnpi.c:95
nnai * nnai_build(delaunay *d, int n, double *x, double *y)
Definition: nnai.c:64
void nnpi_set_point(nnpi *nn, point *p)
Definition: nnpi.c:421
void nnpi_normalize_weights(nnpi *nn)
Definition: nnpi.c:279
int * nnpi_get_vertices(nnpi *nn)
Definition: nnpi.c:437
void nnai_destroy(nnai *nn)
Definition: nnai.c:115
int nnpi_get_nvertices(nnpi *nn)
Definition: nnpi.c:429
void nn_quit(const char *format,...)
Definition: nncommon.c:53
void nnai_interpolate(nnai *nn, double *zin, double *zout)
Definition: nnai.c:141
void nnai_setwmin(nnai *nn, double wmin)
Definition: nnai.c:171
double * nnpi_get_weights(nnpi *nn)
Definition: nnpi.c:445
void nnpi_reset(nnpi *nn)
Definition: nnpi.c:122
void nnpi_calculate_weights(nnpi *nn)
Definition: nnpi.c:259
NN_RULE nn_rule
Definition: nncommon.c:45
static PLCHAR_VECTOR usage
Definition: plargs.c:179
#define NX
static int argc
Definition: qt.cpp:48
static char ** argv
Definition: qt.cpp:49
int nvertices
Definition: nnai.c:34
int * vertices
Definition: nnai.c:35
double * weights
Definition: nnai.c:36
Definition: nnai.c:40
double * x
Definition: nnai.c:44
delaunay * d
Definition: nnai.c:41
double n
Definition: nnai.c:43
nn_weights * weights
Definition: nnai.c:46
double * y
Definition: nnai.c:45
double wmin
Definition: nnai.c:42
Definition: nnpi.c:54
Definition: csa.h:30
double y
Definition: csa.h:32
double x
Definition: csa.h:31
double z
Definition: csa.h:33