PLplot
5.10.0
|
00001 //-------------------------------------------------------------------------- 00002 // 00003 // File: nncommon.c 00004 // 00005 // Created: 04/08/2000 00006 // 00007 // Author: Pavel Sakov 00008 // CSIRO Marine Research 00009 // 00010 // Purpose: Common stuff for NN interpolation library 00011 // 00012 // Description: None 00013 // 00014 // Revisions: 15/11/2002 PS: Changed name from "utils.c" 00015 // 28/02/2003 PS: Modified points_read() to do the job without 00016 // rewinding the file. This allows to read from stdin when 00017 // necessary. 00018 // 09/04/2003 PS: Modified points_read() to read from a 00019 // file specified by name, not by handle. 00020 // Modified: Andrew Ross 20/10/2008 00021 // Change <= comparison in circle_contains() to use EPSILON 00022 // to catch case where the point lies on the circle and there 00023 // is floating point rounding error in the radii. 00024 // 00025 //-------------------------------------------------------------------------- 00026 00027 #include <stdlib.h> 00028 #include <stdio.h> 00029 #include <stdarg.h> 00030 #include <assert.h> 00031 #include <math.h> 00032 #include <limits.h> 00033 #include <float.h> 00034 #include <string.h> 00035 #include <errno.h> 00036 #include "nan.h" 00037 #include "delaunay.h" 00038 00039 #define BUFSIZE 1024 00040 00041 #define EPSILON 1.0e-8 00042 00043 int nn_verbose = 0; 00044 int nn_test_vertice = -1; 00045 NN_RULE nn_rule = SIBSON; 00046 00047 #include "version.h" 00048 00049 void nn_quit( const char* format, ... ); 00050 int circle_build( circle* c, point* p1, point* p2, point* p3 ); 00051 int circle_contains( circle* c, point* p ); 00052 00053 void nn_quit( const char* format, ... ) 00054 { 00055 va_list args; 00056 00057 fflush( stdout ); // just in case, to have the exit message 00058 // last 00059 00060 fprintf( stderr, "error: nn: " ); 00061 va_start( args, format ); 00062 vfprintf( stderr, format, args ); 00063 va_end( args ); 00064 00065 exit( 1 ); 00066 } 00067 00068 int circle_build( circle* c, point* p1, point* p2, point* p3 ) 00069 { 00070 double x1sq = p1->x * p1->x; 00071 double x2sq = p2->x * p2->x; 00072 double x3sq = p3->x * p3->x; 00073 double y1sq = p1->y * p1->y; 00074 double y2sq = p2->y * p2->y; 00075 double y3sq = p3->y * p3->y; 00076 double t1 = x3sq - x2sq + y3sq - y2sq; 00077 double t2 = x1sq - x3sq + y1sq - y3sq; 00078 double t3 = x2sq - x1sq + y2sq - y1sq; 00079 double D = ( p1->x * ( p2->y - p3->y ) + p2->x * ( p3->y - p1->y ) + p3->x * ( p1->y - p2->y ) ) * 2.0; 00080 00081 if ( D == 0.0 ) 00082 return 0; 00083 00084 c->x = ( p1->y * t1 + p2->y * t2 + p3->y * t3 ) / D; 00085 c->y = -( p1->x * t1 + p2->x * t2 + p3->x * t3 ) / D; 00086 c->r = hypot( c->x - p1->x, c->y - p1->y ); 00087 00088 return 1; 00089 } 00090 00091 // This procedure has taken it final shape after a number of tries. The problem 00092 // was to have the calculated and stored radii being the same if (x,y) is 00093 // exactly on the circle border (i.e. not to use FCPU extended precision in 00094 // the radius calculation). This may have little effect in practice but was 00095 // important in some tests when both input and output data were placed 00096 // in rectangular grid nodes. 00097 // 00098 int circle_contains( circle* c, point* p ) 00099 { 00100 return hypot( c->x - p->x, c->y - p->y ) <= c->r * ( 1.0 + EPSILON ); 00101 } 00102 00103 // Smoothes the input point array by averaging the input x,y and z values 00104 // for each cell within virtual rectangular nx by ny grid. The corners of the 00105 // grid are created from min and max values of the input array. It also frees 00106 // the original array and returns results and new dimension via original 00107 // data and size pointers. 00108 // 00109 // @param pn Pointer to number of points (input/output) 00110 // @param ppoints Pointer to array of points (input/output) [*pn] 00111 // @param nx Number of x nodes in decimation 00112 // @param ny Number of y nodes in decimation 00113 // 00114 void points_thin( int* pn, point** ppoints, int nx, int ny ) 00115 { 00116 int n = *pn; 00117 point * points = *ppoints; 00118 double xmin = DBL_MAX; 00119 double xmax = -DBL_MAX; 00120 double ymin = DBL_MAX; 00121 double ymax = -DBL_MAX; 00122 int nxy = nx * ny; 00123 double * sumx = calloc( (size_t) nxy, sizeof ( double ) ); 00124 double * sumy = calloc( (size_t) nxy, sizeof ( double ) ); 00125 double * sumz = calloc( (size_t) nxy, sizeof ( double ) ); 00126 int * count = calloc( (size_t) nxy, sizeof ( int ) ); 00127 double stepx = 0.0; 00128 double stepy = 0.0; 00129 int nnew = 0; 00130 point * pointsnew = NULL; 00131 int i, j, ii; 00132 00133 if ( nn_verbose ) 00134 fprintf( stderr, "thinned: %d points -> ", *pn ); 00135 00136 if ( nx < 1 || ny < 1 ) 00137 { 00138 free( points ); 00139 *ppoints = NULL; 00140 *pn = 0; 00141 if ( nn_verbose ) 00142 fprintf( stderr, "0 points" ); 00143 free( sumx ); 00144 free( sumy ); 00145 free( sumz ); 00146 free( count ); 00147 return; 00148 } 00149 00150 for ( ii = 0; ii < n; ++ii ) 00151 { 00152 point* p = &points[ii]; 00153 00154 if ( p->x < xmin ) 00155 xmin = p->x; 00156 if ( p->x > xmax ) 00157 xmax = p->x; 00158 if ( p->y < ymin ) 00159 ymin = p->y; 00160 if ( p->y > ymax ) 00161 ymax = p->y; 00162 } 00163 00164 stepx = ( nx > 1 ) ? ( xmax - xmin ) / nx : 0.0; 00165 stepy = ( ny > 1 ) ? ( ymax - ymin ) / ny : 0.0; 00166 00167 for ( ii = 0; ii < n; ++ii ) 00168 { 00169 point* p = &points[ii]; 00170 int index; 00171 00172 // 00173 // Following is the portion of the code which really depends on the 00174 // floating point particulars. Do not be surprised if different 00175 // compilers/options give different results here. 00176 // 00177 i = ( nx == 1 ) ? 0 : (int) ( ( p->x - xmin ) / stepx ); 00178 j = ( ny == 1 ) ? 0 : (int) ( ( p->y - ymin ) / stepy ); 00179 00180 if ( i == nx ) 00181 i--; 00182 if ( j == ny ) 00183 j--; 00184 index = i + j * nx; 00185 sumx[index] += p->x; 00186 sumy[index] += p->y; 00187 sumz[index] += p->z; 00188 count[index]++; 00189 } 00190 00191 for ( j = 0; j < ny; ++j ) 00192 { 00193 for ( i = 0; i < nx; ++i ) 00194 { 00195 int index = i + j * nx; 00196 00197 if ( count[index] > 0 ) 00198 nnew++; 00199 } 00200 } 00201 00202 pointsnew = malloc( (size_t) nnew * sizeof ( point ) ); 00203 00204 ii = 0; 00205 for ( j = 0; j < ny; ++j ) 00206 { 00207 for ( i = 0; i < nx; ++i ) 00208 { 00209 int index = i + j * nx; 00210 int nn = count[index]; 00211 00212 if ( nn > 0 ) 00213 { 00214 point* p = &pointsnew[ii]; 00215 00216 p->x = sumx[index] / nn; 00217 p->y = sumy[index] / nn; 00218 p->z = sumz[index] / nn; 00219 ii++; 00220 } 00221 } 00222 } 00223 00224 if ( nn_verbose ) 00225 fprintf( stderr, "%d points\n", nnew ); 00226 00227 free( sumx ); 00228 free( sumy ); 00229 free( sumz ); 00230 free( count ); 00231 00232 free( points ); 00233 *ppoints = pointsnew; 00234 *pn = nnew; 00235 } 00236 00237 // Generates rectangular grid nx by ny using min and max x and y values from 00238 // the input point array. Allocates space for the output point array, be sure 00239 // to free it when necessary! 00240 // 00241 // @param n Number of points 00242 // @param points Array of points [n] 00243 // @param nx Number of x nodes 00244 // @param ny Number of y nodes 00245 // @param zoom Zoom coefficient 00246 // @param nout Pointer to number of output points 00247 // @param pout Pointer to array of output points [*nout] 00248 // 00249 void points_generate1( int nin, point pin[], int nx, int ny, double zoom, int* nout, point** pout ) 00250 { 00251 double xmin = DBL_MAX; 00252 double xmax = -DBL_MAX; 00253 double ymin = DBL_MAX; 00254 double ymax = -DBL_MAX; 00255 double stepx, stepy; 00256 double x0, xx, yy; 00257 int i, j, ii; 00258 00259 if ( nx < 1 || ny < 1 ) 00260 { 00261 *pout = NULL; 00262 *nout = 0; 00263 return; 00264 } 00265 00266 for ( ii = 0; ii < nin; ++ii ) 00267 { 00268 point* p = &pin[ii]; 00269 00270 if ( p->x < xmin ) 00271 xmin = p->x; 00272 if ( p->x > xmax ) 00273 xmax = p->x; 00274 if ( p->y < ymin ) 00275 ymin = p->y; 00276 if ( p->y > ymax ) 00277 ymax = p->y; 00278 } 00279 00280 if ( isnan( zoom ) || zoom <= 0.0 ) 00281 zoom = 1.0; 00282 00283 if ( zoom != 1.0 ) 00284 { 00285 double xdiff2 = ( xmax - xmin ) / 2.0; 00286 double ydiff2 = ( ymax - ymin ) / 2.0; 00287 double xav = ( xmax + xmin ) / 2.0; 00288 double yav = ( ymax + ymin ) / 2.0; 00289 00290 xmin = xav - xdiff2 * zoom; 00291 xmax = xav + xdiff2 * zoom; 00292 ymin = yav - ydiff2 * zoom; 00293 ymax = yav + ydiff2 * zoom; 00294 } 00295 00296 *nout = nx * ny; 00297 *pout = malloc( (size_t) ( *nout ) * sizeof ( point ) ); 00298 00299 stepx = ( nx > 1 ) ? ( xmax - xmin ) / ( nx - 1 ) : 0.0; 00300 stepy = ( ny > 1 ) ? ( ymax - ymin ) / ( ny - 1 ) : 0.0; 00301 x0 = ( nx > 1 ) ? xmin : ( xmin + xmax ) / 2.0; 00302 yy = ( ny > 1 ) ? ymin : ( ymin + ymax ) / 2.0; 00303 00304 ii = 0; 00305 for ( j = 0; j < ny; ++j ) 00306 { 00307 xx = x0; 00308 for ( i = 0; i < nx; ++i ) 00309 { 00310 point* p = &( *pout )[ii]; 00311 00312 p->x = xx; 00313 p->y = yy; 00314 xx += stepx; 00315 ii++; 00316 } 00317 yy += stepy; 00318 } 00319 } 00320 00321 // Generates rectangular grid nx by ny using specified min and max x and y 00322 // values. Allocates space for the output point array, be sure to free it 00323 // when necessary! 00324 // 00325 // @param xmin Min x value 00326 // @param xmax Max x value 00327 // @param ymin Min y value 00328 // @param ymax Max y value 00329 // @param nx Number of x nodes 00330 // @param ny Number of y nodes 00331 // @param nout Pointer to number of output points 00332 // @param pout Pointer to array of output points [*nout] 00333 // 00334 void points_generate2( double xmin, double xmax, double ymin, double ymax, int nx, int ny, int* nout, point** pout ) 00335 { 00336 double stepx, stepy; 00337 double x0, xx, yy; 00338 int i, j, ii; 00339 00340 if ( nx < 1 || ny < 1 ) 00341 { 00342 *pout = NULL; 00343 *nout = 0; 00344 return; 00345 } 00346 00347 *nout = nx * ny; 00348 *pout = malloc( (size_t) ( *nout ) * sizeof ( point ) ); 00349 00350 stepx = ( nx > 1 ) ? ( xmax - xmin ) / ( nx - 1 ) : 0.0; 00351 stepy = ( ny > 1 ) ? ( ymax - ymin ) / ( ny - 1 ) : 0.0; 00352 x0 = ( nx > 1 ) ? xmin : ( xmin + xmax ) / 2.0; 00353 yy = ( ny > 1 ) ? ymin : ( ymin + ymax ) / 2.0; 00354 00355 ii = 0; 00356 for ( j = 0; j < ny; ++j ) 00357 { 00358 xx = x0; 00359 for ( i = 0; i < nx; ++i ) 00360 { 00361 point* p = &( *pout )[ii]; 00362 00363 p->x = xx; 00364 p->y = yy; 00365 xx += stepx; 00366 ii++; 00367 } 00368 yy += stepy; 00369 } 00370 } 00371 00372 static int str2double( char* token, double* value ) 00373 { 00374 char* end = NULL; 00375 00376 if ( token == NULL ) 00377 { 00378 *value = NaN; 00379 return 0; 00380 } 00381 00382 *value = strtod( token, &end ); 00383 00384 if ( end == token ) 00385 { 00386 *value = NaN; 00387 return 0; 00388 } 00389 00390 return 1; 00391 } 00392 00393 #define NALLOCATED_START 1024 00394 00395 // Reads array of points from a columnar file. 00396 // 00397 // @param fname File name (can be "stdin" for standard input) 00398 // @param dim Number of dimensions (must be 2 or 3) 00399 // @param n Pointer to number of points (output) 00400 // @param points Pointer to array of points [*n] (output) (to be freed) 00401 // 00402 void points_read( char* fname, int dim, int* n, point** points ) 00403 { 00404 FILE * f = NULL; 00405 int nallocated = NALLOCATED_START; 00406 char buf[BUFSIZE]; 00407 char seps[] = " ,;\t"; 00408 char * token; 00409 00410 if ( dim < 2 || dim > 3 ) 00411 { 00412 *n = 0; 00413 *points = NULL; 00414 return; 00415 } 00416 00417 if ( fname == NULL ) 00418 f = stdin; 00419 else 00420 { 00421 if ( strcmp( fname, "stdin" ) == 0 || strcmp( fname, "-" ) == 0 ) 00422 f = stdin; 00423 else 00424 { 00425 f = fopen( fname, "r" ); 00426 if ( f == NULL ) 00427 nn_quit( "%s: %s\n", fname, strerror( errno ) ); 00428 } 00429 } 00430 00431 *points = malloc( (size_t) nallocated * sizeof ( point ) ); 00432 *n = 0; 00433 while ( fgets( buf, BUFSIZE, f ) != NULL ) 00434 { 00435 point* p; 00436 00437 if ( *n == nallocated ) 00438 { 00439 nallocated *= 2; 00440 *points = realloc( *points, (size_t) nallocated * sizeof ( point ) ); 00441 } 00442 00443 p = &( *points )[*n]; 00444 00445 if ( buf[0] == '#' ) 00446 continue; 00447 if ( ( token = strtok( buf, seps ) ) == NULL ) 00448 continue; 00449 if ( !str2double( token, &p->x ) ) 00450 continue; 00451 if ( ( token = strtok( NULL, seps ) ) == NULL ) 00452 continue; 00453 if ( !str2double( token, &p->y ) ) 00454 continue; 00455 if ( dim == 2 ) 00456 p->z = NaN; 00457 else 00458 { 00459 if ( ( token = strtok( NULL, seps ) ) == NULL ) 00460 continue; 00461 if ( !str2double( token, &p->z ) ) 00462 continue; 00463 } 00464 ( *n )++; 00465 } 00466 00467 if ( *n == 0 ) 00468 { 00469 free( *points ); 00470 *points = NULL; 00471 } 00472 else 00473 *points = realloc( *points, (size_t) ( *n ) * sizeof ( point ) ); 00474 00475 if ( f != stdin ) 00476 if ( fclose( f ) != 0 ) 00477 nn_quit( "%s: %s\n", fname, strerror( errno ) ); 00478 } 00479 00480 //* Scales Y coordinate so that the resulting set fits into square: 00481 //** xmax - xmin = ymax - ymin 00482 //* 00483 //* @param n Number of points 00484 //* @param points The points to scale 00485 //* @return Y axis compression coefficient 00486 // 00487 double points_scaletosquare( int n, point* points ) 00488 { 00489 double xmin, ymin, xmax, ymax; 00490 double k; 00491 int i; 00492 00493 if ( n <= 0 ) 00494 return NaN; 00495 00496 xmin = xmax = points[0].x; 00497 ymin = ymax = points[0].y; 00498 00499 for ( i = 1; i < n; ++i ) 00500 { 00501 point* p = &points[i]; 00502 00503 if ( p->x < xmin ) 00504 xmin = p->x; 00505 else if ( p->x > xmax ) 00506 xmax = p->x; 00507 if ( p->y < ymin ) 00508 ymin = p->y; 00509 else if ( p->y > ymax ) 00510 ymax = p->y; 00511 } 00512 00513 if ( xmin == xmax || ymin == ymax ) 00514 return NaN; 00515 else 00516 k = ( ymax - ymin ) / ( xmax - xmin ); 00517 00518 for ( i = 0; i < n; ++i ) 00519 points[i].y /= k; 00520 00521 return k; 00522 } 00523 00524 //* Compresses Y domain by a given multiple. 00525 // 00526 // @param n Number of points 00527 // @param points The points to scale 00528 // @param Y axis compression coefficient as returned by points_scaletosquare() 00529 // 00530 void points_scale( int n, point* points, double k ) 00531 { 00532 int i; 00533 00534 for ( i = 0; i < n; ++i ) 00535 points[i].y /= k; 00536 }