PLplot 5.15.0
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
43int nn_verbose = 0;
46
47#include "version.h"
48
49void nn_quit( const char* format, ... );
50int circle_build( circle* c, point* p1, point* p2, point* p3 );
51int circle_contains( circle* c, point* p );
52
53void 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
68int 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//
114void 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//
249void 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//
334void 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
372static 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//
402void 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//
487double 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//
530void 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}
#define NaN
Definition: csa/nan.h:62
NN_RULE
Definition: nn.h:23
@ SIBSON
Definition: nn.h:23
#define EPSILON
Definition: nncommon.c:41
void points_generate2(double xmin, double xmax, double ymin, double ymax, int nx, int ny, int *nout, point **pout)
Definition: nncommon.c:334
int nn_test_vertice
Definition: nncommon.c:44
int nn_verbose
Definition: nncommon.c:43
void points_scale(int n, point *points, double k)
Definition: nncommon.c:530
NN_RULE nn_rule
Definition: nncommon.c:45
double points_scaletosquare(int n, point *points)
Definition: nncommon.c:487
void points_read(char *fname, int dim, int *n, point **points)
Definition: nncommon.c:402
int circle_build(circle *c, point *p1, point *p2, point *p3)
Definition: nncommon.c:68
int circle_contains(circle *c, point *p)
Definition: nncommon.c:98
void points_thin(int *pn, point **ppoints, int nx, int ny)
Definition: nncommon.c:114
#define NALLOCATED_START
Definition: nncommon.c:393
void points_generate1(int nin, point pin[], int nx, int ny, double zoom, int *nout, point **pout)
Definition: nncommon.c:249
static int str2double(char *token, double *value)
Definition: nncommon.c:372
void nn_quit(const char *format,...)
Definition: nncommon.c:53
#define BUFSIZE
Definition: nncommon.c:39
static PLFLT value(double n1, double n2, double hue)
Definition: plctrl.c:1219
#define isnan(x)
Definition: plplotP.h:262
double x
Definition: delaunay.h:35
double r
Definition: delaunay.h:37
double y
Definition: delaunay.h:36
Definition: csa.h:30
double y
Definition: csa.h:32
double x
Definition: csa.h:31
double z
Definition: csa.h:33
static char buf[200]
Definition: tclAPI.c:873