PLplot  5.10.0
nncommon.c
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines