/* ---------------------------------
rbox.c
Generate input points for qhull.
notes:
50 points generated for 'rbox D4'
This code needs a full rewrite. It needs separate procedures for each
distribution with common, helper procedures.
WARNING:
incorrect range if qh_RANDOMmax is defined wrong (user.h)
*/
#include
#include
#include
#include
#include
#include
#include
#include "user.h"
#if __MWERKS__ && __POWERPC__
#include
#include
#include
#include
#endif
#ifdef _MSC_VER /* Microsoft Visual C++ */
#pragma warning( disable : 4244) /* conversion from double to int */
#endif
#define MINVALUE 0.8
#define MAXdim 200
#define PI 3.1415926535897932384
#define DEFAULTzbox 1e6
char prompt[]= "\n\
-rbox- generate various point distributions. Default is random in cube.\n\
\n\
args (any order, space separated): Version: 2001/06/24\n\
3000 number of random points in cube, lens, spiral, sphere or grid\n\
D3 dimension 3-d\n\
c add a unit cube to the output ('c G2.0' sets size)\n\
d add a unit diamond to the output ('d G2.0' sets size)\n\
l generate a regular 3-d spiral\n\
r generate a regular polygon, ('r s Z1 G0.1' makes a cone)\n\
s generate cospherical points\n\
x generate random points in simplex, may use 'r' or 'Wn'\n\
y same as 'x', plus simplex\n\
Pn,m,r add point [n,m,r] first, pads with 0\n\
\n\
Ln lens distribution of radius n. Also 's', 'r', 'G', 'W'.\n\
Mn,m,r lattice (Mesh) rotated by [n,-m,0], [m,n,0], [0,0,r], ...\n\
'27 M1,0,1' is {0,1,2} x {0,1,2} x {0,1,2}. Try 'M3,4 z'.\n\
W0.1 random distribution within 0.1 of the cube's or sphere's surface\n\
Z0.5 s random points in a 0.5 disk projected to a sphere\n\
Z0.5 s G0.6 same as Z0.5 within a 0.6 gap\n\
\n\
Bn bounding box coordinates, default %2.2g\n\
h output as homogeneous coordinates for cdd\n\
n remove command line from the first line of output\n\
On offset coordinates by n\n\
t use time as the random number seed (default is command line)\n\
tn use n as the random number seed\n\
z print integer coordinates, default 'Bn' is %2.2g\n\
";
/* ------------------------------ prototypes ----------------*/
int roundi( double a);
void out1( double a);
void out2n( double a, double b);
void out3n( double a, double b, double c);
int qh_rand( void);
void qh_srand( int seed);
/* ------------------------------ globals -------------------*/
FILE *fp;
int isinteger= 0;
double out_offset= 0.0;
/*--------------------------------------------
-rbox- main procedure of rbox application
*/
int main(int argc, char **argv) {
int i,j,k;
int gendim;
int cubesize, diamondsize, seed=0, count, apex;
int dim=3 , numpoints= 0, totpoints, addpoints=0;
int issphere=0, isaxis=0, iscdd= 0, islens= 0, isregular=0, iswidth=0, addcube=0;
int isgap=0, isspiral=0, NOcommand= 0, adddiamond=0, istime=0;
int isbox=0, issimplex=0, issimplex2=0, ismesh=0;
double width=0.0, gap=0.0, radius= 0.0;
double coord[MAXdim], offset, meshm=3.0, meshn=4.0, meshr=5.0;
double *simplex, *simplexp;
int nthroot, mult[MAXdim];
double norm, factor, randr, rangap, lensangle= 0, lensbase= 1;
double anglediff, angle, x, y, cube= 0.0, diamond= 0.0;
double box= qh_DEFAULTbox; /* scale all numbers before output */
double randmax= qh_RANDOMmax;
char command[200], *s, seedbuf[200];
time_t timedata;
#if __MWERKS__ && __POWERPC__
char inBuf[BUFSIZ], outBuf[BUFSIZ], errBuf[BUFSIZ];
SIOUXSettings.showstatusline= False;
SIOUXSettings.tabspaces= 1;
SIOUXSettings.rows= 40;
if (setvbuf (stdin, inBuf, _IOFBF, sizeof(inBuf)) < 0 /* w/o, SIOUX I/O is slow*/
|| setvbuf (stdout, outBuf, _IOFBF, sizeof(outBuf)) < 0
|| (stdout != stderr && setvbuf (stderr, errBuf, _IOFBF, sizeof(errBuf)) < 0))
fprintf ( stderr, "qhull internal warning (main): could not change stdio to fully buffered.\n");
argc= ccommand(&argv);
#endif
if (argc == 1) {
printf (prompt, box, DEFAULTzbox);
exit(1);
}
if ((s = strrchr( argv[0], '\\'))) /* Borland gives full path */
strcpy (command, s+1);
else
strcpy (command, argv[0]);
if ((s= strstr (command, ".EXE"))
|| (s= strstr (command, ".exe")))
*s= '\0';
/* ============= read flags =============== */
for (i=1; i < argc; i++) {
if (strlen (command) + strlen(argv[i]) + 1 < sizeof(command) ) {
strcat (command, " ");
strcat (command, argv[i]);
}
if (isdigit (argv[i][0])) {
numpoints= atoi (argv[i]);
continue;
}
if (argv[i][0] == '-')
(argv[i])++;
switch (argv[i][0]) {
case 'c':
addcube= 1;
if (i+1 < argc && argv[i+1][0] == 'G')
cube= (double) atof (&argv[++i][1]);
break;
case 'd':
adddiamond= 1;
if (i+1 < argc && argv[i+1][0] == 'G')
diamond= (double) atof (&argv[++i][1]);
break;
case 'h':
iscdd= 1;
break;
case 'l':
isspiral= 1;
break;
case 'n':
NOcommand= 1;
break;
case 'r':
isregular= 1;
break;
case 's':
issphere= 1;
break;
case 't':
istime= 1;
if (isdigit (argv[i][1]))
seed= atoi (&argv[i][1]);
else {
seed= time (&timedata);
sprintf (seedbuf, "%d", seed);
strcat (command, seedbuf);
}
break;
case 'x':
issimplex= 1;
break;
case 'y':
issimplex2= 1;
break;
case 'z':
isinteger= 1;
break;
case 'B':
box= (double) atof (&argv[i][1]);
isbox= 1;
break;
case 'D':
dim= atoi (&argv[i][1]);
if (dim < 1
|| dim > MAXdim) {
fprintf (stderr, "rbox error: dim %d too large or too small\n", dim);
exit (1);
}
break;
case 'G':
if (argv[i][1])
gap= (double) atof (&argv[i][1]);
else
gap= 0.5;
isgap= 1;
break;
case 'L':
if (argv[i][1])
radius= (double) atof (&argv[i][1]);
else
radius= 10;
islens= 1;
break;
case 'M':
ismesh= 1;
s= argv[i]+1;
if (*s)
meshn= strtod (s, &s);
if (*s == ',')
meshm= strtod (++s, &s);
else
meshm= 0.0;
if (*s == ',')
meshr= strtod (++s, &s);
else
meshr= sqrt (meshn*meshn + meshm*meshm);
if (*s) {
fprintf (stderr, "rbox warning: assuming 'M3,4,5' since mesh args are not integers or reals\n");
meshn= 3.0, meshm=4.0, meshr=5.0;
}
break;
case 'O':
out_offset= (double) atof (&argv[i][1]);
break;
case 'P':
addpoints++;
break;
case 'W':
width= (double) atof (&argv[i][1]);
iswidth= 1;
break;
case 'Z':
if (argv[i][1])
radius= (double) atof (&argv[i][1]);
else
radius= 1.0;
isaxis= 1;
break;
default:
fprintf (stderr, "rbox warning: unknown flag %s.\nExecute 'rbox' without arguments for documentation.\n", argv[i]);
}
}
/* ============= defaults, constants, and sizes =============== */
if (isinteger && !isbox)
box= DEFAULTzbox;
if (addcube) {
cubesize= floor(ldexp(1.0,dim)+0.5);
if (cube == 0.0)
cube= box;
}else
cubesize= 0;
if (adddiamond) {
diamondsize= 2*dim;
if (diamond == 0.0)
diamond= box;
}else
diamondsize= 0;
if (islens) {
if (isaxis) {
fprintf (stderr, "rbox error: can not combine 'Ln' with 'Zn'\n");
exit(1);
}
if (radius <= 1.0) {
fprintf (stderr, "rbox error: lens radius %.2g should be greater than 1.0\n",
radius);
exit(1);
}
lensangle= asin (1.0/radius);
lensbase= radius * cos (lensangle);
}
if (!numpoints) {
if (issimplex2)
; /* ok */
else if (isregular + issimplex + islens + issphere + isaxis + isspiral + iswidth + ismesh) {
fprintf (stderr, "rbox error: missing count\n");
exit(1);
}else if (adddiamond + addcube + addpoints)
; /* ok */
else {
numpoints= 50; /* ./rbox D4 is the test case */
issphere= 1;
}
}
if ((issimplex + islens + isspiral + ismesh > 1)
|| (issimplex + issphere + isspiral + ismesh > 1)) {
fprintf (stderr, "rbox error: can only specify one of 'l', 's', 'x', 'Ln', or 'Mn,m,r' ('Ln s' is ok).\n");
exit(1);
}
fp= stdout;
/* ============= print header with total points =============== */
if (issimplex || ismesh)
totpoints= numpoints;
else if (issimplex2)
totpoints= numpoints+dim+1;
else if (isregular) {
totpoints= numpoints;
if (dim == 2) {
if (islens)
totpoints += numpoints - 2;
}else if (dim == 3) {
if (islens)
totpoints += 2 * numpoints;
else if (isgap)
totpoints += 1 + numpoints;
else
totpoints += 2;
}
}else
totpoints= numpoints + isaxis;
totpoints += cubesize + diamondsize + addpoints;
if (iscdd)
fprintf(fp, "%s\nbegin\n %d %d %s\n",
NOcommand ? "" : command,
totpoints, dim+1,
isinteger ? "integer" : "real");
else if (NOcommand)
fprintf(fp, "%d\n%d\n", dim, totpoints);
else
fprintf(fp, "%d %s\n%d\n", dim, command, totpoints);
/* ============= seed randoms =============== */
if (istime == 0) {
for (s=command; *s; s++) {
if (issimplex2 && *s == 'y') /* make 'y' same seed as 'x' */
i= 'x';
else
i= *s;
seed= 11*seed + i;
}
} /* else, seed explicitly set to n or to time */
qh_RANDOMseed_(seed);
/* ============= explicit points =============== */
for (i=1; i < argc; i++) {
if (argv[i][0] == 'P') {
s= argv[i]+1;
count= 0;
if (iscdd)
out1( 1.0);
while (*s) {
out1( strtod (s, &s));
count++;
if (*s) {
if (*s++ != ',') {
fprintf (stderr, "rbox error: missing comma after coordinate in %s\n\n", argv[i]);
exit (1);
}
}
}
if (count < dim) {
for (k= dim-count; k--; )
out1( 0.0);
}else if (count > dim) {
fprintf (stderr, "rbox error: %d coordinates instead of %d coordinates in %s\n\n",
count, dim, argv[i]);
exit (1);
}
fprintf (fp, "\n");
}
}
/* ============= simplex distribution =============== */
if (issimplex+issimplex2) {
if (!(simplex= malloc( dim * (dim+1) * sizeof(double)))) {
fprintf (stderr, "insufficient memory for simplex\n");
exit(0);
}
simplexp= simplex;
if (isregular) {
for (i= 0; i randmax/2)
coord[dim-1]= -coord[dim-1];
/* ============= project 'Wn' point toward boundary =============== */
}else if (iswidth && !issphere) {
j= qh_RANDOMint % gendim;
if (coord[j] < 0)
coord[j]= -1.0 - coord[j] * width;
else
coord[j]= 1.0 - coord[j] * width;
}
/* ============= write point =============== */
if (iscdd)
out1( 1.0);
for (k=0; k < dim; k++)
out1( coord[k] * box);
fprintf (fp, "\n");
}
}
/* ============= write cube vertices =============== */
if (addcube) {
for (j=0; j=0; k--) {
if (j & ( 1 << k))
out1( cube);
else
out1( -cube);
}
fprintf (fp, "\n");
}
}
/* ============= write diamond vertices =============== */
if (adddiamond) {
for (j=0; j=0; k--) {
if (j/2 != k)
out1( 0.0);
else if (j & 0x1)
out1( diamond);
else
out1( -diamond);
}
fprintf (fp, "\n");
}
}
if (iscdd)
fprintf (fp, "end\nhull\n");
return 0;
} /* rbox */
/*------------------------------------------------
-outxxx - output functions
*/
int roundi( double a) {
if (a < 0.0) {
if (a - 0.5 < INT_MIN) {
fprintf(stderr, "rbox input error: coordinate %2.2g is too large. Reduce 'Bn'\n", a);
exit (1);
}
return a - 0.5;
}else {
if (a + 0.5 > INT_MAX) {
fprintf(stderr, "rbox input error: coordinate %2.2g is too large. Reduce 'Bn'\n", a);
exit (1);
}
return a + 0.5;
}
} /* roundi */
void out1(double a) {
if (isinteger)
fprintf(fp, "%d ", roundi( a+out_offset));
else
fprintf(fp, qh_REAL_1, a+out_offset);
} /* out1 */
void out2n( double a, double b) {
if (isinteger)
fprintf(fp, "%d %d\n", roundi(a+out_offset), roundi(b+out_offset));
else
fprintf(fp, qh_REAL_2n, a+out_offset, b+out_offset);
} /* out2n */
void out3n( double a, double b, double c) {
if (isinteger)
fprintf(fp, "%d %d %d\n", roundi(a+out_offset), roundi(b+out_offset), roundi(c+out_offset));
else
fprintf(fp, qh_REAL_3n, a+out_offset, b+out_offset, c+out_offset);
} /* out3n */
/*-------------------------------------------------
-rand & srand- generate pseudo-random number between 1 and 2^31 -2
from Park & Miller's minimimal standard random number generator
Communications of the ACM, 31:1192-1201, 1988.
notes:
does not use 0 or 2^31 -1
this is silently enforced by qh_srand()
copied from geom2.c
*/
static int seed = 1; /* global static */
int qh_rand( void) {
#define qh_rand_a 16807
#define qh_rand_m 2147483647
#define qh_rand_q 127773 /* m div a */
#define qh_rand_r 2836 /* m mod a */
int lo, hi, test;
hi = seed / qh_rand_q; /* seed div q */
lo = seed % qh_rand_q; /* seed mod q */
test = qh_rand_a * lo - qh_rand_r * hi;
if (test > 0)
seed= test;
else
seed= test + qh_rand_m;
return seed;
} /* rand */
void qh_srand( int newseed) {
if (newseed < 1)
seed= 1;
else if (newseed >= qh_rand_m)
seed= qh_rand_m - 1;
else
seed= newseed;
} /* qh_srand */