PLplot  5.10.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
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 
32 typedef struct
33 {
34  int nvertices;
35  int * vertices; // vertex indices [nvertices]
36  double* weights; // vertex weights [nvertices]
37 } nn_weights;
38 
39 struct nnai
40 {
42  double wmin;
43  double n; // number of output points
44  double * x; // [n]
45  double * y; // [n]
47 };
48 
49 void nn_quit( const char* format, ... );
50 void nnpi_calculate_weights( nnpi* nn );
51 int nnpi_get_nvertices( nnpi* nn );
52 int* nnpi_get_vertices( nnpi* nn );
53 double* nnpi_get_weights( nnpi* nn );
54 void nnpi_normalize_weights( nnpi* nn );
55 void nnpi_reset( nnpi* nn );
56 void 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 //
64 nnai* 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 //
115 void nnai_destroy( nnai* nn )
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 //
141 void 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 //
171 void 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 
189 static 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 
199 static 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 
213 int 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