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