1 /*<html><pre>  -<a                             href="index_r.htm#TOC"
2   >-------------------------------</a><a name="TOP">-</a>
3 
4    rboxlib_r.c
5      Generate input points
6 
7    notes:
8      For documentation, see prompt[] of rbox_r.c
9      50 points generated for 'rbox D4'
10 
11    WARNING:
12      incorrect range if qh_RANDOMmax is defined wrong (user_r.h)
13 */
14 
15 #include "libqhull_r.h"  /* First for user_r.h */
16 #include "random_r.h"
17 
18 #include <ctype.h>
19 #include <limits.h>
20 #include <math.h>
21 #include <setjmp.h>
22 #include <string.h>
23 #include <time.h>
24 #include <stdio.h>
25 #include <stdlib.h>
26 
27 #ifdef _MSC_VER  /* Microsoft Visual C++ */
28 #pragma warning( disable : 4706)  /* assignment within conditional expression. */
29 #pragma warning( disable : 4996)  /* this function (strncat,sprintf,strcpy) or variable may be unsafe. */
30 #endif
31 
32 #define MAXdim 200
33 #define PI 3.1415926535897932384
34 
35 /* ------------------------------ prototypes ----------------*/
36 int qh_roundi(qhT *qh, double a);
37 void qh_out1(qhT *qh, double a);
38 void qh_out2n(qhT *qh, double a, double b);
39 void qh_out3n(qhT *qh, double a, double b, double c);
40 void qh_outcoord(qhT *qh, int iscdd, double *coord, int dim);
41 void qh_outcoincident(qhT *qh, int coincidentpoints, double radius, int iscdd, double *coord, int dim);
42 void qh_rboxpoints2(qhT *qh, char* rbox_command, double **simplex);
43 
44 void    qh_fprintf_rbox(qhT *qh, FILE *fp, int msgcode, const char *fmt, ... );
45 void    qh_free(void *mem);
46 void   *qh_malloc(size_t size);
47 int     qh_rand(qhT *qh);
48 void    qh_srand(qhT *qh, int seed);
49 
50 /*-<a                             href="qh-qhull_r.htm#TOC"
51   >-------------------------------</a><a name="rboxpoints">-</a>
52 
53   qh_rboxpoints(qh, rbox_command )
54     Generate points to qh.fout according to rbox options
55     Report errors on qh.ferr
56 
57   returns:
58     0 (qh_ERRnone) on success
59     1 (qh_ERRinput) on input error
60     4 (qh_ERRmem) on memory error
61     5 (qh_ERRqhull) on internal error
62 
63   notes:
64     To avoid using stdio, redefine qh_malloc, qh_free, and qh_fprintf_rbox (user_r.c)
65     Split out qh_rboxpoints2() to avoid -Wclobbered
66 
67   design:
68     Straight line code (consider defining a struct and functions):
69 
70     Parse arguments into variables
71     Determine the number of points
72     Generate the points
73 */
qh_rboxpoints(qhT * qh,char * rbox_command)74 int qh_rboxpoints(qhT *qh, char* rbox_command) {
75   int exitcode;
76   double *simplex;
77 
78   simplex= NULL;
79   exitcode= setjmp(qh->rbox_errexit);
80   if (exitcode) {
81     /* same code for error exit and normal return.  qh.NOerrexit is set */
82     if (simplex)
83       qh_free(simplex);
84     return exitcode;
85   }
86   qh_rboxpoints2(qh, rbox_command, &simplex);
87   /* same code for error exit and normal return */
88   if (simplex)
89     qh_free(simplex);
90   return qh_ERRnone;
91 } /* rboxpoints */
92 
qh_rboxpoints2(qhT * qh,char * rbox_command,double ** simplex)93 void qh_rboxpoints2(qhT *qh, char* rbox_command, double **simplex) {
94   int i,j,k;
95   int gendim;
96   int coincidentcount=0, coincidenttotal=0, coincidentpoints=0;
97   int cubesize, diamondsize, seed=0, count, apex;
98   int dim=3, numpoints=0, totpoints, addpoints=0;
99   int issphere=0, isaxis=0,  iscdd=0, islens=0, isregular=0, iswidth=0, addcube=0;
100   int isgap=0, isspiral=0, NOcommand=0, adddiamond=0;
101   int israndom=0, istime=0;
102   int isbox=0, issimplex=0, issimplex2=0, ismesh=0;
103   double width=0.0, gap=0.0, radius=0.0, coincidentradius=0.0;
104   double coord[MAXdim], offset, meshm=3.0, meshn=4.0, meshr=5.0;
105   double *coordp, *simplexp;
106   int nthroot, mult[MAXdim];
107   double norm, factor, randr, rangap, tempr, lensangle=0, lensbase=1;
108   double anglediff, angle, x, y, cube=0.0, diamond=0.0;
109   double box= qh_DEFAULTbox; /* scale all numbers before output */
110   double randmax= qh_RANDOMmax;
111   char command[250], seedbuf[50];
112   char *s=command, *t, *first_point=NULL;
113   time_t timedata;
114 
115   *command= '\0';
116   strncat(command, rbox_command, sizeof(command)-sizeof(seedbuf)-strlen(command)-1);
117 
118   while (*s && !isspace(*s))  /* skip program name */
119     s++;
120   while (*s) {
121     while (*s && isspace(*s))
122       s++;
123     if (*s == '-')
124       s++;
125     if (!*s)
126       break;
127     if (isdigit(*s)) {
128       numpoints= qh_strtol(s, &s);
129       continue;
130     }
131     /* ============= read flags =============== */
132     switch (*s++) {
133     case 'c':
134       addcube= 1;
135       t= s;
136       while (isspace(*t))
137         t++;
138       if (*t == 'G')
139         cube= qh_strtod(++t, &s);
140       break;
141     case 'd':
142       adddiamond= 1;
143       t= s;
144       while (isspace(*t))
145         t++;
146       if (*t == 'G')
147         diamond= qh_strtod(++t, &s);
148       break;
149     case 'h':
150       iscdd= 1;
151       break;
152     case 'l':
153       isspiral= 1;
154       break;
155     case 'n':
156       NOcommand= 1;
157       break;
158     case 'r':
159       isregular= 1;
160       break;
161     case 's':
162       issphere= 1;
163       break;
164     case 't':
165       istime= 1;
166       if (isdigit(*s)) {
167         seed= qh_strtol(s, &s);
168         israndom= 0;
169       }else
170         israndom= 1;
171       break;
172     case 'x':
173       issimplex= 1;
174       break;
175     case 'y':
176       issimplex2= 1;
177       break;
178     case 'z':
179       qh->rbox_isinteger= 1;
180       break;
181     case 'B':
182       box= qh_strtod(s, &s);
183       isbox= 1;
184       break;
185     case 'C':
186       if (*s)
187         coincidentpoints=  qh_strtol(s, &s);
188       if (*s == ',') {
189         ++s;
190         coincidentradius=  qh_strtod(s, &s);
191       }
192       if (*s == ',') {
193         ++s;
194         coincidenttotal=  qh_strtol(s, &s);
195       }
196       if (*s && !isspace(*s)) {
197         qh_fprintf_rbox(qh, qh->ferr, 7080, "rbox error: arguments for 'Cn,r,m' are not 'int', 'float', and 'int'.  Remaining string is '%s'\n", s);
198         qh_errexit_rbox(qh, qh_ERRinput);
199       }
200       if (coincidentpoints==0){
201         qh_fprintf_rbox(qh, qh->ferr, 6268, "rbox error: missing arguments for 'Cn,r,m' where n is the number of coincident points, r is the radius (default 0.0), and m is the number of points\n");
202         qh_errexit_rbox(qh, qh_ERRinput);
203       }
204       if (coincidentpoints<0 || coincidenttotal<0 || coincidentradius<0.0){
205         qh_fprintf_rbox(qh, qh->ferr, 6269, "rbox error: negative arguments for 'Cn,m,r' where n (%d) is the number of coincident points, m (%d) is the number of points, and r (%.2g) is the radius (default 0.0)\n", coincidentpoints, coincidenttotal, coincidentradius);
206         qh_errexit_rbox(qh, qh_ERRinput);
207       }
208       break;
209     case 'D':
210       dim= qh_strtol(s, &s);
211       if (dim < 1
212       || dim > MAXdim) {
213         qh_fprintf_rbox(qh, qh->ferr, 6189, "rbox error: dimension, D%d, out of bounds (>=%d or <=0)\n", dim, MAXdim);
214         qh_errexit_rbox(qh, qh_ERRinput);
215       }
216       break;
217     case 'G':
218       if (isdigit(*s))
219         gap= qh_strtod(s, &s);
220       else
221         gap= 0.5;
222       isgap= 1;
223       break;
224     case 'L':
225       if (isdigit(*s))
226         radius= qh_strtod(s, &s);
227       else
228         radius= 10;
229       islens= 1;
230       break;
231     case 'M':
232       ismesh= 1;
233       if (*s)
234         meshn= qh_strtod(s, &s);
235       if (*s == ',') {
236         ++s;
237         meshm= qh_strtod(s, &s);
238       }else
239         meshm= 0.0;
240       if (*s == ',') {
241         ++s;
242         meshr= qh_strtod(s, &s);
243       }else
244         meshr= sqrt(meshn*meshn + meshm*meshm);
245       if (*s && !isspace(*s)) {
246         qh_fprintf_rbox(qh, qh->ferr, 7069, "rbox warning: assuming 'M3,4,5' since mesh args are not integers or reals\n");
247         meshn= 3.0, meshm=4.0, meshr=5.0;
248       }
249       break;
250     case 'O':
251       qh->rbox_out_offset= qh_strtod(s, &s);
252       break;
253     case 'P':
254       if (!first_point)
255         first_point= s - 1;
256       addpoints++;
257       while (*s && !isspace(*s))   /* read points later */
258         s++;
259       break;
260     case 'W':
261       width= qh_strtod(s, &s);
262       iswidth= 1;
263       break;
264     case 'Z':
265       if (isdigit(*s))
266         radius= qh_strtod(s, &s);
267       else
268         radius= 1.0;
269       isaxis= 1;
270       break;
271     default:
272       qh_fprintf_rbox(qh, qh->ferr, 6352, "rbox error: unknown flag at '%s'.\nExecute 'rbox' without arguments for documentation.\n", s - 1);
273       qh_errexit_rbox(qh, qh_ERRinput);
274     }
275     if (*s && !isspace(*s)) {
276       qh_fprintf_rbox(qh, qh->ferr, 6353, "rbox error: missing space between flags at %s.\n", s);
277       qh_errexit_rbox(qh, qh_ERRinput);
278     }
279   }
280 
281   /* ============= defaults, constants, and sizes =============== */
282   if (qh->rbox_isinteger && !isbox)
283     box= qh_DEFAULTzbox;
284   if (addcube) {
285     tempr= floor(ldexp(1.0,dim)+0.5);
286     cubesize= (int)tempr;
287     if (cube == 0.0)
288       cube= box;
289   }else
290     cubesize= 0;
291   if (adddiamond) {
292     diamondsize= 2*dim;
293     if (diamond == 0.0)
294       diamond= box;
295   }else
296     diamondsize= 0;
297   if (islens) {
298     if (isaxis) {
299         qh_fprintf_rbox(qh, qh->ferr, 6190, "rbox error: can not combine 'Ln' with 'Zn'\n");
300         qh_errexit_rbox(qh, qh_ERRinput);
301     }
302     if (radius <= 1.0) {
303         qh_fprintf_rbox(qh, qh->ferr, 6191, "rbox error: lens radius %.2g should be greater than 1.0\n",
304                radius);
305         qh_errexit_rbox(qh, qh_ERRinput);
306     }
307     lensangle= asin(1.0/radius);
308     lensbase= radius * cos(lensangle);
309   }
310 
311   if (!numpoints) {
312     if (issimplex2)
313         ; /* ok */
314     else if (isregular + issimplex + islens + issphere + isaxis + isspiral + iswidth + ismesh) {
315         qh_fprintf_rbox(qh, qh->ferr, 6192, "rbox error: missing count\n");
316         qh_errexit_rbox(qh, qh_ERRinput);
317     }else if (adddiamond + addcube + addpoints)
318         ; /* ok */
319     else {
320         numpoints= 50;  /* ./rbox D4 is the test case */
321         issphere= 1;
322     }
323   }
324   if ((issimplex + islens + isspiral + ismesh > 1)
325   || (issimplex + issphere + isspiral + ismesh > 1)) {
326     qh_fprintf_rbox(qh, qh->ferr, 6193, "rbox error: can only specify one of 'l', 's', 'x', 'Ln', or 'Mn,m,r' ('Ln s' is ok).\n");
327     qh_errexit_rbox(qh, qh_ERRinput);
328   }
329   if (coincidentpoints>0 && (numpoints == 0 || coincidenttotal > numpoints)) {
330     qh_fprintf_rbox(qh, qh->ferr, 6270, "rbox error: 'Cn,r,m' requested n coincident points for each of m points.  Either there is no points or m (%d) is greater than the number of points (%d).\n", coincidenttotal, numpoints);
331     qh_errexit_rbox(qh, qh_ERRinput);
332   }
333   if (coincidentpoints > 0 && isregular) {
334     qh_fprintf_rbox(qh, qh->ferr, 6423, "rbox error: 'Cn,r,m' is not implemented for regular points ('r')\n");
335     qh_errexit_rbox(qh, qh_ERRinput);
336   }
337 
338   if (coincidenttotal == 0)
339     coincidenttotal= numpoints;
340 
341   /* ============= print header with total points =============== */
342   if (issimplex || ismesh)
343     totpoints= numpoints;
344   else if (issimplex2)
345     totpoints= numpoints+dim+1;
346   else if (isregular) {
347     totpoints= numpoints;
348     if (dim == 2) {
349         if (islens)
350           totpoints += numpoints - 2;
351     }else if (dim == 3) {
352         if (islens)
353           totpoints += 2 * numpoints;
354       else if (isgap)
355         totpoints += 1 + numpoints;
356       else
357         totpoints += 2;
358     }
359   }else
360     totpoints= numpoints + isaxis;
361   totpoints += cubesize + diamondsize + addpoints;
362   totpoints += coincidentpoints*coincidenttotal;
363 
364   /* ============= seed randoms =============== */
365   if (istime == 0) {
366     for (s=command; *s; s++) {
367       if (issimplex2 && *s == 'y') /* make 'y' same seed as 'x' */
368         i= 'x';
369       else
370         i= *s;
371       seed= 11*seed + i;
372     }
373   }else if (israndom) {
374     seed= (int)time(&timedata);
375     sprintf(seedbuf, " t%d", seed);  /* appends an extra t, not worth removing */
376     strncat(command, seedbuf, sizeof(command) - strlen(command) - 1);
377     t= strstr(command, " t ");
378     if (t)
379       strcpy(t+1, t+3); /* remove " t " */
380   } /* else, seed explicitly set to n */
381   qh_RANDOMseed_(qh, seed);
382 
383   /* ============= print header =============== */
384 
385   if (iscdd)
386       qh_fprintf_rbox(qh, qh->fout, 9391, "%s\nbegin\n        %d %d %s\n",
387       NOcommand ? "" : command,
388       totpoints, dim+1,
389       qh->rbox_isinteger ? "integer" : "real");
390   else if (NOcommand)
391       qh_fprintf_rbox(qh, qh->fout, 9392, "%d\n%d\n", dim, totpoints);
392   else
393       /* qh_fprintf_rbox special cases 9393 to append 'command' to the RboxPoints.comment() */
394       qh_fprintf_rbox(qh, qh->fout, 9393, "%d %s\n%d\n", dim, command, totpoints);
395 
396   /* ============= explicit points =============== */
397   if ((s= first_point)) {
398     while (s && *s) { /* 'P' */
399       count= 0;
400       if (iscdd)
401         qh_out1(qh, 1.0);
402       while (*++s) {
403         qh_out1(qh, qh_strtod(s, &s));
404         count++;
405         if (isspace(*s) || !*s)
406           break;
407         if (*s != ',') {
408           qh_fprintf_rbox(qh, qh->ferr, 6194, "rbox error: missing comma after coordinate in %s\n\n", s);
409           qh_errexit_rbox(qh, qh_ERRinput);
410         }
411       }
412       if (count < dim) {
413         for (k=dim-count; k--; )
414           qh_out1(qh, 0.0);
415       }else if (count > dim) {
416         qh_fprintf_rbox(qh, qh->ferr, 6195, "rbox error: %d coordinates instead of %d coordinates in %s\n\n",
417                   count, dim, s);
418         qh_errexit_rbox(qh, qh_ERRinput);
419       }
420       qh_fprintf_rbox(qh, qh->fout, 9394, "\n");
421       while ((s= strchr(s, 'P'))) {
422         if (isspace(s[-1]))
423           break;
424       }
425     }
426   }
427 
428   /* ============= simplex distribution =============== */
429   if (issimplex+issimplex2) {
430     if (!(*simplex= (double *)qh_malloc( (size_t)(dim * (dim+1)) * sizeof(double)))) {
431       qh_fprintf_rbox(qh, qh->ferr, 6196, "rbox error: insufficient memory for simplex\n");
432       qh_errexit_rbox(qh, qh_ERRmem); /* qh_ERRmem */
433     }
434     simplexp= *simplex;
435     if (isregular) {
436       for (i=0; i<dim; i++) {
437         for (k=0; k<dim; k++)
438           *(simplexp++)= i==k ? 1.0 : 0.0;
439       }
440       for (k=0; k<dim; k++)
441         *(simplexp++)= -1.0;
442     }else {
443       for (i=0; i<dim+1; i++) {
444         for (k=0; k<dim; k++) {
445           randr= qh_RANDOMint;
446           *(simplexp++)= 2.0 * randr/randmax - 1.0;
447         }
448       }
449     }
450     if (issimplex2) {
451         simplexp= *simplex;
452       for (i=0; i<dim+1; i++) {
453         if (iscdd)
454           qh_out1(qh, 1.0);
455         for (k=0; k<dim; k++)
456           qh_out1(qh, *(simplexp++) * box);
457         qh_fprintf_rbox(qh, qh->fout, 9395, "\n");
458       }
459     }
460     for (j=0; j<numpoints; j++) {
461       if (iswidth)
462         apex= qh_RANDOMint % (dim+1);
463       else
464         apex= -1;
465       for (k=0; k<dim; k++)
466         coord[k]= 0.0;
467       norm= 0.0;
468       for (i=0; i<dim+1; i++) {
469         randr= qh_RANDOMint;
470         factor= randr/randmax;
471         if (i == apex)
472           factor *= width;
473         norm += factor;
474         for (k=0; k<dim; k++) {
475           simplexp= *simplex + i*dim + k;
476           coord[k] += factor * (*simplexp);
477         }
478       }
479       for (k=0; k<dim; k++)
480         coord[k] *= box/norm;
481       qh_outcoord(qh, iscdd, coord, dim);
482       if(coincidentcount++ < coincidenttotal)
483         qh_outcoincident(qh, coincidentpoints, coincidentradius, iscdd, coord, dim);
484     }
485     isregular= 0; /* continue with isbox */
486     numpoints= 0;
487   }
488 
489   /* ============= mesh distribution =============== */
490   if (ismesh) {
491     nthroot= (int)(pow((double)numpoints, 1.0/dim) + 0.99999);
492     for (k=dim; k--; )
493       mult[k]= 0;
494     for (i=0; i < numpoints; i++) {
495       coordp= coord;
496       for (k=0; k < dim; k++) {
497         if (k == 0)
498           *(coordp++)= mult[0] * meshn + mult[1] * (-meshm);
499         else if (k == 1)
500           *(coordp++)= mult[0] * meshm + mult[1] * meshn;
501         else
502           *(coordp++)= mult[k] * meshr;
503       }
504       qh_outcoord(qh, iscdd, coord, dim);
505       if(coincidentcount++ < coincidenttotal)
506         qh_outcoincident(qh, coincidentpoints, coincidentradius, iscdd, coord, dim);
507       for (k=0; k < dim; k++) {
508         if (++mult[k] < nthroot)
509           break;
510         mult[k]= 0;
511       }
512     }
513   }
514   /* ============= regular points for 's' =============== */
515   else if (isregular && !islens) {
516     if (dim != 2 && dim != 3) {
517       qh_fprintf_rbox(qh, qh->ferr, 6197, "rbox error: regular points can be used only in 2-d and 3-d\n\n");
518       qh_errexit_rbox(qh, qh_ERRinput);
519     }
520     if (!isaxis || radius == 0.0) {
521       isaxis= 1;
522       radius= 1.0;
523     }
524     if (dim == 3) {
525       if (iscdd)
526         qh_out1(qh, 1.0);
527       qh_out3n(qh, 0.0, 0.0, -box);
528       if (!isgap) {
529         if (iscdd)
530           qh_out1(qh, 1.0);
531         qh_out3n(qh, 0.0, 0.0, box);
532       }
533     }
534     angle= 0.0;
535     anglediff= 2.0 * PI/numpoints;
536     for (i=0; i < numpoints; i++) {
537       angle += anglediff;
538       x= radius * cos(angle);
539       y= radius * sin(angle);
540       if (dim == 2) {
541         if (iscdd)
542           qh_out1(qh, 1.0);
543         qh_out2n(qh, x*box, y*box);
544       }else {
545         norm= sqrt(1.0 + x*x + y*y);
546         if (iscdd)
547           qh_out1(qh, 1.0);
548         qh_out3n(qh, box*x/norm, box*y/norm, box/norm);
549         if (isgap) {
550           x *= 1-gap;
551           y *= 1-gap;
552           norm= sqrt(1.0 + x*x + y*y);
553           if (iscdd)
554             qh_out1(qh, 1.0);
555           qh_out3n(qh, box*x/norm, box*y/norm, box/norm);
556         }
557       }
558     }
559   }
560   /* ============= regular points for 'r Ln D2' =============== */
561   else if (isregular && islens && dim == 2) {
562     double cos_0;
563 
564     angle= lensangle;
565     anglediff= 2 * lensangle/(numpoints - 1);
566     cos_0= cos(lensangle);
567     for (i=0; i < numpoints; i++, angle -= anglediff) {
568       x= radius * sin(angle);
569       y= radius * (cos(angle) - cos_0);
570       if (iscdd)
571         qh_out1(qh, 1.0);
572       qh_out2n(qh, x*box, y*box);
573       if (i != 0 && i != numpoints - 1) {
574         if (iscdd)
575           qh_out1(qh, 1.0);
576         qh_out2n(qh, x*box, -y*box);
577       }
578     }
579   }
580   /* ============= regular points for 'r Ln D3' =============== */
581   else if (isregular && islens && dim != 2) {
582     if (dim != 3) {
583       qh_fprintf_rbox(qh, qh->ferr, 6198, "rbox error: regular points can be used only in 2-d and 3-d\n\n");
584       qh_errexit_rbox(qh, qh_ERRinput);
585     }
586     angle= 0.0;
587     anglediff= 2* PI/numpoints;
588     if (!isgap) {
589       isgap= 1;
590       gap= 0.5;
591     }
592     offset= sqrt(radius * radius - (1-gap)*(1-gap)) - lensbase;
593     for (i=0; i < numpoints; i++, angle += anglediff) {
594       x= cos(angle);
595       y= sin(angle);
596       if (iscdd)
597         qh_out1(qh, 1.0);
598       qh_out3n(qh, box*x, box*y, 0.0);
599       x *= 1-gap;
600       y *= 1-gap;
601       if (iscdd)
602         qh_out1(qh, 1.0);
603       qh_out3n(qh, box*x, box*y, box * offset);
604       if (iscdd)
605         qh_out1(qh, 1.0);
606       qh_out3n(qh, box*x, box*y, -box * offset);
607     }
608   }
609   /* ============= apex of 'Zn' distribution + gendim =============== */
610   else {
611     if (isaxis) {
612       gendim= dim-1;
613       if (iscdd)
614         qh_out1(qh, 1.0);
615       for (j=0; j < gendim; j++)
616         qh_out1(qh, 0.0);
617       qh_out1(qh, -box);
618       qh_fprintf_rbox(qh, qh->fout, 9398, "\n");
619     }else if (islens)
620       gendim= dim-1;
621     else
622       gendim= dim;
623     /* ============= generate random point in unit cube =============== */
624     for (i=0; i < numpoints; i++) {
625       norm= 0.0;
626       for (j=0; j < gendim; j++) {
627         randr= qh_RANDOMint;
628         coord[j]= 2.0 * randr/randmax - 1.0;
629         norm += coord[j] * coord[j];
630       }
631       norm= sqrt(norm);
632       /* ============= dim-1 point of 'Zn' distribution ========== */
633       if (isaxis) {
634         if (!isgap) {
635           isgap= 1;
636           gap= 1.0;
637         }
638         randr= qh_RANDOMint;
639         rangap= 1.0 - gap * randr/randmax;
640         factor= radius * rangap / norm;
641         for (j=0; j<gendim; j++)
642           coord[j]= factor * coord[j];
643       /* ============= dim-1 point of 'Ln s' distribution =========== */
644       }else if (islens && issphere) {
645         if (!isgap) {
646           isgap= 1;
647           gap= 1.0;
648         }
649         randr= qh_RANDOMint;
650         rangap= 1.0 - gap * randr/randmax;
651         factor= rangap / norm;
652         for (j=0; j<gendim; j++)
653           coord[j]= factor * coord[j];
654       /* ============= dim-1 point of 'Ln' distribution ========== */
655       }else if (islens && !issphere) {
656         if (!isgap) {
657           isgap= 1;
658           gap= 1.0;
659         }
660         j= qh_RANDOMint % gendim;
661         if (coord[j] < 0)
662           coord[j]= -1.0 - coord[j] * gap;
663         else
664           coord[j]= 1.0 - coord[j] * gap;
665       /* ============= point of 'l' distribution =============== */
666       }else if (isspiral) {
667         if (dim != 3) {
668           qh_fprintf_rbox(qh, qh->ferr, 6199, "rbox error: spiral distribution is available only in 3d\n\n");
669           qh_errexit_rbox(qh, qh_ERRinput);
670         }
671         coord[0]= cos(2*PI*i/(numpoints - 1));
672         coord[1]= sin(2*PI*i/(numpoints - 1));
673         coord[2]= 2.0*(double)i/(double)(numpoints - 1) - 1.0;
674       /* ============= point of 's' distribution =============== */
675       }else if (issphere) {
676         factor= 1.0/norm;
677         if (iswidth) {
678           randr= qh_RANDOMint;
679           factor *= 1.0 - width * randr/randmax;
680         }
681         for (j=0; j<dim; j++)
682           coord[j]= factor * coord[j];
683       }
684       /* ============= project 'Zn s' point in to sphere =============== */
685       if (isaxis && issphere) {
686         coord[dim-1]= 1.0;
687         norm= 1.0;
688         for (j=0; j<gendim; j++)
689           norm += coord[j] * coord[j];
690         norm= sqrt(norm);
691         for (j=0; j<dim; j++)
692           coord[j]= coord[j] / norm;
693         if (iswidth) {
694           randr= qh_RANDOMint;
695           coord[dim-1] *= 1 - width * randr/randmax;
696         }
697       /* ============= project 'Zn' point onto cube =============== */
698       }else if (isaxis && !issphere) {  /* not very interesting */
699         randr= qh_RANDOMint;
700         coord[dim-1]= 2.0 * randr/randmax - 1.0;
701       /* ============= project 'Ln' point out to sphere =============== */
702       }else if (islens) {
703         coord[dim-1]= lensbase;
704         for (j=0, norm= 0; j<dim; j++)
705           norm += coord[j] * coord[j];
706         norm= sqrt(norm);
707         for (j=0; j<dim; j++)
708           coord[j]= coord[j] * radius/ norm;
709         coord[dim-1] -= lensbase;
710         if (iswidth) {
711           randr= qh_RANDOMint;
712           coord[dim-1] *= 1 - width * randr/randmax;
713         }
714         if (qh_RANDOMint > randmax/2)
715           coord[dim-1]= -coord[dim-1];
716       /* ============= project 'Wn' point toward boundary =============== */
717       }else if (iswidth && !issphere) {
718         j= qh_RANDOMint % gendim;
719         if (coord[j] < 0)
720           coord[j]= -1.0 - coord[j] * width;
721         else
722           coord[j]= 1.0 - coord[j] * width;
723       }
724       /* ============= scale point to box =============== */
725       for (k=0; k<dim; k++)
726           coord[k]= coord[k] * box;
727 
728       /* ============= write output =============== */
729       qh_outcoord(qh, iscdd, coord, dim);
730       if(coincidentcount++ < coincidenttotal)
731         qh_outcoincident(qh, coincidentpoints, coincidentradius, iscdd, coord, dim);
732     }
733   }
734 
735   /* ============= write cube vertices =============== */
736   if (addcube) {
737     for (j=0; j<cubesize; j++) {
738       if (iscdd)
739         qh_out1(qh, 1.0);
740       for (k=dim-1; k>=0; k--) {
741         if (j & ( 1 << k))
742           qh_out1(qh, cube);
743         else
744           qh_out1(qh, -cube);
745       }
746       qh_fprintf_rbox(qh, qh->fout, 9400, "\n");
747     }
748   }
749 
750   /* ============= write diamond vertices =============== */
751   if (adddiamond) {
752     for (j=0; j<diamondsize; j++) {
753       if (iscdd)
754         qh_out1(qh, 1.0);
755       for (k=dim-1; k>=0; k--) {
756         if (j/2 != k)
757           qh_out1(qh, 0.0);
758         else if (j & 0x1)
759           qh_out1(qh, diamond);
760         else
761           qh_out1(qh, -diamond);
762       }
763       qh_fprintf_rbox(qh, qh->fout, 9401, "\n");
764     }
765   }
766 
767   if (iscdd)
768     qh_fprintf_rbox(qh, qh->fout, 9402, "end\nhull\n");
769 } /* rboxpoints2 */
770 
771 /*------------------------------------------------
772 outxxx - output functions for qh_rboxpoints
773 */
qh_roundi(qhT * qh,double a)774 int qh_roundi(qhT *qh, double a) {
775   if (a < 0.0) {
776     if (a - 0.5 < INT_MIN) {
777       qh_fprintf_rbox(qh, qh->ferr, 6200, "rbox input error: negative coordinate %2.2g is too large.  Reduce 'Bn'\n", a);
778       qh_errexit_rbox(qh, qh_ERRinput);
779     }
780     return (int)(a - 0.5);
781   }else {
782     if (a + 0.5 > INT_MAX) {
783       qh_fprintf_rbox(qh, qh->ferr, 6201, "rbox input error: coordinate %2.2g is too large.  Reduce 'Bn'\n", a);
784       qh_errexit_rbox(qh, qh_ERRinput);
785     }
786     return (int)(a + 0.5);
787   }
788 } /* qh_roundi */
789 
qh_out1(qhT * qh,double a)790 void qh_out1(qhT *qh, double a) {
791 
792   if (qh->rbox_isinteger)
793     qh_fprintf_rbox(qh, qh->fout, 9403, "%d ", qh_roundi(qh, a+qh->rbox_out_offset));
794   else
795     qh_fprintf_rbox(qh, qh->fout, 9404, qh_REAL_1, a+qh->rbox_out_offset);
796 } /* qh_out1 */
797 
qh_out2n(qhT * qh,double a,double b)798 void qh_out2n(qhT *qh, double a, double b) {
799 
800   if (qh->rbox_isinteger)
801     qh_fprintf_rbox(qh, qh->fout, 9405, "%d %d\n", qh_roundi(qh, a+qh->rbox_out_offset), qh_roundi(qh, b+qh->rbox_out_offset));
802   else
803     qh_fprintf_rbox(qh, qh->fout, 9406, qh_REAL_2n, a+qh->rbox_out_offset, b+qh->rbox_out_offset);
804 } /* qh_out2n */
805 
qh_out3n(qhT * qh,double a,double b,double c)806 void qh_out3n(qhT *qh, double a, double b, double c) {
807 
808   if (qh->rbox_isinteger)
809     qh_fprintf_rbox(qh, qh->fout, 9407, "%d %d %d\n", qh_roundi(qh, a+qh->rbox_out_offset), qh_roundi(qh, b+qh->rbox_out_offset), qh_roundi(qh, c+qh->rbox_out_offset));
810   else
811     qh_fprintf_rbox(qh, qh->fout, 9408, qh_REAL_3n, a+qh->rbox_out_offset, b+qh->rbox_out_offset, c+qh->rbox_out_offset);
812 } /* qh_out3n */
813 
qh_outcoord(qhT * qh,int iscdd,double * coord,int dim)814 void qh_outcoord(qhT *qh, int iscdd, double *coord, int dim) {
815     double *p= coord;
816     int k;
817 
818     if (iscdd)
819       qh_out1(qh, 1.0);
820     for (k=0; k < dim; k++)
821       qh_out1(qh, *(p++));
822     qh_fprintf_rbox(qh, qh->fout, 9396, "\n");
823 } /* qh_outcoord */
824 
qh_outcoincident(qhT * qh,int coincidentpoints,double radius,int iscdd,double * coord,int dim)825 void qh_outcoincident(qhT *qh, int coincidentpoints, double radius, int iscdd, double *coord, int dim) {
826   double *p;
827   double randr, delta;
828   int i,k;
829   double randmax= qh_RANDOMmax;
830 
831   for (i=0; i<coincidentpoints; i++) {
832     p= coord;
833     if (iscdd)
834       qh_out1(qh, 1.0);
835     for (k=0; k < dim; k++) {
836       randr= qh_RANDOMint;
837       delta= 2.0 * randr/randmax - 1.0; /* -1..+1 */
838       delta *= radius;
839       qh_out1(qh, *(p++) + delta);
840     }
841     qh_fprintf_rbox(qh, qh->fout, 9410, "\n");
842   }
843 } /* qh_outcoincident */
844 
845 /*------------------------------------------------
846    Only called from qh_rboxpoints2 or qh_fprintf_rbox
847    qh_fprintf_rbox is only called from qh_rboxpoints2
848    The largest exitcode is '255' for compatibility with exit()
849 */
qh_errexit_rbox(qhT * qh,int exitcode)850 void qh_errexit_rbox(qhT *qh, int exitcode)
851 {
852     longjmp(qh->rbox_errexit, exitcode);
853 } /* qh_errexit_rbox */
854 
855