43 #define NPASTART 5 // Number of Points Allocated at Start 44 #define SVD_NMAX 30 // Maximal number of iterations in svd() 144 fprintf( stderr,
"error: csa: " );
145 va_start( args, format );
146 vfprintf( stderr, format, args );
158 static void*
alloc2d(
int n1,
int n2,
size_t unitsize )
167 assert( (
double) n1 * (
double) n2 <= (
double) UINT_MAX );
169 size = (size_t) ( n1 * n2 );
170 if ( ( p = calloc( size, unitsize ) ) == NULL )
171 csa_quit(
"alloc2d(): %s\n", strerror( errno ) );
173 assert( (
double) n2 * (
double)
sizeof (
void* ) <= (
double) UINT_MAX );
175 size = (size_t) n2 *
sizeof (
void* );
176 if ( ( pp = malloc( size ) ) == NULL )
177 csa_quit(
"alloc2d(): %s\n", strerror( errno ) );
178 for ( i = 0; i < n2; i++ )
179 pp[i] = &p[(
size_t) ( i * n1 ) * unitsize];
191 assert( pp != NULL );
192 p = ( (
void **) pp )[0];
204 t->
middle.
x = ( vertices[0].
x + vertices[1].
x + vertices[2].
x ) / 3.0;
205 t->
middle.
y = ( vertices[0].
y + vertices[1].
y + vertices[2].
y ) / 3.0;
260 bc[1] = ( dy - dx ) / t->
h;
261 bc[2] = -( dx + dy ) / t->
h;
263 else if ( t->
index == 1 )
265 bc[1] = ( dx + dy ) / t->
h;
266 bc[2] = ( dy - dx ) / t->
h;
268 else if ( t->
index == 2 )
270 bc[1] = ( dx - dy ) / t->
h;
271 bc[2] = ( dx + dy ) / t->
h;
275 bc[1] = -( dx + dy ) / t->
h;
276 bc[2] = ( dx - dy ) / t->
h;
278 bc[0] = 1.0 - bc[1] - bc[2];
286 double h = parent->
h;
301 for ( ii = 0; ii < 4; ++ii )
305 vertices[0].
x = xmin + h / 2.0;
306 vertices[0].
y = ymin + h / 2.0;
307 vertices[1].
x = xmin + h * ( ( ( ii + 1 ) % 4 ) / 2 );
308 vertices[1].
y = ymin + h * ( ( ( ii + 2 ) % 4 ) / 2 );
309 vertices[2].
x = xmin + h * ( ii / 2 );
310 vertices[2].
y = ymin + h * ( ( ( ii + 1 ) % 4 ) / 2 );
314 for ( ii = 0; ii < 25; ++ii )
324 for ( i = 0; i < 4; ++
i )
353 csa* a = malloc(
sizeof (
csa ) );
386 for ( j = 0; j < a->
nj; ++
j )
387 for ( i = 0; i < a->
ni; ++
i )
416 for ( i = 0; i < n; ++
i )
423 if ( p->
x < a->
xmin )
425 if ( p->
x > a->
xmax )
427 if ( p->
y < a->
ymin )
429 if ( p->
y > a->
ymax )
444 for ( j = 1; j < nj1; ++
j )
446 for ( i = 1; i < ni1; ++
i )
448 if ( squares[j][i]->
npoints > 0 )
450 if ( ( i + j ) % 2 == 0 )
453 squares[j - 1][i - 1]->
primary = 1;
454 squares[j + 1][i - 1]->
primary = 1;
455 squares[j - 1][i + 1]->
primary = 1;
456 squares[j + 1][i + 1]->
primary = 1;
474 int nps[7] = { 0, 0, 0, 0, 0, 0 };
484 fprintf( stderr,
"squarizing csa:\n" );
493 h = sqrt( dx * dy * a->
nppc / npoints );
500 a->
ni = (int) ( ceil( dx / h ) + 2 );
501 a->
nj = (int) ( ceil( dy / h ) + 2 );
505 fprintf( stderr,
" %d x %d squares\n", a->
ni, a->
nj );
513 for ( j = 0; j < a->
nj; ++
j )
514 for ( i = 0; i < a->
ni; ++
i )
520 for ( ii = 0; ii <
npoints; ++ii )
524 i = (int) ( floor( ( p->
x - a->
xmin ) / h ) + 1 );
525 j = (int) ( floor( ( p->
y - a->
ymin ) / h ) + 1 );
539 a->
pt = malloc( (
size_t) ( ( a->
ni / 2 + 1 ) * a->
nj ) * sizeof (
triangle* ) );
540 for ( j = 0, ii = 0, nadj = 0; j < a->
nj; ++
j )
542 for ( i = 0; i < a->
ni; ++
i )
568 fprintf( stderr,
" %d non-empty squares\n", ii );
569 fprintf( stderr,
" %d primary squares\n", a->
npt );
570 fprintf( stderr,
" %d primary squares with no data\n", nadj );
571 fprintf( stderr,
" %.2f points per square \n", (
double) a->
npoints / ii );
576 for ( i = 0; i < 6; ++
i )
577 fprintf( stderr,
" %d-%d points -- %d squares\n", i * 5, i * 5 + 4, nps[i] );
578 fprintf( stderr,
" %d or more points -- %d squares\n", i * 5, nps[i] );
583 fprintf( stderr,
" j\\i" );
584 for ( i = 0; i < a->
ni; ++
i )
585 fprintf( stderr,
"%3d ", i );
586 fprintf( stderr,
"\n" );
587 for ( j = a->
nj - 1; j >= 0; --j )
589 fprintf( stderr,
"%3d ", j );
590 for ( i = 0; i < a->
ni; ++
i )
595 fprintf( stderr,
"%3d ", s->
npoints );
597 fprintf( stderr,
" . " );
599 fprintf( stderr,
"\n" );
612 int imin = (int) floor( ( t->
middle.
x - t->
r - a->
xmin ) / t->
h );
613 int imax = (int) ceil( ( t->
middle.
x + t->
r - a->
xmin ) / t->
h );
614 int jmin = (int) floor( ( t->
middle.
y - t->
r - a->
ymin ) / t->
h );
615 int jmax = (int) ceil( ( t->
middle.
y + t->
r - a->
ymin ) / t->
h );
628 ( *squares ) = malloc( (
size_t) ( ( imax - imin + 1 ) * ( jmax - jmin + 1 ) ) *
sizeof (
square* ) );
630 for ( j = jmin; j <= jmax; ++
j )
632 for ( i = imin; i <= imax; ++
i )
647 return hypot( p1->
x - p2->
x, p1->
y - p2->
y );
661 int imax = (int) ceil( sqrt( (
double) ( npmax * 3 / 2 ) ) );
662 square *** squares =
alloc2d( imax, imax,
sizeof (
void* ) );
663 double h = t->
r * 2.0 / imax;
669 for ( j = 0; j < imax; ++
j )
670 for ( i = 0; i < imax; ++
i )
671 squares[j][i] =
square_create( a, xmin + h * i, ymin + h * j, i, j );
673 for ( ii = 0; ii < t->
npoints; ++ii )
677 i = (int) floor( ( p->
x - xmin ) / h );
678 j = (int) floor( ( p->
y - ymin ) / h );
688 pmiddle.
x = xmin + h * i + h2;
689 pmiddle.
y = ymin + h * j + h2;
697 for ( j = 0; j < imax; ++
j )
699 for ( i = 0; i < imax; ++
i )
703 if ( squares[j][i]->
npoints != 0 )
718 int npmin = a->
npmin;
719 int npmax = a->
npmax;
724 assert( a->
npt > 0 );
728 fprintf( stderr,
"pre-processing data points:\n " );
732 for ( i = 0; i < a->
npt; ++
i )
739 fprintf( stderr,
"." );
751 for ( ii = 0; ii < nsquares; ++ii )
756 for ( iii = 0; iii < s->
npoints; ++iii )
801 fprintf( stderr,
"\n %d sets enhanced, %d sets thinned\n", nincreased, nthinned );
828 static void svd(
double** a,
int n,
int m,
double* w,
double** v )
832 double tst1, c, f, g, h, s, scale;
834 assert( m > 0 && n > 0 );
836 rv1 = malloc( (
size_t) n *
sizeof (
double ) );
844 for ( i = 0; i < n; i++ )
853 for ( k = i; k < m; k++ )
854 scale += fabs( a[k][i] );
857 for ( k = i; k < m; k++ )
860 s += a[k][
i] * a[k][
i];
863 g = -copysign( sqrt( s ), f );
868 for ( j = l; j < n; j++ )
871 for ( k = i; k < m; k++ )
872 s += a[k][i] * a[k][j];
874 for ( k = i; k < m; k++ )
875 a[k][j] += f * a[k][i];
878 for ( k = i; k < m; k++ )
886 if ( i < m && i < n - 1 )
888 for ( k = l; k < n; k++ )
889 scale += fabs( a[i][k] );
892 for ( k = l; k < n; k++ )
895 s += a[
i][k] * a[
i][k];
898 g = -copysign( sqrt( s ), f );
901 for ( k = l; k < n; k++ )
902 rv1[k] = a[i][k] / h;
903 for ( j = l; j < m; j++ )
906 for ( k = l; k < n; k++ )
907 s += a[j][k] * a[i][k];
908 for ( k = l; k < n; k++ )
909 a[j][k] += s * rv1[k];
911 for ( k = l; k < n; k++ )
916 double tmp = fabs( w[i] ) + fabs( rv1[i] );
918 tst1 = ( tst1 > tmp ) ? tst1 : tmp;
925 for ( i = n - 1; i >= 0; i-- )
931 for ( j = l; j < n; j++ )
935 v[j][i] = ( a[i][j] / a[i][l] ) / g;
936 for ( j = l; j < n; j++ )
939 for ( k = l; k < n; k++ )
940 s += a[i][k] * v[k][j];
941 for ( k = l; k < n; k++ )
942 v[k][j] += s * v[k][i];
945 for ( j = l; j < n; j++ )
959 for ( i = ( m < n ) ? m - 1 : n - 1; i >= 0; i-- )
964 for ( j = l; j < n; j++ )
968 for ( j = l; j < n; j++ )
971 for ( k = l; k < m; k++ )
972 s += a[k][i] * a[k][j];
976 f = ( s / a[
i][
i] ) / g;
977 for ( k = i; k < m; k++ )
978 a[k][j] += f * a[k][i];
980 for ( j = i; j < m; j++ )
984 for ( j = i; j < m; j++ )
992 for ( k = n - 1; k >= 0; k-- )
999 int docancellation = 1;
1007 for ( l = k; l >= 0; l-- )
1009 double tst2 = fabs( rv1[l] ) + tst1;
1021 tst2 = fabs( w[l - 1] ) + tst1;
1028 if ( docancellation )
1032 for ( i = l; i <= k; i++ )
1035 rv1[
i] = c * rv1[
i];
1036 if ( ( fabs( f ) + tst1 ) == tst1 )
1044 for ( j = 0; j < m; j++ )
1049 a[
j][l1] = y * c + z * s;
1050 a[
j][
i] = z * c - y * s;
1069 f = 0.5 * ( ( ( g + z ) / h ) * ( ( g - z ) / y ) + y / h - h / y );
1070 g = hypot( f, 1.0 );
1071 f = x - ( z / x ) * z + ( h / x ) * ( y / ( f + copysign( g, f ) ) - h );
1077 for ( i1 = l; i1 < k; i1++ )
1092 for ( j = 0; j < n; j++ )
1096 v[
j][i1] = x * c + z * s;
1097 v[
j][
i] = z * c - x * s;
1111 for ( j = 0; j < m; j++ )
1115 a[
j][i1] = y * c + z * s;
1116 a[
j][
i] = z * c - y * s;
1131 for ( j = 0; j < n; j++ )
1144 static void lsq(
double** A,
int ni,
int nj,
double* z,
double* w,
double* sol )
1146 double** V =
alloc2d( ni, ni,
sizeof (
double ) );
1147 double** B =
alloc2d( nj, ni,
sizeof (
double ) );
1150 svd( A, ni, nj, w, V );
1152 for ( j = 0; j < ni; ++
j )
1153 for ( i = 0; i < ni; ++
i )
1155 for ( i = 0; i < ni; ++
i )
1159 for ( j = 0; j < nj; ++
j )
1162 double* b = &B[
i][
j];
1164 for ( ii = 0; ii < ni; ++ii )
1165 *b += v[ii] * a[ii];
1168 for ( i = 0; i < ni; ++
i )
1170 for ( i = 0; i < ni; ++
i )
1171 for ( j = 0; j < nj; ++
j )
1172 sol[i] += B[i][j] * z[j];
1197 int n[4] = { 0, 0, 0, 0 };
1201 fprintf( stderr,
"calculating spline coefficients for primary triangles:\n " );
1203 for ( i = 0; i < a->
npt; ++
i )
1208 double * z = malloc( (
size_t) npoints *
sizeof (
double ) );
1217 fprintf( stderr,
"." );
1221 for ( ii = 0; ii <
npoints; ++ii )
1222 z[ii] = points[ii]->z;
1236 double ** A =
alloc2d( 10, npoints,
sizeof (
double ) );
1239 for ( ii = 0; ii <
npoints; ++ii )
1241 double * aii = A[ii];
1250 tmp = bc[0] * bc[0];
1251 aii[0] = tmp * bc[0];
1253 aii[1] = tmp * bc[1];
1254 aii[2] = tmp * bc[2];
1255 tmp = bc[1] * bc[1];
1256 aii[6] = tmp * bc[1];
1258 aii[3] = tmp * bc[0];
1259 aii[7] = tmp * bc[2];
1260 tmp = bc[2] * bc[2];
1261 aii[9] = tmp * bc[2];
1263 aii[5] = tmp * bc[0];
1264 aii[8] = tmp * bc[1];
1265 aii[4] = bc[0] * bc[1] * bc[2] * 6.0;
1268 lsq( A, 10, npoints, z, w, b );
1272 for ( ii = 1; ii < 10; ++ii )
1276 else if ( w[ii] > wmax )
1279 if ( wmin < wmax / a->k )
1286 double ** A =
alloc2d( 6, npoints,
sizeof (
double ) );
1289 for ( ii = 0; ii <
npoints; ++ii )
1291 double* aii = A[ii];
1300 aii[0] = bc[0] * bc[0];
1301 aii[1] = bc[0] * bc[1] * 2.0;
1302 aii[2] = bc[0] * bc[2] * 2.0;
1303 aii[3] = bc[1] * bc[1];
1304 aii[4] = bc[1] * bc[2] * 2.0;
1305 aii[5] = bc[2] * bc[2];
1308 lsq( A, 6, npoints, z, w, b1 );
1312 for ( ii = 1; ii < 6; ++ii )
1316 else if ( w[ii] > wmax )
1319 if ( wmin < wmax / a->k )
1325 b[1] = ( b1[0] + 2.0 * b1[1] ) / 3.0;
1326 b[2] = ( b1[0] + 2.0 * b1[2] ) / 3.0;
1327 b[3] = ( b1[3] + 2.0 * b1[1] ) / 3.0;
1328 b[4] = ( b1[1] + b1[2] + b1[4] ) / 3.0;
1329 b[5] = ( b1[5] + 2.0 * b1[2] ) / 3.0;
1331 b[7] = ( b1[3] + 2.0 * b1[4] ) / 3.0;
1332 b[8] = ( b1[5] + 2.0 * b1[4] ) / 3.0;
1340 double ** A =
alloc2d( 3, npoints,
sizeof (
double ) );
1343 for ( ii = 0; ii <
npoints; ++ii )
1345 double* aii = A[ii];
1354 lsq( A, 3, npoints, z, w, b1 );
1358 for ( ii = 1; ii < 3; ++ii )
1362 else if ( w[ii] > wmax )
1365 if ( wmin < wmax / a->k )
1371 b[1] = ( 2.0 * b1[0] + b1[1] ) / 3.0;
1372 b[2] = ( 2.0 * b1[0] + b1[2] ) / 3.0;
1373 b[3] = ( 2.0 * b1[1] + b1[0] ) / 3.0;
1374 b[4] = ( b1[0] + b1[1] + b1[2] ) / 3.0;
1375 b[5] = ( 2.0 * b1[2] + b1[0] ) / 3.0;
1377 b[7] = ( 2.0 * b1[1] + b1[2] ) / 3.0;
1378 b[8] = ( 2.0 * b1[2] + b1[1] ) / 3.0;
1386 double ** A =
alloc2d( 1, npoints,
sizeof (
double ) );
1389 for ( ii = 0; ii <
npoints; ++ii )
1392 lsq( A, 1, npoints, z, w, b1 );
1434 fprintf( stderr,
"\n 3rd order -- %d sets\n", n[3] );
1435 fprintf( stderr,
" 2nd order -- %d sets\n", n[2] );
1436 fprintf( stderr,
" 1st order -- %d sets\n", n[1] );
1437 fprintf( stderr,
" 0th order -- %d sets\n", n[0] );
1445 fprintf( stderr,
" j\\i" );
1446 for ( i = 0; i < a->
ni; ++
i )
1447 fprintf( stderr,
"%2d ", i );
1448 fprintf( stderr,
"\n" );
1449 for ( j = a->
nj - 1; j >= 0; --j )
1451 fprintf( stderr,
"%2d ", j );
1452 for ( i = 0; i < a->
ni; ++
i )
1459 fprintf( stderr,
" . " );
1461 fprintf( stderr,
"\n" );
1478 fprintf( stderr,
"propagating spline coefficients to the remaining triangles:\n" );
1485 for ( ii = 0; ii < a->
npt; ++ii )
1492 double * cl = ( i > 0 ) ? squares[j][i - 1]->
coeffs : NULL;
1493 double * cb = ( j > 0 ) ? squares[j - 1][i]->
coeffs : NULL;
1494 double * cbl = ( i > 0 && j > 0 ) ? squares[j - 1][i - 1]->
coeffs : NULL;
1495 double * ca = ( j < nj - 1 ) ? squares[j + 1][i]->
coeffs : NULL;
1496 double * cal = ( j < nj - 1 && i > 0 ) ? squares[j + 1][i - 1]->
coeffs : NULL;
1498 c[7] = 2.0 * c[4] - c[1];
1499 c[11] = 2.0 * c[8] - c[5];
1500 c[15] = 2.0 * c[12] - c[9];
1502 c[10] = 2.0 * c[6] - c[2];
1503 c[13] = 2.0 * c[9] - c[5];
1504 c[16] = 2.0 * c[12] - c[8];
1506 c[19] = 2.0 * c[15] - c[11];
1515 cl[18] = c[0] + c[1] - c[4];
1516 cl[19] = c[1] + c[2] - c[5];
1517 cl[20] = c[2] + c[3] - c[6];
1519 cl[17] = 2.0 * cl[20] - cl[23];
1520 cl[14] = 2.0 * cl[18] - cl[22];
1528 cb[6] = c[0] + c[7] - c[4];
1529 cb[2] = 2.0 * cb[6] - cb[10];
1537 cbl[20] = cb[2] + cb[3] - cb[6];
1546 ca[4] = c[3] + c[10] - c[6];
1547 ca[1] = 2.0 * ca[4] - ca[7];
1555 cal[18] = ca[0] + ca[1] - ca[4];
1563 for ( ii = 0; ii < a->
npt; ++ii )
1570 double * cr = ( i < ni - 1 ) ? squares[j][i + 1]->
coeffs : NULL;
1571 double * car = ( i < ni - 1 && j < nj - 1 ) ? squares[j + 1][i + 1]->
coeffs : NULL;
1572 double * cbr = ( i < ni - 1 && j > 0 ) ? squares[j - 1][i + 1]->
coeffs : NULL;
1575 cr[13] = car[7] + car[14] - car[11];
1578 cr[11] = cbr[10] + cbr[17] - cbr[13];
1581 cr[5] = c[22] + c[23] - c[19];
1587 for ( ii = 0; ii < a->
npt; ++ii )
1593 double * cr = ( i < ni - 1 ) ? squares[j][i + 1]->
coeffs : NULL;
1597 cr[9] = ( cr[5] + cr[13] ) / 2.0;
1598 cr[8] = ( cr[5] + cr[11] ) / 2.0;
1599 cr[15] = ( cr[11] + cr[19] ) / 2.0;
1600 cr[16] = ( cr[13] + cr[19] ) / 2.0;
1601 cr[12] = ( cr[8] + cr[16] ) / 2.0;
1607 fprintf( stderr,
"checking that all coefficients have been set:\n" );
1611 for ( ii = 0; ii < ni * nj; ++ii )
1613 square* s = squares[0][ii];
1619 for ( i = 0; i < 25; ++
i )
1620 if (
isnan( c[i] ) )
1621 fprintf( stderr,
" squares[%d][%d]->coeffs[%d] = NaN\n", s->
j, s->
i, i );
1625 static int i300[] = { 12, 12, 12, 12 };
1626 static int i030[] = { 3, 24, 21, 0 };
1627 static int i003[] = { 0, 3, 24, 21 };
1628 static int i210[] = { 9, 16, 15, 8 };
1629 static int i021[] = { 2, 17, 22, 7 };
1630 static int i102[] = { 4, 6, 20, 18 };
1631 static int i120[] = { 6, 20, 18, 4 };
1632 static int i012[] = { 1, 10, 23, 14 };
1633 static int i201[] = { 8, 9, 16, 15 };
1634 static int i111[] = { 5, 13, 19, 11 };
1642 for ( j = 0; j < a->
nj; ++
j )
1644 for ( i = 0; i < a->
ni; ++
i )
1650 for ( ii = 0; ii < 4; ++ii )
1655 for ( cc = 0; cc < 10; ++cc )
1677 double ii = ( p->
x - a->
xmin ) / h + 1.0;
1678 double jj = ( p->
y - a->
ymin ) / h + 1.0;
1686 if ( ii < 0.0 || jj < 0.0 || ii > (
double) a->
ni - 1.0 || jj > (double) a->
nj - 1.0 )
1692 i = (int) floor( ii );
1693 j = (int) floor( jj );
1700 if ( fi + fj < 1.0 )
1707 if ( fi + fj < 1.0 )
1726 double tmp1 = bc1 * bc1;
1727 double tmp2 = bc2 * bc2;
1728 double tmp3 = bc3 * bc3;
1733 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;
1736 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;
1739 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;
1742 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;
1751 for ( ii = 0; ii < n; ++ii )
1772 a->
nppc = (int) nppc;
1775 #if defined ( STANDALONE ) 1780 #define BUFSIZE 10240 1781 #define STRBUFSIZE 64 1783 static void points_generate(
double xmin,
double xmax,
double ymin,
double ymax,
int nx,
int ny,
int* nout,
point** pout )
1785 double stepx, stepy;
1789 if ( nx < 1 || ny < 1 )
1797 *pout = malloc( *nout *
sizeof (
point ) );
1799 stepx = ( nx > 1 ) ? ( xmax - xmin ) / ( nx - 1 ) : 0.0;
1800 stepy = ( ny > 1 ) ? ( ymax - ymin ) / ( ny - 1 ) : 0.0;
1801 x0 = ( nx > 1 ) ? xmin : ( xmin + xmax ) / 2.0;
1802 yy = ( ny > 1 ) ? ymin : ( ymin + ymax ) / 2.0;
1805 for ( j = 0; j < ny; ++
j )
1808 for ( i = 0; i < nx; ++
i )
1810 point* p = &( *pout )[ii];
1825 if ( token == NULL )
1831 *value = strtod( token, &end );
1842 #define NALLOCATED_START 1024 1856 char seps[] =
" ,;\t";
1859 if ( dim < 2 || dim > 3 )
1866 if ( fname == NULL )
1870 if ( strcmp( fname,
"stdin" ) == 0 || strcmp( fname,
"-" ) == 0 )
1874 f = fopen( fname,
"r" );
1876 csa_quit(
"%s: %s\n", fname, strerror( errno ) );
1880 *points = malloc( nallocated *
sizeof (
point ) );
1882 while ( fgets( buf,
BUFSIZE, f ) != NULL )
1886 if ( *n == nallocated )
1889 *points = realloc( *points, nallocated *
sizeof (
point ) );
1892 p = &( *points )[*n];
1894 if ( buf[0] ==
'#' )
1896 if ( ( token = strtok( buf, seps ) ) == NULL )
1900 if ( ( token = strtok( NULL, seps ) ) == NULL )
1908 if ( ( token = strtok( NULL, seps ) ) == NULL )
1922 *points = realloc( *points, *n *
sizeof (
point ) );
1925 if ( fclose( f ) != 0 )
1926 csa_quit(
"%s: %s\n", fname, strerror( errno ) );
1929 static void points_write(
int n,
point* points )
1934 printf(
"Output:\n" );
1936 for ( i = 0; i < n; ++
i )
1940 printf(
"%.15g %.15g %.15g\n", p->
x, p->
y, p->
z );
1946 double xmin, ymin, xmax, ymax;
1953 xmin = xmax = points[0].
x;
1954 ymin = ymax = points[0].
y;
1956 for ( i = 1; i < n; ++
i )
1962 else if ( p->
x > xmax )
1966 else if ( p->
y > ymax )
1970 if ( xmin == xmax || ymin == ymax )
1973 k = ( ymax - ymin ) / ( xmax - xmin );
1975 for ( i = 0; i < n; ++
i )
1985 for ( i = 0; i < n; ++
i )
1991 printf(
"Usage: csabathy -i <XYZ file> {-o <XY file>|-n <nx>x<ny> [-c|-s] [-z <zoom>]}\n" );
1992 printf(
" [-v|-V] [-P nppc=<value>] [-P k=<value>]\n" );
1993 printf(
"Options:\n" );
1994 printf(
" -c -- scale internally so that the minimal ellipse turns into a\n" );
1995 printf(
" circle (this produces results invariant to affine\n" );
1996 printf(
" transformations)\n" );
1997 printf(
" -i <XYZ file> -- three-column file with points to approximate from (use\n" );
1998 printf(
" \"-i stdin\" or \"-i -\" for standard input)\n" );
1999 printf(
" -n <nx>x<ny> -- generate <nx>x<ny> output rectangular grid\n" );
2000 printf(
" -o <XY file> -- two-column file with points to approximate in (use\n" );
2001 printf(
" \"-o stdin\" or \"-o -\" for standard input)\n" );
2002 printf(
" -s -- scale internally so that Xmax - Xmin = Ymax - Ymin\n" );
2003 printf(
" -v -- verbose / version\n" );
2004 printf(
" -z <zoom> -- zoom in (if <zoom> < 1) or out (<zoom> > 1)\n" );
2005 printf(
" -P nppc=<value> -- set the average number of points per cell (default = 5\n" );
2006 printf(
" increase if the point distribution is strongly non-uniform\n" );
2007 printf(
" to get larger cells)\n" );
2008 printf(
" -P k=<value> -- set the spline sensitivity (default = 140, reduce to get\n" );
2009 printf(
" smoother results)\n" );
2010 printf(
" -V -- very verbose / version\n" );
2011 printf(
"Description:\n" );
2012 printf(
" `csabathy' approximates irregular scalar 2D data in specified points using\n" );
2013 printf(
" C1-continuous bivariate cubic spline. The calculated values are written to\n" );
2014 printf(
" standard output.\n" );
2019 static void version()
2025 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 )
2035 if (
argv[i][0] !=
'-' )
2038 switch (
argv[i][1] )
2049 csa_quit(
"no file name found after -i\n" );
2056 *generate_points = 1;
2058 csa_quit(
"no grid dimensions found after -n\n" );
2059 if ( sscanf(
argv[i],
"%dx%d", nx, ny ) != 2 )
2060 csa_quit(
"could not read grid dimensions after \"-n\"\n" );
2061 if ( *nx <= 0 || *nx > NIMAX || *ny <= 0 || *ny > NIMAX )
2062 csa_quit(
"invalid size for output grid\n" );
2068 csa_quit(
"no file name found after -o\n" );
2085 csa_quit(
"no zoom value found after -z\n" );
2086 *zoom = atof(
argv[i] );
2091 char prmstr[STRBUFSIZE] =
"";
2096 csa_quit(
"no input found after -P\n" );
2098 if ( strlen(
argv[i] ) >= STRBUFSIZE )
2099 csa_quit(
"could not interpret \"%s\" after -P option\n",
argv[i] );
2101 strcpy( prmstr,
argv[i] );
2102 token = strtok( prmstr, delim );
2103 if ( token == NULL )
2104 csa_quit(
"could not interpret \"%s\" after -P option\n",
argv[i] );
2106 if ( strcmp( token,
"nppc" ) == 0 )
2110 token = strtok( NULL, delim );
2111 if ( token == NULL )
2112 csa_quit(
"could not interpret \"%s\" after -P option\n",
argv[i] );
2114 n = strtol( token, NULL, 10 );
2115 if ( n == LONG_MIN || n == LONG_MAX )
2116 csa_quit(
"could not interpret \"%s\" after -P option\n",
argv[i] );
2118 csa_quit(
"non-sensible value for \"nppc\" parameter\n" );
2121 else if ( strcmp( token,
"k" ) == 0 )
2125 token = strtok( NULL, delim );
2126 if ( token == NULL )
2127 csa_quit(
"could not interpret \"%s\" after -P option\n",
argv[i] );
2129 n = strtol( token, NULL, 10 );
2130 if ( n == LONG_MIN || n == LONG_MAX )
2131 csa_quit(
"could not interpret \"%s\" after -P option\n",
argv[i] );
2133 csa_quit(
"non-sensible value for \"k\" parameter\n" );
2158 char * fdata = NULL;
2159 char * fpoints = NULL;
2166 int generate_points = 0;
2167 point * pout = NULL;
2176 parse_commandline(
argc, argv, &fdata, &fpoints, &invariant, &square, &generate_points, &nx, &ny, &nppc, &k, &zoom );
2178 if ( fdata == NULL )
2181 if ( !generate_points && fpoints == NULL )
2182 csa_quit(
"no output grid specified\n" );
2191 me = minell_build( nin, pin );
2192 minell_scalepoints( me, nin, pin );
2205 if ( generate_points )
2207 if (
isnan( zoom ) )
2208 points_generate( a->
xmin - a->
h, a->
xmax + a->
h, a->
ymin - a->
h, a->
ymax + a->
h, nx, ny, &nout, &pout );
2211 double xdiff2 = ( a->
xmax - a->
xmin ) / 2.0;
2212 double ydiff2 = ( a->
ymax - a->
ymin ) / 2.0;
2213 double xav = ( a->
xmax + a->
xmin ) / 2.0;
2214 double yav = ( a->
ymax + a->
ymin ) / 2.0;
2216 points_generate( xav - xdiff2 * zoom, xav + xdiff2 * zoom, yav - ydiff2 * zoom, yav + ydiff2 * zoom, nx, ny, &nout, &pout );
2223 minell_scalepoints( me, nout, pout );
2231 minell_rescalepoints( me, nout, pout );
2235 points_write( nout, pout );
2244 #endif // STANDALONE
void csa_approximate_point(csa *a, point *p)
static void triangle_addpoint(triangle *t, point *p)
static void csa_attachpoints(csa *a)
static square * square_create(csa *parent, double xmin, double ymin, int i, int j)
void points_scale(int n, point *points, double k)
void csa_setk(csa *a, int k)
static PLCHAR_VECTOR usage
static void thindata(triangle *t, int npmax)
static void getsquares(csa *a, triangle *t, int *n, square ***squares)
void csa_calculatespline(csa *a)
static double distance(point *p1, point *p2)
static void csa_findprimarycoeffs(csa *a)
static void triangle_calculatebc(triangle *t, point *p, double bc[])
void points_read(char *fname, int dim, int *n, point **points)
void csa_addpoints(csa *a, int n, point points[])
static int str2double(char *token, double *value)
static void free2d(void *pp)
static void triangle_destroy(triangle *t)
static void csa_quit(const char *format,...)
static void lsq(double **A, int ni, int nj, double *z, double *w, double *sol)
static triangle * triangle_create(square *s, point vertices[], int index)
static void svd(double **a, int n, int m, double *w, double **v)
static PLFLT value(double n1, double n2, double hue)
void csa_setnpmin(csa *a, int npmin)
static void csa_setprimaryflag(csa *a)
static void csa_sethascoeffsflag(csa *a)
static void csa_findsecondarycoeffs(csa *a)
static void * alloc2d(int n1, int n2, size_t unitsize)
static void square_addpoint(square *s, point *p)
void csa_setnpmax(csa *a, int npmax)
static void csa_squarize(csa *a)
void csa_setnppc(csa *a, double nppc)
double points_scaletosquare(int n, point *points)
void csa_approximate_points(csa *a, int n, point *points)
static void square_destroy(square *s)