PLplot  5.10.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
nncommon.c
Go to the documentation of this file.
1 //--------------------------------------------------------------------------
2 //
3 // File: nncommon.c
4 //
5 // Created: 04/08/2000
6 //
7 // Author: Pavel Sakov
8 // CSIRO Marine Research
9 //
10 // Purpose: Common stuff for NN interpolation library
11 //
12 // Description: None
13 //
14 // Revisions: 15/11/2002 PS: Changed name from "utils.c"
15 // 28/02/2003 PS: Modified points_read() to do the job without
16 // rewinding the file. This allows to read from stdin when
17 // necessary.
18 // 09/04/2003 PS: Modified points_read() to read from a
19 // file specified by name, not by handle.
20 // Modified: Andrew Ross 20/10/2008
21 // Change <= comparison in circle_contains() to use EPSILON
22 // to catch case where the point lies on the circle and there
23 // is floating point rounding error in the radii.
24 //
25 //--------------------------------------------------------------------------
26 
27 #include <stdlib.h>
28 #include <stdio.h>
29 #include <stdarg.h>
30 #include <assert.h>
31 #include <math.h>
32 #include <limits.h>
33 #include <float.h>
34 #include <string.h>
35 #include <errno.h>
36 #include "nan.h"
37 #include "delaunay.h"
38 
39 #define BUFSIZE 1024
40 
41 #define EPSILON 1.0e-8
42 
43 int nn_verbose = 0;
44 int nn_test_vertice = -1;
46 
47 #include "version.h"
48 
49 void nn_quit( const char* format, ... );
50 int circle_build( circle* c, point* p1, point* p2, point* p3 );
51 int circle_contains( circle* c, point* p );
52 
53 void nn_quit( const char* format, ... )
54 {
55  va_list args;
56 
57  fflush( stdout ); // just in case, to have the exit message
58  // last
59 
60  fprintf( stderr, "error: nn: " );
61  va_start( args, format );
62  vfprintf( stderr, format, args );
63  va_end( args );
64 
65  exit( 1 );
66 }
67 
68 int circle_build( circle* c, point* p1, point* p2, point* p3 )
69 {
70  double x1sq = p1->x * p1->x;
71  double x2sq = p2->x * p2->x;
72  double x3sq = p3->x * p3->x;
73  double y1sq = p1->y * p1->y;
74  double y2sq = p2->y * p2->y;
75  double y3sq = p3->y * p3->y;
76  double t1 = x3sq - x2sq + y3sq - y2sq;
77  double t2 = x1sq - x3sq + y1sq - y3sq;
78  double t3 = x2sq - x1sq + y2sq - y1sq;
79  double D = ( p1->x * ( p2->y - p3->y ) + p2->x * ( p3->y - p1->y ) + p3->x * ( p1->y - p2->y ) ) * 2.0;
80 
81  if ( D == 0.0 )
82  return 0;
83 
84  c->x = ( p1->y * t1 + p2->y * t2 + p3->y * t3 ) / D;
85  c->y = -( p1->x * t1 + p2->x * t2 + p3->x * t3 ) / D;
86  c->r = hypot( c->x - p1->x, c->y - p1->y );
87 
88  return 1;
89 }
90 
91 // This procedure has taken it final shape after a number of tries. The problem
92 // was to have the calculated and stored radii being the same if (x,y) is
93 // exactly on the circle border (i.e. not to use FCPU extended precision in
94 // the radius calculation). This may have little effect in practice but was
95 // important in some tests when both input and output data were placed
96 // in rectangular grid nodes.
97 //
99 {
100  return hypot( c->x - p->x, c->y - p->y ) <= c->r * ( 1.0 + EPSILON );
101 }
102 
103 // Smoothes the input point array by averaging the input x,y and z values
104 // for each cell within virtual rectangular nx by ny grid. The corners of the
105 // grid are created from min and max values of the input array. It also frees
106 // the original array and returns results and new dimension via original
107 // data and size pointers.
108 //
109 // @param pn Pointer to number of points (input/output)
110 // @param ppoints Pointer to array of points (input/output) [*pn]
111 // @param nx Number of x nodes in decimation
112 // @param ny Number of y nodes in decimation
113 //
114 void points_thin( int* pn, point** ppoints, int nx, int ny )
115 {
116  int n = *pn;
117  point * points = *ppoints;
118  double xmin = DBL_MAX;
119  double xmax = -DBL_MAX;
120  double ymin = DBL_MAX;
121  double ymax = -DBL_MAX;
122  int nxy = nx * ny;
123  double * sumx = calloc( (size_t) nxy, sizeof ( double ) );
124  double * sumy = calloc( (size_t) nxy, sizeof ( double ) );
125  double * sumz = calloc( (size_t) nxy, sizeof ( double ) );
126  int * count = calloc( (size_t) nxy, sizeof ( int ) );
127  double stepx = 0.0;
128  double stepy = 0.0;
129  int nnew = 0;
130  point * pointsnew = NULL;
131  int i, j, ii;
132 
133  if ( nn_verbose )
134  fprintf( stderr, "thinned: %d points -> ", *pn );
135 
136  if ( nx < 1 || ny < 1 )
137  {
138  free( points );
139  *ppoints = NULL;
140  *pn = 0;
141  if ( nn_verbose )
142  fprintf( stderr, "0 points" );
143  free( sumx );
144  free( sumy );
145  free( sumz );
146  free( count );
147  return;
148  }
149 
150  for ( ii = 0; ii < n; ++ii )
151  {
152  point* p = &points[ii];
153 
154  if ( p->x < xmin )
155  xmin = p->x;
156  if ( p->x > xmax )
157  xmax = p->x;
158  if ( p->y < ymin )
159  ymin = p->y;
160  if ( p->y > ymax )
161  ymax = p->y;
162  }
163 
164  stepx = ( nx > 1 ) ? ( xmax - xmin ) / nx : 0.0;
165  stepy = ( ny > 1 ) ? ( ymax - ymin ) / ny : 0.0;
166 
167  for ( ii = 0; ii < n; ++ii )
168  {
169  point* p = &points[ii];
170  int index;
171 
172  //
173  // Following is the portion of the code which really depends on the
174  // floating point particulars. Do not be surprised if different
175  // compilers/options give different results here.
176  //
177  i = ( nx == 1 ) ? 0 : (int) ( ( p->x - xmin ) / stepx );
178  j = ( ny == 1 ) ? 0 : (int) ( ( p->y - ymin ) / stepy );
179 
180  if ( i == nx )
181  i--;
182  if ( j == ny )
183  j--;
184  index = i + j * nx;
185  sumx[index] += p->x;
186  sumy[index] += p->y;
187  sumz[index] += p->z;
188  count[index]++;
189  }
190 
191  for ( j = 0; j < ny; ++j )
192  {
193  for ( i = 0; i < nx; ++i )
194  {
195  int index = i + j * nx;
196 
197  if ( count[index] > 0 )
198  nnew++;
199  }
200  }
201 
202  pointsnew = malloc( (size_t) nnew * sizeof ( point ) );
203 
204  ii = 0;
205  for ( j = 0; j < ny; ++j )
206  {
207  for ( i = 0; i < nx; ++i )
208  {
209  int index = i + j * nx;
210  int nn = count[index];
211 
212  if ( nn > 0 )
213  {
214  point* p = &pointsnew[ii];
215 
216  p->x = sumx[index] / nn;
217  p->y = sumy[index] / nn;
218  p->z = sumz[index] / nn;
219  ii++;
220  }
221  }
222  }
223 
224  if ( nn_verbose )
225  fprintf( stderr, "%d points\n", nnew );
226 
227  free( sumx );
228  free( sumy );
229  free( sumz );
230  free( count );
231 
232  free( points );
233  *ppoints = pointsnew;
234  *pn = nnew;
235 }
236 
237 // Generates rectangular grid nx by ny using min and max x and y values from
238 // the input point array. Allocates space for the output point array, be sure
239 // to free it when necessary!
240 //
241 // @param n Number of points
242 // @param points Array of points [n]
243 // @param nx Number of x nodes
244 // @param ny Number of y nodes
245 // @param zoom Zoom coefficient
246 // @param nout Pointer to number of output points
247 // @param pout Pointer to array of output points [*nout]
248 //
249 void points_generate1( int nin, point pin[], int nx, int ny, double zoom, int* nout, point** pout )
250 {
251  double xmin = DBL_MAX;
252  double xmax = -DBL_MAX;
253  double ymin = DBL_MAX;
254  double ymax = -DBL_MAX;
255  double stepx, stepy;
256  double x0, xx, yy;
257  int i, j, ii;
258 
259  if ( nx < 1 || ny < 1 )
260  {
261  *pout = NULL;
262  *nout = 0;
263  return;
264  }
265 
266  for ( ii = 0; ii < nin; ++ii )
267  {
268  point* p = &pin[ii];
269 
270  if ( p->x < xmin )
271  xmin = p->x;
272  if ( p->x > xmax )
273  xmax = p->x;
274  if ( p->y < ymin )
275  ymin = p->y;
276  if ( p->y > ymax )
277  ymax = p->y;
278  }
279 
280  if ( isnan( zoom ) || zoom <= 0.0 )
281  zoom = 1.0;
282 
283  if ( zoom != 1.0 )
284  {
285  double xdiff2 = ( xmax - xmin ) / 2.0;
286  double ydiff2 = ( ymax - ymin ) / 2.0;
287  double xav = ( xmax + xmin ) / 2.0;
288  double yav = ( ymax + ymin ) / 2.0;
289 
290  xmin = xav - xdiff2 * zoom;
291  xmax = xav + xdiff2 * zoom;
292  ymin = yav - ydiff2 * zoom;
293  ymax = yav + ydiff2 * zoom;
294  }
295 
296  *nout = nx * ny;
297  *pout = malloc( (size_t) ( *nout ) * sizeof ( point ) );
298 
299  stepx = ( nx > 1 ) ? ( xmax - xmin ) / ( nx - 1 ) : 0.0;
300  stepy = ( ny > 1 ) ? ( ymax - ymin ) / ( ny - 1 ) : 0.0;
301  x0 = ( nx > 1 ) ? xmin : ( xmin + xmax ) / 2.0;
302  yy = ( ny > 1 ) ? ymin : ( ymin + ymax ) / 2.0;
303 
304  ii = 0;
305  for ( j = 0; j < ny; ++j )
306  {
307  xx = x0;
308  for ( i = 0; i < nx; ++i )
309  {
310  point* p = &( *pout )[ii];
311 
312  p->x = xx;
313  p->y = yy;
314  xx += stepx;
315  ii++;
316  }
317  yy += stepy;
318  }
319 }
320 
321 // Generates rectangular grid nx by ny using specified min and max x and y
322 // values. Allocates space for the output point array, be sure to free it
323 // when necessary!
324 //
325 // @param xmin Min x value
326 // @param xmax Max x value
327 // @param ymin Min y value
328 // @param ymax Max y value
329 // @param nx Number of x nodes
330 // @param ny Number of y nodes
331 // @param nout Pointer to number of output points
332 // @param pout Pointer to array of output points [*nout]
333 //
334 void points_generate2( double xmin, double xmax, double ymin, double ymax, int nx, int ny, int* nout, point** pout )
335 {
336  double stepx, stepy;
337  double x0, xx, yy;
338  int i, j, ii;
339 
340  if ( nx < 1 || ny < 1 )
341  {
342  *pout = NULL;
343  *nout = 0;
344  return;
345  }
346 
347  *nout = nx * ny;
348  *pout = malloc( (size_t) ( *nout ) * sizeof ( point ) );
349 
350  stepx = ( nx > 1 ) ? ( xmax - xmin ) / ( nx - 1 ) : 0.0;
351  stepy = ( ny > 1 ) ? ( ymax - ymin ) / ( ny - 1 ) : 0.0;
352  x0 = ( nx > 1 ) ? xmin : ( xmin + xmax ) / 2.0;
353  yy = ( ny > 1 ) ? ymin : ( ymin + ymax ) / 2.0;
354 
355  ii = 0;
356  for ( j = 0; j < ny; ++j )
357  {
358  xx = x0;
359  for ( i = 0; i < nx; ++i )
360  {
361  point* p = &( *pout )[ii];
362 
363  p->x = xx;
364  p->y = yy;
365  xx += stepx;
366  ii++;
367  }
368  yy += stepy;
369  }
370 }
371 
372 static int str2double( char* token, double* value )
373 {
374  char* end = NULL;
375 
376  if ( token == NULL )
377  {
378  *value = NaN;
379  return 0;
380  }
381 
382  *value = strtod( token, &end );
383 
384  if ( end == token )
385  {
386  *value = NaN;
387  return 0;
388  }
389 
390  return 1;
391 }
392 
393 #define NALLOCATED_START 1024
394 
395 // Reads array of points from a columnar file.
396 //
397 // @param fname File name (can be "stdin" for standard input)
398 // @param dim Number of dimensions (must be 2 or 3)
399 // @param n Pointer to number of points (output)
400 // @param points Pointer to array of points [*n] (output) (to be freed)
401 //
402 void points_read( char* fname, int dim, int* n, point** points )
403 {
404  FILE * f = NULL;
405  int nallocated = NALLOCATED_START;
406  char buf[BUFSIZE];
407  char seps[] = " ,;\t";
408  char * token;
409 
410  if ( dim < 2 || dim > 3 )
411  {
412  *n = 0;
413  *points = NULL;
414  return;
415  }
416 
417  if ( fname == NULL )
418  f = stdin;
419  else
420  {
421  if ( strcmp( fname, "stdin" ) == 0 || strcmp( fname, "-" ) == 0 )
422  f = stdin;
423  else
424  {
425  f = fopen( fname, "r" );
426  if ( f == NULL )
427  nn_quit( "%s: %s\n", fname, strerror( errno ) );
428  }
429  }
430 
431  *points = malloc( (size_t) nallocated * sizeof ( point ) );
432  *n = 0;
433  while ( fgets( buf, BUFSIZE, f ) != NULL )
434  {
435  point* p;
436 
437  if ( *n == nallocated )
438  {
439  nallocated *= 2;
440  *points = realloc( *points, (size_t) nallocated * sizeof ( point ) );
441  }
442 
443  p = &( *points )[*n];
444 
445  if ( buf[0] == '#' )
446  continue;
447  if ( ( token = strtok( buf, seps ) ) == NULL )
448  continue;
449  if ( !str2double( token, &p->x ) )
450  continue;
451  if ( ( token = strtok( NULL, seps ) ) == NULL )
452  continue;
453  if ( !str2double( token, &p->y ) )
454  continue;
455  if ( dim == 2 )
456  p->z = NaN;
457  else
458  {
459  if ( ( token = strtok( NULL, seps ) ) == NULL )
460  continue;
461  if ( !str2double( token, &p->z ) )
462  continue;
463  }
464  ( *n )++;
465  }
466 
467  if ( *n == 0 )
468  {
469  free( *points );
470  *points = NULL;
471  }
472  else
473  *points = realloc( *points, (size_t) ( *n ) * sizeof ( point ) );
474 
475  if ( f != stdin )
476  if ( fclose( f ) != 0 )
477  nn_quit( "%s: %s\n", fname, strerror( errno ) );
478 }
479 
480 //* Scales Y coordinate so that the resulting set fits into square:
481 //** xmax - xmin = ymax - ymin
482 //*
483 //* @param n Number of points
484 //* @param points The points to scale
485 //* @return Y axis compression coefficient
486 //
487 double points_scaletosquare( int n, point* points )
488 {
489  double xmin, ymin, xmax, ymax;
490  double k;
491  int i;
492 
493  if ( n <= 0 )
494  return NaN;
495 
496  xmin = xmax = points[0].x;
497  ymin = ymax = points[0].y;
498 
499  for ( i = 1; i < n; ++i )
500  {
501  point* p = &points[i];
502 
503  if ( p->x < xmin )
504  xmin = p->x;
505  else if ( p->x > xmax )
506  xmax = p->x;
507  if ( p->y < ymin )
508  ymin = p->y;
509  else if ( p->y > ymax )
510  ymax = p->y;
511  }
512 
513  if ( xmin == xmax || ymin == ymax )
514  return NaN;
515  else
516  k = ( ymax - ymin ) / ( xmax - xmin );
517 
518  for ( i = 0; i < n; ++i )
519  points[i].y /= k;
520 
521  return k;
522 }
523 
524 //* Compresses Y domain by a given multiple.
525 //
526 // @param n Number of points
527 // @param points The points to scale
528 // @param Y axis compression coefficient as returned by points_scaletosquare()
529 //
530 void points_scale( int n, point* points, double k )
531 {
532  int i;
533 
534  for ( i = 0; i < n; ++i )
535  points[i].y /= k;
536 }