Menu
▾
▴
[Plplot-cvs] SF.net SVN: plplot:[10540] trunk
|
From: <ai...@us...> - 2009-10-21 02:21:31
|
Revision: 10540
http://plplot.svn.sourceforge.net/plplot/?rev=10540&view=rev
Author: airwin
Date: 2009-10-21 02:21:21 +0000 (Wed, 21 Oct 2009)
Log Message:
-----------
Style all C source code in the lib subdirectories.
Modified Paths:
--------------
trunk/lib/csa/csa.c
trunk/lib/csa/csa.h
trunk/lib/csa/csadll.h
trunk/lib/csa/nan.h
trunk/lib/csa/version.h
trunk/lib/nistcd/cd.c
trunk/lib/nistcd/cd.h
trunk/lib/nistcd/cddll.h
trunk/lib/nistcd/cdexpert.c
trunk/lib/nistcd/cdmulti.c
trunk/lib/nistcd/cdsimple.c
trunk/lib/nistcd/cdtest.c
trunk/lib/nistcd/cdtext.c
trunk/lib/nistcd/color16.c
trunk/lib/nistcd/defines.h
trunk/lib/nn/delaunay.c
trunk/lib/nn/delaunay.h
trunk/lib/nn/hash.c
trunk/lib/nn/hash.h
trunk/lib/nn/istack.c
trunk/lib/nn/istack.h
trunk/lib/nn/lpi.c
trunk/lib/nn/nan.h
trunk/lib/nn/nn.h
trunk/lib/nn/nnai.c
trunk/lib/nn/nncommon.c
trunk/lib/nn/nndll.h
trunk/lib/nn/nnpi.c
trunk/lib/nn/version.h
trunk/lib/qsastime/bhunt_search_test.c
trunk/lib/qsastime/deltaT-gen.c
trunk/lib/qsastime/deltaT_test.c
trunk/lib/qsastime/dspline.c
trunk/lib/qsastime/dsplint.c
trunk/lib/qsastime/qsastime.c
trunk/lib/qsastime/qsastime.h
trunk/lib/qsastime/qsastimeP.h.in
trunk/lib/qsastime/qsastime_extra.c
trunk/lib/qsastime/qsastime_extra.h
trunk/lib/qsastime/qsastime_test.c
trunk/lib/qsastime/qsastime_testlib.c
trunk/lib/qsastime/qsastimedll.h
trunk/lib/qsastime/tai-utc-gen.c
trunk/scripts/style_source.sh
Modified: trunk/lib/csa/csa.c
===================================================================
--- trunk/lib/csa/csa.c 2009-10-21 01:35:56 UTC (rev 10539)
+++ trunk/lib/csa/csa.c 2009-10-21 02:21:21 UTC (rev 10540)
@@ -14,7 +14,7 @@
* Smooth approximation and rendering of large scattered data
* sets, in ``Proceedings of IEEE Visualization 2001''
* (Th.Ertl, K.Joy and A.Varshney, Eds.), pp.341-347, 571,
- * IEEE Computer Society, 2001.
+ * IEEE Computer Society, 2001.
* http://www.uni-giessen.de/www-Numerische-Mathematik/
* davydov/VIS2001.ps.gz
* http://www.math.uni-mannheim.de/~lsmath4/paper/
@@ -40,86 +40,89 @@
int csa_verbose = 0;
-#define NPASTART 5 /* Number of Points Allocated at Start */
-#define SVD_NMAX 30 /* Maximal number of iterations in svd() */
+#define NPASTART 5 /* Number of Points Allocated at Start */
+#define SVD_NMAX 30 /* Maximal number of iterations in svd() */
/* default algorithm parameters */
-#define NPMIN_DEF 3
-#define NPMAX_DEF 40
-#define K_DEF 140
-#define NPPC_DEF 5
+#define NPMIN_DEF 3
+#define NPMAX_DEF 40
+#define K_DEF 140
+#define NPPC_DEF 5
struct square;
-typedef struct square square;
+typedef struct square square;
-typedef struct {
- square* parent;
- int index; /* index within parent square; 0 <= index <=
+typedef struct
+{
+ square * parent;
+ int index; /* index within parent square; 0 <= index <=
* 3 */
- point vertices[3];
- point middle; /* barycenter */
+ point vertices[3];
+ point middle; /* barycenter */
double h; /* parent square edge length */
double r; /* data visibility radius */
/*
- * points used -- in primary triangles only
+ * points used -- in primary triangles only
*/
- int nallocated;
- int npoints;
+ int nallocated;
+ int npoints;
point** points;
- int primary; /* flag -- whether calculate spline
+ int primary; /* flag -- whether calculate spline
* coefficients directly (by least squares
* method) (primary = 1) or indirectly (from
* C1 smoothness conditions) (primary = 0) */
int hascoeffs; /* flag -- whether there are no NaNs among
* the spline coefficients */
- int order; /* spline order -- for primary triangles only
+ int order; /* spline order -- for primary triangles only
*/
} triangle;
-struct square {
- csa* parent;
- int i, j; /* indices */
+struct square
+{
+ csa * parent;
+ int i, j; /* indices */
- int nallocated;
- int npoints;
- point** points;
+ int nallocated;
+ int npoints;
+ point ** points;
- int primary; /* flag -- whether this square contains a
+ int primary; /* flag -- whether this square contains a
* primary triangle */
triangle* triangles[4];
- double coeffs[25];
+ double coeffs[25];
};
-struct csa {
+struct csa
+{
double xmin;
double xmax;
double ymin;
double ymax;
- int nallocated;
- int npoints;
- point** points;
+ int nallocated;
+ int npoints;
+ point ** points;
/*
- * squarization
+ * squarization
*/
- int ni;
- int nj;
- double h;
- square*** squares; /* square* [j][i] */
+ int ni;
+ int nj;
+ double h;
+ square *** squares; /* square* [j][i] */
- int npt; /* Number of Primary Triangles */
+ int npt; /* Number of Primary Triangles */
triangle** pt; /* Primary Triangles -- triangle* [npt] */
/*
- * algorithm parameters
+ * algorithm parameters
*/
- int npmin; /* minimal number of points locally involved
+ int npmin; /* minimal number of points locally involved
* in * spline calculation (normally = 3) */
- int npmax; /* maximal number of points locally involved
+ int npmax; /* maximal number of points locally involved
* in * spline calculation (required > 10, *
* recommended 20 < npmax < 60) */
double k; /* relative tolerance multiple in fitting
@@ -129,48 +132,48 @@
int nppc; /* average number of points per square */
};
-static void csa_quit(char* format, ...)
+static void csa_quit( char* format, ... )
{
va_list args;
- fflush(stdout); /* just in case -- to have the exit message
+ fflush( stdout ); /* just in case -- to have the exit message
* last */
- fprintf(stderr, "error: csa: ");
- va_start(args, format);
- vfprintf(stderr, format, args);
- va_end(args);
+ fprintf( stderr, "error: csa: " );
+ va_start( args, format );
+ vfprintf( stderr, format, args );
+ va_end( args );
- exit(1);
+ exit( 1 );
}
-/* Allocates n1xn2 matrix of something. Note that it will be accessed as
+/* Allocates n1xn2 matrix of something. Note that it will be accessed as
* [n2][n1].
* @param n1 Number of columns
* @param n2 Number of rows
* @return Matrix
*/
-static void* alloc2d(int n1, int n2, size_t unitsize)
+static void* alloc2d( int n1, int n2, size_t unitsize )
{
unsigned int size;
- char* p;
- char** pp;
- int i;
+ char * p;
+ char ** pp;
+ int i;
- assert(n1 > 0);
- assert(n2 > 0);
- assert((double) n1 * (double) n2 <= (double) UINT_MAX);
+ assert( n1 > 0 );
+ assert( n2 > 0 );
+ assert((double) n1 * (double) n2 <= (double) UINT_MAX );
size = n1 * n2;
- if ((p = calloc(size, unitsize)) == NULL)
- csa_quit("alloc2d(): %s\n", strerror(errno));
+ if (( p = calloc( size, unitsize )) == NULL )
+ csa_quit( "alloc2d(): %s\n", strerror( errno ));
- assert((double) n2 * (double) sizeof(void*) <= (double) UINT_MAX);
+ assert((double) n2 * (double) sizeof ( void* ) <= (double) UINT_MAX );
- size = n2 * sizeof(void*);
- if ((pp = malloc(size)) == NULL)
- csa_quit("alloc2d(): %s\n", strerror(errno));
- for (i = 0; i < n2; i++)
+ size = n2 * sizeof ( void* );
+ if (( pp = malloc( size )) == NULL )
+ csa_quit( "alloc2d(): %s\n", strerror( errno ));
+ for ( i = 0; i < n2; i++ )
pp[i] = &p[i * n1 * unitsize];
return pp;
@@ -179,47 +182,51 @@
/* Destroys n1xn2 matrix.
* @param pp Matrix
*/
-static void free2d(void* pp)
+static void free2d( void* pp )
{
void* p;
- assert(pp != NULL);
- p = ((void**) pp)[0];
- free(pp);
- assert(p != NULL);
- free(p);
+ assert( pp != NULL );
+ p = ((void**) pp )[0];
+ free( pp );
+ assert( p != NULL );
+ free( p );
}
-static triangle* triangle_create(square* s, point vertices[], int index)
+static triangle* triangle_create( square* s, point vertices[], int index )
{
- triangle* t = malloc(sizeof(triangle));
+ triangle* t = malloc( sizeof ( triangle ));
t->parent = s;
- memcpy(t->vertices, vertices, sizeof(point) * 3);
- t->middle.x = (vertices[0].x + vertices[1].x + vertices[2].x) / 3.0;
- t->middle.y = (vertices[0].y + vertices[1].y + vertices[2].y) / 3.0;
- t->h = s->parent->h;
- t->index = index;
+ memcpy( t->vertices, vertices, sizeof ( point ) * 3 );
+ t->middle.x = ( vertices[0].x + vertices[1].x + vertices[2].x ) / 3.0;
+ t->middle.y = ( vertices[0].y + vertices[1].y + vertices[2].y ) / 3.0;
+ t->h = s->parent->h;
+ t->index = index;
- t->r = 0.0;
- t->points = 0;
+ t->r = 0.0;
+ t->points = 0;
t->nallocated = 0;
- t->npoints = 0;
- t->hascoeffs = 0;
- t->primary = 0;
- t->order = -1;
+ t->npoints = 0;
+ t->hascoeffs = 0;
+ t->primary = 0;
+ t->order = -1;
return t;
}
-static void triangle_addpoint(triangle* t, point* p)
+static void triangle_addpoint( triangle* t, point* p )
{
- if (t->nallocated == t->npoints) {
- if (t->nallocated == 0) {
- t->points = malloc(NPASTART * sizeof(point*));
+ if ( t->nallocated == t->npoints )
+ {
+ if ( t->nallocated == 0 )
+ {
+ t->points = malloc( NPASTART * sizeof ( point* ));
t->nallocated = NPASTART;
- } else {
- t->points = realloc(t->points, t->nallocated * 2 * sizeof(point*));
+ }
+ else
+ {
+ t->points = realloc( t->points, t->nallocated * 2 * sizeof ( point* ));
t->nallocated *= 2;
}
}
@@ -228,97 +235,109 @@
t->npoints++;
}
-static void triangle_destroy(triangle* t)
+static void triangle_destroy( triangle* t )
{
- if (t->points != NULL)
- free(t->points);
- free(t);
+ if ( t->points != NULL )
+ free( t->points );
+ free( t );
}
/* Calculates barycentric coordinates of a point.
- * Takes into account that possible triangles are rectangular, with the right
- * angle at t->vertices[0], the vertices[1] vertex being in
+ * Takes into account that possible triangles are rectangular, with the right
+ * angle at t->vertices[0], the vertices[1] vertex being in
* (-3*PI/4) + (PI/2) * t->index direction from vertices[0], and
* vertices[2] being at (-5*PI/4) + (PI/2) * t->index.
*/
-static void triangle_calculatebc(triangle* t, point* p, double bc[])
+static void triangle_calculatebc( triangle* t, point* p, double bc[] )
{
double dx = p->x - t->vertices[0].x;
double dy = p->y - t->vertices[0].y;
- if (t->index == 0) {
- bc[1] = (dy - dx) / t->h;
- bc[2] = -(dx + dy) / t->h;
- } else if (t->index == 1) {
- bc[1] = (dx + dy) / t->h;
- bc[2] = (dy - dx) / t->h;
- } else if (t->index == 2) {
- bc[1] = (dx - dy) / t->h;
- bc[2] = (dx + dy) / t->h;
- } else {
- bc[1] = -(dx + dy) / t->h;
- bc[2] = (dx - dy) / t->h;
+ if ( t->index == 0 )
+ {
+ bc[1] = ( dy - dx ) / t->h;
+ bc[2] = -( dx + dy ) / t->h;
}
+ else if ( t->index == 1 )
+ {
+ bc[1] = ( dx + dy ) / t->h;
+ bc[2] = ( dy - dx ) / t->h;
+ }
+ else if ( t->index == 2 )
+ {
+ bc[1] = ( dx - dy ) / t->h;
+ bc[2] = ( dx + dy ) / t->h;
+ }
+ else
+ {
+ bc[1] = -( dx + dy ) / t->h;
+ bc[2] = ( dx - dy ) / t->h;
+ }
bc[0] = 1.0 - bc[1] - bc[2];
}
-static square* square_create(csa* parent, double xmin, double ymin, int i, int j)
+static square* square_create( csa* parent, double xmin, double ymin, int i, int j )
{
- int ii;
+ int ii;
- square* s = malloc(sizeof(square));
- double h = parent->h;
+ square * s = malloc( sizeof ( square ));
+ double h = parent->h;
s->parent = parent;
- s->i = i;
- s->j = j;
+ s->i = i;
+ s->j = j;
- s->points = NULL;
+ s->points = NULL;
s->nallocated = 0;
- s->npoints = 0;
+ s->npoints = 0;
s->primary = 0;
/*
* create 4 triangles formed by diagonals
*/
- for (ii = 0; ii < 4; ++ii) {
+ for ( ii = 0; ii < 4; ++ii )
+ {
point vertices[3];
- vertices[0].x = xmin + h / 2.0;
- vertices[0].y = ymin + h / 2.0;
- vertices[1].x = xmin + h * (((ii + 1) % 4) / 2); /* 0 1 1 0 */
- vertices[1].y = ymin + h * (((ii + 2) % 4) / 2); /* 1 1 0 0 */
- vertices[2].x = xmin + h * (ii / 2); /* 0 0 1 1 */
- vertices[2].y = ymin + h * (((ii + 1) % 4) / 2); /* 0 1 1 0 */
- s->triangles[ii] = triangle_create(s, vertices, ii);
+ vertices[0].x = xmin + h / 2.0;
+ vertices[0].y = ymin + h / 2.0;
+ vertices[1].x = xmin + h * ((( ii + 1 ) % 4 ) / 2 ); /* 0 1 1 0 */
+ vertices[1].y = ymin + h * ((( ii + 2 ) % 4 ) / 2 ); /* 1 1 0 0 */
+ vertices[2].x = xmin + h * ( ii / 2 ); /* 0 0 1 1 */
+ vertices[2].y = ymin + h * ((( ii + 1 ) % 4 ) / 2 ); /* 0 1 1 0 */
+ s->triangles[ii] = triangle_create( s, vertices, ii );
}
- for (ii = 0; ii < 25; ++ii)
+ for ( ii = 0; ii < 25; ++ii )
s->coeffs[ii] = NaN;
return s;
}
-static void square_destroy(square* s)
+static void square_destroy( square* s )
{
int i;
- for (i = 0; i < 4; ++i)
- triangle_destroy(s->triangles[i]);
- if (s->points != NULL)
- free(s->points);
- free(s);
+ for ( i = 0; i < 4; ++i )
+ triangle_destroy( s->triangles[i] );
+ if ( s->points != NULL )
+ free( s->points );
+ free( s );
}
-static void square_addpoint(square* s, point* p)
+static void square_addpoint( square* s, point* p )
{
- if (s->nallocated == s->npoints) {
- if (s->nallocated == 0) {
- s->points = malloc(NPASTART * sizeof(point*));
+ if ( s->nallocated == s->npoints )
+ {
+ if ( s->nallocated == 0 )
+ {
+ s->points = malloc( NPASTART * sizeof ( point* ));
s->nallocated = NPASTART;
- } else {
- s->points = realloc(s->points, s->nallocated * 2 * sizeof(point*));
+ }
+ else
+ {
+ s->points = realloc( s->points, s->nallocated * 2 * sizeof ( point* ));
s->nallocated *= 2;
}
}
@@ -329,104 +348,113 @@
csa* csa_create()
{
- csa* a = malloc(sizeof(csa));
+ csa* a = malloc( sizeof ( csa ));
a->xmin = DBL_MAX;
a->xmax = -DBL_MAX;
a->ymin = DBL_MAX;
a->ymax = -DBL_MAX;
- a->points = malloc(NPASTART * sizeof(point*));
+ a->points = malloc( NPASTART * sizeof ( point* ));
a->nallocated = NPASTART;
- a->npoints = 0;
+ a->npoints = 0;
- a->ni = 0;
- a->nj = 0;
- a->h = NaN;
+ a->ni = 0;
+ a->nj = 0;
+ a->h = NaN;
a->squares = NULL;
a->npt = 0;
- a->pt = NULL;
+ a->pt = NULL;
a->npmin = NPMIN_DEF;
a->npmax = NPMAX_DEF;
- a->k = K_DEF;
- a->nppc = NPPC_DEF;
+ a->k = K_DEF;
+ a->nppc = NPPC_DEF;
return a;
}
-void csa_destroy(csa* a)
+void csa_destroy( csa* a )
{
int i, j;
- if (a->squares != NULL) {
- for (j = 0; j < a->nj; ++j)
- for (i = 0; i < a->ni; ++i)
- square_destroy(a->squares[j][i]);
- free2d(a->squares);
+ if ( a->squares != NULL )
+ {
+ for ( j = 0; j < a->nj; ++j )
+ for ( i = 0; i < a->ni; ++i )
+ square_destroy( a->squares[j][i] );
+ free2d( a->squares );
}
- if (a->pt != NULL)
- free(a->pt);
- if (a->points != NULL)
- free(a->points);
- free(a);
+ if ( a->pt != NULL )
+ free( a->pt );
+ if ( a->points != NULL )
+ free( a->points );
+ free( a );
}
-void csa_addpoints(csa* a, int n, point points[])
+void csa_addpoints( csa* a, int n, point points[] )
{
int na = a->nallocated;
int i;
- assert(a->squares == NULL);
+ assert( a->squares == NULL );
/*
- * (can be called prior to squarization only)
+ * (can be called prior to squarization only)
*/
- while (na < a->npoints + n)
+ while ( na < a->npoints + n )
na *= 2;
- if (na != a->nallocated) {
- a->points = realloc(a->points, na * sizeof(point*));
+ if ( na != a->nallocated )
+ {
+ a->points = realloc( a->points, na * sizeof ( point* ));
a->nallocated = na;
}
- for (i = 0; i < n; ++i) {
+ for ( i = 0; i < n; ++i )
+ {
point* p = &points[i];
a->points[a->npoints] = p;
a->npoints++;
- if (p->x < a->xmin)
+ if ( p->x < a->xmin )
a->xmin = p->x;
- if (p->x > a->xmax)
+ if ( p->x > a->xmax )
a->xmax = p->x;
- if (p->y < a->ymin)
+ if ( p->y < a->ymin )
a->ymin = p->y;
- if (p->y > a->ymax)
+ if ( p->y > a->ymax )
a->ymax = p->y;
}
}
-/* Marks the squares containing "primary" triangles by setting "primary" flag
+/* Marks the squares containing "primary" triangles by setting "primary" flag
* to 1.
*/
-static void csa_setprimaryflag(csa* a)
+static void csa_setprimaryflag( csa* a )
{
square*** squares = a->squares;
- int nj1 = a->nj - 1;
- int ni1 = a->ni - 1;
- int i, j;
+ int nj1 = a->nj - 1;
+ int ni1 = a->ni - 1;
+ int i, j;
- for (j = 1; j < nj1; ++j) {
- for (i = 1; i < ni1; ++i) {
- if (squares[j][i]->npoints > 0) {
- if ((i + j) % 2 == 0) {
- squares[j][i]->primary = 1;
+ for ( j = 1; j < nj1; ++j )
+ {
+ for ( i = 1; i < ni1; ++i )
+ {
+ if ( squares[j][i]->npoints > 0 )
+ {
+ if (( i + j ) % 2 == 0 )
+ {
+ squares[j][i]->primary = 1;
squares[j - 1][i - 1]->primary = 1;
squares[j + 1][i - 1]->primary = 1;
squares[j - 1][i + 1]->primary = 1;
squares[j + 1][i + 1]->primary = 1;
- } else {
+ }
+ else
+ {
squares[j - 1][i]->primary = 1;
squares[j + 1][i]->primary = 1;
squares[j][i - 1]->primary = 1;
@@ -439,167 +467,182 @@
/* Splits the data domain in a number of squares.
*/
-static void csa_squarize(csa* a)
+static void csa_squarize( csa* a )
{
- int nps[7] = { 0, 0, 0, 0, 0, 0 }; /* stats on number of points per
+ int nps[7] = { 0, 0, 0, 0, 0, 0 }; /* stats on number of points per
* square */
- double dx = a->xmax - a->xmin;
- double dy = a->ymax - a->ymin;
- int npoints = a->npoints;
+ double dx = a->xmax - a->xmin;
+ double dy = a->ymax - a->ymin;
+ int npoints = a->npoints;
double h;
- int i, j, ii, nadj;
+ int i, j, ii, nadj;
- if (csa_verbose) {
- fprintf(stderr, "squarizing csa:\n");
- fflush(stderr);
+ if ( csa_verbose )
+ {
+ fprintf( stderr, "squarizing csa:\n" );
+ fflush( stderr );
}
- assert(a->squares == NULL);
+ assert( a->squares == NULL );
/*
- * (can be done only once)
+ * (can be done only once)
*/
- h = sqrt(dx * dy * a->nppc / npoints); /* square edge size */
- if (dx < h)
+ h = sqrt( dx * dy * a->nppc / npoints ); /* square edge size */
+ if ( dx < h )
h = dy * a->nppc / npoints;
- if (dy < h)
+ if ( dy < h )
h = dx * a->nppc / npoints;
a->h = h;
- a->ni = (int)(ceil(dx / h) + 2);
- a->nj = (int)(ceil(dy / h) + 2);
+ a->ni = (int) ( ceil( dx / h ) + 2 );
+ a->nj = (int) ( ceil( dy / h ) + 2 );
- if (csa_verbose) {
- fprintf(stderr, " %d x %d squares\n", a->ni, a->nj);
- fflush(stderr);
+ if ( csa_verbose )
+ {
+ fprintf( stderr, " %d x %d squares\n", a->ni, a->nj );
+ fflush( stderr );
}
/*
- * create squares
+ * create squares
*/
- a->squares = alloc2d(a->ni, a->nj, sizeof(void*));
- for (j = 0; j < a->nj; ++j)
- for (i = 0; i < a->ni; ++i)
- a->squares[j][i] = square_create(a, a->xmin + h * (i - 1), a->ymin + h * (j - 1), i, j);
+ a->squares = alloc2d( a->ni, a->nj, sizeof ( void* ));
+ for ( j = 0; j < a->nj; ++j )
+ for ( i = 0; i < a->ni; ++i )
+ a->squares[j][i] = square_create( a, a->xmin + h * ( i - 1 ), a->ymin + h * ( j - 1 ), i, j );
/*
- * map points to squares
+ * map points to squares
*/
- for (ii = 0; ii < npoints; ++ii) {
+ for ( ii = 0; ii < npoints; ++ii )
+ {
point* p = a->points[ii];
- i = (int)(floor((p->x - a->xmin) / h) + 1);
- j = (int)(floor((p->y - a->ymin) / h) + 1);
- square_addpoint(a->squares[j][i], p);
+ i = (int) ( floor(( p->x - a->xmin ) / h ) + 1 );
+ j = (int) ( floor(( p->y - a->ymin ) / h ) + 1 );
+ square_addpoint( a->squares[j][i], p );
}
/*
- * mark relevant squares with no points
+ * mark relevant squares with no points
*/
- csa_setprimaryflag(a);
+ csa_setprimaryflag( a );
/*
* Create a list of "primary" triangles, for which spline coefficients
* will be calculated directy (by least squares method), without using
- * C1 smoothness conditions.
+ * C1 smoothness conditions.
*/
- a->pt = malloc((a->ni / 2 + 1) * a->nj * sizeof(triangle*));
- for (j = 0, ii = 0, nadj = 0; j < a->nj; ++j) {
- for (i = 0; i < a->ni; ++i) {
+ a->pt = malloc(( a->ni / 2 + 1 ) * a->nj * sizeof ( triangle* ));
+ for ( j = 0, ii = 0, nadj = 0; j < a->nj; ++j )
+ {
+ for ( i = 0; i < a->ni; ++i )
+ {
square* s = a->squares[j][i];
- if (s->npoints > 0) {
+ if ( s->npoints > 0 )
+ {
int nn = s->npoints / 5;
- if (nn > 6)
+ if ( nn > 6 )
nn = 6;
nps[nn]++;
ii++;
}
- if (s->primary && s->npoints == 0)
+ if ( s->primary && s->npoints == 0 )
nadj++;
- if (s->primary) {
- a->pt[a->npt] = s->triangles[0];
+ if ( s->primary )
+ {
+ a->pt[a->npt] = s->triangles[0];
s->triangles[0]->primary = 1;
a->npt++;
}
}
}
- if (csa_verbose) {
- fprintf(stderr, " %d non-empty squares\n", ii);
- fprintf(stderr, " %d primary squares\n", a->npt);
- fprintf(stderr, " %d primary squares with no data\n", nadj);
- fprintf(stderr, " %.2f points per square \n", (double) a->npoints / ii);
+ if ( csa_verbose )
+ {
+ fprintf( stderr, " %d non-empty squares\n", ii );
+ fprintf( stderr, " %d primary squares\n", a->npt );
+ fprintf( stderr, " %d primary squares with no data\n", nadj );
+ fprintf( stderr, " %.2f points per square \n", (double) a->npoints / ii );
}
- if (csa_verbose == 2) {
- for (i = 0; i < 6; ++i)
- fprintf(stderr, " %d-%d points -- %d squares\n", i * 5, i * 5 + 4, nps[i]);
- fprintf(stderr, " %d or more points -- %d squares\n", i * 5, nps[i]);
+ if ( csa_verbose == 2 )
+ {
+ for ( i = 0; i < 6; ++i )
+ fprintf( stderr, " %d-%d points -- %d squares\n", i * 5, i * 5 + 4, nps[i] );
+ fprintf( stderr, " %d or more points -- %d squares\n", i * 5, nps[i] );
}
- if (csa_verbose == 2) {
- fprintf(stderr, " j\\i");
- for (i = 0; i < a->ni; ++i)
- fprintf(stderr, "%3d ", i);
- fprintf(stderr, "\n");
- for (j = a->nj - 1; j >= 0; --j) {
- fprintf(stderr, "%3d ", j);
- for (i = 0; i < a->ni; ++i) {
+ if ( csa_verbose == 2 )
+ {
+ fprintf( stderr, " j\\i" );
+ for ( i = 0; i < a->ni; ++i )
+ fprintf( stderr, "%3d ", i );
+ fprintf( stderr, "\n" );
+ for ( j = a->nj - 1; j >= 0; --j )
+ {
+ fprintf( stderr, "%3d ", j );
+ for ( i = 0; i < a->ni; ++i )
+ {
square* s = a->squares[j][i];
- if (s->npoints > 0)
- fprintf(stderr, "%3d ", s->npoints);
+ if ( s->npoints > 0 )
+ fprintf( stderr, "%3d ", s->npoints );
else
- fprintf(stderr, " . ");
+ fprintf( stderr, " . " );
}
- fprintf(stderr, "\n");
+ fprintf( stderr, "\n" );
}
}
- if (csa_verbose)
- fflush(stderr);
+ if ( csa_verbose )
+ fflush( stderr );
}
/* Returns all squares intersecting with a square with center in t->middle
* and edges of length 2*t->r parallel to X and Y axes.
*/
-static void getsquares(csa* a, triangle* t, int* n, square*** squares)
+static void getsquares( csa* a, triangle* t, int* n, square*** squares )
{
- int imin = (int)floor((t->middle.x - t->r - a->xmin) / t->h);
- int imax = (int)ceil((t->middle.x + t->r - a->xmin) / t->h);
- int jmin = (int)floor((t->middle.y - t->r - a->ymin) / t->h);
- int jmax = (int)ceil((t->middle.y + t->r - a->ymin) / t->h);
+ int imin = (int) floor(( t->middle.x - t->r - a->xmin ) / t->h );
+ int imax = (int) ceil(( t->middle.x + t->r - a->xmin ) / t->h );
+ int jmin = (int) floor(( t->middle.y - t->r - a->ymin ) / t->h );
+ int jmax = (int) ceil(( t->middle.y + t->r - a->ymin ) / t->h );
int i, j;
- if (imin < 0)
+ if ( imin < 0 )
imin = 0;
- if (imax >= a->ni)
+ if ( imax >= a->ni )
imax = a->ni - 1;
- if (jmin < 0)
+ if ( jmin < 0 )
jmin = 0;
- if (jmax >= a->nj)
+ if ( jmax >= a->nj )
jmax = a->nj - 1;
- *n = 0;
- (*squares) = malloc((imax - imin + 1) * (jmax - jmin + 1) * sizeof(square*));
+ *n = 0;
+ ( *squares ) = malloc(( imax - imin + 1 ) * ( jmax - jmin + 1 ) * sizeof ( square* ));
- for (j = jmin; j <= jmax; ++j) {
- for (i = imin; i <= imax; ++i) {
+ for ( j = jmin; j <= jmax; ++j )
+ {
+ for ( i = imin; i <= imax; ++i )
+ {
square* s = a->squares[j][i];
- if (s->npoints > 0) {
- (*squares)[*n] = a->squares[j][i];
- (*n)++;
+ if ( s->npoints > 0 )
+ {
+ ( *squares )[*n] = a->squares[j][i];
+ ( *n )++;
}
}
}
}
-static double distance(point* p1, point* p2)
+static double distance( point* p1, point* p2 )
{
- return hypot(p1->x - p2->x, p1->y - p2->y);
+ return hypot( p1->x - p2->x, p1->y - p2->y );
}
/* Thins data by creating an auxiliary regular grid and for leaving only
@@ -610,143 +653,160 @@
* this would require quite a bit of structural changes. So, leaving it as is
* for now.)
*/
-static void thindata(triangle* t, int npmax)
+static void thindata( triangle* t, int npmax )
{
- csa* a = t->parent->parent;
- int imax = (int)ceil(sqrt((double) (npmax * 3 / 2)));
- square*** squares = alloc2d(imax, imax, sizeof(void*));
- double h = t->r * 2.0 / imax;
- double h2 = h / 2.0;
- double xmin = t->middle.x - t->r;
- double ymin = t->middle.y - t->r;
- int i, j, ii;
+ csa * a = t->parent->parent;
+ int imax = (int) ceil( sqrt((double) ( npmax * 3 / 2 )));
+ square *** squares = alloc2d( imax, imax, sizeof ( void* ));
+ double h = t->r * 2.0 / imax;
+ double h2 = h / 2.0;
+ double xmin = t->middle.x - t->r;
+ double ymin = t->middle.y - t->r;
+ int i, j, ii;
- for (j = 0; j < imax; ++j)
- for (i = 0; i < imax; ++i)
- squares[j][i] = square_create(a, xmin + h * i, ymin + h * j, i, j);
+ for ( j = 0; j < imax; ++j )
+ for ( i = 0; i < imax; ++i )
+ squares[j][i] = square_create( a, xmin + h * i, ymin + h * j, i, j );
- for (ii = 0; ii < t->npoints; ++ii) {
- point* p = t->points[ii];
- int i = (int)floor((p->x - xmin) / h);
- int j = (int)floor((p->y - ymin) / h);
+ for ( ii = 0; ii < t->npoints; ++ii )
+ {
+ point * p = t->points[ii];
+ int i = (int) floor(( p->x - xmin ) / h );
+ int j = (int) floor(( p->y - ymin ) / h );
square* s = squares[j][i];
- if (s->npoints == 0)
- square_addpoint(s, p);
- else { /* npoints == 1 */
+ if ( s->npoints == 0 )
+ square_addpoint( s, p );
+ else /* npoints == 1 */
+ {
point pmiddle;
pmiddle.x = xmin + h * i + h2;
pmiddle.y = ymin + h * j + h2;
- if (distance(s->points[0], &pmiddle) > distance(p, &pmiddle))
+ if ( distance( s->points[0], &pmiddle ) > distance( p, &pmiddle ))
s->points[0] = p;
}
}
t->npoints = 0;
- for (j = 0; j < imax; ++j) {
- for (i = 0; i < imax; ++i) {
+ for ( j = 0; j < imax; ++j )
+ {
+ for ( i = 0; i < imax; ++i )
+ {
square* s = squares[j][i];
- if (squares[j][i]->npoints != 0)
- triangle_addpoint(t, s->points[0]);
- square_destroy(s);
+ if ( squares[j][i]->npoints != 0 )
+ triangle_addpoint( t, s->points[0] );
+ square_destroy( s );
}
}
- free2d(squares);
+ free2d( squares );
imax++;
}
-/* Finds data points to be used in calculating spline coefficients for each
+/* Finds data points to be used in calculating spline coefficients for each
* primary triangle.
*/
-static void csa_attachpoints(csa* a)
+static void csa_attachpoints( csa* a )
{
- int npmin = a->npmin;
- int npmax = a->npmax;
+ int npmin = a->npmin;
+ int npmax = a->npmax;
int nincreased = 0;
- int nthinned = 0;
+ int nthinned = 0;
int i;
- assert(a->npt > 0);
+ assert( a->npt > 0 );
- if (csa_verbose) {
- fprintf(stderr, "pre-processing data points:\n ");
- fflush(stderr);
+ if ( csa_verbose )
+ {
+ fprintf( stderr, "pre-processing data points:\n " );
+ fflush( stderr );
}
- for (i = 0; i < a->npt; ++i) {
- triangle* t = a->pt[i];
- int increased = 0;
+ for ( i = 0; i < a->npt; ++i )
+ {
+ triangle* t = a->pt[i];
+ int increased = 0;
- if (csa_verbose) {
- fprintf(stderr, ".");
- fflush(stderr);
+ if ( csa_verbose )
+ {
+ fprintf( stderr, "." );
+ fflush( stderr );
}
t->r = t->h * 1.25;
- while (1) {
- int nsquares = 0;
+ while ( 1 )
+ {
+ int nsquares = 0;
square** squares = NULL;
- int ii;
+ int ii;
- getsquares(a, t, &nsquares, &squares);
- for (ii = 0; ii < nsquares; ++ii) {
+ getsquares( a, t, &nsquares, &squares );
+ for ( ii = 0; ii < nsquares; ++ii )
+ {
square* s = squares[ii];
- int iii;
+ int iii;
- for (iii = 0; iii < s->npoints; ++iii) {
+ for ( iii = 0; iii < s->npoints; ++iii )
+ {
point* p = s->points[iii];
- if (distance(p, &t->middle) <= t->r)
- triangle_addpoint(t, p);
+ if ( distance( p, &t->middle ) <= t->r )
+ triangle_addpoint( t, p );
}
}
- free(squares);
+ free( squares );
- if (t->npoints < npmin) {
- if (!increased) {
+ if ( t->npoints < npmin )
+ {
+ if ( !increased )
+ {
increased = 1;
nincreased++;
}
- t->r *= 1.25;
+ t->r *= 1.25;
t->npoints = 0;
- } else if (t->npoints > npmax) {
+ }
+ else if ( t->npoints > npmax )
+ {
nthinned++;
- thindata(t, npmax);
- if (t->npoints > npmin)
+ thindata( t, npmax );
+ if ( t->npoints > npmin )
break;
- else {
+ else
+ {
/*
- * Sometimes you have too much data, you thin it and --
+ * Sometimes you have too much data, you thin it and --
* oops -- you have too little. This is not a frequent
- * event, so let us not bother to put a new subdivision.
+ * event, so let us not bother to put a new subdivision.
*/
- t->r *= 1.25;
+ t->r *= 1.25;
t->npoints = 0;
}
- } else
+ }
+ else
break;
}
}
- if (csa_verbose) {
- fprintf(stderr, "\n %d sets enhanced, %d sets thinned\n", nincreased, nthinned);
- fflush(stderr);
+ if ( csa_verbose )
+ {
+ fprintf( stderr, "\n %d sets enhanced, %d sets thinned\n", nincreased, nthinned );
+ fflush( stderr );
}
}
-static int n2q(int n)
+static int n2q( int n )
{
- assert(n >= 3);
+ assert( n >= 3 );
- if (n >= 10)
+ if ( n >= 10 )
return 3;
- else if (n >= 6)
+ else if ( n >= 6 )
return 2;
else /* n == 3 */
return 1;
@@ -762,144 +822,163 @@
* @param w Ouput vector that presents diagonal matrix W
* @param V output matrix V
*/
-static void svd(double** a, int n, int m, double* w, double** v)
+static void svd( double** a, int n, int m, double* w, double** v )
{
- double* rv1;
- int i, j, k, l = -1;
+ double * rv1;
+ int i, j, k, l = -1;
double tst1, c, f, g, h, s, scale;
- assert(m > 0 && n > 0);
+ assert( m > 0 && n > 0 );
- rv1 = malloc(n * sizeof(double));
+ rv1 = malloc( n * sizeof ( double ));
/*
- * householder reduction to bidiagonal form
+ * householder reduction to bidiagonal form
*/
- g = 0.0;
+ g = 0.0;
scale = 0.0;
- tst1 = 0.0;
- for (i = 0; i < n; i++) {
- l = i + 1;
+ tst1 = 0.0;
+ for ( i = 0; i < n; i++ )
+ {
+ l = i + 1;
rv1[i] = scale * g;
- g = 0.0;
- s = 0.0;
- scale = 0.0;
- if (i < m) {
- for (k = i; k < m; k++)
- scale += fabs(a[k][i]);
- if (scale != 0.0) {
- for (k = i; k < m; k++) {
+ g = 0.0;
+ s = 0.0;
+ scale = 0.0;
+ if ( i < m )
+ {
+ for ( k = i; k < m; k++ )
+ scale += fabs( a[k][i] );
+ if ( scale != 0.0 )
+ {
+ for ( k = i; k < m; k++ )
+ {
a[k][i] /= scale;
- s += a[k][i] * a[k][i];
+ s += a[k][i] * a[k][i];
}
- f = a[i][i];
- g = -copysign(sqrt(s), f);
- h = f * g - s;
+ f = a[i][i];
+ g = -copysign( sqrt( s ), f );
+ h = f * g - s;
a[i][i] = f - g;
- if (i < n - 1) {
- for (j = l; j < n; j++) {
- s = 0.0;
- for (k = i; k < m; k++)
- s += a[k][i] * a[k][j];
- f = s / h;
- for (k = i; k < m; k++)
- a[k][j] += f * a[k][i];
- }
- }
- for (k = i; k < m; k++)
+ if ( i < n - 1 )
+ {
+ for ( j = l; j < n; j++ )
+ {
+ s = 0.0;
+ for ( k = i; k < m; k++ )
+ s += a[k][i] * a[k][j];
+ f = s / h;
+ for ( k = i; k < m; k++ )
+ a[k][j] += f * a[k][i];
+ }
+ }
+ for ( k = i; k < m; k++ )
a[k][i] *= scale;
}
}
- w[i] = scale * g;
- g = 0.0;
- s = 0.0;
- scale = 0.0;
- if (i < m && i < n - 1) {
- for (k = l; k < n; k++)
- scale += fabs(a[i][k]);
- if (scale != 0.0) {
- for (k = l; k < n; k++) {
+ w[i] = scale * g;
+ g = 0.0;
+ s = 0.0;
+ scale = 0.0;
+ if ( i < m && i < n - 1 )
+ {
+ for ( k = l; k < n; k++ )
+ scale += fabs( a[i][k] );
+ if ( scale != 0.0 )
+ {
+ for ( k = l; k < n; k++ )
+ {
a[i][k] /= scale;
- s += a[i][k] * a[i][k];
+ s += a[i][k] * a[i][k];
}
- f = a[i][l];
- g = -copysign(sqrt(s), f);
- h = f * g - s;
+ f = a[i][l];
+ g = -copysign( sqrt( s ), f );
+ h = f * g - s;
a[i][l] = f - g;
- for (k = l; k < n; k++)
+ for ( k = l; k < n; k++ )
rv1[k] = a[i][k] / h;
- for (j = l; j < m; j++) {
- s = 0.0;
- for (k = l; k < n; k++)
+ for ( j = l; j < m; j++ )
+ {
+ s = 0.0;
+ for ( k = l; k < n; k++ )
s += a[j][k] * a[i][k];
- for (k = l; k < n; k++)
+ for ( k = l; k < n; k++ )
a[j][k] += s * rv1[k];
}
- for (k = l; k < n; k++)
+ for ( k = l; k < n; k++ )
a[i][k] *= scale;
}
}
{
- double tmp = fabs(w[i]) + fabs(rv1[i]);
+ double tmp = fabs( w[i] ) + fabs( rv1[i] );
- tst1 = (tst1 > tmp) ? tst1 : tmp;
+ tst1 = ( tst1 > tmp ) ? tst1 : tmp;
}
}
/*
- * accumulation of right-hand transformations
+ * accumulation of right-hand transformations
*/
- for (i = n - 1; i >= 0; i--) {
- if (i < n - 1) {
- if (g != 0.0) {
- for (j = l; j < n; j++)
+ for ( i = n - 1; i >= 0; i-- )
+ {
+ if ( i < n - 1 )
+ {
+ if ( g != 0.0 )
+ {
+ for ( j = l; j < n; j++ )
/*
- * double division avoids possible underflow
+ * double division avoids possible underflow
*/
- v[j][i] = (a[i][j] / a[i][l]) / g;
- for (j = l; j < n; j++) {
- s = 0.0;
- for (k = l; k < n; k++)
+ v[j][i] = ( a[i][j] / a[i][l] ) / g;
+ for ( j = l; j < n; j++ )
+ {
+ s = 0.0;
+ for ( k = l; k < n; k++ )
s += a[i][k] * v[k][j];
- for (k = l; k < n; k++)
+ for ( k = l; k < n; k++ )
v[k][j] += s * v[k][i];
}
}
- for (j = l; j < n; j++) {
+ for ( j = l; j < n; j++ )
+ {
v[i][j] = 0.0;
- v[j][i] = 0.0;
- }
+ v[j][i] = 0.0;
+ }
}
v[i][i] = 1.0;
- g = rv1[i];
- l = i;
+ g = rv1[i];
+ l = i;
}
/*
- * accumulation of left-hand transformations
+ * accumulation of left-hand transformations
*/
- for (i = (m < n) ? m - 1 : n - 1; i >= 0; i--) {
+ for ( i = ( m < n ) ? m - 1 : n - 1; i >= 0; i-- )
+ {
l = i + 1;
g = w[i];
- if (i != n - 1)
- for (j = l; j < n; j++)
- a[i][j] = 0.0;
- if (g != 0.0) {
- for (j = l; j < n; j++) {
- s = 0.0;
- for (k = l; k < m; k++)
+ if ( i != n - 1 )
+ for ( j = l; j < n; j++ )
+ a[i][j] = 0.0;
+ if ( g != 0.0 )
+ {
+ for ( j = l; j < n; j++ )
+ {
+ s = 0.0;
+ for ( k = l; k < m; k++ )
s += a[k][i] * a[k][j];
- /*
- * double division avoids possible underflow
- */
- f = (s / a[i][i]) / g;
- for (k = i; k < m; k++)
+ /*
+ * double division avoids possible underflow
+ */
+ f = ( s / a[i][i] ) / g;
+ for ( k = i; k < m; k++ )
a[k][j] += f * a[k][i];
}
- for (j = i; j < m; j++)
+ for ( j = i; j < m; j++ )
a[j][i] /= g;
- } else
- for (j = i; j < m; j++)
+ }
+ else
+ for ( j = i; j < m; j++ )
a[j][i] = 0.0;
a[i][i] += 1.0;
}
@@ -907,131 +986,146 @@
/*
* diagonalization of the bidiagonal form
*/
- for (k = n - 1; k >= 0; k--) {
- int k1 = k - 1;
+ for ( k = n - 1; k >= 0; k-- )
+ {
+ int k1 = k - 1;
int its = 0;
- while (1) {
- int docancellation = 1;
- double x, y, z;
- int l1 = -1;
+ while ( 1 )
+ {
+ int docancellation = 1;
+ double x, y, z;
+ int l1 = -1;
- its++;
- if (its > SVD_NMAX)
- csa_quit("svd(): no convergence in %d iterations", SVD_NMAX);
+ its++;
+ if ( its > SVD_NMAX )
+ csa_quit( "svd(): no convergence in %d iterations", SVD_NMAX );
- for (l = k; l >= 0; l--) { /* test for splitting */
- double tst2 = fabs(rv1[l]) + tst1;
+ for ( l = k; l >= 0; l-- ) /* test for splitting */
+ {
+ double tst2 = fabs( rv1[l] ) + tst1;
- if (tst2 == tst1) {
+ if ( tst2 == tst1 )
+ {
docancellation = 0;
break;
}
l1 = l - 1;
- /*
- * rv1(1) is always zero, so there is no exit through the
- * bottom of the loop
- */
- tst2 = fabs(w[l - 1]) + tst1;
- if (tst2 == tst1)
+ /*
+ * rv1(1) is always zero, so there is no exit through the
+ * bottom of the loop
+ */
+ tst2 = fabs( w[l - 1] ) + tst1;
+ if ( tst2 == tst1 )
break;
}
- /*
- * cancellation of rv1[l] if l > 1
- */
- if (docancellation) {
+ /*
+ * cancellation of rv1[l] if l > 1
+ */
+ if ( docancellation )
+ {
c = 0.0;
s = 1.0;
- for (i = l; i <= k; i++) {
- f = s * rv1[i];
+ for ( i = l; i <= k; i++ )
+ {
+ f = s * rv1[i];
rv1[i] = c * rv1[i];
- if ((fabs(f) + tst1) == tst1)
+ if (( fabs( f ) + tst1 ) == tst1 )
break;
- g = w[i];
- h = hypot(f, g);
+ g = w[i];
+ h = hypot( f, g );
w[i] = h;
- h = 1.0 / h;
- c = g * h;
- s = -f * h;
- for (j = 0; j < m; j++) {
+ h = 1.0 / h;
+ c = g * h;
+ s = -f * h;
+ for ( j = 0; j < m; j++ )
+ {
double y = a[j][l1];
double z = a[j][i];
a[j][l1] = y * c + z * s;
- a[j][i] = z * c - y * s;
+ a[j][i] = z * c - y * s;
}
}
}
- /*
- * test for convergence
- */
+ /*
+ * test for convergence
+ */
z = w[k];
- if (l != k) {
- int i1;
+ if ( l != k )
+ {
+ int i1;
- /*
- * shift from bottom 2 by 2 minor
- */
- x = w[l];
- y = w[k1];
- g = rv1[k1];
- h = rv1[k];
- f = 0.5 * (((g + z) / h) * ((g - z) / y) + y / h - h / y);
- g = hypot(f, 1.0);
- f = x - (z / x) * z + (h / x) * (y / (f + copysign(g, f)) - h);
- /*
- * next qr transformation
- */
- c = 1.0;
- s = 1.0;
- for (i1 = l; i1 < k; i1++) {
- i = i1 + 1;
- g = rv1[i];
- y = w[i];
- h = s * g;
- g = c * g;
- z = hypot(f, h);
- rv1[i1] = z;
- c = f / z;
- s = h / z;
- f = x * c + g * s;
- g = g * c - x * s;
- h = y * s;
- y *= c;
- for (j = 0; j < n; j++) {
- x = v[j][i1];
- z = v[j][i];
- v[j][i1] = x * c + z * s;
- v[j][i] = z * c - x * s;
- }
- z = hypot(f, h);
- w[i1] = z;
- /*
- * rotation can be arbitrary if z = 0
- */
- if (z != 0.0) {
- c = f / z;
- s = h / z;
- }
- f = c * g + s * y;
- x = c * y - s * g;
- for (j = 0; j < m; j++) {
- y = a[j][i1];
- z = a[j][i];
- a[j][i1] = y * c + z * s;
- a[j][i] = z * c - y * s;
- }
- }
- rv1[l] = 0.0;
- rv1[k] = f;
- w[k] = x;
- } else {
- /*
- * w[k] is made non-negative
- */
- if (z < 0.0) {
+ /*
+ * shift from bottom 2 by 2 minor
+ */
+ x = w[l];
+ y = w[k1];
+ g = rv1[k1];
+ h = rv1[k];
+ f = 0.5 * ((( g + z ) / h ) * (( g - z ) / y ) + y / h - h / y );
+ g = hypot( f, 1.0 );
+ f = x - ( z / x ) * z + ( h / x ) * ( y / ( f + copysign( g, f )) - h );
+ /*
+ * next qr transformation
+ */
+ c = 1.0;
+ s = 1.0;
+ for ( i1 = l; i1 < k; i1++ )
+ {
+ i = i1 + 1;
+ g = rv1[i];
+ y = w[i];
+ h = s * g;
+ g = c * g;
+ z = hypot( f, h );
+ rv1[i1] = z;
+ c = f / z;
+ s = h / z;
+ f = x * c + g * s;
+ g = g * c - x * s;
+ h = y * s;
+ y *= c;
+ for ( j = 0; j < n; j++ )
+ {
+ x = v[j][i1];
+ z = v[j][i];
+ v[j][i1] = x * c + z * s;
+ v[j][i] = z * c - x * s;
+ }
+ z = hypot( f, h );
+ w[i1] = z;
+ /*
+ * rotation can be arbitrary if z = 0
+ */
+ if ( z != 0.0 )
+ {
+ c = f / z;
+ s = h / z;
+ }
+ f = c * g + s * y;
+ x = c * y - s * g;
+ for ( j = 0; j < m; j++ )
+ {
+ y = a[j][i1];
+ z = a[j][i];
+ a[j][i1] = y * c + z * s;
+ a[j][i] = z * c - y * s;
+ }
+ }
+ rv1[l] = 0.0;
+ rv1[k] = f;
+ w[k] = x;
+ }
+ else
+ {
+ /*
+ * w[k] is made non-negative
+ */
+ if ( z < 0.0 )
+ {
w[k] = -z;
- for (j = 0; j < n; j++)
+ for ( j = 0; j < n; j++ )
v[j][k] = -v[j][k];
}
break;
@@ -1039,41 +1133,43 @@
}
}
- free(rv1);
+ free( rv1 );
}
/* Least squares fitting via singular value decomposition.
*/
-static void lsq(double** A, int ni, int nj, double* z, double* w, double* sol)
+static void lsq( double** A, int ni, int nj, double* z, double* w, double* sol )
{
- double** V = alloc2d(ni, ni, sizeof(double));
- double** B = alloc2d(nj, ni, sizeof(double));
- int i, j, ii;
+ double** V = alloc2d( ni, ni, sizeof ( double ));
+ double** B = alloc2d( nj, ni, sizeof ( double ));
+ int i, j, ii;
- svd(A, ni, nj, w, V);
+ svd( A, ni, nj, w, V );
- for (j = 0; j < ni; ++j)
- for (i = 0; i < ni; ++i)
+ for ( j = 0; j < ni; ++j )
+ for ( i = 0; i < ni; ++i )
V[j][i] /= w[i];
- for (i = 0; i < ni; ++i) {
+ for ( i = 0; i < ni; ++i )
+ {
double* v = V[i];
- for (j = 0; j < nj; ++j) {
+ for ( j = 0; j < nj; ++j )
+ {
double* a = A[j];
double* b = &B[i][j];
- for (ii = 0; ii < ni; ++ii)
+ for ( ii = 0; ii < ni; ++ii )
*b += v[ii] * a[ii];
}
}
- for (i = 0; i < ni; ++i)
+ for ( i = 0; i < ni; ++i )
sol[i] = 0.0;
- for (i = 0; i < ni; ++i)
- for (j = 0; j < nj; ++j)
+ for ( i = 0; i < ni; ++i )
+ for ( j = 0; j < nj; ++j )
sol[i] += B[i][j] * z[j];
- free2d(B);
- free2d(V);
+ free2d( B );
+ free2d( V );
}
/*
@@ -1093,101 +1189,109 @@
/* Calculates spline coefficients in each primary triangle by least squares
* fitting to data attached by csa_attachpoints().
*/
-static void csa_findprimarycoeffs(csa* a)
+static void csa_findprimarycoeffs( csa* a )
{
int n[4] = { 0, 0, 0, 0 };
int i;
- if (csa_verbose)
- fprintf(stderr, "calculating spline coefficients for primary triangles:\n ");
+ if ( csa_verbose )
+ fprintf( stderr, "calculating spline coefficients for primary triangles:\n " );
- for (i = 0; i < a->npt; ++i) {
- triangle* t = a->pt[i];
- int npoints = t->npoints;
- point** points = t->points;
- double* z = malloc(npoints * sizeof(double));
- int q = n2q(t->npoints);
- int ok = 1;
- double b[10];
- double b1[6];
- int ii;
+ for ( i = 0; i < a->npt; ++i )
+ {
+ triangle* t = a->pt[i];
+ int npoints = t->npoints;
+ point ** points = t->points;
+ double * z = malloc( npoints * sizeof ( double ));
+ int q = n2q( t->npoints );
+ int ok = 1;
+ double b[10];
+ double b1[6];
+ int ii;
- if (csa_verbose) {
- fprintf(stderr, ".");
- fflush(stderr);
+ if ( csa_verbose )
+ {
+ fprintf( stderr, "." );
+ fflush( stderr );
}
- for (ii = 0; ii < npoints; ++ii)
+ for ( ii = 0; ii < npoints; ++ii )
z[ii] = points[ii]->z;
- do {
+ do
+ {
double bc[3];
double wmin, wmax;
- if (!ok)
+ if ( !ok )
q--;
- assert(q >= 0);
+ assert( q >= 0 );
- if (q == 3) {
- double** A = alloc2d(10, npoints, sizeof(double));
+ if ( q == 3 )
+ {
+ double ** A = alloc2d( 10, npoints, sizeof ( double ));
double w[10];
- for (ii = 0; ii < npoints; ++ii) {
- double* aii = A[ii];
+ for ( ii = 0; ii < npoints; ++ii )
+ {
+ double * aii = A[ii];
double tmp;
- triangle_calculatebc(t, points[ii], bc);
+ triangle_calculatebc( t, points[ii], bc );
/*
- * 0 1 2 3 4 5 6 7 8 9
- * 300 210 201 120 111 102 030 021 012 003
+ * 0 1 2 3 4 5 6 7 8 9
+ * 300 210 201 120 111 102 030 021 012 003
*/
- tmp = bc[0] * bc[0];
+ tmp = bc[0] * bc[0];
aii[0] = tmp * bc[0];
- tmp *= 3.0;
+ tmp *= 3.0;
aii[1] = tmp * bc[1];
aii[2] = tmp * bc[2];
- tmp = bc[1] * bc[1];
+ tmp = bc[1] * bc[1];
aii[6] = tmp * bc[1];
- tmp *= 3.0;
+ tmp *= 3.0;
aii[3] = tmp * bc[0];
aii[7] = tmp * bc[2];
- tmp = bc[2] * bc[2];
+ tmp = bc[2] * bc[2];
aii[9] = tmp * bc[2];
- tmp *= 3.0;
+ tmp *= 3.0;
aii[5] = tmp * bc[0];
aii[8] = tmp * bc[1];
aii[4] = bc[0] * bc[1] * bc[2] * 6.0;
}
- lsq(A, 10, npoints, z, w, b);
+ lsq( A, 10, npoints, z, w, b );
wmin = w[0];
wmax = w[0];
- for (ii = 1; ii < 10; ++ii) {
- if (w[ii] < wmin)
+ for ( ii = 1; ii < 10; ++ii )
+ {
+ if ( w[ii] < wmin )
wmin = w[ii];
- else if (w[ii] > wmax)
+ else if ( w[ii] > wmax )
wmax = w[ii];
}
- if (wmin < wmax / a->k)
+ if ( wmin < wmax / a->k )
ok = 0;
- free2d(A);
-
- } else if (q == 2) {
- double** A = alloc2d(6, npoints, sizeof(double));
+ free2d( A );
+ }
+ else if ( q == 2 )
+ {
+ double ** A = alloc2d( 6, npoints, sizeof ( double ));
double w[6];
- for (ii = 0; ii < npoints; ++ii) {
+ for ( ii = 0; ii < npoints; ++ii )
+ {
double* aii = A[ii];
- triangle_calculatebc(t, points[ii], bc);
+ triangle_calculatebc( t, points[ii], bc );
/*
- * 0 1 2 3 4 5
- * 200 110 101 020 011 002
+ * 0 1 2 3 4 5
+ * 200 110 101 020 011 002
*/
aii[0] = bc[0] * bc[0];
@@ -1198,85 +1302,93 @@
aii[5] = bc[2] * bc[2];
}
- lsq(A, 6, npoints, z, w, b1);
+ lsq( A, 6, npoints, z, w, b1 );
wmin = w[0];
wmax = w[0];
- for (ii = 1; ii < 6; ++ii) {
- if (w[ii] < wmin)
+ for ( ii = 1; ii < 6; ++ii )
+ {
+ if ( w[ii] < wmin )
wmin = w[ii];
- else if (w[ii] > wmax)
+ else if ( w[ii] > wmax )
wmax = w[ii];
}
- if (wmin < wmax / a->k)
+ if ( wmin < wmax / a->k )
ok = 0;
- else { /* degree raising */
- ok = 1;
+ else /* degree raising */
+ {
+ ok = 1;
b[0] = b1[0];
- b[1] = (b1[0] + 2.0 * b1[1]) / 3.0;
- b[2] = (b1[0] + 2.0 * b1[2]) / 3.0;
- b[3] = (b1[3] + 2.0 * b1[1]) / 3.0;
- b[4] = (b1[1] + b1[2] + b1[4]) / 3.0;
- b[5] = (b1[5] + 2.0 * b1[2]) / 3.0;
+ b[1] = ( b1[0] + 2.0 * b1[1] ) / 3.0;
+ b[2] = ( b1[0] + 2.0 * b1[2] ) / 3.0;
+ b[3] = ( b1[3] + 2.0 * b1[1] ) / 3.0;
+ b[4] = ( b1[1] + b1[2] + b1[4] ) / 3.0;
+ b[5] = ( b1[5] + 2.0 * b1[2] ) / 3.0;
b[6] = b1[3];
- b[7] = (b1[3] + 2.0 * b1[4]) / 3.0;
- b[8] = (b1[5] + 2.0 * b1[4]) / 3.0;
+ b[7] = ( b1[3] + 2.0 * b1[4] ) / 3.0;
+ b[8] = ( b1[5] + 2.0 * b1[4] ) / 3.0;
b[9] = b1[5];
}
- free2d(A);
-
- } else if (q == 1) {
- double** A = alloc2d(3, npoints, sizeof(double));
+ free2d( A );
+ }
+ else if ( q == 1 )
+ {
+ double ** A = alloc2d( 3, npoints, sizeof ( double ));
double w[3];
- for (ii = 0; ii < npoints; ++ii) {
+ for ( ii = 0; ii < npoints; ++ii )
+ {
double* aii = A[ii];
- triangle_calculatebc(t, points[ii], bc);
+ triangle_calculatebc( t, points[ii], bc );
aii[0] = bc[0];
aii[1] = bc[1];
aii[2] = bc[2];
}
- lsq(A, 3, npoints, z, w, b1);
+ lsq( A, 3, npoints, z, w, b1 );
wmin = w[0];
wmax = w[0];
- for (ii = 1; ii < 3; ++ii) {
- if (w[ii] < wmin)
+ for ( ii = 1; ii < 3; ++ii )
+ {
+ if ( w[ii] < wmin )
wmin = w[ii];
- else if (w[ii] > wmax)
+ else if ( w[ii] > wmax )
wmax = w[ii];
}
- if (wmin < wmax / a->k)
+ if ( wmin < wmax / a->k )
ok = 0;
- else { /* degree raising */
- ok = 1;
+ else /* degree raising */
+ {
+ ok = 1;
b[0] = b1[0];
- b[1] = (2.0 * b1[0] + b1[1]) / 3.0;
- b[2] = (2.0 * b1[0] + b1[2]) / 3.0;
- b[3] = (2.0 * b1[1] + b1[0]) / 3.0;
- b[4] = (b1[0] + b1[1] + b1[2]) /...
[truncated message content] |
