PLplot  5.10.0
csa.c
Go to the documentation of this file.
00001 //--------------------------------------------------------------------------
00002 //
00003 // File:           csa.c
00004 //
00005 // Created:        16/10/2002
00006 //
00007 // Author:         Pavel Sakov
00008 //                 CSIRO Marine Research
00009 //
00010 // Purpose:        2D data approximation with bivariate C1 cubic spline.
00011 //                 A set of library functions + standalone utility.
00012 //
00013 // Description:    See J. Haber, F. Zeilfelder, O.Davydov and H.-P. Seidel,
00014 //                 Smooth approximation and rendering of large scattered data
00015 //                 sets, in  ``Proceedings of IEEE Visualization 2001''
00016 //                 (Th.Ertl, K.Joy and A.Varshney, Eds.), pp.341-347, 571,
00017 //                 IEEE Computer Society, 2001.
00018 //                 http://www.uni-giessen.de/www-Numerische-Mathematik/
00019 //                        davydov/VIS2001.ps.gz
00020 //                 http://www.math.uni-mannheim.de/~lsmath4/paper/
00021 //                        VIS2001.pdf.gz
00022 //
00023 // Revisions:      09/04/2003 PS: Modified points_read() to read from a
00024 //                   file specified by name, not by handle.
00025 //
00026 //--------------------------------------------------------------------------
00027 
00028 #include <stdlib.h>
00029 #include <stdio.h>
00030 #include <stdarg.h>
00031 #include <limits.h>
00032 #include <float.h>
00033 #include <math.h>
00034 #include <assert.h>
00035 #include <string.h>
00036 #include <errno.h>
00037 #include "version.h"
00038 #include "nan.h"
00039 #include "csa.h"
00040 
00041 int csa_verbose = 0;
00042 
00043 #define NPASTART     5          // Number of Points Allocated at Start
00044 #define SVD_NMAX     30         // Maximal number of iterations in svd()
00045 
00046 // default algorithm parameters
00047 #define NPMIN_DEF    3
00048 #define NPMAX_DEF    40
00049 #define K_DEF        140
00050 #define NPPC_DEF     5
00051 
00052 struct square;
00053 typedef struct square   square;
00054 
00055 typedef struct
00056 {
00057     square * parent;
00058     int    index;               // index within parent square; 0 <= index <=
00059                                 // 3
00060     point  vertices[3];
00061     point  middle;              // barycenter
00062     double h;                   // parent square edge length
00063     double r;                   // data visibility radius
00064 
00065     //
00066     // points used -- in primary triangles only
00067     //
00068     int  nallocated;
00069     int  npoints;
00070     point** points;
00071 
00072     int  primary;               // flag -- whether calculate spline
00073                                 // coefficients directly (by least squares
00074                                 // method) (primary = 1) or indirectly (from
00075                                 // C1 smoothness conditions) (primary = 0)
00076     int hascoeffs;              // flag -- whether there are no NaNs among
00077                                 // the spline coefficients
00078     int order;                  // spline order -- for primary triangles only
00079                                 //
00080 } triangle;
00081 
00082 struct square
00083 {
00084     csa     * parent;
00085     int     i, j;               // indices
00086 
00087     int     nallocated;
00088     int     npoints;
00089     point   ** points;
00090 
00091     int     primary;            // flag -- whether this square contains a
00092                                 // primary triangle
00093     triangle* triangles[4];
00094 
00095     double  coeffs[25];
00096 };
00097 
00098 struct csa
00099 {
00100     double xmin;
00101     double xmax;
00102     double ymin;
00103     double ymax;
00104 
00105     int    nallocated;
00106     int    npoints;
00107     point  ** points;
00108 
00109     //
00110     // squarization
00111     //
00112     int     ni;
00113     int     nj;
00114     double  h;
00115     square  *** squares;        // square* [j][i]
00116 
00117     int     npt;                // Number of Primary Triangles
00118     triangle** pt;              // Primary Triangles -- triangle* [npt]
00119 
00120     //
00121     // algorithm parameters
00122     //
00123     int    npmin;               // minimal number of points locally involved
00124                                 // in * spline calculation (normally = 3)
00125     int    npmax;               // maximal number of points locally involved
00126                                 // in * spline calculation (required > 10, *
00127                                 // recommended 20 < npmax < 60)
00128     double k;                   // relative tolerance multiple in fitting
00129                                 // spline coefficients: the higher this
00130                                 // value, the higher degree of the locally
00131                                 // fitted spline (recommended 80 < k < 200)
00132     int nppc;                   // average number of points per square
00133 };
00134 
00135 void csa_setnppc( csa* a, double nppc );
00136 
00137 static void csa_quit( const char* format, ... )
00138 {
00139     va_list args;
00140 
00141     fflush( stdout );             // just in case -- to have the exit message
00142                                   // last
00143 
00144     fprintf( stderr, "error: csa: " );
00145     va_start( args, format );
00146     vfprintf( stderr, format, args );
00147     va_end( args );
00148 
00149     exit( 1 );
00150 }
00151 
00152 // Allocates n1xn2 matrix of something. Note that it will be accessed as
00153 // [n2][n1].
00154 // @param n1 Number of columns
00155 // @param n2 Number of rows
00156 // @return Matrix
00157 //
00158 static void* alloc2d( int n1, int n2, size_t unitsize )
00159 {
00160     size_t size;
00161     char   * p;
00162     char   ** pp;
00163     int    i;
00164 
00165     assert( n1 > 0 );
00166     assert( n2 > 0 );
00167     assert( (double) n1 * (double) n2 <= (double) UINT_MAX );
00168 
00169     size = (size_t) ( n1 * n2 );
00170     if ( ( p = calloc( size, unitsize ) ) == NULL )
00171         csa_quit( "alloc2d(): %s\n", strerror( errno ) );
00172 
00173     assert( (double) n2 * (double) sizeof ( void* ) <= (double) UINT_MAX );
00174 
00175     size = (size_t) n2 * sizeof ( void* );
00176     if ( ( pp = malloc( size ) ) == NULL )
00177         csa_quit( "alloc2d(): %s\n", strerror( errno ) );
00178     for ( i = 0; i < n2; i++ )
00179         pp[i] = &p[(size_t) ( i * n1 ) * unitsize];
00180 
00181     return pp;
00182 }
00183 
00184 // Destroys n1xn2 matrix.
00185 // @param pp Matrix
00186 //
00187 static void free2d( void* pp )
00188 {
00189     void* p;
00190 
00191     assert( pp != NULL );
00192     p = ( (void **) pp )[0];
00193     free( pp );
00194     assert( p != NULL );
00195     free( p );
00196 }
00197 
00198 static triangle* triangle_create( square* s, point vertices[], int index )
00199 {
00200     triangle* t = malloc( sizeof ( triangle ) );
00201 
00202     t->parent = s;
00203     memcpy( t->vertices, vertices, sizeof ( point ) * 3 );
00204     t->middle.x = ( vertices[0].x + vertices[1].x + vertices[2].x ) / 3.0;
00205     t->middle.y = ( vertices[0].y + vertices[1].y + vertices[2].y ) / 3.0;
00206     t->h        = s->parent->h;
00207     t->index    = index;
00208 
00209     t->r          = 0.0;
00210     t->points     = 0;
00211     t->nallocated = 0;
00212     t->npoints    = 0;
00213     t->hascoeffs  = 0;
00214     t->primary    = 0;
00215     t->order      = -1;
00216 
00217     return t;
00218 }
00219 
00220 static void triangle_addpoint( triangle* t, point* p )
00221 {
00222     if ( t->nallocated == t->npoints )
00223     {
00224         if ( t->nallocated == 0 )
00225         {
00226             t->points     = malloc( NPASTART * sizeof ( point* ) );
00227             t->nallocated = NPASTART;
00228         }
00229         else
00230         {
00231             t->points      = realloc( t->points, (size_t) ( t->nallocated * 2 ) * sizeof ( point* ) );
00232             t->nallocated *= 2;
00233         }
00234     }
00235 
00236     t->points[t->npoints] = p;
00237     t->npoints++;
00238 }
00239 
00240 static void triangle_destroy( triangle* t )
00241 {
00242     if ( t->points != NULL )
00243         free( t->points );
00244     free( t );
00245 }
00246 
00247 // Calculates barycentric coordinates of a point.
00248 // Takes into account that possible triangles are rectangular, with the right
00249 // angle at t->vertices[0], the vertices[1] vertex being in
00250 // (-3*PI/4) + (PI/2) * t->index direction from vertices[0], and
00251 // vertices[2] being at (-5*PI/4) + (PI/2) * t->index.
00252 //
00253 static void triangle_calculatebc( triangle* t, point* p, double bc[] )
00254 {
00255     double dx = p->x - t->vertices[0].x;
00256     double dy = p->y - t->vertices[0].y;
00257 
00258     if ( t->index == 0 )
00259     {
00260         bc[1] = ( dy - dx ) / t->h;
00261         bc[2] = -( dx + dy ) / t->h;
00262     }
00263     else if ( t->index == 1 )
00264     {
00265         bc[1] = ( dx + dy ) / t->h;
00266         bc[2] = ( dy - dx ) / t->h;
00267     }
00268     else if ( t->index == 2 )
00269     {
00270         bc[1] = ( dx - dy ) / t->h;
00271         bc[2] = ( dx + dy ) / t->h;
00272     }
00273     else
00274     {
00275         bc[1] = -( dx + dy ) / t->h;
00276         bc[2] = ( dx - dy ) / t->h;
00277     }
00278     bc[0] = 1.0 - bc[1] - bc[2];
00279 }
00280 
00281 static square* square_create( csa* parent, double xmin, double ymin, int i, int j )
00282 {
00283     int    ii;
00284 
00285     square * s = malloc( sizeof ( square ) );
00286     double h   = parent->h;
00287 
00288     s->parent = parent;
00289     s->i      = i;
00290     s->j      = j;
00291 
00292     s->points     = NULL;
00293     s->nallocated = 0;
00294     s->npoints    = 0;
00295 
00296     s->primary = 0;
00297 
00298     //
00299     // create 4 triangles formed by diagonals
00300     //
00301     for ( ii = 0; ii < 4; ++ii )
00302     {
00303         point vertices[3];
00304 
00305         vertices[0].x    = xmin + h / 2.0;
00306         vertices[0].y    = ymin + h / 2.0;
00307         vertices[1].x    = xmin + h * ( ( ( ii + 1 ) % 4 ) / 2 ); // 0 1 1 0
00308         vertices[1].y    = ymin + h * ( ( ( ii + 2 ) % 4 ) / 2 ); // 1 1 0 0
00309         vertices[2].x    = xmin + h * ( ii / 2 );                 // 0 0 1 1
00310         vertices[2].y    = ymin + h * ( ( ( ii + 1 ) % 4 ) / 2 ); // 0 1 1 0
00311         s->triangles[ii] = triangle_create( s, vertices, ii );
00312     }
00313 
00314     for ( ii = 0; ii < 25; ++ii )
00315         s->coeffs[ii] = NaN;
00316 
00317     return s;
00318 }
00319 
00320 static void square_destroy( square* s )
00321 {
00322     int i;
00323 
00324     for ( i = 0; i < 4; ++i )
00325         triangle_destroy( s->triangles[i] );
00326     if ( s->points != NULL )
00327         free( s->points );
00328     free( s );
00329 }
00330 
00331 static void square_addpoint( square* s, point* p )
00332 {
00333     if ( s->nallocated == s->npoints )
00334     {
00335         if ( s->nallocated == 0 )
00336         {
00337             s->points     = malloc( NPASTART * sizeof ( point* ) );
00338             s->nallocated = NPASTART;
00339         }
00340         else
00341         {
00342             s->points      = realloc( s->points, (size_t) ( s->nallocated * 2 ) * sizeof ( point* ) );
00343             s->nallocated *= 2;
00344         }
00345     }
00346 
00347     s->points[s->npoints] = p;
00348     s->npoints++;
00349 }
00350 
00351 csa* csa_create()
00352 {
00353     csa* a = malloc( sizeof ( csa ) );
00354 
00355     a->xmin = DBL_MAX;
00356     a->xmax = -DBL_MAX;
00357     a->ymin = DBL_MAX;
00358     a->ymax = -DBL_MAX;
00359 
00360     a->points     = malloc( NPASTART * sizeof ( point* ) );
00361     a->nallocated = NPASTART;
00362     a->npoints    = 0;
00363 
00364     a->ni      = 0;
00365     a->nj      = 0;
00366     a->h       = NaN;
00367     a->squares = NULL;
00368 
00369     a->npt = 0;
00370     a->pt  = NULL;
00371 
00372     a->npmin = NPMIN_DEF;
00373     a->npmax = NPMAX_DEF;
00374     a->k     = K_DEF;
00375     a->nppc  = NPPC_DEF;
00376 
00377     return a;
00378 }
00379 
00380 void csa_destroy( csa* a )
00381 {
00382     int i, j;
00383 
00384     if ( a->squares != NULL )
00385     {
00386         for ( j = 0; j < a->nj; ++j )
00387             for ( i = 0; i < a->ni; ++i )
00388                 square_destroy( a->squares[j][i] );
00389         free2d( a->squares );
00390     }
00391     if ( a->pt != NULL )
00392         free( a->pt );
00393     if ( a->points != NULL )
00394         free( a->points );
00395     free( a );
00396 }
00397 
00398 void csa_addpoints( csa* a, int n, point points[] )
00399 {
00400     int na = a->nallocated;
00401     int i;
00402 
00403     assert( a->squares == NULL );
00404     //
00405     // (can be called prior to squarization only)
00406     //
00407 
00408     while ( na < a->npoints + n )
00409         na *= 2;
00410     if ( na != a->nallocated )
00411     {
00412         a->points     = realloc( a->points, (size_t) na * sizeof ( point* ) );
00413         a->nallocated = na;
00414     }
00415 
00416     for ( i = 0; i < n; ++i )
00417     {
00418         point* p = &points[i];
00419 
00420         a->points[a->npoints] = p;
00421         a->npoints++;
00422 
00423         if ( p->x < a->xmin )
00424             a->xmin = p->x;
00425         if ( p->x > a->xmax )
00426             a->xmax = p->x;
00427         if ( p->y < a->ymin )
00428             a->ymin = p->y;
00429         if ( p->y > a->ymax )
00430             a->ymax = p->y;
00431     }
00432 }
00433 
00434 // Marks the squares containing "primary" triangles by setting "primary" flag
00435 // to 1.
00436 //
00437 static void csa_setprimaryflag( csa* a )
00438 {
00439     square*** squares = a->squares;
00440     int   nj1         = a->nj - 1;
00441     int   ni1         = a->ni - 1;
00442     int   i, j;
00443 
00444     for ( j = 1; j < nj1; ++j )
00445     {
00446         for ( i = 1; i < ni1; ++i )
00447         {
00448             if ( squares[j][i]->npoints > 0 )
00449             {
00450                 if ( ( i + j ) % 2 == 0 )
00451                 {
00452                     squares[j][i]->primary         = 1;
00453                     squares[j - 1][i - 1]->primary = 1;
00454                     squares[j + 1][i - 1]->primary = 1;
00455                     squares[j - 1][i + 1]->primary = 1;
00456                     squares[j + 1][i + 1]->primary = 1;
00457                 }
00458                 else
00459                 {
00460                     squares[j - 1][i]->primary = 1;
00461                     squares[j + 1][i]->primary = 1;
00462                     squares[j][i - 1]->primary = 1;
00463                     squares[j][i + 1]->primary = 1;
00464                 }
00465             }
00466         }
00467     }
00468 }
00469 
00470 // Splits the data domain in a number of squares.
00471 //
00472 static void csa_squarize( csa* a )
00473 {
00474     int    nps[7] = { 0, 0, 0, 0, 0, 0 };  // stats on number of points per
00475                                            // square
00476     double dx      = a->xmax - a->xmin;
00477     double dy      = a->ymax - a->ymin;
00478     int    npoints = a->npoints;
00479     double h;
00480     int    i, j, ii, nadj;
00481 
00482     if ( csa_verbose )
00483     {
00484         fprintf( stderr, "squarizing csa:\n" );
00485         fflush( stderr );
00486     }
00487 
00488     assert( a->squares == NULL );
00489     //
00490     // (can be done only once)
00491     //
00492 
00493     h = sqrt( dx * dy * a->nppc / npoints );      // square edge size
00494     if ( dx < h )
00495         h = dy * a->nppc / npoints;
00496     if ( dy < h )
00497         h = dx * a->nppc / npoints;
00498     a->h = h;
00499 
00500     a->ni = (int) ( ceil( dx / h ) + 2 );
00501     a->nj = (int) ( ceil( dy / h ) + 2 );
00502 
00503     if ( csa_verbose )
00504     {
00505         fprintf( stderr, "  %d x %d squares\n", a->ni, a->nj );
00506         fflush( stderr );
00507     }
00508 
00509     //
00510     // create squares
00511     //
00512     a->squares = alloc2d( a->ni, a->nj, sizeof ( void* ) );
00513     for ( j = 0; j < a->nj; ++j )
00514         for ( i = 0; i < a->ni; ++i )
00515             a->squares[j][i] = square_create( a, a->xmin + h * ( i - 1 ), a->ymin + h * ( j - 1 ), i, j );
00516 
00517     //
00518     // map points to squares
00519     //
00520     for ( ii = 0; ii < npoints; ++ii )
00521     {
00522         point* p = a->points[ii];
00523 
00524         i = (int) ( floor( ( p->x - a->xmin ) / h ) + 1 );
00525         j = (int) ( floor( ( p->y - a->ymin ) / h ) + 1 );
00526         square_addpoint( a->squares[j][i], p );
00527     }
00528 
00529     //
00530     // mark relevant squares with no points
00531     //
00532     csa_setprimaryflag( a );
00533 
00534     //
00535     // Create a list of "primary" triangles, for which spline coefficients
00536     // will be calculated directy (by least squares method), without using
00537     // C1 smoothness conditions.
00538     //
00539     a->pt = malloc( (size_t) ( ( a->ni / 2 + 1 ) * a->nj ) * sizeof ( triangle* ) );
00540     for ( j = 0, ii = 0, nadj = 0; j < a->nj; ++j )
00541     {
00542         for ( i = 0; i < a->ni; ++i )
00543         {
00544             square* s = a->squares[j][i];
00545 
00546             if ( s->npoints > 0 )
00547             {
00548                 int nn = s->npoints / 5;
00549 
00550                 if ( nn > 6 )
00551                     nn = 6;
00552                 nps[nn]++;
00553                 ii++;
00554             }
00555             if ( s->primary && s->npoints == 0 )
00556                 nadj++;
00557             if ( s->primary )
00558             {
00559                 a->pt[a->npt]            = s->triangles[0];
00560                 s->triangles[0]->primary = 1;
00561                 a->npt++;
00562             }
00563         }
00564     }
00565 
00566     if ( csa_verbose )
00567     {
00568         fprintf( stderr, "  %d non-empty squares\n", ii );
00569         fprintf( stderr, "  %d primary squares\n", a->npt );
00570         fprintf( stderr, "  %d primary squares with no data\n", nadj );
00571         fprintf( stderr, "  %.2f points per square \n", (double) a->npoints / ii );
00572     }
00573 
00574     if ( csa_verbose == 2 )
00575     {
00576         for ( i = 0; i < 6; ++i )
00577             fprintf( stderr, "  %d-%d points -- %d squares\n", i * 5, i * 5 + 4, nps[i] );
00578         fprintf( stderr, "  %d or more points -- %d squares\n", i * 5, nps[i] );
00579     }
00580 
00581     if ( csa_verbose == 2 )
00582     {
00583         fprintf( stderr, " j\\i" );
00584         for ( i = 0; i < a->ni; ++i )
00585             fprintf( stderr, "%3d ", i );
00586         fprintf( stderr, "\n" );
00587         for ( j = a->nj - 1; j >= 0; --j )
00588         {
00589             fprintf( stderr, "%3d ", j );
00590             for ( i = 0; i < a->ni; ++i )
00591             {
00592                 square* s = a->squares[j][i];
00593 
00594                 if ( s->npoints > 0 )
00595                     fprintf( stderr, "%3d ", s->npoints );
00596                 else
00597                     fprintf( stderr, "  . " );
00598             }
00599             fprintf( stderr, "\n" );
00600         }
00601     }
00602 
00603     if ( csa_verbose )
00604         fflush( stderr );
00605 }
00606 
00607 // Returns all squares intersecting with a square with center in t->middle
00608 // and edges of length 2*t->r parallel to X and Y axes.
00609 //
00610 static void getsquares( csa* a, triangle* t, int* n, square*** squares )
00611 {
00612     int imin = (int) floor( ( t->middle.x - t->r - a->xmin ) / t->h );
00613     int imax = (int) ceil( ( t->middle.x + t->r - a->xmin ) / t->h );
00614     int jmin = (int) floor( ( t->middle.y - t->r - a->ymin ) / t->h );
00615     int jmax = (int) ceil( ( t->middle.y + t->r - a->ymin ) / t->h );
00616     int i, j;
00617 
00618     if ( imin < 0 )
00619         imin = 0;
00620     if ( imax >= a->ni )
00621         imax = a->ni - 1;
00622     if ( jmin < 0 )
00623         jmin = 0;
00624     if ( jmax >= a->nj )
00625         jmax = a->nj - 1;
00626 
00627     *n           = 0;
00628     ( *squares ) = malloc( (size_t) ( ( imax - imin + 1 ) * ( jmax - jmin + 1 ) ) * sizeof ( square* ) );
00629 
00630     for ( j = jmin; j <= jmax; ++j )
00631     {
00632         for ( i = imin; i <= imax; ++i )
00633         {
00634             square* s = a->squares[j][i];
00635 
00636             if ( s->npoints > 0 )
00637             {
00638                 ( *squares )[*n] = a->squares[j][i];
00639                 ( *n )++;
00640             }
00641         }
00642     }
00643 }
00644 
00645 static double distance( point* p1, point* p2 )
00646 {
00647     return hypot( p1->x - p2->x, p1->y - p2->y );
00648 }
00649 
00650 // Thins data by creating an auxiliary regular grid and for leaving only
00651 // the most central point within each grid cell.
00652 // (I follow the paper here. It is possible that taking average -- in terms of
00653 // both value and position -- of all points within a cell would be a bit more
00654 // robust. However, because of keeping only shallow copies of input points,
00655 // this would require quite a bit of structural changes. So, leaving it as is
00656 // for now.)
00657 //
00658 static void thindata( triangle* t, int npmax )
00659 {
00660     csa    * a         = t->parent->parent;
00661     int    imax        = (int) ceil( sqrt( (double) ( npmax * 3 / 2 ) ) );
00662     square *** squares = alloc2d( imax, imax, sizeof ( void* ) );
00663     double h           = t->r * 2.0 / imax;
00664     double h2          = h / 2.0;
00665     double xmin        = t->middle.x - t->r;
00666     double ymin        = t->middle.y - t->r;
00667     int    i, j, ii;
00668 
00669     for ( j = 0; j < imax; ++j )
00670         for ( i = 0; i < imax; ++i )
00671             squares[j][i] = square_create( a, xmin + h * i, ymin + h * j, i, j );
00672 
00673     for ( ii = 0; ii < t->npoints; ++ii )
00674     {
00675         square* s;
00676         point * p = t->points[ii];
00677         i = (int) floor( ( p->x - xmin ) / h );
00678         j = (int) floor( ( p->y - ymin ) / h );
00679         s = squares[j][i];
00680 
00681         if ( s->npoints == 0 )
00682             square_addpoint( s, p );
00683         else                    // npoints == 1
00684 
00685         {
00686             point pmiddle;
00687 
00688             pmiddle.x = xmin + h * i + h2;
00689             pmiddle.y = ymin + h * j + h2;
00690 
00691             if ( distance( s->points[0], &pmiddle ) > distance( p, &pmiddle ) )
00692                 s->points[0] = p;
00693         }
00694     }
00695 
00696     t->npoints = 0;
00697     for ( j = 0; j < imax; ++j )
00698     {
00699         for ( i = 0; i < imax; ++i )
00700         {
00701             square* s = squares[j][i];
00702 
00703             if ( squares[j][i]->npoints != 0 )
00704                 triangle_addpoint( t, s->points[0] );
00705             square_destroy( s );
00706         }
00707     }
00708 
00709     free2d( squares );
00710     imax++;
00711 }
00712 
00713 // Finds data points to be used in calculating spline coefficients for each
00714 // primary triangle.
00715 //
00716 static void csa_attachpoints( csa* a )
00717 {
00718     int npmin      = a->npmin;
00719     int npmax      = a->npmax;
00720     int nincreased = 0;
00721     int nthinned   = 0;
00722     int i;
00723 
00724     assert( a->npt > 0 );
00725 
00726     if ( csa_verbose )
00727     {
00728         fprintf( stderr, "pre-processing data points:\n  " );
00729         fflush( stderr );
00730     }
00731 
00732     for ( i = 0; i < a->npt; ++i )
00733     {
00734         triangle* t       = a->pt[i];
00735         int     increased = 0;
00736 
00737         if ( csa_verbose )
00738         {
00739             fprintf( stderr, "." );
00740             fflush( stderr );
00741         }
00742 
00743         t->r = t->h * 1.25;
00744         while ( 1 )
00745         {
00746             int   nsquares   = 0;
00747             square** squares = NULL;
00748             int   ii;
00749 
00750             getsquares( a, t, &nsquares, &squares );
00751             for ( ii = 0; ii < nsquares; ++ii )
00752             {
00753                 square* s = squares[ii];
00754                 int   iii;
00755 
00756                 for ( iii = 0; iii < s->npoints; ++iii )
00757                 {
00758                     point* p = s->points[iii];
00759 
00760                     if ( distance( p, &t->middle ) <= t->r )
00761                         triangle_addpoint( t, p );
00762                 }
00763             }
00764 
00765             free( squares );
00766 
00767             if ( t->npoints < npmin )
00768             {
00769                 if ( !increased )
00770                 {
00771                     increased = 1;
00772                     nincreased++;
00773                 }
00774                 t->r      *= 1.25;
00775                 t->npoints = 0;
00776             }
00777             else if ( t->npoints > npmax )
00778             {
00779                 nthinned++;
00780                 thindata( t, npmax );
00781                 if ( t->npoints > npmin )
00782                     break;
00783                 else
00784                 {
00785                     //
00786                     // Sometimes you have too much data, you thin it and --
00787                     // oops -- you have too little. This is not a frequent
00788                     // event, so let us not bother to put a new subdivision.
00789                     //
00790                     t->r      *= 1.25;
00791                     t->npoints = 0;
00792                 }
00793             }
00794             else
00795                 break;
00796         }
00797     }
00798 
00799     if ( csa_verbose )
00800     {
00801         fprintf( stderr, "\n  %d sets enhanced, %d sets thinned\n", nincreased, nthinned );
00802         fflush( stderr );
00803     }
00804 }
00805 
00806 static int n2q( int n )
00807 {
00808     assert( n >= 3 );
00809 
00810     if ( n >= 10 )
00811         return 3;
00812     else if ( n >= 6 )
00813         return 2;
00814     else                        // n == 3
00815         return 1;
00816 }
00817 
00818 //* Singular value decomposition.
00819 // Borrowed from EISPACK (1972-1973).
00820 // Presents input matrix A as  A = U.W.V'.
00821 //
00822 // @param a Input matrix A = U.W[0..m-1][0..n-1]; output matrix U
00823 // @param n Number of columns
00824 // @param m Number of rows
00825 // @param w Ouput vector that presents diagonal matrix W
00826 // @param V output matrix V
00827 //
00828 static void svd( double** a, int n, int m, double* w, double** v )
00829 {
00830     double * rv1;
00831     int    i, j, k, l = -1;
00832     double tst1, c, f, g, h, s, scale;
00833 
00834     assert( m > 0 && n > 0 );
00835 
00836     rv1 = malloc( (size_t) n * sizeof ( double ) );
00837 
00838     //
00839     // householder reduction to bidiagonal form
00840     //
00841     g     = 0.0;
00842     scale = 0.0;
00843     tst1  = 0.0;
00844     for ( i = 0; i < n; i++ )
00845     {
00846         l      = i + 1;
00847         rv1[i] = scale * g;
00848         g      = 0.0;
00849         s      = 0.0;
00850         scale  = 0.0;
00851         if ( i < m )
00852         {
00853             for ( k = i; k < m; k++ )
00854                 scale += fabs( a[k][i] );
00855             if ( scale != 0.0 )
00856             {
00857                 for ( k = i; k < m; k++ )
00858                 {
00859                     a[k][i] /= scale;
00860                     s       += a[k][i] * a[k][i];
00861                 }
00862                 f       = a[i][i];
00863                 g       = -copysign( sqrt( s ), f );
00864                 h       = f * g - s;
00865                 a[i][i] = f - g;
00866                 if ( i < n - 1 )
00867                 {
00868                     for ( j = l; j < n; j++ )
00869                     {
00870                         s = 0.0;
00871                         for ( k = i; k < m; k++ )
00872                             s += a[k][i] * a[k][j];
00873                         f = s / h;
00874                         for ( k = i; k < m; k++ )
00875                             a[k][j] += f * a[k][i];
00876                     }
00877                 }
00878                 for ( k = i; k < m; k++ )
00879                     a[k][i] *= scale;
00880             }
00881         }
00882         w[i]  = scale * g;
00883         g     = 0.0;
00884         s     = 0.0;
00885         scale = 0.0;
00886         if ( i < m && i < n - 1 )
00887         {
00888             for ( k = l; k < n; k++ )
00889                 scale += fabs( a[i][k] );
00890             if ( scale != 0.0 )
00891             {
00892                 for ( k = l; k < n; k++ )
00893                 {
00894                     a[i][k] /= scale;
00895                     s       += a[i][k] * a[i][k];
00896                 }
00897                 f       = a[i][l];
00898                 g       = -copysign( sqrt( s ), f );
00899                 h       = f * g - s;
00900                 a[i][l] = f - g;
00901                 for ( k = l; k < n; k++ )
00902                     rv1[k] = a[i][k] / h;
00903                 for ( j = l; j < m; j++ )
00904                 {
00905                     s = 0.0;
00906                     for ( k = l; k < n; k++ )
00907                         s += a[j][k] * a[i][k];
00908                     for ( k = l; k < n; k++ )
00909                         a[j][k] += s * rv1[k];
00910                 }
00911                 for ( k = l; k < n; k++ )
00912                     a[i][k] *= scale;
00913             }
00914         }
00915         {
00916             double tmp = fabs( w[i] ) + fabs( rv1[i] );
00917 
00918             tst1 = ( tst1 > tmp ) ? tst1 : tmp;
00919         }
00920     }
00921 
00922     //
00923     // accumulation of right-hand transformations
00924     //
00925     for ( i = n - 1; i >= 0; i-- )
00926     {
00927         if ( i < n - 1 )
00928         {
00929             if ( g != 0.0 )
00930             {
00931                 for ( j = l; j < n; j++ )
00932                     //
00933                     // double division avoids possible underflow
00934                     //
00935                     v[j][i] = ( a[i][j] / a[i][l] ) / g;
00936                 for ( j = l; j < n; j++ )
00937                 {
00938                     s = 0.0;
00939                     for ( k = l; k < n; k++ )
00940                         s += a[i][k] * v[k][j];
00941                     for ( k = l; k < n; k++ )
00942                         v[k][j] += s * v[k][i];
00943                 }
00944             }
00945             for ( j = l; j < n; j++ )
00946             {
00947                 v[i][j] = 0.0;
00948                 v[j][i] = 0.0;
00949             }
00950         }
00951         v[i][i] = 1.0;
00952         g       = rv1[i];
00953         l       = i;
00954     }
00955 
00956     //
00957     // accumulation of left-hand transformations
00958     //
00959     for ( i = ( m < n ) ? m - 1 : n - 1; i >= 0; i-- )
00960     {
00961         l = i + 1;
00962         g = w[i];
00963         if ( i != n - 1 )
00964             for ( j = l; j < n; j++ )
00965                 a[i][j] = 0.0;
00966         if ( g != 0.0 )
00967         {
00968             for ( j = l; j < n; j++ )
00969             {
00970                 s = 0.0;
00971                 for ( k = l; k < m; k++ )
00972                     s += a[k][i] * a[k][j];
00973                 //
00974                 // double division avoids possible underflow
00975                 //
00976                 f = ( s / a[i][i] ) / g;
00977                 for ( k = i; k < m; k++ )
00978                     a[k][j] += f * a[k][i];
00979             }
00980             for ( j = i; j < m; j++ )
00981                 a[j][i] /= g;
00982         }
00983         else
00984             for ( j = i; j < m; j++ )
00985                 a[j][i] = 0.0;
00986         a[i][i] += 1.0;
00987     }
00988 
00989     //
00990     // diagonalization of the bidiagonal form
00991     //
00992     for ( k = n - 1; k >= 0; k-- )
00993     {
00994         int k1  = k - 1;
00995         int its = 0;
00996 
00997         while ( 1 )
00998         {
00999             int    docancellation = 1;
01000             double x, y, z;
01001             int    l1 = -1;
01002 
01003             its++;
01004             if ( its > SVD_NMAX )
01005                 csa_quit( "svd(): no convergence in %d iterations", SVD_NMAX );
01006 
01007             for ( l = k; l >= 0; l-- )    // test for splitting
01008             {
01009                 double tst2 = fabs( rv1[l] ) + tst1;
01010 
01011                 if ( tst2 == tst1 )
01012                 {
01013                     docancellation = 0;
01014                     break;
01015                 }
01016                 l1 = l - 1;
01017                 //
01018                 // rv1(1) is always zero, so there is no exit through the
01019                 // bottom of the loop
01020                 //
01021                 tst2 = fabs( w[l - 1] ) + tst1;
01022                 if ( tst2 == tst1 )
01023                     break;
01024             }
01025             //
01026             // cancellation of rv1[l] if l > 1
01027             //
01028             if ( docancellation )
01029             {
01030                 c = 0.0;
01031                 s = 1.0;
01032                 for ( i = l; i <= k; i++ )
01033                 {
01034                     f      = s * rv1[i];
01035                     rv1[i] = c * rv1[i];
01036                     if ( ( fabs( f ) + tst1 ) == tst1 )
01037                         break;
01038                     g    = w[i];
01039                     h    = hypot( f, g );
01040                     w[i] = h;
01041                     h    = 1.0 / h;
01042                     c    = g * h;
01043                     s    = -f * h;
01044                     for ( j = 0; j < m; j++ )
01045                     {
01046                         y = a[j][l1];
01047                         z = a[j][i];
01048 
01049                         a[j][l1] = y * c + z * s;
01050                         a[j][i]  = z * c - y * s;
01051                     }
01052                 }
01053             }
01054             //
01055             // test for convergence
01056             //
01057             z = w[k];
01058             if ( l != k )
01059             {
01060                 int i1;
01061 
01062                 //
01063                 // shift from bottom 2 by 2 minor
01064                 //
01065                 x = w[l];
01066                 y = w[k1];
01067                 g = rv1[k1];
01068                 h = rv1[k];
01069                 f = 0.5 * ( ( ( g + z ) / h ) * ( ( g - z ) / y ) + y / h - h / y );
01070                 g = hypot( f, 1.0 );
01071                 f = x - ( z / x ) * z + ( h / x ) * ( y / ( f + copysign( g, f ) ) - h );
01072                 //
01073                 // next qr transformation
01074                 //
01075                 c = 1.0;
01076                 s = 1.0;
01077                 for ( i1 = l; i1 < k; i1++ )
01078                 {
01079                     i       = i1 + 1;
01080                     g       = rv1[i];
01081                     y       = w[i];
01082                     h       = s * g;
01083                     g       = c * g;
01084                     z       = hypot( f, h );
01085                     rv1[i1] = z;
01086                     c       = f / z;
01087                     s       = h / z;
01088                     f       = x * c + g * s;
01089                     g       = g * c - x * s;
01090                     h       = y * s;
01091                     y      *= c;
01092                     for ( j = 0; j < n; j++ )
01093                     {
01094                         x        = v[j][i1];
01095                         z        = v[j][i];
01096                         v[j][i1] = x * c + z * s;
01097                         v[j][i]  = z * c - x * s;
01098                     }
01099                     z     = hypot( f, h );
01100                     w[i1] = z;
01101                     //
01102                     // rotation can be arbitrary if z = 0
01103                     //
01104                     if ( z != 0.0 )
01105                     {
01106                         c = f / z;
01107                         s = h / z;
01108                     }
01109                     f = c * g + s * y;
01110                     x = c * y - s * g;
01111                     for ( j = 0; j < m; j++ )
01112                     {
01113                         y        = a[j][i1];
01114                         z        = a[j][i];
01115                         a[j][i1] = y * c + z * s;
01116                         a[j][i]  = z * c - y * s;
01117                     }
01118                 }
01119                 rv1[l] = 0.0;
01120                 rv1[k] = f;
01121                 w[k]   = x;
01122             }
01123             else
01124             {
01125                 //
01126                 // w[k] is made non-negative
01127                 //
01128                 if ( z < 0.0 )
01129                 {
01130                     w[k] = -z;
01131                     for ( j = 0; j < n; j++ )
01132                         v[j][k] = -v[j][k];
01133                 }
01134                 break;
01135             }
01136         }
01137     }
01138 
01139     free( rv1 );
01140 }
01141 
01142 // Least squares fitting via singular value decomposition.
01143 //
01144 static void lsq( double** A, int ni, int nj, double* z, double* w, double* sol )
01145 {
01146     double** V = alloc2d( ni, ni, sizeof ( double ) );
01147     double** B = alloc2d( nj, ni, sizeof ( double ) );
01148     int   i, j, ii;
01149 
01150     svd( A, ni, nj, w, V );
01151 
01152     for ( j = 0; j < ni; ++j )
01153         for ( i = 0; i < ni; ++i )
01154             V[j][i] /= w[i];
01155     for ( i = 0; i < ni; ++i )
01156     {
01157         double* v = V[i];
01158 
01159         for ( j = 0; j < nj; ++j )
01160         {
01161             double* a = A[j];
01162             double* b = &B[i][j];
01163 
01164             for ( ii = 0; ii < ni; ++ii )
01165                 *b += v[ii] * a[ii];
01166         }
01167     }
01168     for ( i = 0; i < ni; ++i )
01169         sol[i] = 0.0;
01170     for ( i = 0; i < ni; ++i )
01171         for ( j = 0; j < nj; ++j )
01172             sol[i] += B[i][j] * z[j];
01173 
01174     free2d( B );
01175     free2d( V );
01176 }
01177 
01178 //
01179 //  square->coeffs[]:
01180 //
01181 //   ---------------------
01182 //  | 3    10    17    24 |
01183 //  |    6    13    20    |
01184 //  | 2     9    16    23 |
01185 //  |    5    12    19    |
01186 //  | 1     8    15    22 |
01187 //  |    4    11    18    |
01188 //  | 0     7    14    21 |
01189 //   ---------------------
01190 //
01191 
01192 // Calculates spline coefficients in each primary triangle by least squares
01193 // fitting to data attached by csa_attachpoints().
01194 //
01195 static void csa_findprimarycoeffs( csa* a )
01196 {
01197     int n[4] = { 0, 0, 0, 0 };
01198     int i;
01199 
01200     if ( csa_verbose )
01201         fprintf( stderr, "calculating spline coefficients for primary triangles:\n  " );
01202 
01203     for ( i = 0; i < a->npt; ++i )
01204     {
01205         triangle* t       = a->pt[i];
01206         int     npoints   = t->npoints;
01207         point   ** points = t->points;
01208         double  * z       = malloc( (size_t) npoints * sizeof ( double ) );
01209         int     q         = n2q( t->npoints );
01210         int     ok        = 1;
01211         double  b[10];
01212         double  b1[6];
01213         int     ii;
01214 
01215         if ( csa_verbose )
01216         {
01217             fprintf( stderr, "." );
01218             fflush( stderr );
01219         }
01220 
01221         for ( ii = 0; ii < npoints; ++ii )
01222             z[ii] = points[ii]->z;
01223 
01224         do
01225         {
01226             double bc[3];
01227             double wmin, wmax;
01228 
01229             if ( !ok )
01230                 q--;
01231 
01232             assert( q >= 0 );
01233 
01234             if ( q == 3 )
01235             {
01236                 double ** A = alloc2d( 10, npoints, sizeof ( double ) );
01237                 double w[10];
01238 
01239                 for ( ii = 0; ii < npoints; ++ii )
01240                 {
01241                     double * aii = A[ii];
01242                     double tmp;
01243 
01244                     triangle_calculatebc( t, points[ii], bc );
01245 
01246                     //
01247                     //  0   1   2   3   4   5   6   7   8   9
01248                     // 300 210 201 120 111 102 030 021 012 003
01249                     //
01250                     tmp    = bc[0] * bc[0];
01251                     aii[0] = tmp * bc[0];
01252                     tmp   *= 3.0;
01253                     aii[1] = tmp * bc[1];
01254                     aii[2] = tmp * bc[2];
01255                     tmp    = bc[1] * bc[1];
01256                     aii[6] = tmp * bc[1];
01257                     tmp   *= 3.0;
01258                     aii[3] = tmp * bc[0];
01259                     aii[7] = tmp * bc[2];
01260                     tmp    = bc[2] * bc[2];
01261                     aii[9] = tmp * bc[2];
01262                     tmp   *= 3.0;
01263                     aii[5] = tmp * bc[0];
01264                     aii[8] = tmp * bc[1];
01265                     aii[4] = bc[0] * bc[1] * bc[2] * 6.0;
01266                 }
01267 
01268                 lsq( A, 10, npoints, z, w, b );
01269 
01270                 wmin = w[0];
01271                 wmax = w[0];
01272                 for ( ii = 1; ii < 10; ++ii )
01273                 {
01274                     if ( w[ii] < wmin )
01275                         wmin = w[ii];
01276                     else if ( w[ii] > wmax )
01277                         wmax = w[ii];
01278                 }
01279                 if ( wmin < wmax / a->k )
01280                     ok = 0;
01281 
01282                 free2d( A );
01283             }
01284             else if ( q == 2 )
01285             {
01286                 double ** A = alloc2d( 6, npoints, sizeof ( double ) );
01287                 double w[6];
01288 
01289                 for ( ii = 0; ii < npoints; ++ii )
01290                 {
01291                     double* aii = A[ii];
01292 
01293                     triangle_calculatebc( t, points[ii], bc );
01294 
01295                     //
01296                     //  0   1   2   3   4   5
01297                     // 200 110 101 020 011 002
01298                     //
01299 
01300                     aii[0] = bc[0] * bc[0];
01301                     aii[1] = bc[0] * bc[1] * 2.0;
01302                     aii[2] = bc[0] * bc[2] * 2.0;
01303                     aii[3] = bc[1] * bc[1];
01304                     aii[4] = bc[1] * bc[2] * 2.0;
01305                     aii[5] = bc[2] * bc[2];
01306                 }
01307 
01308                 lsq( A, 6, npoints, z, w, b1 );
01309 
01310                 wmin = w[0];
01311                 wmax = w[0];
01312                 for ( ii = 1; ii < 6; ++ii )
01313                 {
01314                     if ( w[ii] < wmin )
01315                         wmin = w[ii];
01316                     else if ( w[ii] > wmax )
01317                         wmax = w[ii];
01318                 }
01319                 if ( wmin < wmax / a->k )
01320                     ok = 0;
01321                 else            // degree raising
01322                 {
01323                     ok   = 1;
01324                     b[0] = b1[0];
01325                     b[1] = ( b1[0] + 2.0 * b1[1] ) / 3.0;
01326                     b[2] = ( b1[0] + 2.0 * b1[2] ) / 3.0;
01327                     b[3] = ( b1[3] + 2.0 * b1[1] ) / 3.0;
01328                     b[4] = ( b1[1] + b1[2] + b1[4] ) / 3.0;
01329                     b[5] = ( b1[5] + 2.0 * b1[2] ) / 3.0;
01330                     b[6] = b1[3];
01331                     b[7] = ( b1[3] + 2.0 * b1[4] ) / 3.0;
01332                     b[8] = ( b1[5] + 2.0 * b1[4] ) / 3.0;
01333                     b[9] = b1[5];
01334                 }
01335 
01336                 free2d( A );
01337             }
01338             else if ( q == 1 )
01339             {
01340                 double ** A = alloc2d( 3, npoints, sizeof ( double ) );
01341                 double w[3];
01342 
01343                 for ( ii = 0; ii < npoints; ++ii )
01344                 {
01345                     double* aii = A[ii];
01346 
01347                     triangle_calculatebc( t, points[ii], bc );
01348 
01349                     aii[0] = bc[0];
01350                     aii[1] = bc[1];
01351                     aii[2] = bc[2];
01352                 }
01353 
01354                 lsq( A, 3, npoints, z, w, b1 );
01355 
01356                 wmin = w[0];
01357                 wmax = w[0];
01358                 for ( ii = 1; ii < 3; ++ii )
01359                 {
01360                     if ( w[ii] < wmin )
01361                         wmin = w[ii];
01362                     else if ( w[ii] > wmax )
01363                         wmax = w[ii];
01364                 }
01365                 if ( wmin < wmax / a->k )
01366                     ok = 0;
01367                 else            // degree raising
01368                 {
01369                     ok   = 1;
01370                     b[0] = b1[0];
01371                     b[1] = ( 2.0 * b1[0] + b1[1] ) / 3.0;
01372                     b[2] = ( 2.0 * b1[0] + b1[2] ) / 3.0;
01373                     b[3] = ( 2.0 * b1[1] + b1[0] ) / 3.0;
01374                     b[4] = ( b1[0] + b1[1] + b1[2] ) / 3.0;
01375                     b[5] = ( 2.0 * b1[2] + b1[0] ) / 3.0;
01376                     b[6] = b1[1];
01377                     b[7] = ( 2.0 * b1[1] + b1[2] ) / 3.0;
01378                     b[8] = ( 2.0 * b1[2] + b1[1] ) / 3.0;
01379                     b[9] = b1[2];
01380                 }
01381 
01382                 free2d( A );
01383             }
01384             else if ( q == 0 )
01385             {
01386                 double ** A = alloc2d( 1, npoints, sizeof ( double ) );
01387                 double w[1];
01388 
01389                 for ( ii = 0; ii < npoints; ++ii )
01390                     A[ii][0] = 1.0;
01391 
01392                 lsq( A, 1, npoints, z, w, b1 );
01393 
01394                 ok   = 1;
01395                 b[0] = b1[0];
01396                 b[1] = b1[0];
01397                 b[2] = b1[0];
01398                 b[3] = b1[0];
01399                 b[4] = b1[0];
01400                 b[5] = b1[0];
01401                 b[6] = b1[0];
01402                 b[7] = b1[0];
01403                 b[8] = b1[0];
01404                 b[9] = b1[0];
01405 
01406                 free2d( A );
01407             }
01408         } while ( !ok );
01409 
01410         n[q]++;
01411         t->order = q;
01412 
01413         {
01414             square* s      = t->parent;
01415             double* coeffs = s->coeffs;
01416 
01417             coeffs[12] = b[0];
01418             coeffs[9]  = b[1];
01419             coeffs[6]  = b[3];
01420             coeffs[3]  = b[6];
01421             coeffs[2]  = b[7];
01422             coeffs[1]  = b[8];
01423             coeffs[0]  = b[9];
01424             coeffs[4]  = b[5];
01425             coeffs[8]  = b[2];
01426             coeffs[5]  = b[4];
01427         }
01428 
01429         free( z );
01430     }
01431 
01432     if ( csa_verbose )
01433     {
01434         fprintf( stderr, "\n  3rd order -- %d sets\n", n[3] );
01435         fprintf( stderr, "  2nd order -- %d sets\n", n[2] );
01436         fprintf( stderr, "  1st order -- %d sets\n", n[1] );
01437         fprintf( stderr, "  0th order -- %d sets\n", n[0] );
01438         fflush( stderr );
01439     }
01440 
01441     if ( csa_verbose == 2 )
01442     {
01443         int j;
01444 
01445         fprintf( stderr, " j\\i" );
01446         for ( i = 0; i < a->ni; ++i )
01447             fprintf( stderr, "%2d ", i );
01448         fprintf( stderr, "\n" );
01449         for ( j = a->nj - 1; j >= 0; --j )
01450         {
01451             fprintf( stderr, "%2d  ", j );
01452             for ( i = 0; i < a->ni; ++i )
01453             {
01454                 square* s = a->squares[j][i];
01455 
01456                 if ( s->triangles[0]->primary )
01457                     fprintf( stderr, "%2d ", s->triangles[0]->order );
01458                 else
01459                     fprintf( stderr, " . " );
01460             }
01461             fprintf( stderr, "\n" );
01462         }
01463     }
01464 }
01465 
01466 // Finds spline coefficients in (adjacent to primary triangles) secondary
01467 // triangles from C1 smoothness conditions.
01468 //
01469 static void csa_findsecondarycoeffs( csa* a )
01470 {
01471     square*** squares = a->squares;
01472     int   ni          = a->ni;
01473     int   nj          = a->nj;
01474     int   ii;
01475 
01476     if ( csa_verbose )
01477     {
01478         fprintf( stderr, "propagating spline coefficients to the remaining triangles:\n" );
01479         fflush( stderr );
01480     }
01481 
01482     //
01483     // red
01484     //
01485     for ( ii = 0; ii < a->npt; ++ii )
01486     {
01487         triangle* t   = a->pt[ii];
01488         square  * s   = t->parent;
01489         int     i     = s->i;
01490         int     j     = s->j;
01491         double  * c   = s->coeffs;
01492         double  * cl  = ( i > 0 ) ? squares[j][i - 1]->coeffs : NULL;
01493         double  * cb  = ( j > 0 ) ? squares[j - 1][i]->coeffs : NULL;
01494         double  * cbl = ( i > 0 && j > 0 ) ? squares[j - 1][i - 1]->coeffs : NULL;
01495         double  * ca  = ( j < nj - 1 ) ? squares[j + 1][i]->coeffs : NULL;
01496         double  * cal = ( j < nj - 1 && i > 0 ) ? squares[j + 1][i - 1]->coeffs : NULL;
01497 
01498         c[7]  = 2.0 * c[4] - c[1];
01499         c[11] = 2.0 * c[8] - c[5];
01500         c[15] = 2.0 * c[12] - c[9];
01501 
01502         c[10] = 2.0 * c[6] - c[2];
01503         c[13] = 2.0 * c[9] - c[5];
01504         c[16] = 2.0 * c[12] - c[8];
01505 
01506         c[19] = 2.0 * c[15] - c[11];
01507 
01508         if ( cl != NULL )
01509         {
01510             cl[21] = c[0];
01511             cl[22] = c[1];
01512             cl[23] = c[2];
01513             cl[24] = c[3];
01514 
01515             cl[18] = c[0] + c[1] - c[4];
01516             cl[19] = c[1] + c[2] - c[5];
01517             cl[20] = c[2] + c[3] - c[6];
01518 
01519             cl[17] = 2.0 * cl[20] - cl[23];
01520             cl[14] = 2.0 * cl[18] - cl[22];
01521         }
01522 
01523         if ( cb != NULL )
01524         {
01525             cb[3]  = c[0];
01526             cb[10] = c[7];
01527 
01528             cb[6] = c[0] + c[7] - c[4];
01529             cb[2] = 2.0 * cb[6] - cb[10];
01530         }
01531 
01532         if ( cbl != NULL )
01533         {
01534             cbl[23] = cb[2];
01535             cbl[24] = cb[3];
01536 
01537             cbl[20] = cb[2] + cb[3] - cb[6];
01538             cbl[17] = cl[14];
01539         }
01540 
01541         if ( ca != NULL )
01542         {
01543             ca[0] = c[3];
01544             ca[7] = c[10];
01545 
01546             ca[4] = c[3] + c[10] - c[6];
01547             ca[1] = 2.0 * ca[4] - ca[7];
01548         }
01549 
01550         if ( cal != NULL )
01551         {
01552             cal[21] = c[3];
01553             cal[22] = ca[1];
01554 
01555             cal[18] = ca[0] + ca[1] - ca[4];
01556             cal[14] = cl[17];
01557         }
01558     }
01559 
01560     //
01561     // blue
01562     //
01563     for ( ii = 0; ii < a->npt; ++ii )
01564     {
01565         triangle* t   = a->pt[ii];
01566         square  * s   = t->parent;
01567         int     i     = s->i;
01568         int     j     = s->j;
01569         double  * c   = s->coeffs;
01570         double  * cr  = ( i < ni - 1 ) ? squares[j][i + 1]->coeffs : NULL;
01571         double  * car = ( i < ni - 1 && j < nj - 1 ) ? squares[j + 1][i + 1]->coeffs : NULL;
01572         double  * cbr = ( i < ni - 1 && j > 0 ) ? squares[j - 1][i + 1]->coeffs : NULL;
01573 
01574         if ( car != NULL )
01575             cr[13] = car[7] + car[14] - car[11];
01576 
01577         if ( cbr != NULL )
01578             cr[11] = cbr[10] + cbr[17] - cbr[13];
01579 
01580         if ( cr != NULL )
01581             cr[5] = c[22] + c[23] - c[19];
01582     }
01583 
01584     //
01585     // green & yellow
01586     //
01587     for ( ii = 0; ii < a->npt; ++ii )
01588     {
01589         triangle* t  = a->pt[ii];
01590         square  * s  = t->parent;
01591         int     i    = s->i;
01592         int     j    = s->j;
01593         double  * cr = ( i < ni - 1 ) ? squares[j][i + 1]->coeffs : NULL;
01594 
01595         if ( cr != NULL )
01596         {
01597             cr[9]  = ( cr[5] + cr[13] ) / 2.0;
01598             cr[8]  = ( cr[5] + cr[11] ) / 2.0;
01599             cr[15] = ( cr[11] + cr[19] ) / 2.0;
01600             cr[16] = ( cr[13] + cr[19] ) / 2.0;
01601             cr[12] = ( cr[8] + cr[16] ) / 2.0;
01602         }
01603     }
01604 
01605     if ( csa_verbose )
01606     {
01607         fprintf( stderr, "checking that all coefficients have been set:\n" );
01608         fflush( stderr );
01609     }
01610 
01611     for ( ii = 0; ii < ni * nj; ++ii )
01612     {
01613         square* s = squares[0][ii];
01614         double* c = s->coeffs;
01615         int   i;
01616 
01617         if ( s->npoints == 0 )
01618             continue;
01619         for ( i = 0; i < 25; ++i )
01620             if ( isnan( c[i] ) )
01621                 fprintf( stderr, "  squares[%d][%d]->coeffs[%d] = NaN\n", s->j, s->i, i );
01622     }
01623 }
01624 
01625 static int i300[] = { 12, 12, 12, 12 };
01626 static int i030[] = { 3, 24, 21, 0 };
01627 static int i003[] = { 0, 3, 24, 21 };
01628 static int i210[] = { 9, 16, 15, 8 };
01629 static int i021[] = { 2, 17, 22, 7 };
01630 static int i102[] = { 4, 6, 20, 18 };
01631 static int i120[] = { 6, 20, 18, 4 };
01632 static int i012[] = { 1, 10, 23, 14 };
01633 static int i201[] = { 8, 9, 16, 15 };
01634 static int i111[] = { 5, 13, 19, 11 };
01635 
01636 static int * iall[] = { i300, i030, i003, i210, i021, i102, i120, i012, i201, i111 };
01637 
01638 static void csa_sethascoeffsflag( csa* a )
01639 {
01640     int i, j;
01641 
01642     for ( j = 0; j < a->nj; ++j )
01643     {
01644         for ( i = 0; i < a->ni; ++i )
01645         {
01646             square* s      = a->squares[j][i];
01647             double* coeffs = s->coeffs;
01648             int   ii;
01649 
01650             for ( ii = 0; ii < 4; ++ii )
01651             {
01652                 triangle* t = s->triangles[ii];
01653                 int     cc;
01654 
01655                 for ( cc = 0; cc < 10; ++cc )
01656                     if ( isnan( coeffs[iall[cc][ii]] ) )
01657                         break;
01658                 if ( cc == 10 )
01659                     t->hascoeffs = 1;
01660             }
01661         }
01662     }
01663 }
01664 
01665 void csa_calculatespline( csa* a )
01666 {
01667     csa_squarize( a );
01668     csa_attachpoints( a );
01669     csa_findprimarycoeffs( a );
01670     csa_findsecondarycoeffs( a );
01671     csa_sethascoeffsflag( a );
01672 }
01673 
01674 void csa_approximate_point( csa* a, point* p )
01675 {
01676     double  h  = a->h;
01677     double  ii = ( p->x - a->xmin ) / h + 1.0;
01678     double  jj = ( p->y - a->ymin ) / h + 1.0;
01679     int     i, j;
01680     square  * s;
01681     double  fi, fj;
01682     int     ti;
01683     triangle* t;
01684     double  bc[3];
01685 
01686     if ( ii < 0.0 || jj < 0.0 || ii > (double) a->ni - 1.0 || jj > (double) a->nj - 1.0 )
01687     {
01688         p->z = NaN;
01689         return;
01690     }
01691 
01692     i  = (int) floor( ii );
01693     j  = (int) floor( jj );
01694     s  = a->squares[j][i];
01695     fi = ii - i;
01696     fj = jj - j;
01697 
01698     if ( fj < fi )
01699     {
01700         if ( fi + fj < 1.0 )
01701             ti = 3;
01702         else
01703             ti = 2;
01704     }
01705     else
01706     {
01707         if ( fi + fj < 1.0 )
01708             ti = 0;
01709         else
01710             ti = 1;
01711     }
01712 
01713     t = s->triangles[ti];
01714     if ( !t->hascoeffs )
01715     {
01716         p->z = NaN;
01717         return;
01718     }
01719     triangle_calculatebc( t, p, bc );
01720 
01721     {
01722         double * c  = s->coeffs;
01723         double bc1  = bc[0];
01724         double bc2  = bc[1];
01725         double bc3  = bc[2];
01726         double tmp1 = bc1 * bc1;
01727         double tmp2 = bc2 * bc2;
01728         double tmp3 = bc3 * bc3;
01729 
01730         switch ( ti )
01731         {
01732         case 0:
01733             p->z = c[12] * bc1 * tmp1 + c[3] * bc2 * tmp2 + c[0] * bc3 * tmp3 + 3.0 * ( c[9] * tmp1 * bc2 + c[2] * tmp2 * bc3 + c[4] * tmp3 * bc1 + c[6] * bc1 * tmp2 + c[1] * bc2 * tmp3 + c[8] * tmp1 * bc3 ) + 6.0 * c[5] * bc1 * bc2 * bc3;
01734             break;
01735         case 1:
01736             p->z = c[12] * bc1 * tmp1 + c[24] * bc2 * tmp2 + c[3] * bc3 * tmp3 + 3.0 * ( c[16] * tmp1 * bc2 + c[17] * tmp2 * bc3 + c[6] * tmp3 * bc1 + c[20] * bc1 * tmp2 + c[10] * bc2 * tmp3 + c[9] * tmp1 * bc3 ) + 6.0 * c[13] * bc1 * bc2 * bc3;
01737             break;
01738         case 2:
01739             p->z = c[12] * bc1 * tmp1 + c[21] * bc2 * tmp2 + c[24] * bc3 * tmp3 + 3.0 * ( c[15] * tmp1 * bc2 + c[22] * tmp2 * bc3 + c[20] * tmp3 * bc1 + c[18] * bc1 * tmp2 + c[23] * bc2 * tmp3 + c[16] * tmp1 * bc3 ) + 6.0 * c[19] * bc1 * bc2 * bc3;
01740             break;
01741         default:               // 3
01742             p->z = c[12] * bc1 * tmp1 + c[0] * bc2 * tmp2 + c[21] * bc3 * tmp3 + 3.0 * ( c[8] * tmp1 * bc2 + c[7] * tmp2 * bc3 + c[18] * tmp3 * bc1 + c[4] * bc1 * tmp2 + c[14] * bc2 * tmp3 + c[15] * tmp1 * bc3 ) + 6.0 * c[11] * bc1 * bc2 * bc3;
01743         }
01744     }
01745 }
01746 
01747 void csa_approximate_points( csa* a, int n, point* points )
01748 {
01749     int ii;
01750 
01751     for ( ii = 0; ii < n; ++ii )
01752         csa_approximate_point( a, &points[ii] );
01753 }
01754 
01755 void csa_setnpmin( csa* a, int npmin )
01756 {
01757     a->npmin = npmin;
01758 }
01759 
01760 void csa_setnpmax( csa* a, int npmax )
01761 {
01762     a->npmax = npmax;
01763 }
01764 
01765 void csa_setk( csa* a, int k )
01766 {
01767     a->k = k;
01768 }
01769 
01770 void csa_setnppc( csa* a, double nppc )
01771 {
01772     a->nppc = (int) nppc;
01773 }
01774 
01775 #if defined ( STANDALONE )
01776 
01777 #include "minell.h"
01778 
01779 #define NIMAX         2048
01780 #define BUFSIZE       10240
01781 #define STRBUFSIZE    64
01782 
01783 static void points_generate( double xmin, double xmax, double ymin, double ymax, int nx, int ny, int* nout, point** pout )
01784 {
01785     double stepx, stepy;
01786     double x0, xx, yy;
01787     int    i, j, ii;
01788 
01789     if ( nx < 1 || ny < 1 )
01790     {
01791         *pout = NULL;
01792         *nout = 0;
01793         return;
01794     }
01795 
01796     *nout = nx * ny;
01797     *pout = malloc( *nout * sizeof ( point ) );
01798 
01799     stepx = ( nx > 1 ) ? ( xmax - xmin ) / ( nx - 1 ) : 0.0;
01800     stepy = ( ny > 1 ) ? ( ymax - ymin ) / ( ny - 1 ) : 0.0;
01801     x0    = ( nx > 1 ) ? xmin : ( xmin + xmax ) / 2.0;
01802     yy    = ( ny > 1 ) ? ymin : ( ymin + ymax ) / 2.0;
01803 
01804     ii = 0;
01805     for ( j = 0; j < ny; ++j )
01806     {
01807         xx = x0;
01808         for ( i = 0; i < nx; ++i )
01809         {
01810             point* p = &( *pout )[ii];
01811 
01812             p->x = xx;
01813             p->y = yy;
01814             xx  += stepx;
01815             ii++;
01816         }
01817         yy += stepy;
01818     }
01819 }
01820 
01821 static int str2double( char* token, double* value )
01822 {
01823     char* end = NULL;
01824 
01825     if ( token == NULL )
01826     {
01827         *value = NaN;
01828         return 0;
01829     }
01830 
01831     *value = strtod( token, &end );
01832 
01833     if ( end == token )
01834     {
01835         *value = NaN;
01836         return 0;
01837     }
01838 
01839     return 1;
01840 }
01841 
01842 #define NALLOCATED_START    1024
01843 
01844 // Reads array of points from a columnar file.
01845 //
01846 // @param fname File name (can be "stdin" or "-" for standard input)
01847 // @param dim Number of dimensions (must be 2 or 3)
01848 // @param n Pointer to number of points (output)
01849 // @param points Pointer to array of points [*n] (output) (to be freed)
01850 //
01851 void points_read( char* fname, int dim, int* n, point** points )
01852 {
01853     FILE * f        = NULL;
01854     int  nallocated = NALLOCATED_START;
01855     char buf[BUFSIZE];
01856     char seps[] = " ,;\t";
01857     char * token;
01858 
01859     if ( dim < 2 || dim > 3 )
01860     {
01861         *n      = 0;
01862         *points = NULL;
01863         return;
01864     }
01865 
01866     if ( fname == NULL )
01867         f = stdin;
01868     else
01869     {
01870         if ( strcmp( fname, "stdin" ) == 0 || strcmp( fname, "-" ) == 0 )
01871             f = stdin;
01872         else
01873         {
01874             f = fopen( fname, "r" );
01875             if ( f == NULL )
01876                 csa_quit( "%s: %s\n", fname, strerror( errno ) );
01877         }
01878     }
01879 
01880     *points = malloc( nallocated * sizeof ( point ) );
01881     *n      = 0;
01882     while ( fgets( buf, BUFSIZE, f ) != NULL )
01883     {
01884         point* p;
01885 
01886         if ( *n == nallocated )
01887         {
01888             nallocated *= 2;
01889             *points     = realloc( *points, nallocated * sizeof ( point ) );
01890         }
01891 
01892         p = &( *points )[*n];
01893 
01894         if ( buf[0] == '#' )
01895             continue;
01896         if ( ( token = strtok( buf, seps ) ) == NULL )
01897             continue;
01898         if ( !str2double( token, &p->x ) )
01899             continue;
01900         if ( ( token = strtok( NULL, seps ) ) == NULL )
01901             continue;
01902         if ( !str2double( token, &p->y ) )
01903             continue;
01904         if ( dim == 2 )
01905             p->z = NaN;
01906         else
01907         {
01908             if ( ( token = strtok( NULL, seps ) ) == NULL )
01909                 continue;
01910             if ( !str2double( token, &p->z ) )
01911                 continue;
01912         }
01913         ( *n )++;
01914     }
01915 
01916     if ( *n == 0 )
01917     {
01918         free( *points );
01919         *points = NULL;
01920     }
01921     else
01922         *points = realloc( *points, *n * sizeof ( point ) );
01923 
01924     if ( f != stdin )
01925         if ( fclose( f ) != 0 )
01926             csa_quit( "%s: %s\n", fname, strerror( errno ) );
01927 }
01928 
01929 static void points_write( int n, point* points )
01930 {
01931     int i;
01932 
01933     if ( csa_verbose )
01934         printf( "Output:\n" );
01935 
01936     for ( i = 0; i < n; ++i )
01937     {
01938         point* p = &points[i];
01939 
01940         printf( "%.15g %.15g %.15g\n", p->x, p->y, p->z );
01941     }
01942 }
01943 
01944 static double points_scaletosquare( int n, point* points )
01945 {
01946     double xmin, ymin, xmax, ymax;
01947     double k;
01948     int    i;
01949 
01950     if ( n <= 0 )
01951         return NaN;
01952 
01953     xmin = xmax = points[0].x;
01954     ymin = ymax = points[0].y;
01955 
01956     for ( i = 1; i < n; ++i )
01957     {
01958         point* p = &points[i];
01959 
01960         if ( p->x < xmin )
01961             xmin = p->x;
01962         else if ( p->x > xmax )
01963             xmax = p->x;
01964         if ( p->y < ymin )
01965             ymin = p->y;
01966         else if ( p->y > ymax )
01967             ymax = p->y;
01968     }
01969 
01970     if ( xmin == xmax || ymin == ymax )
01971         return NaN;
01972     else
01973         k = ( ymax - ymin ) / ( xmax - xmin );
01974 
01975     for ( i = 0; i < n; ++i )
01976         points[i].y /= k;
01977 
01978     return k;
01979 }
01980 
01981 static void points_scale( int n, point* points, double k )
01982 {
01983     int i;
01984 
01985     for ( i = 0; i < n; ++i )
01986         points[i].y /= k;
01987 }
01988 
01989 static void usage()
01990 {
01991     printf( "Usage: csabathy -i <XYZ file> {-o <XY file>|-n <nx>x<ny> [-c|-s] [-z <zoom>]}\n" );
01992     printf( "       [-v|-V] [-P nppc=<value>] [-P k=<value>]\n" );
01993     printf( "Options:\n" );
01994     printf( "  -c              -- scale internally so that the minimal ellipse turns into a\n" );
01995     printf( "                     circle (this produces results invariant to affine\n" );
01996     printf( "                     transformations)\n" );
01997     printf( "  -i <XYZ file>   -- three-column file with points to approximate from (use\n" );
01998     printf( "                     \"-i stdin\" or \"-i -\" for standard input)\n" );
01999     printf( "  -n <nx>x<ny>    -- generate <nx>x<ny> output rectangular grid\n" );
02000     printf( "  -o <XY file>    -- two-column file with points to approximate in (use\n" );
02001     printf( "                     \"-o stdin\" or \"-o -\" for standard input)\n" );
02002     printf( "  -s              -- scale internally so that Xmax - Xmin = Ymax - Ymin\n" );
02003     printf( "  -v              -- verbose / version\n" );
02004     printf( "  -z <zoom>       -- zoom in (if <zoom> < 1) or out (<zoom> > 1)\n" );
02005     printf( "  -P nppc=<value> -- set the average number of points per cell (default = 5\n" );
02006     printf( "                     increase if the point distribution is strongly non-uniform\n" );
02007     printf( "                     to get larger cells)\n" );
02008     printf( "  -P k=<value>    -- set the spline sensitivity (default = 140, reduce to get\n" );
02009     printf( "                     smoother results)\n" );
02010     printf( "  -V              -- very verbose / version\n" );
02011     printf( "Description:\n" );
02012     printf( "  `csabathy' approximates irregular scalar 2D data in specified points using\n" );
02013     printf( "  C1-continuous bivariate cubic spline. The calculated values are written to\n" );
02014     printf( "  standard output.\n" );
02015 
02016     exit( 0 );
02017 }
02018 
02019 static void version()
02020 {
02021     printf( "csa version %s\n", csa_version );
02022     exit( 0 );
02023 }
02024 
02025 static void parse_commandline( int argc, char* argv[], char** fdata, char** fpoints, int* invariant, int* square, int* generate_points, int* nx, int* ny, int* nppc, int* k, double* zoom )
02026 {
02027     int i;
02028 
02029     if ( argc < 2 )
02030         usage();
02031 
02032     i = 1;
02033     while ( i < argc )
02034     {
02035         if ( argv[i][0] != '-' )
02036             usage();
02037 
02038         switch ( argv[i][1] )
02039         {
02040         case 'c':
02041             i++;
02042             *invariant = 1;
02043             *square    = 0;
02044 
02045             break;
02046         case 'i':
02047             i++;
02048             if ( i >= argc )
02049                 csa_quit( "no file name found after -i\n" );
02050             *fdata = argv[i];
02051             i++;
02052             break;
02053         case 'n':
02054             i++;
02055             *fpoints         = NULL;
02056             *generate_points = 1;
02057             if ( i >= argc )
02058                 csa_quit( "no grid dimensions found after -n\n" );
02059             if ( sscanf( argv[i], "%dx%d", nx, ny ) != 2 )
02060                 csa_quit( "could not read grid dimensions after \"-n\"\n" );
02061             if ( *nx <= 0 || *nx > NIMAX || *ny <= 0 || *ny > NIMAX )
02062                 csa_quit( "invalid size for output grid\n" );
02063             i++;
02064             break;
02065         case 'o':
02066             i++;
02067             if ( i >= argc )
02068                 csa_quit( "no file name found after -o\n" );
02069             *fpoints = argv[i];
02070             i++;
02071             break;
02072         case 's':
02073             i++;
02074             *square = 1;
02075 
02076             *invariant = 0;
02077             break;
02078         case 'v':
02079             i++;
02080             csa_verbose = 1;
02081             break;
02082         case 'z':
02083             i++;
02084             if ( i >= argc )
02085                 csa_quit( "no zoom value found after -z\n" );
02086             *zoom = atof( argv[i] );
02087             i++;
02088             break;
02089         case 'P': {
02090             char delim[]            = "=";
02091             char prmstr[STRBUFSIZE] = "";
02092             char * token;
02093 
02094             i++;
02095             if ( i >= argc )
02096                 csa_quit( "no input found after -P\n" );
02097 
02098             if ( strlen( argv[i] ) >= STRBUFSIZE )
02099                 csa_quit( "could not interpret \"%s\" after -P option\n", argv[i] );
02100 
02101             strcpy( prmstr, argv[i] );
02102             token = strtok( prmstr, delim );
02103             if ( token == NULL )
02104                 csa_quit( "could not interpret \"%s\" after -P option\n", argv[i] );
02105 
02106             if ( strcmp( token, "nppc" ) == 0 )
02107             {
02108                 long int n;
02109 
02110                 token = strtok( NULL, delim );
02111                 if ( token == NULL )
02112                     csa_quit( "could not interpret \"%s\" after -P option\n", argv[i] );
02113 
02114                 n = strtol( token, NULL, 10 );
02115                 if ( n == LONG_MIN || n == LONG_MAX )
02116                     csa_quit( "could not interpret \"%s\" after -P option\n", argv[i] );
02117                 else if ( n <= 0 )
02118                     csa_quit( "non-sensible value for \"nppc\" parameter\n" );
02119                 *nppc = (int) n;
02120             }
02121             else if ( strcmp( token, "k" ) == 0 )
02122             {
02123                 long int n;
02124 
02125                 token = strtok( NULL, delim );
02126                 if ( token == NULL )
02127                     csa_quit( "could not interpret \"%s\" after -P option\n", argv[i] );
02128 
02129                 n = strtol( token, NULL, 10 );
02130                 if ( n == LONG_MIN || n == LONG_MAX )
02131                     csa_quit( "could not interpret \"%s\" after -P option\n", argv[i] );
02132                 else if ( n <= 0 )
02133                     csa_quit( "non-sensible value for \"k\" parameter\n" );
02134                 *k = (int) n;
02135             }
02136             else
02137                 usage();
02138 
02139             i++;
02140             break;
02141         }
02142         case 'V':
02143             i++;
02144             csa_verbose = 2;
02145             break;
02146         default:
02147             usage();
02148             break;
02149         }
02150     }
02151 
02152     if ( csa_verbose && argc == 2 )
02153         version();
02154 }
02155 
02156 int main( int argc, char* argv[] )
02157 {
02158     char   * fdata         = NULL;
02159     char   * fpoints       = NULL;
02160     int    nin             = 0;
02161     point  * pin           = NULL;
02162     int    invariant       = 0;
02163     minell * me            = NULL;
02164     int    square          = 0;
02165     int    nout            = 0;
02166     int    generate_points = 0;
02167     point  * pout          = NULL;
02168     int    nx   = -1;
02169     int    ny   = -1;
02170     csa    * a  = NULL;
02171     int    nppc = -1;
02172     int    k    = -1;
02173     double ks   = NaN;
02174     double zoom = NaN;
02175 
02176     parse_commandline( argc, argv, &fdata, &fpoints, &invariant, &square, &generate_points, &nx, &ny, &nppc, &k, &zoom );
02177 
02178     if ( fdata == NULL )
02179         csa_quit( "no input data\n" );
02180 
02181     if ( !generate_points && fpoints == NULL )
02182         csa_quit( "no output grid specified\n" );
02183 
02184     points_read( fdata, 3, &nin, &pin );
02185 
02186     if ( nin < 3 )
02187         return 0;
02188 
02189     if ( invariant )
02190     {
02191         me = minell_build( nin, pin );
02192         minell_scalepoints( me, nin, pin );
02193     }
02194     else if ( square )
02195         ks = points_scaletosquare( nin, pin );
02196 
02197     a = csa_create();
02198     csa_addpoints( a, nin, pin );
02199     if ( nppc > 0 )
02200         csa_setnppc( a, nppc );
02201     if ( k > 0 )
02202         csa_setk( a, k );
02203     csa_calculatespline( a );
02204 
02205     if ( generate_points )
02206     {
02207         if ( isnan( zoom ) )
02208             points_generate( a->xmin - a->h, a->xmax + a->h, a->ymin - a->h, a->ymax + a->h, nx, ny, &nout, &pout );
02209         else
02210         {
02211             double xdiff2 = ( a->xmax - a->xmin ) / 2.0;
02212             double ydiff2 = ( a->ymax - a->ymin ) / 2.0;
02213             double xav    = ( a->xmax + a->xmin ) / 2.0;
02214             double yav    = ( a->ymax + a->ymin ) / 2.0;
02215 
02216             points_generate( xav - xdiff2 * zoom, xav + xdiff2 * zoom, yav - ydiff2 * zoom, yav + ydiff2 * zoom, nx, ny, &nout, &pout );
02217         }
02218     }
02219     else
02220     {
02221         points_read( fpoints, 2, &nout, &pout );
02222         if ( invariant )
02223             minell_scalepoints( me, nout, pout );
02224         else if ( square )
02225             points_scale( nout, pout, ks );
02226     }
02227 
02228     csa_approximate_points( a, nout, pout );
02229 
02230     if ( invariant )
02231         minell_rescalepoints( me, nout, pout );
02232     else if ( square )
02233         points_scale( nout, pout, 1.0 / ks );
02234 
02235     points_write( nout, pout );
02236 
02237     csa_destroy( a );
02238     free( pin );
02239     free( pout );
02240 
02241     return 0;
02242 }
02243 
02244 #endif                          // STANDALONE
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines