PLplot  5.10.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
csa.c
Go to the documentation of this file.
1 //--------------------------------------------------------------------------
2 //
3 // File: csa.c
4 //
5 // Created: 16/10/2002
6 //
7 // Author: Pavel Sakov
8 // CSIRO Marine Research
9 //
10 // Purpose: 2D data approximation with bivariate C1 cubic spline.
11 // A set of library functions + standalone utility.
12 //
13 // Description: See J. Haber, F. Zeilfelder, O.Davydov and H.-P. Seidel,
14 // Smooth approximation and rendering of large scattered data
15 // sets, in ``Proceedings of IEEE Visualization 2001''
16 // (Th.Ertl, K.Joy and A.Varshney, Eds.), pp.341-347, 571,
17 // IEEE Computer Society, 2001.
18 // http://www.uni-giessen.de/www-Numerische-Mathematik/
19 // davydov/VIS2001.ps.gz
20 // http://www.math.uni-mannheim.de/~lsmath4/paper/
21 // VIS2001.pdf.gz
22 //
23 // Revisions: 09/04/2003 PS: Modified points_read() to read from a
24 // file specified by name, not by handle.
25 //
26 //--------------------------------------------------------------------------
27 
28 #include <stdlib.h>
29 #include <stdio.h>
30 #include <stdarg.h>
31 #include <limits.h>
32 #include <float.h>
33 #include <math.h>
34 #include <assert.h>
35 #include <string.h>
36 #include <errno.h>
37 #include "version.h"
38 #include "nan.h"
39 #include "csa.h"
40 
41 int csa_verbose = 0;
42 
43 #define NPASTART 5 // Number of Points Allocated at Start
44 #define SVD_NMAX 30 // Maximal number of iterations in svd()
45 
46 // default algorithm parameters
47 #define NPMIN_DEF 3
48 #define NPMAX_DEF 40
49 #define K_DEF 140
50 #define NPPC_DEF 5
51 
52 struct square;
53 typedef struct square square;
54 
55 typedef struct
56 {
58  int index; // index within parent square; 0 <= index <=
59  // 3
60  point vertices[3];
61  point middle; // barycenter
62  double h; // parent square edge length
63  double r; // data visibility radius
64 
65  //
66  // points used -- in primary triangles only
67  //
69  int npoints;
71 
72  int primary; // flag -- whether calculate spline
73  // coefficients directly (by least squares
74  // method) (primary = 1) or indirectly (from
75  // C1 smoothness conditions) (primary = 0)
76  int hascoeffs; // flag -- whether there are no NaNs among
77  // the spline coefficients
78  int order; // spline order -- for primary triangles only
79  //
80 } triangle;
81 
82 struct square
83 {
85  int i, j; // indices
86 
88  int npoints;
90 
91  int primary; // flag -- whether this square contains a
92  // primary triangle
94 
95  double coeffs[25];
96 };
97 
98 struct csa
99 {
100  double xmin;
101  double xmax;
102  double ymin;
103  double ymax;
104 
106  int npoints;
108 
109  //
110  // squarization
111  //
112  int ni;
113  int nj;
114  double h;
115  square *** squares; // square* [j][i]
116 
117  int npt; // Number of Primary Triangles
118  triangle** pt; // Primary Triangles -- triangle* [npt]
119 
120  //
121  // algorithm parameters
122  //
123  int npmin; // minimal number of points locally involved
124  // in * spline calculation (normally = 3)
125  int npmax; // maximal number of points locally involved
126  // in * spline calculation (required > 10, *
127  // recommended 20 < npmax < 60)
128  double k; // relative tolerance multiple in fitting
129  // spline coefficients: the higher this
130  // value, the higher degree of the locally
131  // fitted spline (recommended 80 < k < 200)
132  int nppc; // average number of points per square
133 };
134 
135 void csa_setnppc( csa* a, double nppc );
136 
137 static void csa_quit( const char* format, ... )
138 {
139  va_list args;
140 
141  fflush( stdout ); // just in case -- to have the exit message
142  // last
143 
144  fprintf( stderr, "error: csa: " );
145  va_start( args, format );
146  vfprintf( stderr, format, args );
147  va_end( args );
148 
149  exit( 1 );
150 }
151 
152 // Allocates n1xn2 matrix of something. Note that it will be accessed as
153 // [n2][n1].
154 // @param n1 Number of columns
155 // @param n2 Number of rows
156 // @return Matrix
157 //
158 static void* alloc2d( int n1, int n2, size_t unitsize )
159 {
160  size_t size;
161  char * p;
162  char ** pp;
163  int i;
164 
165  assert( n1 > 0 );
166  assert( n2 > 0 );
167  assert( (double) n1 * (double) n2 <= (double) UINT_MAX );
168 
169  size = (size_t) ( n1 * n2 );
170  if ( ( p = calloc( size, unitsize ) ) == NULL )
171  csa_quit( "alloc2d(): %s\n", strerror( errno ) );
172 
173  assert( (double) n2 * (double) sizeof ( void* ) <= (double) UINT_MAX );
174 
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];
180 
181  return pp;
182 }
183 
184 // Destroys n1xn2 matrix.
185 // @param pp Matrix
186 //
187 static void free2d( void* pp )
188 {
189  void* p;
190 
191  assert( pp != NULL );
192  p = ( (void **) pp )[0];
193  free( pp );
194  assert( p != NULL );
195  free( p );
196 }
197 
198 static triangle* triangle_create( square* s, point vertices[], int index )
199 {
200  triangle* t = malloc( sizeof ( triangle ) );
201 
202  t->parent = s;
203  memcpy( t->vertices, vertices, sizeof ( point ) * 3 );
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;
206  t->h = s->parent->h;
207  t->index = index;
208 
209  t->r = 0.0;
210  t->points = 0;
211  t->nallocated = 0;
212  t->npoints = 0;
213  t->hascoeffs = 0;
214  t->primary = 0;
215  t->order = -1;
216 
217  return t;
218 }
219 
220 static void triangle_addpoint( triangle* t, point* p )
221 {
222  if ( t->nallocated == t->npoints )
223  {
224  if ( t->nallocated == 0 )
225  {
226  t->points = malloc( NPASTART * sizeof ( point* ) );
227  t->nallocated = NPASTART;
228  }
229  else
230  {
231  t->points = realloc( t->points, (size_t) ( t->nallocated * 2 ) * sizeof ( point* ) );
232  t->nallocated *= 2;
233  }
234  }
235 
236  t->points[t->npoints] = p;
237  t->npoints++;
238 }
239 
240 static void triangle_destroy( triangle* t )
241 {
242  if ( t->points != NULL )
243  free( t->points );
244  free( t );
245 }
246 
247 // Calculates barycentric coordinates of a point.
248 // Takes into account that possible triangles are rectangular, with the right
249 // angle at t->vertices[0], the vertices[1] vertex being in
250 // (-3*PI/4) + (PI/2) * t->index direction from vertices[0], and
251 // vertices[2] being at (-5*PI/4) + (PI/2) * t->index.
252 //
253 static void triangle_calculatebc( triangle* t, point* p, double bc[] )
254 {
255  double dx = p->x - t->vertices[0].x;
256  double dy = p->y - t->vertices[0].y;
257 
258  if ( t->index == 0 )
259  {
260  bc[1] = ( dy - dx ) / t->h;
261  bc[2] = -( dx + dy ) / t->h;
262  }
263  else if ( t->index == 1 )
264  {
265  bc[1] = ( dx + dy ) / t->h;
266  bc[2] = ( dy - dx ) / t->h;
267  }
268  else if ( t->index == 2 )
269  {
270  bc[1] = ( dx - dy ) / t->h;
271  bc[2] = ( dx + dy ) / t->h;
272  }
273  else
274  {
275  bc[1] = -( dx + dy ) / t->h;
276  bc[2] = ( dx - dy ) / t->h;
277  }
278  bc[0] = 1.0 - bc[1] - bc[2];
279 }
280 
281 static square* square_create( csa* parent, double xmin, double ymin, int i, int j )
282 {
283  int ii;
284 
285  square * s = malloc( sizeof ( square ) );
286  double h = parent->h;
287 
288  s->parent = parent;
289  s->i = i;
290  s->j = j;
291 
292  s->points = NULL;
293  s->nallocated = 0;
294  s->npoints = 0;
295 
296  s->primary = 0;
297 
298  //
299  // create 4 triangles formed by diagonals
300  //
301  for ( ii = 0; ii < 4; ++ii )
302  {
303  point vertices[3];
304 
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 ); // 0 1 1 0
308  vertices[1].y = ymin + h * ( ( ( ii + 2 ) % 4 ) / 2 ); // 1 1 0 0
309  vertices[2].x = xmin + h * ( ii / 2 ); // 0 0 1 1
310  vertices[2].y = ymin + h * ( ( ( ii + 1 ) % 4 ) / 2 ); // 0 1 1 0
311  s->triangles[ii] = triangle_create( s, vertices, ii );
312  }
313 
314  for ( ii = 0; ii < 25; ++ii )
315  s->coeffs[ii] = NaN;
316 
317  return s;
318 }
319 
320 static void square_destroy( square* s )
321 {
322  int i;
323 
324  for ( i = 0; i < 4; ++i )
325  triangle_destroy( s->triangles[i] );
326  if ( s->points != NULL )
327  free( s->points );
328  free( s );
329 }
330 
331 static void square_addpoint( square* s, point* p )
332 {
333  if ( s->nallocated == s->npoints )
334  {
335  if ( s->nallocated == 0 )
336  {
337  s->points = malloc( NPASTART * sizeof ( point* ) );
338  s->nallocated = NPASTART;
339  }
340  else
341  {
342  s->points = realloc( s->points, (size_t) ( s->nallocated * 2 ) * sizeof ( point* ) );
343  s->nallocated *= 2;
344  }
345  }
346 
347  s->points[s->npoints] = p;
348  s->npoints++;
349 }
350 
352 {
353  csa* a = malloc( sizeof ( csa ) );
354 
355  a->xmin = DBL_MAX;
356  a->xmax = -DBL_MAX;
357  a->ymin = DBL_MAX;
358  a->ymax = -DBL_MAX;
359 
360  a->points = malloc( NPASTART * sizeof ( point* ) );
361  a->nallocated = NPASTART;
362  a->npoints = 0;
363 
364  a->ni = 0;
365  a->nj = 0;
366  a->h = NaN;
367  a->squares = NULL;
368 
369  a->npt = 0;
370  a->pt = NULL;
371 
372  a->npmin = NPMIN_DEF;
373  a->npmax = NPMAX_DEF;
374  a->k = K_DEF;
375  a->nppc = NPPC_DEF;
376 
377  return a;
378 }
379 
380 void csa_destroy( csa* a )
381 {
382  int i, j;
383 
384  if ( a->squares != NULL )
385  {
386  for ( j = 0; j < a->nj; ++j )
387  for ( i = 0; i < a->ni; ++i )
388  square_destroy( a->squares[j][i] );
389  free2d( a->squares );
390  }
391  if ( a->pt != NULL )
392  free( a->pt );
393  if ( a->points != NULL )
394  free( a->points );
395  free( a );
396 }
397 
398 void csa_addpoints( csa* a, int n, point points[] )
399 {
400  int na = a->nallocated;
401  int i;
402 
403  assert( a->squares == NULL );
404  //
405  // (can be called prior to squarization only)
406  //
407 
408  while ( na < a->npoints + n )
409  na *= 2;
410  if ( na != a->nallocated )
411  {
412  a->points = realloc( a->points, (size_t) na * sizeof ( point* ) );
413  a->nallocated = na;
414  }
415 
416  for ( i = 0; i < n; ++i )
417  {
418  point* p = &points[i];
419 
420  a->points[a->npoints] = p;
421  a->npoints++;
422 
423  if ( p->x < a->xmin )
424  a->xmin = p->x;
425  if ( p->x > a->xmax )
426  a->xmax = p->x;
427  if ( p->y < a->ymin )
428  a->ymin = p->y;
429  if ( p->y > a->ymax )
430  a->ymax = p->y;
431  }
432 }
433 
434 // Marks the squares containing "primary" triangles by setting "primary" flag
435 // to 1.
436 //
437 static void csa_setprimaryflag( csa* a )
438 {
439  square*** squares = a->squares;
440  int nj1 = a->nj - 1;
441  int ni1 = a->ni - 1;
442  int i, j;
443 
444  for ( j = 1; j < nj1; ++j )
445  {
446  for ( i = 1; i < ni1; ++i )
447  {
448  if ( squares[j][i]->npoints > 0 )
449  {
450  if ( ( i + j ) % 2 == 0 )
451  {
452  squares[j][i]->primary = 1;
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;
457  }
458  else
459  {
460  squares[j - 1][i]->primary = 1;
461  squares[j + 1][i]->primary = 1;
462  squares[j][i - 1]->primary = 1;
463  squares[j][i + 1]->primary = 1;
464  }
465  }
466  }
467  }
468 }
469 
470 // Splits the data domain in a number of squares.
471 //
472 static void csa_squarize( csa* a )
473 {
474  int nps[7] = { 0, 0, 0, 0, 0, 0 }; // stats on number of points per
475  // square
476  double dx = a->xmax - a->xmin;
477  double dy = a->ymax - a->ymin;
478  int npoints = a->npoints;
479  double h;
480  int i, j, ii, nadj;
481 
482  if ( csa_verbose )
483  {
484  fprintf( stderr, "squarizing csa:\n" );
485  fflush( stderr );
486  }
487 
488  assert( a->squares == NULL );
489  //
490  // (can be done only once)
491  //
492 
493  h = sqrt( dx * dy * a->nppc / npoints ); // square edge size
494  if ( dx < h )
495  h = dy * a->nppc / npoints;
496  if ( dy < h )
497  h = dx * a->nppc / npoints;
498  a->h = h;
499 
500  a->ni = (int) ( ceil( dx / h ) + 2 );
501  a->nj = (int) ( ceil( dy / h ) + 2 );
502 
503  if ( csa_verbose )
504  {
505  fprintf( stderr, " %d x %d squares\n", a->ni, a->nj );
506  fflush( stderr );
507  }
508 
509  //
510  // create squares
511  //
512  a->squares = alloc2d( a->ni, a->nj, sizeof ( void* ) );
513  for ( j = 0; j < a->nj; ++j )
514  for ( i = 0; i < a->ni; ++i )
515  a->squares[j][i] = square_create( a, a->xmin + h * ( i - 1 ), a->ymin + h * ( j - 1 ), i, j );
516 
517  //
518  // map points to squares
519  //
520  for ( ii = 0; ii < npoints; ++ii )
521  {
522  point* p = a->points[ii];
523 
524  i = (int) ( floor( ( p->x - a->xmin ) / h ) + 1 );
525  j = (int) ( floor( ( p->y - a->ymin ) / h ) + 1 );
526  square_addpoint( a->squares[j][i], p );
527  }
528 
529  //
530  // mark relevant squares with no points
531  //
532  csa_setprimaryflag( a );
533 
534  //
535  // Create a list of "primary" triangles, for which spline coefficients
536  // will be calculated directy (by least squares method), without using
537  // C1 smoothness conditions.
538  //
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 )
541  {
542  for ( i = 0; i < a->ni; ++i )
543  {
544  square* s = a->squares[j][i];
545 
546  if ( s->npoints > 0 )
547  {
548  int nn = s->npoints / 5;
549 
550  if ( nn > 6 )
551  nn = 6;
552  nps[nn]++;
553  ii++;
554  }
555  if ( s->primary && s->npoints == 0 )
556  nadj++;
557  if ( s->primary )
558  {
559  a->pt[a->npt] = s->triangles[0];
560  s->triangles[0]->primary = 1;
561  a->npt++;
562  }
563  }
564  }
565 
566  if ( csa_verbose )
567  {
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 );
572  }
573 
574  if ( csa_verbose == 2 )
575  {
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] );
579  }
580 
581  if ( csa_verbose == 2 )
582  {
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 )
588  {
589  fprintf( stderr, "%3d ", j );
590  for ( i = 0; i < a->ni; ++i )
591  {
592  square* s = a->squares[j][i];
593 
594  if ( s->npoints > 0 )
595  fprintf( stderr, "%3d ", s->npoints );
596  else
597  fprintf( stderr, " . " );
598  }
599  fprintf( stderr, "\n" );
600  }
601  }
602 
603  if ( csa_verbose )
604  fflush( stderr );
605 }
606 
607 // Returns all squares intersecting with a square with center in t->middle
608 // and edges of length 2*t->r parallel to X and Y axes.
609 //
610 static void getsquares( csa* a, triangle* t, int* n, square*** squares )
611 {
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 );
616  int i, j;
617 
618  if ( imin < 0 )
619  imin = 0;
620  if ( imax >= a->ni )
621  imax = a->ni - 1;
622  if ( jmin < 0 )
623  jmin = 0;
624  if ( jmax >= a->nj )
625  jmax = a->nj - 1;
626 
627  *n = 0;
628  ( *squares ) = malloc( (size_t) ( ( imax - imin + 1 ) * ( jmax - jmin + 1 ) ) * sizeof ( square* ) );
629 
630  for ( j = jmin; j <= jmax; ++j )
631  {
632  for ( i = imin; i <= imax; ++i )
633  {
634  square* s = a->squares[j][i];
635 
636  if ( s->npoints > 0 )
637  {
638  ( *squares )[*n] = a->squares[j][i];
639  ( *n )++;
640  }
641  }
642  }
643 }
644 
645 static double distance( point* p1, point* p2 )
646 {
647  return hypot( p1->x - p2->x, p1->y - p2->y );
648 }
649 
650 // Thins data by creating an auxiliary regular grid and for leaving only
651 // the most central point within each grid cell.
652 // (I follow the paper here. It is possible that taking average -- in terms of
653 // both value and position -- of all points within a cell would be a bit more
654 // robust. However, because of keeping only shallow copies of input points,
655 // this would require quite a bit of structural changes. So, leaving it as is
656 // for now.)
657 //
658 static void thindata( triangle* t, int npmax )
659 {
660  csa * a = t->parent->parent;
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;
664  double h2 = h / 2.0;
665  double xmin = t->middle.x - t->r;
666  double ymin = t->middle.y - t->r;
667  int i, j, ii;
668 
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 );
672 
673  for ( ii = 0; ii < t->npoints; ++ii )
674  {
675  square* s;
676  point * p = t->points[ii];
677  i = (int) floor( ( p->x - xmin ) / h );
678  j = (int) floor( ( p->y - ymin ) / h );
679  s = squares[j][i];
680 
681  if ( s->npoints == 0 )
682  square_addpoint( s, p );
683  else // npoints == 1
684 
685  {
686  point pmiddle;
687 
688  pmiddle.x = xmin + h * i + h2;
689  pmiddle.y = ymin + h * j + h2;
690 
691  if ( distance( s->points[0], &pmiddle ) > distance( p, &pmiddle ) )
692  s->points[0] = p;
693  }
694  }
695 
696  t->npoints = 0;
697  for ( j = 0; j < imax; ++j )
698  {
699  for ( i = 0; i < imax; ++i )
700  {
701  square* s = squares[j][i];
702 
703  if ( squares[j][i]->npoints != 0 )
704  triangle_addpoint( t, s->points[0] );
705  square_destroy( s );
706  }
707  }
708 
709  free2d( squares );
710  imax++;
711 }
712 
713 // Finds data points to be used in calculating spline coefficients for each
714 // primary triangle.
715 //
716 static void csa_attachpoints( csa* a )
717 {
718  int npmin = a->npmin;
719  int npmax = a->npmax;
720  int nincreased = 0;
721  int nthinned = 0;
722  int i;
723 
724  assert( a->npt > 0 );
725 
726  if ( csa_verbose )
727  {
728  fprintf( stderr, "pre-processing data points:\n " );
729  fflush( stderr );
730  }
731 
732  for ( i = 0; i < a->npt; ++i )
733  {
734  triangle* t = a->pt[i];
735  int increased = 0;
736 
737  if ( csa_verbose )
738  {
739  fprintf( stderr, "." );
740  fflush( stderr );
741  }
742 
743  t->r = t->h * 1.25;
744  while ( 1 )
745  {
746  int nsquares = 0;
747  square** squares = NULL;
748  int ii;
749 
750  getsquares( a, t, &nsquares, &squares );
751  for ( ii = 0; ii < nsquares; ++ii )
752  {
753  square* s = squares[ii];
754  int iii;
755 
756  for ( iii = 0; iii < s->npoints; ++iii )
757  {
758  point* p = s->points[iii];
759 
760  if ( distance( p, &t->middle ) <= t->r )
761  triangle_addpoint( t, p );
762  }
763  }
764 
765  free( squares );
766 
767  if ( t->npoints < npmin )
768  {
769  if ( !increased )
770  {
771  increased = 1;
772  nincreased++;
773  }
774  t->r *= 1.25;
775  t->npoints = 0;
776  }
777  else if ( t->npoints > npmax )
778  {
779  nthinned++;
780  thindata( t, npmax );
781  if ( t->npoints > npmin )
782  break;
783  else
784  {
785  //
786  // Sometimes you have too much data, you thin it and --
787  // oops -- you have too little. This is not a frequent
788  // event, so let us not bother to put a new subdivision.
789  //
790  t->r *= 1.25;
791  t->npoints = 0;
792  }
793  }
794  else
795  break;
796  }
797  }
798 
799  if ( csa_verbose )
800  {
801  fprintf( stderr, "\n %d sets enhanced, %d sets thinned\n", nincreased, nthinned );
802  fflush( stderr );
803  }
804 }
805 
806 static int n2q( int n )
807 {
808  assert( n >= 3 );
809 
810  if ( n >= 10 )
811  return 3;
812  else if ( n >= 6 )
813  return 2;
814  else // n == 3
815  return 1;
816 }
817 
818 //* Singular value decomposition.
819 // Borrowed from EISPACK (1972-1973).
820 // Presents input matrix A as A = U.W.V'.
821 //
822 // @param a Input matrix A = U.W[0..m-1][0..n-1]; output matrix U
823 // @param n Number of columns
824 // @param m Number of rows
825 // @param w Ouput vector that presents diagonal matrix W
826 // @param V output matrix V
827 //
828 static void svd( double** a, int n, int m, double* w, double** v )
829 {
830  double * rv1;
831  int i, j, k, l = -1;
832  double tst1, c, f, g, h, s, scale;
833 
834  assert( m > 0 && n > 0 );
835 
836  rv1 = malloc( (size_t) n * sizeof ( double ) );
837 
838  //
839  // householder reduction to bidiagonal form
840  //
841  g = 0.0;
842  scale = 0.0;
843  tst1 = 0.0;
844  for ( i = 0; i < n; i++ )
845  {
846  l = i + 1;
847  rv1[i] = scale * g;
848  g = 0.0;
849  s = 0.0;
850  scale = 0.0;
851  if ( i < m )
852  {
853  for ( k = i; k < m; k++ )
854  scale += fabs( a[k][i] );
855  if ( scale != 0.0 )
856  {
857  for ( k = i; k < m; k++ )
858  {
859  a[k][i] /= scale;
860  s += a[k][i] * a[k][i];
861  }
862  f = a[i][i];
863  g = -copysign( sqrt( s ), f );
864  h = f * g - s;
865  a[i][i] = f - g;
866  if ( i < n - 1 )
867  {
868  for ( j = l; j < n; j++ )
869  {
870  s = 0.0;
871  for ( k = i; k < m; k++ )
872  s += a[k][i] * a[k][j];
873  f = s / h;
874  for ( k = i; k < m; k++ )
875  a[k][j] += f * a[k][i];
876  }
877  }
878  for ( k = i; k < m; k++ )
879  a[k][i] *= scale;
880  }
881  }
882  w[i] = scale * g;
883  g = 0.0;
884  s = 0.0;
885  scale = 0.0;
886  if ( i < m && i < n - 1 )
887  {
888  for ( k = l; k < n; k++ )
889  scale += fabs( a[i][k] );
890  if ( scale != 0.0 )
891  {
892  for ( k = l; k < n; k++ )
893  {
894  a[i][k] /= scale;
895  s += a[i][k] * a[i][k];
896  }
897  f = a[i][l];
898  g = -copysign( sqrt( s ), f );
899  h = f * g - s;
900  a[i][l] = f - g;
901  for ( k = l; k < n; k++ )
902  rv1[k] = a[i][k] / h;
903  for ( j = l; j < m; j++ )
904  {
905  s = 0.0;
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];
910  }
911  for ( k = l; k < n; k++ )
912  a[i][k] *= scale;
913  }
914  }
915  {
916  double tmp = fabs( w[i] ) + fabs( rv1[i] );
917 
918  tst1 = ( tst1 > tmp ) ? tst1 : tmp;
919  }
920  }
921 
922  //
923  // accumulation of right-hand transformations
924  //
925  for ( i = n - 1; i >= 0; i-- )
926  {
927  if ( i < n - 1 )
928  {
929  if ( g != 0.0 )
930  {
931  for ( j = l; j < n; j++ )
932  //
933  // double division avoids possible underflow
934  //
935  v[j][i] = ( a[i][j] / a[i][l] ) / g;
936  for ( j = l; j < n; j++ )
937  {
938  s = 0.0;
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];
943  }
944  }
945  for ( j = l; j < n; j++ )
946  {
947  v[i][j] = 0.0;
948  v[j][i] = 0.0;
949  }
950  }
951  v[i][i] = 1.0;
952  g = rv1[i];
953  l = i;
954  }
955 
956  //
957  // accumulation of left-hand transformations
958  //
959  for ( i = ( m < n ) ? m - 1 : n - 1; i >= 0; i-- )
960  {
961  l = i + 1;
962  g = w[i];
963  if ( i != n - 1 )
964  for ( j = l; j < n; j++ )
965  a[i][j] = 0.0;
966  if ( g != 0.0 )
967  {
968  for ( j = l; j < n; j++ )
969  {
970  s = 0.0;
971  for ( k = l; k < m; k++ )
972  s += a[k][i] * a[k][j];
973  //
974  // double division avoids possible underflow
975  //
976  f = ( s / a[i][i] ) / g;
977  for ( k = i; k < m; k++ )
978  a[k][j] += f * a[k][i];
979  }
980  for ( j = i; j < m; j++ )
981  a[j][i] /= g;
982  }
983  else
984  for ( j = i; j < m; j++ )
985  a[j][i] = 0.0;
986  a[i][i] += 1.0;
987  }
988 
989  //
990  // diagonalization of the bidiagonal form
991  //
992  for ( k = n - 1; k >= 0; k-- )
993  {
994  int k1 = k - 1;
995  int its = 0;
996 
997  while ( 1 )
998  {
999  int docancellation = 1;
1000  double x, y, z;
1001  int l1 = -1;
1002 
1003  its++;
1004  if ( its > SVD_NMAX )
1005  csa_quit( "svd(): no convergence in %d iterations", SVD_NMAX );
1006 
1007  for ( l = k; l >= 0; l-- ) // test for splitting
1008  {
1009  double tst2 = fabs( rv1[l] ) + tst1;
1010 
1011  if ( tst2 == tst1 )
1012  {
1013  docancellation = 0;
1014  break;
1015  }
1016  l1 = l - 1;
1017  //
1018  // rv1(1) is always zero, so there is no exit through the
1019  // bottom of the loop
1020  //
1021  tst2 = fabs( w[l - 1] ) + tst1;
1022  if ( tst2 == tst1 )
1023  break;
1024  }
1025  //
1026  // cancellation of rv1[l] if l > 1
1027  //
1028  if ( docancellation )
1029  {
1030  c = 0.0;
1031  s = 1.0;
1032  for ( i = l; i <= k; i++ )
1033  {
1034  f = s * rv1[i];
1035  rv1[i] = c * rv1[i];
1036  if ( ( fabs( f ) + tst1 ) == tst1 )
1037  break;
1038  g = w[i];
1039  h = hypot( f, g );
1040  w[i] = h;
1041  h = 1.0 / h;
1042  c = g * h;
1043  s = -f * h;
1044  for ( j = 0; j < m; j++ )
1045  {
1046  y = a[j][l1];
1047  z = a[j][i];
1048 
1049  a[j][l1] = y * c + z * s;
1050  a[j][i] = z * c - y * s;
1051  }
1052  }
1053  }
1054  //
1055  // test for convergence
1056  //
1057  z = w[k];
1058  if ( l != k )
1059  {
1060  int i1;
1061 
1062  //
1063  // shift from bottom 2 by 2 minor
1064  //
1065  x = w[l];
1066  y = w[k1];
1067  g = rv1[k1];
1068  h = rv1[k];
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 );
1072  //
1073  // next qr transformation
1074  //
1075  c = 1.0;
1076  s = 1.0;
1077  for ( i1 = l; i1 < k; i1++ )
1078  {
1079  i = i1 + 1;
1080  g = rv1[i];
1081  y = w[i];
1082  h = s * g;
1083  g = c * g;
1084  z = hypot( f, h );
1085  rv1[i1] = z;
1086  c = f / z;
1087  s = h / z;
1088  f = x * c + g * s;
1089  g = g * c - x * s;
1090  h = y * s;
1091  y *= c;
1092  for ( j = 0; j < n; j++ )
1093  {
1094  x = v[j][i1];
1095  z = v[j][i];
1096  v[j][i1] = x * c + z * s;
1097  v[j][i] = z * c - x * s;
1098  }
1099  z = hypot( f, h );
1100  w[i1] = z;
1101  //
1102  // rotation can be arbitrary if z = 0
1103  //
1104  if ( z != 0.0 )
1105  {
1106  c = f / z;
1107  s = h / z;
1108  }
1109  f = c * g + s * y;
1110  x = c * y - s * g;
1111  for ( j = 0; j < m; j++ )
1112  {
1113  y = a[j][i1];
1114  z = a[j][i];
1115  a[j][i1] = y * c + z * s;
1116  a[j][i] = z * c - y * s;
1117  }
1118  }
1119  rv1[l] = 0.0;
1120  rv1[k] = f;
1121  w[k] = x;
1122  }
1123  else
1124  {
1125  //
1126  // w[k] is made non-negative
1127  //
1128  if ( z < 0.0 )
1129  {
1130  w[k] = -z;
1131  for ( j = 0; j < n; j++ )
1132  v[j][k] = -v[j][k];
1133  }
1134  break;
1135  }
1136  }
1137  }
1138 
1139  free( rv1 );
1140 }
1141 
1142 // Least squares fitting via singular value decomposition.
1143 //
1144 static void lsq( double** A, int ni, int nj, double* z, double* w, double* sol )
1145 {
1146  double** V = alloc2d( ni, ni, sizeof ( double ) );
1147  double** B = alloc2d( nj, ni, sizeof ( double ) );
1148  int i, j, ii;
1149 
1150  svd( A, ni, nj, w, V );
1151 
1152  for ( j = 0; j < ni; ++j )
1153  for ( i = 0; i < ni; ++i )
1154  V[j][i] /= w[i];
1155  for ( i = 0; i < ni; ++i )
1156  {
1157  double* v = V[i];
1158 
1159  for ( j = 0; j < nj; ++j )
1160  {
1161  double* a = A[j];
1162  double* b = &B[i][j];
1163 
1164  for ( ii = 0; ii < ni; ++ii )
1165  *b += v[ii] * a[ii];
1166  }
1167  }
1168  for ( i = 0; i < ni; ++i )
1169  sol[i] = 0.0;
1170  for ( i = 0; i < ni; ++i )
1171  for ( j = 0; j < nj; ++j )
1172  sol[i] += B[i][j] * z[j];
1173 
1174  free2d( B );
1175  free2d( V );
1176 }
1177 
1178 //
1179 // square->coeffs[]:
1180 //
1181 // ---------------------
1182 // | 3 10 17 24 |
1183 // | 6 13 20 |
1184 // | 2 9 16 23 |
1185 // | 5 12 19 |
1186 // | 1 8 15 22 |
1187 // | 4 11 18 |
1188 // | 0 7 14 21 |
1189 // ---------------------
1190 //
1191 
1192 // Calculates spline coefficients in each primary triangle by least squares
1193 // fitting to data attached by csa_attachpoints().
1194 //
1195 static void csa_findprimarycoeffs( csa* a )
1196 {
1197  int n[4] = { 0, 0, 0, 0 };
1198  int i;
1199 
1200  if ( csa_verbose )
1201  fprintf( stderr, "calculating spline coefficients for primary triangles:\n " );
1202 
1203  for ( i = 0; i < a->npt; ++i )
1204  {
1205  triangle* t = a->pt[i];
1206  int npoints = t->npoints;
1207  point ** points = t->points;
1208  double * z = malloc( (size_t) npoints * sizeof ( double ) );
1209  int q = n2q( t->npoints );
1210  int ok = 1;
1211  double b[10];
1212  double b1[6];
1213  int ii;
1214 
1215  if ( csa_verbose )
1216  {
1217  fprintf( stderr, "." );
1218  fflush( stderr );
1219  }
1220 
1221  for ( ii = 0; ii < npoints; ++ii )
1222  z[ii] = points[ii]->z;
1223 
1224  do
1225  {
1226  double bc[3];
1227  double wmin, wmax;
1228 
1229  if ( !ok )
1230  q--;
1231 
1232  assert( q >= 0 );
1233 
1234  if ( q == 3 )
1235  {
1236  double ** A = alloc2d( 10, npoints, sizeof ( double ) );
1237  double w[10];
1238 
1239  for ( ii = 0; ii < npoints; ++ii )
1240  {
1241  double * aii = A[ii];
1242  double tmp;
1243 
1244  triangle_calculatebc( t, points[ii], bc );
1245 
1246  //
1247  // 0 1 2 3 4 5 6 7 8 9
1248  // 300 210 201 120 111 102 030 021 012 003
1249  //
1250  tmp = bc[0] * bc[0];
1251  aii[0] = tmp * bc[0];
1252  tmp *= 3.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];
1257  tmp *= 3.0;
1258  aii[3] = tmp * bc[0];
1259  aii[7] = tmp * bc[2];
1260  tmp = bc[2] * bc[2];
1261  aii[9] = tmp * bc[2];
1262  tmp *= 3.0;
1263  aii[5] = tmp * bc[0];
1264  aii[8] = tmp * bc[1];
1265  aii[4] = bc[0] * bc[1] * bc[2] * 6.0;
1266  }
1267 
1268  lsq( A, 10, npoints, z, w, b );
1269 
1270  wmin = w[0];
1271  wmax = w[0];
1272  for ( ii = 1; ii < 10; ++ii )
1273  {
1274  if ( w[ii] < wmin )
1275  wmin = w[ii];
1276  else if ( w[ii] > wmax )
1277  wmax = w[ii];
1278  }
1279  if ( wmin < wmax / a->k )
1280  ok = 0;
1281 
1282  free2d( A );
1283  }
1284  else if ( q == 2 )
1285  {
1286  double ** A = alloc2d( 6, npoints, sizeof ( double ) );
1287  double w[6];
1288 
1289  for ( ii = 0; ii < npoints; ++ii )
1290  {
1291  double* aii = A[ii];
1292 
1293  triangle_calculatebc( t, points[ii], bc );
1294 
1295  //
1296  // 0 1 2 3 4 5
1297  // 200 110 101 020 011 002
1298  //
1299 
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];
1306  }
1307 
1308  lsq( A, 6, npoints, z, w, b1 );
1309 
1310  wmin = w[0];
1311  wmax = w[0];
1312  for ( ii = 1; ii < 6; ++ii )
1313  {
1314  if ( w[ii] < wmin )
1315  wmin = w[ii];
1316  else if ( w[ii] > wmax )
1317  wmax = w[ii];
1318  }
1319  if ( wmin < wmax / a->k )
1320  ok = 0;
1321  else // degree raising
1322  {
1323  ok = 1;
1324  b[0] = b1[0];
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;
1330  b[6] = b1[3];
1331  b[7] = ( b1[3] + 2.0 * b1[4] ) / 3.0;
1332  b[8] = ( b1[5] + 2.0 * b1[4] ) / 3.0;
1333  b[9] = b1[5];
1334  }
1335 
1336  free2d( A );
1337  }
1338  else if ( q == 1 )
1339  {
1340  double ** A = alloc2d( 3, npoints, sizeof ( double ) );
1341  double w[3];
1342 
1343  for ( ii = 0; ii < npoints; ++ii )
1344  {
1345  double* aii = A[ii];
1346 
1347  triangle_calculatebc( t, points[ii], bc );
1348 
1349  aii[0] = bc[0];
1350  aii[1] = bc[1];
1351  aii[2] = bc[2];
1352  }
1353 
1354  lsq( A, 3, npoints, z, w, b1 );
1355 
1356  wmin = w[0];
1357  wmax = w[0];
1358  for ( ii = 1; ii < 3; ++ii )
1359  {
1360  if ( w[ii] < wmin )
1361  wmin = w[ii];
1362  else if ( w[ii] > wmax )
1363  wmax = w[ii];
1364  }
1365  if ( wmin < wmax / a->k )
1366  ok = 0;
1367  else // degree raising
1368  {
1369  ok = 1;
1370  b[0] = b1[0];
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;
1376  b[6] = b1[1];
1377  b[7] = ( 2.0 * b1[1] + b1[2] ) / 3.0;
1378  b[8] = ( 2.0 * b1[2] + b1[1] ) / 3.0;
1379  b[9] = b1[2];
1380  }
1381 
1382  free2d( A );
1383  }
1384  else if ( q == 0 )
1385  {
1386  double ** A = alloc2d( 1, npoints, sizeof ( double ) );
1387  double w[1];
1388 
1389  for ( ii = 0; ii < npoints; ++ii )
1390  A[ii][0] = 1.0;
1391 
1392  lsq( A, 1, npoints, z, w, b1 );
1393 
1394  ok = 1;
1395  b[0] = b1[0];
1396  b[1] = b1[0];
1397  b[2] = b1[0];
1398  b[3] = b1[0];
1399  b[4] = b1[0];
1400  b[5] = b1[0];
1401  b[6] = b1[0];
1402  b[7] = b1[0];
1403  b[8] = b1[0];
1404  b[9] = b1[0];
1405 
1406  free2d( A );
1407  }
1408  } while ( !ok );
1409 
1410  n[q]++;
1411  t->order = q;
1412 
1413  {
1414  square* s = t->parent;
1415  double* coeffs = s->coeffs;
1416 
1417  coeffs[12] = b[0];
1418  coeffs[9] = b[1];
1419  coeffs[6] = b[3];
1420  coeffs[3] = b[6];
1421  coeffs[2] = b[7];
1422  coeffs[1] = b[8];
1423  coeffs[0] = b[9];
1424  coeffs[4] = b[5];
1425  coeffs[8] = b[2];
1426  coeffs[5] = b[4];
1427  }
1428 
1429  free( z );
1430  }
1431 
1432  if ( csa_verbose )
1433  {
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] );
1438  fflush( stderr );
1439  }
1440 
1441  if ( csa_verbose == 2 )
1442  {
1443  int j;
1444 
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 )
1450  {
1451  fprintf( stderr, "%2d ", j );
1452  for ( i = 0; i < a->ni; ++i )
1453  {
1454  square* s = a->squares[j][i];
1455 
1456  if ( s->triangles[0]->primary )
1457  fprintf( stderr, "%2d ", s->triangles[0]->order );
1458  else
1459  fprintf( stderr, " . " );
1460  }
1461  fprintf( stderr, "\n" );
1462  }
1463  }
1464 }
1465 
1466 // Finds spline coefficients in (adjacent to primary triangles) secondary
1467 // triangles from C1 smoothness conditions.
1468 //
1469 static void csa_findsecondarycoeffs( csa* a )
1470 {
1471  square*** squares = a->squares;
1472  int ni = a->ni;
1473  int nj = a->nj;
1474  int ii;
1475 
1476  if ( csa_verbose )
1477  {
1478  fprintf( stderr, "propagating spline coefficients to the remaining triangles:\n" );
1479  fflush( stderr );
1480  }
1481 
1482  //
1483  // red
1484  //
1485  for ( ii = 0; ii < a->npt; ++ii )
1486  {
1487  triangle* t = a->pt[ii];
1488  square * s = t->parent;
1489  int i = s->i;
1490  int j = s->j;
1491  double * c = s->coeffs;
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;
1497 
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];
1501 
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];
1505 
1506  c[19] = 2.0 * c[15] - c[11];
1507 
1508  if ( cl != NULL )
1509  {
1510  cl[21] = c[0];
1511  cl[22] = c[1];
1512  cl[23] = c[2];
1513  cl[24] = c[3];
1514 
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];
1518 
1519  cl[17] = 2.0 * cl[20] - cl[23];
1520  cl[14] = 2.0 * cl[18] - cl[22];
1521  }
1522 
1523  if ( cb != NULL )
1524  {
1525  cb[3] = c[0];
1526  cb[10] = c[7];
1527 
1528  cb[6] = c[0] + c[7] - c[4];
1529  cb[2] = 2.0 * cb[6] - cb[10];
1530  }
1531 
1532  if ( cbl != NULL )
1533  {
1534  cbl[23] = cb[2];
1535  cbl[24] = cb[3];
1536 
1537  cbl[20] = cb[2] + cb[3] - cb[6];
1538  cbl[17] = cl[14];
1539  }
1540 
1541  if ( ca != NULL )
1542  {
1543  ca[0] = c[3];
1544  ca[7] = c[10];
1545 
1546  ca[4] = c[3] + c[10] - c[6];
1547  ca[1] = 2.0 * ca[4] - ca[7];
1548  }
1549 
1550  if ( cal != NULL )
1551  {
1552  cal[21] = c[3];
1553  cal[22] = ca[1];
1554 
1555  cal[18] = ca[0] + ca[1] - ca[4];
1556  cal[14] = cl[17];
1557  }
1558  }
1559 
1560  //
1561  // blue
1562  //
1563  for ( ii = 0; ii < a->npt; ++ii )
1564  {
1565  triangle* t = a->pt[ii];
1566  square * s = t->parent;
1567  int i = s->i;
1568  int j = s->j;
1569  double * c = s->coeffs;
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;
1573 
1574  if ( car != NULL )
1575  cr[13] = car[7] + car[14] - car[11];
1576 
1577  if ( cbr != NULL )
1578  cr[11] = cbr[10] + cbr[17] - cbr[13];
1579 
1580  if ( cr != NULL )
1581  cr[5] = c[22] + c[23] - c[19];
1582  }
1583 
1584  //
1585  // green & yellow
1586  //
1587  for ( ii = 0; ii < a->npt; ++ii )
1588  {
1589  triangle* t = a->pt[ii];
1590  square * s = t->parent;
1591  int i = s->i;
1592  int j = s->j;
1593  double * cr = ( i < ni - 1 ) ? squares[j][i + 1]->coeffs : NULL;
1594 
1595  if ( cr != NULL )
1596  {
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;
1602  }
1603  }
1604 
1605  if ( csa_verbose )
1606  {
1607  fprintf( stderr, "checking that all coefficients have been set:\n" );
1608  fflush( stderr );
1609  }
1610 
1611  for ( ii = 0; ii < ni * nj; ++ii )
1612  {
1613  square* s = squares[0][ii];
1614  double* c = s->coeffs;
1615  int i;
1616 
1617  if ( s->npoints == 0 )
1618  continue;
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 );
1622  }
1623 }
1624 
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 };
1635 
1636 static int * iall[] = { i300, i030, i003, i210, i021, i102, i120, i012, i201, i111 };
1637 
1638 static void csa_sethascoeffsflag( csa* a )
1639 {
1640  int i, j;
1641 
1642  for ( j = 0; j < a->nj; ++j )
1643  {
1644  for ( i = 0; i < a->ni; ++i )
1645  {
1646  square* s = a->squares[j][i];
1647  double* coeffs = s->coeffs;
1648  int ii;
1649 
1650  for ( ii = 0; ii < 4; ++ii )
1651  {
1652  triangle* t = s->triangles[ii];
1653  int cc;
1654 
1655  for ( cc = 0; cc < 10; ++cc )
1656  if ( isnan( coeffs[iall[cc][ii]] ) )
1657  break;
1658  if ( cc == 10 )
1659  t->hascoeffs = 1;
1660  }
1661  }
1662  }
1663 }
1664 
1666 {
1667  csa_squarize( a );
1668  csa_attachpoints( a );
1669  csa_findprimarycoeffs( a );
1671  csa_sethascoeffsflag( a );
1672 }
1673 
1675 {
1676  double h = a->h;
1677  double ii = ( p->x - a->xmin ) / h + 1.0;
1678  double jj = ( p->y - a->ymin ) / h + 1.0;
1679  int i, j;
1680  square * s;
1681  double fi, fj;
1682  int ti;
1683  triangle* t;
1684  double bc[3];
1685 
1686  if ( ii < 0.0 || jj < 0.0 || ii > (double) a->ni - 1.0 || jj > (double) a->nj - 1.0 )
1687  {
1688  p->z = NaN;
1689  return;
1690  }
1691 
1692  i = (int) floor( ii );
1693  j = (int) floor( jj );
1694  s = a->squares[j][i];
1695  fi = ii - i;
1696  fj = jj - j;
1697 
1698  if ( fj < fi )
1699  {
1700  if ( fi + fj < 1.0 )
1701  ti = 3;
1702  else
1703  ti = 2;
1704  }
1705  else
1706  {
1707  if ( fi + fj < 1.0 )
1708  ti = 0;
1709  else
1710  ti = 1;
1711  }
1712 
1713  t = s->triangles[ti];
1714  if ( !t->hascoeffs )
1715  {
1716  p->z = NaN;
1717  return;
1718  }
1719  triangle_calculatebc( t, p, bc );
1720 
1721  {
1722  double * c = s->coeffs;
1723  double bc1 = bc[0];
1724  double bc2 = bc[1];
1725  double bc3 = bc[2];
1726  double tmp1 = bc1 * bc1;
1727  double tmp2 = bc2 * bc2;
1728  double tmp3 = bc3 * bc3;
1729 
1730  switch ( ti )
1731  {
1732  case 0:
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;
1734  break;
1735  case 1:
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;
1737  break;
1738  case 2:
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;
1740  break;
1741  default: // 3
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;
1743  }
1744  }
1745 }
1746 
1747 void csa_approximate_points( csa* a, int n, point* points )
1748 {
1749  int ii;
1750 
1751  for ( ii = 0; ii < n; ++ii )
1752  csa_approximate_point( a, &points[ii] );
1753 }
1754 
1755 void csa_setnpmin( csa* a, int npmin )
1756 {
1757  a->npmin = npmin;
1758 }
1759 
1760 void csa_setnpmax( csa* a, int npmax )
1761 {
1762  a->npmax = npmax;
1763 }
1764 
1765 void csa_setk( csa* a, int k )
1766 {
1767  a->k = k;
1768 }
1769 
1770 void csa_setnppc( csa* a, double nppc )
1771 {
1772  a->nppc = (int) nppc;
1773 }
1774 
1775 #if defined ( STANDALONE )
1776 
1777 #include "minell.h"
1778 
1779 #define NIMAX 2048
1780 #define BUFSIZE 10240
1781 #define STRBUFSIZE 64
1782 
1783 static void points_generate( double xmin, double xmax, double ymin, double ymax, int nx, int ny, int* nout, point** pout )
1784 {
1785  double stepx, stepy;
1786  double x0, xx, yy;
1787  int i, j, ii;
1788 
1789  if ( nx < 1 || ny < 1 )
1790  {
1791  *pout = NULL;
1792  *nout = 0;
1793  return;
1794  }
1795 
1796  *nout = nx * ny;
1797  *pout = malloc( *nout * sizeof ( point ) );
1798 
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;
1803 
1804  ii = 0;
1805  for ( j = 0; j < ny; ++j )
1806  {
1807  xx = x0;
1808  for ( i = 0; i < nx; ++i )
1809  {
1810  point* p = &( *pout )[ii];
1811 
1812  p->x = xx;
1813  p->y = yy;
1814  xx += stepx;
1815  ii++;
1816  }
1817  yy += stepy;
1818  }
1819 }
1820 
1821 static int str2double( char* token, double* value )
1822 {
1823  char* end = NULL;
1824 
1825  if ( token == NULL )
1826  {
1827  *value = NaN;
1828  return 0;
1829  }
1830 
1831  *value = strtod( token, &end );
1832 
1833  if ( end == token )
1834  {
1835  *value = NaN;
1836  return 0;
1837  }
1838 
1839  return 1;
1840 }
1841 
1842 #define NALLOCATED_START 1024
1843 
1844 // Reads array of points from a columnar file.
1845 //
1846 // @param fname File name (can be "stdin" or "-" for standard input)
1847 // @param dim Number of dimensions (must be 2 or 3)
1848 // @param n Pointer to number of points (output)
1849 // @param points Pointer to array of points [*n] (output) (to be freed)
1850 //
1851 void points_read( char* fname, int dim, int* n, point** points )
1852 {
1853  FILE * f = NULL;
1854  int nallocated = NALLOCATED_START;
1855  char buf[BUFSIZE];
1856  char seps[] = " ,;\t";
1857  char * token;
1858 
1859  if ( dim < 2 || dim > 3 )
1860  {
1861  *n = 0;
1862  *points = NULL;
1863  return;
1864  }
1865 
1866  if ( fname == NULL )
1867  f = stdin;
1868  else
1869  {
1870  if ( strcmp( fname, "stdin" ) == 0 || strcmp( fname, "-" ) == 0 )
1871  f = stdin;
1872  else
1873  {
1874  f = fopen( fname, "r" );
1875  if ( f == NULL )
1876  csa_quit( "%s: %s\n", fname, strerror( errno ) );
1877  }
1878  }
1879 
1880  *points = malloc( nallocated * sizeof ( point ) );
1881  *n = 0;
1882  while ( fgets( buf, BUFSIZE, f ) != NULL )
1883  {
1884  point* p;
1885 
1886  if ( *n == nallocated )
1887  {
1888  nallocated *= 2;
1889  *points = realloc( *points, nallocated * sizeof ( point ) );
1890  }
1891 
1892  p = &( *points )[*n];
1893 
1894  if ( buf[0] == '#' )
1895  continue;
1896  if ( ( token = strtok( buf, seps ) ) == NULL )
1897  continue;
1898  if ( !str2double( token, &p->x ) )
1899  continue;
1900  if ( ( token = strtok( NULL, seps ) ) == NULL )
1901  continue;
1902  if ( !str2double( token, &p->y ) )
1903  continue;
1904  if ( dim == 2 )
1905  p->z = NaN;
1906  else
1907  {
1908  if ( ( token = strtok( NULL, seps ) ) == NULL )
1909  continue;
1910  if ( !str2double( token, &p->z ) )
1911  continue;
1912  }
1913  ( *n )++;
1914  }
1915 
1916  if ( *n == 0 )
1917  {
1918  free( *points );
1919  *points = NULL;
1920  }
1921  else
1922  *points = realloc( *points, *n * sizeof ( point ) );
1923 
1924  if ( f != stdin )
1925  if ( fclose( f ) != 0 )
1926  csa_quit( "%s: %s\n", fname, strerror( errno ) );
1927 }
1928 
1929 static void points_write( int n, point* points )
1930 {
1931  int i;
1932 
1933  if ( csa_verbose )
1934  printf( "Output:\n" );
1935 
1936  for ( i = 0; i < n; ++i )
1937  {
1938  point* p = &points[i];
1939 
1940  printf( "%.15g %.15g %.15g\n", p->x, p->y, p->z );
1941  }
1942 }
1943 
1944 static double points_scaletosquare( int n, point* points )
1945 {
1946  double xmin, ymin, xmax, ymax;
1947  double k;
1948  int i;
1949 
1950  if ( n <= 0 )
1951  return NaN;
1952 
1953  xmin = xmax = points[0].x;
1954  ymin = ymax = points[0].y;
1955 
1956  for ( i = 1; i < n; ++i )
1957  {
1958  point* p = &points[i];
1959 
1960  if ( p->x < xmin )
1961  xmin = p->x;
1962  else if ( p->x > xmax )
1963  xmax = p->x;
1964  if ( p->y < ymin )
1965  ymin = p->y;
1966  else if ( p->y > ymax )
1967  ymax = p->y;
1968  }
1969 
1970  if ( xmin == xmax || ymin == ymax )
1971  return NaN;
1972  else
1973  k = ( ymax - ymin ) / ( xmax - xmin );
1974 
1975  for ( i = 0; i < n; ++i )
1976  points[i].y /= k;
1977 
1978  return k;
1979 }
1980 
1981 static void points_scale( int n, point* points, double k )
1982 {
1983  int i;
1984 
1985  for ( i = 0; i < n; ++i )
1986  points[i].y /= k;
1987 }
1988 
1989 static void usage()
1990 {
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" );
2015 
2016  exit( 0 );
2017 }
2018 
2019 static void version()
2020 {
2021  printf( "csa version %s\n", csa_version );
2022  exit( 0 );
2023 }
2024 
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 )
2026 {
2027  int i;
2028 
2029  if ( argc < 2 )
2030  usage();
2031 
2032  i = 1;
2033  while ( i < argc )
2034  {
2035  if ( argv[i][0] != '-' )
2036  usage();
2037 
2038  switch ( argv[i][1] )
2039  {
2040  case 'c':
2041  i++;
2042  *invariant = 1;
2043  *square = 0;
2044 
2045  break;
2046  case 'i':
2047  i++;
2048  if ( i >= argc )
2049  csa_quit( "no file name found after -i\n" );
2050  *fdata = argv[i];
2051  i++;
2052  break;
2053  case 'n':
2054  i++;
2055  *fpoints = NULL;
2056  *generate_points = 1;
2057  if ( i >= argc )
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" );
2063  i++;
2064  break;
2065  case 'o':
2066  i++;
2067  if ( i >= argc )
2068  csa_quit( "no file name found after -o\n" );
2069  *fpoints = argv[i];
2070  i++;
2071  break;
2072  case 's':
2073  i++;
2074  *square = 1;
2075 
2076  *invariant = 0;
2077  break;
2078  case 'v':
2079  i++;
2080  csa_verbose = 1;
2081  break;
2082  case 'z':
2083  i++;
2084  if ( i >= argc )
2085  csa_quit( "no zoom value found after -z\n" );
2086  *zoom = atof( argv[i] );
2087  i++;
2088  break;
2089  case 'P': {
2090  char delim[] = "=";
2091  char prmstr[STRBUFSIZE] = "";
2092  char * token;
2093 
2094  i++;
2095  if ( i >= argc )
2096  csa_quit( "no input found after -P\n" );
2097 
2098  if ( strlen( argv[i] ) >= STRBUFSIZE )
2099  csa_quit( "could not interpret \"%s\" after -P option\n", argv[i] );
2100 
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] );
2105 
2106  if ( strcmp( token, "nppc" ) == 0 )
2107  {
2108  long int n;
2109 
2110  token = strtok( NULL, delim );
2111  if ( token == NULL )
2112  csa_quit( "could not interpret \"%s\" after -P option\n", argv[i] );
2113 
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] );
2117  else if ( n <= 0 )
2118  csa_quit( "non-sensible value for \"nppc\" parameter\n" );
2119  *nppc = (int) n;
2120  }
2121  else if ( strcmp( token, "k" ) == 0 )
2122  {
2123  long int n;
2124 
2125  token = strtok( NULL, delim );
2126  if ( token == NULL )
2127  csa_quit( "could not interpret \"%s\" after -P option\n", argv[i] );
2128 
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] );
2132  else if ( n <= 0 )
2133  csa_quit( "non-sensible value for \"k\" parameter\n" );
2134  *k = (int) n;
2135  }
2136  else
2137  usage();
2138 
2139  i++;
2140  break;
2141  }
2142  case 'V':
2143  i++;
2144  csa_verbose = 2;
2145  break;
2146  default:
2147  usage();
2148  break;
2149  }
2150  }
2151 
2152  if ( csa_verbose && argc == 2 )
2153  version();
2154 }
2155 
2156 int main( int argc, char* argv[] )
2157 {
2158  char * fdata = NULL;
2159  char * fpoints = NULL;
2160  int nin = 0;
2161  point * pin = NULL;
2162  int invariant = 0;
2163  minell * me = NULL;
2164  int square = 0;
2165  int nout = 0;
2166  int generate_points = 0;
2167  point * pout = NULL;
2168  int nx = -1;
2169  int ny = -1;
2170  csa * a = NULL;
2171  int nppc = -1;
2172  int k = -1;
2173  double ks = NaN;
2174  double zoom = NaN;
2175 
2176  parse_commandline( argc, argv, &fdata, &fpoints, &invariant, &square, &generate_points, &nx, &ny, &nppc, &k, &zoom );
2177 
2178  if ( fdata == NULL )
2179  csa_quit( "no input data\n" );
2180 
2181  if ( !generate_points && fpoints == NULL )
2182  csa_quit( "no output grid specified\n" );
2183 
2184  points_read( fdata, 3, &nin, &pin );
2185 
2186  if ( nin < 3 )
2187  return 0;
2188 
2189  if ( invariant )
2190  {
2191  me = minell_build( nin, pin );
2192  minell_scalepoints( me, nin, pin );
2193  }
2194  else if ( square )
2195  ks = points_scaletosquare( nin, pin );
2196 
2197  a = csa_create();
2198  csa_addpoints( a, nin, pin );
2199  if ( nppc > 0 )
2200  csa_setnppc( a, nppc );
2201  if ( k > 0 )
2202  csa_setk( a, k );
2203  csa_calculatespline( a );
2204 
2205  if ( generate_points )
2206  {
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 );
2209  else
2210  {
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;
2215 
2216  points_generate( xav - xdiff2 * zoom, xav + xdiff2 * zoom, yav - ydiff2 * zoom, yav + ydiff2 * zoom, nx, ny, &nout, &pout );
2217  }
2218  }
2219  else
2220  {
2221  points_read( fpoints, 2, &nout, &pout );
2222  if ( invariant )
2223  minell_scalepoints( me, nout, pout );
2224  else if ( square )
2225  points_scale( nout, pout, ks );
2226  }
2227 
2228  csa_approximate_points( a, nout, pout );
2229 
2230  if ( invariant )
2231  minell_rescalepoints( me, nout, pout );
2232  else if ( square )
2233  points_scale( nout, pout, 1.0 / ks );
2234 
2235  points_write( nout, pout );
2236 
2237  csa_destroy( a );
2238  free( pin );
2239  free( pout );
2240 
2241  return 0;
2242 }
2243 
2244 #endif // STANDALONE