rboxlib.c
Go to the documentation of this file.
00001 /*<html><pre>  -<a                             href="index.htm#TOC"
00002   >-------------------------------</a><a name="TOP">-</a>
00003 
00004    rboxlib.c
00005      Generate input points
00006 
00007    notes:
00008      For documentation, see prompt[] of rbox.c
00009      50 points generated for 'rbox D4'
00010 
00011    WARNING:
00012      incorrect range if qh_RANDOMmax is defined wrong (user.h)
00013 */
00014 
00015 #include "random.h"
00016 #include "libqhull.h"
00017 
00018 #include <ctype.h>
00019 #include <limits.h>
00020 #include <math.h>
00021 #include <setjmp.h>
00022 #include <string.h>
00023 #include <time.h>
00024 #include <stdio.h>
00025 #include <stdlib.h>
00026 
00027 #ifdef _MSC_VER  /* Microsoft Visual C++ */
00028 #pragma warning( disable : 4706)  /* assignment within conditional expression. */
00029 #pragma warning( disable : 4996)  /* this function (strncat,sprintf,strcpy) or variable may be unsafe. */
00030 #endif
00031 
00032 #define MAXdim 200
00033 #define PI 3.1415926535897932384
00034 
00035 /* ------------------------------ prototypes ----------------*/
00036 int roundi( double a);
00037 void out1( double a);
00038 void out2n( double a, double b);
00039 void out3n( double a, double b, double c);
00040 
00041 void    qh_fprintf_rbox(FILE *fp, int msgcode, const char *fmt, ... );
00042 void    qh_free(void *mem);
00043 void   *qh_malloc(size_t size);
00044 int     qh_rand( void);
00045 void    qh_srand( int seed);
00046 
00047 
00048 /* ------------------------------ globals -------------------*/
00049 
00050 /* No state is carried between rbox requests */
00051 typedef struct rboxT rboxT;
00052 struct rboxT {
00053   FILE *fout;
00054   FILE *ferr;
00055   int isinteger;
00056   double out_offset;
00057   jmp_buf errexit;        /* exit label for rboxpoints, defined by setjmp(), called by qh_errexit_rbox() */
00058 };
00059 
00060 
00061 int rbox_inuse= 0;
00062 rboxT rbox;
00063 
00064 /*-<a                             href="qh-qhull.htm#TOC"
00065   >-------------------------------</a><a name="rboxpoints">-</a>
00066 
00067   qh_rboxpoints( fout, ferr, rbox_command )
00068     Generate points to fout according to rbox options
00069     Report errors on ferr
00070 
00071   returns:
00072     0 (qh_ERRnone) on success
00073     1 (qh_ERRinput) on input error
00074     4 (qh_ERRmem) on memory error
00075     5 (qh_ERRqhull) on internal error
00076 
00077   notes:
00078     To avoid stdio, redefine qh_malloc, qh_free, and qh_fprintf_rbox (user.c)
00079     Rbox is not multithreaded.
00080 
00081   design:
00082     Straight line code (consider defining a struct and functions):
00083 
00084     Parse arguments into variables
00085     Determine the number of points
00086     Generate the points
00087 */
00088 int qh_rboxpoints(FILE* fout, FILE* ferr, char* rbox_command) {
00089   int i,j,k;
00090   int gendim;
00091   int cubesize, diamondsize, seed=0, count, apex;
00092   int dim=3 , numpoints= 0, totpoints, addpoints=0;
00093   int issphere=0, isaxis=0,  iscdd= 0, islens= 0, isregular=0, iswidth=0, addcube=0;
00094   int isgap=0, isspiral=0, NOcommand= 0, adddiamond=0;
00095   int israndom=0, istime=0;
00096   int isbox=0, issimplex=0, issimplex2=0, ismesh=0;
00097   double width=0.0, gap=0.0, radius= 0.0;
00098   double coord[MAXdim], offset, meshm=3.0, meshn=4.0, meshr=5.0;
00099   double *simplex= NULL, *simplexp;
00100   int nthroot, mult[MAXdim];
00101   double norm, factor, randr, rangap, lensangle= 0, lensbase= 1;
00102   double anglediff, angle, x, y, cube= 0.0, diamond= 0.0;
00103   double box= qh_DEFAULTbox; /* scale all numbers before output */
00104   double randmax= qh_RANDOMmax;
00105   char command[200], seedbuf[200];
00106   char *s= command, *t, *first_point= NULL;
00107   time_t timedata;
00108   int exitcode;
00109 
00110   if (rbox_inuse) {
00111     qh_fprintf_rbox(rbox.ferr, 6188, "rbox error: rbox in use by another process.  Please lock calls to rbox.\n");
00112     return qh_ERRqhull;
00113   }
00114   rbox_inuse= True;
00115   rbox.ferr= ferr;
00116   rbox.fout= fout;
00117 
00118   exitcode= setjmp(rbox.errexit);
00119   if (exitcode) {
00120     /* same code for error exit and normal return */
00121     if (simplex)
00122         qh_free(simplex);
00123     rbox_inuse= False;
00124     return exitcode;
00125   }
00126 
00127   *command= '\0';
00128   strncat(command, rbox_command, sizeof(command)-strlen(command)-1);
00129 
00130   while (*s && !isspace(*s))  /* skip program name */
00131     s++;
00132   while (*s) {
00133     while (*s && isspace(*s))
00134       s++;
00135     if (*s == '-')
00136       s++;
00137     if (!*s)
00138       break;
00139     if (isdigit(*s)) {
00140       numpoints= qh_strtol(s, &s);
00141       continue;
00142     }
00143     /* ============= read flags =============== */
00144     switch (*s++) {
00145     case 'c':
00146       addcube= 1;
00147       t= s;
00148       while (isspace(*t))
00149         t++;
00150       if (*t == 'G')
00151         cube= qh_strtod(++t, &s);
00152       break;
00153     case 'd':
00154       adddiamond= 1;
00155       t= s;
00156       while (isspace(*t))
00157         t++;
00158       if (*t == 'G')
00159         diamond= qh_strtod(++t, &s);
00160       break;
00161     case 'h':
00162       iscdd= 1;
00163       break;
00164     case 'l':
00165       isspiral= 1;
00166       break;
00167     case 'n':
00168       NOcommand= 1;
00169       break;
00170     case 'r':
00171       isregular= 1;
00172       break;
00173     case 's':
00174       issphere= 1;
00175       break;
00176     case 't':
00177       istime= 1;
00178       if (isdigit(*s)) {
00179         seed= qh_strtol(s, &s);
00180         israndom= 0;
00181       }else
00182         israndom= 1;
00183       break;
00184     case 'x':
00185       issimplex= 1;
00186       break;
00187     case 'y':
00188       issimplex2= 1;
00189       break;
00190     case 'z':
00191       rbox.isinteger= 1;
00192       break;
00193     case 'B':
00194       box= qh_strtod(s, &s);
00195       isbox= 1;
00196       break;
00197     case 'D':
00198       dim= qh_strtol(s, &s);
00199       if (dim < 1
00200       || dim > MAXdim) {
00201         qh_fprintf_rbox(rbox.ferr, 6189, "rbox error: dimension, D%d, out of bounds (>=%d or <=0)", dim, MAXdim);
00202         qh_errexit_rbox(qh_ERRinput);
00203       }
00204       break;
00205     case 'G':
00206       if (isdigit(*s))
00207         gap= qh_strtod(s, &s);
00208       else
00209         gap= 0.5;
00210       isgap= 1;
00211       break;
00212     case 'L':
00213       if (isdigit(*s))
00214         radius= qh_strtod(s, &s);
00215       else
00216         radius= 10;
00217       islens= 1;
00218       break;
00219     case 'M':
00220       ismesh= 1;
00221       if (*s)
00222         meshn= qh_strtod(s, &s);
00223       if (*s == ',') {
00224         ++s;
00225         meshm= qh_strtod(s, &s);
00226       }else
00227         meshm= 0.0;
00228       if (*s == ',') {
00229         ++s;
00230         meshr= qh_strtod(s, &s);
00231       }else
00232         meshr= sqrt(meshn*meshn + meshm*meshm);
00233       if (*s && !isspace(*s)) {
00234         qh_fprintf_rbox(rbox.ferr, 7069, "rbox warning: assuming 'M3,4,5' since mesh args are not integers or reals\n");
00235         meshn= 3.0, meshm=4.0, meshr=5.0;
00236       }
00237       break;
00238     case 'O':
00239       rbox.out_offset= qh_strtod(s, &s);
00240       break;
00241     case 'P':
00242       if (!first_point)
00243         first_point= s-1;
00244       addpoints++;
00245       while (*s && !isspace(*s))   /* read points later */
00246         s++;
00247       break;
00248     case 'W':
00249       width= qh_strtod(s, &s);
00250       iswidth= 1;
00251       break;
00252     case 'Z':
00253       if (isdigit(*s))
00254         radius= qh_strtod(s, &s);
00255       else
00256         radius= 1.0;
00257       isaxis= 1;
00258       break;
00259     default:
00260       qh_fprintf_rbox(rbox.ferr, 7070, "rbox error: unknown flag at %s.\nExecute 'rbox' without arguments for documentation.\n", s);
00261       qh_errexit_rbox(qh_ERRinput);
00262     }
00263     if (*s && !isspace(*s)) {
00264       qh_fprintf_rbox(rbox.ferr, 7071, "rbox error: missing space between flags at %s.\n", s);
00265       qh_errexit_rbox(qh_ERRinput);
00266     }
00267   }
00268 
00269   /* ============= defaults, constants, and sizes =============== */
00270   if (rbox.isinteger && !isbox)
00271     box= qh_DEFAULTzbox;
00272   if (addcube) {
00273     cubesize= (int)floor(ldexp(1.0,dim)+0.5);
00274     if (cube == 0.0)
00275       cube= box;
00276   }else
00277     cubesize= 0;
00278   if (adddiamond) {
00279     diamondsize= 2*dim;
00280     if (diamond == 0.0)
00281       diamond= box;
00282   }else
00283     diamondsize= 0;
00284   if (islens) {
00285     if (isaxis) {
00286         qh_fprintf_rbox(rbox.ferr, 6190, "rbox error: can not combine 'Ln' with 'Zn'\n");
00287         qh_errexit_rbox(qh_ERRinput);
00288     }
00289     if (radius <= 1.0) {
00290         qh_fprintf_rbox(rbox.ferr, 6191, "rbox error: lens radius %.2g should be greater than 1.0\n",
00291                radius);
00292         qh_errexit_rbox(qh_ERRinput);
00293     }
00294     lensangle= asin(1.0/radius);
00295     lensbase= radius * cos(lensangle);
00296   }
00297 
00298   if (!numpoints) {
00299     if (issimplex2)
00300         ; /* ok */
00301     else if (isregular + issimplex + islens + issphere + isaxis + isspiral + iswidth + ismesh) {
00302         qh_fprintf_rbox(rbox.ferr, 6192, "rbox error: missing count\n");
00303         qh_errexit_rbox(qh_ERRinput);
00304     }else if (adddiamond + addcube + addpoints)
00305         ; /* ok */
00306     else {
00307         numpoints= 50;  /* ./rbox D4 is the test case */
00308         issphere= 1;
00309     }
00310   }
00311   if ((issimplex + islens + isspiral + ismesh > 1)
00312   || (issimplex + issphere + isspiral + ismesh > 1)) {
00313     qh_fprintf_rbox(rbox.ferr, 6193, "rbox error: can only specify one of 'l', 's', 'x', 'Ln', or 'Mn,m,r' ('Ln s' is ok).\n");
00314     qh_errexit_rbox(qh_ERRinput);
00315   }
00316 
00317   /* ============= print header with total points =============== */
00318   if (issimplex || ismesh)
00319     totpoints= numpoints;
00320   else if (issimplex2)
00321     totpoints= numpoints+dim+1;
00322   else if (isregular) {
00323     totpoints= numpoints;
00324     if (dim == 2) {
00325         if (islens)
00326           totpoints += numpoints - 2;
00327     }else if (dim == 3) {
00328         if (islens)
00329           totpoints += 2 * numpoints;
00330       else if (isgap)
00331         totpoints += 1 + numpoints;
00332       else
00333         totpoints += 2;
00334     }
00335   }else
00336     totpoints= numpoints + isaxis;
00337   totpoints += cubesize + diamondsize + addpoints;
00338 
00339   /* ============= seed randoms =============== */
00340   if (istime == 0) {
00341     for (s=command; *s; s++) {
00342       if (issimplex2 && *s == 'y') /* make 'y' same seed as 'x' */
00343         i= 'x';
00344       else
00345         i= *s;
00346       seed= 11*seed + i;
00347     }
00348   }else if (israndom) {
00349     seed= (int)time(&timedata);
00350     sprintf(seedbuf, " t%d", seed);  /* appends an extra t, not worth removing */
00351     strncat(command, seedbuf, sizeof(command)-strlen(command)-1);
00352     t= strstr(command, " t ");
00353     if (t)
00354       strcpy(t+1, t+3); /* remove " t " */
00355   } /* else, seed explicitly set to n */
00356   qh_RANDOMseed_(seed);
00357 
00358   /* ============= print header =============== */
00359 
00360   if (iscdd)
00361       qh_fprintf_rbox(rbox.fout, 9391, "%s\nbegin\n        %d %d %s\n",
00362       NOcommand ? "" : command,
00363       totpoints, dim+1,
00364       rbox.isinteger ? "integer" : "real");
00365   else if (NOcommand)
00366       qh_fprintf_rbox(rbox.fout, 9392, "%d\n%d\n", dim, totpoints);
00367   else
00368       qh_fprintf_rbox(rbox.fout, 9393, "%d %s\n%d\n", dim, command, totpoints);
00369 
00370   /* ============= explicit points =============== */
00371   if ((s= first_point)) {
00372     while (s && *s) { /* 'P' */
00373       count= 0;
00374       if (iscdd)
00375         out1( 1.0);
00376       while (*++s) {
00377         out1( qh_strtod(s, &s));
00378         count++;
00379         if (isspace(*s) || !*s)
00380           break;
00381         if (*s != ',') {
00382           qh_fprintf_rbox(rbox.ferr, 6194, "rbox error: missing comma after coordinate in %s\n\n", s);
00383           qh_errexit_rbox(qh_ERRinput);
00384         }
00385       }
00386       if (count < dim) {
00387         for (k=dim-count; k--; )
00388           out1( 0.0);
00389       }else if (count > dim) {
00390         qh_fprintf_rbox(rbox.ferr, 6195, "rbox error: %d coordinates instead of %d coordinates in %s\n\n",
00391                   count, dim, s);
00392         qh_errexit_rbox(qh_ERRinput);
00393       }
00394       qh_fprintf_rbox(rbox.fout, 9394, "\n");
00395       while ((s= strchr(s, 'P'))) {
00396         if (isspace(s[-1]))
00397           break;
00398       }
00399     }
00400   }
00401 
00402   /* ============= simplex distribution =============== */
00403   if (issimplex+issimplex2) {
00404     if (!(simplex= (double*)qh_malloc( dim * (dim+1) * sizeof(double)))) {
00405       qh_fprintf_rbox(rbox.ferr, 6196, "rbox error: insufficient memory for simplex\n");
00406       qh_errexit_rbox(qh_ERRmem); /* qh_ERRmem */
00407     }
00408     simplexp= simplex;
00409     if (isregular) {
00410       for (i=0; i<dim; i++) {
00411         for (k=0; k<dim; k++)
00412           *(simplexp++)= i==k ? 1.0 : 0.0;
00413       }
00414       for (k=0; k<dim; k++)
00415         *(simplexp++)= -1.0;
00416     }else {
00417       for (i=0; i<dim+1; i++) {
00418         for (k=0; k<dim; k++) {
00419           randr= qh_RANDOMint;
00420           *(simplexp++)= 2.0 * randr/randmax - 1.0;
00421         }
00422       }
00423     }
00424     if (issimplex2) {
00425         simplexp= simplex;
00426       for (i=0; i<dim+1; i++) {
00427         if (iscdd)
00428           out1( 1.0);
00429         for (k=0; k<dim; k++)
00430           out1( *(simplexp++) * box);
00431         qh_fprintf_rbox(rbox.fout, 9395, "\n");
00432       }
00433     }
00434     for (j=0; j<numpoints; j++) {
00435       if (iswidth)
00436         apex= qh_RANDOMint % (dim+1);
00437       else
00438         apex= -1;
00439       for (k=0; k<dim; k++)
00440         coord[k]= 0.0;
00441       norm= 0.0;
00442       for (i=0; i<dim+1; i++) {
00443         randr= qh_RANDOMint;
00444         factor= randr/randmax;
00445         if (i == apex)
00446           factor *= width;
00447         norm += factor;
00448         for (k=0; k<dim; k++) {
00449           simplexp= simplex + i*dim + k;
00450           coord[k] += factor * (*simplexp);
00451         }
00452       }
00453       for (k=0; k<dim; k++)
00454         coord[k] /= norm;
00455       if (iscdd)
00456         out1( 1.0);
00457       for (k=0; k < dim; k++)
00458         out1( coord[k] * box);
00459       qh_fprintf_rbox(rbox.fout, 9396, "\n");
00460     }
00461     isregular= 0; /* continue with isbox */
00462     numpoints= 0;
00463   }
00464 
00465   /* ============= mesh distribution =============== */
00466   if (ismesh) {
00467     nthroot= (int)(pow((double)numpoints, 1.0/dim) + 0.99999);
00468     for (k=dim; k--; )
00469       mult[k]= 0;
00470     for (i=0; i < numpoints; i++) {
00471       for (k=0; k < dim; k++) {
00472         if (k == 0)
00473           out1( mult[0] * meshn + mult[1] * (-meshm));
00474         else if (k == 1)
00475           out1( mult[0] * meshm + mult[1] * meshn);
00476         else
00477           out1( mult[k] * meshr );
00478       }
00479       qh_fprintf_rbox(rbox.fout, 9397, "\n");
00480       for (k=0; k < dim; k++) {
00481         if (++mult[k] < nthroot)
00482           break;
00483         mult[k]= 0;
00484       }
00485     }
00486   }
00487   /* ============= regular points for 's' =============== */
00488   else if (isregular && !islens) {
00489     if (dim != 2 && dim != 3) {
00490       qh_fprintf_rbox(rbox.ferr, 6197, "rbox error: regular points can be used only in 2-d and 3-d\n\n");
00491       qh_errexit_rbox(qh_ERRinput);
00492     }
00493     if (!isaxis || radius == 0.0) {
00494       isaxis= 1;
00495       radius= 1.0;
00496     }
00497     if (dim == 3) {
00498       if (iscdd)
00499         out1( 1.0);
00500       out3n( 0.0, 0.0, -box);
00501       if (!isgap) {
00502         if (iscdd)
00503           out1( 1.0);
00504         out3n( 0.0, 0.0, box);
00505       }
00506     }
00507     angle= 0.0;
00508     anglediff= 2.0 * PI/numpoints;
00509     for (i=0; i < numpoints; i++) {
00510       angle += anglediff;
00511       x= radius * cos(angle);
00512       y= radius * sin(angle);
00513       if (dim == 2) {
00514         if (iscdd)
00515           out1( 1.0);
00516         out2n( x*box, y*box);
00517       }else {
00518         norm= sqrt(1.0 + x*x + y*y);
00519         if (iscdd)
00520           out1( 1.0);
00521         out3n( box*x/norm, box*y/norm, box/norm);
00522         if (isgap) {
00523           x *= 1-gap;
00524           y *= 1-gap;
00525           norm= sqrt(1.0 + x*x + y*y);
00526           if (iscdd)
00527             out1( 1.0);
00528           out3n( box*x/norm, box*y/norm, box/norm);
00529         }
00530       }
00531     }
00532   }
00533   /* ============= regular points for 'r Ln D2' =============== */
00534   else if (isregular && islens && dim == 2) {
00535     double cos_0;
00536 
00537     angle= lensangle;
00538     anglediff= 2 * lensangle/(numpoints - 1);
00539     cos_0= cos(lensangle);
00540     for (i=0; i < numpoints; i++, angle -= anglediff) {
00541       x= radius * sin(angle);
00542       y= radius * (cos(angle) - cos_0);
00543       if (iscdd)
00544         out1( 1.0);
00545       out2n( x*box, y*box);
00546       if (i != 0 && i != numpoints - 1) {
00547         if (iscdd)
00548           out1( 1.0);
00549         out2n( x*box, -y*box);
00550       }
00551     }
00552   }
00553   /* ============= regular points for 'r Ln D3' =============== */
00554   else if (isregular && islens && dim != 2) {
00555     if (dim != 3) {
00556       qh_fprintf_rbox(rbox.ferr, 6198, "rbox error: regular points can be used only in 2-d and 3-d\n\n");
00557       qh_errexit_rbox(qh_ERRinput);
00558     }
00559     angle= 0.0;
00560     anglediff= 2* PI/numpoints;
00561     if (!isgap) {
00562       isgap= 1;
00563       gap= 0.5;
00564     }
00565     offset= sqrt(radius * radius - (1-gap)*(1-gap)) - lensbase;
00566     for (i=0; i < numpoints; i++, angle += anglediff) {
00567       x= cos(angle);
00568       y= sin(angle);
00569       if (iscdd)
00570         out1( 1.0);
00571       out3n( box*x, box*y, 0.0);
00572       x *= 1-gap;
00573       y *= 1-gap;
00574       if (iscdd)
00575         out1( 1.0);
00576       out3n( box*x, box*y, box * offset);
00577       if (iscdd)
00578         out1( 1.0);
00579       out3n( box*x, box*y, -box * offset);
00580     }
00581   }
00582   /* ============= apex of 'Zn' distribution + gendim =============== */
00583   else {
00584     if (isaxis) {
00585       gendim= dim-1;
00586       if (iscdd)
00587         out1( 1.0);
00588       for (j=0; j < gendim; j++)
00589         out1( 0.0);
00590       out1( -box);
00591       qh_fprintf_rbox(rbox.fout, 9398, "\n");
00592     }else if (islens)
00593       gendim= dim-1;
00594     else
00595       gendim= dim;
00596     /* ============= generate random point in unit cube =============== */
00597     for (i=0; i < numpoints; i++) {
00598       norm= 0.0;
00599       for (j=0; j < gendim; j++) {
00600         randr= qh_RANDOMint;
00601         coord[j]= 2.0 * randr/randmax - 1.0;
00602         norm += coord[j] * coord[j];
00603       }
00604       norm= sqrt(norm);
00605       /* ============= dim-1 point of 'Zn' distribution ========== */
00606       if (isaxis) {
00607         if (!isgap) {
00608           isgap= 1;
00609           gap= 1.0;
00610         }
00611         randr= qh_RANDOMint;
00612         rangap= 1.0 - gap * randr/randmax;
00613         factor= radius * rangap / norm;
00614         for (j=0; j<gendim; j++)
00615           coord[j]= factor * coord[j];
00616       /* ============= dim-1 point of 'Ln s' distribution =========== */
00617       }else if (islens && issphere) {
00618         if (!isgap) {
00619           isgap= 1;
00620           gap= 1.0;
00621         }
00622         randr= qh_RANDOMint;
00623         rangap= 1.0 - gap * randr/randmax;
00624         factor= rangap / norm;
00625         for (j=0; j<gendim; j++)
00626           coord[j]= factor * coord[j];
00627       /* ============= dim-1 point of 'Ln' distribution ========== */
00628       }else if (islens && !issphere) {
00629         if (!isgap) {
00630           isgap= 1;
00631           gap= 1.0;
00632         }
00633         j= qh_RANDOMint % gendim;
00634         if (coord[j] < 0)
00635           coord[j]= -1.0 - coord[j] * gap;
00636         else
00637           coord[j]= 1.0 - coord[j] * gap;
00638       /* ============= point of 'l' distribution =============== */
00639       }else if (isspiral) {
00640         if (dim != 3) {
00641           qh_fprintf_rbox(rbox.ferr, 6199, "rbox error: spiral distribution is available only in 3d\n\n");
00642           longjmp(rbox.errexit,qh_ERRinput);
00643         }
00644         coord[0]= cos(2*PI*i/(numpoints - 1));
00645         coord[1]= sin(2*PI*i/(numpoints - 1));
00646         coord[2]= 2.0*(double)i/(double)(numpoints-1) - 1.0;
00647       /* ============= point of 's' distribution =============== */
00648       }else if (issphere) {
00649         factor= 1.0/norm;
00650         if (iswidth) {
00651           randr= qh_RANDOMint;
00652           factor *= 1.0 - width * randr/randmax;
00653         }
00654         for (j=0; j<dim; j++)
00655           coord[j]= factor * coord[j];
00656       }
00657       /* ============= project 'Zn s' point in to sphere =============== */
00658       if (isaxis && issphere) {
00659         coord[dim-1]= 1.0;
00660         norm= 1.0;
00661         for (j=0; j<gendim; j++)
00662           norm += coord[j] * coord[j];
00663         norm= sqrt(norm);
00664         for (j=0; j<dim; j++)
00665           coord[j]= coord[j] / norm;
00666         if (iswidth) {
00667           randr= qh_RANDOMint;
00668           coord[dim-1] *= 1 - width * randr/randmax;
00669         }
00670       /* ============= project 'Zn' point onto cube =============== */
00671       }else if (isaxis && !issphere) {  /* not very interesting */
00672         randr= qh_RANDOMint;
00673         coord[dim-1]= 2.0 * randr/randmax - 1.0;
00674       /* ============= project 'Ln' point out to sphere =============== */
00675       }else if (islens) {
00676         coord[dim-1]= lensbase;
00677         for (j=0, norm= 0; j<dim; j++)
00678           norm += coord[j] * coord[j];
00679         norm= sqrt(norm);
00680         for (j=0; j<dim; j++)
00681           coord[j]= coord[j] * radius/ norm;
00682         coord[dim-1] -= lensbase;
00683         if (iswidth) {
00684           randr= qh_RANDOMint;
00685           coord[dim-1] *= 1 - width * randr/randmax;
00686         }
00687         if (qh_RANDOMint > randmax/2)
00688           coord[dim-1]= -coord[dim-1];
00689       /* ============= project 'Wn' point toward boundary =============== */
00690       }else if (iswidth && !issphere) {
00691         j= qh_RANDOMint % gendim;
00692         if (coord[j] < 0)
00693           coord[j]= -1.0 - coord[j] * width;
00694         else
00695           coord[j]= 1.0 - coord[j] * width;
00696       }
00697       /* ============= write point =============== */
00698       if (iscdd)
00699         out1( 1.0);
00700       for (k=0; k < dim; k++)
00701         out1( coord[k] * box);
00702       qh_fprintf_rbox(rbox.fout, 9399, "\n");
00703     }
00704   }
00705 
00706   /* ============= write cube vertices =============== */
00707   if (addcube) {
00708     for (j=0; j<cubesize; j++) {
00709       if (iscdd)
00710         out1( 1.0);
00711       for (k=dim-1; k>=0; k--) {
00712         if (j & ( 1 << k))
00713           out1( cube);
00714         else
00715           out1( -cube);
00716       }
00717       qh_fprintf_rbox(rbox.fout, 9400, "\n");
00718     }
00719   }
00720 
00721   /* ============= write diamond vertices =============== */
00722   if (adddiamond) {
00723     for (j=0; j<diamondsize; j++) {
00724       if (iscdd)
00725         out1( 1.0);
00726       for (k=dim-1; k>=0; k--) {
00727         if (j/2 != k)
00728           out1( 0.0);
00729         else if (j & 0x1)
00730           out1( diamond);
00731         else
00732           out1( -diamond);
00733       }
00734       qh_fprintf_rbox(rbox.fout, 9401, "\n");
00735     }
00736   }
00737 
00738   if (iscdd)
00739     qh_fprintf_rbox(rbox.fout, 9402, "end\nhull\n");
00740 
00741   /* same code for error exit and normal return */
00742   if (simplex)
00743     qh_free(simplex);
00744   rbox_inuse= False;
00745   return qh_ERRnone;
00746 } /* rboxpoints */
00747 
00748 /*------------------------------------------------
00749 outxxx - output functions
00750 */
00751 int roundi( double a) {
00752   if (a < 0.0) {
00753     if (a - 0.5 < INT_MIN) {
00754       qh_fprintf_rbox(rbox.ferr, 6200, "rbox input error: negative coordinate %2.2g is too large.  Reduce 'Bn'\n", a);
00755       qh_errexit_rbox(qh_ERRinput);
00756     }
00757     return (int)(a - 0.5);
00758   }else {
00759     if (a + 0.5 > INT_MAX) {
00760       qh_fprintf_rbox(rbox.ferr, 6201, "rbox input error: coordinate %2.2g is too large.  Reduce 'Bn'\n", a);
00761       qh_errexit_rbox(qh_ERRinput);
00762     }
00763     return (int)(a + 0.5);
00764   }
00765 } /* roundi */
00766 
00767 void out1(double a) {
00768 
00769   if (rbox.isinteger)
00770     qh_fprintf_rbox(rbox.fout, 9403, "%d ", roundi( a+rbox.out_offset));
00771   else
00772     qh_fprintf_rbox(rbox.fout, 9404, qh_REAL_1, a+rbox.out_offset);
00773 } /* out1 */
00774 
00775 void out2n( double a, double b) {
00776 
00777   if (rbox.isinteger)
00778     qh_fprintf_rbox(rbox.fout, 9405, "%d %d\n", roundi(a+rbox.out_offset), roundi(b+rbox.out_offset));
00779   else
00780     qh_fprintf_rbox(rbox.fout, 9406, qh_REAL_2n, a+rbox.out_offset, b+rbox.out_offset);
00781 } /* out2n */
00782 
00783 void out3n( double a, double b, double c) {
00784 
00785   if (rbox.isinteger)
00786     qh_fprintf_rbox(rbox.fout, 9407, "%d %d %d\n", roundi(a+rbox.out_offset), roundi(b+rbox.out_offset), roundi(c+rbox.out_offset));
00787   else
00788     qh_fprintf_rbox(rbox.fout, 9408, qh_REAL_3n, a+rbox.out_offset, b+rbox.out_offset, c+rbox.out_offset);
00789 } /* out3n */
00790 
00791 void qh_errexit_rbox(int exitcode)
00792 {
00793     longjmp(rbox.errexit, exitcode);
00794 } /* rbox_errexit */
00795 


libqhull-ours
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:11