PLplot
5.10.0
|
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